featureCounts 软件说明

featuresCounts 软件用于定量，不仅可以支持gene的定量，也支持exon, gene bodies, genomic bins, chromsomal locations的定量；

官网 ： http://bioinf.wehi.edu.au/featureCounts/

只需要输入reads的比对情况，就是BAM 文件，再输入一个你感兴趣的区间的注释（通常是基因或者转录本的注释gtf 文件，就可以了），所以不论是DNA seq 还是RNA seq, 这个软件都是可以定量的。

featureCounts 集成在subreads 软件中，类似 word 和 office 的关系，subreads 这个软件也有对应的 R包（Rsubreads）.

featureCounts 需要两个输入文件：

1）reads的比对情况，这种信息通常都用BAM/ SAM文件来存储

2）区间注释文件，支持两种格式

最常见的gtf 格式

simplified annotation format(SAF) 格式， 示例如下

GeneID	Chr	Start	End	Strand
497097	chr1	3204563	3207049	-
497097	chr1	3411783	3411982	-
497097	chr1	3660633	3661579	-

在featureCounts 软件中，有两个核心概念：

1） feature , 类似 exon 这种

2) metafeature， 可以看做是一组 feature, 比如属于同一个gene 的外显子的组合

在定量的时候，支持对单个feature 定量（对外显子定量），也支持对meta-feature 进行定量（对基因进行定量）， meta-feature 的定量是属于同一meta-features 下的所有features 的总和；

当reads 比对到2个或者以上的features 时，默认情况下，featureCounts在统计时会忽略到这部分reads, 如果你想要统计上这部分reads, 可以添加-O 参数，此时一条reads 比对到多个feature， 每个feature 定量时，都会加1，对于meta-features 来说，如果比对到多个features 属于同一个 meta-features(比如一条reads比对到了exon, 但这些exon 属于同一个gene), 则对于这个gene 而言，只会计数1次；

总之，不管对于feature 还是meta-feature, 只有比对多个不同的区间时，才会分别计数；

定量：

features 支持对单个样本定量，还支持对多个样本进行归一化

单个样本定量的用法示例

featureCounts -T 5 -t exon -g gene_id -a annotation.gtf -o counts.txt mapping.sam

多个样本归一化的用法示例

featureCounts -t exon -g gene_id -a annotation.gtf -o counts.txt library1.bam library2.bam library3.bam

下面对几个常用的选项详细解释一下：

-a  :  指定的区间注释文件，默认是gtf格式

-T  :  线程数，默认是1

-t   :  想要统计的feature 的名称， 取值范围是gtf 文件中的第3列的值，默认是exon

-g  :  想要统计的meta-feature 的名称，取值范围参考gtf 第9列注释信息，gtf 的第9列为 key=value 的格式， -g 参数可能的取值就是所有的key, 默认值是gene_id

其他的一些参数可以根据自己的目的，实际做调整。

输出结果的解读：

featuresCount 会输出两个文件，如果-o 指定的是gene, 则会产生gene 和 gene.summary 两个文件

gene 文件的部分内容如下

# Program:featureCounts v1.6.0; Command:"./featureCounts" "-T" "20" "-t" "exon" "-g" "gene_id" "-a" "hg19.gtf" "-o" "gene" "accepted_hits.bam" 						
Geneid	Chr	Start	End	Strand	Length	accepted_hits.bam
LOC102725121	chr1;chr1;chr1;chr15;chr15;chr15	11869;12613;13221;102516808;102518449;102518943	12227;12721;14362;102517949;102518557;102519301	+;+;+;-;-;-	3220	0
DDX11L1	chr1;chr1;chr1	11874;12613;13221	12227;12721;14409	+;+;+	1652	0
WASH7P	chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1	14362;14970;15796;16607;16858;17233;17606;17915;18268;24738;29321	14829;15038;15947;16765;17055;17368;17742;18061;18366;24891;29370	-;-;-;-;-;-;-;-;-;-;-	1769	88

# 号开头的注释行，记录了产生这个结果文件所用的命令，（感觉这个思路特别好，在输出的文件中记录当时的命令，便于核对）

Geneid 开头的行是表头，Geneid 代表统计的meta-features 的名称，Chr ， Start, End 染色体上的位置，Strand 正负链，Length 该区间的长度，最后一列的表头是你的输入文件的名称，代表的是这个meta-feature 的count 值，即表达量

接下来看下正文部分，以第一行为例，在gtf 文件中共gene_id 为 LOC102725121 的行如下

chr1	refGene	transcript	11869	14362	.	+	.	gene_id "LOC102725121"; transcript_id "NR_148357";  gene_name "LOC102725121";
chr1	refGene	exon	11869	12227	.	+	.	gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "1"; exon_id "NR_148357.1"; gene_name "LOC102725121";
chr1	refGene	exon	12613	12721	.	+	.	gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "2"; exon_id "NR_148357.2"; gene_name "LOC102725121";
chr1	refGene	exon	13221	14362	.	+	.	gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "3"; exon_id "NR_148357.3"; gene_name "LOC102725121";
chr15	refGene	transcript	102516808	102519301	.	-	.	gene_id "LOC102725121"; transcript_id "NR_148357_2";  gene_name "LOC102725121";
chr15	refGene	exon	102516808	102517949	.	-	.	gene_id "LOC102725121"; transcript_id "NR_148357_2"; exon_number "1"; exon_id "NR_148357_2.1"; gene_name "LOC102725121";
chr15	refGene	exon	102518449	102518557	.	-	.	gene_id "LOC102725121"; transcript_id "NR_148357_2"; exon_number "2"; exon_id "NR_148357_2.2"; gene_name "LOC102725121";
chr15	refGene	exon	102518943	102519301	.	-	.	gene_id "LOC102725121"; transcript_id "NR_148357_2"; exon_number "3"; exon_id "NR_148357_2.3"; gene_name "LOC102725121";

对于 LOC102725121 这个meta-features 而言，在gtf 文件中有6个exon的记录，就是说有6个features , 所以可以看到对应的Chr, Start, End, Strand 这些列都有；分号分隔的6个值，Length 则是这6个exon 区间的的长度的总和，最后一列就是LOC102725121的表达量

这个结果文件有1个问题，就是同一个gene_id 会有多个染色体编号，这是因为gtf 文件中的gene_id 不是唯一标识符导致的，这样和我们想要的定量结果是不一样的，所以在实际分析中，应该挑选gtf 文件中的唯一标识符；

总结： 

这个软件最大的特点就是运以行的非常快，几分钟就可以运行完1个人类基因组样本的定量；但是准备gtf 文件时，要确保-g 参数指定的值都是唯一标识符，才能达到预期的效果；