087-问题：转录本ID转换？

Feb 25, 2019 4 min read cnposts

刘小泽写于19.2.25
好久没写咯，最近一直在写毕业论文，就给自己一周的时间写完，然后开始学习新知识😎期待！到时不断更的我又会出现
重回正题~Trinity确实是个非常好的软件，官方文档写的十分详细，而且主流的无参转录组分析都推荐使用这款，它可以预测出更长的亚型、更多的基因和转录本，另外它还有几十个小工具，比如不会做热图，不会做火山图都可以借助它来进行可视化，完全不需要本人需要会R知识，直接一个shell命令搞定。去年2月我刚开始接触生信不久，就做了这个内容，因此记得特别清楚那种感觉——不明不白整出来一个自己看不懂的pdf，导出以后，哇这么高大上！
但随着时间的推移，知识是不断更迭的，现在如果还是用老一套程序自带的小工具，就有点落伍了

最近就遇到这样的一个问题

不想用Trinity自带的作图脚本，因为做出来有点nan kan🤥并且不能按照我们自己的需求去定制，这一点是综合性工具的通病。

例如：自带的火山图【想改颜色？想加基因名？想加辅助线？对不起，Trinity不可以】

Trinity火山图

因此，如果对R语言还算熟悉的话，直接可以从定量结束后，自己去做，无论用limma、deseq、edger，只要你喜欢，都可以实现，并且有太多的可视化包，比如ComplexHeatmap（https://github.com/jokergoo/ComplexHeatmap）做出来的图都很好看。这里为了快速演示学习，我就用了简单的pheatmap功能

Trinity热图

可以看到，我这里截取了前30基因进行不同样本间的热图绘制，结果看右侧全是Trinity，我们自己可能知道是Trinity拼接的转录本（或者也勉强可以叫基因），但是对于别人不清楚Trinity的，就不知道我们做了啥。

别人想看到的，可能仅仅是几十个基因名

例如，在利用DESeq2得到差异基因后，几十个Trinity ID，如何让它们变成基因名，方便下游的功能注释【因为使用clusterProfiler对非模式生物进行注释，需要ID的转换，并且转换列表里肯定没有Trinity ID选项】

肯定不能让我们自己一条一条去找序列然后blast吧~想想都不行！

于是，我们之前对Trinity.fasta做的各种注释派上了用场：nr注释、uniprot注释，其中都有基因信息啊

再说一个经验：不要想着一次就实现批处理，先对一条序列或一个基因做通了，批处理也不再是难事

举个例子：我现在在热图中看到20条基因在样本间差异表达还比较明显，我想想看看它们，先从中抽出来一条让我们写写脚本练练手，就拿TRINITY_DN1458_c0_g1_i1来说吧

我们之前跑完了Uniprot的注释，结果文件是blastx.outfmt6，重点关注第二列

UniProt注释

第二列记录了数据库比对到的Uniprot编号与物种，以第一行为例，Q14203是UniProt编号，可以直接查询；DCTN1_HUMAN是比对到的名称和物种

我们看看TRINITY_DN1458_c0_g1_i1这个转录本对应的信息是什么吧：

grep "TRINITY_DN1458_c0_g1_i1" blastx.outfmt6

# TRINITY_DN1458_c0_g1_i1	sp|P31417|FABP2_MANSE	54.89	133	59	1	63	461	1	132	1e-45	153

然后我们需要导出Uniprot ID号

grep "TRINITY_DN1458_c0_g1_i1" blastx.outfmt6 | cut -f 2 | cut -d "|" -f2
# P31417

这样一条基因的就搞定啦！

看看clusterProfiler支不支持吧~

library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
## [25] "UNIGENE"      "UNIPROT"

看到了"UNIPROT"的选项，于是可以愉快地使用bitr函数去转换啦

只不过需要使用自己研究物种的orgdb

或者，使用Uniprot在线工具也行（https://www.uniprot.org/uploadlists/）

如果有多条转录本呢？
一条的会了，多条的就是一个批处理而已

还是首先从Rstudio中导出那些需要转换的，首先这是一个字符型向量，怎么导出为txt呢？由于好久不接触R，有点生疏了😅，问了花花，2行代码搞定【最后花花还附加了一句：一切能放到R里的问题都难不倒我】

top30 <- as.data.frame(top30_gene)
write.table(top30, file="test.txt", row.names = F, quote = F)

多条trinity ID

然后，去linux中放到一个文本

# 这里截取一部分
cat >test.txt
TRINITY_DN1458_c0_g1
TRINITY_DN22728_c0_g1
TRINITY_DN841_c0_g1
TRINITY_DN4036_c1_g1
TRINITY_DN5293_c0_g1
TRINITY_DN2511_c0_g1
TRINITY_DN5135_c1_g1
TRINITY_DN63171_c0_g1
TRINITY_DN28151_c1_g1
TRINITY_DN2177_c0_g1
TRINITY_DN1250_c11_g1

然后批处理

需要考虑到有的转录本没有对应的注释，此时就需要指明，否则顺序就乱了；另外还有一种情况需要注意，就是一个转录本存在好几个注释，如下图

多个对应

cat top30.txt | while read id;do if grep -q -w $id blastx.outfmt6 ;then (i=${id};grep -w ${i} blastx.outfmt6|cut -f 2 | cut -d "|" -f2);else echo "${id}";fi ;done

cat trans.txt
# 截取一部分
P31417
P31416
P31416
P31416
Q26365
O61387
P0C220
P0C220
Q9CWK8
Q9CWK8
Q2TBW7

当然，这只是最粗略的方法，还有一些问题需要解决，比如多个比对选分值最高；完全相同的比对需要去重复
最后得到了对应的基因名，可以看到，虽然有的转录本可以转换成Uniprot ID，但依然存在对应不到基因的情况，这种在后续注释中可能就要抛弃了

网页版结果

RNA-Seq

087-问题：转录本ID转换？

最近就遇到这样的一个问题

别人想看到的，可能仅仅是几十个基因名

Yunze Liu

Bioinformatics Sharer

Related