featurecounts_featurecounts定量

本文带来了【featurecounts】内容供参考阅读，并对相关内容featurecounts定量进行了分析，下面就跟随币王网小编一起了解featurecountsfeaturecounts定量。

featureCounts 计算基因 reads 数目

建议先阅读以前文章 《StringTie DESeq2 进行 RNA-seq 差异基因分析》

再阅读本文。在前文中使用 StringTie 取得基因 reads 数目，如果你无计划做其他的分析，仅仅进行表达谱分析，可以使用 featureCounts 更直接取得 read counts 矩阵。

featureCounts 非常的灵活，默认情况下是以外显子(exon)为单位 feature , 基因(gene)由多个外显子构成是 meta-feature。 featureCounts 会自动汇总外显子计数给基因。这个默认可以根据自己需要进行修改，分别是 -t 和 -g 参数。

featureCounts默认不计算比对到多个 feature 上的 reads，如果是双端测序那就是 fragment(插入片段)。对于 meta-feature 如果 reads 比对到多个 feature 那只会被计算一次，不会重复计算。

featureCounts 使用很简单自己看看文档即可，可以一次给多个样本给出合并的结果，不用自己再去合并了，只要手动提取一下数据就可以用 DESeq2 进行下游分析了。

featureCounts统计counts后的cpm和tpm计算

该脚本适用于featureCounts数counts后的cpm和tpm计算

setwd("~/Desktop")#设置工作目录

countdata-read.table("counts.txt",skip = 1,sep="\t",header = T,row.names = 1)

metadata - countdata[,1:5]#提取基因信息count数据前的几列

countdata - countdata[,6:ncol(countdata)]#提取counts数，counts数据主题部分

prefix-"couts"#设置输出文件前缀名

cpm - t(t(countdata)/colSums(countdata) * 1000000)#参考cpm定义

avg_cpm - data.frame(avg_cpm=rowMeans(cpm))

#-----TPM Calculation------

kb - metadata$Length / 1000

rpk - countdata / kb

tpm - t(t(rpk)/colSums(rpk) * 1000000)

avg_tpm - data.frame(avg_tpm=rowMeans(tpm))

write.csv(avg_tpm,paste0(prefix,"_avg_tpm.csv"))

write.csv(avg_cpm,paste0(prefix,"_avg_cpm.csv"))

write.csv(tpm,paste0(prefix,"_tpm.csv"))

write.csv(cpm,paste0(prefix,"_cpm.csv"))

转录组定量工具-featureCounts安装及使用

        计算表达量可以用 StringTie、Htseq-count或featureCount ，第一次做转录组分析时，参照了一篇Cell的子刊文章的分析方法，里面使用的STAR featureCount，就直接用了这个软件，也就没再使用别的，回头看第一次使用时，发现好多细节没有注意到，温故而知新。featureCount是subread软件包里的一个命令，所以安装subread即可。而subread又有命令行版和R版，有服务器，自然选择命令行版了。

featureCounts ，有两个核心概念：

       Feature: 指的是基因组区间的最小单位，比如exon;

       Metafeature: 可以看做是许多的feature构成的区间，比如属于同一个gene的外显子的组合。

       在定量的时候，支持对单个feature 定量(对外显子定量), 也支持对meta-feature进行定量(对基因进行定量)。当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, 只有比对多个不同的区间时，才会分别计数。

  首先是官方网站

       

       

-a 指定注释文件

-o 指定结果输出目录及文件名

-p 能用在paired-end的情况中，会统计fragment而不统计read

-t 指定feature的类型，默认是exon，当然gtf里面还有gene、CDS或者直接以feature命名的分类方式。

其它参数：

 -f 参数   该参数设置后统计的是 feature 层面（默认是 exon ）的参数，如果不设置则是直接统计 meta-feature 参数（即一个 gene 中的多个 exon ）

这时按exon分类进行统计，但是由于没有设置-f，在同一个gene内的exon会被统计成一个meta-feature，但是每个exon仍然会被显示出来，遇到一个gene有多个exon的时候看着就很乱。

第二种： 然后我加上-f，这样设置-t exon -f , 看一下结果：

        我现在还不确定-f参数及-t参数对后面差异表达会不会有影响，初步判断不会的，但我注意到，-t gene -f设置后，count计数基于gene 层面，就不会出现相同基因的不同外显子count值，也就是第一列不会出现重复，并且可以直接得到基因信息，避免了注释、删除重复这个过程，我们做转录组测序，不就是想看基因水平的变化吗，我觉得这是很好的一个参数设置，不知道为什么网上一堆的帖子都没有这样设置，官网上示例也只是-t exon。希望未来有人和我讨论一下这个问题。

最终：我基于自己的理解，加上-t gene -f参数了。

1 、运行过程情况：

        Successfully assignedalignments: 14212190 (32.7%), 说明只有32.7的paired reads 定量到了基因上，如果想知道那些没有分配上的reads是出于什么原因，则可看下图，输出中的summary文件。

     Unmapped: 没有比对上；     

     MultiMapping:多个序列比对在有限的序列区域上，即参考组上有多个匹配点； 

     NoFeatures: 其比对与任何基因都不重叠； 

     ambiguous: 其比对与多个基因重叠。

2. 合并不同样本的 count 文件：

        join count1.txt count2.txt count_12.txt

        或者先提取出来每个样本的第一列和第七列信息，再通过join合并

        cut -f 1,7 count1.txt | grep -v ‘^#’ count1_cut.txt

        这样就能得到所有的样本的Count矩阵了。

总结：使用这个工具时要根据不同的项目，不同的目的，参数也要进行适当的调整，尤其是模式生物和非模式生物研究时，一定要想想参数设置合适不合适，我不认为写好了一个流程，就可以用来做所有课题的转录组分析了。这也是自己会和交给公司来做最大的好处了，自己的课题，只有自己才能对数据结果负责。

附：

STAR有一个参数-quantMode，可以指定--quantMode GeneCounts输出STAR计算出的reads计数结果，如果是比对完之后未做转录本拼装，直接对已知基因（构建基因组索引时GTF中囊括的基因）进行定量时，完全不需要再次用featureCounts或HTSeq再计算reads count。以后试试。

参考：

有关featurecountsfeaturecounts定量分享到这里，想要阅读更多相关内容请关注币王网。

本文 原创，转载保留链接！网址：https://licai.bangqike.com/lzs/154485.html