ID转换那点事

在进行测序数据下游分析的时候常常需要用到不同的数据库，而这些数据库的分析的输入文件经常是各有区别，因此我们常常需要将各类ID（ensamble ID、Entrez ID、gene symbol等）进行转换。这里我们将简单介绍使用编程进行ID转换的方法。

芯片ID注释

芯片注释主要有两种方法：使用R包进行注释或在数据库中下载芯片平台信息然后手动注释，这两种方法其实基本都是相同的，都是利用芯片平台信息进行注释，区别只在于一个是利用R包API而另一个则是自己写代码。

载入数据：

> library('GEOquery')
> gse20986 <- getGEO("GSE20986", destdir=".")  # 下载处理好的表达矩阵
> exprmt <- exprs(gse20986[[1]]);head(exprmt);class(exprmt)
          GSM524662 GSM524663 GSM524664 GSM524665 GSM524666 GSM524667 GSM524668 GSM524669 GSM524670 GSM524671 GSM524672 GSM524673
1007_s_at  2.867331  3.941245  4.003509  3.230621  3.546867  2.507474  3.040758  2.456892  2.404267  3.568920  3.307409  3.225964
1053_at   10.503022  8.320770  9.523183 10.645904  8.140238 10.228986 11.019461 10.447158 10.799832  9.765222 10.069229 10.089496
117_at     2.702517  2.749144  3.226281  2.788124  4.664960  2.722543  2.745330  4.500872  2.702614  3.015139  2.768683  2.853809
121_at     3.052316  3.057630  3.056844  3.697705  3.057630  3.057630  3.057630  3.057630  3.072568  3.057630  3.057630  3.057630
1255_g_at  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998
1294_at    5.360226  4.465648  5.914568  5.558741  5.241448  5.037804  4.784192  4.978118  4.713441  4.303196  4.452989  4.766019
[1] "matrix"
# 将矩阵转换为数据框方便后面注释操作，并将行名称装换为ID
> exprmt <- as.data.frame(exprmt); exprmt$ID <- rownames(exprmt); head(exprmt)
          GSM524662 GSM524663 GSM524664 GSM524665 GSM524666 GSM524667 GSM524668 GSM524669 GSM524670 GSM524671 GSM524672 GSM524673        ID
1007_s_at  2.867331  3.941245  4.003509  3.230621  3.546867  2.507474  3.040758  2.456892  2.404267  3.568920  3.307409  3.225964 1007_s_at
1053_at   10.503022  8.320770  9.523183 10.645904  8.140238 10.228986 11.019461 10.447158 10.799832  9.765222 10.069229 10.089496   1053_at
117_at     2.702517  2.749144  3.226281  2.788124  4.664960  2.722543  2.745330  4.500872  2.702614  3.015139  2.768683  2.853809    117_at
121_at     3.052316  3.057630  3.056844  3.697705  3.057630  3.057630  3.057630  3.057630  3.072568  3.057630  3.057630  3.057630    121_at
1255_g_at  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998  2.278998 1255_g_at
1294_at    5.360226  4.465648  5.914568  5.558741  5.241448  5.037804  4.784192  4.978118  4.713441  4.303196  4.452989  4.766019   1294_at

使用annotate包注释芯片

在数据集GSE20986页面可以看到芯片平台信息：

芯片平台信息

当然也可以直接通过代码获得芯片平台信息：

> gse20986
$GSE20986_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 54675 features, 12 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM524662 GSM524663 ... GSM524673 (12 total)
  varLabels: title geo_accession ... tissue:ch1 (34 total)
  varMetadata: labelDescription
featureData
  featureNames: 1007_s_at 1053_at ... AFFX-TrpnX-M_at (54675 total)
  fvarLabels: ID GB_ACC ... Gene Ontology Molecular Function (16 total)
  fvarMetadata: Column Description labelDescription
experimentData: use 'experimentData(object)'
Annotation: GPL570

知道了芯片平台信息，接下来我们就可以使用annotate包进行注释了。
先查看一下有没有重复的芯片ID:

