Как применить нормализацию rma () к уникальному файлу CEL? - PullRequest
Новая
       28

Как применить нормализацию rma () к уникальному файлу CEL?

0 голосов
/ 05 марта 2019

Я реализовал R-скрипт, который выполняет пакетную коррекцию набора данных генной экспрессии.Чтобы выполнить пакетную коррекцию, мне сначала нужно нормализовать данные в каждом файле CEL с помощью функции Affy rma () в Bioconductor.

Если я запускаю ее на наборе данных GSE59867 получено от GEO, все работает.Я определяю пакет как дату сбора данных: я помещаю все файлы CEL с одинаковой датой в определенную папку, а затем рассматриваю эту дату / папку как конкретный пакет.В наборе данных GSE59867 пакет / папка содержит только 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 объектами).

( EDIT: Я тщательно отредактировал вопрос и добавил фрагмент кода R, который вы сможете запустить на своем компьютере.)

1 Ответ

0 голосов
/ 06 марта 2019

Хм.Это загадочная проблема.функция oligo::rma() может содержать ошибки для класса GeneFeatureSet с отдельными выборками.Я заставил его работать с одним образцом, используя функции более низкого уровня, но это означает, что мне также пришлось создавать набор выражений с нуля, указав слоты: 

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

Надеюсь, что приведенный выше подход будет работать вваш код.

...