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

# The class entry points are guarded: SummarizedExperiment and GenomicRanges
# are Suggests, so the vignette still builds where they are absent.
has_bioc <- all(vapply(
  c("SummarizedExperiment", "GenomicRanges", "IRanges", "S4Vectors"),
  requireNamespace, logical(1), quietly = TRUE))
```

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

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

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

The package takes a plain matrix of M-values with probes in rows plus parallel
vectors for the design. Nothing here needs the IDATs or an array manifest: the
harmonisation stage that builds this matrix from both array generations is a
separate step in the pipeline, and `read_f64()` is what loads the matrix it
writes.

The matrix is simulated here so the vignette runs anywhere; on real data it
is the harmonised matrix written by stage 01 of the pipeline, read back with
`read_f64()`.

# Probe-level effects, fitted within subject

`fit_within()` does the whole transform: it demeans the design and the
M-values by subject, drops design columns with no within-subject variation,
solves once for all probes, and moderates the residual variances towards their
common trend.

```{r fit-within}
fit <- fit_within(M, exposure, subject, covars)
str(fit)
```

The reported degrees of freedom are what the within transform leaves, not the
sample count. A subject with a single visit contributes nothing to this
contrast and is carried but not informative.

```{r fit-within-recovery}
# effects are recovered where they were planted, and not elsewhere
where <- ifelse(beta_true > 0.1, "region",
                ifelse(beta_true > 0, "block", "null"))
round(tapply(abs(fit$z), where, median), 3)
```

`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: maintained, well behaved in the
degenerate cases, and already a dependency of the harmonisation stage. `"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 rule to follow is the one the pipeline follows. The stage drivers default
to `"limma"` and the published analysis was run with it, so report with the
default. Choose `"mom"` when you are reproducing the Python path, or deploying
where the Bioconductor stack is not installed, and name it in the methods if a
reported result uses it.

# Co-methylation clusters

A region should not be allowed to run across probes that are not correlated,
so probes are grouped by their **residual** co-methylation, after the exposure
effect has been projected out. This is the idiom the stage-04 driver uses, and
`pinv()` is exported so it can be reproduced here.

```{r clusters}
resid <- demean_rows_by_group(M, subject)
Xw    <- demean_by_group(cbind(exposure, covars), subject)
Xw    <- Xw[, apply(abs(Xw), 2, max) > 1e-9, drop = FALSE]
resid <- resid - (resid %*% Xw %*% pinv(crossprod(Xw))) %*% t(Xw)

cl <- comethylation_clusters(rep(21L, n_probes), pos, resid,
                             max_gap = 1000L, rho_min = 0.30)
n_cl <- if (max(cl) >= 0) max(cl) + 1L else 0L
c(clusters = n_cl, clustered_probes = sum(cl >= 0), singletons = sum(cl < 0))
```

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

# Regions

Two steps: a distance-decaying total-variation penalty, then segmentation of
the denoised track. The penalty falls off with the gap between probes, so a
wide gap is cheap to jump and a tight pair is not, and the region boundaries
are the jumps in the piecewise-constant solution rather than crossings of a
threshold.

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

cl_arg <- ifelse(cl >= 0, cl, -seq_along(cl))
obs <- call_regions(chrom, pos, fit$beta, fit$se, theta,
                    min_probes = 3L, min_effect = 0.05, cluster = cl_arg)
reg <- obs$table
c(regions = nrow(reg), max_abs_z = round(obs$max_abs_z, 2))
head(reg[order(-abs(reg$z)),
         c("chr", "start", "end", "n_probes", "effect_M", "z")], 5)
```

`min_effect` is on the M-value scale and is a scientific choice, not a tuning
knob. In the clozapine analysis `min_effect` was 0.05 on the M-value scale 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.

`lam0` is cross-validated over held-out **subjects** in the pipeline rather
than fixed as it is here. `fold_assign()` is the split it uses: deterministic
and language-independent, so the R and Python drivers agree on it.

```{r folds}
folds <- fold_assign(subject, n_folds = 5L, seed = 0)
lengths(folds)          # subjects per fold, not samples
```

`within_df()` is the cheap test for whether a fold or a subsample is estimable
before paying for the fit: it returns the residual degrees of freedom the
within transform would leave, and zero when the exposure is annihilated
altogether.

```{r within-df}
held_out <- subject %in% folds[[1]]
c(all = within_df(exposure, subject, covars),
  train = within_df(exposure[!held_out], subject[!held_out],
                    covars[!held_out, , drop = FALSE]))
```

# Working from a SummarizedExperiment

Everything above is matrices and parallel vectors, which is what the pipeline
stages pass around. If your data is already a `SummarizedExperiment` — or a
`GenomicRatioSet`, which extends one — four functions take you in and out
without unpacking it by hand: `sort_probes()`, `probe_coords()`,
`fit_within_se()` and `as_granges()`.

```{r se-in, eval = has_bioc}
se <- SummarizedExperiment::SummarizedExperiment(
  assays = list(M = M),
  rowRanges = GenomicRanges::GRanges(chrom, IRanges::IRanges(pos, width = 2L),
                                     cluster = cl),
  colData = data.frame(exposure = exposure, subject = subject, covars,
                       array_type = array_type))

se  <- sort_probes(se)          # genomic order, once, before the fit
co  <- probe_coords(se)
fse <- fit_within_se(se, "exposure", "subject",
                     covars = c("cell_fraction", "smoking_score"))
identical(round(fse$z[names(fit$z)], 10), round(fit$z, 10))
```

