R生信分析流程R语言与科研

如何让基因名称在多个数据库间随意转换?

2019-07-09  本文已影响46人  9d760c7ce737

这个帖子会成为一个经典。
在GEO芯片的分析中,探针ID转换是个限速环节。这个限速环节迫使我在早期分析的时候,常常在多个在线网站中跳来跳去,我觉得不够优雅,在学会了R语言之后,我就一直尝试去解决这个问题。
简单说来,GEO芯片的探针ID转换,包括三个层面:

第一,R包注释

如果有平台对应的R包,我们就直接下载对应的R包去转换
平台和R包的对应关系我们对应了一个platformMap文件,在果子学生信微信公众号回复“果子学生信”即可获取

第二,平台获取

如果平台没有对应的R包,我们可以下载平台的注释文件,自己提取。
所有的探针ID转换,我们至少要获取两列数据,第一列是已有的探针,第二列对应的基因ID。
有以下的帖子可供参考
skr!GEO芯片数据的探针ID转换
有些GEO平台的探针转换比较麻烦
正则表达式是我们认识世界的哲学
学习正则表达式-stringr这个包的使用-极简入门
GEO芯片中的NM_,NR_开头的识别号如何转换成基因名称?

第三,序列比对

非编码GEO芯片的探针ID转换常常平台信息给出的是序列。
这个稍微有点困难,我们也写了教程,甚至提供了常见平台转换好的文件。在果子学生信微信公众号回复“果子非编码”即可自助获取。
GEO芯片分析的倒数第2个关卡被没有了

这三个方法可以搞定几乎所有GEO芯片的探针ID转换,接下来的问题就是,如果和获得的基因名称变换为各种数据库的名称,比如ensemble ID, entrez ID , 你看我连这些名称都写不清楚。

当然有很多在线网站可以做这个事情,R语言里面也可以完成,今天介绍的是biomaRt,介绍的风格跟学习正则表达式-stringr这个包的使用-极简入门这个帖子一样,用最小的精力,学习最重要的功能。

使用biomaRt来进行ID转换

这个包我们之前介绍过,当时是处理NM_,NR_开头的基因名称
GEO芯片中的NM_,NR_开头的识别号如何转换成基因名称?
今天我们重新讲一下,作为课程的补充。

1.安装及加载R包

options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
if(!require("biomaRt")) BiocManager::install("biomaRt",update = F,ask = F)
library(biomaRt)

biomaRt包ID转换三部曲

我们现在假设一个需求,把三个基因名称转换为entrez ID

第一,要有自己的名称

your_gene_ids <- c("TP53","SMAD4","BRCA1")

第二,要选择所用的数据库

ensembl=useMart(biomart="ensembl",dataset="hsapiens_gene_ensembl")

第三,根据目的直接转换

annotations1 <- getBM(values = your_gene_ids,
                     filters = 'hgnc_symbol', 
                     mart = ensembl,
                     attributes = c('hgnc_symbol', 'entrezgene_id'))

最终得到这样的结果


这里面的核心是第三步getBM转换函数,里面包括四个参数,我来解释一下。

扩展需求来加深理解

1. 如何同时返回多个数据库的结果?

这个很简单,只要在getBM函数中,修改attributes的参数内容即可。难点是,不知道如何表示不同数据库的基因名称。
这个可以通过以下代码实现:

listAttributes <- listAttributes(mart = ensembl)

第一列表示可以转换的类型,比如ensemble数据库的名称用ensembl_gene_id来表示,我们就可以修改代码如下
annotations2 <- getBM(values = your_gene_ids,
                      filters = 'hgnc_symbol', 
                      mart = ensembl,
                      attributes = c('hgnc_symbol',
                                     'entrezgene_id',
                                     "ensembl_gene_id"))

如果是假如转录本的名称,结果会自动扩展行数,因为一个基因会对应多个转录本。

annotations3 <- getBM(values = your_gene_ids,
                     filters = 'hgnc_symbol', 
                     mart = ensembl,
                     attributes = c('hgnc_symbol',
                                    'entrezgene_id',
                                    "ensembl_gene_id",
                                    "ensembl_transcript_id",
                                    "refseq_mrna"))

2. 如果我不是人,怎么办?

这个R包支持多物种的注释,只要在第二步中在useMart里面限定dataset参数即可,首先我们需要知道有哪些物种的数据库可用。

listDatasets= listDatasets(mart = useMart('ensembl'))

给出了190个动物的注释数据,图中我能认识的是熊猫和鸭子。
可以写一个简单的表达式来返回数据库名称,比如人类的

index <- grep("human",listDatasets$description,ignore.case = T)
listDatasets$dataset[index]

"hsapiens_gene_ensembl"

尝试一下兔子的

index <- grep("rabbit",listDatasets$description,ignore.case = T)
listDatasets$dataset[index]

"ocuniculus_gene_ensembl"

现在只要修改useMart函数里面的dataset参数即可。

3.混杂基因名称怎么办?

读入一个数据

dd <- data.table::fread("Diffgenes_JDS.csv")

这个数据是兔子的转录组测序数据,里面可以看到基因名称是混杂的,如何转换为entrez ID呢?


我实际上没有一步实现的方法,要是我就分别转换,先把ensemble的ID抽取出来转换,再把剩下的gene symbol转换一下
先把ensemble ID转换一下,用grep来提取
以下全是三部曲。

library("biomaRt")
index <- grep("ENSOCUG",dd$GeneName)
your_gene_ids <- dd$GeneName[index]
ensembl = useMart("ensembl",dataset="ocuniculus_gene_ensembl")
annotations4 <- getBM(values = your_gene_ids, 
                     mart = ensembl,
                     filters = 'ensembl_gene_id', 
                     attributes=c('ensembl_gene_id', 'entrezgene_id'))

再来转换剩下的gene symbol,他有专业的名称,叫做hgnc_symbol.

library("biomaRt")
index <- grep("ENSOCUG",dd$GeneName)
your_gene_ids <- dd$GeneName[-index]
ensembl = useMart("ensembl",dataset="ocuniculus_gene_ensembl")
annotations5 <- getBM(values = your_gene_ids, 
                      mart = ensembl,
                      filters = 'hgnc_symbol', 
                      attributes=c('hgnc_symbol', 'entrezgene_id'))

最后用rbind函数合并以下就可以了。

好了,我是果子,明天见。

上一篇下一篇

猜你喜欢

热点阅读