一個基因在不同的數(shù)據(jù)庫有不同的名字：
1.Entrez gene ID：我們一般說的Gnen ID即Entrez gene ID，是用一串?dāng)?shù)字表示的（在NCBI里面用）
2.Gene Symbol：可以理解為基因的官方名稱，如TP53
3.Ensembl ID：Ensembl ID形式：ENSG00000223972
思路：實現(xiàn) 兩者的轉(zhuǎn)換，首先要找對原ID和你要注釋的ID的對應(yīng)關(guān)系
怎么找對對應(yīng)關(guān)系呢？
可以用R包，如org.Hs.eg.db
也可以加載相關(guān)文件，如gene2ensembl，gene2accession，gene_info
下載地址:[ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/]
理解這個原理后我們再實際操作一下：
方法一：利用encode數(shù)據(jù)庫的gtf注釋轉(zhuǎn)換ensembl_id
加載需要ID轉(zhuǎn)換的文件
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)
因為表達(dá)矩陣的行名點(diǎn)前面的部分就是ensembl_id，所以用str_split(rownames(exprSet),'[.]',simplify = T)[,1]取出前面部分
median=apply(exprSet,1,median)對每一行計算中位數(shù)
ids
加載gencode數(shù)據(jù)庫的gtf注釋文件
load('human_geneInfo_genecode_v25.rda')#gencode數(shù)據(jù)庫的gtf注釋
head(human_geneInfo_genecode_v25)#可以看到有symbol和ensembl的對應(yīng)關(guān)系
human_geneInfo_genecode_v25
s2e=human_geneInfo_genecode_v25[,c(4,6)]#取出symbol和ensembl的對應(yīng)關(guān)系
head(s2e)
s2e
table(ids$ensembl_id %in% s2e$ensembl)
QQ截圖20190719155942.jpg
可以看到我們需要的表達(dá)矩陣的ensembl_id在gencode數(shù)據(jù)庫的gtf注釋能找到25103個，只有31個找不到
利用數(shù)據(jù)庫的gtf注釋文件轉(zhuǎn)換ensembl_id
ids=ids[ids$ensembl_id %in% s2e$ensembl,]
#取出在gencode數(shù)據(jù)庫的gtf注釋中能找到的ensembl_id
ids$symbol=s2e[match(ids$ensembl_id,s2e$ensembl),2]
#match返回其第二個參數(shù)中第一個參數(shù)匹配的位置
# 把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),]#取出不重復(fù)的ids$symbol
dim(ids)
exprSet= exprSet[rownames(ids),]#取出表達(dá)矩陣中ids有的行
微信截圖_20190719161238.png
rownames(exprSet)=ids$symbol#把ids$symbol變?yōu)閑xprSet的行名
exprSet[1:4,1:4]
dim(exprSet)
進(jìn)行轉(zhuǎn)換成功
方法二：org.Hs.eg.db包來注釋ensembl_id
加載需要轉(zhuǎn)換的數(shù)據(jù)，提取出其中的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)
ids
這里用的是org.Hs.eg.db包來注釋ensembl_id，可以把org.Hs.eg.db包理解為一個數(shù)據(jù)存儲包，里面含著對應(yīng)探針I(yè)D關(guān)系
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')
整合兩個數(shù)據(jù)，得到ENSEMBL和SYMBOL的對應(yīng)關(guān)系
table(ids$ensembl_id %in% s2e$ensembl)
微信截圖_20190719163248.png
這里可以看到，用這個方法，又一半的基因在注釋包中找不到
得到對應(yīng)關(guān)系后進(jìn)行ensembl_id轉(zhuǎn)換，方法和上面的就是一樣的了
ids=ids[ids$ensembl_id %in% s2e$ensembl,]
#取出在gencode數(shù)據(jù)庫的gtf注釋中能找到的ensembl_id
ids$symbol=s2e[match(ids$ensembl_id,s2e$ensembl),2]
#match返回其第二個參數(shù)中第一個參數(shù)匹配的位置
# 把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),]#取出不重復(fù)的ids$symbol
dim(ids)
exprSet= exprSet[rownames(ids),]#取出表達(dá)矩陣中ids有的行
rownames(exprSet)=ids$symbol#把ids$symbol變?yōu)閑xprSet的行名
exprSet[1:4,1:4]
dim(exprSet)
微信截圖_20190719163534.png
可以看到只有14224個基因注釋成功
