GEOquery 下載 GEO 數(shù)據(jù)

前言

NCBI Gene Expression Omnibus(基因表達(dá)綜合數(shù)據(jù)庫(kù),GEO)公開了很多高通量基因組學(xué)數(shù)據(jù)集。盡管名字很好聽,但這個(gè)數(shù)據(jù)庫(kù)并不是專門針對(duì)基因表達(dá)數(shù)據(jù)的。

NCBI GEO

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

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

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

他們對(duì)兩種不同的因子(NFkB-IIPol II)進(jìn)行 RNA-seqChIP-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)行了匯總。

最后,GEOSeries 標(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 numberGPLxxx)。一個(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 numberGSExxx)。系列記錄有幾種不同的格式,均由 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 入門

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 類),表示 GEOGDS507 條目。

您會(huì)注意到用于存儲(chǔ)下載的文件名被輸出到屏幕上(但沒有保存到任何地方),以便以后調(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)下載(如果沒有指定文件名)并將 GEO SOFT 格式文件解析為 R 數(shù)據(jù)結(jié)構(gòu),該結(jié)構(gòu)是專門設(shè)計(jì)的,便于訪問(wèn) GEO SOFT 格式的每個(gè)重要部分。

  1. 使用
getGEO(GEO = NULL, filename = NULL, destdir = tempdir(),
       GSElimits = NULL, GSEMatrix = TRUE, 
       AnnotGPL = FALSE, getGPL = TRUE,
       parseCharacteristics = TRUE)
  1. 參數(shù)
image.png

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

GEOquery 的數(shù)據(jù)結(jié)構(gòu)其實(shí)有兩種形式。第一種,由 GDSGPLGSM 組成,它們的行為都是相似的,訪問(wèn)器對(duì)每一種都有相似的作用。

第四種 GEOquery 數(shù)據(jù)結(jié)構(gòu) GSE,是由 GSMGPL 對(duì)象組合而成的復(fù)合數(shù)據(jù)類型。我先把前三個(gè)一起解釋一下。

2.1 GDS、GSM 和 GPL 類

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

GEODataTable 有兩個(gè)簡(jiǎn)單的部分: ColumnsTable,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 類

GSEGEO 中最容易混淆的。一個(gè) GSE 條目可以代表在任意數(shù)量的平臺(tái)上運(yùn)行的任意數(shù)量的樣本。

GSE 類和其他類一樣,有一個(gè)元數(shù)據(jù)部分,但是它沒有 GEODataTable。

相反,它包含了兩個(gè)列表,可以使用 GPLListGSMList 方法進(jìn)行訪問(wèn),這兩個(gè)列表分別是 GPLGSM 對(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ù)集與 limmaMAList 數(shù)據(jù)結(jié)構(gòu)和 BiobaseExpressionSet 數(shù)據(jù)結(jié)構(gòu)非常相似

存在兩個(gè)函數(shù) GDS2MAGDS2eSet,能夠?qū)?GEOquery 數(shù)據(jù)結(jié)構(gòu)分別轉(zhuǎn)換為 MAListExpressionSet

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 沒有包含基因信息,所以沒有檢索到任何注釋信息(也稱為平臺(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)自 GEOGPL97 的信息。

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"

MAMAList 類型,不僅包含數(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 的底層 GSMGPL 對(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) GPL96GPL97。我們可以對(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è) GSMTable 的前幾行信息

> 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 文件,我們可以使用 affyoligo 包進(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"
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡(jiǎn)書系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容