library(BatchVaria)
library(SummarizedExperiment)
library(BiocStyle)

1 Overview

High-throughput transcriptomic datasets frequently contain unwanted technical variation arising from batch effects. Numerous batch correction methods exist, but selecting an appropriate method is challenging because correction can remove biological signal or distort sample relationships. BatchVaria provides a framework for quantifying variance components, tracking variance changes after correction, and evaluating correction performance using multiple metrics, enabling evidence-based selection of batch correction strategies.

While sva, limma and BatchQC provide batch correction and diagnostic visualisations, BatchVaria focuses specifically on quantitative evaluation of correction performance across multiple assays and correction strategies within a unified framework. BatchVaria unifies variance decomposition, tracking, and comparison into a single coherent workflow, built on SummarizedExperiment and using variancePartition as one of its two variance engines.

If you use BatchVaria in published work, please cite Moir N, Sherwood K, Simpson TI (2026), BatchVaria: a variance-aware framework for evaluating batch correction in high-dimensional omics data, bioRxiv, doi:10.64898/2026.05.07.721996. citation("BatchVaria") prints the same entry.

2 Installation

Install the release version with BiocManager:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("BatchVaria")

For the development version:

BiocManager::install(version = "devel")
BiocManager::install("BatchVaria")

The development source is at BatchVaria.

3 Example Dataset

We begin with a synthetic dataset with partially confounded batch and biology. This creates a situation where effective batch correction can be challenging. Within this dataset, most control samples sit in one batch and most treated samples in the other, so effects are not cleanly separable.

## exampleBatchVaria() draws random data; seed for a reproducible vignette
set.seed(1)
bv <- exampleBatchVaria(confounding = 0.6)

class(bv)
## [1] "SummarizedExperiment"
## attr(,"package")
## [1] "SummarizedExperiment"
table(
    batch = bv$batch,
    group = bv$group
)
##      group
## batch Ctrl Trt
##    B1    8   2
##    B2    2   8

BatchVaria leverages SummarizedExperiment for storage. Expression matrices are assays, sample annotation is colData(), and everything the package records goes in metadata(). See ?"BatchVaria-requirements" for further detail.

Correction methods are registered rather than hard-coded, so the set available is extensible in the same way the variance engines are:

availableCorrectionMethods()
## [1] "combat" "limma"

This vignette uses combat throughout, so that the comparison below isolates the effect of the design rather than differences between methods.

ComBat is applied twice. The first call corrects for batch alone. The second additionally names group as a variable to preserve, so the batch effect is estimated adjusted for the named preservation variable.

bv <- runCorrection(
    bv,
    method = "combat", batch = "batch",
    newAssayName = "naive"
)

bv <- runCorrection(
    bv,
    method = "combat", batch = "batch", preserve = "group",
    newAssayName = "protected"
)

assayNames(bv)
## [1] "raw"        "raw_center" "raw_scale"  "raw_noise"  "naive"     
## [6] "protected"

raw_noise is a deliberately distorted assay supplied with the example data, kept here as a reference point for what a genuinely destructive transformation looks like. raw_center and raw_scale are no-op controls: centring and scaling change no variance fraction, so a diagnostic that reports a change on them is reporting an artefact of the diagnostic rather than of the data.

4 Variance Profiling

We explicitly quantify the contribution of known covariates.

assaysUsed <- c("raw", "naive", "protected", "raw_noise")

bv <- profileVariance(
    bv,
    formula = ~ batch + group,
    assays = assaysUsed,
    methods = c("anova", "variancePartition")
)
## Warning: The shared variance component is negative (-0.035), meaning the model
## terms explain more jointly than separately. This indicates a suppression effect
## in the design

This step fits the specified model to each assay and decomposes variance into batch (unwanted), group (biological) and residual (unexplained). Two engines are run: anova, which fits linear models per feature and partitions variance with Type II sums of squares, and variancePartition, which fits mixed models. Every result is stored in an internal variance ledger keyed on assay, method, formula and term.

Because the ledger now holds two decompositions per assay, summaries must say which one they mean. Asking for a table without choosing returns an error rather than a silent mixture, since the engines’ terms are not interchangeable:

