基因芯片数据分析（二）：读取芯片数据

2019 年 12 月 13 日
笔记

在microarray的处理中，第一步就是读取数据。无论是自己的保存在本地的数据，还是在线保存的数据，对于不同公司的芯片可以使用不同的软件包读取。在这里，我们说的在线数据，主要是指保存在GEO (Gene Expression Omnibus) 数据库中的数据，当然GEO的数据可先下载后再读入。

读取本地数据

我们可以使用affy包中的ReadAffy函数读取cel文件。

library(affy) ##加载库文件  Data <- ReadAffy("/path/of/CELs") ##读取工作目录下的CEL文件

ReadAffy()格式如下：

ReadAffy(..., filenames=character(0),                widget=getOption("BioC")$affy$use.widgets,                compress=getOption("BioC")$affy$compress.cel,                celfile.path=NULL,                sampleNames=NULL,                phenoData=NULL,                description=NULL,                notes="",                rm.mask=FALSE, rm.outliers=FALSE, rm.extra=FALSE,                verbose=FALSE,sd=FALSE, cdfname = NULL)

参数	说明
…	需要读入的文件名，可以用逗号间隔，输入多个CEL文件
filenames	文件名列表构成的字符向量(character vector)
phenoData	an AnnotatedDataFrame object, a character of length one, or a data.frame.
description	MIAME对象，它是对microarray实验的完整描述，包括名称，实验室，通信方式，实验摘要，网址，样品，杂交信息，对照，预处理，PubMedId，等等。除继承的方法外，还提供了提取相应信息的方法。
notes	注释
compress	CEL文件是否为压缩文件，支持zip和gzip
rm.mask	should the spots marked as ‘MASKS’ set to NA?
rm.outliers	should the spots marked as ‘OUTLIERS’ set to NA?
rm.extra	if TRUE, then overrides what is in rm.mask and rm.oultiers.
verbose	verbosity flag.
widget	a logical specifying if widgets should be used.
celfile.path	文件所在的目录，缺省时为R的工作目录
sampleNames	样品名列表构成的字符向量(character vector)
sd	是否读入CEL文件中的标准差？默认不读入，可以节省大量的内存。
cdfname	指定CDF库的文件名。如果设置为NULL，程序会自动从标准库中下载。

对于Affymetrix Exon/Gene ST Arrays，我们不能使用affy包来读取，我们需要使用oligo或者xps来进行分析。这里介绍oligo包。

library(oligo)  geneCELs <- list.celfiles("path/to/cels", full.names=TRUE)  affyGeneFS <- read.celfiles(geneCELs)

对于NimbleGen数据（XYS文件）,可以使用oligo读取。

library(oligo)  xysFiles <- list.celfiles('myXYSs', full.names=TRUE)  rawData <- read.xysfiles(xysFiles)

对于Agilent的gpr文件或者txt文件，可以使用limma的read.maimages来读取。需要注意的是，对于不同的文件source参数有多种选择："generic", "agilent", "agilent.median", "agilent.mean", "arrayvision", "arrayvision.ARM", "arrayvision.MTM", "bluefuse", "genepix", "genepix.custom", "genepix.median", "imagene", "imagene9", "quantarray", "scanarrayexpress", "smd.old", "smd", "spot" or "spot.close.open"。对于单色芯片，注意将green.only设置成TRUE.

library(limma)  data <- read.maimages(files=filelist, source = "agilent")

读取在线数据

在GEO数据库中保存有大量的microarray的原始数据。许多文章在发表之前，作者为了提高文章的可重复性，都会将高通量的数据提交至GEO数据库当中，以方便审稿人以及公从读者调验。

本文以GSE46106数据为例，讲述如何从GEO上下载数据。分为两种方式，第一，直接从GEO上下载表达数据，第二，直接从GEO上下载CEL文件，然后以读取本地数据的方式读入。

首先我们示例如何下载表达数据。

library(GEOquery)  gset <- getGEO("GSE46106", GSEMatrix =TRUE)  length(gset)

gset只是从geo从抓回的数据，它可能是多个数据，所以返回结果保存在了一个list当中。为了以后操作的方便，对于返回长度为1的list，我将其结果从list中抽取出来。就是类似：a <- list(a=c(1,2,3)) 以后每次访问c(1,2,3)，我都要写成a[[1]]这样，感觉不方便，于是 a <- a[[1]] 这样以后访问c(1,2,3)就只需要写成a就可以了。

gset <- gset[[1]]  head(pData(gset)[,1:5])  # load NCBI platform annotation  gpl <- annotation(gset)  platf <- getGEO(gpl, AnnotGPL=TRUE)  ncbifd <- data.frame(attr(dataTable(platf), "table"))  eset <- exprs(gset)  head(eset[,1:5])  head(ncbifd[,1:5])

其次我们示例如何下载原始的CEL文件。

getGEOSuppFiles("GSE46106")  setwd("GSE46106/")  dir()  untar("GSE46106_RAW.tar")  files <- dir(pattern="gz$")  sapply(files, gunzip)  filelist <- dir(pattern="CEL$")  library(affy)  library(annotate)  data <- ReadAffy(filenames=filelist)  affydb<-annPkgName(data@annotation,type="db")  require(affydb, character.only=TRUE)  eset<-rma(data,verbose=FALSE)  eset.e <- exprs(eset)  library(annaffy)  symbols<-as.character(aafSymbol(as.character(rownames(eset)),affydb))  genes<-as.character(aafUniGene(as.character(rownames(eset)),affydb))