GEO SRA ArrayExpress ENA 数据挖掘ArrayExpress

从 EBI ArrayExpress 获取芯片数据

2019-12-16  本文已影响0人  BeeBee生信

要说寻找公共芯片测序数据大家都知道上 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

上一篇 下一篇

猜你喜欢

热点阅读