htseq-count中怎么实现定量分析操作,很多新手对此不是很清楚,为了帮助大家解决这个难题,下面小编将为大家详细讲解,有这方面需求的人可以来学习下,希望你能有所收获。
创新互联公司始终坚持【策划先行,效果至上】的经营理念,通过多达十余年累计超上千家客户的网站建设总结了一套系统有效的全网整合营销推广解决方案,现已广泛运用于各行各业的客户,其中包括:成都三维植被网等企业,备受客户称扬。
和featurecounts一样,htseq-count也是一款进行raw count定量的软件。该软件采用python语言进行开发,集成在HTseq这个包中。
对于python的包,通过pip
可以方便的进行安装,代码如下
pip install HTSeq
HTSeq提供了许多处理NGS数据的功能,htseq-count只是其中进行定量分析的一个模块。
htseq-count的设计思想和featurecounts非常类似,也包含了feature和meta-features两个概念。对于转录组数据而言,feature指的是exon
, 而meta-feature
可以是gene
, 也可以是transcript
。
进行定量分析需要以下两个文件
比对的BAM/SAM文件
基因组的GTF文件
对于双端数据,要求输入sort之后的BAM文件。
由于序列读长的限制和基因组的同源性,一条reads可能比对到多个基因上,而且基因之间也存在overlap, 在对这些特殊情况进行处理时,htseq-count内置了以下3种模式
union
intersection-strict
intersection-nonempty
通过--mode
参数指定某种模式,默认值为union
。这3种模式在判断一条reads是否属于某个feature时,有不同的判别标准,示意图如下
在BAM文件,包含了比对上的reads和没有比对上的reads, 只有比对上的reads 会用来计数,htseq-count默认会根据mapping的质量值对BAM文件进行过滤,默认值为10
, 意味着只有mapping quality > 10的reads才会用来计数,当然可以通过-a
参数来修改这个阈值。
能够明确reads属于一个featurer时,比如示意图种的第一种情况,reads完全是gene_A的一个片段,该feature的计数就加1;能明确reads不属于一个feature时,称之为no_feature, 比如示意图种的第二种情况,reads的一部分比对上了gene_A, 在intersection_strict模式下,判定该reads不属于gene_A, 就是no_feature。
当不明确一条reads是否属于某个feature时,通常是由于reads在两个gene的overlap区,比如示意图中的第六和第七种情况,这样的reads被标记为ambigous。
当一条reads比对上了两个feature时,会被标记为alignment_not_unique。
在统计属于某个基因的reads数时,需要重点关注对 ambiguous 和 alignment_not_unique 的reads的处理, 通过--nonunique
参数来指定,取值有以下两种
none
all
默认值为none
时,这两种reads被忽略掉,不进行任何的计数;取值为all
时,对应的所有feature的计数都会加1。
除了--mode
和--nonunique
两个参数外,还需要关注--stranded
参数,这个参数指定文库的类型,默认值为yes
, 代表文库为链特异性文库,no
代表为非链特异性文库。对于非链特异性文库文库,在判断一条reads是否属于一个基因时,只需要关注比对位置,而链特异性文库还需要关注比对的正负链和基因的正负链是否一致,只有一致时,才会计数。
理解了以上3个参数,就能够正确的使用htseq-count了。对于非链特异性的数据,常规用法如下
htseq-count \ -f bam \ -r name \ -s no \ -a 10 \ -t exon \ -i gene_id \ -m union \ --nonunique=none \ -o htseq.count \ align.sorted.bam \ hg19.gtf
在运行速度上,featurecounts比htseq-count快很多倍,而且feature-count不仅支持基因/转录本的定量,也支持exon等单个feature的定量。所以更加推荐使用featurecounts来定量。
看完上述内容是否对您有帮助呢?如果还想对相关知识有进一步的了解或阅读更多相关文章,请关注创新互联行业资讯频道,感谢您对创新互联的支持。