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

## ----Load_Packages, message = FALSE--------------------------------------
library(MEAL)
library(MultiDataSet)
library(MEALData)
library(GenomicRanges)

## ----Methylation_Data----------------------------------------------------
betavals[1:5, 1:5]
dim(betavals)
summary(pheno)

## ----Expression_Data-----------------------------------------------------
eset
summary(pData(eset))

## ----Snps_Data-----------------------------------------------------------
str(snps)

## ------------------------------------------------------------------------
snps$genotypes[1:5, 1:5]

## ------------------------------------------------------------------------
head(snps$map)

## ----Meth_Object, message = FALSE----------------------------------------
methSet <- prepareMethylationSet(matrix = betavals, phenotypes = pheno)
table(fData(methSet)$chr)

## ----Filtering_Meth------------------------------------------------------
annot <- fData(methSet)
autosomiccpgs <- rownames(annot)[!annot$chr %in% c("chrX", "chrY")]
methSet <- methSet[autosomiccpgs, ]
methSet

## ----Remove NAs----------------------------------------------------------
methSet <- methSet[!is.na(fData(methSet)$position), ]
methSet

## ----Meth_Analysis-------------------------------------------------------
methRes <- DAPipeline(set = methSet, variable_names = "inv", probe_method = "ls")
methRes

## ----Meth_Manhattan, fig.cap = fn("Manhattan plot of methylation results. Some differentially methylated probes are detected in the 17q21.31 region.")----
plotEWAS(methRes)

## ----Meth_QQplot, fig.cap = fn("QQplot of methylation analysis. All the p-values are on the theoretical line but the three probes differentially methylated.")----
plotQQ(methRes)

## ----Exp show------------------------------------------------------------
eset
fvarLabels(eset)

## ----Exp Analysis--------------------------------------------------------
colnames(fData(eset))[colnames(fData(eset)) == "chr"] <- "chromosome"
expRes <- DAPipeline(eset, variable_names = "inv", probe_method = "ls")
expRes

## ----Expr_Manhattan, fig.cap = fn("Manhattan plot of expression results. No probe is differentially expressed after multiple testing correction.")----
plotEWAS(expRes)

## ----ExprQQplot, fig.cap = fn("QQplot of expression results. Most of the p-values are below the theoretical line so there is a great deflation.")----
plotQQ(expRes)

## ----Exp Analysis SVA----------------------------------------------------
expResSVA <- DAPipeline(eset, variable_names = "inv", probe_method = "ls", sva = TRUE)
expResSVA

## ----ExprManhattanSVA, fig.cap = fn("Manhattan plot of expression results using SVA. Again, no probe is differentially expressed after multiple testing correction.")----
plotEWAS(expResSVA)

## ----ExprQQSVA, fig.cap = fn("QQ plot of expression results using SVA.")----
plotQQ(expResSVA)

## ----ExprVolcanoSVA, fig.cap = fn("Volcano of expression results using SVA. All the probes are inside the non significant region.")----
plotVolcano(expResSVA)

## ------------------------------------------------------------------------
head(probeResults(expResSVA)[[1]])

## ----Exprs Region Analysis-----------------------------------------------
range <- GRanges(seqnames = Rle("chr17"),
                 ranges = IRanges(43000000, end = 45000000))
exprsRegion <- DARegionAnalysis(set = eset, range = range, variable_names = "inv", probe_method = "ls")
exprsRegion

## ----ExprsRDA, fig.cap = fn("RDA plot for the range analysis for expression data. The different genotypes are very closed but separated."), fig.width=8, fig.height=5----
plotRDA(exprsRegion)

## ----RDA hits expression-------------------------------------------------
rdahits <- topRDAhits(exprsRegion)
rdahits

## ----New Multi Meth Exp--------------------------------------------------
multi2 <- new("MultiDataSet")
multi2 <- add_genexp(multi2, eset)
multi2 <- add_methy(multi2, methSet)

## ----Corr Meth Exp-------------------------------------------------------
topcpgs <- feats(methRes)
methExprs <- correlationMethExprs(multi2, sel_cpgs = topcpgs)
head(methExprs)