Design-aware region and block detection for cross-array longitudinal EWAS
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.
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.
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"))
#> subjects samples on_450K
#> 24 66 38Probe 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.
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)
#> probes planted_clusters
#> 2225 640M-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.
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)
#> [1] 2225 66The 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.
fit <- fit_within(M, exposure, subject, covars)
str(fit)
#> List of 6
#> $ beta : Named num [1:2225] -0.0113 -0.00754 -0.00465 0.00369 -0.01431 ...
#> ..- attr(*, "names")= chr [1:2225] "cg0000001" "cg0000002" "cg0000003" "cg0000004" ...
#> $ se : Named num [1:2225] 0.0197 0.0197 0.0197 0.0197 0.0197 ...
#> ..- attr(*, "names")= chr [1:2225] "cg0000001" "cg0000002" "cg0000003" "cg0000004" ...
#> $ z : Named num [1:2225] -0.572 -0.382 -0.236 0.187 -0.725 ...
#> ..- attr(*, "names")= chr [1:2225] "cg0000001" "cg0000002" "cg0000003" "cg0000004" ...
#> $ df : int 39
#> $ n_subjects: int 24
#> $ n_samples : int 66The 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.
# 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)
#> block null region
#> 2.967 0.629 15.450var_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.
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))
#> clusters clustered_probes singletons
#> 642 2220 5Cluster 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.
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))
#> regions max_abs_z
#> 26.00 30.67
head(reg[order(-abs(reg$z)),
c("chr", "start", "end", "n_probes", "effect_M", "z")], 5)
#> chr start end n_probes effect_M z
#> 303 chr21 4727790 4729167 4 0.30290881 30.669435
#> 304 chr21 4735476 4736527 3 0.30357855 26.654351
#> 456 chr21 6486497 6488011 5 0.09585164 10.875439
#> 443 chr21 6315208 6315667 3 0.10278703 9.027847
#> 433 chr21 6200642 6201611 4 0.08380458 8.472760min_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.
folds <- fold_assign(subject, n_folds = 5L, seed = 0)
lengths(folds) # subjects per fold, not samples
#> 1 2 3 4 5
#> 5 5 5 5 4within_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.
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().
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))
#> [1] TRUEThe 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.
gr <- as_granges(obs)
gr[order(-abs(GenomicRanges::mcols(gr)$z))][1:3]
#> GRanges object with 3 ranges and 8 metadata columns:
#> seqnames ranges strand | n_probes effect_M se
#> <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric>
#> 303 chr21 4727790-4729167 * | 4 0.3029088 0.00987657
#> 304 chr21 4735476-4736527 * | 3 0.3035786 0.01138946
#> 456 chr21 6486497-6488011 * | 5 0.0958516 0.00881359
#> z p_pointwise theta probe_start probe_end
#> <numeric> <numeric> <numeric> <integer> <integer>
#> 303 30.6694 1.45539e-206 0.3029088 1042 1045
#> 304 26.6544 1.59346e-156 0.3035786 1046 1048
#> 456 10.8754 1.50927e-27 0.0958516 1575 1579
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengths
S4Vectors::metadata(gr)$max_abs_z
#> [1] 30.66944SummarizedExperiment 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.
perm <- within_subject_permutation(exposure, subject)
identical(sort(exposure[subject == "S01"]), sort(perm[subject == "S01"]))
#> [1] TRUE
fit_perm <- fit_within(M, perm, subject, covars)
c(observed = round(max(abs(fit$z)), 2),
permuted = round(max(abs(fit_perm$z)), 2))
#> observed permuted
#> 15.91 3.60stability_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.
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)
#> chr start end n_probes z stability
#> 303 chr21 4727790 4729167 4 30.669435 1
#> 304 chr21 4735476 4736527 3 26.654351 1
#> 456 chr21 6486497 6488011 5 10.875439 1
#> 443 chr21 6315208 6315667 3 9.027847 1
#> 433 chr21 6200642 6201611 4 8.472760 1The 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.
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)
#> [1] 642 66fit_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))
#> state mu sigma pi
#> 1 hypo -0.0052 0.0068 0.000
#> 2 neutral 0.0000 0.0001 0.915
#> 3 hyper 0.0681 0.0466 0.085hs$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.
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
#> chr start end width n_clusters direction posterior
#> 1 chr21 6078741 6696728 617988 51 hyper 0.997824planted <- 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)
Cluster-level effects along the chromosome. Shaded span: the planted block. Bars: the blocks called.
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, which generates this package from its numerical core, and the wrappers are published as a Tool Shed-installable suite at bioinfbrad/ewas-crossarray-harmonise; both extend the EWASGalaxy tool suite (Murat et al., bioRxiv 2019, doi:10.1101/553784).
Session info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] crossarrayEWAS_0.99.0 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.7-6 limma_3.99.0
#> [3] jsonlite_2.0.0 compiler_4.6.1
#> [5] BiocManager_1.30.27 SummarizedExperiment_1.43.0
#> [7] Biobase_2.73.2 GenomicRanges_1.65.4
#> [9] jquerylib_0.1.4 IRanges_2.47.5
#> [11] Seqinfo_1.3.2 statmod_1.5.2
#> [13] yaml_2.3.12 fastmap_1.2.0
#> [15] lattice_0.23-1 XVector_0.53.0
#> [17] R6_2.6.1 S4Arrays_1.13.0
#> [19] generics_0.1.4 knitr_1.52
#> [21] BiocGenerics_0.59.12 DelayedArray_0.39.6
#> [23] MatrixGenerics_1.25.0 maketools_1.3.2
#> [25] bslib_0.12.0 rlang_1.3.0
#> [27] cachem_1.1.0 xfun_0.61
#> [29] sass_0.4.10 sys_3.4.3
#> [31] otel_0.2.0 SparseArray_1.13.2
#> [33] cli_3.6.6 digest_0.6.39
#> [35] grid_4.6.1 lifecycle_1.0.5
#> [37] S4Vectors_0.51.10 evaluate_1.0.5
#> [39] buildtools_1.0.0 abind_1.4-8
#> [41] stats4_4.6.1 rmarkdown_2.32
#> [43] matrixStats_1.5.0 tools_4.6.1
#> [45] htmltools_0.5.9