强强联合的基因ID转换

⼤神⼀句话，菜鸟跑半年。我不是⼤神，但我可以缩短你⾛弯路的半年~ 就像歌⼉唱的那样，如果你不知道该往哪⼉⾛，就留在这学点⽣信好不好~ 这⾥有⾖⾖和花花的学习历程，从新⼿到进阶，⽣信路上有你有我！

⾖⾖写于2020.5.28

Ensemble转Symbol其实不是这么简单，问题百出，需要留意

感谢花花提供的tidyverse⽀持，说实话我不太会⽤，但不怕，我有花花啊！

0 今天遇到⼀个需求，先来看看

有这样的四个⽂件，每个⽂件中都有两列Ensembl ID，从命名可以看到是⼈类的ID

如何辨别Ensembl ID？

G表⽰：这个ID指向⼀个基因；E指向Exon；FM指向 protein family；GT指向gene tree；P指向protein；R指向regulary feature；T指向transcript

后⾯11位数字部分如00000279928 表⽰基因真正的编号

⼩数点后不⼀定每个都有，表⽰这个ID在数据库中做了⼏次变更，⽐如.1做了1次变更，在分析时需要去除⾄于怎么去掉⼩数点后⾯的信息也很简单，之前花花介绍过（去掉ensembl id的version信息）然后每⼀⾏表⽰：⼀个基因（第⼀列）与另⼀个基因（第⼆列）编码的蛋⽩存在互作关系

于是可以看到很多这样的情况：⼀个基因（例如ENSG00000000005）与其他⼗⼏个基因都有关系。我们需要保证转换后的顺序和原来的顺序⼀⼀对应

1 先处理⼀个⽂件的转换

既然有两列，那思路就是：先对第⼀列的id进⾏转换，然后按照原来的顺序添加到新的⼀列，即第三列；再对第⼆列的id进⾏转换，然后按照原来的顺序添加到新的⼀列，即第四列1.1 ⾸先读取数据看看

rm(list=ls())

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

df=read.csv('test1.tsv',sep='\',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 然后对第⼀列进⾏转换

# drop参数：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时，存在⼀对多的关系，即⼀个Ensemble 可能对应多个Symbol，检查⼀下conv[which(duplicated(conv$ENSEMBL)),]

# 如果⼀个ID出现了2次或者多次，那就会返回这些重复的位置，再取出来看看

1.3 特殊情况⼀：存在⼀对多

看到这4个ENSG00000215269，实际上总共出现了5次，在6862这个位置是第⼀次出现，之后的6863-6866这⼏个symbol就被认定为了重复

> conv[6862,]

ENSEMBL SYMBOL6862 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去搜索，会发现并不是第⼀个，⽽应该是最后⼀个GAGE12G：再次检索⼀下第⼀个GAGE4基因：对⽐⼀下⼆者：

这时会想，那我不⽤第⼀个了，⽤最后⼀个呗。事情依然没有这么简单：

1.3.2 第⼆个例⼦：ENSG00000063587

可以看到这个ENSG00000063587基因存在两个symbol：ZNF275和LOC105373378

> conv[389:392,]

ENSEMBL SYMBOL389 ENSG00000063515 GSC2390 ENSG00000063587 ZNF275391 ENSG00000063587 LOC105373378392 ENSG00000063854 HAGH

搜索⼀下就知道，这第⼀个symbol：ZNF275 是可靠的

好，⽤第⼀个也不对，⽤最后⼀个也不对，那应该怎么办呢？这时要考虑换⼀下注释数据库

ID转换⽅法不⽌⼀种，详见之前我写的：ID转换失败，爬起来继续high尤其是针对Ensemble ID时，⾃家的转换能获得更多的结果

library('biomaRt')

# ⽤useMart函数链接到⼈类的数据库

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【补充】特殊情况⼆：存在多对⼀这⼀⼩部分只是为了⽂章结构的完整性，这⾥的数据没有体现那么只存在⼀对多吗？并不是！还有多个Ensemble 对应⼀个symbol但其实看⼀下Ensemble ID在染⾊体上的位置就能知道：只有下图中加粗的那个ID才位于染⾊体的主体（chr19），其他的都来⾃haplotypic regions名词解释：Haplotyptic regions are defined by the Genome Reference Consortium (GRC) and are visible on all their genomeassemblies for human, mouse, and zebrafish. They can also be called ‘alternate locus’, ‘partial chromosomes’,and ‘alternate alleles’ or ‘assembly 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的结果，ids2是biomart的结果，df是原来的ID数据框conv2=conv[!duplicated(conv$ENSEMBL),] # 简单去重，只保留第⼀个df2=left_join(df,conv2,by=c('V1'='ENSEMBL')) # bitr结果df3=left_join(df,ids2,by=c('V1'='ensembl_gene_id')) # biomart结果> dim(df);dim(df2);dim(df3)[1] 64006 2[1] 64006 3[1] 64006 3# 上⾯的第⼀个例⼦：ENSG00000215269> df2[df2$V1=='ENSG00000215269',] V1 V2 SYMBOL63685 ENSG00000215269 ENSG00000283632 GAGE4> df3[df3$V1=='ENSG00000215269',] V1 V2 hgnc_symbol63685 ENSG00000215269 ENSG00000283632 GAGE12G看看df2中⼀些NA的情况：例如ENSG00000064489在bitr中没有，但在biomart中有df2[df2$V1=='ENSG00000064489',]df3[df3$V1=='ENSG00000064489',]image但能说biomart就⼀定⽐bitr好吗？显然不能这么说，例如ENSG00000036549就在bitr结果中有，⽽在biomart结果没有# 再看⼀个例⼦：df2[df2$V1=='ENSG00000036549',]df3[df3$V1=='ENSG00000036549',]并且搜索验证⼀下，ENSG00000036549 还真的有这个Symbol ID对应不过这⾥是GRCh37 (hg19)基因组版本如果是GRCh38（hg38）的话，结果就是：3 代码整合版上⾯的内容作为探索过程，理解做了什么就好df=read.csv(files[i],sep='\',header = F)## ⽤biomart转换id# ⽤useMart函数链接到⼈类的数据库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 IDif(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 IDif(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 为何不强强联合呢？

既然看到有少部分在bitr有注释，但biomart没有

那么，在biomart的基础上加上bitr的补充结果，岂不是更好？

# 就是看看hgnc_symbol这⼀列中显⽰ENSG0*****的基因在bitr中能不能转换成功，能的话就取代原来的结果## 整合bitr

# 先得到org.Hs.eg.db中的所有Ensemble IDorg_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 ⼩结

当转换Ensemble时，最好还是⽤⾃家的biomart，可以注释到更多的Ensemble ID，准确性也会更好⼀些，但速度可能会稍微

慢⼀点点；

如果基因ID数量不多，只是想快速看看转换后的ID，⽤bitr 也是⾜够的

最后，⽤了⼀个上午重新探索了⼀下ID转换，最后将bitr和biomart合⼆为⼀，增加了ID转换的准确性和丰富性，并且可以批处理多个⽂件，还增加了⼀些⼩修饰（⽐如每个⽂件转换的时间计算、转换结果输出前进⼀步检查，保证转换前后的结果每⾏的Ensemble ID⼀致）