(16)芯片探针与基因的对应关系-生信菜鸟团博客2周年精选文章集

这个我非常喜欢，目录如下：

用R获取芯片探针与基因的对应关系三部曲-bioconductor 用R获取芯片探针与基因的对应关系三部曲-NCBI下载对应关系 gene的各种ID转换终结者-bioconductor系列包

现有的基因芯片种类不要太多了！

但是重要而且常用的芯片并不多！

一般分析芯片数据都需要把探针的ID切换成基因的ID，我一般喜欢用基因的entrez ID。

一般有三种方法可以得到芯片探针与gene的对应关系。

金标准当然是去基因芯片的厂商的官网直接去下载啦！！！

一种是直接用bioconductor的包

一种是从NCBI里面下载文件来解析好！

首先，我们说官网，肯定可以找到，不然这种芯片出来就没有意义了！

然后，我们看看NCBI下载的，会比较大

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6947

这两种方法都比较麻烦，需要一个个的来！

所以我接下来要讲的是用R的bioconductor包来批量得到芯片探针与gene的对应关系！

一般重要的芯片在R的bioconductor里面都是有包的，用一个R包可以批量获取有注释信息的芯片平台，我选取了常见的物种，如下：

gpl organism bioc_package 1 GPL32 Mus musculus mgu74a 2 GPL33 Mus musculus mgu74b 3 GPL34 Mus musculus mgu74c 6 GPL74 Homo sapiens hcg110 7 GPL75 Mus musculus mu11ksuba 8 GPL76 Mus musculus mu11ksubb 9 GPL77 Mus musculus mu19ksuba 10 GPL78 Mus musculus mu19ksubb 11 GPL79 Mus musculus mu19ksubc 12 GPL80 Homo sapiens hu6800 13 GPL81 Mus musculus mgu74av2 14 GPL82 Mus musculus mgu74bv2 15 GPL83 Mus musculus mgu74cv2 16 GPL85 Rattus norvegicus rgu34a 17 GPL86 Rattus norvegicus rgu34b 18 GPL87 Rattus norvegicus rgu34c 19 GPL88 Rattus norvegicus rnu34 20 GPL89 Rattus norvegicus rtu34 22 GPL91 Homo sapiens hgu95av2 23 GPL92 Homo sapiens hgu95b 24 GPL93 Homo sapiens hgu95c 25 GPL94 Homo sapiens hgu95d 26 GPL95 Homo sapiens hgu95e 27 GPL96 Homo sapiens hgu133a 28 GPL97 Homo sapiens hgu133b 29 GPL98 Homo sapiens hu35ksuba 30 GPL99 Homo sapiens hu35ksubb 31 GPL100 Homo sapiens hu35ksubc 32 GPL101 Homo sapiens hu35ksubd 36 GPL201 Homo sapiens hgfocus 37 GPL339 Mus musculus moe430a 38 GPL340 Mus musculus mouse4302 39 GPL341 Rattus norvegicus rae230a 40 GPL342 Rattus norvegicus rae230b 41 GPL570 Homo sapiens hgu133plus2 42 GPL571 Homo sapiens hgu133a2 43 GPL886 Homo sapiens hgug4111a 44 GPL887 Homo sapiens hgug4110b 45 GPL1261 Mus musculus mouse430a2 49 GPL1352 Homo sapiens u133x3p 50 GPL1355 Rattus norvegicus rat2302 51 GPL1708 Homo sapiens hgug4112a 54 GPL2891 Homo sapiens h20kcod 55 GPL2898 Rattus norvegicus adme16cod 60 GPL3921 Homo sapiens hthgu133a 63 GPL4191 Homo sapiens h10kcod 64 GPL5689 Homo sapiens hgug4100a 65 GPL6097 Homo sapiens illuminaHumanv1 66 GPL6102 Homo sapiens illuminaHumanv2 67 GPL6244 Homo sapiens hugene10sttranscriptcluster 68 GPL6947 Homo sapiens illuminaHumanv3 69 GPL8300 Homo sapiens hgu95av2 70 GPL8490 Homo sapiens IlluminaHumanMethylation27k 71 GPL10558 Homo sapiens illuminaHumanv4 72 GPL11532 Homo sapiens hugene11sttranscriptcluster 73 GPL13497 Homo sapiens HsAgilentDesign026652 74 GPL13534 Homo sapiens IlluminaHumanMethylation450k 75 GPL13667 Homo sapiens hgu219 76 GPL15380 Homo sapiens GGHumanMethCancerPanelv1 77 GPL15396 Homo sapiens hthgu133b 78 GPL17897 Homo sapiens hthgu133a