The fit is the same code: `fit_within_se()` only resolves the object, and the
identity above is the point of it. Two things it does that a hand-written
extraction usually forgets — it orders nothing silently, which is why
`sort_probes()` is a separate call, and it warns when the assay looks like
beta values, since the region and block callers assume the M scale.

Region and block calls come back as data frames because that is what the
pipeline writes to CSV. `as_granges()` turns either into a `GRanges`, keeping
`max_abs_z` in `metadata()` because that is the value a permutation null is
compared against.

```{r se-out, eval = has_bioc}
gr <- as_granges(obs)
gr[order(-abs(GenomicRanges::mcols(gr)$z))][1:3]
S4Vectors::metadata(gr)$max_abs_z
```

`SummarizedExperiment` and `GenomicRanges` are `Suggests`, not `Imports`: the
numerical core is base R plus limma so the pipeline can source it in a conda
environment with no Bioconductor stack, and the conda build of this package
inherits that. That is also why these are plain functions rather than S4
methods on `fit_within()` — `setMethod()` would need the generic's package at
install time.

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

`stability_selection()` refits on subject-level subsamples and reports how
often each region survives. The driver default is 100 resamples; this runs 8 to
keep the vignette quick.

The callback is handed a set of **subject identifiers**, not column indices,
and returns the keys it re-discovered. Boundaries move between subsamples, so
an observed region counts as re-discovered when the intervals overlap rather
than when the endpoints match, the same convention the stage-04 driver uses.

```{r stability}
keys <- sprintf("%s:%d-%d", reg$chr, reg$start, reg$end)

fit_subset <- function(sub) {
  sel <- subject %in% sub
  if (within_df(exposure[sel], subject[sel],
                covars[sel, , drop = FALSE]) <= 0)
    return(character(0))    # this subsample cannot estimate the contrast
  f  <- fit_within(M[, sel, drop = FALSE], exposure[sel], subject[sel],
                   covars[sel, , drop = FALSE])
  th <- tv_denoise(f$beta, 1 / f$se^2, lam, n_iter = 500L)
  r  <- call_regions(chrom, pos, f$beta, f$se, th, min_probes = 3L,
                     min_effect = 0.05, cluster = cl_arg)$table
  if (!nrow(r)) return(character(0))
  keys[vapply(seq_len(nrow(reg)), function(i)
    any(r$chr == reg$chr[i] & r$start <= reg$end[i] & reg$start[i] <= r$end),
    logical(1))]
}

set.seed(1L)
ss <- stability_selection(fit_subset, subject, n_boot = 8L, frac = 0.5)
reg$stability <- ifelse(keys %in% names(ss$freq), ss$freq[keys], 0)
reg$stability[is.na(reg$stability)] <- 0
head(reg[order(-reg$stability, -abs(reg$z)),
         c("chr", "start", "end", "n_probes", "z", "stability")], 5)
```

The published run resampled 50 subject subsamples at half the cohort; the
eight used above only keep the example fast. Selection frequency is carried as
a column 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, which is conventional in
the stability-selection literature: regions selected in at least 60% of
subsamples are described as stable, and the column annotates the region table
rather than removing rows from it.

# Long-range blocks

Blocks are called at a coarser resolution: clusters are collapsed to their mean
profile, 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 collapse}
runs <- split(which(cl >= 0), cl[cl >= 0])

Mc     <- t(vapply(runs, function(rows) colMeans(M[rows, , drop = FALSE]),
                   numeric(ncol(M))))
cstart <- vapply(runs, function(rows) min(pos[rows]), numeric(1))
cend   <- vapply(runs, function(rows) max(pos[rows]), numeric(1))
cmid   <- (cstart + cend) / 2
cchr   <- rep("chr21", nrow(Mc))

o <- order(cmid)
Mc <- Mc[o, , drop = FALSE]
cstart <- cstart[o]; cend <- cend[o]; cmid <- cmid[o]
dim(Mc)
```

```{r hsmm}
fit_cl <- fit_within(Mc, exposure, subject, covars)
hs <- fit_block_hsmm(fit_cl$beta, fit_cl$se, cmid, cchr,
                     length_scale = 250000, seed = 1L)
labs <- state_labels(hs$neutral)
data.frame(state = labs, mu = round(hs$mu, 4), sigma = round(hs$sigma, 4),
           pi = round(hs$pi, 3))
```

`hs$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 then be called. Pass it to `call_blocks()` rather than assuming the middle
column, or every block comes back mislabelled.

```{r blocks}
bl <- call_blocks(cchr, cstart, cend, hs$posterior, min_post = 0.80,
                  min_clusters = 3L, neutral = hs$neutral)
if (nrow(bl)) bl[order(-bl$width), ][seq_len(min(5L, nrow(bl))), ] else bl
```

```{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 <- range(pos[block_probes])
plot(cmid / 1e6, fit_cl$beta, pch = 16, cex = 0.5, col = "grey40",
     xlab = "position (Mb)", ylab = "effect (M-value per exposure unit)")
rect(planted[1] / 1e6, par("usr")[3], planted[2] / 1e6, par("usr")[4],
     col = "#1b9e7733", border = NA)
abline(h = 0, col = "grey70")
if (nrow(bl))
  segments(bl$start / 1e6, 0, bl$end / 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.

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