## ---- 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 <- 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)