空中不产楼阁Linux与生物信息2020生物信息学

关于stringtie定量基因的时候,最后输出很多MSTRG样式

2020-03-20  本文已影响0人  Caspare

    相信大家在用hisat2-stringtie-DESeq2这一套流程做差异表达基因分析的时候的时候,最后会输出很多带MSTRG字样的geneid。

我一开始搜索这个问题,网上有些人给出的答案是这个带MSTRG样id的是新发现的新转录本,但是我做出来的结果几乎一半都是带这个id的基因,我觉得不太可能。于是我在外网搜索到了一篇关于这个问题的开发者的答案,大概的意思就是用stringtie在run脚本的时候因为是多线程的,所以每一个线程分开运行,当接收到一个gene_id的时候会先给他一个MSTRG id,方便之后在合并的时候不会乱,于是下面有人回答可以修改官网给出的prepDE.py的脚本将第26行识别的gene_id改成ref_gene_id,经过我的尝试我发现将RE_GENE_ID=re.compile('gene_id "([^"]+)"')修改为RE_GENE_ID=re.compile('gene_name "([^"]+)"'),之后就仅仅出现了大概一千多个MSTRG。

这比较符合新发现的基因的说法。

我以为这件事告一段落,但是就在我做RNA_seq的时候,我发现了一个问题,在一个基因内我发现了一个snp/indel,但是在gene_count.csv文件中没有搜到这个基因,我觉得不太可能,如果这个基因没有表达又是如何找到这个snp/indel差异的呢,于是我看了Stringtie产生的gtf文件,在里面发现了这个基因,但是他的gene_id是MSTRG id,而后面的gene_name是我要找到基因名

于是我以为是prepDE.py的问题,我继续修改,改成了ref_gene_id,结果还是没有用。于是今天我又继续寻找问题的答案,就在我看到一篇文章解决 Stringtie 基因重复定量的意外收获,让我恍然大悟,大概就是这个merge gtf这一步的时候,生成的,同一个基因没有重叠的转录本会分割成两个基因,所以会赋予MSTRG的id ,后来我重新搜索基于stringtie-DESeq2的分析流程,hisat2+stringtie+deseq2分析RNA-SEQ数据,发现可以不用merge gtf这一步,所以我直接用自己的注释基因组,最后做出来果然没有MSTRG  id了。也就是说如果你对新发现的转录本不感兴趣,就可以不用merge gtf这一步。其实如果你用bollown包来做差异分析,理论上不会有MSTRG的id了,因为他后面会根据GFF文件注释你的MSTRG对应的gene_name。但是我没有做过,所以有人做过可以分享一下。


      因为看到很多人有这个问题,但是没有人分享,不知道我理解的对不对,生信路上感恩前人的经验,分享一下自己的经验。

上一篇下一篇

猜你喜欢

热点阅读