varianceTable(bv, assays = assaysUsed)
## Error in `.resolveMethod()`:
## ! The variance ledger holds results for 2 methods (anova, variancePartition). Their terms are not comparable, so specify method = to choose one

4.1 What the object is carrying

provenance() reports what has been done, reading both ledgers:

p <- provenance(bv)
## BatchVaria provenance
## ---------------------
## 
## Dimensions:
##   features: 800 
##   samples : 20 
## 
## Assays (6):
## [1] "raw"        "raw_center" "raw_scale"  "raw_noise"  "naive"     
## [6] "protected" 
## 
## Correction history:
##  step method input    output batch preserve no_op                  timestamp
##     1 combat   raw     naive batch        - FALSE 2026-10-13 16:26:15.980271
##     2 combat   raw protected batch    group FALSE 2026-10-13 16:26:16.060592
## 
## Variance profiling history:
##  step     assay        formula           methods                  timestamp
##     1       raw ~batch + group             anova 2026-10-13 16:26:16.133166
##     2       raw ~batch + group variancePartition 2026-10-13 16:26:20.021419
##     3     naive ~batch + group             anova 2026-10-13 16:26:20.038455
##     4     naive ~batch + group variancePartition 2026-10-13 16:26:23.032021
##     5 protected ~batch + group             anova 2026-10-13 16:26:23.049358
##     6 protected ~batch + group variancePartition 2026-10-13 16:26:26.120103
##     7 raw_noise ~batch + group             anova 2026-10-13 16:26:26.138851
##     8 raw_noise ~batch + group variancePartition 2026-10-13 16:26:29.179335

Two things there are worth pointing at. The correction table’s preserve column is what distinguishes the two ComBat runs - without it the rows are identical but for the assay names, and the comparison this vignette is built on would not be reconstructable from the object. The no_op column flags a correction that returned an assay indistinguishable from its input; a method can succeed and change nothing, and recording it on the object means the condition outlives the console session that produced it.

It prints that report and returns the same content as a data.frame, so the history can be filtered or joined rather than only read:

p$corrections[, c("method", "output", "preserve", "no_op")]
##   method    output preserve no_op
## 1 combat     naive        - FALSE
## 2 combat protected    group FALSE

5 Variance Summary Tables

A core output of BatchVaria is a pair of summary tables.

vt <- varianceTable(bv, assays = assaysUsed, method = "anova")

6 Variance Composition Table

This table shows the percentage of variance explained by each component for each assay.

vt$percent
## # A tibble: 4 × 5
##   component   raw naive protected raw_noise
##   <chr>     <dbl> <dbl>     <dbl>     <dbl>
## 1 batch      11.0   5.6      1.07      7.86
## 2 group      11.3  15.2     13.7       7.89
## 3 shared     11.4  -3.5      6.96      4.94
## 4 residual   66.3  82.7     78.3      79.3

7 Variance Change Relative to Raw

This table shows how each correction changed variance relative to the raw data. Every entry is a change in a fraction, which constrains what it can be read to mean; the next section sets out what that constraint is.

varianceTable(
    bv,
    assays = assaysUsed,
    method = "anova",
    formatDelta = TRUE
)$delta
## # A tibble: 4 × 4
##   component naive  protected raw_noise
##   <chr>     <chr>  <chr>     <chr>    
## 1 batch     -5.41  -9.94     -3.15    
## 2 group     +3.92  +2.41     -3.4     
## 3 shared    -14.91 -4.45     -6.48    
## 4 residual  +16.4  +11.98    +13.02

8 Total Variance: Reading Fractions Correctly

Every number above is a fraction. Within an assay, the fractions sum to one. This has a consequence that must be accounted for before any of these tables can be interpreted: removing variance associated with one term mechanically inflates the fraction attributed to every remaining term, even when the absolute variance of those terms has not increased at all.

assayVariance() supplies the missing denominator.

av <- assayVariance(bv, assays = assaysUsed)
av
##       assay total_variance n_features n_samples
## 1       raw      1281.1431        800        20
## 2     naive       828.6556        800        20
## 3 protected       927.4081        800        20
## 4 raw_noise      3077.1481        800        20

Multiplying a fraction by its assay’s total variance converts it to an absolute quantity, which is directly comparable across assays.

