ID转换过程中会有基因丢失这件事情，在部分做干实验的人看来，是非常正常的。但不做干实验的湿实验人知道这件事，心里可能会有疙瘩。为什么要丢失呢？为什么丢失了不补回来？
在分析结果中查阅重要的基因并绘图时发现竟然无此基因的数据，或者是绘图时row.name的部分丢失，会让我们的强迫症爆发。真心希望有一个完整的ID转换方法。在这之前我使用的是ENSEMBL的官方工具BiomaRt。目前BiomaRt囊括了多达208种物种的ID数据，可以说是又大由全。然而即使我做的是人类基因ID转换，仍然还是发现有基因丢失，是可忍孰不可忍。后经过多方查阅，决定使用参考基因组的gtf文件进行完整的ID转换。下面我们就来对比一下，BiomaRt和gtf的ID转换结果。
（注：使用BiomaRt进行ID转化出现ID丢失的原因不一定是因为R包不好，有可能是数据的版本不同）
数据准备
在这里使用本人的转录组数据用于测试，数据经过上游处理后简单整理数据如下：
0
> dim(count) # 检查数据框大小
[1] 58884 9
> table(duplicated(count$geneid)) # 检查是否有重复的ensemble ID
FALSE
58884
由上信息可见，我的转录组数据经过比对后得到58884个“基因”，因此gene symbol转换要越接近这个数值越好。
1. BiomaRt
1.1 install BiomaRt
安装代码可以从bioconductor http://www.bioconductor.org/packages/release/bioc/html/biomaRt.html上查阅
if (!requireNamespace("biomaRt", quietly = TRUE) ){
BiocManager::install("biomaRt")
}
library(biomaRt) # 激活R包
1.2 选定人类数据集
listMarts() ## 查看目前提供的数据库
# formal class mart
my_mart<- useMart("ENSEMBL_MART_ENSEMBL") # 选定数据集
## 查看数据集
datasets<- listDatasets(my_mart)
datasets
dim(datasets) # [1] 202 3 目前有208个数据集（物种ID信息）
# 设定人类ID数据集
human_dataset<- useDataset("hsapiens_gene_ensembl",mart = my_mart) # 约1.3 M
1
1.3 简单查看人类ID数据集
human_dataset@attributes$name[1:20] ## 查看一下都有什么名字
2
可以看到数据集中包含：ensemble ID和gene id的转换，基因转录本ID等等内容。简单查阅一遍，可知“ensembl_gene_id”是我们想要的内容。
1.4 设定attributes参数
attributes参数定义了四个输出项ensembl_gene_id，chromosome_name， hgnc_smbol以及hgnc_id。
count_value<- count$geneid # 设定需要转换的ID
attr1<- c("ensembl_gene_id","chromosome_name","hgnc_symbol","hgnc_id") # 设定参数
count_ID<- getBM(attributes = attr1,
filters = "ensembl_gene_id",
values = count_value,
mart = human_dataset)
输出结果如下：可见attribute定义的四个输出项。
ID conversion
1.5 查看ID转换的完整度
> dim(count_ID) # 查看ID转换结果
[1] 58666 4 # 有58666个ensemble ID完成了比对
可见有58666个ensemble ID完成了转换，这比原始数据中的58884少了218个ensemble ID。但这仅仅是ensemble ID的转换结果，我们还需要查看gene symbol（也就是表格中的hgnc_symbol）的完成度。
table(is.na(count_ID$ensembl_gene_id))
table(is.na(count_ID$hgnc_symbol))
table( count_ID$hgnc_symbol == "")
table(duplicated(count_ID$hgnc_symbol))
biomaRt
小结：biomaRt ID转换会出现大量的gene symbol的丢失，具体原因可能是已经将重复的gene symbol去除（有多个ensemble ID对应gene symbol的情况）。
2. 使用基因组的gtf注释文件
在自己用于比对的参考基因组那里可以找到相应的注释文件，我使用的是hg38版本
hg38
2.1 配置R包
BiocManager::install("rtracklayer")
library(rtracklayer)
进行相应的文件设置
# input
gff <- readGFF("Homo_sapiens.GRCh38.96.gtf")
head(gff)
mapid <- gff[gff$type == "gene", c("gene_id", "gene_name")]
head(mapid)
#gff文件很大，用掉就删掉
# 判断我们要转换的基因是不是都在
table(count$geneid %in%mapid$gene_id)
mapid <- read.csv("ensemble2symbol.csv")
head(mapid)
whether all
至此，mapid文件则是存储ID转换信息的文件。
2.2 合并
查看mapid文件并与需要进行ID转换的数据合并
# 查看gene symbol是否有空值
table( mapid$gene_name == "")
# 用merge进行合并
df <- merge(count, mapid, by.x="geneid", by.y="gene_id")
# 先判断是不是存在重复的基因名，如果存在重复，先考虑去重
table(duplicated(df$gene_name))
df <- df[!duplicated(df$gene_name),]
table(duplicated(df$gene_name))
image
由上可知，使用gtf文件进行ID转换会得到最全的结果。但最全的结果是否最适合后续分析，我们还需要进一步考察。
2.3 查看gtf转换文件
image
由上图可见，有一些ensemble ID其实是对应同一个gene symbol，只不过这些gene symbol有不同的版本号，可见名字后面的小数点。于是我们会提出疑问，该如何处理不同版本的gene symbol呢？暂未明确，须待考察。
3. 小结
总体来说，使用gtf注释文件进行ID转换是最完整的ID转换。关于ID转换的事情困扰我已久，或许这本不是什么大问题，但仍是一个问题。
网友评论