(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”))

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-06

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏杨建荣的学习笔记

闪回数据库不是“万金油”(r11笔记第73天)

闪回数据库这个特性在很多Oracle DBA眼里就是鸡肋特性,因为谁会因为恢复数据而需要在主库闪回,最后可能丢掉更多的数据,这个观点没错。 但是...

3256
来自专栏企鹅号快讯

一文教会你数据库性能调优

前言 微软工程师的一个工程师曾经对性能调优有一个非常形象的比喻:剥洋葱 。我也非常认可,让我们来一层一层拨开外面它神秘的面纱。 ? 六大因素 下面祭出的是我们在...

1819
来自专栏数据和云

藏在表分区统计信息背后的小秘密

作者介绍 ? 曾令军 云和恩墨技术专家,8年数据库运维经验。思维敏捷,擅长于数据库开发、解决棘手的数据库故障和性能问题,在数据库故障诊断、运维监控、性能优化方面...

2695
来自专栏牛客网

链家,阿里面经链家:阿里:

今天下午面的北京链家现场面,虽然凉凉还是总结下面经吧~ 链家: 一面: 拿出手机问我笔试做错的一道笔试怎么分析,提醒了半天我也没想到(实际是拆装箱相关的知识)...

7369
来自专栏数据和云

经典文档:Oracle Database 12.2新特性概览解读下载

在2017 OOW大会上,关于Oracle Database 12.2 数据库的新特性介绍仍然引人瞩目,会后公布了 Oracle VP Swonger的文档,我...

3117
来自专栏程序员的SOD蜜

PDF.NET 数据开发框架 许可限制 框架源码的获取

欢迎使用 PDF.NET 数据开发框架 (Ver 4.0) 关于框架的名字由来          在我设计www.pwmis.cn 站点(原域名已经过期,现在正...

2016
来自专栏腾讯Bugly的专栏

【Dev Club 分享】微信 iOS SQLite 源码优化实践

Dev Club 是一个交流移动开发技术,结交朋友,扩展人脉的社群,成员都是经过审核的移动开发工程师。每周都会举行嘉宾分享,话题讨论等活动。 本期,我们邀请了腾...

3518
来自专栏搜云库

手把手教你 MongoDB 的安装与详细使用(一)

MongoDB 是由C++语言编写的,是一个基于分布式文件存储的开源数据库系统。

4867
来自专栏Python中文社区

使用Python自动生成报表以邮件发送

数据分析师肯定每天都被各种各样的数据数据报表搞得焦头烂额,老板的,运营的、产品的等等。而且大部分报表都是重复性的工作,这篇文章就是帮助大家如何用Python来实...

3515
来自专栏架构之美

MongoDB在58同城的应用实践

1413

扫码关注云+社区