使用htseq-count进行定量分析

使用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种模式

unionintersection-strictintersection-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​​参数来指定，取值有以下两种

noneall

默认值为​​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来定量。

·end·

—如果喜欢，快分享给你的朋友们吧—

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。