解读GEO数据存放规律及下载,一文就够
做生物信息学分析的朋友应该是对GEO数据库耳熟能详了,总会用到公共数据库的,而GEO数据库则是首选,起先只是为表达芯片数据准备的,后期纳入了各种NGS组学数据,文章里面会给出数据地址,GSE ID号,由此我们就可以进入GEO数据库,进而了解它!
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84498
GEO数据库基本介绍
其实只需要理解下面的4个概念。
GEO Platform (GPL)
GEO Sample (GSM)
GEO Series (GSE)
GEO Dataset (GDS)
理解起来也很容易。一篇文章可以有一个或者多个GSE数据集,一个GSE里面可以有一个或者多个GSM样本。多个研究的GSM样本可以根据研究目的整合为一个GDS,不过GDS本身用的很少。而每个数据集都有着自己对应的芯片平台,就是GPL。
使用GEOquery包从GEO数据库下载数据
官网详细用法:https://warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/geo/ 而且我写了这个软件包的中文用法:http://www.bio-info-trainee.com/bioconductor_China/software/GEOquery.html
而且需要理解所有bioconductor支持的芯片平台对应关系:通过bioconductor包来获取所有的芯片探针与gene的对应关系
首先是下载函数的使用
包的本质就是getGEO
函数,用法列举如下:
1、根据GDS号来下载数据,下载soft文件
gds858 <- getGEO(‘GDS858’, destdir=“.”)
2、根据GPL号下载的是芯片设计的信息!
gpl96 <- getGEO(‘GPL96’, destdir=“.”) ```
3、根据GSE号下载数据,下载_series_matrix.txt.gz
gse1009 <- getGEO(‘GSE1009’, destdir=“.”)
上面的代码下载的文件都会保存在本地,destdir参数指定下载地址。比较重要的三个参数是:GSEMatrix=TRUE,AnnotGPL=FALSE,getGPL=TRUE
然后是了解下载函数返回的对象
1、查看下载GDS后返回的对象
用Table(gds858)
可以得到表达矩阵!
用Meta(gds858)
可以得到描述信息!
具体代码如下:
options(warn=-1)
suppressMessages(library(GEOquery))
gds858 <- getGEO('GDS858', destdir=".")
names(Meta(gds858))
Table(gds858)[1:5,1:5]
然后还可以用 GDS2eSet函数把它转变为 expression set 对象
, expression set这个对象
非常之重要,后续会详细讲解:
eset <- GDS2eSet(gds858, do.log2=TRUE)
2、GSE
直接根据GSE号返回的对象:gse1009,就是expression set 对象
我们的处理函数有:geneNames/sampleNames/pData/exprs
(这个是重点,对expression set 对象的操作函数)
3、GPL
根据GPL号下载返回的对象跟GDS一样,也是用Table/Meta处理!
options(warn=-1)
suppressMessages(library(GEOquery))
gpl96 <- getGEO('GPL96', destdir=".")
names(Meta(gpl96))
Table(gpl96)[1:10,1:4]
一般来说我们下载GPL是为了得到芯片的探针对应基因ID的关系列表,下面这个代码就是芯片ID的基因注释信息
#Table(gpl96)[1:10,c("ID","GB_LIST","Gene.Title","Gene.Symbol","Entrez.Gene")]
Table(gpl96)[1:10,c("ID","Gene Title","Gene Symbol","ENTREZ_GENE_ID")]
包装成函数
downGSE <- function(studyID = "GSE1009", destdir = ".") {
library(GEOquery)
eSet <- getGEO(studyID, destdir = destdir, getGPL = F)
exprSet = exprs(eSet[[1]])
pdata = pData(eSet[[1]])
write.csv(exprSet, paste0(studyID, "_exprSet.csv"))
write.csv(pdata, paste0(studyID, "_metadata.csv"))
return(eSet)
}
参考链接:http://www.bio-info-trainee.com/1085.html
番外
也可以用GEOmetadb包来获取对应GEO数据的实验信息(得到metadata数据),可以批量得到多个GSE数据集的信息,如下:
GSE1009
GSE10785
GSE1133
GSE11975
GSE121
GSE12409
执行下面代码即可:
library(GEOmetadb)
if(!file.exists('GEOmetadb.sqlite')) getSQLiteFile()
## 取决于网速哦
file.info('/path/GEOmetadb.sqlite')
con <- dbConnect(SQLite(),'/path/GEOmetadb.sqlite'))
#dbListTables(con2)
#dbListFields(con2,'gse')
GeoList = read.table("diabetes.GEO.list")
query = paste("select + from gsm where series_id in
( ' ", gsub(", ", " ', ' ", paste(Geolist[,1], collapse=",")," ')", seq=" ")
query
tmp = dbGetQuery(con2, query)
write.csv(tmp, "diabetes.GEO.meta.csv")
作业
大家可以使用上面介绍的R包及GEO数据库知识,把文中列出的几个GSE数据集探索一下。
或者下面这些:
GSE11072 2009-gastric cancer SBC Human 16K cDNA Microarray
GSE42872 2015-melanoma-vemurafenib HuGene-1_0-st
GSE24673 2015-hub-gene-mcode-retinoblastoma HuGene-1_0-st
GSE22863 2011-NSCLC HuGene-1_0-st
GSE622221, GSE4180414, GSE5140122 A total of 117 samples (54 cases and 63 controls) Affymetrix Human Genome U133 Plus 2.0 Array 2015-HCC
GSE21815 2016-CRC Agilent-014850 Whole Human Genome Microarray 4x44K