## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4)

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("crossarrayEWAS")

## ----load---------------------------------------------------------------------
library(crossarrayEWAS)

## ----simulate-design----------------------------------------------------------
set.seed(20240607)

n_sub   <- 24L
visits  <- sample(2:4, n_sub, replace = TRUE)
subject <- rep(sprintf("S%02d", seq_len(n_sub)), visits)
n       <- length(subject)

# array generation is a property of the subject, not of the visit: this is the
# confounding the within-subject transform absorbs
array_type <- rep(ifelse(seq_len(n_sub) <= 14L, "450K", "EPIC"), visits)

# exposure varies between visits of the same subject, which is what makes the
# longitudinal contrast identified at all
exposure <- round(runif(n, 0, 9), 3)
covars   <- cbind(cell_fraction = rnorm(n, 0.20, 0.05),
                  smoking_score = rnorm(n, -1.0, 0.40))

c(subjects = n_sub, samples = n, on_450K = sum(array_type == "450K"))

## ----simulate-probes----------------------------------------------------------
n_clust <- 640L
sizes   <- sample(2:5, n_clust, replace = TRUE)
pos     <- integer(0)
p       <- 1000000L
for (k in seq_len(n_clust)) {
  within <- c(0L, cumsum(sample(80:600, sizes[k] - 1L, replace = TRUE)))
  pos <- c(pos, p + within)
  p   <- p + max(within) + sample(3000:20000, 1L)
}
n_probes   <- length(pos)
chrom      <- rep("chr21", n_probes)
cluster_of <- rep(seq_len(n_clust), sizes)   # the truth, for comparison below
c(probes = n_probes, planted_clusters = n_clust)

## ----simulate-mvalues---------------------------------------------------------
sub_idx <- match(subject, unique(subject))

baseline    <- rnorm(n_probes, 0, 1.5)
subject_eff <- matrix(rnorm(n_probes * n_sub, 0, 0.30), n_probes, n_sub)
array_off   <- rnorm(n_probes, 0, 0.40)
latent      <- matrix(rnorm(n_clust * n, 0, 0.25), n_clust, n)[cluster_of, ]

# planted signal: one short region, one long shallow block
beta_true <- numeric(n_probes)
region_probes <- which(cluster_of %in% 300:301)
block_probes  <- which(cluster_of %in% 420:470)
beta_true[region_probes] <- 0.30
beta_true[block_probes]  <- 0.06

M <- baseline +
  subject_eff[, sub_idx] +
  outer(array_off, as.numeric(array_type == "EPIC")) +
  latent +
  matrix(rnorm(n_probes * n, 0, 0.20), n_probes, n) +
  outer(beta_true, exposure - mean(exposure))

rownames(M) <- sprintf("cg%07d", seq_len(n_probes))
colnames(M) <- sprintf("%s_v%d", subject,
                       unlist(lapply(visits, seq_len), use.names = FALSE))
dim(M)

## ----assemble-----------------------------------------------------------------
se <- SummarizedExperiment::SummarizedExperiment(
  assays = list(M = M),
  rowRanges = GenomicRanges::GRanges(chrom, IRanges::IRanges(pos, width = 2L)),
  colData = data.frame(exposure = exposure, subject = subject, covars,
                       array_type = array_type))
se

## ----regions------------------------------------------------------------------
regions <- ewas_regions(se, exposure = "exposure", subject = "subject",
                        covars = c("cell_fraction", "smoking_score"),
                        lambda = 0.5, min_probes = 3L, min_effect = 0.05,
                        n_perm = 20L, n_boot = 8L, perm_iter = 400L)
regions[order(-abs(regions$z))][seq_len(min(5L, length(regions)))]

## ----regions-metadata---------------------------------------------------------
md <- S4Vectors::metadata(regions)
c(max_abs_z = round(md$max_abs_z, 2), lambda = md$lambda,
  clusters = md$n_clusters, df = md$fit$df)

## ----regions-recovery---------------------------------------------------------
planted <- range(pos[region_probes])
hit <- GenomicRanges::start(regions) <= planted[2] &
  planted[1] <= GenomicRanges::end(regions)
c(regions = length(regions), overlapping_the_planted_region = sum(hit))

## ----blocks-------------------------------------------------------------------
blocks <- ewas_blocks(se, exposure = "exposure", subject = "subject",
                      covars = c("cell_fraction", "smoking_score"),
                      max_gap = 1500L, rho_min = 0.20,
                      min_fit_clusters = 100L)
blocks

