inst/tests/prepareData.R

###############################
# Try preprocessing cel files #
###############################

# Load packages
#require(BiocInstaller)
#biocLite("oligo")

# require(ff) # it will change the exprs matrix format
suppressMessages(require(foreach))
suppressMessages(require(doMC))
suppressMessages(require(oligo))
# require(parallel) # this package has already loaded by doMC package
# celFiles <- list.celfiles(data.dir, listGzipped = TRUE, full.names = TRUE)
# rawData  <- read.celfiles(celFiles)
# rmaRes   <- rma(rawData)
# sampleNames(rmaRes) <- sub("(^GSM.*)\\.CEL.gz","\\1",sampleNames(rmaRes))
#
# rawData1 <- affy::ReadAffy(filenames=celFiles)
# rmaRes2  <- affy::rma(rawData1)
# hist(rmaRes2)
# gcrmaRes2 <- gcrma::gcrma(rawData1)
# hist(gcrmaRes2)

Ncores <- parallel::detectCores()
registerDoMC(cores = Ncores)
data.dir <- "./inst/tests/iGCC/cel_files/GSE7670/"
rmaRes <- affyPreprocess.oligo(data.dir)
t1 <- exprs(rmaRes)

rmaRes2 <- affyPreprocess.affy(data.dir)
t2 <- exprs(rmaRes2)
gcrmaRes2 <- affyPreprocess.affy(data.dir, method="gcRMA")
t3 <- exprs(gcrmaRes2)

# compare
hist(rmaRes)
hist(rmaRes2)
hist(gcrmaRes2)

# get directory of cel files
# GSE72094 has no annotation package, remove it
dir_list <- dir("./inst/tests/iGCC/cel_files", pattern = "^GSE[0-9]{4,5}$",
                full.names = TRUE)
dir_list <- setdiff(dir_list, "./inst/tests/iGCC/cel_files/GSE72094")
Studynames <- basename(dir_list)
Study_Res  <- list()
for (i in 1:length(dir_list)){
    Study_Res[[Studynames[i]]] <- affyPreprocess.oligo(dir_list[i])
}

# make annotation package for GSE72094
# library(pdInfoBuilder)
# cdf <- "./inst/tests/iGCC/makePd/GPL15048_HuRSTA_2a520709.CDF"
# cel <- oligoClasses::list.celfiles("./inst/tests/iGCC/cel_files/GSE72094", listGzipped = TRUE, full.names = TRUE)[]
# tab <- "./inst/tests/iGCC/makePd/GPL15048_probe_sequences.txt"
#
# new("AffyExpressionPDInfoPkgSeed",
#     cdfFile = cdf, celFile = cel,
#     tabSeqFile = tab, author = "ShixiangWang",
#     email = "wangshx@shanghaitech.edu.cn",
#     biocViews = "AnnotationData",
#     genomebuild = "NCBI Build 36",
#     organism = "Human", species = "Homo Sapiens",
#     url = "http://www.biostat.jhsph.edu/~bcarvalh")
ShixiangWang/iProfile documentation built on May 11, 2019, 6:25 p.m.