这些包首先需要都下载

gpl_info=read.csv(“GPL_info.csv”,stringsAsFactors = F) ### first download all of the annotation packages from bioconductor for (i in 1:nrow(gpl_info)){ print(i) platform=gpl_info[i,4] platform=gsub(‘^ ‘,””,platform) ##主要是因为我处理包的字符串前面有空格 #platformDB=’hgu95av2.db’ platformDB=paste(platform,”.db”,sep=””) if( platformDB %in% rownames(installed.packages()) == FALSE) { BiocInstaller::biocLite(platformDB) #source(“http://bioconductor.org/biocLite.R“); #biocLite(platformDB ) } }

下载完了所有的包， 就可以进行批量导出芯片探针与gene的对应关系！

for (i in 1:nrow(gpl_info)){ print(i) platform=gpl_info[i,4] platform=gsub(‘^ ‘,””,platform) #platformDB=’hgu95av2.db’ platformDB=paste(platform,”.db”,sep=””) if( platformDB %in% rownames(installed.packages()) != FALSE) { library(platformDB,character.only = T) #tmp=paste(‘head(mappedkeys(‘,platform,’ENTREZID))’,sep=”) #eval(parse(text = tmp)) ###重点在这里，把字符串当做命令运行 all_probe=eval(parse(text = paste(‘mappedkeys(‘,platform,’ENTREZID)’,sep=”))) EGID <- as.numeric(lookUp(all_probe, platformDB, “ENTREZID”)) ##自己把内容写出来即可 } }

参考：http://blog.sina.com.cn/s/blog_62b37bfe0101jbuq.html

这是系列文章，请先看：

用R获取芯片探针与基因的对应关系三部曲-bioconductor

ncbi现有的GPL已经过万了，但是bioconductor的芯片注释包不到一千，虽然bioconductor可以解决我们大部分的需要，比如affymetrix的95,133系列，深圳1.0st系列，HTA2.0系列，但是如果碰到比较生僻的芯片，bioconductor也不会刻意为之制作一个bioconductor的包，这时候就需要自行下载NCBI的GPL信息了，也可以通过R来解决：

##本质上是下载一个文件，读进R里面，然后解析行列式，得到芯片探针与基因的对应关系，看下面的代码，你就能理解了。

## A-AGIL-28 – Agilent Whole Human Genome Microarray 4x44K 014850 G4112F (85 cols x 532 rows) library(Biobase) library(GEOquery) #Download GPL file, put it in the current directory, and load it: gpl <- getGEO(‘GPL6480′, destdir=”.”) colnames(Table(gpl)) ## [1] 41108 17 head(Table(gpl)[,c(1,6,7)]) ## you need to check this , which column do you need write.csv(Table(gpl)[,c(1,6,7)],”GPL6400.csv”) #platformDB=’hgu133plus2.db’ #library(platformDB, character.only=TRUE) probeset <- featureNames(GSE32575[[1]]) library(Biobase) library(GEOquery) #Download GPL file, put it in the current directory, and load it: gpl <- getGEO(‘GPL6102′, destdir=”.”) colnames(Table(gpl)) ## [1] 41108 17 head(Table(gpl)[,c(1,10,13)]) ## you need to check this , which column do you need probe2symbol=Table(gpl)[,c(1,13)] ## GPL15207 [PrimeView] Affymetrix Human Gene Expression Array probeset <- featureNames(GSE58979[[1]]) library(Biobase) library(GEOquery) #Download GPL file, put it in the current directory, and load it: gpl <- getGEO(‘GPL15207′, destdir=”.”) colnames(Table(gpl)) ## [1] 49395 24 head(Table(gpl)[,c(1,15,19)]) ## you need to check this , which column do you need probe2symbol=Table(gpl)[,c(1,15)]

## GPL10558 Illumina HumanHT-12 V4.0 expression beadchip library(Biobase) library(GEOquery) #Download GPL file, put it in the current directory, and load it: gpl <- getGEO(‘GPL10558′, destdir=”.”) colnames(Table(gpl)) ## [1] 41108 17 head(Table(gpl)[,c(1,10,13)]) ## you need to check this , which column do you need probe2symbol=Table(gpl)[,c(1,13)]

经常会有人问这样的问题I have list of 10,000 Entrez IDs and i want to convert the multiple Entrez IDs into the respective gene names. Could someone suggest me the way to do this?等等类似的基因转换，能做的基因转换的方法非常多，以前不懂编程的时候，都是用各种网站，而最常用的就是ensembl的biomart了，它支持的ID非常多，高达几百种，好多ID我到现在都不知道是什么意思。

现在学会编程了，我比较喜欢的是R的一些包，是bioconductor系列，一般来说，其中有biomart，org.Hs.eg.db，annotate，等等。关于biomart我就不再讲了，我前面的博客至少有七八篇都提到了它。本次我们讲讲简单的， 我就以把gene entrez ID转换为gene symbol 为例子把。

当然，首先要安装这些包，并且加载。

if(“org.Hs.eg.db” %in% rownames(installed.packages()) == FALSE) {source(“http://bioconductor.org/biocLite.R”);biocLite(“org.Hs.eg.db”)} suppressMessages(library(org.Hs.eg.db)) 我比较喜欢这样加载包

library(annotate) #一般都是这样加载包

如果是用org.Hs.eg.db包，首先你只需要读取你的待转换ID文件，构造成一个向量，tmp，然后只需要symbols <- org.Hs.egSYMBOL[as.character(tmp)]就可以得到结果了，返回的symbols是一个对象，需要用toTable这个函数变成数据框。但是这样转换容易出一些问题，比如如果你的输入数据tmp，里面含有一些无法转换的gene entrez ID，就会报错。

而且它支持的ID转换很有限，具体看看它的说明书即可：https://www.bioconductor.org/packages/release/data/annotation/manuals/org.Hs.eg.db/man/org.Hs.eg.db.pdf

org.Hs.eg.db org.Hs.eg_dbconn org.Hs.egACCNUM org.Hs.egALIAS2EG org.Hs.egCHR org.Hs.egCHRLENGTHS org.Hs.egCHRLOC org.Hs.egENSEMBL org.Hs.egENSEMBLPROT org.Hs.egENSEMBLTRANS org.Hs.egENZYME org.Hs.egGENENAME org.Hs.egGO org.Hs.egMAP org.Hs.egMAPCOUNTS org.Hs.egOMIM org.Hs.egORGANISM org.Hs.egPATH org.Hs.egPMID org.Hs.egREFSEQ org.Hs.egSYMBOL org.Hs.egUCSCKG org.Hs.egUNIGENE org.Hs.egUNIPROT

如果是用annotate包，首先你还是需要读取你的待转换ID文件，构造成一个向量，tmp，然后用getSYMBOL(as.character(tmp), data=’org.Hs.eg’)这样直接就返回的还是以向量，只是在原来向量的基础上面加上了names属性。说明书：http://www.bioconductor.org/packages/3.3/bioc/manuals/annotate/man/annotate.pdf

然后你可以把转换好的向量写出去，如下：

1 A1BG 2 A2M 3 A2MP1 9 NAT1 10 NAT2 12 SERPINA3 13 AADAC 14 AAMP 15 AANAT 16 AARS

PS：如果是芯片数据，需要把探针的ID转换成gene，那么一般还需要加载特定芯片的数据包才行：

platformDB <- paste(eset.mas5@annotation, “.db”, sep=””) #这里需要确定你用的是什么芯片 cat(“the annotation is “,platformDB,”\n”) if(platformDB %in% rownames(installed.packages()) == FALSE) {source(“http://bioconductor.org/biocLite.R”);tmp=try(biocLite(platformDB))} library(platformDB, character.only=TRUE) probeset <- featureNames(eset.mas5) rowMeans <- rowMeans(exprSet)

library(annotate) # lookUp函数是属于annotate这个包的 EGID <- as.numeric(lookUp(probeset, platformDB, “ENTREZID”))