## ----blocks-states------------------------------------------------------------
hs <- S4Vectors::metadata(blocks)$hsmm
data.frame(state = hs$states, mu = round(hs$mu, 4),
           sigma = round(hs$sigma, 4), pi = round(hs$pi, 3))
clusters <- S4Vectors::metadata(blocks)$clusters
length(clusters)

## ----caption, include = FALSE-------------------------------------------------
# chunk options are evaluated before the chunk body, so the caption is set here
cap <- paste("Cluster-level effects along the chromosome. Shaded span:",
             "the planted block. Bars: the blocks called.")

## ----block-plot, fig.cap = cap------------------------------------------------
planted_block <- range(pos[block_probes])
cmid <- (GenomicRanges::start(clusters) + GenomicRanges::end(clusters)) / 2
plot(cmid / 1e6, clusters$effect_M, pch = 16, cex = 0.5, col = "grey40",
     xlab = "position (Mb)", ylab = "effect (M-value per exposure unit)")
rect(planted_block[1] / 1e6, par("usr")[3], planted_block[2] / 1e6,
     par("usr")[4], col = "#1b9e7733", border = NA)
abline(h = 0, col = "grey70")
if (length(blocks))
  segments(GenomicRanges::start(blocks) / 1e6, 0,
           GenomicRanges::end(blocks) / 1e6, 0, col = "#d95f02", lwd = 4)

## ----harmonise, eval = FALSE--------------------------------------------------
# library(minfi)
# 
# rg450 <- read.metharray.exp(targets = targets_450k)
# rgepic <- read.metharray.exp(targets = targets_epic)
# 
# # combineArrays() keeps the probes shared by both designs and sets the
# # resulting object's array to the older of the two
# rg <- combineArrays(rg450, rgepic, outType = "IlluminaHumanMethylation450k")
# 
# grset <- mapToGenome(ratioConvert(preprocessQuantile(rg)))
# SummarizedExperiment::colData(grset)$subject <- targets$Subject_ID
# SummarizedExperiment::colData(grset)$exposure <- targets$days_on_drug / 100
# 
# regions <- ewas_regions(grset, "exposure", "subject")

## ----real-data----------------------------------------------------------------
path <- system.file("extdata", "GSE237561_chr21_subset.rds",
                    package = "crossarrayEWAS")
grset <- readRDS(path)
grset
table(SummarizedExperiment::colData(grset)$array_type)

## ----real-regions-------------------------------------------------------------
real <- ewas_regions(grset, exposure = "exposure", subject = "subject",
                     lambda = 0.5, min_probes = 3L, min_effect = 0.05)
length(real)
if (length(real))
  real[order(-abs(real$z))][seq_len(min(3L, length(real)))]

## ----building-blocks----------------------------------------------------------
design <- stats::model.matrix(
  ~ subject + exposure + cell_fraction + smoking_score,
  data = as.data.frame(SummarizedExperiment::colData(se)))
resid <- t(stats::lm.fit(design, t(M))$residuals)

cl <- comethylation_clusters(chrom, pos, resid, max_gap = 1000L,
                             rho_min = 0.30)
c(clusters = max(cl) + 1L, clustered = sum(cl >= 0), singletons = sum(cl < 0))

## ----blocks-manual------------------------------------------------------------
fit <- fit_within_se(se, "exposure", "subject",
                     covars = c("cell_fraction", "smoking_score"))
w <- 1 / fit$se^2
lam <- distance_penalty(pos, lam0 = 0.5, decay_bp = 1000, w = w)
theta <- tv_denoise(fit$beta, w, lam, n_iter = 3000L)
obs <- call_regions(chrom, pos, fit$beta, fit$se, theta, min_probes = 3L,
                    min_effect = 0.05, cluster = ifelse(cl >= 0, cl, -seq_along(cl)))
identical(nrow(obs$table), length(regions))

## ----permutation--------------------------------------------------------------
perm <- within_subject_permutation(exposure, subject)
identical(sort(exposure[subject == "S01"]), sort(perm[subject == "S01"]))

fit_perm <- fit_within(M, perm, subject, covars)
c(observed = round(max(abs(fit$z)), 2),
  permuted = round(max(abs(fit_perm$z)), 2))

## ----perm-table---------------------------------------------------------------
GenomicRanges::mcols(regions)[order(-abs(regions$z)),
                              c("z", "p_fwer", "p_fwer_mcse")][1:3, ]
S4Vectors::metadata(regions)$perm_scheme

## ----stability----------------------------------------------------------------
GenomicRanges::mcols(regions)[order(-regions$stability, -abs(regions$z)),
                              c("n_probes", "z", "stability")][1:3, ]

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