基因ID转换方法总结

2021-06-19 本文已影响0人 Hayley笔记

1. clusterProfiler包，bitr()转换ID函数

clusterProfiler 是Y叔写的一个功能强大的R包，可以用来做各种富集分析，如GO、KEGG、DO（Disease Ontology analysis）、Reactome pathway analysis、GSEA富集分析等。还具有非常优秀的富集分析结果可视化功能。

1.1 bitr()的使用方法：

bitr(geneID, fromType, toType, OrgDb, drop = TRUE)

geneID：一个含有基因名的向量
orgDb：人类的注释包是org.Hs.eg.db，小鼠是org.Mm.eg.db
fromType：输入的基因名的类型
toType：需要转换成的类型，可以是多种类型，用大写character类型向量表示

1.2 查看org.Hs.eg.db中可以被选择/使用的类型

keytypes(org.Hs.eg.db)
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
 [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
[11] "GO"           "GOALL"        "IPI"          "MAP"          "OMIM"        
[16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
[21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"       "UNIGENE"     
[26] "UNIPROT"

1.3 示例

library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)

deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL")) #将转换后的结果加入矩阵，变成一列

需要注意的是，一些数据比如从TCGA上下载的数据，行名是带有版本号的，如下：

View(exp)
a = rownames(exp)
head(a)
# [1] "ENSG00000000003.13"
# [2] "ENSG00000000419.11"
# [3] "ENSG00000000457.12"
# [4] "ENSG00000000460.15"
# [5] "ENSG00000000938.11"
# [6] "ENSG00000000971.14"

这样的ENTREZID是不能被bitr()识别的，需要按点号分隔开，保留点号前面的部分，才能被bitr()识别。

library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
a = str_split(a,"\\.",simplify = T)[,1] #不带\\直接写"."是正则表达式任意字符的意思，需要\\来转义。
id = bitr(a,
          fromType = "ENSEMBL",
          toType = "SYMBOL",
          OrgDb = "org.Hs.eg.db")

2. 探针ID转化为gene symbol

2.1 获得探针和gene symbol的对应关系

方法一：根据gpl_number，从http://www.bio-info-trainee.com/1399.html中寻找对应的bioc_package，下载对应的包。
比如最常用的GPL570，对应hgu133plus2

if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db") #注意加上.db
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
#  [1] "hgu133plus2"             
#  [2] "hgu133plus2_dbconn"      
#  [3] "hgu133plus2_dbfile"      
#  [4] "hgu133plus2_dbInfo"      
#  [5] "hgu133plus2_dbschema"    
#  [6] "hgu133plus2.db"          
#  [7] "hgu133plus2ACCNUM"       
#  [8] "hgu133plus2ALIAS2PROBE"  
#  [9] "hgu133plus2CHR"          
# [10] "hgu133plus2CHRLENGTHS"   
# [11] "hgu133plus2CHRLOC"       
# [12] "hgu133plus2CHRLOCEND"    
# [13] "hgu133plus2ENSEMBL"      
# [14] "hgu133plus2ENSEMBL2PROBE"
# [15] "hgu133plus2ENTREZID"     
# [16] "hgu133plus2ENZYME"       
# [17] "hgu133plus2ENZYME2PROBE" 
# [18] "hgu133plus2GENENAME"     
# [19] "hgu133plus2GO"           
# [20] "hgu133plus2GO2ALLPROBES" 
# [21] "hgu133plus2GO2PROBE"     
# [22] "hgu133plus2MAP"          
# [23] "hgu133plus2MAPCOUNTS"    
# [24] "hgu133plus2OMIM"         
# [25] "hgu133plus2ORGANISM"     
# [26] "hgu133plus2ORGPKG"       
# [27] "hgu133plus2PATH"         
# [28] "hgu133plus2PATH2PROBE"   
# [29] "hgu133plus2PFAM"         
# [30] "hgu133plus2PMID"         
# [31] "hgu133plus2PMID2PROBE"   
# [32] "hgu133plus2PROSITE"      
# [33] "hgu133plus2REFSEQ"       
# [34] "hgu133plus2SYMBOL"       
# [35] "hgu133plus2UNIGENE"      
# [36] "hgu133plus2UNIPROT"  
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
#    probe_id symbol
# 1   1053_at   RFC2
# 2    117_at  HSPA6
# 3    121_at   PAX8
# 4 1255_g_at GUCA1A
# 5   1316_at   THRA
# 6   1320_at PTPN21

方法二：读取GPL平台的soft文件，按列取子集
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570

a = getGEO(gpl_number,destdir = ".")
b = a@dataTable@table
colnames(b)
ids2 = b[,c("ID","Gene Symbol")]
colnames(ids2) = c("probe_id","symbol")
ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"///"),]

方法三：官网下载，文件读取
参考：http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus
方法四：自主注释
参考：https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA

2.2 加上探针注释

#为deg数据框添加几列（deg是做完差异分析后，行为probe_id，列为logFC等差异分析结果的矩阵）
#1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)
#2.加上探针注释
#多个探针对应一个基因：按照基因去重复
table(!duplicated(ids$probe_id)) 
table(!duplicated(ids$symbol)) #结果显示很多基因对应了多个探针
#按symbol列去重，常见标准有3个：最大值/平均值/随机去重
#随机去重，另两个见zz.去重方式.R
ids = ids[!duplicated(ids$symbol),] #保留重复值中第一个出现的

