前言
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),有 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 入門
從 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ǔ)下載的文件名被輸出到屏幕上(但沒有保存到任何地方),以便以后調(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
- 描述
該函數(shù)是 GEOquery 包中主要的用戶級(jí)函數(shù)。它引導(dǎo)下載(如果沒有指定文件名)并將 GEO SOFT 格式文件解析為 R 數(shù)據(jù)結(jié)構(gòu),該結(jié)構(gòu)是專門設(shè)計(jì)的,便于訪問(wèn) GEO SOFT 格式的每個(gè)重要部分。
- 使用
getGEO(GEO = NULL, filename = NULL, destdir = tempdir(),
GSElimits = NULL, GSEMatrix = TRUE,
AnnotGPL = FALSE, getGPL = TRUE,
parseCharacteristics = TRUE)
- 參數(shù)

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ù)部分,但是它沒有 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 沒有包含基因信息,所以沒有檢索到任何注釋信息(也稱為平臺(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"