Installation

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("exposomeSC")

# Development version from GitHub:
BiocManager::install("CuiweiG/exposomeSC")

Introduction

The network module estimates cell-type-resolved conditional-association structures that combine transcript pseudobulk with donor-level molecular measurements such as metabolomics. It is exploratory. An edge is a regularised conditional association under the fitted model—not a biochemical interaction or a causal relation.

Network analysis answers a different question from sc-ExWAS. sc-ExWAS tests a prespecified marginal exposure coefficient for each feature and cell type. Network analysis describes multivariate dependence after preprocessing and covariate adjustment. Selection and tuning uncertainty, sample size and external stability must be reported separately.

The fitted network object records its matrices, feature annotations, method and provenance-related metadata in an explicit S4 contract.

library(exposomeSC)
methods::getSlots("CelltypeNetworkResult")
#> precision_matrix adjacency_matrix stability_scores        node_info 
#>         "matrix"         "matrix"         "matrix"      "DataFrame" 
#>         celltype           method         metadata 
#>      "character"      "character"           "list"

Input alignment

The transcript and metabolite data must refer to the same donors. Raw single-cell counts are aggregated within donor and cell type before network estimation; donor-level metabolite rows are then aligned by identifier. Do not match rows by position.

The examples below use a small simulated data set so that the vignette runs in seconds: 24 donors, monocytes and B cells, 12 genes and 4 donor-level metabolites. In monocytes four genes move with the exposure, and one metabolite tracks it.

has_coglasso <- requireNamespace("coglasso", quietly = TRUE)
has_plot_packages <- all(vapply(c("ggplot2", "igraph", "ggraph", "tidygraph"),
    requireNamespace, logical(1), quietly = TRUE))

set.seed(2026)
donors <- sprintf("D%02d", 1:24)
celltypes <- c("Monocyte", "B cell")
cell_donor <- rep(rep(donors, each = 20), times = length(celltypes))
cell_type <- rep(celltypes, each = length(donors) * 20)
pm25 <- stats::setNames(stats::rnorm(length(donors)), donors)

counts <- matrix(stats::rpois(12 * length(cell_donor), 15), nrow = 12,
    dimnames = list(paste0("G", 1:12), paste0("c", seq_along(cell_donor))))
in_mono <- cell_type == "Monocyte"
for (gene in 1:4) {
    effect <- c(0.6, 0.6, -0.6, -0.6)[gene]
    counts[gene, in_mono] <- stats::rpois(sum(in_mono),
        15 * exp(effect * pm25[cell_donor[in_mono]]))
}
sce <- SingleCellExperiment::SingleCellExperiment(
    assays = list(counts = counts),
    colData = S4Vectors::DataFrame(cell_id = colnames(counts),
        donor_id = cell_donor, cell_type = cell_type))

donor_data <- cbind(PM2.5 = pm25,
    age = round(stats::runif(length(donors), 30, 70)),
    sex = rep(0:1, length.out = length(donors)))
rownames(donor_data) <- donors
scee <- build_scee(sce, donor_data, sample_col = "donor_id")

metabolites <- matrix(stats::rnorm(length(donors) * 4), ncol = 4,
    dimnames = list(donors, paste0("M", 1:4)))
metabolites[, "M1"] <- metabolites[, "M1"] + 1.2 * pm25

mono_net <- run_celltype_network(
    scee,
    metabolites = metabolites,
    celltype = "Monocyte",
    sample_col = "donor_id",
    covariates = c("age", "sex"),
    method = "coglasso",
    stability = TRUE,
    nlambda_w = 3L,
    nlambda_b = 3L,
    rep_num = 5L
)

The exposure is deliberately left out of this model. The four exposure genes and M1 share its signal, so the network is expected to connect them; any other edge is a false positive, which 24 donors can produce even under stability selection. Adding exposure = "PM2.5" regresses the exposure out of both matrices before estimation, and that module is no longer expected: the network then describes dependence the exposure does not explain. Which of the two is the scientific question has to be settled before fitting.

The tuning grid and the number of XStARS subsamples are deliberately small here; the defaults are a 15 by 15 grid and 20 subsamples. These calls do not reproduce the HLCA primary model. In a multistudy application, study and other prespecified design terms must be included consistently, and the feature-selection rule must be recorded before network fitting.

Collaborative graphical lasso

The coglasso backend permits different penalties for within-omic and between-omic precision-matrix blocks. Penalty selection is based on resampling stability. The complete tuning grid, resample identifiers, edge-selection frequencies and final penalty must be retained; a network plot alone is not an adequate result record.

