基因芯片數據分析（二）：讀取芯片數據

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))