htseq-count的使用

转：https://blog.youkuaiyun.com/herokoking/article/details/78257714

没找到原文地址，上面的博主也是转的别人的，但是没有标出原文出处

htseq-count 是一款用于reads计数的软件，他能对位于基因组上的一些单位的reads数进行统计，这里所说的单位主要是指染色体上的一组位置区间（我们常见的就是gene exon）。

＊教程基于0.6.1p2版本

基本用法（软件安装方法省略）：

htseq-count [options] <alignment_file> <gff_file>

输入：

1. <alignment_file>:   比对到基因组后得到的SAM文件 （SAMtools 包含一些perl 脚本可以将大多数的比对文件转换成SAM格式 ），注意 基因组mapping时一定要用支持剪接的比对软件（splicing-aware aligner）进行比对软件如TopHat.  HTSeq-count 需要用到SAM格式中的 CIGAR 区域的信息。
   a.  想要通过标准输入来传入 基因组mapping得到SAM文件，用 – 替换 <alignment_file> 即可
   b.  如果你是双端测序，必须要对SAM进行排序（单端可不必排序，但这里我也推荐对单端测序结果排序已减少内存消耗并提高软件效率），对read name或 染色体位置 进行排序皆可（这里我推荐按read name排序，因为通过位置排序我遇到过错误）。具体需要通过 -r 参数指定，所以排序请详见参数 -r
2. <gff_file>:  包含单位信息的gff/GTF文件（gff文件格式），大多数情况下就是指注释文件;   由于GTF文件其实就是gff文件格式的变形，在这里同样可以传入GTF格式文件。

模型：

如何判断一个reads属于某个基因， htseq-count 提供了union, intersection_strict, intersection_nonempty 3 种模型，如图（大多数情况下作者推荐用union模型）：

如果这3种模型不能满足您的需求，您还可以通过htseq-count 自定义模型，方法详见  A tour through HTSeq 。

参数：

-f <format>, –format=<format> ：

指定输入文件的格式，<format> 可以是 sam (for text SAM files) 或 bam (for binary BAM files)格式，默认是sam.

-r <order>, –order=<order>

对于双端测序数据，必须要对SAM文件进行排序，对read name或 位置 进行排序皆可，通过 -r 可以指定您的数据是以什么方式排序的： <order> 可以是 name 或 pos ， 默认是name.
如果您的数据没有排序，可以通过samtools sort(推荐)或msort 排序， 命令：

samtools sort -n accepted_hits_unique.bam -o accepted_hits_unique_name_sorted.bam   # -n 按read name 排序 ，如果不指定则按染色体位置排序

如果你指定name（按read name 排序）, htseq-count 期望输入的文件中的read pair是紧邻的2行；而指定pos (按位置排序) ，则不需要这样，进一步说，如果指定pos , 可以不排序，但这样会耗费更多的内存和效率。

-s <yes/no/reverse>, –stranded=<yes/no/reverse>

是否这个数据是来自链特异性建库（默认 yes)
For stranded=no, a read is considered overlapping with a feature regardless of whether it is mapped to the same or the opposite strand as the feature. For stranded=yesand single-end reads, the read has to be mapped to the same strand as the feature. For paired-end reads, the first read has to be on the same strand and the second read on the opposite strand. For stranded=reverse, these rules are reversed.

-a <minaqual>, –a=<minaqual>

指定一个最低 read mapping质量值，低于<minaqual>值会被过滤掉（默认是10，0.5.4版本以前默认值是 0）

-t <feature type>, --type=<feature type>

指定最小计数单位类型（gff文件的第3列中的类型如： exon）, 指定后其他单位类型将被忽略（默认情况下，对于rna-seq分析采用 Ensembl GTF 文件类型时,默认值是：exon）。

-i <id attribute>, --idattr=<id attribute>

GFF文件的一类属性，最终的技术单位，默认情况下，对于rna-seq分析采用 Ensembl GTF 文件类型时,默认值是：gene_id）。

-m <mode>, --mode=<mode>

判断一个reads属于某个基因的模型，用来判断统计reads的时候对一些比较特殊的reads定义是否计入。 <mode> 包括：默认的union和intersection-strict、 intersection-nonempty  （默认：union）。

-o <samout>, --samout=<samout>输出所有alignment的reads到一个叫 <samout>的sam文件中，通过一个可选的sam标签 ‘ XF ’来标注每一行对应的单位和计数，可以不设置。 -q, --quiet

屏蔽程序报告和警告

-h, --help

打印帮助

输出：

文件类似于：

ENSG00000000003.14      21
ENSG00000000005.5       8
ENSG00000000419.12      13
ENSG00000000457.13      65
......
ENSG00000282815.1 25
__no_feature 42987809     #不能对应到任何单位类型的reads数
__ambiguous 183025        #不能判断落在那个单位类型的reads数
__too_low_aQual 0         #低于-a设定的reads mapping质量的reads数
__not_aligned 0           #存在于SAM文件，但没有比对上的reads数
__alignment_not_unique 0   #比对到多个位置的reads数

 

下面是我自己的运行代码以供借鉴：

tophat得到accepted_hits.bam后：
samtools view -h -q 20 accepted_hits.bam |grep -E "^@|NH:i:1$|NH:i:1[^0-9]" > accepted_hits_unique.sam  #筛选unique mapping 
samtools view -bS accepted_hits_unique.sam > accepted_hits_unique.bam  #sam转到bam
samtools sort -n accepted_hits_unique.bam -o accepted_hits_unique_name_sorted.bam  #按read name排序
htseq-count -f bam accepted_hits_unique_name_sorted.bam /leofs/chenf_group/yangl/pub/gencode/annotation/gencode.v23.annotation.gff3 > htseq_count_number.txt   #计数

 

拓展：
软件manual(常见问题可以查看FAQ): http://www-huber.embl.de/users/anders/HTSeq/doc/count.html