exposomeSC estimates donor-level exposure associations
with relative gene expression within a cell type. Cells are measurements
nested within a donor, not independent replicates. The implemented
default workflow aggregates raw counts by donor and cell type and fits a
robust edgeR quasi-likelihood model.
The current real-data application is ever- versus never-smoking in the HLCA. It is an adjusted observational contrast, not a causal effect.
sce must contain a raw integer counts assay
and colData columns identifying donor and cell type.
donor_design must be numeric, have one row per donor and
use the same donor identifiers as its row names.
library(exposomeSC)
# A small simulated demonstration object. Cells within a donor share a
# donor-level offset, which is the dependence donor pseudobulk exists to
# handle; Gene1 is given a 2.5-fold increase in exposed AT2 cells. The full
# analysis uses the pinned HLCA inputs described below, not these counts.
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
)
)
donor_design <- cbind(smoking_ever = smoking_ever)
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"For a multistudy atlas, construct the study indicators and nonlinear age terms before building the object. Never derive those terms, the cell-type universe or the feature family from observed exposure effects.
fit <- run_sc_exwas(
scee,
exposure = "smoking_ever",
celltype_col = "ann_level_3",
sample_col = "donor_id",
min_cells = 10L,
min_donors = 10L,
min_group_donors = 5L,
method = "edgeR",
filter_method = "fixed_support",
filter_min_cpm = 1,
filter_min_donors = 5L,
filter_min_total_count = 15L,
robust = TRUE
)The model retains the complete target-cell-type library size for offsets and TMM normalisation even when the reported feature family is restricted. The fixed-support filter depends on counts and donor support, not smoking labels.
required_columns <- c(
"gene", "celltype", "log2FC", "se", "statistic", "pvalue",
"padj", "padj_global", "n_donors"
)
stopifnot(all(required_columns %in% colnames(fit)))
# Sorting is descriptive; the inferential family was fixed before this step.
head(fit[order(fit$padj_global, fit$pvalue), required_columns])
#> DataFrame with 6 rows and 9 columns
#> gene celltype log2FC se statistic pvalue padj
#> <character> <character> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 Gene1 AT2 1.010578 NA 4.06064 8.76098e-05 0.00105132
#> 2 Gene3 AT2 -0.564568 NA -2.09356 3.84046e-02 0.23042741
#> 3 Gene10 Macrophage 0.469246 NA 1.65573 1.00390e-01 0.68464830
#> 4 Gene2 Macrophage -0.456994 NA -1.57063 1.18901e-01 0.68464830
#> 5 Gene11 Macrophage 0.367593 NA 1.37672 1.71162e-01 0.68464830
#> 6 Gene2 AT2 -0.354718 NA -1.31321 1.91618e-01 0.61964129
#> padj_global n_donors
#> <numeric> <integer>
#> 1 0.00210264 12
#> 2 0.46085482 12
#> 3 0.71340858 12
#> 4 0.71340858 12
#> 5 0.76647236 12
#> 6 0.76647236 12
# What the backend reports in place of a standard error
S4Vectors::metadata(fit)$parameters[c("se_method", "statistic_type")]
#> $se_method
#> [1] "not_available_edgeR_QL"
#>
#> $statistic_type
#> [1] "signed_sqrt_qlf"For the default edgeR backend, se is deliberately
numeric NA because the quasi-likelihood F test does not
supply a coefficient standard error. statistic is the
signed square root of that QL F statistic—a directional ranking
quantity, not a Wald z statistic and not a basis for a coefficient
confidence interval. The last line of the chunk above prints the two
parameters that record this; check them before any downstream operation
that uses uncertainty.
In the full HLCA analysis, feature annotation defines a primary host-expression family and separate receptor/residual families after fitting. The reporting adjustment is recomputed across each complete prespecified gene-by-cell-type family, so do not infer that family from this abbreviated example.
Gene1 in AT2 cells heads the table, which is the one effect the simulation plants. The counts are simulated and the cohort is a dozen donors, so this shows that the workflow recovers a planted effect and carries no biological conclusion.
Reasonable analyses should examine annotation resolution, target-cell support, age specification, missing BMI, within-study exposure overlap, sex strata, current and former smoking separately, and leave-one-study-out fits. Donors with repeated biological samples require a sample-pseudobulk sensitivity with donor blocking. Influence is assessed by refitting after donor or study deletion, not by selectively removing observations.
Relative cell-type composition should be modelled jointly at donor level. A within-type expression coefficient does not estimate absolute abundance, and a captured-cell proportion is affected by tissue sampling, dissociation, capture, quality control and annotation.
run_decomposed_exwas() has been removed, and
run_erd_interventional() is an unexported developmental
function that is not for confirmatory interpretation. The former
direct/indirect labels and E-values do not establish a causal pathway or
an abundance contribution. A replacement estimand is under
development.
The exact HLCA asset, /raw/X layer and SHA-256 value are
listed in the package README. Reportable outputs must include matching
input and code manifests, parameters, design diagnostics, session
information and output hashes.
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