前言

NCBI Gene Expression Omnibus（基因表達(dá)綜合數(shù)據(jù)庫(kù)，GEO）公開(kāi)了很多高通量基因組學(xué)數(shù)據(jù)集榨汤。盡管名字很好聽(tīng)，但這個(gè)數(shù)據(jù)庫(kù)并不是專門(mén)針對(duì)基因表達(dá)數(shù)據(jù)的怎茫。

NCBI GEO

NCBI GEO 是以樣本的形式組織起來(lái)的收壕，這些樣本被歸入 series。對(duì)于較大的實(shí)驗(yàn)轨蛤，有 SubSeries 和 SuperSeries 兩種

一個(gè) SuperSeries 就是一篇論文的所有實(shí)驗(yàn)蜜宪，一個(gè) SuperSeries 可以分解為不同技術(shù)的 SubSeries。

舉個(gè)例子祥山，看看 SuperSeries GSE19486圃验。在這篇論文中，他們使用了 2 個(gè)不同的平臺(tái)（這是個(gè)奇怪的名字缝呕，平臺(tái)是一種技術(shù)和物種的組合）

他們對(duì)兩種不同的因子（NFkB-II 和 Pol II）進(jìn)行 RNA-seq 和 ChIP-seq澳窑。這導(dǎo)致 4 個(gè)子系列（2 個(gè)用于 RNA-seq，2 個(gè)用于 ChIP-seq）

一個(gè)更簡(jiǎn)單的例子供常，如 GSE994摊聋，使用 Affymetrix 微陣列比較來(lái)自當(dāng)前和以前吸煙者的樣本

提交給 NCBI GEO 的數(shù)據(jù)既可以是 "raw"，也可以是 "proceseed"栈暇。我們暫且把重點(diǎn)放在基因表達(dá)數(shù)據(jù)上麻裁。"proceseed" 的數(shù)據(jù)是經(jīng)過(guò)標(biāo)準(zhǔn)化和量化的，通常是在基因水平上，以基因?qū)?yīng)樣本的矩陣形式提供煎源。

"raw" 數(shù)據(jù)可以是任何東西色迂，從測(cè)序讀數(shù)到微陣列圖像文件。甚至可以有不同狀態(tài)的 "raw" 數(shù)據(jù)薪夕，例如對(duì)于一個(gè) RNA-seq 數(shù)據(jù)集脚草，可能包含：

• FASTQ 文件（原始 read）
• BAM 文件（比對(duì) read）
• 基因樣本表達(dá)矩陣(非標(biāo)準(zhǔn)化)
• 基因樣本表達(dá)矩陣(標(biāo)準(zhǔn)化)

所以赫悄，什么是 "raw"原献，什么是 "proceseed"，可能非常依賴于背景信息埂淮。

在某些情況下姑隅，該領(lǐng)域有一個(gè)共識(shí)。對(duì)于 Affymetrix 基因表達(dá)芯片來(lái)說(shuō)倔撞，"raw" 文件是所謂的 CEL 文件（Affymetrix 發(fā)明的一種文件格式）讲仰，"proceseed" 數(shù)據(jù)是經(jīng)過(guò)歸一化和量化的數(shù)據(jù)，在探針集層面進(jìn)行了匯總痪蝇。

最后，GEO 有 Series 標(biāo)識(shí)符（如 GSE19486）和樣本標(biāo)識(shí)符（GSM486297）躏啰。用戶幾乎總是對(duì)一個(gè)給定 Series 中的所有樣品感興趣给僵，所以系列標(biāo)識(shí)符毫捣，也叫 accession number。

Platforms(平臺(tái))

平臺(tái)記錄了陣列上的元素列表（如 cDNA蔓同、寡核苷酸探針組蹲诀、ORF斑粱、抗體）或該實(shí)驗(yàn)中可能被檢測(cè)和量化的元素列表（如 SAGE 標(biāo)簽、肽）则北。

每個(gè)平臺(tái)記錄都有一個(gè)唯一的 GEO accession number（GPLxxx）披粟。一個(gè)平臺(tái)可以包含多個(gè)提交者提交的樣品守屉。

Samples(樣本)

一份樣品記錄描述了處理單個(gè)樣品的條件、所涉及的操作以及從中獲得的每種元素的豐度測(cè)量滨巴。

每個(gè)樣品記錄都有一個(gè)唯一的 GEO accession number(GSMxxx)恭取。一個(gè)樣品只能屬于一個(gè)平臺(tái)，但是可以被加入多個(gè)系列

Series(系列)