vr <- varianceResults(bv, method = "anova")

absVar <- merge(
    vr[, c("assay", "term", "variance_fraction")],
    av[, c("assay", "total_variance")],
    by = "assay"
)
absVar$absolute <- round(absVar$variance_fraction * absVar$total_variance, 1)

## as.data.frame(), because pivot_wider() returns a tibble and tibble
## prints three significant figure.
as.data.frame(tidyr::pivot_wider(
    absVar[, c("assay", "term", "absolute")],
    names_from = "term",
    values_from = "absolute"
))
##       assay batch group shared residual
## 1     naive  46.4 126.0  -29.0    685.3
## 2 protected   9.9 127.0   64.6    725.9
## 3       raw 141.0 144.6  146.2    849.4
## 4 raw_noise 241.9 242.7  151.9   2440.6

Read down the group column: 144.6 in raw, 126.0 after the naive correction and 127.0 after the protected one. In absolute terms group variance barely moves. Its fraction meanwhile climbs from 11.3 to 15.2 per cent, because the denominator fell from 1281 to 829. The fraction rose because something else was removed, not because biological signal grew.

9 Two Engines, One Answer - Until the Design Stops Cooperating

BatchVaria harmonises independent variance-profiling engines onto a common schema: the same columns, and the same names for the two terms that mean the same thing in any engine - residual and shared. Both runs above reach their decompositions by different routes: anova fits a linear model per feature and partitions with Type II sums of squares, while variancePartition fits mixed models via lme4.

cmp <- varianceResults(bv, assayName = "raw")[
    , c("method", "term", "variance_fraction")
]

as.data.frame(
    tidyr::pivot_wider(
        cmp,
        names_from = "method",
        values_from = "variance_fraction"
    )
)
##       term     anova variancePartition
## 1    batch 0.1100799         0.1698180
## 2    group 0.1128310         0.1751827
## 3   shared 0.1141206                NA
## 4 residual 0.6629684         0.6549993

Here they disagree, and the disagreement is the point. On an orthogonal design they do not:

orth <- profileVariance(
    exampleBatchVaria(confounding = 0),
    ~ batch + group,
    assays = "raw",
    methods = c("anova", "variancePartition")
)

## Printed at full precision deliberately: at three significant figures
## both columns read 0.163 and the agreement below cannot be checked.
as.data.frame(
    tidyr::pivot_wider(
        varianceResults(orth, assayName = "raw")[
            , c("method", "term", "variance_fraction")
        ],
        names_from = "method",
        values_from = "variance_fraction"
    )
)
##       term     anova variancePartition
## 1    batch 0.1695154         0.1695154
## 2    group 0.1598141         0.1598141
## 3   shared 0.0000000                NA
## 4 residual 0.6706704         0.6706704

Every term the two engines share agrees to floating-point precision. shared is exactly zero, because an orthogonal design leaves nothing unattributable, and variancePartition reports NA for it because a mixed model has no such parameter to report.

The divergence above is therefore a property of the design, not of either implementation. The two engines estimate different quantities: variancePartition reports how much variance is associated with a term, anova with Type II sums of squares reports how much is uniquely attributable to it. When a design is orthogonal these coincide, which is why the two tables differ as they do. When batch and group are partially aliased they cannot coincide, because roughly 11 per cent of the variance here is consistent with either term and identifiable to neither.

anova reports that share separately, as shared. A per-gene mixed model has no corresponding parameter - and could not, since the effects are not separately identifiable in this design - so the ambiguous variance is reflected in both terms, each rising from about 11 per cent to about 17. Neither number is wrong for the quantity it estimates; they simply stop being interchangeable once the design cannot separate the effects.

The gap between them is therefore a measurement of the design rather than a disagreement between tools: it says how much of the variance the experiment permits you to attribute at all. Running both, and seeing them coincide on a clean design and separate on a compromised one, is more informative than trusting either alone.

10 Interpretation

The naive correction looks like a success. Read the fraction table on its own: batch variance halves, from 11.0 per cent to 5.6, and group variance rises, from 11.3 per cent to 15.2. Unwanted variation down, biological signal up. By the standard before-and-after comparison - decompose variance, correct, decompose again - this is a correction working exactly as intended.