mono_net
#> CelltypeNetworkResult for 'Monocyte'
#>   Method: coglasso
#>   Nodes: 16 (12 transcripts, 4 metabolites)
#>   Edges: 5 (5 cross-omic)
#>   Stability selection: yes
#>   Donors: 24
plot_stability_surface(mono_net)

plot_celltype_network(
    mono_net,
    layout = "fr",
    color_by = "omic_layer",
    highlight_cross_omic = TRUE
)

Comparing cell types

Cell-type comparisons require aligned features, compatible preprocessing and adequate donor support in every network. Apparent differences can result from unequal measurement error or selection stability rather than a true change in conditional association.

bcell_net <- run_celltype_network(
    scee,
    metabolites = metabolites,
    celltype = "B cell",
    sample_col = "donor_id",
    covariates = c("age", "sex"),
    method = "coglasso",
    stability = TRUE,
    nlambda_w = 3L,
    nlambda_b = 3L,
    rep_num = 5L
)

comparison <- run_comparative_network(
    mono_net,
    bcell_net,
    method = "descriptive"
)
comparison
#> NetworkComparison
#>   Cell types compared: Monocyte, B cell
#>   Shared edges: 0
#>   Differential edges: 7
#>   Unique to Monocyte: 5
#>   Unique to B cell: 2

plot_network_comparison(comparison)

This comparison is descriptive and reports no p-values. The Fisher-Z mode is experimental: networks estimated from the same donors are not independent samples, and penalised partial correlations do not follow the distribution the test assumes. Edge-presence or weight differences are not evidence of differential conditional association.

Exposure-restricted networks

run_exposure_network() uses an exposure-association stage followed by network estimation. This is a two-stage exploratory procedure. Selecting and estimating on the same donors can exaggerate stability and effect size; nested resampling, sample splitting or independent validation is required for confirmatory use.

exposure_net <- run_exposure_network(
    scee,
    metabolites,
    celltype = "Monocyte",
    exposure = "PM2.5",
    sample_col = "donor_id",
    covariates = c("age", "sex"),
    selection_method = "univariate",
    fdr_threshold = 0.05,
    network_method = "block_glasso"
)
exposure_net
#> CelltypeNetworkResult for 'Monocyte'
#>   Method: exposure_driven
#>   Nodes: 5 (4 transcripts, 1 metabolites)
#>   Edges: 6 (0 cross-omic)
#>   Donors: 24

Longitudinal networks

Time-point-specific networks can be organised with run_temporal_network(). Repeated measurements must remain linked to their donor, and an apparent emerging or disappearing edge is descriptive unless its uncertainty has been assessed under an appropriate longitudinal resampling scheme.

metabolites_follow_up <- metabolites
metabolites_follow_up[, "M2"] <- metabolites_follow_up[, "M2"] + 1.5 * pm25

temporal <- run_temporal_network(
    scee_list = list(baseline = scee, follow_up = scee),
    metabolites_list = list(
        baseline = metabolites,
        follow_up = metabolites_follow_up
    ),
    celltype = "Monocyte",
    sample_col = "donor_id",
    covariates = c("age", "sex"),
    method = "block_glasso",
    stability = FALSE
)
temporal
#> TemporalNetwork
#>   Time points: baseline -> follow_up
#>   Networks: 2
#>   Edge dynamics tracked: 6 edges

plot_temporal_dynamics(temporal)

The follow-up here reuses the baseline cells with one changed metabolite, so the example shows the mechanics of edge tracking rather than a longitudinal design.

Mediation terminology

run_network_mediation() is a model-based exploratory utility. Network topology does not establish temporal ordering, eliminate unmeasured confounding or identify a causal indirect effect. Results must be described as conditional model associations unless the study design and identification assumptions independently support stronger language.

Composition and provenance

Building a network within a cell type removes variation in broad cell-type mixture from the expression outcome, but it does not eliminate subtype mixture, sampling, dissociation, capture, annotation or residual confounding. The network model does not supply evidence for a composition-mediated pathway.

Every reported network should retain input hashes, donor inclusion, feature definitions, preprocessing, covariate design, tuning grid, random seeds, resampling output, software versions and the exact object used for plotting. The examples above run on simulated data and report no result from a real cohort.

References

Albanese A, Kohlen W, Behrouzi P (2024). Collaborative graphical lasso. arXiv:2403.18602.
Liu H, Roeder K, Wasserman L (2010). Stability Approach to Regularization Selection (StARS) for High Dimensional Graphical Models. Advances in Neural Information Processing Systems 23.
Meinshausen N, Bühlmann P (2010). Stability selection. Journal of the Royal Statistical Society: Series B 72:417–473.

Session information

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] 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

Exploratory cross-omic network inference with exposomeSC