> dupID <- exprmt[duplicated(exprmt$ID),]; head(dupID); dim(dupID)
 [1] GSM524662 GSM524663 GSM524664 GSM524665 GSM524666 GSM524667 GSM524668 GSM524669 GSM524670 GSM524671 GSM524672 GSM524673 ID       
<0 行> (或0-长度的row.names)
[1]  0 13

没有重复，直接进行注释

> biocLite("hgu133plus2.db")
> library(hgu133plus2.db)
> library(annotate)

# 利用getSYMBOL函数进行注释，这里只选择后几列额，方便观察
> exprmt$symbol <- getSYMBOL(exprmt$ID, 'hgu133plus2');exprmt <- exprmt[,8:14]; head(exprmt)
          GSM524669 GSM524670 GSM524671 GSM524672 GSM524673        ID symbol
1007_s_at  2.456892  2.404267  3.568920  3.307409  3.225964 1007_s_at   <NA>
1053_at   10.447158 10.799832  9.765222 10.069229 10.089496   1053_at   RFC2
117_at     4.500872  2.702614  3.015139  2.768683  2.853809    117_at  HSPA6
121_at     3.057630  3.072568  3.057630  3.057630  3.057630    121_at   PAX8
1255_g_at  2.278998  2.278998  2.278998  2.278998  2.278998 1255_g_at GUCA1A
1294_at    4.978118  4.713441  4.303196  4.452989  4.766019   1294_at   <NA>

可以看到有一些芯片ID是没有注释上的，接下来对注释情况进行检查，主要包括三个方面：（1）未注释上--NA；（2）未注释上--空字符串；（3）注释上了，但有重复。

# 查看一共多少行
> dim(exprmt)
[1] 54675     7
# symbol非NA的行
> exprmt.nona <- exprmt[!is.na(exprmt$symbol),]; dim(exprmt.nona)
[1] 41941     7
> dim(exprmt[is.na(exprmt$symbol),])
[1] 12734     7
# symbol重复的行数
> symbol.dup <- exprmt.nona[duplicated(exprmt.nona$symbol),]; head(symbol.dup)
             GSM524669 GSM524670 GSM524671 GSM524672 GSM524673           ID   symbol
1552264_a_at  8.016109  7.672156  7.661139  7.835451  7.225329 1552264_a_at    MAPK1
1552272_a_at  2.279247  2.279247  2.279247  2.279247  2.279247 1552272_a_at    PRR22
1552275_s_at  4.862782  5.539661  4.672476  4.709687  3.931607 1552275_s_at      PXK
1552279_a_at  2.436048  2.731219  2.436048  2.436048  2.464805 1552279_a_at  SLC46A1
1552289_a_at  2.278998  2.278998  2.278998  2.278998  2.278998 1552289_a_at    CILP2
1552303_a_at  4.421032  4.601013  4.614977  4.628660  4.628942 1552303_a_at TMEM106A
> dim(symbol.dup)
[1] 21749     7
> exprmt.nona[exprmt.nona$symbol=='PXK',]
             GSM524669 GSM524670 GSM524671 GSM524672 GSM524673           ID symbol
1552274_at    5.290347  5.539950  5.081292  5.016670  5.032546   1552274_at    PXK
1552275_s_at  4.862782  5.539661  4.672476  4.709687  3.931607 1552275_s_at    PXK
225796_at     7.140168  7.106552  7.232615  6.650670  7.167075    225796_at    PXK
238223_at     2.281464  2.281464  2.394397  2.281464  2.281464    238223_at    PXK
# 最后看看有没有空字符串
> exprmt.nona[exprmt.nona$symbol=='',]
[1] GSM524669 GSM524670 GSM524671 GSM524672 GSM524673 ID        symbol   
<0 行> (或0-长度的row.names)

# 最后得到非重复的，无NA的数据
> exprmt.final <- exprmt.nona[!duplicated(exprmt.nona$symbol),]; dim(exprmt.final)
[1] 20192     7
> exprmt.final[exprmt.final$symbol=='PXK',]
           GSM524669 GSM524670 GSM524671 GSM524672 GSM524673         ID symbol
