強強聯(lián)合的基因ID轉(zhuǎn)換

劉小澤寫于2020.5.28
Ensemble轉(zhuǎn)Symbol其實不是這么簡單阅茶，問題百出，需要留意

0 今天遇到一個需求谅海，先來看看

有這樣的四個文件脸哀，每個文件中都有兩列Ensembl ID，從命名可以看到是人類的ID

如何辨別Ensembl ID扭吁？

例如ENSG00000279928.1中撞蜂，ENS是固定字符，表示它是屬于Ensembl數(shù)據(jù)庫的智末。默認物種是人谅摄，如果是小鼠就要用ENSMUS開頭，關(guān)于物種代碼：http://www.ensembl.org/info/genome/stable_ids/index.html

G表示：這個ID指向一個基因系馆；E指向Exon送漠；FM指向 protein family；GT指向gene tree由蘑；P指向protein闽寡；R指向regulary feature；T指向transcript

后面11位數(shù)字部分如00000279928 表示基因真正的編號

小數(shù)點后不一定每個都有尼酿，表示這個ID在數(shù)據(jù)庫中做了幾次變更爷狈，比如.1做了1次變更板熊，在分析時需要去除

至于怎么去掉小數(shù)點后面的信息也很簡單丢胚，之前花花介紹過（去掉ensembl id的version信息）

然后每一行表示：一個基因（第一列）與另一個基因（第二列）編碼的蛋白存在互作關(guān)系

于是可以看到很多這樣的情況：一個基因（例如ENSG00000000005）與其他十幾個基因都有關(guān)系景用。我們需要保證轉(zhuǎn)換后的順序和原來的順序一一對應

1 先處理一個文件的轉(zhuǎn)換

既然有兩列蔚晨，那思路就是：先對第一列的id進行轉(zhuǎn)換，然后按照原來的順序添加到新的一列兼贡，即第三列夕春；再對第二列的id進行轉(zhuǎn)換投慈，然后按照原來的順序添加到新的一列惶我，即第四列

1.1 首先讀取數(shù)據(jù)看看

rm(list=ls())
options(stringsAsFactors = F)
library(tidyverse)
library(org.Hs.eg.db)
library(clusterProfiler)
library(stringr)

df=read.csv('test1.tsv',sep="\t",header = F)
> dim(df)
[1] 64006     2
# 這里的6萬多行不是真實的6萬個基因妈倔，因為存在大量的重復ID情況
> length(unique(df$V1))
[1] 6919 # 第一列實際有6919個基因
> length(unique(df$V2))
[1] 7162

1.2 然后對第一列進行轉(zhuǎn)換

# drop參數(shù)：drop NA or not
conv=bitr(df$V1,fromType = "ENSEMBL",
           toType = "SYMBOL",OrgDb = org.Hs.eg.db, drop = F)
> dim(conv)
[1] 6980    2
# 看到相比原來的6919個基因，少了一些绸贡，就要想到Ensemble ID對應Symbol時盯蝴，存在一對多的關(guān)系，即一個Ensemble 可能對應多個Symbol听怕，檢查一下
conv[which(duplicated(conv$ENSEMBL)),]
# 如果一個ID出現(xiàn)了2次或者多次捧挺，那就會返回這些重復的位置，再取出來看看

1.3 特殊情況一：存在一對多

看到這4個ENSG00000215269尿瞭，實際上總共出現(xiàn)了5次松忍，在6862這個位置是第一次出現(xiàn)，之后的6863-6866這幾個symbol就被認定為了重復

> conv[6862,]
             ENSEMBL SYMBOL
6862 ENSG00000215269  GAGE4

既然存在一對多筷厘，那么就應該保留唯一的一個symbol

一般為了代碼的便利性鸣峭，會選擇第一個，但第一個對不對呢酥艳？這里舉兩個例子：

1.3.1 第一個例子：ENSG00000215269

library("AnnotationDbi")
# 第一個例子：ENSG00000215269
geneSymbols <- mapIds(org.Hs.eg.db, keys='ENSG00000215269', column="SYMBOL", keytype="ENSEMBL", multiVals="CharacterList")

> unlist(geneSymbols)
ENSG00000215269 ENSG00000215269 ENSG00000215269 ENSG00000215269 ENSG00000215269 
        "GAGE4"         "GAGE6"         "GAGE7"       "GAGE12I"       "GAGE12G" 
# 第一個是GAGE4摊溶，最后一個是GAGE12G

但是當把這個Ensemble ID去搜索，會發(fā)現(xiàn)并不是第一個充石，而應該是最后一個GAGE12G：

再次檢索一下第一個GAGE4基因：

對比一下二者：

這時會想莫换，那我不用第一個了，用最后一個唄骤铃。事情依然沒有這么簡單：

1.3.2 第二個例子：ENSG00000063587

