如何将 rma() 标准化应用于唯一的 CEL 文件？

Posted 2023-03-23

【中文标题】如何将 rma() 标准化应用于唯一的 CEL 文件？

【英文标题】：How to apply rma() normalization to a unique CEL file?

【发布时间】：2019-07-26 08:04:34

【问题描述】：

我已经实现了一个对基因表达数据集执行批量校正的 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 代码。）

【问题讨论】：

首先，这是一个非常具体的生物信息学/生物导体问题；尽我所能与您的问题相关，我不确定 SO 是否是此类问题的正确论坛。其次，您应该始终清楚地说明您正在使用哪些非基础 R 库。例如，read.celfiles 来自oligo； rma 来自 affy。请在您的代码中包含相关调用，即library(oligo); library(affy)。

我建议在Bioconductor support site问这个问题。

第三，您链接到的 GEO 文件非常大（分别为 1.8 Gb 和 6.0 Gb）。我认为您不会很幸运地找到愿意下载近 8 Gb 数据来重现您的问题的人。有没有办法让您共享数据的最小子集？最后：ExpressionSet 是 rma 返回的 post 探测摘要。询问如何将 GeneFeatureSet 转换为 ExpressionSet 以在 rma 中使用是没有意义的。 RMA 需要*FeatureSet（强调“功能”）。

我快速浏览了一下：affy 和 oligo 都提供了 rma 函数； oligo::rma 接受 ExpressionFeatureSet 和 GeneFeatureSet。确保不要同时加载affy 和oligo；或将rma 明确指定为oligo::rma。

@MauritsEvers 感谢您的 cmets。我的问题实际上被问得很糟糕；我完全重写了它并添加了一段应该在你的电脑上运行的代码。如果你能再次帮助我，请告诉我。谢谢！

【参考方案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")

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

【讨论】：

感谢@Paggles01 的回答。我应该在我的脚本中的哪里插入你的代码？那么我应该如何处理gse3_expset 变量呢？

从 CEL 文件创建 GeneFeatureSet 对象（在我的示例中为 gse3_efs）后插入它。我假设如果您正在分析您想要创建一个ExpressionSet 的数据，这就是本示例中的gse3_expset。如果您想要的只是 rma 数据，那么您可以停止使用 basicRMA() 创建 exprsData。例如，您可以将 rma 数据作为.RData 文件保存到包含芯片数据的文件夹中。如果解决方案在您的代码中仍然不起作用，请更新。

嗨@Paggles01，非常感谢！它显然解决了:-) 当我使用 combine() 函数将rma() 的输出和你的gse3_expset 合并在一起时，它也产生了这些警告：Warning messages: 1: In alleq(levels(x[[nm]]), levels(y[[nm]])) : Lengths (2, 1) differ (string compare on first 1)1 string mismatch 2: data frame column 'exprs' levels not all.equal 3: In alleq(levels(x[[nm]]), levels(y[[nm]])) : Lengths (2, 1) differ (string compare on first 1)1 string mismatch 4: data frame column 'dates' levels not all.equal 你知道如何删除它们吗？谢谢！

一个警告就是——一个警告。但它提出了一个问题，即为什么要组合rma()（我假设形成另一个 GEO 集？）和gse3_expset 的输出。这两组的 CEL 文件来自两个不同的芯片平台，大部分数据不能直接比较。在不深入研究设计的细节（不是真正的 SO 主题）的情况下，请小心您尝试分析的内容。此外，我并不精通 RMA 文献，但如果您要基于不同芯片对样本进行 RMA 标准化，您可能需要验证它们是否应该一起分析。