normalizeAndCreateTables.R

library(minfi)
source("config.R")

targetDir = paste(base,"normalized",sep="/")
if(!file.exists(targetDir)) { dir.create(file.path(targetDir)) }

rgSet = read.metharray(basenames = unique(sampleInfo$Basename))

detM = detectionP(rgSet)
failed_detection = detM > det_p
failed_probes = which(rowMeans(failed_detection) > sample_threshold)

if(exists("qc_file")) {
	qcReport(rgSet = rgSet, sampNames = sampleInfo$Sample_Name, sampGroups = sampleInfo$Sample_Group, pdf = qc_file)
}

if(normalization == "fn") {
	mSet = preprocessFunnorm(rgSet, bgCorr = bg_noob, dyeCorr = dye_corr)
	write.table(getM(mSet)[-failed_probes,], file = paste(targetDir, "mValues_fn.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
	write.table(getBeta(mSet)[-failed_probes,], file = paste(targetDir, "betaValues_fn.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
}

if(normalization == "swan") {
	mSet = preprocessSWAN(rgSet)
	write.table(getM(mSet)[-failed_probes,], file = paste(targetDir, "mValues_swan.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
	write.table(getBeta(mSet)[-failed_probes,], file = paste(targetDir, "betaValues_swan.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
}

if(normalization == "raw") {
	mSet = preprocessRaw(rgSet)
	write.table(getM(mSet)[-failed_probes,], file = paste(targetDir, "mValues_raw.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
	write.table(getBeta(mSet)[-failed_probes,], file = paste(targetDir, "betaValues_raw.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
}

if(normalization == "noob") {
	mSet = preprocessNoob(rgSet, dyeCorr = dye_corr)
	write.table(getM(mSet)[-failed_probes,], file = paste(targetDir, "mValues_noob.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
	write.table(getBeta(mSet)[-failed_probes,], file = paste(targetDir, "betaValues_noob.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
}

if(normalization == "illumina") {
	mSet = preprocessIllumina(rgSet)
	write.table(getM(mSet)[-failed_probes,], file = paste(targetDir, "mValues_illumina.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
	write.table(getBeta(mSet)[-failed_probes,], file = paste(targetDir, "betaValues_illumina.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
}

if(normalization == "quantile") {
	mSet = preprocessQuantile(rgSet, fixOutliers=fixOutliers, removeBadSamples=removeBadSamples, badSampleCutoff=badSampleCutoff)
	write.table(getM(mSet)[-failed_probes,], file = paste(targetDir, "mValues_quantile.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
	write.table(getBeta(mSet)[-failed_probes,], file = paste(targetDir, "betaValues_qunatile.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
}

if(normalization == "qn.bmiq") {
	mSet = preprocessQuantile(rgSet, fixOutliers=fixOutliers, removeBadSamples=removeBadSamples, badSampleCutoff=badSampleCutoff)
	betas = getBeta(mSet)[-failed_probes,]
	design = ifelse(rownames(betas) %in% getProbeInfo(rgSet, type="I")$Name, 1, 2)
	bnorm = lapply(1:ncol(betas), function(i) BMIQ(betas[,i],design=design)$all)
	betaValues = do.call("cbind", bnorm)
	colnames(betaValues) = colnames(betas)
	mValues = apply(betaValues, 2, logit2)
	write.table(mValues, file=paste(targetDir, "mValues_qn.bmiq.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
	write.table(betaValues, file=paste(targetDir, "betaValues_qn.bmiq.txt", sep="/"), quote = F, sep = "\t", row.names = T, col.names = T)
}