系列記錄定義了一組相關(guān)的樣品耗跛，這些樣品是如何相關(guān)的调塌，以及它們是否排序或如何排序的惠猿。

系列記錄提供了整體實(shí)驗(yàn)的關(guān)注點(diǎn)和描述信息偶妖。還可以包含所提取的數(shù)據(jù)的描述趾访、簡(jiǎn)要結(jié)論或分析腹缩。

每個(gè)系列記錄都有一個(gè)唯一的 GEO accession number（GSExxx）。系列記錄有幾種不同的格式润讥，均由 GEOquery 獨(dú)立處理楚殿。

小而新的 GSEMatrix 文件的解析速度相當(dāng)快竿痰；GEOquery 使用一個(gè)參數(shù)來(lái)選擇是否使用 GSEMatrix 文件

Datasets(數(shù)據(jù)集)

GEO 數(shù)據(jù)集（GDSxxx）是 GEO 樣本數(shù)據(jù)的精選集。一個(gè) GDS 記錄代表了一個(gè)在生物學(xué)和統(tǒng)計(jì)學(xué)上可比較的 GEO 樣本的集合影涉，并構(gòu)成了 GEO 數(shù)據(jù)可視化和分析工具的基礎(chǔ)蟹倾。

GDS 內(nèi)的樣本來(lái)自同一個(gè)平臺(tái)，它們共享一組共同的探針集培慌。

假設(shè) GDS 內(nèi)每個(gè)樣品的測(cè)量值是以同樣方式計(jì)算的吵护，即考慮背景處理和歸一化等因素表鳍。

GEOquery

今天我們的重點(diǎn)并不是要介紹上面的內(nèi)容进胯，而是要和大家聊聊該怎么用 R 來(lái)下載 GEO 的數(shù)據(jù)胁镐。

我們使用的是 GEOquery 包盯漂，下來(lái)就來(lái)詳細(xì)介紹一下

安裝

從 Bioconductor 安裝

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

從 GitHub 安裝

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

使用

1 GEOquery 入門(mén)

從 GEO 中獲取數(shù)據(jù)真的很容易就缆。只需要一個(gè)命令: getGEO竭宰。這一個(gè)函數(shù)對(duì)其輸入進(jìn)行解釋切揭，以確定如何從 GEO 中獲取數(shù)據(jù)锁摔，然后將數(shù)據(jù)解析成有用的 R 數(shù)據(jù)結(jié)構(gòu)谐腰。

用法很簡(jiǎn)單十气。首先加載 GEOquery 庫(kù)

library(GEOquery)

現(xiàn)在砸西，我們可以自由訪問(wèn)任何 GEO accession number。請(qǐng)注意竟闪，在下文中炼蛤，我使用了一個(gè)與 GEOquery 包一起打包的文件理朋。

一般來(lái)說(shuō)嗽上，你只要傳入正確的 GEO accession 就行兽愤，正如代碼注釋中所指出的那樣

# 更常用的方式是傳入 GEO accession浅萧，從 GEO 數(shù)據(jù)庫(kù)中下載數(shù)據(jù):
# gds <- getGEO("GDS507")
gds <- getGEO(filename=system.file("extdata/GDS507.soft.gz",package="GEOquery"))

現(xiàn)在洼畅，gds 包含了 R 數(shù)據(jù)結(jié)構(gòu)(GDS 類)，表示 GEO 的 GDS507 條目徘郭。

您會(huì)注意到用于存儲(chǔ)下載的文件名被輸出到屏幕上（但沒(méi)有保存到任何地方）残揉，以便以后調(diào)用 getGEO(filename=...) 時(shí)使用冲甘。

# 更常用的方式是直接傳入 GSM 號(hào)江醇，從數(shù)據(jù)庫(kù)中下載：
# gds <- getGEO("GSM11805")
gsm <- getGEO(filename=system.file("extdata/GSM11805.txt.gz",package="GEOquery"))

1.1 getGEO

1. 描述

該函數(shù)是 GEOquery 包中主要的用戶級(jí)函數(shù)陶夜。它引導(dǎo)下載（如果沒(méi)有指定文件名）并將 GEO SOFT 格式文件解析為 R 數(shù)據(jù)結(jié)構(gòu)条辟，該結(jié)構(gòu)是專門(mén)設(shè)計(jì)的羽嫡，便于訪問(wèn) GEO SOFT 格式的每個(gè)重要部分杭棵。

2. 使用

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

3. 參數(shù)

image.png

2 GEOquery 數(shù)據(jù)結(jié)構(gòu)

