htseq-count使用说明

Posted on | In | | Visitors |

htseq-count是一款用于reads计数的轻便软件，可以用于多种mapping软件（tophat、HISAT2、BWA等）的输出结果进行计数。

一、htseq-count参数简介

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# 用法概述
usage: htseq-count [options] alignment_file gff_file

# -h参数可显示帮助列表
  -h, --help            show this help message and exit

# -f参数指定输入文件格式类型，默认文件类型为sam
  -f {sam,bam}, --format {sam,bam}

# -r参数指定文件的排序方式，pos:按照染色体位置排序，name:按照read名称进行排序。双端测序数据必选参数，默认值为name。对于单端测序数据，该选项可以忽略
  -r {pos,name}, --order {pos,name}
  
# 当-r参数设定为pos，该选项可以选择最大内存，该参数对单端测序数据无效
  --max-reads-in-buffer MAX_BUFFER_SIZE

# -s参数指定数据建库的链特异性情况，默认值为yes。对于双端测序数据，大多数为非链特异性建库
  -s {yes,no,reverse}, --stranded {yes,no,reverse}
                       
# -a参数指定最低read mapping质量值，低于<minaqual>值会被过滤掉（默认值为10）
  -a MINAQUAL, --minaqual MINAQUAL
                  
# -t参数指定指定最小计数单位类型，（GFF文件中的第三列：如exon），当RNA-seq分析采用 Ensembl GTF 文件类型时，默认值是exon           
  -t FEATURETYPE, --type FEATURETYPE

# -i参数指定最终作为特征id的值，当分析采用Ensembl GTF文件类型是，默认值是gene_id
  -i IDATTR, --idattr IDATTR

# 
  --additional-attr ADDITIONAL_ATTR [ADDITIONAL_ATTR ...]
                        Additional feature attributes (default: none, suitable
                        for Ensembl GTF files: gene_name) 

# -m参数指定判断一个reads属于某个基因的模型，用来判断统计reads的时候对一些比较特殊的reads定义是否计入。<mode> 包括：默认的union和intersection-strict、 intersection-nonempty  （默认：union） 
  -m {union,intersection-strict,intersection-nonempty}, --mode {union,intersection-strict,intersection-nonempty}
                        mode to handle reads overlapping more than one feature
                        (choices: union, intersection-strict, intersection-
                        nonempty; default: union)
  --nonunique {none,all}
                        Whether to score reads that are not uniquely aligned
                        or ambiguously assigned to features
  --secondary-alignments {score,ignore}
                        Whether to score secondary alignments (0x100 flag)
  --supplementary-alignments {score,ignore}
                        Whether to score supplementary alignments (0x800 flag)
  -o SAMOUTS [SAMOUTS ...], --samout SAMOUTS [SAMOUTS ...]
                        write out all SAM alignment records into an output SAM
                        file called SAMOUT, annotating each line with its
                        feature assignment (as an optional field with tag
                        'XF')
  -q, --quiet           suppress progress report

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

对于双端测序数据，必须要对SAM文件进行排序，对read name或 染色体位置 进行排序皆可，但是推荐使用read name进行排序，将会大大节约时间及服务器资源。通过 -r 可以指定数据是以什么方式排序的： 可以是 name 或 pos ， 默认是name.
如果数据没有排序，可以通过samtools sort(推荐)或msort 排序。

二、htseq-count的输入文件

输入为sam或bam格式文件，需要注意的是，如果是paired-end数据必须按照reads名称排序（sort by name），否则运行时间很长，且内存占用高。

1
samtools sort -n file.bam #sort bam by name

三、htseq-count的使用和参数

1
Usage：htseq-count [options] <sam_file> <gff_file>

主要参数说明：

-m 计数模型，统计reads的时候对一些比较特殊的reads定义是否计入。包括：默认的union和intersection-strict、 intersection-nonempty具体说明如图所示。
htseq-count-s reads是否匹配到同一条链上，默认：yes，可以设置no 、 reverse。对于一般的RNA-seq数据，一般选择no
-t feature type 计数单位类型，在gtf或者gff文件中，外显子为最小的定义单位，可选参数有：exon, gene, transcript。默认值：exon
-i 最终的计数单位，一般为基因。 默认为：gene_id 也可以设置转录本，但由于模型问题，计数效果不佳。
-o 输出所有alignment的reads到一个sam文件中。可以不设置。

四、htseq-count的结果文件

结果文件如下图所示：

![image_1bni91970eac128789jmstl4o9.png-56.4kB][1]
[1]: http://static.zybuluo.com/oxygenjing/cneczjpj4ruad3puui5231zp/image_1bni91970eac128789jmstl4o9.png

可以看到是两列的文件，第一列是ensemble id，第二列则是该基因的count值。文件末尾会对重复值或没有匹配上的值进行统计。

多个样本通过htseq-count得到的count文件可以通过脚本进行可并为count值矩阵，进而使用DEseq2等R包进行下游分析。