Design-aware region and block detection for cross-array longitudinal EWAS

Katarzyna Kamieniecka, Samantha McLean, Krzysztof Poterlowicz

Package

crossarrayEWAS 0.99.1

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

2 Installation

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

library(crossarrayEWAS)

3 The two functions most users need

ewas_regions() and ewas_blocks() take a SummarizedExperiment — including a GenomicRatioSet from 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.

4 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       38

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.

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              640

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.

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   66

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.

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
#> class: RangedSummarizedExperiment 
#> dim: 2225 66 
#> metadata(0):
#> assays(1): M
#> rownames(2225): cg0000001 cg0000002 ... cg0002224 cg0002225
#> rowData names(0):
#> colnames(66): S01_v1 S01_v2 ... S24_v2 S24_v3
#> colData names(5): exposure subject cell_fraction smoking_score
#>   array_type

5 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)))]
#> GRanges object with 5 ranges and 11 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
#>   443    chr21 6315208-6315667      * |         3 0.1027870 0.01138555
#>   433    chr21 6200642-6201611      * |         4 0.0838046 0.00989106
#>               z  p_pointwise     theta probe_start probe_end    p_fwer
#>       <numeric>    <numeric> <numeric>   <integer> <integer> <numeric>
#>   303  30.66944 1.45539e-206 0.3029088        1042      1045 0.0476190
#>   304  26.65435 1.59346e-156 0.3035786        1046      1048 0.0476190
#>   456  10.87544  1.50927e-27 0.0958516        1575      1579 0.0952381
#>   443   9.02785  1.75082e-19 0.1027870        1535      1537 0.1428571
#>   433   8.47276  2.39646e-17 0.0838046        1506      1509 0.1428571
#>       p_fwer_mcse stability
#>         <numeric> <numeric>
#>   303   0.0476190         1
#>   304   0.0476190         1
#>   456   0.0656383         1
#>   443   0.0782461         1
#>   433   0.0782461         1
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

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.

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)
#> max_abs_z    lambda  clusters        df 
#>     30.67      0.50    642.00     39.00

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

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))
#>                        regions overlapping_the_planted_region 
#>                             26                              2

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.

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

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
#> GRanges object with 1 range and 3 metadata columns:
#>       seqnames          ranges strand | n_clusters   direction posterior
#>          <Rle>       <IRanges>  <Rle> |  <integer> <character> <numeric>
#>   [1]    chr21 6078741-6696728      * |         51       hyper  0.997761
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

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.

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))
#>     state      mu  sigma    pi
#> 1    hypo -0.0052 0.0069 0.000
#> 2 neutral  0.0000 0.0001 0.914
#> 3   hyper  0.0681 0.0466 0.086
clusters <- S4Vectors::metadata(blocks)$clusters
length(clusters)
#> [1] 641

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)

Figure 1: 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.

7 From 450K and EPIC data

ewas_regions() needs one object spanning both array generations. Building it is minfi’s job, and it is three calls: combine the two RGChannelSets 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.

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.

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

path <- system.file("extdata", "GSE237561_chr21_subset.rds",
                    package = "crossarrayEWAS")
grset <- readRDS(path)
grset
#> Loading required namespace: minfi
#> Setting options('download.file.method.GEOquery'='auto')
#> Setting options('GEOquery.inmemory.gpl'=FALSE)
#> class: GenomicRatioSet 
#> dim: 2000 126 
#> metadata(0):
#> assays(1): Beta
#> rownames(2000): cg17035109 cg12459059 ... cg14296488 cg15214567
#> rowData names(0):
#> colnames(126): 3998567005_R01C01 3998567005_R02C01 ...
#>   204026590037_R02C01 204026590037_R03C01
#> colData names(4): subject visit exposure array_type
#> Annotation
#>   array: IlluminaHumanMethylation450k
#>   annotation: ilmn12.hg19
#> Preprocessing
#>   Method: GEO GSE237561 normalised beta matrix
#>   minfi version: NA
#>   Manifest version: NA
table(SummarizedExperiment::colData(grset)$array_type)
#> 
#> 450K EPIC 
#>   92   34

real <- ewas_regions(grset, exposure = "exposure", subject = "subject",
                     lambda = 0.5, min_probes = 3L, min_effect = 0.05)
#> assay "Beta" holds beta values; converting to M-values with log2(beta / (1 - beta))
length(real)
#> [1] 1
if (length(real))
  real[order(-abs(real$z))][seq_len(min(3L, length(real)))]
#> GRanges object with 1 range and 8 metadata columns:
#>        seqnames            ranges strand |  n_probes   effect_M        se
#>           <Rle>         <IRanges>  <Rle> | <integer>  <numeric> <numeric>
#>   1398    chr21 42218932-42219139      * |         3 -0.0510731 0.0211092
#>                z p_pointwise       theta probe_start probe_end
#>        <numeric>   <numeric>   <numeric>   <integer> <integer>
#>   1398  -2.41947   0.0155432 -0.00633606        1493      1495
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

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.

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

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))
#>   clusters  clustered singletons 
#>        642       2220          5

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

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))
#> [1] TRUE

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.

10 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.23

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 bumphunter uses, and exists to quantify its mis-calibration in a longitudinal design rather than to be used.

GenomicRanges::mcols(regions)[order(-abs(regions$z)),
                              c("z", "p_fwer", "p_fwer_mcse")][1:3, ]
