--- title: "Introduction to exposomeSC: a smoking-association workflow" author: "Cuiwei Gao" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true vignette: > %\VignetteIndexEntry{Introduction to exposomeSC: a smoking-association workflow} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Installation ```{r install, eval=FALSE} if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("exposomeSC") # Development version from GitHub: BiocManager::install("CuiweiG/exposomeSC") ``` # Introduction `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. # Unit of analysis and estimand 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. # Required data contract A defensible input has: - a raw, non-negative integer count assay; - unique cell and feature identifiers; - a donor identifier for every cell; - donor-constant exposure and baseline-covariate values; - cell-type annotations fixed before exposure testing; and - an auditable map from biological samples to donors when donors have repeated samples. 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. ```{r build} 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) exposureVariables(scee) ``` `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. # Canonical donor-pseudobulk model 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. ```{r primary-model} 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), ]) ``` The 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. # Multiplicity and output interpretation 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. # Repeated samples and composition 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. # Developmental decomposition interfaces `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. # Public source and reproducibility The HLCA source is the CC BY 4.0 CZ CELLxGENE collection associated with Sikkema *et al.* (2023): - collection: ; - exact H5AD: ; - publication: ; - count layer: `/raw/X`; and - H5AD SHA-256: `1cbdde1e513a31bdb6a3f70296e2d0612760d2af4032e239c94596ea87b44007`. 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. # Session information ```{r session-info} sessionInfo() ```