library(IMA) library("minfi") library("sva") library("IlluminaHumanMethylation450kmanifest") rm(list = ls()) #####Load data from the .idat files pd <- read.table("SamplesTableFinalReport.txt", sep="\t", header=TRUE, skip=8) pd$X <- NULL # noticed that the names have to change to proper R name, as read in above rownames(pd) <- make.names(pd$"Sample.ID") colnames(pd) ## all idat files baseDir <- file.path("../data","idat") list.files(baseDir) # CRITICAL ## create "target" file manually ## which need a column "Basename" for path to the file ## make sure the filename and pheno info are correct Basename = with(pd, file.path(baseDir, paste(Sentrix.Barcode, "_", Sample.Section, sep=""))) target <- cbind(pd,Basename) RGset <- read.450k.exp(target=target, verbose=TRUE) save(RGset, file="../Rdata/RGset") ####Preprocess and Swan MSet.raw <- preprocessRaw(RGset) ##MethylSet MSet.SWAN <- preprocessSWAN(RGset, MSet.raw)##MethylSet ####Convert minfi to IMA data format load("annot.Rdata") detP <- detectionP(RGset) grpInfo <- pData(MSet.SWAN) convert.minfi.to.IMA <- function(x, annot, detectP=NULL, groupinfo=NULL){ if (!is(x, "RGChannelSet") & !is(x, "MethylSet") ) stop("needs to be a 'RGChannelSet' or 'MethylSet' object") if (is(x, "RGChannelSet") ){ bmatrix <- getBeta(preprocessRaw(x) ) detectP <- detectionP(x) } else { if(is.null(detectP) ) stop("needs to provide detectP") bmatrix <- getBeta(x) } probes <- rownames(bmatrix) extraProbes <- probes[which(!probes %in% rownames(detectP) )] if (length(extraProbes) >0){ e <- matrix(rep(NA, ncol(detectP)*length(extraProbes)), ncol=ncol(detectP) ) colnames(e) <- colnames(detectP) rownames(e) <- extraProbes detectP <- rbind(detectP, e) } detectP <- detectP[probes,] if(! ( identical(probes, rownames(detectP)) & identical(colnames(bmatrix), colnames(detectP)) ) ) stop("different probes/samples in beta and detectP matrix") annot <- annot[probes,] #subset and sort if(!identical(probes, rownames(annot) )) stop("Annotation file does not include probe info for all probes.\n") y <- new("exprmethy450", bmatrix = bmatrix, detectP = detectP, annot = as.matrix(annot) ) if(is.null(groupinfo)) { warning("Please add groupinfo\n") } else{ y@groupinfo <- groupinfo } y } IMA.swan <- convert.minfi.to.IMA(MSet.SWAN,anno,detectP = detP,groupinfo = grpInfo) IMA.raw <- convert.minfi.to.IMA(RGset,anno,detectP = detP,groupinfo = grpInfo)