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