一个基因在不同的数据库有不同的名字：
1.Entrez gene ID：我们一般说的Gnen ID即Entrez gene ID，是用一串数字表示的（在NCBI里面用）
2.Gene Symbol：可以理解为基因的官方名称，如TP53
3.Ensembl ID：Ensembl ID形式：ENSG00000223972
思路：实现 两者的转换，首先要找对原ID和你要注释的ID的对应关系
怎么找对对应关系呢？
可以用R包，如org.Hs.eg.db
也可以加载相关文件，如gene2ensembl，gene2accession，gene_info
下载地址:[ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/]
理解这个原理后我们再实际操作一下：
方法一：利用encode数据库的gtf注释转换ensembl_id
加载需要ID转换的文件
load('input_exprSet.Rdata')
head(exprSet)
exprSet
library(stringr)
ids=data.frame(ensembl_id=str_split(rownames(exprSet),
'[.]',simplify = T)[,1],
median=apply(exprSet,1,median)
)
head(ids)
head(ids$ensembl_id)
因为表达矩阵的行名点前面的部分就是ensembl_id，所以用str_split(rownames(exprSet),'[.]',simplify = T)[,1]取出前面部分
median=apply(exprSet,1,median)对每一行计算中位数
ids
加载gencode数据库的gtf注释文件
load('human_geneInfo_genecode_v25.rda')#gencode数据库的gtf注释
head(human_geneInfo_genecode_v25)#可以看到有symbol和ensembl的对应关系
human_geneInfo_genecode_v25
s2e=human_geneInfo_genecode_v25[,c(4,6)]#取出symbol和ensembl的对应关系
head(s2e)
s2e
table(ids$ensembl_id %in% s2e$ensembl)
QQ截图20190719155942.jpg
可以看到我们需要的表达矩阵的ensembl_id在gencode数据库的gtf注释能找到25103个，只有31个找不到
利用数据库的gtf注释文件转换ensembl_id
ids=ids[ids$ensembl_id %in% s2e$ensembl,]
#取出在gencode数据库的gtf注释中能找到的ensembl_id
ids$symbol=s2e[match(ids$ensembl_id,s2e$ensembl),2]
#match返回其第二个参数中第一个参数匹配的位置
# 把s2e的ensembl按照ids$ensembl的顺序一个个取出来，从而得到ids$symbol这一列
QQ截图20190719160550.jpg
length(unique(ids$symbol))
head(ids)
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#把ids$symbol按照ids$median排序
QQ截图20190719160807.jpg
可以看到整个ids的行都按照median的大小排列了
ids=ids[!duplicated(ids$symbol),]#取出不重复的ids$symbol
dim(ids)
exprSet= exprSet[rownames(ids),]#取出表达矩阵中ids有的行
微信截图_20190719161238.png
rownames(exprSet)=ids$symbol#把ids$symbol变为exprSet的行名
exprSet[1:4,1:4]
dim(exprSet)
进行转换成功
方法二：org.Hs.eg.db包来注释ensembl_id
加载需要转换的数据，提取出其中的ensembl_id（和前面是一样的）
rm(list = ls()) ## 魔幻操作，一键清空~
load('input_exprSet.Rdata')
library(stringr)
ids=data.frame(ensembl_id=str_split(rownames(exprSet),
'[.]',simplify = T)[,1],
median=apply(exprSet,1,median)
)
head(ids)
head(ids$ensembl_id)
这里用的是org.Hs.eg.db包来注释ensembl_id，可以把org.Hs.eg.db包理解为一个数据存储包，里面含着对应探针ID关系
library(org.Hs.eg.db)
g2s=unique(toTable(org.Hs.egSYMBOL))
head(g2s)
提取出SYMBOL信息
g2e=unique(toTable(org.Hs.egENSEMBL))
head(g2e)
提取出ENSEMBL信息
s2e=merge(g2e,g2s,by='gene_id')
整合两个数据，得到ENSEMBL和SYMBOL的对应关系
table(ids$ensembl_id %in% s2e$ensembl)
微信截图_20190719163248.png
这里可以看到，用这个方法，又一半的基因在注释包中找不到
得到对应关系后进行ensembl_id转换，方法和上面的就是一样的了
ids=ids[ids$ensembl_id %in% s2e$ensembl,]
#取出在gencode数据库的gtf注释中能找到的ensembl_id
ids$symbol=s2e[match(ids$ensembl_id,s2e$ensembl),2]
#match返回其第二个参数中第一个参数匹配的位置
# 把s2e的ensembl按照ids$ensembl的顺序一个个取出来，从而得到ids$symbol这一列
length(unique(ids$symbol))
head(ids)
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#把ids$symbol按照ids$median排序
ids=ids[!duplicated(ids$symbol),]#取出不重复的ids$symbol
dim(ids)
exprSet= exprSet[rownames(ids),]#取出表达矩阵中ids有的行
rownames(exprSet)=ids$symbol#把ids$symbol变为exprSet的行名
exprSet[1:4,1:4]
dim(exprSet)
微信截图_20190719163534.png
可以看到只有14224个基因注释成功
网友评论