---
title: "Design-aware region and block detection for cross-array longitudinal EWAS"
author: "Katarzyna Kamieniecka, Samantha McLean, Krzysztof Poterlowicz"
date: "`r Sys.Date()`"
package: crossarrayEWAS
output:
  BiocStyle::html_document:
    toc_float: true
vignette: >
  %\VignetteIndexEntry{Design-aware region and block detection for cross-array longitudinal EWAS}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# What this package is for

A cohort assembled over a decade is usually assembled over two array
generations as well: earlier visits on HumanMethylation450, later ones on
MethylationEPIC. Array type is then not a nuisance variable that happens to be
present, it is a variable that moves with calendar time and therefore with
follow-up, treatment duration and often with the exposure itself. Regressing
out an array indicator does not repair that, because the indicator and the
contrast of interest are partly the same column.

`crossarrayEWAS` takes the other route: the probe-level effect is identified
**within subject**. Anything constant within a subject's visit series is
absorbed by the transformation rather than modelled: subject identity, and in
a cohort where each subject sits on one array generation, array type and chip
as well. What is left is estimated once, grouped by observed
co-methylation, and turned into regions and long-range blocks by two
estimators that place boundaries rather than threshold a smooth curve.

The design this was written for is a longitudinal cohort, sampled repeatedly
within subject over a recruitment window long enough that the array generation
changed mid-study: earlier subjects assayed on 450K, later subjects on EPIC,
and no subject spanning the two. The validation cohort is GSE237561, a
clozapine treatment series of 126 arrays from 38 subjects, 92 on 450K and 34
on EPIC. What makes that split non-exchangeable with exposure is that one
variable decides both: subjects recruited later sit on EPIC and also have
longer treatment duration, so an array indicator and the contrast of interest
are partly the same column. In GSE237561 neither a subject nor a chip spans
the two generations, so cohort is nested inside array generation rather than
crossed with it, and a batch correction taken on cohort is not identifiable
there.

# Installation

```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("crossarrayEWAS")
```

```{r load}
library(crossarrayEWAS)
```

# The two functions most users need

`ewas_regions()` and `ewas_blocks()` take a `SummarizedExperiment` — including
a `GenomicRatioSet` from `r BiocStyle::Biocpkg("minfi")`, which extends one —
and return a `GRanges`. Each runs the analysis the pipeline's corresponding
stage runs: assay and design extraction, within-subject probe effects,
co-methylation clustering, smoothness chosen by held-out subjects, region or
block calling, and, when asked for, a permutation null and subject-level
selection stability. Everything else exported here is a building block one of
them assembles, and the later sections take those apart.

# A simulated two-array cohort

No study data ships with the package, so this vignette simulates a cohort with
the structure the estimators are designed for: subjects followed over an
unbalanced number of visits, a continuous exposure that varies **within**
subject, two array generations nested in subject, probes spaced in
co-methylated clusters, and two planted signals: a short region and a long
low-amplitude block.

```{r 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"))
```

Probe positions are drawn as tight clusters separated by wide gaps, which is
what `comethylation_clusters()` looks for and roughly what a methylation array
gives you.

```{r 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)
```

M-values are assembled from the terms the model assumes: a probe baseline, a
subject intercept, an array offset, a per-cluster latent component that makes
neighbouring probes co-methylated, independent noise, and the exposure effect.

```{r 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)
```

Assembled into a `SummarizedExperiment`, this is the shape the package works
on: an M-value assay, probe coordinates in `rowRanges()`, and the design in
`colData()`. A `GenomicRatioSet` carries the same three things, which is why
one can be passed straight in.