It is not. Three things the fraction table cannot show say otherwise.

Total variance fell by 35 per cent, from 1281 to 829. The group fraction rose because the denominator collapsed, not because biology was preserved. In absolute terms group variance fell from 145 to 126. This is the compositional artefact: within an assay the fractions must sum to one, so removing variance from any term mechanically inflates the share of every other. A rising fraction is not evidence of preserved signal, and reading it as such is the easiest way to certify a bad correction as a good one.

The batch effect was not removed. 5.6 per cent remains, half of what was there.

shared went negative, from +11.4 per cent to -3.5. The correction did not merely fail to resolve the confounding; it distorted the design’s attributability, leaving the model terms explaining more jointly than they do separately.

Protecting the covariate addresses the second of these. Telling ComBat that group matters lets it estimate the batch effect adjusted for biology, and batch variance drops to 1.1 per cent - essentially removed - while total variance falls less far, to 927. It does not rescue biological signal on this data: absolute group variance is 126 for the naive correction and 127 for the protected one. The case for protecting the covariate here is that the correction works, not that the biology survives better.

That balance is specific to this level of confounding. This vignette uses a single value, confounding = 0.6; the trade-off moves with it, because as batch and group approach aliasing, protecting the covariate increasingly means protecting the batch effect alongside it. Raising or lowering the parameter will change which correction looks better on absolute group variance, so treat that comparison as a property of the design rather than a general ranking of the two calls.

The sign of shared is the more dependable signal. Wherever there is any confounding to resolve it is consistently negative after naive correction and positive after protected correction, so a correction that drives it below zero has distorted the design’s attributability regardless of what the other fractions report. Being a sign rather than a magnitude, that reading needs no threshold and no reference distribution - but note it is demonstrated here on simulated data with one correction method and one design shape.

For contrast, raw_noise shows what a destructive transformation looks like under the same instruments. Its batch and group fractions both fall, which reads as damage. But its total variance more than doubles, to 3077, so in absolute terms every component grew: nothing was removed, the batch effect was diluted. basisRetention() is unambiguous about the difference.

basisRetention(bv, assays = assaysUsed)
##       assay reference retention  off_basis n_features
## 1       raw       raw 1.0000000 0.00000000        800
## 2     naive       raw 0.9891828 0.01081725        800
## 3 protected       raw 0.9951290 0.00487105        800
## 4 raw_noise       raw 0.3699325 0.63006753        800

It retains 0.37 of the original variance structure against 0.99 for both ComBat runs, because added noise occupies directions the raw data did not. Note what this measure does and does not separate: it distinguishes a correction from a distortion, not a good correction from a bad one. Both ComBat runs stay on the original axes because ComBat adjusts each feature’s location and scale within a batch, which rescales the existing structure rather than introducing new directions.

The general rule: a variance fraction is a statement about composition, not about quantity. Use fractions to compare terms within an assay, and totals or absolute variances to compare the same term across assays.

This is also the answer to what BatchVaria offers over calling the underlying tools twice. Decomposing variance before and after correction is exactly what certified the naive correction as successful. Recording each decomposition alongside the total variance it is relative to, the share the design cannot attribute, and the structure the correction preserved is what catches it.

11 Summary Figure

evaluateCorrections() runs the diagnostics above in one call and returns them in a form the plotting layer reads directly.

ev <- evaluateCorrections(
    bv,
    baseline = "raw",
    assays = assaysUsed,
    diagnostics = "variance"
)

plotVarianceDelta(ev)

One group of bars per model term, one bar per assay, and height is the signed change in that term’s fraction relative to raw - as a fraction here, where the delta table above reported the same quantities as percentages. Terms are drawn in alphabetical order: batch, group, residual, shared.

Both corrections cut batch, the protected run roughly twice as hard, and both raise residual. The bar that separates them is shared. Both push it down, and by different amounts - 0.149 for the naive run against 0.044 for the protected one - but the amount is not the point. What distinguishes them is where they land: the naive run takes shared from +0.114 to -0.035, while the protected run takes it from +0.114 to +0.070. As above, a negative shared means the terms explain more jointly than they do separately, so it is the side of zero that carries the finding, not the size of the drop. raw_noise cuts batch too, by dilution rather than correction, which is why no single bar here is sufficient on its own.