1552274_at  5.290347   5.53995  5.081292   5.01667  5.032546 1552274_at    PXK

从上面可以看到，我们的芯片数据共有54675个features，其中注释到gene symbol的为41941个，在这些注释到的features中，有21749个是重复的。无空字符串gene symbol。最后得到非重复数据20192个（重复的则取其第一个出现的）。

下面我们下载芯片平台对应信息文件手动进行注释：

# 下载芯片对应信息
> gpl570 <- getGEO('GPL570', destdir=".")
> names(Meta(gpl570))
 [1] "contact_address"         "contact_city"            "contact_country"         "contact_email"           "contact_institute"       "contact_name"           
 [7] "contact_phone"           "contact_state"           "contact_web_link"        "contact_zip/postal_code" "data_row_count"          "description"            
[13] "distribution"            "geo_accession"           "last_update_date"        "manufacture_protocol"    "manufacturer"            "organism"               
[19] "relation"                "sample_id"               "series_id"               "status"                  "submission_date"         "taxid"                  
[25] "technology"              "title"                   "web_link"     
> colnames(Table(gpl570))
 [1] "ID"                               "GB_ACC"                           "SPOT_ID"                          "Species Scientific Name"         
 [5] "Annotation Date"                  "Sequence Type"                    "Sequence Source"                  "Target Description"              
 [9] "Representative Public ID"         "Gene Title"                       "Gene Symbol"                      "ENTREZ_GENE_ID"                  
[13] "RefSeq Transcript ID"             "Gene Ontology Biological Process" "Gene Ontology Cellular Component" "Gene Ontology Molecular Function"

# 选取其中的几列用来转换
> iddata <- Table(gpl570)[c('ID','Gene Symbol')]
> head(iddata)
         ID      Gene Symbol
1 1007_s_at DDR1 /// MIR4640
2   1053_at             RFC2
3    117_at            HSPA6
4    121_at             PAX8
5 1255_g_at           GUCA1A
6   1294_at MIR5193 /// UBA7

我们选取gpl信息中的ID、Gene Symbol用于转换。

> mann.anno <- merge(x=exprmt, y=iddata, by='ID',all.x=T,all.y=F); head(mann.anno); dim(mann.anno)
         ID GSM524669 GSM524670 GSM524671 GSM524672 GSM524673 symbol.RP      Gene Symbol
1 1007_s_at  2.456892  2.404267  3.568920  3.307409  3.225964      <NA> DDR1 /// MIR4640
2   1053_at 10.447158 10.799832  9.765222 10.069229 10.089496      RFC2             RFC2
3    117_at  4.500872  2.702614  3.015139  2.768683  2.853809     HSPA6            HSPA6
4    121_at  3.057630  3.072568  3.057630  3.057630  3.057630      PAX8             PAX8
5 1255_g_at  2.278998  2.278998  2.278998  2.278998  2.278998    GUCA1A           GUCA1A
6   1294_at  4.978118  4.713441  4.303196  4.452989  4.766019      <NA> MIR5193 /// UBA7

可以看到用R包注释的（symbol.RP列）和直接下载平台信息注释的有些许差别，R包未能注释到的，手动注释结果是对应多个信息。下面我们进一步查看一下其区别：

> library(dplyr)
> consistent <- filter(mann.anno, symbol.RP==`Gene Symbol`); head(consistent)
         ID GSM524669 GSM524670 GSM524671 GSM524672 GSM524673 symbol.RP Gene Symbol