可以看到這個ENSG00000063587基因存在兩個symbol：ZNF275和LOC105373378

> conv[389:392,]
            ENSEMBL       SYMBOL
389 ENSG00000063515         GSC2
390 ENSG00000063587       ZNF275
391 ENSG00000063587 LOC105373378
392 ENSG00000063854         HAGH

搜索一下就知道拉岁，這第一個symbol：ZNF275 是可靠的

好，用第一個也不對惰爬，用最后一個也不對喊暖，那應該怎么辦呢？
這時要考慮換一下注釋數(shù)據(jù)庫

ID轉(zhuǎn)換方法不止一種撕瞧，詳見之前我寫的：gene的symbol與entrez ID并不是絕對的一一對應的

尤其是針對Ensemble ID時陵叽，自家的轉(zhuǎn)換能獲得更多的結(jié)果

library("biomaRt")
# 用useMart函數(shù)鏈接到人類的數(shù)據(jù)庫
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
# 除以以外還有很多：使用 listDatasets(ensembl) 查看

attributes <- listAttributes(ensembl)
View(attributes)

value <- df$V1
attr <- c("ensembl_gene_id","hgnc_symbol")
ids <- getBM(attributes = attr,
             filters = "ensembl_gene_id",
             values = value,
             mart = ensembl)

# 如果有匹配不上的，還是用原來的Ensembl ID
ids$hgnc_symbol[ids$hgnc_symbol == ""] <- ids$ensembl_gene_id[ids$hgnc_symbol == ""]
ids

1.4【補充】特殊情況二：存在多對一

這一小部分只是為了文章結(jié)構(gòu)的完整性丛版，這里的數(shù)據(jù)沒有體現(xiàn)
那么只存在一對多嗎巩掺？并不是！還有多個Ensemble 對應一個symbol

看：https://www.researchgate.net/post/How_to_deal_with_multiple_ensemble_IDs_mapping_to_one_gene_symbol_in_a_RNA-Seq_dataset

但其實看一下Ensemble ID在染色體上的位置就能知道：

只有下圖中加粗的那個ID才位于染色體的主體（chr19）页畦，其他的都來自haplotypic regions

名詞解釋：
Haplotyptic regions are defined by the Genome Reference Consortium (GRC) and are visible on all their genome assemblies for human, mouse, and zebrafish. They can also be called ‘alternate locus’, ‘partial chromosomes’, and ‘a(chǎn)lternate alleles’ or ‘a(chǎn)ssembly exceptions’
These regions can have small sequence differences, contain different gene structures or different genes entirely, or contain genes in a different order compared to the reference genome assembly.

2 下面對比bitr和biomart

下面的conv是bitr的結(jié)果胖替，ids2是biomart的結(jié)果，df是原來的ID數(shù)據(jù)框

conv2=conv[!duplicated(conv$ENSEMBL),] # 簡單去重豫缨，只保留第一個
df2=left_join(df,conv2,by=c("V1"="ENSEMBL")) # bitr結(jié)果
df3=left_join(df,ids2,by=c("V1"="ensembl_gene_id")) # biomart結(jié)果
> dim(df);dim(df2);dim(df3)
[1] 64006     2
[1] 64006     3
[1] 64006     3

# 上面的第一個例子：ENSG00000215269
> df2[df2$V1=='ENSG00000215269',]
                   V1              V2 SYMBOL
63685 ENSG00000215269 ENSG00000283632  GAGE4
> df3[df3$V1=='ENSG00000215269',] 
                   V1              V2 hgnc_symbol
63685 ENSG00000215269 ENSG00000283632     GAGE12G

看看df2中一些NA的情況：例如ENSG00000064489在bitr中沒有独令，但在biomart中有

https://asia.ensembl.org/Homo_sapiens/Gene/Summary?g=ENSG00000064489;r=19:19145569-19192158

df2[df2$V1=='ENSG00000064489',]
df3[df3$V1=='ENSG00000064489',]

image

但能說biomart就一定比bitr好嗎？顯然不能這么說州胳，例如ENSG00000036549就在bitr結(jié)果中有记焊，而在biomart結(jié)果沒有

# 再看一個例子：
df2[df2$V1=='ENSG00000036549',]
df3[df3$V1=='ENSG00000036549',]

并且搜索驗證一下，ENSG00000036549 還真的有這個Symbol ID對應

不過這里是GRCh37 (hg19)基因組版本

http://grch37.ensembl.org/Homo_sapiens/Gene/Summary?g=ENSG00000036549;r=1:78028101-78149104

如果是GRCh38（hg38）的話栓撞，結(jié)果就是：

http://asia.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000036549;r=1:77562416-77683419

3 代碼整合版

上面的內(nèi)容作為探索過程遍膜，理解做了什么就好

df=read.csv(files[i],sep="\t",header = F)

