## ---- echo = FALSE-------------------------------------------------------
fn = local({
  i = 0
  function(x) {
    i <<- i + 1
    paste('Figure ', i, ': ', x, sep = '')
  }
})

## ---- message = FALSE----------------------------------------------------
library(MEAL)
library(MultiDataSet)
library(minfiData)
library(GenomicRanges)

## ------------------------------------------------------------------------
MsetExFilt <- dropMethylationLoci(MsetEx)

## ------------------------------------------------------------------------
set.seed(0)
betas <- getBeta(MsetExFilt)
betas <- betas[sample(1:nrow(betas), 30000), ]
phenotypes <- data.frame(pData(MsetExFilt))

## ----Raw_Set-------------------------------------------------------------
set <- prepareMethylationSet(matrix = betas, phenotypes = phenotypes, 
                                    annotation = "IlluminaHumanMethylation450kanno.ilmn12.hg19")
set

## ------------------------------------------------------------------------
pData(set)
summary(pData(set))

## ----Probe_Analysis------------------------------------------------------
mod <- model.matrix(~as.factor(status), data = pData(set))
proberes <- DAProbe(set = set, model = mod, method = "ls", coefficient = 2)
head(proberes)

## ----Region_Analysis, message=FALSE, warning=FALSE-----------------------
regionres <- DARegion(set = set, model = mod, coefficient = 2)
names(regionres)
head(regionres$bumphunter)
head(regionres$blockFinder)
head(regionres$DMRcate)

## ----Pipeline,  warning=FALSE--------------------------------------------
res <- DAPipeline(set = set, variable_names = "status", variable_types = "categorical", 
                      probe_method = "ls", verbose = FALSE)
res

## ----Pipeline_Equation, warning = FALSE----------------------------------
complexres <- DAPipeline(set = set, variable_names = c("status", "sex"),
                             variable_types = c("categorical", "categorical"), 
                             probe_method = "ls", num_var = 3, verbose = FALSE,
                             equation = "~ status:sex + status + sex")
complexres

## ------------------------------------------------------------------------
#Bumphunter
head(bumps(res)[[1]])
#BlockFinder
head(blocks(res)[[1]])
#DMRcate
head(dmrCate(res)[[1]])
#Probe
head(probeResults(res)[[1]])
#Region
names(regionResults(res))

## ---- eval = FALSE-------------------------------------------------------
#  exportResults(res, dir = "./results")

## ----Plot_Features, fig.cap = fn("Beta values of cg25937714 split by cancer status")----
plotFeature(res, 1)

## ----QQplot, fig.cap = fn("QQplot of the analysis.")---------------------
plotQQ(res)

## ----EWAS_plot, fig.cap = fn("Manhattan plot with the CpGs of the range highlighted")----
range <- GRanges(seqnames = Rle("chr1"), 
                                ranges = IRanges(1000000, end = 10000000))
plotEWAS(res, range = range)

## ----Range_Analysis, warning=FALSE---------------------------------------
range <- GRanges(seqnames = Rle("chr12"), 
                                ranges = IRanges(70000000, end = 80000000))
region <- DARegionAnalysis(set = set, variable_names = "status", 
                               variable_types = "categorical", 
                               covariable_names = "age", range = range, 
                               verbose = FALSE)
region

## ------------------------------------------------------------------------
#Bumphunter
head(bumps(region)[[1]])
#BlockFinder
head(blocks(region)[[1]])
#DMRcate
head(dmrCate(region)[[1]])
#Probe
head(probeResults(region)[[1]])
#Region
names(regionResults(region))

## ----Region_plot, fig.cap = fn("Plot of the differential methylation for each CpG of the region.")----
plotRegion(region)

## ----plot_RDA, fig.cap = fn("RDA plot of the region."), fig.height = 5, fig.width=8----
plotRDA(region)