exposomeSC relates donor-level exposures to expression
within annotated single-cell populations. The package stores counts,
cell metadata, donor exposures and exposure metadata in a validated
SingleCellExposomeExperiment, then fits models at the
biological-replicate level.
The motivating analysis contrasts ever- with never-smoking in the integrated Human Lung Cell Atlas (HLCA). It is an observational association analysis. Smoking history is not a randomised intervention, and the available ever-smoking indicator combines active and former exposure. We therefore avoid causal language.
Cells inherit their donor’s exposure. Treating those cells as independent replicates would underestimate uncertainty. For each cell type, the canonical workflow sums raw counts within donor, producing one pseudobulk column per eligible donor. The target is the exposure-associated rate ratio for relative expression within that cell type.
This target differs from cell-type abundance. If smoking changes only the number of captured cells of a type, both its gene counts and complete library size may change together while its within-type expression rate remains stable. Relative composition is consequently analysed as a separate multivariate outcome.
A defensible input has:
The rows of the donor design matrix must match donor identifiers exactly. A categorical study variable is represented by prespecified numeric indicator columns before the object is built.
library(exposomeSC)
# A small simulated object makes the vignette executable. Cells within a
# donor share a donor-level offset, and Gene1 is raised 2.5-fold in exposed
# AT2 cells. It illustrates the data contract only; it is not used for the
# HLCA results.
set.seed(20260916L)
donor_ids <- paste0("D", seq_len(12L))
smoking_ever <- rep(c(0, 1), each = 6L)
cell_donor <- rep(donor_ids, each = 24L)
cell_type <- rep(
rep(c("AT2", "Macrophage"), each = 12L),
times = length(donor_ids)
)
gene_ids <- paste0("Gene", seq_len(12L))
cell_ids <- paste0("cell", seq_along(cell_donor))
donor_offset <- matrix(
stats::rnorm(length(gene_ids) * length(donor_ids), sd = 0.35),
nrow = length(gene_ids),
dimnames = list(gene_ids, donor_ids)
)
exposed_at2 <- smoking_ever[match(cell_donor, donor_ids)] == 1 &
cell_type == "AT2"
cell_mean <- exp(log(25) + donor_offset[, cell_donor])
cell_mean[1L, ] <- cell_mean[1L, ] * ifelse(exposed_at2, 2.5, 1)
counts <- matrix(
stats::rnbinom(length(cell_mean), mu = as.vector(cell_mean), size = 8),
nrow = length(gene_ids),
dimnames = list(gene_ids, cell_ids)
)
storage.mode(counts) <- "integer"
sce <- SingleCellExperiment::SingleCellExperiment(
assays = list(counts = counts),
colData = S4Vectors::DataFrame(
cell_id = cell_ids,
donor_id = cell_donor,
ann_level_3 = cell_type
)
)
age_scaled <- as.numeric(scale(c(
45, 52, 61, 48, 57, 66, 50, 59, 68, 54, 63, 72
)))
donor_design <- cbind(
smoking_ever = smoking_ever,
age_scaled = age_scaled
)
rownames(donor_design) <- donor_ids
scee <- build_scee(
sce,
donor_design,
sample_col = "donor_id"
)
methods::validObject(scee)
#> [1] TRUE
exposureVariables(scee)
#> [1] "smoking_ever" "age_scaled"exposureVariables() lists every column of the donor
design, so the age term appears alongside the exposure. Which column is
the exposure, and which are adjustments, is set in the call to
run_sc_exwas() below.
The implemented primary backend uses robust edgeR quasi-likelihood
inference. Expression support is fixed independently of the smoking
labels, and the complete target-cell-type transcriptome remains
available for library-size and normalisation calculations after
filtering. Adjustment terms are supplied through
covariates, as columns of the donor design; a nonlinear age
term or a multistudy indicator set is built before the object is
created, never derived from the observed exposure effects.
fit <- run_sc_exwas(
scee,
exposure = "smoking_ever",
celltype_col = "ann_level_3",
sample_col = "donor_id",
covariates = "age_scaled",
min_cells = 10L,
min_donors = 10L,
min_group_donors = 5L,
method = "edgeR",
filter_genes = TRUE,
filter_method = "fixed_support",
filter_min_cpm = 1,
filter_min_donors = 5L,
filter_min_total_count = 15L,
robust = TRUE
)
head(fit[order(fit$padj_global, fit$pvalue), ])
#> DataFrame with 6 rows and 17 columns
#> gene celltype log2FC se statistic pvalue
#> <character> <character> <numeric> <numeric> <numeric> <numeric>
#> 1 Gene1 AT2 0.960521 NA 3.244107 0.00156839
#> 2 Gene3 AT2 -0.575408 NA -1.984446 0.04974257
#> 3 Gene2 Macrophage -0.483770 NA -1.494824 0.13787667
#> 4 Gene11 Macrophage 0.457719 NA 1.474802 0.14317513
#> 5 Gene10 Macrophage 0.381676 NA 1.179911 0.24062813
#> 6 Gene12 AT2 -0.240536 NA -0.804872 0.42266176
#> pvalue_underflow_clamped padj logCPM exposure n_donors
#> <logical> <numeric> <numeric> <character> <integer>
#> 1 FALSE 0.0188207 17.0812 smoking_ever 12
#> 2 FALSE 0.2984554 16.3350 smoking_ever 12
#> 3 FALSE 0.8590508 16.4468 smoking_ever 12
#> 4 FALSE 0.8590508 16.3278 smoking_ever 12
#> 5 FALSE 0.9625125 16.4211 smoking_ever 12
#> 6 FALSE 0.8797329 16.0188 smoking_ever 12
#> n_unexposed n_exposed min_cells median_cells method padj_global
#> <integer> <integer> <integer> <numeric> <character> <numeric>
#> 1 6 6 12 12 edgeR_robust_QL 0.0376414
#> 2 6 6 12 12 edgeR_robust_QL 0.5969108
#> 3 6 6 12 12 edgeR_robust_QL 0.8590508
#> 4 6 6 12 12 edgeR_robust_QL 0.8590508
#> 5 6 6 12 12 edgeR_robust_QL 0.9798704
#> 6 6 6 12 12 edgeR_robust_QL 0.9798704The call above is deliberately short. An analysis meant for reporting also builds its adjustment terms, records the input and code versions it ran on, inspects the design for rank and separation problems, and states the feature family before the model is fitted.
The primary feature family comprises the annotation-defined, exposure-blind support-filtered host-expression transcriptome. Benjamini–Hochberg adjustment is applied across its complete realised tested gene-by-cell-type family. The host category is externally annotated; the realised tested universe also depends on the donor-support rule, which is fixed before fitting. Immune-receptor and residual-annotation features are retained as separate secondary or exploratory audits; their adjusted values are not merged into the primary family after looking at the results.
The returned data include the effect estimate, backend-specific test
statistic, unadjusted p-value, within-cell-type adjustment,
cross-cell-type adjustment and donor-support diagnostics. edgeR
quasi-likelihood fits report the coefficient standard-error field as
numeric NA: glmQLFTest() provides a QL F test,
not a coefficient standard error. The signed square-root QL statistic
preserves direction for ranking but is neither a Wald statistic nor a
basis for a coefficient confidence interval. DESeq2 and voom-dream
retain their own valid Wald and moderated-t uncertainty semantics. A
small p-value is insufficient on its own. We also require a stable
feature definition, adequate support, an identifiable design,
sensitivity to annotation and sampling choices, empirical-null
calibration, and influence assessment.
The object built above is synthetic, so the fit illustrates the interface and carries no biological conclusion.
The donor-aggregated model is the primary expression estimand. A sensitivity analysis that keeps repeated biological samples apart, for example with donor blocking in voom/limma, is not part of the package; it would be run separately and is not a second discovery screen.
The composition analysis operates on the supported
ann_level_3 simplex. It uses Jeffreys-regularised log
ratios, donor-level aggregation and a study-adjusted multivariate test.
The resulting coefficients concern relative captured-cell
composition—not absolute lung-cell abundance—and no mediation
interpretation is made.
run_decomposed_exwas() has been removed.
run_erd_interventional() is kept as an unexported
developmental function and is not for confirmatory interpretation. The
earlier within-type-rate outcome cannot identify the mechanical
contribution of cell-type abundance to a tissue mixture, so former
direct/indirect labels and ERD E-values must not be used for causal or
abundance claims.
A replacement estimand remains developmental until it is implemented, calibrated under appropriate data-generating mechanisms and assessed independently.
The HLCA source is the CC BY 4.0 CZ CELLxGENE collection associated with Sikkema et al. (2023):
/raw/X; and1cbdde1e513a31bdb6a3f70296e2d0612760d2af4032e239c94596ea87b44007.Every reportable result must retain the verified input hashes, ordered code manifest, model parameters, design diagnostics, session information and output hashes. Figures and tables should read the stored result objects rather than have numerical values transcribed into display code.
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] exposomeSC_0.99.1 rmarkdown_2.32
#>
#> 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] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [7] GenomicRanges_1.65.4 jquerylib_0.1.4
#> [9] splines_4.6.1 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 R6_2.6.1
#> [17] XVector_0.53.0 S4Arrays_1.13.1
#> [19] generics_0.1.4 knitr_1.52
#> [21] BiocGenerics_0.59.12 DelayedArray_0.39.7
#> [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.3
#> [33] cli_3.6.6 locfit_1.5-9.12
#> [35] digest_0.6.39 grid_4.6.1
#> [37] edgeR_4.99.6 lifecycle_1.0.5
#> [39] S4Vectors_0.51.10 SingleCellExperiment_1.35.2
#> [41] evaluate_1.0.5 buildtools_1.0.0
#> [43] abind_1.4-8 stats4_4.6.1
#> [45] matrixStats_1.5.0 tools_4.6.1
#> [47] htmltools_0.5.9