转录组完整的ID转换：biomaRt和gtf

2021-05-08 本文已影响0人冻春卷

ID转换过程中会有基因丢失这件事情，在部分做干实验的人看来，是非常正常的。但不做干实验的湿实验人知道这件事，心里可能会有疙瘩。为什么要丢失呢？为什么丢失了不补回来？

在分析结果中查阅重要的基因并绘图时发现竟然无此基因的数据，或者是绘图时row.name的部分丢失，会让我们的强迫症爆发。真心希望有一个完整的ID转换方法。在这之前我使用的是ENSEMBL的官方工具BiomaRt。目前BiomaRt囊括了多达208种物种的ID数据，可以说是又大由全。然而即使我做的是人类基因ID转换，仍然还是发现有基因丢失，是可忍孰不可忍。后经过多方查阅，决定使用参考基因组的gtf文件进行完整的ID转换。下面我们就来对比一下，BiomaRt和gtf的ID转换结果。

（注：使用BiomaRt进行ID转化出现ID丢失的原因不一定是因为R包不好，有可能是数据的版本不同）

数据准备

在这里使用本人的转录组数据用于测试，数据经过上游处理后简单整理数据如下：

> 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.3 简单查看人类ID数据集

human_dataset@attributes$name[1:20] ## 查看一下都有什么名字

可以看到数据集中包含：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转换文件