轉(zhuǎn)錄組完整的ID轉(zhuǎn)換：biomaRt和gtf

ID轉(zhuǎn)換過(guò)程中會(huì)有基因丟失這件事情李破，在部分做干實(shí)驗(yàn)的人看來(lái)宠哄，是非常正常的。但不做干實(shí)驗(yàn)的濕實(shí)驗(yàn)人知道這件事嗤攻，心里可能會(huì)有疙瘩琳拨。為什么要丟失呢？為什么丟失了不補(bǔ)回來(lái)屯曹？

在分析結(jié)果中查閱重要的基因并繪圖時(shí)發(fā)現(xiàn)竟然無(wú)此基因的數(shù)據(jù)狱庇，或者是繪圖時(shí)row.name的部分丟失，會(huì)讓我們的強(qiáng)迫癥爆發(fā)恶耽。真心希望有一個(gè)完整的ID轉(zhuǎn)換方法密任。在這之前我使用的是ENSEMBL的官方工具BiomaRt。目前BiomaRt囊括了多達(dá)208種物種的ID數(shù)據(jù)偷俭，可以說(shuō)是又大由全浪讳。然而即使我做的是人類基因ID轉(zhuǎn)換，仍然還是發(fā)現(xiàn)有基因丟失涌萤，是可忍孰不可忍淹遵。后經(jīng)過(guò)多方查閱，決定使用參考基因組的gtf文件進(jìn)行完整的ID轉(zhuǎn)換负溪。下面我們就來(lái)對(duì)比一下，BiomaRt和gtf的ID轉(zhuǎn)換結(jié)果川抡。

（注：使用BiomaRt進(jìn)行ID轉(zhuǎn)化出現(xiàn)ID丟失的原因不一定是因?yàn)镽包不好辐真，有可能是數(shù)據(jù)的版本不同）

數(shù)據(jù)準(zhǔn)備

在這里使用本人的轉(zhuǎn)錄組數(shù)據(jù)用于測(cè)試崖堤，數(shù)據(jù)經(jīng)過(guò)上游處理后簡(jiǎn)單整理數(shù)據(jù)如下：

> dim(count) # 檢查數(shù)據(jù)框大小
[1] 58884     9

> table(duplicated(count$geneid)) # 檢查是否有重復(fù)的ensemble ID
FALSE 
58884

由上信息可見(jiàn)，我的轉(zhuǎn)錄組數(shù)據(jù)經(jīng)過(guò)比對(duì)后得到58884個(gè)“基因”，因此gene symbol轉(zhuǎn)換要越接近這個(gè)數(shù)值越好胯甩。

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 選定人類數(shù)據(jù)集

listMarts()  ## 查看目前提供的數(shù)據(jù)庫(kù)
# formal class mart 
my_mart<- useMart("ENSEMBL_MART_ENSEMBL") # 選定數(shù)據(jù)集
## 查看數(shù)據(jù)集
datasets<- listDatasets(my_mart)
datasets
dim(datasets)  # [1] 202   3  目前有208個(gè)數(shù)據(jù)集（物種ID信息）
# 設(shè)定人類ID數(shù)據(jù)集
human_dataset<- useDataset("hsapiens_gene_ensembl",mart = my_mart) # 約1.3 M

1.3 簡(jiǎn)單查看人類ID數(shù)據(jù)集

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

可以看到數(shù)據(jù)集中包含：ensemble ID和gene id的轉(zhuǎn)換昧廷，基因轉(zhuǎn)錄本ID等等內(nèi)容。簡(jiǎn)單查閱一遍蜡豹，可知“ensembl_gene_id”是我們想要的內(nèi)容麸粮。

1.4 設(shè)定attributes參數(shù)

attributes參數(shù)定義了四個(gè)輸出項(xiàng)ensembl_gene_id，chromosome_name镜廉， hgnc_smbol以及hgnc_id弄诲。

count_value<- count$geneid  # 設(shè)定需要轉(zhuǎn)換的ID
attr1<- c("ensembl_gene_id","chromosome_name","hgnc_symbol","hgnc_id") # 設(shè)定參數(shù)
count_ID<- getBM(attributes = attr1,
          filters = "ensembl_gene_id",
          values = count_value, 
          mart = human_dataset)

輸出結(jié)果如下：可見(jiàn)attribute定義的四個(gè)輸出項(xiàng)。

ID conversion

1.5 查看ID轉(zhuǎn)換的完整度

> dim(count_ID)  # 查看ID轉(zhuǎn)換結(jié)果
[1] 58666     4  # 有58666個(gè)ensemble ID完成了比對(duì)

可見(jiàn)有58666個(gè)ensemble ID完成了轉(zhuǎn)換，這比原始數(shù)據(jù)中的58884少了218個(gè)ensemble ID齐遵。但這僅僅是ensemble ID的轉(zhuǎn)換結(jié)果寂玲，我們還需要查看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

小結(jié)：biomaRt ID轉(zhuǎn)換會(huì)出現(xiàn)大量的gene symbol的丟失梗摇，具體原因可能是已經(jīng)將重復(fù)的gene symbol去除（有多個(gè)ensemble ID對(duì)應(yīng)gene symbol的情況）拓哟。

2. 使用基因組的gtf注釋文件

在自己用于比對(duì)的參考基因組那里可以找到相應(yīng)的注釋文件，我使用的是hg38版本

hg38

2.1 配置R包

BiocManager::install("rtracklayer")
library(rtracklayer)

進(jìn)行相應(yīng)的文件設(shè)置

# input 
gff <- readGFF("Homo_sapiens.GRCh38.96.gtf")
head(gff)

mapid <- gff[gff$type == "gene", c("gene_id", "gene_name")]
head(mapid)
#gff文件很大伶授，用掉就刪掉

# 判斷我們要轉(zhuǎn)換的基因是不是都在
table(count$geneid %in%mapid$gene_id)
mapid <- read.csv("ensemble2symbol.csv")
head(mapid)

whether all

至此断序，mapid文件則是存儲(chǔ)ID轉(zhuǎn)換信息的文件。

2.2 合并

查看mapid文件并與需要進(jìn)行ID轉(zhuǎn)換的數(shù)據(jù)合并

# 查看gene symbol是否有空值
table( mapid$gene_name == "")

# 用merge進(jìn)行合并
df <- merge(count, mapid, by.x="geneid", by.y="gene_id")

# 先判斷是不是存在重復(fù)的基因名糜烹，如果存在重復(fù)违诗，先考慮去重
table(duplicated(df$gene_name))

df <- df[!duplicated(df$gene_name),]
table(duplicated(df$gene_name))

image

由上可知，使用gtf文件進(jìn)行ID轉(zhuǎn)換會(huì)得到最全的結(jié)果疮蹦。但最全的結(jié)果是否最適合后續(xù)分析诸迟，我們還需要進(jìn)一步考察。

2.3 查看gtf轉(zhuǎn)換文件