12 Diagnostic Plots

The argument above is made with tables, deliberately: every claim in it is a number that can be checked. The plotting layer draws those same quantities, and is collected here rather than interleaved so that the two are not confused. Each plot below answers one question.

evaluateCorrections() computes all three diagnostic families in one call and returns them in the form the plotting functions read:

ev <- evaluateCorrections(
    bv,
    baseline = "raw",
    assays = assaysUsed,
    varianceMethod = "anova"
)

names(ev)
## [1] "variance_delta"     "pca_comparison"     "correlation_change"

12.1 Where the variance sits

plotVarianceComposition() draws each assay’s decomposition as a composition: one stacked bar per assay, summing to 100 per cent. It answers “what is this assay made of”, and needs no baseline - an assay profiled on its own can be drawn.

plotVarianceComposition(bv, assays = assaysUsed, method = "anova")
## Warning: Negative variance percentages for term(s) shared. The stacked bars
## extend below zero and above 100 per cent; read the values from varianceTable()
## rather than the chart

Read it against the totals from assayVariance() above, not on its own: the bars are shares, so raw_noise shows a smaller batch share than raw while its absolute batch variance has risen, from 141 to 242. The share fell only because the total more than doubled underneath it.

plotVarianceRadar() shows the same deltas as the summary figure on polar axes, one spoke per term. It is the quicker read when several corrections are being compared at once and the question is the shape of the change rather than its size - two corrections that move the same terms in the same proportions trace the same outline.

plotVarianceRadar(bv, assays = assaysUsed, method = "anova")

12.2 Where the samples sit

Every assay is projected onto a single basis fitted on the reference assay, rather than each being fitted separately. This is what makes the panels comparable: a per-assay fit gives PC1 a different meaning in each, so samples appear to move when only the axes have.

The coordinates are available directly, not only through a plot:

coords <- basisProjection(
    bv,
    assays = c("raw", "naive"),
    reference = "raw"
)

head(coords)
##   assay sample reference       PC1        PC2
## 1   raw     S1       raw -15.28668  0.6697611
## 2   raw     S2       raw -16.21937 -1.0138277
## 3   raw     S3       raw -16.29337  2.1242652
## 4   raw     S4       raw -15.99028 -0.0769103
## 5   raw     S5       raw -15.06171 -0.9763562
## 6   raw     S6       raw -16.15882 -1.0309512

plotBasisProjection() draws exactly what that returns, one panel per assay:

for (p in plotBasisProjection(
    bv,
    assays = c("raw", "naive", "protected"),
    reference = "raw"
)) {
    print(p)
}

The naive panel repays a careful look, because it does not show what a reader expects of a corrected assay: the batches still appear separated. They have in fact been brought together - the distance between batch means along PC1 falls from 1.73 to 0.54 standard deviations. What remains on the panel is group, whose separation rises over the same step, from 1.75 to 1.81. At confounding = 0.6 the batch and group labels agree for sixteen of the twenty samples, so a panel coloured by batch draws surviving group structure in batch colours, and the four samples whose labels disagree are the only ones that give it away.

Colouring by the other variable shows what is actually there:

plotBasisProjection(
    bv,
    assays = "naive",
    reference = "raw",
    colourBy = "group"
)[[1]]

This is a limit of the panel rather than of the correction, and it is not peculiar to this function: wherever batch and biology are confounded, a projection coloured by batch cannot distinguish a batch effect that survived from biological signal that was preserved. Here it ranks the two corrections backwards: set the two corrected panels side by side and protected looks more batch-separated than naive - 1.31 standard deviations against 0.54 - while carrying a fifth of the batch variance, 1.1 per cent against 5.6. The variance table settles what the panel cannot.

There is a second thing these panels cannot show. raw_noise, drawn on this basis, is nearly indistinguishable from raw: the noise it added lies almost entirely off the reference basis, so it never enters these coordinates, despite the assay carrying more than twice the variance. Projecting onto a basis shows what happened within the retained structure and is silent about structure that left it, which is why basisRetention() is reported alongside it - scoring raw_noise at 0.37 against 0.99 for both ComBat runs.