#> DataFrame with 3 rows and 3 columns
#>             z    p_fwer p_fwer_mcse
#>     <numeric> <numeric>   <numeric>
#> 303   30.6694 0.0476190   0.0476190
#> 304   26.6544 0.0476190   0.0476190
#> 456   10.8754 0.0952381   0.0656383
S4Vectors::metadata(regions)$perm_scheme
#> [1] "within"

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.

GenomicRanges::mcols(regions)[order(-regions$stability, -abs(regions$z)),
                              c("n_probes", "z", "stability")][1:3, ]
#> DataFrame with 3 rows and 3 columns
#>      n_probes         z stability
#>     <integer> <numeric> <numeric>
#> 303         4   30.6694         1
#> 304         3   26.6544         1
#> 456         5   10.8754         1

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.

11 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. bumphunter, DMRcate, ENmix (combp(), ipdmr()) and 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 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.

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

13 Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [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: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] crossarrayEWAS_0.99.1 BiocStyle_2.41.0     
#> 
#> loaded via a namespace (and not attached):
#>   [1] bitops_1.1-0                beanplot_1.3.1             
#>   [3] DBI_1.3.0                   rlang_1.3.0                
#>   [5] magrittr_2.0.5              scrime_1.3.7               
#>   [7] otel_0.3.0                  matrixStats_1.5.0          
#>   [9] compiler_4.6.1              RSQLite_3.53.3             
#>  [11] GenomicFeatures_1.65.0      DelayedMatrixStats_1.35.0  
#>  [13] png_0.1-9                   vctrs_0.7.3                
#>  [15] quadprog_1.5-8              minfi_1.59.1               
#>  [17] pkgconfig_2.0.3             crayon_1.5.3               
#>  [19] fastmap_1.2.0               magick_2.9.1               
#>  [21] XVector_0.53.0              Rsamtools_2.29.0           
#>  [23] rmarkdown_2.32              tzdb_0.5.0                 
#>  [25] preprocessCore_1.75.3       purrr_1.2.2                
#>  [27] tinytex_0.61                bit_4.6.0                  
#>  [29] xfun_0.62                   cachem_1.1.0               
#>  [31] cigarillo_1.3.1             jsonlite_2.0.0             
#>  [33] blob_1.3.0                  rhdf5filters_1.25.5        
#>  [35] DelayedArray_0.39.8         reshape_0.8.10             
#>  [37] Rhdf5lib_2.1.0              BiocParallel_1.47.0        
#>  [39] parallel_4.6.1              R6_2.6.1                   
#>  [41] bslib_0.12.0                RColorBrewer_1.1-3         
#>  [43] limma_3.99.4                rtracklayer_1.73.0         
#>  [45] genefilter_1.95.0           GenomicRanges_1.65.4       
#>  [47] jquerylib_0.1.4             Rcpp_1.1.2                 
#>  [49] Seqinfo_1.3.2               bookdown_0.48              
#>  [51] SummarizedExperiment_1.43.0 iterators_1.0.14           
#>  [53] knitr_1.52                  readr_2.2.0                
#>  [55] IRanges_2.47.5              BiocBaseUtils_1.15.1       
#>  [57] rentrez_1.2.4               tidyselect_1.2.1           
#>  [59] illuminaio_0.55.0           Matrix_1.7-6               
#>  [61] splines_4.6.1               abind_1.4-8                
#>  [63] yaml_2.3.12                 siggenes_1.87.0            
#>  [65] codetools_0.2-20            curl_8.1.0                 
#>  [67] doRNG_1.8.6.3               tibble_3.3.1               
#>  [69] lattice_0.23-1              plyr_1.8.9                 
#>  [71] Biobase_2.73.2              KEGGREST_1.53.6            
#>  [73] askpass_1.2.1               evaluate_1.0.5             
#>  [75] survival_3.8-12             xml2_1.6.0                 
#>  [77] mclust_6.1.3                Biostrings_2.81.9          
#>  [79] pillar_1.11.1               BiocManager_1.30.27        
#>  [81] MatrixGenerics_1.25.0       rngtools_1.5.2             
#>  [83] foreach_1.5.2               stats4_4.6.1               
#>  [85] generics_0.1.4              RCurl_1.98-1.20            
#>  [87] hms_1.1.4                   S4Vectors_0.51.10          
#>  [89] sparseMatrixStats_1.25.0    bumphunter_1.55.1          
#>  [91] xtable_1.8-8                glue_1.8.1                 
#>  [93] tools_4.6.1                 BiocIO_1.23.3              
#>  [95] data.table_1.18.6.1         GenomicAlignments_1.49.2   
#>  [97] annotate_1.91.0             locfit_1.5-9.12            
#>  [99] GEOquery_2.81.36            XML_3.99-0.25              
#> [101] rhdf5_2.57.20               grid_4.6.1                 
#> [103] tidyr_1.3.2                 AnnotationDbi_1.75.2       
#> [105] base64_2.0.2                nlme_3.1-172               
#> [107] nor1mix_1.3-3               HDF5Array_1.41.3           
#> [109] restfulr_0.0.17             cli_3.6.6                  
#> [111] S4Arrays_1.13.2             dplyr_1.2.1                
#> [113] sass_0.4.10                 digest_0.6.39              
#> [115] BiocGenerics_0.59.12        SparseArray_1.13.4         
#> [117] rjson_0.2.23                memoise_2.0.1              
#> [119] htmltools_0.5.9             multtest_2.69.0            
#> [121] lifecycle_1.0.5             h5mread_1.5.3              
#> [123] httr_1.4.9                  statmod_1.5.2              
#> [125] openssl_2.5.0               bit64_4.8.6                
#> [127] MASS_7.3-66