计算表达量可以用 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,只有比对多个不同的区间时,才会分别计数。
一、软件下载及安装:
首先是官方网站
https://sourceforge.net/projects/subread/
http://subread.sourceforge.net
二、参数说明:
-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。以后试试。
参考:
//www.greatytc.com/p/9cc4e8657d62
http://www.360doc.com/content/21/0714/12/76149697_986499746.shtml
https://pubmed.ncbi.nlm.nih.gov/24227677/
http://subread.sourceforge.net/featureCounts.html
http://subread.sourceforge.net/RNAseqCaseStudy.html
本文使用 文章同步助手 同步