如何将rma（）规范化应用于唯一的CEL文件？

Question

我已经实现了一个R脚本，它对基因表达数据集执行批量修正。要进行批量校正，我首先需要通过Bioconductor的Affy rma() function对每个CEL文件中的数据进行标准化。

如果我在从GEO获得的GSE59867 dataset上运行它，一切正常。我将批处理定义为数据收集日期：我将具有相同日期的所有CEL文件放入特定文件夹，然后将该日期/文件夹视为特定批处理。在GSE59867 dataset上，批处理/文件夹只包含1个CEL文件。尽管如此，rma()功能完美地适用于它。

但是，相反，如果我尝试在另一个数据集（GSE36809）上运行我的脚本，我有一些麻烦：如果我尝试将rma()函数应用于仅包含1个文件的批处理/文件夹，我会收到以下错误：

Error in `colnames<-`(`*tmp*`, value = "GSM901376_c23583161.CEL.gz") : 
  attempt to set 'colnames' on an object with less than two dimensions

这是我特定的R代码，让您了解。你首先要下载文件GSM901376_c23583161.CEL.gz：

setwd(".")
options(stringsAsFactors = FALSE)

fileURL <- "ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM901nnn/GSM901376/suppl/GSM901376%5Fc23583161%2ECEL%2Egz"
fileDownloadCommand <- paste("wget ", fileURL, " ", sep="")
system(fileDownloadCommand)

图书馆安装：

source("https://bioconductor.org/biocLite.R")    
list.of.packages <- c("easypackages")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)    
listOfBiocPackages <- c("oligo", "affyio","BiocParallel")

bioCpackagesNotInstalled <- which( !listOfBiocPackages %in% rownames(installed.packages()) )
cat("package missing listOfBiocPackages[", bioCpackagesNotInstalled, "]: ", listOfBiocPackages[bioCpackagesNotInstalled], "\n", sep="")

if( length(bioCpackagesNotInstalled) ) {
    biocLite(listOfBiocPackages[bioCpackagesNotInstalled])
}

library("easypackages")
libraries(list.of.packages)
libraries(listOfBiocPackages)

rma（）的应用

thisFileDate <- "GSM901376_c23583161.CEL.gz"
thisDateRawData <- read.celfiles(thisDateCelFiles)
thisDateNormData <- rma(thisDateRawData)

在调用rma()之后，我收到错误。我怎么解决这个问题？

我还尝试通过直接保存thisDateRawData对象来跳过此规范化。但后来我遇到的问题是我无法将这个thisDateRawData（即ExpressionFeatureSet）与rma（）的输出结合在一起（即ExpressionSet对象）。

（编辑：我对这个问题进行了大量编辑，并添加了一段你应该可以在你的电脑上运行的R代码。）

Answer 1

嗯。这是一个令人费解的问题。对于具有单个样本的类GeneFeatureSet，oligo::rma()函数可能是错误的。我通过使用较低级别的函数使用单个样本，但这意味着我还必须通过指定插槽从头创建表达式集：

# source("https://bioconductor.org/biocLite.R")
# biocLite("GEOquery")
# biocLite("pd.hg.u133.plus.2")
# biocLite("pd.hugene.1.0.st.v1")
library(GEOquery)
library(oligo)

# # Instead of using .gz files, I extracted the actual CELs.
# # This is just to illustrate how I read in the files; your usage will differ.
# projectDir <- ""  # Path to .tar files here
# setwd(projectDir)
# untar("GSE36809_RAW.tar", exdir = "GSE36809")
# untar("GSE59867_RAW.tar", exdir = "GSE59867")
# setwd("GSE36809"); gse3_cels <- dir()
# sapply(paste(gse3_cels, sep = "/"), gunzip); setwd(projectDir)
# setwd("GSE59867"); gse5_cels <- dir()
# sapply(paste(gse5_cels, sep = "/"), gunzip); setwd(projectDir)
#
# Read in CEL
#
# setwd("GSE36809"); gse3_cels <- dir()
# gse3_efs <- read.celfiles(gse3_cels[1])

# # Assuming you've read in the CEL files as a GeneFeatureSet or 
# # ExpressionFeatureSet object (i.e. gse3_efs in this example),
# # you can now fit the RMA and create an ExpressionSet object with it:

exprsData <- basicRMA(exprs(gse3_efs), pnVec = featureNames(gse3_efs))
gse3_expset <- new("ExpressionSet")
slot(gse3_expset, "assayData") <- assayDataNew(exprs = exprsData)
slot(gse3_expset, "phenoData") <- phenoData(gse3_efs)
slot(gse3_expset, "featureData") <- annotatedDataFrameFrom(attr(gse3_expset, 
  'assayData'), byrow = TRUE)
slot(gse3_expset, "protocolData") <- protocolData(gse3_efs)
slot(gse3_expset, "annotation") <- slot(gse3_efs, "annotation")

希望上述方法适用于您的代码。