plotPCATrajectory() makes the same comparison as a single panel, drawing an arrow from each sample’s position before correction to its position after. Use it when the question is which samples moved rather than whether the groups separate.

plotPCATrajectory(bv, assayBefore = "raw", assayAfter = "naive")

It inherits the caveat above: an arrow shows that a sample moved, not which of the two confounded variables moved it.

12.3 Whether the sample relationships survived

A correction can leave every variance fraction looking reasonable and still rearrange which samples resemble which. Two diagnostics check that directly.

plotSampleDistance() draws the sample-sample distance matrix, so blocks of mutually similar samples are visible as blocks:

plotSampleDistance(bv, assayName = "raw")

plotCorrelationChange() reduces that to one number per assay: the mean absolute change in sample-sample correlation relative to the baseline. Large values mean the correction rewrote the relationships between samples, which is a different failure from removing too much variance and is not visible in any of the variance tables.

plotCorrelationChange(ev)

Note which assay scores best here: raw_noise, at 0.100 against 0.187 for the naive correction. The assay this vignette has used throughout as its example of a destructive transformation is the one that disturbs sample relationships least. That is not a defect in the diagnostic, it is what the diagnostic measures - adding independent noise to every feature shrinks all sample correlations towards zero fairly evenly, which is a small change in the correlation structure, whereas shifting feature means per batch rearranges it more. The two diagnostics agree that the protected run is gentler than the naive one; where they part company is raw_noise, which this measure ranks first and basisRetention() ranks last. That single disagreement is the clearest illustration in this vignette of why no one diagnostic settles the question.

12.4 Whether the batches actually mixed

plotBatchEntropy() measures local batch mixing: for each sample it looks at its k nearest neighbours and computes the entropy of their batch labels. High entropy means a sample’s neighbourhood is batch-mixed; low entropy means it still sits among its own batch. This is the neighbourhood-level question the global variance fractions cannot answer.

It draws one assay at a time, so the comparison is made by drawing two:

plotBatchEntropy(bv, assayName = "raw", batchVar = "batch", k = 10)
plotBatchEntropy(bv, assayName = "naive", batchVar = "batch", k = 10)

The distribution moves right after correction, which is the batch effect leaving the local structure.

Under partial confounding this should be read with care, and this design shows why. With two batches the entropy is bounded above by log(2), about 0.693; after the naive correction the mean neighbourhood entropy is 0.65, close to that ceiling. But batch and group overlap here, so neighbourhoods that are perfectly batch-mixed are neighbourhoods in which the biological groups have been mixed too. Entropy approaching its maximum is evidence of overcorrection in this design rather than of a job well done - which is the same verdict the shared term reached, arrived at independently.

13 Why This Matters

Batch correction is frequently applied as a routine preprocessing step, yet it is often performed without systematic or quantitative validation. This can lead to suboptimal outcomes, where overcorrection removes meaningful biological signal or undercorrection leaves residual technical artefacts that confound downstream analyses. BatchVaria addresses this gap by providing a transparent and reproducible framework for evaluating correction strategies. By enabling direct comparison of variance components across multiple corrected assays, the package supports quantitative validation of correction outcomes and facilitates more robust, interpretable analytical pipelines.

14 BatchVaria Workflow Summary

The typical BatchVaria workflow is:

  1. Load expression data into a SummarizedExperiment
  2. Apply batch correction methods with runCorrection(), naming any variables to preserve
  3. Profile variance using profileVariance()
  4. Summarise variance using varianceTable(method = …), naming the engine
  5. Convert fractions to absolute variance using assayVariance()
  6. Compare corrections using delta variance and basisRetention()
  7. Visualise variance composition

15 Summary

BatchVaria provides a coherent and extensible framework for evaluating batch correction in high-dimensional expression data. By combining variance decomposition, structured tracking of results across multiple assays, and intuitive summary outputs, the package enables users to quantitatively assess how different correction strategies impact both unwanted technical variation and biological signal. The integration of variance composition tables, relative change metrics, and visual summaries supports transparent comparison of correction methods within a single workflow. Together, these features facilitate reproducible and interpretable evaluation of batch correction, addressing a key gap in many transcriptomic analysis pipelines.

