GEO数据库之芯片数据下载

风雨都停了 2018-06-17

展开全文

GEO数据库是NCBI开发的基因表达数据库，主要接收通过高通量测序、基因芯片等方法获得的基因表达数据——这就方便大家利用他人数据发文章了。

而进行GEO数据挖掘的第一步就是进行数据下载，但是进入网站点点来下载，查询搜索工作就不少，下载下来的数据还不一定能看懂。有没有什么方法可以解决这个问题呢？当然有——R包GEOquery！下面就针对芯片数据，教大家用GEOquery包完成下载工作。

GEO数据

在下载之前要先了解GEO数据库具体存放的四类数据：GSE、GDS、GSM、和GPL。

一个GSE Accession对应的是整个研究项目的系列的数据，可能涉及不同平台；

一个GDS Accession对应的一个同一平台的数据集；

一个GSM Accession对应单一样品的数据信息，它只能是单一平台的数据，往往，GSE 和GDS中会包含多个GSM数据；

一个GPL Accession，则对应一个platform信息。

R包安装与加载

GEOquery

## try http:// if https:// URLs are not supported
source('https:///biocLite.R')
biocLite('GEOquery')

Biobase

## try http:// if https:// URLs are not supported
source('https:///biocLite.R')
biocLite('Biobase')

正确加载

library('Biobase')
library('GEOquery')
setwd('F:/GEO') ############有需要可以设置路径

利用GSE Accession

通过阅读文献查找感兴趣的GSE Accession，下载对应的表达数据和平台信息等，可以利用GEOquery中的getGEO()函数下载series_matrix.txt。例如GSE70213：

> gse = getGEO('GSE70213', GSEMatrix =TRUE, destdir = '.', getGPL = T, AnnotGPL = T)
###destdir设置当前目录，getGPL 和AnnotGPL都设置TRUE，可以下载和获得平台的注释文件

gse为列表数据，对应的GSM是单平台，则length为1，之后分别利用Biobase包中的exprs()、pData()和fData()获得表达数据、样品处理分组等信息、芯片平台的设计注释信息，还可以利用annotation()函数了解对应GPL Accession，譬如exprs()函数：

> exprSet=exprs(gse[[1]])
> head(exprSet,2)
GSM1720833 GSM1720834 GSM1720835 GSM1720836 GSM1720837 GSM1720838 GSM1720839 GSM1720840 GSM1720841 GSM1720842
10338001 2041.40800 2200.86100 2323.7600 3216.26300 2362.77500 2195.31800 2013.35900 2146.25800 1785.9460 2067.04100
10338002 63.78059 65.08438 58.3082 75.86145 66.95605 43.81526 49.11361 51.29279 48.9604 42.14286
GSM1720843 GSM1720844 GSM1720845 GSM1720846 GSM1720847 GSM1720848 GSM1720849 GSM1720850 GSM1720851 GSM1720852
10338001 1769.1150 1720.77400 1847.42900 2214.69800 2279.51500 2530.45600 2303.26400 2358.83400 1701.40000 1970.92400
10338002 42.5472 43.48373 64.34628 59.75188 57.48852 60.26423 54.81179 53.70885 57.86877 57.02808
GSM1720853 GSM1720854 GSM1720855 GSM1720856
10338001 1822.78600 2014.26000 1737.84200 2001.73400
10338002 59.26121 55.27306 54.36722 49.43959

注释信息的获取可以进行探针和基因的对应，方便后续分析。通过exprs()、pData()和fData()获得的数据都可以利用write.table等进行文件保存。

利用GDS Accession

GDS数据同样可以利用getGEO()函数下载soft文件。例如GDS5881：

> gds = getGEO('GDS5881', GSEMatrix =TRUE, destdir = '.', getGPL = T, AnnotGPL = T)
###destdir设置当前目录

gds可以利用GEOquery包中的Table()获取表达数据，并利用Meta()获得描述信息，其中Meta(gds)$platform可以获得GPL Accession。

> exprSet=Table(gds)
> head(exprSet,1)
ID_REF IDENTIFIER GSM1720845 GSM1720846 GSM1720847 GSM1720848 GSM1720849 GSM1720850 GSM1720851 GSM1720852 GSM1720853
1 10344614 Gm2889 48.4971 47.252 39.3331 49.9048 36.8313 41.9501 37.5569 38.1924 46.0668
GSM1720854 GSM1720855 GSM1720856
1 34.689 38.5762 32.2618

> Meta(gset)$platform
[1] 'GPL6246'

针对getGEO返回的gds——GDS数据，可以利用GEOquery包中GDS2Set()和GDS2MA()转变为为ExpressionSets 和limma MALists。

> gds2eSet=GDS2eSet(gds)
> MA=GDS2MA(gds)

再针对返回的gds2eSet,利用exprs()、pData()和fData()同样可以获得表达数据、样品处理分组信息、芯片平台的设计注释信息。返回的MA中涉及大量的描述信息，其中MA$tragets也是样品处理信息。

利用GSM Accession

利用GSM Accession下载的是单样本的表达数据，例如GSM1720833：

> gsm = getGEO('GSM1720833', GSEMatrix =TRUE, destdir = '.', getGPL = T, AnnotGPL = T)

针对gsm，同样是利用GEOquery包中的Table()获取表达数据，并利用Meta()获得描述信息，而获取对应的GSE Accession 和GPL Accsesion利用Meta(gsm)$series_id和Meta(gsm)$platform_id。

利用GPL Accession

针对芯片平台，利用GPL Accession下载得到的数据是芯片的设计和注释信息，可以获得探针组和基因的对应关系，利用Table()函数可以显示annotation，例如GPL6246：

> gpl = getGEO('GPL6246', GSEMatrix =TRUE, destdir = '.', getGPL = T, AnnotGPL = T)
> ann=Table(gpl)
> head(ann,2)
       ID          Gene title Gene symbol   Gene ID UniGene title UniGene symbol UniGene ID
1 10344614 predicted gene 2889      Gm2889 100040658                                        
2 10344616                                                                                  
                                                                                                                                                                                                                                                                                                                                                                                                                                                            Nucleotide Title
1 Mus musculus blastocyst blastocyst cDNA, RIKEN full-length enriched library, clone:I1C0009C06 product:hypothetical DeoxyUTP pyrophosphatase/Aspartyl protease, retroviral-type family profile/Retrovirus capsid, C-terminal/Peptidase aspartic/Peptidase aspartic, active site containing protein, full insert sequence///Mus musculus blastocyst blastocyst cDNA, RIKEN full-length enriched library, clone:I1C0042P10 product:hypothetical protein, full insert sequence
2                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
                   GI   GenBank Accession Platform_CLONEID Platform_ORF      Platform_SPOTID Chromosome location
1 74211482///74217103 AK145513///AK145782                               chr1:3054233-3054733                  18
2                                                                       chr1:3102016-3102125                    
  Chromosome annotation GO:Function GO:Process GO:Component GO:Function ID GO:Process ID GO:Component ID
1         Chromosome 18                                                                                 
2