GEOquery 的數(shù)據(jù)結(jié)構(gòu)其實(shí)有兩種形式魂爪。第一種滓侍，由 GDS撩笆、GPL 和 GSM 組成浇衬，它們的行為都是相似的耘擂，訪問(wèn)器對(duì)每一種都有相似的作用醉冤。

第四種 GEOquery 數(shù)據(jù)結(jié)構(gòu) GSE篙悯，是由 GSM 和 GPL 對(duì)象組合而成的復(fù)合數(shù)據(jù)類型鸽照。我先把前三個(gè)一起解釋一下矮燎。

2.1 GDS诞外、GSM 和 GPL 類

這些類中的每一個(gè)類都由一個(gè)元數(shù)據(jù)頭（與 SOFT 格式頭部信息幾乎一致）和一個(gè) GEODataTable 組成峡谊。

GEODataTable 有兩個(gè)簡(jiǎn)單的部分: Columns 和 Table，Columns 部分是對(duì) Table 部分的列的描述正什。此外埠忘，每個(gè)類還有一個(gè) show 方法莹妒。

例如旨怠，對(duì)于上面的 gsm鉴腻，查看元數(shù)據(jù)

> 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 相關(guān)的數(shù)據(jù)

> 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 中最容易混淆的。一個(gè) GSE 條目可以代表在任意數(shù)量的平臺(tái)上運(yùn)行的任意數(shù)量的樣本固棚。

GSE 類和其他類一樣此洲，有一個(gè)元數(shù)據(jù)部分黍翎，但是它沒(méi)有 GEODataTable匣掸。

相反，它包含了兩個(gè)列表霎匈，可以使用 GPLList 和 GSMList 方法進(jìn)行訪問(wèn)铛嘱，這兩個(gè)列表分別是 GPL 和 GSM 對(duì)象的列表墨吓。

舉個(gè)例子帖烘，我們從本地讀取 soft 文件

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

查看元數(shù)據(jù)

> 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 對(duì)象的名稱

> 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"

獲得列表中的第一個(gè) GSM 對(duì)象

> 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 轉(zhuǎn)換為表達(dá)譜

GEO 數(shù)據(jù)集與 limma 的 MAList 數(shù)據(jù)結(jié)構(gòu)和 Biobase 的 ExpressionSet 數(shù)據(jù)結(jié)構(gòu)非常相似

存在兩個(gè)函數(shù) GDS2MA 和 GDS2eSet式矫，能夠?qū)?GEOquery 數(shù)據(jù)結(jié)構(gòu)分別轉(zhuǎn)換為 MAList 和 ExpressionSet

3.1 獲取 GSE Series Matrix

除了解析較大的 soft 文件采转，我們也可以直接獲取處理后的 Series matrix 文件氏义。

getGEO 能夠快速解析該類型的文件惯悠，解析返回的數(shù)據(jù)結(jié)構(gòu)是 ExpressionSet 的列表

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

顯示文件結(jié)構(gòu)信息

> 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 包的函數(shù)克婶，獲取樣本的表型數(shù)據(jù)

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

使用 exprs 獲取處理后的表達(dá)譜數(shù)據(jù)

> head(exprs(gset[[1]]))
          GSM296630 GSM296635 GSM296636 GSM296637  GSM296638 GSM296639
1000_at   103.88700  138.0500 133.32200 242.86600  175.73400 165.28400
1001_at     3.39447   31.1682  98.97760  69.78960   56.16700 101.31200
1002_f_at  10.56290    3.4511   1.63724   5.92232    3.26876   3.33937
1003_s_at   8.42204   20.7501  18.65330  12.75090    5.53738   6.18855
1004_at     5.55025   14.6198  17.52940  14.25250   16.76820  16.43160
1005_at   802.71900  846.3760 761.27500 911.98200 2650.83000 801.71700

3.2 將 GDS 轉(zhuǎn)換為 ExpressionSet

對(duì)于上面的 gds 對(duì)象，我們可以簡(jiǎn)單地進(jìn)行轉(zhuǎn)換

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

現(xiàn)在筋岛，eset 是一個(gè) ExpressionSet 類型睁宰，它包含了 GEO 數(shù)據(jù)集信息柒傻，以及樣本信息

> 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 轉(zhuǎn)換為 MAList

由于 ExpressionSet 沒(méi)有包含基因信息青柄，所以沒(méi)有檢索到任何注釋信息（也稱為平臺(tái)信息）

不過(guò)致开，要獲得這些信息也很容易双戳。首先拣技，我們需要知道這個(gè) GDS 使用的是什么平臺(tái)膏斤。然后莫辨，再調(diào)用 getGEO 就可以得到我們需要的信息沮榜。

