要说寻找公共芯片测序数据大家都知道上 GEO 查找，其实 EBI 也有个叫 ArrayExpress 的网站(网址ArrayExpress < EMBL-EBI)托管了大量的芯片数据。同时还提供了同名 ArrayExpress R包在 Bioconductor 上，像 GEOquery 下载和整理来自 GEO 的芯片数据一样，下载和整理 ArrayExpress 上数据。
这篇教程用 E-MTAB-1940 数据集为例展示 ArrayExpress 数据操作过程，具体的代码解释看注释。
安装R包并导入
# 安装 ArrayExpress 包
BiocManager::install("ArrayExpress")
# 导入R包
library(ArrayExpress, quietly=TRUE)
library(tidyverse, quietly=TRUE)
library(affy, quietly=TRUE)
一键获取数据使用 ArrayExpress 函数。
# 一定要记得指定路径path, save参数表示运行结束后是否保留下载文件
eset <- ArrayExpress(accession="E-MTAB-1940", path=".", save=TRUE)
从我自身体验感觉这方式不好，下载速度太慢了。所以我自己去网站找到对应数据集，使用浏览器或者复制链接后在服务器用 wget 下载好。
使用 getAE 函数来下载数据，也可以用来读取手动下载好的本地数据。用 local=TRUE 表示读取本地下载好数据， sourcedir 是本地存储位置。
# 如果是下载数据，使用 type 参数控制下载哪些
# 默认 full 下载所有数据
# raw 下载原始数据
# processed 下载处理好数据
> 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 函数读取数据，用 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 进行 RMA normalize得到ExpressionSet对象。
> rmae <- rma(expr)
Background correcting
Normalizing
Calculating Expression
像处理GEO芯片数据的 GEOquery 一样用 exprs 函数取得表达矩阵, pData 取得其他实验信息。
> 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>,
# 省略后续输出
pdata <- pData(rmae) %>% as_tibble()
芯片可能以后用得会越来越少，但是如果有人经常进行这些数据挖掘的，可以更深入去学习，想要的建议把 affy 包的文档读一遍。像刚刚用到 rma 就是来自于 affy 包。
欢迎关注我的微信公众号 Hello BioInfo
网友评论