1   1053_at 10.447158 10.799832  9.765222 10.069229 10.089496      RFC2        RFC2
2    117_at  4.500872  2.702614  3.015139  2.768683  2.853809     HSPA6       HSPA6
3    121_at  3.057630  3.072568  3.057630  3.057630  3.057630      PAX8        PAX8
4 1255_g_at  2.278998  2.278998  2.278998  2.278998  2.278998    GUCA1A      GUCA1A
5   1316_at  5.246463  5.547843  5.269052  4.794183  4.620719      THRA        THRA
6   1320_at  4.578203  4.446357  4.492585  4.502381  4.446357    PTPN21      PTPN21
> dim(mann.anno); dim(consistent)
[1] 54675     8
[1] 39145     8
# 可以看到一共有39145行注释结果都是相同的，那么这些中有多少是都未注释到的呢（即NA）
> consistent.na <- filter(consistent, is.na(symbol.RP));head(consistent.na)
[1] ID          GSM524669   GSM524670   GSM524671   GSM524672   GSM524673   symbol.RP   Gene Symbol
<0 行> (或0-长度的row.names)
# 结果是这些相同的中没有NA
# 再来看看注释到重复gene的情况吧
> dup <- consistent[duplicated(consistent$symbol.RP),]
> head(dup);dim(dup)
             ID GSM524669 GSM524670 GSM524671 GSM524672 GSM524673 symbol.RP Gene Symbol
16 1552264_a_at  8.016109  7.672156  7.661139  7.835451  7.225329     MAPK1       MAPK1
20 1552272_a_at  2.279247  2.279247  2.279247  2.279247  2.279247     PRR22       PRR22
22 1552275_s_at  4.862782  5.539661  4.672476  4.709687  3.931607       PXK         PXK
26 1552279_a_at  2.436048  2.731219  2.436048  2.436048  2.464805   SLC46A1     SLC46A1
32 1552289_a_at  2.278998  2.278998  2.278998  2.278998  2.278998     CILP2       CILP2
40 1552303_a_at  4.421032  4.601013  4.614977  4.628660  4.628942  TMEM106A    TMEM106A
[1] 20267     8
# 因此最后非重复的
> dim(consistent[!duplicated(consistent$symbol.RP),])
[1] 18878     8

对两种注释结果相同的 gene symbol 进行分析，一共得到了39145个注释上的gene，这其中有20267个是重复的，因此最后得到了18878个非重复注释的gene。
我们最后再对手动注释的结果进行分析一下。

# 看看为注释到的情况（NA）
> mann.anno[is.na(mann.anno$`Gene Symbol`),]
[1] ID          GSM524669   GSM524670   GSM524671   GSM524672   GSM524673   symbol.RP   Gene Symbol
<0 行> (或0-长度的row.names)
# 没有注释不上的，看来R包注释不上的这里都注释到多个gene symbol了（如DDR1 /// MIR4640）
> library(stringr)
> multi.anno <- mann.anno[str_detect(mann.anno$`Gene Symbol`, '///'),]; head(multi.anno); dim(multi.anno)
              ID GSM524669 GSM524670 GSM524671 GSM524672 GSM524673 symbol.RP                Gene Symbol
1      1007_s_at  2.456892  2.404267  3.568920  3.307409  3.225964      <NA>           DDR1 /// MIR4640
6        1294_at  4.978118  4.713441  4.303196  4.452989  4.766019      <NA>           MIR5193 /// UBA7
16    1552258_at  2.437344  2.458223  2.437344  2.437344  2.429044     CYTOR LINC00152 /// LOC101930489
32  1552283_s_at  2.278998  2.278998  2.278998  2.281762  2.278998      <NA>       ZDHHC11 /// ZDHHC11B
97    1552388_at  2.550713  2.389557  2.389557  2.389557  2.389557      <NA>        FLJ30901 /// SCUBE1
108 1552401_a_at  2.278998  2.278998  2.278998  2.294001  2.278998      <NA>        BACH1 /// GRIK1-AS2
[1] 2796    8
# 这里2796个features注释到多gene symbol，但意外的是，有的多gene symbol位点R包也注释到了
# 看看dup情况怎么样
> mann.anno.rmdup <- mann.anno[!duplicated(mann.anno$`Gene Symbol`),]; head(mann.anno.rmdup); dim(mann.anno.rmdup)
         ID GSM524669 GSM524670 GSM524671 GSM524672 GSM524673 symbol.RP      Gene Symbol
