This vignette sets out the protocol for evaluating sc-ExWAS inference by simulation: which methods can be compared, how data must be generated, which metrics are reported and what a frozen benchmark release has to contain. It runs two small examples, a complete-null calibration check and a power approximation, so that the mechanics can be seen; it reports no benchmark result, and the package makes no quantitative calibration, power, speed or superiority claim until a frozen release does.
The primary sc-ExWAS estimand is the exposure association with relative expression within a cell type. Differential-abundance methods target relative composition and are not expression-model competitors. A future molecular-contribution decomposition would target the product of composition and within-type expression and requires a distinct simulation truth.
Comparators are evaluated only when they share:
A cell-level test that treats cells as independent can be included as a negative control, but not as a plausible population-inference method.
The principal calibration protocol is a fixed-design negative-binomial parametric bootstrap. The observed exposure and nuisance design, donor and sample support, offsets, expression-support family and feature-by-cell-type keys remain fixed; counts are generated under a reduced model with a zero exposure coefficient, after which the same edgeR pipeline is rerun. The result is a conditional plug-in diagnostic for that exact generator and family, not proof of strong FDR control under arbitrary mixtures of true and false null hypotheses.
Permuting the exposure labels within a study is a stress audit, not a calibration. An observational exposure was not randomised and its conditional exchangeability cannot be assumed, so no amount of auditing of the accepted permutations makes the resulting rate a guarantee.
A calibration run is worth keeping only if it can be repeated: each replicate should retain its seed, its fit diagnostics, the call summary over the whole family, the package and input versions it ran on, and its run time. A failed fit stays a failure and cannot be recoded as no call.
The simulation below has the shape of a null replicate: donor-level biological variation in the counts, an exposure drawn independently of them, so that every null hypothesis is true, and the default pipeline run on the result. Five replicates of 20 donors, 60 genes and two cell types run in a few seconds; a formal run uses far more of each and keeps every replicate’s seed and receipt.
library(exposomeSC)
simulate_null <- function(seed, n_donors = 20L, n_genes = 60L,
cells_per_donor = 30L,
celltypes = c("Monocyte", "NK")) {
set.seed(seed)
donors <- sprintf("D%02d", seq_len(n_donors))
cell_donor <- rep(rep(donors, each = cells_per_donor),
times = length(celltypes))
cell_type <- rep(celltypes, each = n_donors * cells_per_donor)
donor_effect <- matrix(stats::rnorm(n_genes * n_donors, sd = 0.3),
n_genes, dimnames = list(NULL, donors))
mu <- 8 * exp(donor_effect[, cell_donor])
counts <- matrix(stats::rnbinom(length(mu), mu = mu, size = 10),
nrow = n_genes,
dimnames = list(paste0("G", seq_len(n_genes)),
paste0("c", seq_along(cell_donor))))
sce <- SingleCellExperiment::SingleCellExperiment(
assays = list(counts = counts),
colData = S4Vectors::DataFrame(cell_id = colnames(counts),
donor_id = cell_donor, cell_type = cell_type))
## drawn independently of the counts: every null is true
design <- cbind(exposure = stats::rnorm(n_donors),
age = round(stats::runif(n_donors, 30, 70)))
rownames(design) <- donors
build_scee(sce, design, sample_col = "donor_id")
}
null_replicate <- function(seed) {
fit <- run_sc_exwas(simulate_null(seed), exposure = "exposure",
covariates = "age", celltype_col = "cell_type",
sample_col = "donor_id", min_cells = 5L)
c(tests = nrow(fit),
calls = sum(fit$padj_global < 0.05, na.rm = TRUE),
raw_p_below_0.05 = mean(fit$pvalue < 0.05, na.rm = TRUE))
}
replicates <- t(vapply(1001:1005, null_replicate, numeric(3)))
replicates## tests calls raw_p_below_0.05
## [1,] 120 0 0.09166667
## [2,] 120 0 0.04166667
## [3,] 120 0 0.03333333
## [4,] 120 0 0.03333333
## [5,] 120 0 0.04166667Under the complete null the family-wise any-call indicator is
any(calls > 0) per replicate, and the fraction of raw
p-values below 0.05 should sit near 0.05:
c(any_call_rate = mean(replicates[, "calls"] > 0),
mean_raw_p_below_0.05 = mean(replicates[, "raw_p_below_0.05"]))## any_call_rate mean_raw_p_below_0.05
## 0.00000000 0.04833333Five replicates cannot bound a rejection rate; the protocol asks for independent replicates in the hundreds, with Wilson intervals over them.
Simulation must generate the hierarchy explicitly:
At minimum, the design should cross null and non-null expression effects with balanced and imbalanced exposure, heterogeneous study effects, repeated samples, rare cell types, unequal library sizes, influential donors and misspecified age trends. Parameter values must be justified from held-out or external summaries, not tuned to favour one implementation.
A pure composition scenario changes only the composition-generating process. It must not insert the composition score directly into within-type expression. A contextual-expression scenario may do so, but it is labelled and scored as a different mechanism.
For each independent simulation replicate, report:
False-discovery rate is the mean of replicate-level false-discovery proportions. Gene-by-replicate rows are correlated and must not be treated as independent observations when calculating uncertainty. Empirical-null rejection-rate intervals are calculated over independent simulation or bootstrap replicates. Failed fits remain failures; they are not recoded as no call.
estimate_power() gives the large-sample power of the
pseudobulk Wald test for one gene in one cell type at the Bonferroni
per-test level, from the donor count, the log2 effect per unit exposure,
the negative-binomial dispersion and the mean pseudobulk count. It
treats the dispersion as known and is optimistic for small donor
numbers; it is a planning aid, and the simulation above is the
check.
estimate_power(n_donors = c(20, 40, 80, 160), effect_size = 0.5,
n_genes = 5000, n_celltypes = 5,
dispersion = 0.1, base_mean = 100)## DataFrame with 4 rows and 7 columns
## n_donors power effect_size n_tests alpha_adjusted dispersion base_mean
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## 1 20 0.468028 0.5 25000 2e-06 0.1 100
## 2 40 0.968236 0.5 25000 2e-06 0.1 100
## 3 80 0.999998 0.5 25000 2e-06 0.1 100
## 4 160 1.000000 0.5 25000 2e-06 0.1 100Donors, not cells, are the independent units: more cells per donor
raise the mean count and shrink its term in the standard error, but they
do not add to n_donors.
run_sc_exwas() offers four backends, and a comparison is
only meaningful between backends that target the same estimand.
"edgeR" is the default and fits a robust quasi-likelihood
model to donor pseudobulk; it reports a signed square root of the F
statistic and no coefficient standard error. "DESeq2", with
"pseudobulk" as its alias, gives a Wald test on the same
donor pseudobulk. "dreamlet" applies variancePartition’s
voomWithDreamWeights() and dream() without a
random effect, which is a moderated t test; it does not call
the dreamlet package, so the label voom-dream describes it better.
Two exported functions leave donor pseudobulk.
run_sc_exwas_glmm() fits a cell-level mixed model with a
donor random intercept, either a linear model with Satterthwaite or
Kenward-Roger degrees of freedom, or a negative-binomial
lme4::glmer.nb fit tested by likelihood ratio.
run_sc_exwas_pb_offset() fits a donor-level count model
with an offset. Both target a conditional cell-level mean, not the
donor-aggregated contrast, so their coefficients are not interchangeable
with the pseudobulk ones.
A cell-level test that ignores the donor is a negative control, not a comparator. Method labels should describe the model that was fitted rather than the wrapper package it came from, and composition methods are assessed only against composition truth.
For a targeted model evaluation, test_features fixes the
reported family while the full feature matrix remains available for
library sizes, normalisation and dispersion estimation. A missing or
unsupported requested feature is reported as not estimable; the family
must not shrink silently.
Filtering and feature-family definitions are fixed without exposure labels. The primary adjusted values use the complete host-expression-by-cell-type family. Separate receptor and residual-annotation audits do not alter that denominator. Same-data screening followed by unadjusted decomposition or pathway inference is not a valid confirmatory procedure.
A frozen benchmark release must include:
Figures and tables read these objects directly. No metric is typed into display code. Until this record is complete, the package makes no quantitative calibration, power, speed or superiority claim.
The removed decomposition interfaces and the developmental
run_erd_interventional() are excluded from confirmatory
benchmarking. Their within-type-rate outcome did not identify the
abundance contribution described by the former direct/indirect labels.
The redesigned molecular-contribution estimand must first be
implemented, then tested under pure abundance, within-type expression,
combined, contextual, null and misspecified scenarios. No mediational
E-value is reported without a derivation valid for that exact
multivariate stochastic-composition estimand.
## 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] exposomeSC_0.99.1 rmarkdown_2.32
##
## loaded via a namespace (and not attached):
## [1] SummarizedExperiment_1.43.0 gtable_0.3.6
## [3] xfun_0.61 bslib_0.12.0
## [5] ggplot2_4.0.3 ggrepel_0.9.8
## [7] Biobase_2.73.2 lattice_0.23-1
## [9] vctrs_0.7.3 tools_4.6.1
## [11] generics_0.1.4 parallel_4.6.1
## [13] stats4_4.6.1 tibble_3.3.1
## [15] pkgconfig_2.0.3 Matrix_1.7-6
## [17] RColorBrewer_1.1-3 S7_0.2.2
## [19] S4Vectors_0.51.10 lifecycle_1.0.5
## [21] compiler_4.6.1 farver_2.1.2
## [23] statmod_1.5.2 ggforce_0.5.0
## [25] DESeq2_1.53.3 codetools_0.2-20
## [27] graphlayouts_1.2.5 Seqinfo_1.3.2
## [29] htmltools_0.5.9 sys_3.4.3
## [31] buildtools_1.0.0 sass_0.4.10
## [33] yaml_2.3.12 pillar_1.11.1
## [35] jquerylib_0.1.4 tidyr_1.3.2
## [37] MASS_7.3-66 BiocParallel_1.47.0
## [39] SingleCellExperiment_1.35.2 DelayedArray_0.39.7
## [41] cachem_1.1.0 limma_3.99.0
## [43] viridis_0.6.5 abind_1.4-8
## [45] tidyselect_1.2.1 locfit_1.5-9.12
## [47] digest_0.6.39 purrr_1.2.2
## [49] dplyr_1.2.1 labeling_0.4.3
## [51] maketools_1.3.2 splines_4.6.1
## [53] polyclip_1.10-7 fastmap_1.2.0
## [55] grid_4.6.1 cli_3.6.6
## [57] SparseArray_1.13.3 magrittr_2.0.5
## [59] S4Arrays_1.13.1 ggraph_2.2.2
## [61] tidygraph_1.3.1 edgeR_4.99.6
## [63] withr_3.0.3 scales_1.4.0
## [65] XVector_0.53.0 matrixStats_1.5.0
## [67] igraph_2.3.3 otel_0.2.0
## [69] gridExtra_2.3.1 memoise_2.0.1
## [71] evaluate_1.0.5 knitr_1.52
## [73] GenomicRanges_1.65.4 IRanges_2.47.5
## [75] viridisLite_0.4.3 rlang_1.3.0
## [77] Rcpp_1.1.2 glue_1.8.1
## [79] tweenr_2.0.3 BiocGenerics_0.59.12
## [81] jsonlite_2.0.0 R6_2.6.1
## [83] MatrixGenerics_1.25.0 coglasso_1.1.1