要說(shuō)尋找公共芯片測(cè)序數(shù)據(jù)大家都知道上 GEO 查找派歌，其實(shí) EBI 也有個(gè)叫 ArrayExpress 的網(wǎng)站(網(wǎng)址ArrayExpress < EMBL-EBI)托管了大量的芯片數(shù)據(jù)鸥咖。同時(shí)還提供了同名 ArrayExpress R包在 Bioconductor 上衅鹿，像 GEOquery 下載和整理來(lái)自 GEO 的芯片數(shù)據(jù)一樣铅乡，下載和整理 ArrayExpress 上數(shù)據(jù)延柠。
這篇教程用 E-MTAB-1940 數(shù)據(jù)集為例展示 ArrayExpress 數(shù)據(jù)操作過(guò)程爬橡，具體的代碼解釋看注釋忘闻。
安裝R包并導(dǎo)入
# 安裝 ArrayExpress 包
BiocManager::install("ArrayExpress")
# 導(dǎo)入R包
library(ArrayExpress, quietly=TRUE)
library(tidyverse, quietly=TRUE)
library(affy, quietly=TRUE)
一鍵獲取數(shù)據(jù)使用 ArrayExpress 函數(shù)。
# 一定要記得指定路徑path, save參數(shù)表示運(yùn)行結(jié)束后是否保留下載文件
eset <- ArrayExpress(accession="E-MTAB-1940", path=".", save=TRUE)
從我自身體驗(yàn)感覺(jué)這方式不好酿炸，下載速度太慢了瘫絮。所以我自己去網(wǎng)站找到對(duì)應(yīng)數(shù)據(jù)集，使用瀏覽器或者復(fù)制鏈接后在服務(wù)器用 wget 下載好填硕。
使用 getAE 函數(shù)來(lái)下載數(shù)據(jù)麦萤，也可以用來(lái)讀取手動(dòng)下載好的本地?cái)?shù)據(jù)。用 local=TRUE 表示讀取本地下載好數(shù)據(jù)扁眯， sourcedir 是本地存儲(chǔ)位置频鉴。
# 如果是下載數(shù)據(jù)，使用 type 參數(shù)控制下載哪些
# 默認(rèn) full 下載所有數(shù)據(jù)
# raw 下載原始數(shù)據(jù)
# processed 下載處理好數(shù)據(jù)
> ae <- getAE(accession="E-MTAB-1940", path=".", type="raw", local=TRUE, sourcedir=".")
Unpacking data files
Warning message:
In getAE(accession = "E-MTAB-1940", path = ".", type = "raw", local = TRUE, :
No processed data files found in directory .
> str(ae)
List of 8
$ path : chr "."
$ rawFiles : chr [1:86] "FR_196_U133_2.CEL" "FR_327_U133_2.CEL" "FRI_328GRI_b_U133_2.CEL" "FR_329_U133_2.CEL" ...
$ rawArchive : chr [1:6] "E-MTAB-1940.raw.1.zip" "E-MTAB-1940.raw.2.zip" "E-MTAB-1940.raw.3.zip" "E-MTAB-1940.raw.4.zip" ...
$ processedFiles : NULL
$ processedArchive: chr(0)
$ sdrf : chr "E-MTAB-1940.sdrf.txt"
$ idf : chr "E-MTAB-1940.idf.txt"
$ adf : chr "A-AFFY-44.adf.txt"
下載后用 ae2bioc 函數(shù)讀取數(shù)據(jù)恋拍，用 methods 查看可用的方法。
> expr <- ae2bioc(ae)
> class(expr)
[1] "ExpressionFeatureSet"
attr(,"package")
[1] "oligoClasses"
> methods(class="ExpressionFeatureSet")
[1] [ [[ [[<- $
[5] $<- abstract annotation annotation<-
[9] assayData assayData<- backgroundCorrect bg
[13] bg<- bgSequence boxplot channel
[17] channelNames channelNames<- classVersion classVersion<-
[21] coerce combine db description
[25] description<- dim dimnames dimnames<-
[29] dims experimentData experimentData<- exprs
[33] exprs<- fData fData<- featureData
[37] featureData<- featureNames featureNames<- fvarLabels
[41] fvarLabels<- fvarMetadata fvarMetadata<- genomeBuild
[45] geometry getPlatformDesign getX getY
[49] hist image initialize intensity
[53] isCurrent isVersioned kind manufacturer
[57] MAplot mm mm<- mmindex
[61] mmSequence normalize notes notes<-
[65] paCalls pData pData<- phenoData
[69] phenoData<- pm pm<- pmChr
[73] pmindex pmPosition pmSequence preproc
[77] preproc<- probeNames probesetNames protocolData
[81] protocolData<- pubMedIds pubMedIds<- rma
[85] runDate sampleNames sampleNames<- selectChannels
[89] show storageMode storageMode<- updateObject
[93] updateObjectTo varLabels varLabels<- varMetadata
[97] varMetadata<-
用 rma 進(jìn)行 RMA normalize得到ExpressionSet對(duì)象藕甩。
> rmae <- rma(expr)
Background correcting
Normalizing
Calculating Expression
像處理GEO芯片數(shù)據(jù)的 GEOquery 一樣用 exprs 函數(shù)取得表達(dá)矩陣, pData 取得其他實(shí)驗(yàn)信息施敢。
> probe_expr <- exprs(rmae) %>% as_tibble(rownames="probe_id")
> head(probe_expr, n=2)
# A tibble: 2 x 87
probe_id FR_196_U133_2.C… FR_327_U133_2.C… FRI_328GRI_b_U1… FR_329_U133_2.C…
<chr> <dbl> <dbl> <dbl> <dbl>
1 1007_s_… 9.99 10.3 10.4 9.92
2 1053_at 5.65 6.33 6.22 5.76
# … with 82 more variables: FR_46_U133_2.CEL <dbl>, FRI_47BEN_U133_2.CEL <dbl>,
# 省略后續(xù)輸出
pdata <- pData(rmae) %>% as_tibble()
芯片可能以后用得會(huì)越來(lái)越少，但是如果有人經(jīng)常進(jìn)行這些數(shù)據(jù)挖掘的狭莱，可以更深入去學(xué)習(xí)僵娃，想要的建議把 affy 包的文檔讀一遍。像剛剛用到 rma 就是來(lái)自于 affy 包腋妙。
歡迎關(guān)注我的微信公眾號(hào) Hello BioInfo