先獲取 GDS 的平臺(tái)

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

再使用 getGEO 獲取該平臺(tái)的信息

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

gpl 變量現(xiàn)在包含了來(lái)自 GEO 的 GPL97 的信息蟆融。

與 ExpressionSet 不同的是型酥，limma MAList 是存儲(chǔ)了基因的注釋信息弥喉，我們可以在調(diào)用 GDS2MA 時(shí)將我們新創(chuàng)建的 GPL 結(jié)構(gòu)對(duì)象 gpl 傳入進(jìn)去

MA <- GDS2MA(gds,GPL=gpl)

查看 MA 類型信息

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

MA 是 MAList 類型由境，不僅包含數(shù)據(jù)虏杰，還包含與 GDS507 相關(guān)的樣本信息和基因信息

3.4 將 GSE 轉(zhuǎn)換為 ExpressionSet

首先嘹屯，請(qǐng)確保使用 3.1 的方法是無(wú)法滿足分析需求的州弟。因?yàn)樵摲椒焖偾液?jiǎn)便婆翔。

如果確實(shí)無(wú)法滿足分析需求啃奴，那么可以使用下面的方法

將 GSE 對(duì)象轉(zhuǎn)換為 ExpressionSet 對(duì)象需要一點(diǎn)點(diǎn) R 基礎(chǔ)數(shù)據(jù)操作最蕾，因?yàn)?GSE 的底層 GSM 和 GPL 對(duì)象中可以存儲(chǔ)各種數(shù)據(jù)瘟则。

首先醋拧，我們需要確保所有的 GSM 都來(lái)自同一個(gè)平臺(tái)丹壕。

> 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"

這個(gè)例子包含兩個(gè)平臺(tái) GPL96 和 GPL97菌赖。我們可以對(duì) GSMlist 的結(jié)果進(jìn)行過(guò)濾琉用，提取出 GPL96 的樣本

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

現(xiàn)在我們想知道哪一列數(shù)據(jù)是我們想要的，我們可以先查看某一個(gè) 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

然后獲取對(duì)應(yīng)的列描述信息

> head(Columns(gsmlist[[1]]))
    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

現(xiàn)在辕羽，我們知道了 VALUE 列是我們想要的，然后使用下面的代碼來(lái)將這些數(shù)據(jù)轉(zhuǎn)換為矩陣形式

# 先獲取探針集
probesets <- Table(GPLList(gse)[[1]])$ID
# 將每個(gè)樣本的 VALUE 列按列拼接起來(lái)
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)

查看結(jié)果

> 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

最后垄惧，我們將其轉(zhuǎn)換為 ExpressionSet 對(duì)象

require(Biobase)
# 構(gòu)建 ExpressionSet 對(duì)象
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)

查看對(duì)象

> 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 獲取原始數(shù)據(jù)

NCBI GEO 會(huì)保存有一些原始數(shù)據(jù)刁愿，如 .CEL、.CDF 等文件到逊。有時(shí)我們想快速的獲取這個(gè)文件滤钱，可以使用 getGEOSuppFiles 函數(shù)

getGEOSuppFiles 函數(shù)接受一個(gè) GEO accession number，然后下載與其相關(guān)的所有原始數(shù)據(jù)脑题。

默認(rèn)情況下件缸，該函數(shù)會(huì)在當(dāng)前工作目錄下創(chuàng)建一個(gè)文件夾來(lái)存儲(chǔ)這些數(shù)據(jù)。

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 包進(jìn)行分析他炊。有空我們?cè)僭敿?xì)介紹其用法

5 使用示例

GEOquery 對(duì)于快速收集大量數(shù)據(jù)來(lái)說(shuō)是相當(dāng)強(qiáng)大的，下面展示一個(gè)例子

5.1 獲取給定平臺(tái)的所有 Series

有時(shí)我們想對(duì)某個(gè)平臺(tái)的所有 GSE 數(shù)據(jù)進(jìn)行數(shù)據(jù)挖掘已艰，使用 GEOquery 可以讓這一切變得非常簡(jiǎn)單痊末。

但在使用之前，我們需要對(duì) 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

獲取樣本號(hào)

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

> length(info$sample_id)
[1] 7917

該平臺(tái)共包含了 165 個(gè) Series凿叠，樣本總量為 7917。獲取到了這些信息嚼吞，就可以進(jìn)行批量下載了盒件。

我們以下載前 5 個(gè)樣本為例，進(jìn)行說(shuō)明

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