16 Session Information

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] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] BatchVaria_0.99.2           SummarizedExperiment_1.43.0
##  [3] Biobase_2.73.2              GenomicRanges_1.65.4       
##  [5] Seqinfo_1.3.2               IRanges_2.47.5             
##  [7] S4Vectors_0.51.10           BiocGenerics_0.59.12       
##  [9] generics_0.1.4              MatrixGenerics_1.25.0      
## [11] matrixStats_1.5.0           BiocStyle_2.41.0           
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3       jsonlite_2.0.0           magrittr_2.0.5          
##   [4] magick_2.9.1             farver_2.1.2             nloptr_2.2.1            
##   [7] rmarkdown_2.32           vctrs_0.7.3              memoise_2.0.1           
##  [10] minqa_1.2.8              tinytex_0.61             htmltools_0.5.9         
##  [13] S4Arrays_1.13.2          BiocBaseUtils_1.15.1     broom_1.0.13            
##  [16] SparseArray_1.13.4       variancePartition_1.43.1 sass_0.4.10             
##  [19] KernSmooth_2.23-27       bslib_0.12.0             pbkrtest_0.5.5          
##  [22] plyr_1.8.9               cachem_1.1.0             lifecycle_1.0.5         
##  [25] iterators_1.0.14         pkgconfig_2.0.3          Matrix_1.7-6            
##  [28] R6_2.6.1                 fastmap_1.2.0            rbibutils_2.4.1         
##  [31] digest_0.6.39            numDeriv_2016.8-1.1      AnnotationDbi_1.75.2    
##  [34] RSQLite_3.53.3           labeling_0.4.3           httr_1.4.9              
##  [37] abind_1.4-8              mgcv_1.9-4               compiler_4.6.1          
##  [40] bit64_4.8.6              aod_1.3.3                withr_3.0.3             
##  [43] S7_0.2.2                 backports_1.5.1          BiocParallel_1.47.0     
##  [46] DBI_1.3.0                gplots_3.3.0             MASS_7.3-66             
##  [49] DelayedArray_0.39.8      corpcor_1.6.10           gtools_3.9.5            
##  [52] caTools_1.18.4           tools_4.6.1              otel_0.2.0              
##  [55] remaCor_0.0.20           glue_1.8.1               nlme_3.1-171            
##  [58] grid_4.6.1               reshape2_1.4.5           sva_3.61.0              
##  [61] gtable_0.3.6             tidyr_1.3.2              utf8_1.2.6              
##  [64] XVector_0.53.0           pillar_1.11.1            stringr_1.6.0           
##  [67] limma_3.99.2             genefilter_1.95.0        splines_4.6.1           
##  [70] dplyr_1.2.1              lattice_0.23-1           survival_3.8-12         
##  [73] FNN_1.1.4.1              bit_4.6.0                annotate_1.91.0         
##  [76] tidyselect_1.2.1         locfit_1.5-9.12          Biostrings_2.81.9       
##  [79] knitr_1.52               reformulas_0.4.4         bookdown_0.48           
##  [82] edgeR_4.99.7             RhpcBLASctl_0.23-42      xfun_0.61               
##  [85] statmod_1.5.2            stringi_1.8.9            yaml_2.3.12             
##  [88] boot_1.3-32              evaluate_1.0.5           codetools_0.2-20        
##  [91] tibble_3.3.1             BiocManager_1.30.27      cli_3.6.6               
##  [94] xtable_1.8-8             Rdpack_2.6.6             jquerylib_0.1.4         
##  [97] dichromat_2.0-1          Rcpp_1.1.2               EnvStats_3.1.0          
## [100] png_0.1-9                XML_3.99-0.25            parallel_4.6.1          
## [103] ggplot2_4.0.3            blob_1.3.0               bitops_1.1-0            
## [106] lme4_2.0-6               viridisLite_0.4.3        mvtnorm_1.4-2           
## [109] lmerTest_3.2-1           scales_1.4.0             purrr_1.2.2             
## [112] crayon_1.5.3             fANCOVA_0.6-1            rlang_1.3.0             
## [115] KEGGREST_1.53.6