HTseq使用python编写的一款用于基因Count表达量分析的软件，可根据SAM/BAM 以及gff文件得到基因水平的Count表达量
安装
wget https://pypi.python.org/packages/46/f7/6105848893b1d280692eac4f4f3c08ed7f424cec636aeda66b50bbcf217e/HTSeq-0.7.2.tar.gz
tar zxf HTSeq-0.7.2.tar.gz
cd HTSeq-0.7.2
/opt/sysoft/Python-2.7.11/bin/python setup.py build
/opt/sysoft/Python-2.7.11/bin/python setup.py install
cd ../ && rm HTSeq-0.7.2 -rf
#或者
conda install htseq
HTseq统计Count的模式
共有3种，union, intersection_strict, intersection_nonempty
使用流程
htseq-count -f bam -r name -s no -a 10 \
-t exon -i ID -m union .bam GFF \
> count.txt
#参数
-f 输入格式bam/sam
-r 对sam或者bam的排序方式。默认为name。HTseq推介使用name，一般比对软件也都是按name进行排序
-s stranded；默认为yes。一般选择no， yes是成对的reads数量
-a 默认10；忽略比对质量低于此值的比对结果
-t 指定gff第三列的特征进行统计；默认为exon；mNRA选择exon
-i 指定gff第九列的特征用于输出；默认为gene_id，根据gff进行修改
-m 比对模式；对于原核生物，推荐使用intersection-strict模式；对于真核生物，推荐使用union模式。
-o 输入一个sam文件；该sam文件的比对结果中多了一个XF标签，表示该read比对到了某个feature上。
-q 不输出程序运行的状态信息和警告信息。
结果：得到一txt的reads统计文件共两列
# 第一列基因名；第二列reads数量
g1 0
g10 1311
g100 523
g1000 0
g1001 242
g1002 0
g1003 41
g1004 8453
若有多个bam文件，可以使用如下命令将不同样本的reads数量进行合并
paste *.txt |awk '{printf $1"\t";for(i=2;i<=ALL;i+=2)printf $i"\t";printf $i"\n"}' >Count_all.txt
# ALL 为一数值，其值等于txt文件的2倍
==小提示==
1. HTSeq是对有参考基因组的转录组测序数据进行表达量分析的，其输入文件必须有SAM和GTF文件。
2. 一般情况下HTSeq得到的Counts结果会用于下一步不同样品间的基因表达量差异分析，而不是一个样品内部基因的表达量比较。因此，HTSeq设置了-a参数的默认值10，来忽略掉比对到多个位置的reads信息，其结果有利于后续的差异分析。
3. 输入的GTF文件中不能包含可变剪接信息，否则HTSeq会认为每个可变剪接都是单独的基因，导致能比对到多个可变剪接转录本上的reads的计算结果是ambiguous，从而不能计算到基因的count中。即使设置-i参数的值为transcript_id，其结果一样是不准确的，只是得到transcripts的表达量。
网友评论