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 參數指定的值都是唯一標識符,才能達到預期的效果;
