等位基因的特定表达（allele specific expression，ASE）数据是后续分多分析的基础。利用它可以识别印记基因、研究蛋白截断变异以及用于量化顺式（cis）调控遗传变异的影响。虽然下游的分析具有多种形式，但是第一步总是要获得ASE数据。在最基本的层面上，ASE数据由RNA-seq读到的杂合子SNP位点生成，根据它们的序列可以对reads进行分型。但是一般情况下，不使用indel型变异，因为indel周围干扰数据的正确比对。

    这里整理工具phASER（phasing and ASE fromRNA-seq）的具体用法，这个软件可以获得等位基因全部单倍型的可靠表达水平。软件的作者认识 这是产生ASE数据的黄金标准，并且获得的数据可以用于任何感兴趣的下游分析。

1. 下载和安装phASER

软件下载自 https://github.com/secastel/phaser.git

phASER 运行于Python 2.7.x 并依赖下列包: SciPy, NumPy, samtools, tabix, bedtools, Cython

随后对软件进行解压

unzip phaser-master.zip

cd phaser-master/phaser

python setup.py build_ext --inplace

cd ../..

2. 准备分析数据

phASER 需要VCF文件和比对的bam文件。

注意：phASER是利用VCF文件中的杂合位点计算等位基因表达水平。

首先VCF需要压缩和建立索引

bgzip input.vcf

tabix -p vcf input.vcf.gz

bam文件需要建立索引

samtools index reads.bam

3. 运行phASER

python phaser-master/phaser/phaser.py --vcf input.vcf.gz --bam reads.bam --paired_end 1 --mapq 255/60 --baseq 10 --sample sampleID --threads N --temp_dir ./tmp --id_separator ';' --o prefix_output --output_read_ids 1

如果bam文件是STAR产生的，mapq选255；如果bam文件是hisat2产生的，mapq选60。sampleID由call SNP之前的AddOrReplaceReadGroups 决定。

--threads N 线程数，但是phASER对多线程优化不好，使用多线程以后，很容易出报错。

prefix_output是输出文件的前缀，需要注意的是输出文件较多，为了防止文件混乱，最好每个样本输出到一个目录下。

--temp_dir ./tmp 为了防止系统没有空间储存临时文件，最好指定临时文件的位置。

--id_separator默认的分隔符是下划线 ‘_’ ， 如果基因组scaffold的命子中包含下划线，这里要重新指定。

--output_read_ids 输出单倍型的reads ID，这个参数默认是0。reads ID的结果在haplotypic_counts.txt 这个文件里。

如果软件运行成功，最后会有类似信息

COMPLETED using 1591165 reads in 595 seconds using 8 threads

4. 定量单倍型表达

python phasr-master/phaser_gene_ae/phaser_gene_ae.py --haplotypic_counts prefix.haplotypic_counts.txt --features gene.bed --id_separator ';' --o prefix.gene_ae.txt

--haplotypic_counts 入读的文件是上步产生的 prefix.haplotypic_counts.txt

--features 为基因bed注释，这里要求要有基因名，染色体名和之前bam、vcf中的一定要一致。

5. 单倍型表达文件说明

结果文件如下：

contig – Gene chromosome.

start – Gene start (0 based).

stop – Gene stop (0 based).

name – Gene name.

aCount – Total allelic count for haplotype A.

bCount – Total allelic count for haplotype B.

totalCount – Total allelic coverage of this feature (aCount + bCount).

log2_aFC – Effect size for the allelic imbalance reported as allelic fold change (log2(aCount/bCount)) defined in our paper.

n_variants – Number of variants with allelic data in this feature.

variants – List of variants with allelic data in this feature (contig_position_ref_alt).

aCount 和 bCount 列是唯一比对到单倍型的reads数，可用于下游分析。

aCount 对应vcf中ref单倍型，bCount对应vcf中alt单倍型。