在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))
本文分享自 MedBioInfoCloud 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!