deg <- inner_join(deg,ids,by="probe_id")
head(deg)
nrow(deg)

3. gtf文件获得id转换信息

对于TCGA上下载的数据，其使用的geneid是ensembl id，做完差异分析后需要把结果中的每个ensembl id转换为对应的symbol和类型(mRNA/lncRNA或其它)。而gtf文件中就包含了我们需要的gene symbol和ensembl id的对应关系。

思路：
1. 找到TCGA数据对应的参考基因组注释版本。
2. 下载该版本的参考基因组注释文件，提取ensembl id 与symbol的对应关系及每个基因的gene type信息。
3. 可以将symbol和gene type 用merge添加到差异分析结果中，也可以在差异分析前先转换矩阵的行名。

3.1 找参考基因组版本

可以看到，使用的参考基因组版本是genecode的v22。（版本很多，这个是14年的版本了）

在gtf文件里并不是直接分出了lncRNA，需要找gtf文件里对biotype的说明，不看不知道，一看发现这是一个很长的表格。
进入https://www.gencodegenes.org，点击Documentation下的Biotypes，
其中对lncRNA的说明是：Generic long non-coding RNA biotype that replaced the following biotypes: 3prime_overlapping_ncRNA, antisense, bidirectional_promoter_lncRNA, lincRNA, macro_lncRNA, non_coding, processed_transcript, sense_intronic and sense_overlapping.
所以需要将genetype里这些类型对应的行挑出来，就是lncRNA了。然后与表达矩阵行名进行匹配替换，就可以分别得到mRNA和lncRNA的矩阵了。

3.2 转换

#step1:读取并探索gtf文件----
#BiocManager::install("rtracklayer")
library(rtracklayer)
gtf = rtracklayer::import("gencode.v22.annotation.gtf")
class(gtf)
gtf = as.data.frame(gtf);dim(gtf)
colnames(gtf)
table(gtf$type)
#step2:先筛选出gene对应的行
gtf_gene = gtf[gtf$type=="gene",]
save(gtf_gene,file = "gtf_gene.Rdata")
table(rownames(deg) %in% gtf_gene$gene_id) #deg是需要进行id转换的矩阵
# FALSE  TRUE 
#     3 30345

an = gtf_gene[,c("gene_name","gene_id","gene_type")] #取出gtf文件中这三列
head(an)
#       gene_name           gene_id                          gene_type
# 1       DDX11L1 ENSG00000223972.5 transcribed_unprocessed_pseudogene
# 13       WASH7P ENSG00000227232.5             unprocessed_pseudogene
# 26    MIR6859-3 ENSG00000278267.1                              miRNA
# 29 RP11-34P13.3 ENSG00000243485.3                            lincRNA
# 37    MIR1302-9 ENSG00000274890.1                              miRNA
# 40      FAM138A ENSG00000237613.2                            lincRNA
deg = merge(deg,an,by.x = "row.names",by.y = "gene_id")

# mRNA和lncRNA总共有多少个？
lnc = c("3prime_overlapping_ncRNA", "antisense", "bidirectional_promoter_lncRNA", "lincRNA", "macro_lncRNA", "non_coding", "processed_transcript", "sense_intronic" , "sense_overlapping")
k1 = gtf_gene$gene_type %in% lnc;table(k1) #gtf中共14826个lnc
# k1
# FALSE  TRUE 
# 45657 14826
k2 = gtf_gene$gene_type == "protein_coding";table(k2) #gtf中共19814个mRNA
# k2
# FALSE  TRUE 
# 40669 19814

# deg中有多少mRNA和lncRNA?
k3 = deg$gene_type %in% lnc;table(k3) #deg中共7501个lnc
# k3
# FALSE  TRUE 
# 22844  7501
k4 = deg$gene_type =="protein_coding";table(k4) #deg中共17464个mRNA
# k4
# FALSE  TRUE 
# 12881 17464

# 差异的mRNA和lncRNA 各有多少
k5 = deg$change !="NOT"
table(k3&k5)
# FALSE  TRUE   #有396个差异性lnc
# 29949   396
table(k4&k5) 
# FALSE  TRUE  #有1084个差异性mRNA
# 29261  1084

表达矩阵的行名id转换

load("gtf_gene.Rdata")
an = gtf_gene[,c("gene_name","gene_id","gene_type")]
exp = exp[rownames(exp) %in% an$gene_id,] #从表达矩阵中去除不存在于gtf中的gene名
an = an[match(rownames(exp),an$gene_id),] #将an的顺序调整成和表达矩阵行名一致
identical(an$gene_id,rownames(exp)) #检查是否一致
# [1] TRUE
k = !duplicated(an$gene_name);table(k) #因为要用gene_symbol做行名，而行名是不能重复的，所以要先去重
# k
# FALSE  TRUE 
#   193 30152
an = an[k,]
exp = exp[k,]
rownames(exp) = an$gene_name

常用的转换需求总结：

GO分析：需要将gene symbol转化为ENTREZID

分析芯片数据：需要将探针ID转化为gene symbol

将TCGA上下载的矩阵中的ENTREZID转化为gene symbol（1和3实质上是一样的），可以直接转换矩阵，也可以在做完差异分析后，把gene symbol作为一列连接在差异表达矩阵上。

代码来自2021生信技能树数据挖掘课