## -----------------------------------------------------------------------------
#| label: init
#| include: false
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
options(width = 80)
library(lattice)


## -----------------------------------------------------------------------------
#| label: load
suppressPackageStartupMessages(library(methylumi))
samps <- read.table(system.file("extdata/samples.txt", package = "methylumi"),
                    sep = "\t", header = TRUE)
mldat <- methylumiR(
  system.file("extdata/exampledata.samples.txt", package = "methylumi"),
  qcfile = system.file("extdata/exampledata.controls.txt", package = "methylumi"),
  sampleDescriptions = samps
)


## -----------------------------------------------------------------------------
mldat


## -----------------------------------------------------------------------------
getAssayDataNameSubstitutions()


## -----------------------------------------------------------------------------
#| label: mds
md <- cmdscale(dist(t(exprs(mldat)[fData(mldat)$CHROMOSOME == "X", ])), 2)
plot(md,
     pch = c("F", "M")[pData(mldat)$Gender],
     col = c("red", "blue")[pData(mldat)$Gender])


## -----------------------------------------------------------------------------
#| label: pvals
avgPval <- colMeans(pvals(mldat))
par(las = 2)
barplot(avgPval, ylab = "Average P-Value")


## -----------------------------------------------------------------------------
controlTypes(mldat)


## -----------------------------------------------------------------------------
#| label: qcplot
qcplot(mldat, "FIRST HYBRIDIZATION")


## -----------------------------------------------------------------------------
#| label: normalize
toKeep <- (avgPval < 0.05)
pData(mldat)$Gender[9] <- "F"
mldat.norm <- normalizeMethyLumiSet(mldat[, toKeep])


## -----------------------------------------------------------------------------
#| label: limma
library(limma)
dm <- model.matrix(~ 1 + Gender, data = pData(mldat.norm))
colnames(dm)
fit1 <- lmFit(exprs(mldat.norm), dm)
fit2 <- eBayes(fit1)
tt <- topTable(fit2, coef = 2,
               genelist = fData(mldat.norm)[, c("SYMBOL", "CHROMOSOME")],
               number = 1536)
x <- aggregate(tt$adj.P.Val, by = list(tt$CHROMOSOME), median)
colnames(x) <- c("Chromosome", "Median adjusted P-value")


## -----------------------------------------------------------------------------
#| label: tbl-chromosomepvals
#| tbl-cap: "The median adjusted p-value for each chromosome, showing that the
#|   X chromosome is highly significantly different between males and females."
knitr::kable(x, digits = 6)


## -----------------------------------------------------------------------------
#| label: fig-genderProbesByChrom
#| fig-cap: "Probes differentially methylated, plotted by chromosome. Note that
#|   the p-values plotted here are based on a linear model. Since the underlying
#|   data are not normally distributed, the p-values representing the outcomes
#|   of the linear models are not exact."
print(xyplot(-log10(adj.P.Val) ~ CHROMOSOME, tt,
             ylab = "-log10(Adjusted P-value)",
             main = "P-values for probes\ndistinguishing males from females"))


## -----------------------------------------------------------------------------
sessionInfo()