```{r 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

```{r 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)))]
```

The ranges come back in genomic order with the per-region statistics as
metadata columns, and the run's own summary in `metadata()`: the observed
maximum |z| the null is compared against, the penalty used, the number of
co-methylation clusters, and the degrees of freedom the within-subject
transform left.

```{r 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)
```

The planted region is recovered and the rest of the genome is not:

```{r 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))
```

Three arguments carry most of the scientific content of a run.

`min_effect` is on the M-value scale and is a choice about what counts as an
effect, not a tuning knob. In the clozapine analysis it was 0.05 M-value units
per 100 days of treatment: a segment has to move at least that much to be
called, whatever its z. At mid-range methylation that is about 0.9 percentage
points of beta per 100 days, so roughly 3 points over a year of treatment and 6
over two, falling to about 2 and 4 points at beta 0.2 or 0.8. The floor is set
there because a slope below it accumulates to less than the change a single
array can resolve over a realistic follow-up window, so a segment that clears
it is reportable on the beta scale as well as on the z scale.

`lambda` is fixed above to keep the vignette quick. Left at its default it is a
vector, and the penalty is then chosen by cross-validation over folds of whole
**subjects**, with the score table returned in `metadata()$lambda_cv`.
Held-out subjects rather than held-out samples are required: a subject's
visits are the unit the within-subject effect is estimated from, so splitting
them across train and test would leak the quantity being validated.

`var_method` chooses how the per-probe residual variance is moderated before
the z statistics are formed. The default `"limma"` calls `limma::squeezeVar()`,
the reference implementation of that model. `"mom"` is the method-of-moments
fit ported from the Python twin of this core, and it exists so the two
implementations can be compared arithmetically on identical input. The stage
drivers default to `"limma"` and the published analysis was run with it, so
report with the default and name `"mom"` in the methods if a reported result
uses it.

# Long-range blocks

Blocks are a coarser statement: co-methylated probes are collapsed to one
profile per cluster, effects are refitted at cluster level, and a three-state
hidden semi-Markov model runs along the chromosome with transition
probabilities that decay with genomic distance.

```{r 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
```

The fitted states, the collapsed clusters and their posteriors are in
`metadata()`. `hsmm$neutral` is the column of the posterior holding the state
pinned at `mu = 0`, and it is **not** always the middle one: if every cluster
effect falls on one side of zero that column sorts to an end and only one
direction can be called, which `hsmm$one_sided` reports.

```{r 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)
```

```{r 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.")
```

```{r 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)
```

A region and a block are not two attempts at the same quantity. Regions are
called on every harmonised probe at the resolution of a few probes, median 84
bp over 3 probes in the published run, and locate a focal effect. Blocks are
called on open-sea probes only, at cluster resolution, and in the same run
returned 38 blocks of median width 134 kb spanning a median of 9 clusters,
against 207 regions. A block is therefore evidence of a shallow shift across a
domain that is largely unmeasured between clusters, and a region is evidence
of a local one; the two can coexist, and neither subsumes the other.

When they disagree, they are reported as what they are: two independent lines
of evidence, each on its own scale. A block with no region inside it is not
demoted, and a region inside a neutral block is not discounted; the region
table and the block table are read side by side rather than one filtering the
other.

# From 450K and EPIC data

`ewas_regions()` needs one object spanning both array generations. Building it
is `r BiocStyle::Biocpkg("minfi")`'s job, and it is three calls: combine the
two `RGChannelSet`s on the probes they share, convert to ratios, and map to
the genome. The pipeline's stage 01 does exactly this, with quality control
and probe filtering around it, and writes the matrix `read_f64()` reads back.

```{r 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")
```

Two things to check before the fit, because neither is detectable afterwards.
The exposure has to vary between visits of the same subject, or there is
nothing within-subject to estimate and `ewas_regions()` stops. And a covariate
that is constant within subject — sex, genotype, array type in a cohort where
no subject spans both generations — is absorbed by the transform and cannot be
passed as a covariate; this is the point of the design rather than a
limitation of it.

# A real cohort

A 2,000-probe subset of GSE237561 ships with the package as a
`GenomicRatioSet`, with the script that built it in `inst/scripts/`. It is the
clozapine treatment series the method was developed on: 126 arrays from 38
subjects, 92 on 450K and 34 on EPIC, with the exposure being days on treatment
and no subject spanning the two array generations.

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

```{r 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)))]
```

The subset is a thousandth of the harmonised panel and is here so the
interface can be exercised on real data, not to reproduce the published
result: regions called on 2,000 probes of one chromosome are not the regions
called on 422,898 probes genome-wide, and the permutation null is not run here
because its cost is what the pipeline exists for.

# Building blocks

The two functions above are assembled from exported pieces, and an analysis
that needs to depart from them can use those directly. The clustering step is
representative: probes are grouped by their **residual** co-methylation, after
subject means and the design have been projected out.

```{r 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))
```

Cluster labels are 0-based and `-1` marks a probe that joined nothing, the
convention the Python implementation of the same estimator uses.

```{r 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))
```

`as_granges()` converts either caller's table, keeping `max_abs_z` in
`metadata()`; `sort_probes()` and `probe_coords()` are the object accessors the
high-level functions use, and are exported because an analysis that calls
`fit_within()` directly needs the same ordering guarantee.

# Resampling subjects, not samples

Both null constructions resample at the subject level, so a subject's visit
series is never split across a fold or a permutation.
`within_subject_permutation()` permutes exposure *inside* each subject, which
leaves the between-subject distribution intact and destroys only the
within-subject contrast, which is the thing being tested.

```{r 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))
```

`n_perm` runs that null inside `ewas_regions()`, and `p_fwer` is the fraction
of null replicates whose maximum |z| reaches a region's own, so it is a
family-wise statement across the genome rather than a pointwise one.
`p_fwer_mcse` is the Monte Carlo standard error of that fraction, reported
because the decision at 0.05 is a step function of the null tail: a p-value
whose error straddles the threshold is evidence about `n_perm`, not about the
region. `perm_scheme = "naive"` shuffles across subjects instead, the scheme
`r BiocStyle::Biocpkg("bumphunter")` uses, and exists to quantify its
mis-calibration in a longitudinal design rather than to be used.

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

`n_boot` adds selection stability: regions are re-called on repeated
half-samples of subjects, and `stability` is the fraction of subsamples in
which an overlapping region was called again. Boundaries move between
subsamples, so overlap rather than matching endpoints is the criterion.
`stability_selection()` is exported for analyses that want the same resampling
around a different estimator; the callback is handed subject identifiers, not
column indices, and returns the keys it re-discovered.

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

The published run used 200 permutations and 50 half-samples; the small numbers
here keep the vignette quick. Selection frequency is carried beside z rather
than used as a filter, so a region that survives the permutation null but is
unstable across subsamples stays visible and is read as exactly that. The
threshold reported on is 0.6, conventional in the stability-selection
literature: regions selected in at least 60% of subsamples are described as
stable, and the column annotates the table rather than removing rows from it.

# Comparison with related packages

What `crossarrayEWAS` adds is a design, not a new smoother. Three differences
matter when choosing between it and the established tools.

*Where the design sits.* `r BiocStyle::Biocpkg("bumphunter")`,
`r BiocStyle::Biocpkg("DMRcate")`, `r BiocStyle::Biocpkg("ENmix")`
(`combp()`, `ipdmr()`) and `r BiocStyle::Biocpkg("coMethDMR")` all take a
per-probe statistic from a model the user fits, and a repeated-measures design
has to be handled in that model: `limma::duplicateCorrelation()` with a
subject block, or `r BiocStyle::Biocpkg("variancePartition")`'s `dream()` with
a random intercept. Both estimate a between-subject component and a
within-subject one and then combine them, which is exactly what a cohort with
array generation nested in subject cannot support: the between-subject
contrast is partly the array contrast. Here the probe-level effect is
differenced within subject before any region calling, so the confounded
component is absorbed instead of modelled, and no array indicator is needed.

*How the null is built.* `bumphunter` permutes the outcome labels across
samples, which in a longitudinal cohort breaks the within-subject pairing and
the between-subject structure at once. `within_subject_permutation()` permutes
exposure among the visits of one subject, leaving subject, array, chip, visit
count and the within-subject covariance intact; `perm_scheme = "naive"`
reproduces the other scheme so the difference can be measured on the data at
hand rather than assumed.

*How boundaries are placed.* `bumphunter` and `DMRcate` smooth the statistic
and threshold the result; `combp()` and `ipdmr()` combine spatially
autocorrelated p-values; `coMethDMR` tests pre-defined co-methylated regions.
`ewas_regions()` solves a total-variation problem whose penalty decays with
the gap between probes, so the solution is piecewise constant and a boundary
is a jump in it: a wide gap is cheap to cross, a tight pair is not, and no
kernel width has to be chosen. For blocks, `minfi::blockFinder()` collapses
probes into fixed-width bins and smooths; `ewas_blocks()` collapses into
clusters defined by observed within-subject co-methylation and calls blocks
from a three-state model whose transition probabilities decay with distance,
so a block's extent is a model statement with a posterior rather than a
thresholded curve.

Where the design is not longitudinal, or a single array generation is in play,
the established tools are the better choice: they are mature, widely used, and
their smoothing assumptions are well characterised. This package is for the
case they are not built for.

# Relationship to the pipeline and the Galaxy suite

These functions are the numerical core of a Nextflow pipeline whose earlier
stages harmonise 450K and EPIC intensities, and of a Galaxy tool suite that
wraps the same stage drivers. The pipeline is
[kkamieniecka/ewas-crossarray-harmonise](https://github.com/kkamieniecka/ewas-crossarray-harmonise),
which generates this package from its numerical core, and the wrappers are
published as a Tool Shed-installable suite at
[bioinfbrad/ewas-crossarray-harmonise](https://github.com/bioinfbrad/ewas-crossarray-harmonise);
both extend the EWASGalaxy tool suite (Murat *et al.*, *bioRxiv* 2019,
doi:10.1101/553784).

# Session info

```{r session}
sessionInfo()
```