## 用biomart轉(zhuǎn)換id
# 用useMart函數(shù)鏈接到人類的數(shù)據(jù)庫
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
# 除以以外還有很多：使用 listDatasets(ensembl) 查看
# attributes <- listAttributes(ensembl)
# View(attributes)

## 對第一列操作
value <- df$V1
attr <- c("ensembl_gene_id","hgnc_symbol")
conv <- getBM(attributes = attr,
              filters = "ensembl_gene_id",
              values = value,
              mart = ensembl)

# 如果有匹配不上的，還是用原來的Ensembl ID
if(sum(conv$hgnc_symbol == "") != 0){
    conv$hgnc_symbol[conv$hgnc_symbol == ""] <- conv$ensembl_gene_id[conv$hgnc_symbol == ""]
}

# 存在Ensemble與Symbol一對多的情況瓤湘，去重
conv2=conv[!duplicated(conv$ensembl_gene_id),]
# dim(conv);sum(duplicated(conv$ensembl_gene_id));dim(conv2)
# 沒有對應存為NA
df2=left_join(df,conv2,by=c("V1"="ensembl_gene_id"))

## 對第二列操作
value <- df$V2
attr <- c("ensembl_gene_id","hgnc_symbol")
conv3 <- getBM(attributes = attr,
               filters = "ensembl_gene_id",
               values = value,
               mart = ensembl)
# 如果有匹配不上的瓢颅，還是用原來的Ensembl ID
if(sum(conv3$hgnc_symbol == "") != 0){
    conv3$hgnc_symbol[conv3$hgnc_symbol == ""] <- conv3$ensembl_gene_id[conv3$hgnc_symbol == ""]
}

conv4=conv3[!duplicated(conv3$ensembl_gene_id),]
# dim(conv3);sum(duplicated(conv3$ensembl_gene_id));dim(conv4)

df3=left_join(df2,conv4,by=c("V2"="ensembl_gene_id"))
df3[,3][which(is.na(df3[,3]))]=df3[,1][which(is.na(df3[,3]))]
df3[,4][which(is.na(df3[,4]))]=df3[,2][which(is.na(df3[,4]))]

4 為何不強強聯(lián)合呢？

既然看到有少部分在bitr有注釋弛说，但biomart沒有
那么挽懦，在biomart的基礎(chǔ)上加上bitr的補充結(jié)果，豈不是更好木人？

# 就是看看hgnc_symbol這一列中顯示ENSG0*****的基因在bitr中能不能轉(zhuǎn)換成功信柿，能的話就取代原來的結(jié)果
## 整合bitr
# 先得到org.Hs.eg.db中的所有Ensemble ID
org_ensembl=toTable(org.Hs.egENSEMBL)[,2]

# 對第一列操作
for(i in 1:length(df3$hgnc_symbol.x)){
        name=df3$hgnc_symbol.x[i]
        # 判斷這里的name是不是在org.Hs.eg.db中(org_ensembl)
        if(str_detect(name,"ENSG0") & name %in% org_ensembl){
            df3$hgnc_symbol.x[i]=bitr(name,fromType = "ENSEMBL",
                                    toType = "SYMBOL",
                                    OrgDb = org.Hs.eg.db, drop = F)[,2]
    }
}

# 對第二列操作
for(i in 1:length(df3$hgnc_symbol.y)){
        # i=1
        name2=df3$hgnc_symbol.y[i]
        if(str_detect(name,"ENSG0") & name2 %in% org_ensembl ){
            df3$hgnc_symbol.y[i]=bitr(name,fromType = "ENSEMBL",
                                      toType = "SYMBOL",
                                      OrgDb = org.Hs.eg.db, drop = F)[,2]
    }
}

5 小結(jié)

當轉(zhuǎn)換Ensemble時冀偶，最好還是用自家的biomart，可以注釋到更多的Ensemble ID渔嚷，準確性也會更好一些进鸠，但速度可能會稍微慢一點點；

如果基因ID數(shù)量不多形病，只是想快速看看轉(zhuǎn)換后的ID客年，用bitr 也是足夠的

最后，用了一個上午重新探索了一下ID轉(zhuǎn)換漠吻，最后將bitr和biomart合二為一量瓜，增加了ID轉(zhuǎn)換的準確性和豐富性，并且可以批處理多個文件途乃，還增加了一些小修飾（比如每個文件轉(zhuǎn)換的時間計算绍傲、轉(zhuǎn)換結(jié)果輸出前進一步檢查，保證轉(zhuǎn)換前后的結(jié)果每行的Ensemble ID一致）

歡迎關(guān)注我們的公眾號～_～　　
我們是兩個農(nóng)轉(zhuǎn)生信的小碩欺劳，打造生信星球唧取，想讓它成為一個不拽術(shù)語、通俗易懂的生信知識平臺划提。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!