152-ID转换失败,爬起来继续high

刘小泽写于19.12.9 昨天遇到一个基因ID转换的问题,才发现原来ID转换不是这么简单…

0 前言

拿到一个接近60000 x 6000的原始表达矩阵,行是基因,列是样本。并且行名都是Ensembl ID,还带着版本号。

我的第一直觉就是先把ID转换成Symbol,保证后续操作的易读性

下面👇是我对这个原始矩阵的操作

1 第一步 对原始表达矩阵进行过滤

很简单粗暴,直接去掉那些在所有样本不表达的基因

> dim(origin_count)
[1] 58347  6000
# 过滤后的矩阵df
df <- origin_count[rowSums(origin_count)>0,]
> dim(df)
[1] 43523  6295

2 第二步 修改原始行名的Ensembl ID

这里要先说一下Ensembl 的ID

2.1 我们常见的Ensembl ID比如:

ENSG00000279928.1 ENSG00000279457.2 详细说明:https://asia.ensembl.org/Help/Faq?id=488

  • 其中,ENS是固定字符,表示它是属于Ensembl数据库的。默认物种是人,如果是小鼠就要用ENSMUS开头,关于物种代码: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位数字部分如00000279928 表示基因真正的编号
  • 小数点后不一定每个都有,表示这个ID在数据库中做了几次变更,比如.1做了1次变更,在分析时需要去除

至于怎么去掉小数点后面的信息也很简单,之前花花介绍过( 去掉ensembl id的version信息

2.2 小心两个坑

🤡 第一坑

不过这里,仅仅是简单去掉还不够,因为我们的目的是让新的Ensembl ID当做行名,但是去掉小数点后会出现一个问题:重复值 ,而行名是不允许重复值出现的

因此,需要保证新的Ensembl ID中如果出现重复值,那么就只保留第一个

🤡 第二坑

之前原始的Ensembl ID是和行名的长度一致的,现在得到了去重复的Ensembl ID,那么由于长度不等,因此不能一一对应到行名上。

要根据之前Ensembl ID的重复情况,去掉原始数据框对应的行(不过不要担心,重复的值没有多少,所以也去不掉多少行,对结果没有太大影响;并且既然是重复的行,它们的表达量应该也是类似的,所以取其中一个即可)

2.3 那么整体的操作如下:

head(rownames(df),3)
# "ENSG00000000003.14" "ENSG00000000005.5"  "ENSG00000000419.12"

# 去掉小数点后的gene name(gn)
gn <- stringr::str_split(rownames(df),"\\.",simplify = T)[,1]
# 根据gn的重复情况,df也进行对应的去除
df <- df[!duplicated(gn),]
# 此时df的行数才和gn相等
rownames(df) <- gn[!duplicated(gn)]

> dim(df) # 看到相比之前(43523行),过滤掉了39个重复行
[1] 43484  6295
> head(rownames(df),3)
[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419"

3 第三步 基因ID转换

3.1 使用org.Hs.eg.db + mapIds

suppressMessages(library(org.Hs.eg.db))
symb <- mapIds(org.Hs.eg.db, keys=rownames(df), 
               keytype="ENSEMBL", column="SYMBOL")

> length(rownames(df))
[1] 43484
> length(na.omit(symb))
[1] 22965

可以看到基本包含了一半的NA值

3.2 使用org.Hs.eg.db + bitr

其实结果一样,因为都是使用的org.Hs.eg.db数据库,这里只是列一下方法

library(clusterProfiler)
library(org.Hs.eg.db)
gene_tr <- bitr(rownames(df), fromType = "ENSEMBL",
                toType = c("SYMBOL"),
                OrgDb = org.Hs.eg.db)
# 47.19% of input gene IDs are fail to map...

> length(na.omit(symb_2$SYMBOL))
[1] 23196

可以看到比mapIds能多得到一些基因

疑问

**这么多的NA基因,难道真的不存在对应的Symbol基因名吗?**这个问题值得思考🤔

很多时候,使用了两种工具,得到差不多的结果,可能就认为事已至此了。但是换种角度思考问题,如果问题不是出在工具呢?如果是数据库呢?

好,接下来就看看没有对应Symbol的那些Ensembl ID
not_mapped_ensembl <- names(symb[is.na(symb)])

> length(not_mapped_ensembl)
[1] 20519

> head(not_mapped_ensembl,3)
[1] "ENSG00000002079" "ENSG00000018607" "ENSG00000020219"

那么随便挑一个(例如第一个)直接复制、搜索,看看到底有没有基因名

搜索的结果让人大吃一惊,竟然它真的有基因名,而且还在HGNC、Ensembl、Genecards等数据库有收录

看到Ensembl 列出的该基因的信息很全,不由得想到,基因ID转换不止这两种工具,还有一种:Ensembl自家的biomaRt,尝试一下

3.3 使用hsapiens_gene_ensembl + biomaRt

看它的帮助文档就很清楚它怎么操作,过程也很简单

library("biomaRt")
# 用useMart函数链接到人类的数据库
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
# 除以以外还有很多:使用 listDatasets(ensembl) 查看

# 依然是要列出查询的gene ID,以及要转换的格式
# 【biomaRt强大之处就在于它的转换格式种类非常之丰富】
attributes <- listAttributes(ensembl)
View(attributes) # 查看转换格式

# 列出要转换的ID(比如前面匹配失败的3个基因)
value <- head(not_mapped_ensembl,3)
# 列出当前格式和要比对的格式(当前是ensembl_gene_id;我要比对hgnc_symbol和chromosome_name)
attr <- c("ensembl_gene_id","hgnc_symbol","chromosome_name")

# 最后运行就好了(其中filters指定当前基因 ID类型)
ids <- getBM(attributes = attr,
             filters = "ensembl_gene_id",
             values = value,
             mart = ensembl)

最后看一下结果,而且这个结果扩展性极强:

> ids
  ensembl_gene_id  hgnc_symbol chromosome_name
1 ENSG00000002079        MYH16               7
2 ENSG00000018607       ZNF806               2
3 ENSG00000020219      CCT8L1P               7
biomaRt的小问题

当提供的基因数量太多时,容易断线;它的检索速度非常慢!【不过用来应急还是可以的】

心得

当使用某个或某几个工具得到的结果感觉不满意时,不要只考虑工具的优劣,还要思考它们在后台帮我们做了什么。

比如这里基因ID的转换,靠的其实就是背后的数据库,工具虽有差异但没那么大(比如bitrmapIds多了几个基因)。有时一条路走不通,换个工具“背后的”数据库未尝不是个好的解决办法!

Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related