基因id转换方法大盘点

0.R包的安装和加载

if(!require(mygene)) BiocManager::install("mygene")
library(mygene)
if(!require(biomaRt))BiocManager::install("biomaRt")
library(biomaRt)
if(!require(AnnoProbe))devtools::install_github("jmzeng1314/AnnoProbe")
if(!require(rtracklayer))BiocManager::install("rtracklayer")
if(!require(clusterProfiler))BiocManager::install("clusterProfiler")
if(!require(org.Hs.eg.db))BiocManager::install("org.Hs.eg.db")

1.mygene ⭐⭐⭐⭐⭐

优点：一下搞定多种基因id

这是本文的开端，本来想介绍这一个方法。结果和豆豆讨论着，又揪出了后面的一大串方法。为了避免大家选择困难症，打个星级吧。

xli <-  c('DDX26B','CCDC83',  'MAST3',
          'RPL11', 'ZDHHC20',  'LUC7L3',  
          'SNORD49A',  'CTSH', 'ACOT8')
tmp = queryMany(xli, scopes="symbol", 
          fields=c("uniprot", "ensembl.gene", "reporter"), 
          species="human")

## Finished
## Pass returnall=TRUE to return lists of duplicate or missing query terms.

tmp = as.data.frame(tmp)
tmp[1:3,1:5]

##    query notfound   X_id  X_score    ensembl.gene
## 1 DDX26B     TRUE   <NA>       NA            <NA>
## 2 CCDC83       NA 220047 89.64873 ENSG00000150676
## 3  MAST3       NA  23031 89.68144 ENSG00000099308

colnames(tmp)

##  [1] "query"                   "notfound"               
##  [3] "X_id"                    "X_score"                
##  [5] "ensembl.gene"            "reporter.GNF1H"         
##  [7] "reporter.HG.U133_Plus_2" "reporter.HTA.2_0"       
##  [9] "reporter.HuEx.1_0"       "reporter.HuGene.1_1"    
## [11] "reporter.HuGene.2_1"     "reporter.HG.U95Av2"     
## [13] "reporter.HG.U95B"        "uniprot.Swiss.Prot"     
## [15] "uniprot.TrEMBL"

2.bitr ⭐⭐⭐⭐

优点：方便快捷。

缺点：比起上面那个，一次只能转换出一种id哦。

library(clusterProfiler)
library(org.Hs.eg.db)
x = bitr(xli,fromType = "SYMBOL",
         toType = "ENSEMBL",OrgDb = "org.Hs.eg.db")
head(x)

##     SYMBOL         ENSEMBL
## 2   CCDC83 ENSG00000150676
## 3    MAST3 ENSG00000099308
## 4    RPL11 ENSG00000142676
## 5  ZDHHC20 ENSG00000180776
## 6   LUC7L3 ENSG00000108848
## 7 SNORD49A ENSG00000277370

3.AnnoProbe ⭐⭐⭐⭐

优点:方便快捷，代码简单，风格简洁，此处省略一万字。

IDs <- c("DDX11L1", "MIR6859-1", "OR4G4P", "OR4F5")
ID_type = "SYMBOL"
annoGene(IDs, ID_type)

##    SYMBOL                           biotypes         ENSEMBL  chr start   end
## 1 DDX11L1 transcribed_unprocessed_pseudogene ENSG00000223972 chr1 11869 14409
## 5  OR4G4P             unprocessed_pseudogene ENSG00000268020 chr1 52473 53312
## 7   OR4F5                     protein_coding ENSG00000186092 chr1 65419 71585

缺点？这是我老板写的包，不能写缺点，哦不，没有缺点。

4.biomaRt ⭐⭐

也是个常用方法，但我从来不用。豆豆让我加上哦。优点：id也是比较齐全 缺点：代码有点啰嗦。我jio的可以被bitr或者mygene无缝替代。

library(biomaRt)

attr <- c("hgnc_symbol","uniprotswissprot",'chromosome_name','start_position','end_position','strand')

# human use: hsapiens_gene_ensembl
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl") 

ids <- getBM(attributes = attr,
             filters = "hgnc_symbol",
             values = xli,
             mart = ensembl,
             useCache = F)
head(ids)

##   hgnc_symbol uniprotswissprot chromosome_name start_position end_position
## 1     ZDHHC20                               13       21372571     21459370
## 2     ZDHHC20           Q5W0Z9              13       21372571     21459370
## 3       ACOT8           O14734              20       45841721     45857405
## 4       ACOT8                               20       45841721     45857405
## 5      LUC7L3           O95232              17       50719565     50756219
## 6      LUC7L3                               17       50719565     50756219
##   strand
## 1     -1
## 2     -1
## 3     -1
## 4     -1
## 5      1
## 6      1

5.NCBI有数据 ⭐⭐

这个我也不怎么用，需要下载数据。但他是官方出品。

https://ftp.ncbi.nih.gov/gene/DATA/

举个栗子

6.gtf文件⭐⭐⭐⭐⭐

如果知道比对基因时使用的gtf版本，那用gtf来找注释是最好的。例如tcga的转录组数据用的genecodev22版本ensemblid

缺点：版本十分具体，普适性稍差一点。

library(rtracklayer)
g = import("gencode.v22.annotation.gtf.gz")
colnames(as.data.frame(g))

##  [1] "seqnames"                 "start"                   
##  [3] "end"                      "width"                   
##  [5] "strand"                   "source"                  
##  [7] "type"                     "score"                   
##  [9] "phase"                    "gene_id"                 
## [11] "gene_type"                "gene_status"             
## [13] "gene_name"                "level"                   
## [15] "havana_gene"              "transcript_id"           
## [17] "transcript_type"          "transcript_status"       
## [19] "transcript_name"          "tag"                     
## [21] "transcript_support_level" "havana_transcript"       
## [23] "exon_number"              "exon_id"                 
## [25] "ont"                      "protein_id"              
## [27] "ccdsid"

g = as.data.frame(g)
g = g[,c("gene_name","gene_id","gene_type")]
head(g)

##   gene_name           gene_id                          gene_type
## 1   DDX11L1 ENSG00000223972.5 transcribed_unprocessed_pseudogene
## 2   DDX11L1 ENSG00000223972.5 transcribed_unprocessed_pseudogene
## 3   DDX11L1 ENSG00000223972.5 transcribed_unprocessed_pseudogene
## 4   DDX11L1 ENSG00000223972.5 transcribed_unprocessed_pseudogene
## 5   DDX11L1 ENSG00000223972.5 transcribed_unprocessed_pseudogene
## 6   DDX11L1 ENSG00000223972.5 transcribed_unprocessed_pseudogene

b = g[g$gene_name %in% xli,]
head(b)

##       gene_name            gene_id      gene_type
## 35253     RPL11 ENSG00000142676.11 protein_coding
## 35254     RPL11 ENSG00000142676.11 protein_coding
## 35255     RPL11 ENSG00000142676.11 protein_coding
## 35256     RPL11 ENSG00000142676.11 protein_coding
## 35257     RPL11 ENSG00000142676.11 protein_coding
## 35258     RPL11 ENSG00000142676.11 protein_coding