1 1007_s_at  2.456892  2.404267  3.568920  3.307409  3.225964      <NA> DDR1 /// MIR4640
2   1053_at 10.447158 10.799832  9.765222 10.069229 10.089496      RFC2             RFC2
3    117_at  4.500872  2.702614  3.015139  2.768683  2.853809     HSPA6            HSPA6
4    121_at  3.057630  3.072568  3.057630  3.057630  3.057630      PAX8             PAX8
5 1255_g_at  2.278998  2.278998  2.278998  2.278998  2.278998    GUCA1A           GUCA1A
6   1294_at  4.978118  4.713441  4.303196  4.452989  4.766019      <NA> MIR5193 /// UBA7
[1] 23521     8

从上面可以看到所有芯片ID都注释到了，没有NA，其中没有重复的有23521个（其中2796个是多gene symbol，因此单个gene symbol 20725个），较R包注释到的（20276）要多一些。

ID转换

接下来我们将上面注释到的Gene symbol转换为Entrez ID，使用的是HGNC提供的ID对应文件。

> hgnc <- read.csv('D:/TCGA/hgnc_complete_set.txt',header = T,stringsAsFactors = F,sep='\t');colnames(hgnc)
 [1] "hgnc_id"                  "symbol"                   "name"                     "locus_group"              "locus_type"               "status"                  
 [7] "location"                 "location_sortable"        "alias_symbol"             "alias_name"               "prev_symbol"              "prev_name"               
[13] "gene_family"              "gene_family_id"           "date_approved_reserved"   "date_symbol_changed"      "date_name_changed"        "date_modified"           
[19] "entrez_id"                "ensembl_gene_id"          "vega_id"                  "ucsc_id"                  "ena"                      "refseq_accession"        
[25] "ccds_id"                  "uniprot_ids"              "pubmed_id"                "mgd_id"                   "rgd_id"                   "lsdb"                    
[31] "cosmic"                   "omim_id"                  "mirbase"                  "homeodb"                  "snornabase"               "bioparadigms_slc"        
[37] "orphanet"                 "pseudogene.org"           "horde_id"                 "merops"                   "imgt"                     "iuphar"                  
[43] "kznf_gene_catalog"        "mamit.trnadb"             "cd"                       "lncrnadb"                 "enzyme_id"                "intermediate_filament_db"
[49] "rna_central_ids"     

> destid <- hgnc[c('symbol','entrez_id','ensembl_gene_id')];head(destid)
    symbol entrez_id ensembl_gene_id
1     A1BG         1 ENSG00000121410
2 A1BG-AS1    503538 ENSG00000268895
3     A1CF     29974 ENSG00000148584
4      A2M         2 ENSG00000175899
5  A2M-AS1    144571 ENSG00000245105
6    A2ML1    144568 ENSG00000166535

# ID 转换
> trans <- merge(x=exprmt[!duplicated(exprmt$symbol.RP),], y=destid, by.x='symbol.RP',by.y='symbol');cat('exprmt.rmdup\n');dim(exprmt[!duplicated(exprmt$symbol.RP),]);cat('destid\n');dim(destid);cat('trans\n');dim(trans)
exprmt.rmdup
[1] 20193     7
destid
[1] 42508     3
trans
[1] 19062     9

> trans.rmdup <- trans[!duplicated(trans$symbol.RP),]; dim(trans.rmdup)
[1] 19061     9

> trans.rmdup <- trans[!duplicated(trans$entrez_id),]; dim(trans.rmdup)
[1] 19060     9

可以看到，使用HGNC提供的ID转换文件，我们共得到了19062个无重复ID，比上面注释到的少了快2000个，对于这些差异基因有机会在探究一下。

最后编辑于：2018.05.12 15:33:00

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 204,684评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 87,143评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 151,214评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,788评论 1赞 277
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,796评论 5赞 368
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,665评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,027评论 3赞 399
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,679评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 41,346评论 1赞 299
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,664评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,766评论 1赞 331
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,412评论 4赞 321
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,015评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,974评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,203评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,073评论 2赞 350
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,501评论 2赞 343