### R code from vignette source 'ENmix.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: options
###################################################
options(width=70)


###################################################
### code chunk number 2: UnevaluatedCode (eval = FALSE)
###################################################
## library(ENmix)
## require(minfi)
## #see minfi user guide for the format of sample_sheet.txt file
## targets <- read.table("./sample_sheet.txt", header=T)
## rgSet <- read.450k.exp( targets = targets, extended = TRUE)
## # or read in all idat files under a directory
## rgSet <- read.450k.exp(base = "path_to_directory_idat_files", 
## targets = NULL, extended = TRUE, recursive=TRUE)


###################################################
### code chunk number 3: UnevaluatedCode (eval = FALSE)
###################################################
## M<-matrix_for_methylated_intensity
## U<-matrix_for_unmethylated_intensity
## pheno<-as.data.frame(cbind(colnames(M), colnames(M)))
## names(pheno)<-c("Basename","filenames")
## rownames(pheno)<-pheno$Basename
## pheno<-AnnotatedDataFrame(data=pheno)
## anno<-c("IlluminaHumanMethylation450k", "ilmn12.hg19")
## names(anno)<-c("array", "annotation")
## mdat<-MethylSet(Meth = M, Unmeth = U, annotation=anno, 
## phenoData=pheno)


###################################################
### code chunk number 4: load (eval = FALSE)
###################################################
## library(ENmix)
## require(minfi)
## require(minfiData)
## sheet <- read.450k.sheet(file.path(find.package("minfiData"),
## "extdata"), pattern = "csv$")
## rgSet <- read.450k.exp(targets = sheet, extended = TRUE)


###################################################
### code chunk number 5: ctrlplot (eval = FALSE)
###################################################
## plotCtrl(rgSet)


###################################################
### code chunk number 6: filter (eval = FALSE)
###################################################
## qc<-QCinfo(rgSet)
## mraw <- preprocessRaw(rgSet)
## mraw<-QCfilter(mraw, qcinfo=qc, samplethre = 0.01, CpGthre = 0.05
##     ,bisulthre=5000,plot=TRUE, outid=NULL, outCpG=NULL)


###################################################
### code chunk number 7: preprocessENmix (eval = FALSE)
###################################################
## mdat<-preprocessENmix(rgSet, bgParaEst="oob", nCores=6)
## mdat<-QCfilter(mdat, qcinfo=qc, samplethre = 0.01, CpGthre = 0.05
##     ,bisulthre=5000,plot=TRUE, outid=NULL, outCpG=NULL)


###################################################
### code chunk number 8: normalize.quantile.450k (eval = FALSE)
###################################################
## mdat<-normalize.quantile.450k(mdat, method="quantile1")


###################################################
### code chunk number 9: bmiq.mc (eval = FALSE)
###################################################
## beta<-bmiq.mc(mdat, nCores=6)


###################################################
### code chunk number 10: pcrplot (eval = FALSE)
###################################################
## cov<-data.frame(group=pData(mdat)$Sample_Group,
##     slide=factor(pData(mdat)$Slide))
## pcrplot(beta, cov, npc=6)


###################################################
### code chunk number 11: combat.mc (eval = FALSE)
###################################################
## batch<-factor(pData(mdat)$Slide)
## betaC<-ComBat.mc(beta, batch, nCores=6, mod=NULL)


###################################################
### code chunk number 12: nmode.mc (eval = FALSE)
###################################################
## nmode<- nmode.mc(beta, minN = 3, modedist=0.2, nCores = 5)


###################################################
### code chunk number 13: example (eval = FALSE)
###################################################
## library(ENmix)
## #read in raw intensity data
## sheet <- read.450k.sheet(file.path(find.package("minfiData"),
##     "extdata"), pattern = "csv$")
## rgSet <- read.450k.exp(targets = sheet, extended = TRUE)
## #Control plots
## plotCtrl(rgSet)
## #QC info
## qc<-QCinfo(rgSet)
## #background correction
## mdat<-preprocessENmix(rgSet, bgParaEst="oob", nCores=6)
## #filter out low quality samples and probes
## mdat<-QCfilter(mdat, qcinfo=qc, samplethre = 0.01, CpGthre = 0.05
##     ,bisulthre=5000,plot=TRUE, outid=NULL, outCpG=NULL)
## #Search for multimodal CpGs 
## #sample size in this example data is too small for this purpose!
## #should not use beta matrix after ComBat analysis for this purpose!
## beta<-getBeta(mdat, "Illumina")
## nmode<-nmode.mc(beta, minN = 3, modedist=0.2, nCores = 6)
## #between-array normalization
## mdat<-normalize.quantile.450k(mdat, method="quantile1")
## #probe type bias correction
## beta<-bmiq.mc(mdat, nCores=6)
## # Principal component regression analysis plot
## cov<-data.frame(group=pData(mdat)$Sample_Group,
##     slide=factor(pData(mdat)$Slide))
## pcrplot(beta, cov, npc=6)
## #batch correction 
## batch<-factor(pData(mdat)$Slide)
## betaC<-ComBat.mc(beta, batch, nCores=6, mod=NULL)


###################################################
### code chunk number 14: sessionInfo
###################################################
toLatex(sessionInfo())