GEO数据下载方法

港控/mmm° 2024-03-26 08:09 67阅读 0赞

### **方法一：直接从**[**官网**][Link 1]**下载，手动** ###

GEO（GENE EXPRESSION OMNIBUS），[https://www.ncbi.nlm.nih.gov/geo/][https_www.ncbi.nlm.nih.gov_geo]，由美国国立生物技术信息中心NCBI创建维护的，是个公开的基因数据库，包含了测序和芯片数据。在前面，我们介绍过利用GEO数据库进行芯片数据检索，今天我们再进一步细化，如何利用GEO数据库下载信息。GEO数据包括五种，platforms、samples、series、datasets和profiles。

其中原始数据包括GPL、GSM和GSE。

GPL（GEOPlatform）:平台信息，由芯片或测序公司提供，含有芯片或测序平台的描述信息，芯片还包含了其注释信息，每个平台列出了使用该平台的所有样本和系列。

GSM（GEO Sample）:样本信息，记录单个样本的生物学信息，处理流程及该样本的原始数据（芯片或测序）。注意的是每个样本数据仅对应一个平台。

GSE（GSESeries）：系列，将一项研究中所关联的GSM信息集合在一起，含研究名称、设计、概要信息。与GSM不同的是，GSE可以包含多个平台以及子系列，一个样本可以出现在多个系列中。

GDS（GEO DataSet）：经挑选整理的数据集记录，如进行背景校正、均一化处理。每个GDS对应一个平台。

GEO profiles：来源与GDS数据，可以展现单个基因表达水平。

GEO测序文件存储形式包含SOFT、MINiML、Series Matrix files以及Supplementary files。

SOFT和MINiML存储的内容相同，但格式不同。SOFT为ASCII格式，MINiML为XML格式。Series Matrix files则以制表符为分隔的包含每个样本具体数值的文本文件，包含GSM和GSE。Supplementary files列出GSM原始数据或一些样本临床信息。

接下来我们到了我们演示环节，首先进入GEO官网。官网首页含有GEO数据库的概要，一些工具菜单，GEO数据库概况以及对于上传者的操作指南。

![format_png][]

![format_png 1][]

我们这一步需要下载两个文件：

（1）表达矩阵“series matrix”文件和（2）GPL平台注释文件。

![a7bad55fb2111f52f32332a0c181249a.png][]

![f63b4b7e09936271a2c735f7e5efdb16.png][]

之后可以用R进行表达谱和平台数据的合并。

**注意！到这里我们获得的只是初步的表达数据，还没有经过预处理，需要用R处理多个探针对应一个表达值，无对应symbol，以及合并多个探针对应一个symbol的情况后才可进行后续分析。**

### **方法二：代码下载，自动化** ###

我们使用的是 GEOquery 包，下来就来详细介绍一下

#### 安装GEOquery ####

从 Bioconductor 安装

if (!requireNamespace("BiocManager", quietly=TRUE))
        install.packages("BiocManager")
    BiocManager::install("GEOquery")

从 GitHub 安装

library(devtools)
    install_github('GEOquery','seandavi')

#### **使用** ####

##### 1 GEOquery 入门 #####

从 GEO 中获取数据真的很容易。只需要一个命令: getGEO。这一个函数对其输入进行解释，以确定如何从 GEO 中获取数据，然后将数据解析成有用的 R 数据结构。

用法很简单。首先加载 GEOquery 库

library(GEOquery)

现在，我们可以自由访问任何 GEO accession number。请注意，在下文中，我使用了一个与 GEOquery 包一起打包的文件。

一般来说，你只要传入正确的 GEO accession 就行，正如代码注释中所指出的那样

# 更常用的方式是传入 GEO accession，从 GEO 数据库中下载数据:
    # gds <- getGEO("GDS507")
    gds <- getGEO(filename=system.file("extdata/GDS507.soft.gz",package="GEOquery"))

现在，gds 包含了 R 数据结构(GDS 类)，表示 GEO 的 GDS507 条目。

你会注意到用于存储下载的文件名被输出到屏幕上（但没有保存到任何地方），以便以后调用 getGEO(filename=...) 时使用。

# 更常用的方式是直接传入 GSM 号，从数据库中下载：
    # gds <- getGEO("GSM11805")
    gsm <- getGEO(filename=system.file("extdata/GSM11805.txt.gz",package="GEOquery"))

###### 1.1 getGEO ######

1.  描述

该函数是 GEOquery 包中主要的用户级函数。它引导下载（如果没有指定文件名）并将 GEO SOFT 格式文件解析为 R 数据结构，该结构是专门设计的，便于访问 GEO SOFT 格式的每个重要部分。

1.  使用

getGEO(GEO =NULL, filename =NULL, destdir = tempdir(),
           GSElimits =NULL, GSEMatrix =TRUE, 
           AnnotGPL =FALSE, getGPL =TRUE,
           parseCharacteristics =TRUE)

1.  参数

![format_png 2][]

##### 2 GEOquery 数据结构 #####

GEOquery 的数据结构其实有两种形式。第一种，由 GDS、GPL 和 GSM 组成，它们的行为都是相似的，访问器对每一种都有相似的作用。

第四种 GEOquery 数据结构 GSE，是由 GSM 和 GPL 对象组合而成的复合数据类型。我先把前三个一起解释一下。

###### 2.1 GDS、GSM 和 GPL 类 ######

这些类中的每一个类都由一个元数据头（与 SOFT 格式头部信息几乎一致）和一个 GEODataTable 组成。

GEODataTable 有两个简单的部分: Columns 和 Table，Columns 部分是对 Table 部分的列的描述。此外，每个类还有一个 show 方法。

例如，对于上面的 gsm，查看元数据

> head(Meta(gsm))
    $channel_count
    [1] "1"
    
    $comment
    [1] "Raw data provided as supplementary file"
    
    $contact_address
    [1] "715 Albany Street, E613B"
    
    $contact_city
    [1] "Boston"
    
    $contact_country
    [1] "USA"
    
    $contact_department
    [1] "Genetics and Genomics"

查看与 GSM 相关的数据

> Table(gsm)[1:5,]
              ID_REF  VALUE ABS_CALL
    1 AFFX-BioB-5_at  953.9        P
    2 AFFX-BioB-M_at 2982.8        P
    3 AFFX-BioB-3_at 1657.9        P
    4 AFFX-BioC-5_at 2652.7        P
    5 AFFX-BioC-3_at 2019.5        P

查看列描述

> Columns(gsm)
        Column                                                                Description
    1   ID_REF                                                                           
    2    VALUE                         MAS 5.0 Statistical Algorithm (mean scaled to 500)
    3 ABS_CALL MAS 5.0 Absent, Marginal, Present call  with Alpha1 =0.05, Alpha2 =0.065

GPL 类的行为与 GSM 类完全相同。然而，GDS 类的 Columns 方法保存了更多的信息。

> Columns(gds)[,1:3]
         sample disease.state individual
    1  GSM11815           RCC        035
    2  GSM11832           RCC        023
    3  GSM12069           RCC        001
    4  GSM12083           RCC        005
    5  GSM12101           RCC        011
    6  GSM12106           RCC        032
    7  GSM12274           RCC          2
    8  GSM12299           RCC          3
    9  GSM12412           RCC          4
    10 GSM11810        normal        035
    11 GSM11827        normal        023
    12 GSM12078        normal        001
    13 GSM12099        normal        005
    14 GSM12269        normal          1
    15 GSM12287        normal          2
    16 GSM12301        normal          3
    17 GSM12448        normal          4

###### 2.2 GSE 类 ######

GSE 是 GEO 中最容易混淆的。一个 GSE 条目可以代表在任意数量的平台上运行的任意数量的样本。

GSE 类和其他类一样，有一个元数据部分，但是它没有 GEODataTable。

相反，它包含了两个列表，可以使用 GPLList 和 GSMList 方法进行访问，这两个列表分别是 GPL 和 GSM 对象的列表。

举个例子，我们从本地读取 soft 文件

gse <- getGEO(filename=system.file("extdata/GSE781_family.soft.gz",package="GEOquery"))

查看元数据

> head(Meta(gse))
    $contact_address
    [1] "715 Albany Street, E613B"
    
    $contact_city
    [1] "Boston"
    
    $contact_country
    [1] "USA"
    
    $contact_department
    [1] "Genetics and Genomics"
    
    $contact_email
    [1] "mlenburg@bu.edu"
    
    $contact_fax
    [1] "617-414-1646"

获取 GSE 中包含的所有 GSM 对象的名称

> names(GSMList(gse))
     [1] "GSM11805" "GSM11810" "GSM11814" "GSM11815" "GSM11823" "GSM11827" "GSM11830" "GSM11832" "GSM12067"
    [10] "GSM12069" "GSM12075" "GSM12078" "GSM12079" "GSM12083" "GSM12098" "GSM12099" "GSM12100" "GSM12101"
    [19] "GSM12105" "GSM12106" "GSM12268" "GSM12269" "GSM12270" "GSM12274" "GSM12283" "GSM12287" "GSM12298"
    [28] "GSM12299" "GSM12300" "GSM12301" "GSM12399" "GSM12412" "GSM12444" "GSM12448"

获得列表中的第一个 GSM 对象

> GSMList(gse)[[1]]
    An object of class "GSM"
    channel_count 
    [1] "1"
    comment 
    [1] "Raw data provided as supplementary file"
    contact_address 
    [1] "715 Albany Street, E613B"
    contact_city 
    [1] "Boston"
    contact_country 
    [1] "USA"
    contact_department 
    [1] "Genetics and Genomics"
    contact_email 
    [1] "mlenburg@bu.edu"
    ...
    ****** Column Descriptions ******
        Column                                                                Description
    1   ID_REF                                                                           
    2    VALUE                         MAS 5.0 Statistical Algorithm (mean scaled to 500)
    3 ABS_CALL MAS 5.0 Absent, Marginal, Present call  with Alpha1 = 0.05, Alpha2 = 0.065
    ****** Data Table ******
              ID_REF  VALUE ABS_CALL
    1 AFFX-BioB-5_at  953.9        P
    2 AFFX-BioB-M_at 2982.8        P
    3 AFFX-BioB-3_at 1657.9        P
    4 AFFX-BioC-5_at 2652.7        P
    5 AFFX-BioC-3_at 2019.5        P
    22278 more rows ...

获取所有 GPL 名称

> names(GPLList(gse))
    [1] "GPL96" "GPL97"

##### 3 转换为表达谱 #####

GEO 数据集与 limma 的 MAList 数据结构和 Biobase 的 ExpressionSet 数据结构非常相似

存在两个函数 GDS2MA 和 GDS2eSet，能够将 GEOquery 数据结构分别转换为 MAList 和 ExpressionSet

###### 3.1 获取 GSE Series Matrix ######

除了解析较大的 soft 文件，我们也可以直接获取处理后的 Series matrix 文件。

getGEO 能够快速解析该类型的文件，解析返回的数据结构是 ExpressionSet 的列表

gset <- getGEO("GSE11675",destdir ='./')

> show(gset)
    $GSE11675_series_matrix.txt.gz
    ExpressionSet (storageMode: lockedEnvironment)
    assayData:12625 features,6 samples 
      element names: exprs 
    protocolData: none
    phenoData
      sampleNames: GSM296630 GSM296635 ... GSM296639 (6 total)
      varLabels: title geo_accession ... data_row_count (34 total)
      varMetadata: labelDescription
    featureData
      featureNames:1000_at 1001_at ... AFFX-YEL024w/RIP1_at (12625 total)
      fvarLabels: ID GB_ACC ... Gene Ontology Molecular Function (16 total)
      fvarMetadata: Column Description labelDescription
    experimentData: use 'experimentData(object)'
      pubMedIds:19855080 
    Annotation: GPL8300

使用 Biobase 包的函数，获取样本的表型数据

> show(pData(phenoData(gset[[1]]))[,1:5])
                      title geo_accession                status submission_date last_update_date
    GSM296630     Lin-CD34-     GSM296630 Public on Jan 052010     Jun 042008      Jan 052010
    GSM296635 CML Lin-CD34-     GSM296635 Public on Jan 052010     Jun 042008      Jan 052010
    GSM296636 CML Lin-CD34+     GSM296636 Public on Jan 052010     Jun 042008      Jan 052010
    GSM296637     Lin-CD34+     GSM296637 Public on Jan 052010     Jun 042008      Jan 052010
    GSM296638     Lin+CD34+     GSM296638 Public on Jan 052010     Jun 042008      Jan 052010
    GSM296639 CML Lin+CD34+     GSM296639 Public on Jan 052010     Jun 042008      Jan 052010

###### **3.2 将 GDS 转换为 ExpressionSet** ######

对于上面的 gds 对象，我们可以简单地进行转换

eset <- GDS2eSet(gds,do.log2=TRUE)

现在，eset 是一个 ExpressionSet 类型，它包含了 GEO 数据集信息，以及样本信息

> eset
    ExpressionSet (storageMode: lockedEnvironment)
    assayData:22645 features,17 samples 
      element names: exprs 
    protocolData: none
    phenoData
      sampleNames: GSM11815 GSM11832 ... GSM12448 (17 total)
      varLabels: sample disease.state individual description
      varMetadata: labelDescription
    featureData
      featureNames:200000_s_at 200001_at ... AFFX-TrpnX-M_at (22645 total)
      fvarLabels: ID Gene title ... GO:Component ID (21 total)
      fvarMetadata: Column labelDescription
    experimentData: use 'experimentData(object)'
      pubMedIds:14641932 
    Annotation:

获取样本信息

> pData(eset)[,1:3]
               sample disease.state individual
    GSM11815 GSM11815           RCC        035
    GSM11832 GSM11832           RCC        023
    GSM12069 GSM12069           RCC        001
    GSM12083 GSM12083           RCC        005
    GSM12101 GSM12101           RCC        011
    GSM12106 GSM12106           RCC        032
    GSM12274 GSM12274           RCC          2
    GSM12299 GSM12299           RCC          3
    GSM12412 GSM12412           RCC          4
    GSM11810 GSM11810        normal        035
    GSM11827 GSM11827        normal        023
    GSM12078 GSM12078        normal        001
    GSM12099 GSM12099        normal        005
    GSM12269 GSM12269        normal          1
    GSM12287 GSM12287        normal          2
    GSM12301 GSM12301        normal          3
    GSM12448 GSM12448        normal          4

###### 3.3 将 GDS 转换为 MAList ######

由于 ExpressionSet 没有包含基因信息，所以没有检索到任何注释信息（也称为平台信息）

不过，要获得这些信息也很容易。首先，我们需要知道这个 GDS 使用的是什么平台。然后，再调用 getGEO 就可以得到我们需要的信息。

先获取 GDS 的平台

> Meta(gds)$platform
    [1]"GPL97"

再使用 getGEO 获取该平台的信息

gpl <- getGEO(filename=system.file("extdata/GPL97.annot.gz",package="GEOquery"))

gpl 变量现在包含了来自 GEO 的 GPL97 的信息。

与 ExpressionSet 不同的是，limma MAList 是存储了基因的注释信息，我们可以在调用 GDS2MA 时将我们新创建的 GPL 结构对象 gpl 传入进去

MA <- GDS2MA(gds,GPL=gpl)

查看 MA 类型信息

> class(MA)[1]"MAList"
    attr(,"package")
    [1]"limma"

MA 是 MAList 类型，不仅包含数据，还包含与 GDS507 相关的样本信息和基因信息

###### 3.4 将 GSE 转换为 ExpressionSet ######

首先，请确保使用 3.1 的方法是无法满足分析需求的。因为该方法快速且简便。

如果确实无法满足分析需求，那么可以使用下面的方法

将 GSE 对象转换为 ExpressionSet 对象需要一点点 R 基础数据操作，因为 GSE 的底层 GSM 和 GPL 对象中可以存储各种数据。

首先，我们需要确保所有的 GSM 都来自同一个平台。

> gsmplatforms <- lapply(GSMList(gse),function(x) {Meta(x)$platform_id})
    > head(gsmplatforms)
    $GSM11805
    [1] "GPL96"
    
    $GSM11810
    [1] "GPL97"
    
    $GSM11814
    [1] "GPL96"
    
    $GSM11815
    [1] "GPL97"
    
    $GSM11823
    [1] "GPL96"
    
    $GSM11827
    [1] "GPL97"

这个例子包含两个平台 GPL96 和 GPL97。我们可以对 GSMlist 的结果进行过滤，提取出 GPL96 的样本

> gsmlist <- Filter(function(gsm){Meta(gsm)$platform_id=='GPL96'},GSMList(gse))
    >length(gsmlist)
    [1]17

现在我们想知道哪一列数据是我们想要的，我们可以先查看某一个 GSM 的 Table 的前几行信息

> head(Table(gsmlist[[1]]))
               ID_REF  VALUE ABS_CALL
    1  AFFX-BioB-5_at  953.9        P
    2  AFFX-BioB-M_at 2982.8        P
    3  AFFX-BioB-3_at 1657.9        P
    4  AFFX-BioC-5_at 2652.7        P
    5  AFFX-BioC-3_at 2019.5        P
    6 AFFX-BioDn-5_at 3531.5        P

现在，我们知道了 VALUE 列是我们想要的，然后使用下面的代码来将这些数据转换为矩阵形式

# 先获取探针集
    probesets <- Table(GPLList(gse)[[1]])$ID
    # 将每个样本的 VALUE 列按列拼接起来
    data.matrix <- do.call('cbind',
                            lapply(gsmlist,function(x) {
                            tab <- Table(x)
                            # 我们使用 match 保证每列的探针顺序一致
                            mymatch <- match(probesets,tab$ID_REF)
                            return(tab$VALUE[mymatch])
                            })
                          )
    data.matrix <- apply(data.matrix,2,function(x) {as.numeric(as.character(x))})
    data.matrix <- log2(data.matrix)

查看结果

> head(data.matrix)
          GSM11805  GSM11814  GSM11823  GSM11830  GSM12067  GSM12075  GSM12079  GSM12098  GSM12100
    [1,] 10.926963 11.105254 11.275019 11.438636 11.424376 11.222795 11.469845 10.823367 10.835971
    [2,]  5.749534  7.908092  7.093814  7.514122  7.901470  6.407693  5.165912  6.556123  8.207014
    [3,]  7.066089  7.750205  7.244126  7.962896  7.337176  6.569856  7.477354  7.708739  7.428779
    [4,] 12.660353 12.479755 12.215897 11.458355 11.397568 12.529870 12.240046 12.336534 11.762839
    [5,]  6.195741  6.061776  6.565293  6.583459  6.877744  6.652486  3.981853  5.501439  6.247928
    [6,]  9.251956  9.480790  8.774458  9.878817  9.321252  9.275892  9.802355  9.046578  9.414474
          GSM12105  GSM12268  GSM12270  GSM12283  GSM12298  GSM12300  GSM12399  GSM12444
    [1,] 10.810893 11.062653 10.323055 11.181028 11.566387 11.078151 11.535178 11.105450
    [2,]  6.816344  6.563768  7.353147  5.770829  6.912889  4.812498  7.471675  7.488644
    [3,]  7.754888  7.126188  8.742815  7.339850  7.602142  7.383704  7.432959  7.381110
    [4,] 11.237509 12.412490 11.213408 12.678380 12.232901 12.090939 11.421802 12.172834
    [5,]  6.017922  6.525129  6.683696  5.918863  5.837943  6.281698  5.419539  5.469235
    [6,]  9.030115  9.252665  9.631359  8.656782  8.920948  8.629357  9.526695  9.494656

最后，我们将其转换为 ExpressionSet 对象

require(Biobase)
    # 构建 ExpressionSet 对象
    rownames(data.matrix)<- probesets
    colnames(data.matrix)<- names(gsmlist)
    pdata <- data.frame(samples=names(gsmlist))
    rownames(pdata)<- names(gsmlist)
    pheno <- as(pdata,"AnnotatedDataFrame")
    eset2 <- new('ExpressionSet',exprs=data.matrix,phenoData=pheno)

查看对象

> eset2
    ExpressionSet (storageMode: lockedEnvironment)
    assayData:22283 features,17 samples 
      element names: exprs 
    protocolData: none
    phenoData
      sampleNames: GSM11805 GSM11814 ... GSM12444 (17 total)
      varLabels: samples
      varMetadata: labelDescription
    featureData: none
    experimentData: use 'experimentData(object)'
    Annotation:

##### 4 获取原始数据 #####

NCBI GEO 会保存有一些原始数据，如 .CEL、.CDF 等文件。有时我们想快速的获取这个文件，可以使用 getGEOSuppFiles 函数

getGEOSuppFiles 函数接受一个 GEO accession number，然后下载与其相关的所有原始数据。

默认情况下，该函数会在当前工作目录下创建一个文件夹来存储这些数据。

getGEOSuppFiles('GSE11675')

获取下载的原始文件的信息

> eList <- getGEOSuppFiles('GSE11675', fetch_files =FALSE)
    > eList  fname url
    1 GSE11675_RAW.tar https://ftp.ncbi.nlm.nih.gov/geo/series/GSE11nnn/GSE11675/suppl//GSE11675_RAW.tar

下载的原始文件中包含 .CEL 或 .CDF 文件，我们可以使用 affy 或 oligo 包进行分析。

##### 5 使用示例 #####

GEOquery 对于快速收集大量数据来说是相当强大的，下面展示一个例子

###### 5.1 获取给定平台的所有 Series ######

有时我们想对某个平台的所有 GSE 数据进行数据挖掘，使用 GEOquery 可以让这一切变得非常简单。

但在使用之前，我们需要对 GPL 记录有一些了解。

gpl97 <- getGEO('GPL97')
    info <- Meta(gpl97)

获取 title

> info$title
    [1]"[HG-U133B] Affymetrix Human Genome U133B Array"

获取 Series ID

> head(info$series_id)[1]"GSE362""GSE473""GSE620""GSE674""GSE781""GSE907"> length(info$series_id)[1]165

获取样本号

> head(info$sample_id)[1]"GSM3922""GSM3924""GSM3926""GSM3928""GSM3930""GSM3932"> length(info$sample_id)[1]7917

该平台共包含了 165 个 Series，样本总量为 7917。获取到了这些信息，就可以进行批量下载了。

我们以下载前 5 个样本为例，进行说明

> gsmids <- Meta(gpl97)$sample_id
    > gsmlist <- sapply(gsmids[1:5], getGEO)> names(gsmlist)[1]"GSM3922""GSM3924""GSM3926""GSM3928""GSM3930"

[Link 1]: https://www.ncbi.nlm.nih.gov/
[https_www.ncbi.nlm.nih.gov_geo]: https://www.ncbi.nlm.nih.gov/geo/
[format_png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/26/e68b9fca6b38459aa107c98393d824b2.png
[format_png 1]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/26/a5a7730529ae438c923324a187ec0a0e.png
[a7bad55fb2111f52f32332a0c181249a.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/26/749100d764cf4ff191965783aa49e67f.png
[f63b4b7e09936271a2c735f7e5efdb16.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/26/4d8044989ded444e80e9359d2535618c.png
[format_png 2]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/26/0048870d23bf490c8b7889cbcfcfe4b0.png