# Modeling continuous cell-level covariates

# Introduction

Since read counts are summed across cells in a pseudobulk approach, modeling continuous cell-level covariates also requires a collapsing step. Here we summarize the values of a variable from a set of cells using the mean, and store the value for each cell type. Including these variables in a regression formula uses the summarized values from the corresponding cell type.

We demonstrate this feature on a lightly modified analysis of PBMCs from 8 individuals stimulated with interferon-β (Kang, et al, 2018, Nature Biotech).

# Standard processing

Here is the code from the main vignette:

    library(dreamlet)
    library(muscat)
    library(ExperimentHub)
    library(scater)

    # Download data, specifying EH2259 for the Kang, et al study
    eh <- ExperimentHub()

    sce <- eh[["EH2259"]]

    # only keep singlet cells with sufficient reads
    sce <- sce[rowSums(counts(sce) > 0) > 0, ]
    sce <- sce[, colData(sce)$multiplets == "singlet"]

    # compute QC metrics
    qc <- perCellQCMetrics(sce)

    # remove cells with few or many detected genes
    ol <- isOutlier(metric = qc$detected, nmads = 2, log = TRUE)
    sce <- sce[, !ol]

    # set variable indicating stimulated (stim) or control (ctrl)
    sce$StimStatus <- sce$stim

In many datasets, continuous cell-level variables could be mapped reads, gene count, mitochondrial rate, etc. There are no continuous cell-level variables in this dataset, so we can simulate two from a normal distribution:

    sce$value1 <- rnorm(ncol(sce))
    sce$value2 <- rnorm(ncol(sce))

# Pseudobulk

Now compute the pseudobulk using standard code:

    sce$id <- paste0(sce$StimStatus, sce$ind)

    # Create pseudobulk
    pb <- aggregateToPseudoBulk(sce,
      assay = "counts",
      cluster_id = "cell",
      sample_id = "id",
      verbose = FALSE
    )

The means per variable, cell type, and sample are stored in the pseudobulk `SingleCellExperiment` object:

    metadata(pb)$aggr_means

    ## # A tibble: 128 × 5
    ## # Groups:   cell [8]
    ##    cell    id       cluster  value1   value2
    ##    <fct>   <fct>      <dbl>   <dbl>    <dbl>
    ##  1 B cells ctrl101     3.96 -0.0165 -0.0391
    ##  2 B cells ctrl1015    4.00  0.0177  0.0228
    ##  3 B cells ctrl1016    4    -0.0332 -0.00106
    ##  4 B cells ctrl1039    4.04 -0.175  -0.267
    ##  5 B cells ctrl107     4     0.162   0.120
    ##  6 B cells ctrl1244    4     0.0890 -0.0366
    ##  7 B cells ctrl1256    4.01  0.0272  0.121
    ##  8 B cells ctrl1488    4.02  0.0663  0.0294
    ##  9 B cells stim101     4.09  0.0172  0.224
    ## 10 B cells stim1015    4.06  0.0270  0.0192
    ## # ℹ 118 more rows

# Analysis

Including these variables in a regression formula uses the summarized values from the corresponding cell type. This happens behind the scenes, so the user doesn’t need to distinguish bewteen sample-level variables stored in `colData(pb)` and cell-level variables stored in `metadata(pb)$aggr_means`.

Variance partition and hypothesis testing proceeds as ususal:

    form <- ~ StimStatus + value1 + value2

    # Normalize and apply voom/voomWithDreamWeights
    res.proc <- processAssays(pb, form, min.count = 5)

    # run variance partitioning analysis
    vp.lst <- fitVarPart(res.proc, form)

    # Summarize variance fractions genome-wide for each cell type
    plotVarPart(vp.lst, label.angle = 60)

    # Differential expression analysis within each assay
    res.dl <- dreamlet(res.proc, form)

    # dreamlet results include coefficients for value1 and value2
    res.dl

    ## class: dreamletResult
    ## assays(8): B cells CD14+ Monocytes ... Megakaryocytes NK cells
    ## Genes:
    ##  min: 164
    ##  max: 5262
    ## details(7): assay n_retain ... n_errors error_initial
    ## coefNames(4): (Intercept) StimStatusstim value1 value2

# Details

A variable in `colData(sce)` is handled according to if the variable is

- continuous: the mean per donor/cell type is stored in `metadata(pb)$aggr_means`
- discrete
  - [constant within each donor/cell type] it is stored in `colData(pb)`
  - [varies within each donor/cell type] there is no good way to summarize it. The variable is dropped.

# Session Info

    ## R version 4.5.2 (2025-10-31)
    ## Platform: x86_64-pc-linux-gnu
    ## Running under: Ubuntu 24.04.4 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.26.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=C
    ##  [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] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] muscData_1.25.0             scater_1.39.2
    ##  [3] scuttle_1.21.0              ExperimentHub_3.1.0
    ##  [5] AnnotationHub_4.1.0         BiocFileCache_3.1.0
    ##  [7] dbplyr_2.5.2                muscat_1.25.2
    ##  [9] dreamlet_1.9.1              SingleCellExperiment_1.33.0
    ## [11] SummarizedExperiment_1.41.1 Biobase_2.71.0
    ## [13] GenomicRanges_1.63.1        Seqinfo_1.1.0
    ## [15] IRanges_2.45.0              S4Vectors_0.49.0
    ## [17] BiocGenerics_0.57.0         generics_0.1.4
    ## [19] MatrixGenerics_1.23.0       matrixStats_1.5.0
    ## [21] variancePartition_1.41.3    BiocParallel_1.45.0
    ## [23] limma_3.67.0                ggplot2_4.0.2
    ## [25] BiocStyle_2.39.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] bitops_1.0-9              httr_1.4.8
    ##   [3] RColorBrewer_1.1-3        doParallel_1.0.17
    ##   [5] Rgraphviz_2.55.0          numDeriv_2016.8-1.1
    ##   [7] tools_4.5.2               sctransform_0.4.3
    ##   [9] backports_1.5.0           utf8_1.2.6
    ##  [11] R6_2.6.1                  metafor_4.8-0
    ##  [13] mgcv_1.9-4                GetoptLong_1.1.0
    ##  [15] withr_3.0.2               prettyunits_1.2.0
    ##  [17] gridExtra_2.3             fdrtool_1.2.18
    ##  [19] cli_3.6.5                 sandwich_3.1-1
    ##  [21] labeling_0.4.3            slam_0.1-55
    ##  [23] sass_0.4.10               KEGGgraph_1.71.0
    ##  [25] SQUAREM_2021.1            mvtnorm_1.3-3
    ##  [27] S7_0.2.1                  blme_1.0-7
    ##  [29] mixsqp_0.3-54             zenith_1.13.0
    ##  [31] parallelly_1.46.1         invgamma_1.2
    ##  [33] RSQLite_2.4.6             shape_1.4.6.1
    ##  [35] gtools_3.9.5              dplyr_1.2.0
    ##  [37] Matrix_1.7-4              metadat_1.4-0
    ##  [39] ggbeeswarm_0.7.3          abind_1.4-8
    ##  [41] lifecycle_1.0.5           yaml_2.3.12
    ##  [43] edgeR_4.9.2               mathjaxr_2.0-0
    ##  [45] gplots_3.3.0              SparseArray_1.11.10
    ##  [47] grid_4.5.2                blob_1.3.0
    ##  [49] crayon_1.5.3              lattice_0.22-9
    ##  [51] beachmat_2.27.2           msigdbr_25.1.1
    ##  [53] annotate_1.89.0           KEGGREST_1.51.1
    ##  [55] sys_3.4.3                 maketools_1.3.2
    ##  [57] pillar_1.11.1             knitr_1.51
    ##  [59] ComplexHeatmap_2.27.1     rjson_0.2.23
    ##  [61] boot_1.3-32               corpcor_1.6.10
    ##  [63] future.apply_1.20.2       codetools_0.2-20
    ##  [65] glue_1.8.0                data.table_1.18.2.1
    ##  [67] vctrs_0.7.1               png_0.1-8
    ##  [69] Rdpack_2.6.6              gtable_0.3.6
    ##  [71] assertthat_0.2.1          cachem_1.1.0
    ##  [73] zigg_0.0.2                xfun_0.56
    ##  [75] rbibutils_2.4.1           S4Arrays_1.11.1
    ##  [77] Rfast_2.1.5.2             reformulas_0.4.4
    ##  [79] iterators_1.0.14          statmod_1.5.1
    ##  [81] nlme_3.1-168              pbkrtest_0.5.5
    ##  [83] bit64_4.6.0-1             filelock_1.0.3
    ##  [85] progress_1.2.3            EnvStats_3.1.0
    ##  [87] bslib_0.10.0              TMB_1.9.19
    ##  [89] irlba_2.3.7               vipor_0.4.7
    ##  [91] KernSmooth_2.23-26        colorspace_2.1-2
    ##  [93] rmeta_3.0                 DBI_1.2.3
    ##  [95] DESeq2_1.51.6             tidyselect_1.2.1
    ##  [97] bit_4.6.0                 compiler_4.5.2
    ##  [99] curl_7.0.0                httr2_1.2.2
    ## [101] graph_1.89.1              BiocNeighbors_2.5.4
    ## [103] DelayedArray_0.37.0       scales_1.4.0
    ## [105] caTools_1.18.3            remaCor_0.0.20
    ## [107] rappdirs_0.3.4            stringr_1.6.0
    ## [109] digest_0.6.39             minqa_1.2.8
    ## [111] rmarkdown_2.30            aod_1.3.3
    ## [113] XVector_0.51.0            RhpcBLASctl_0.23-42
    ## [115] htmltools_0.5.9           pkgconfig_2.0.3
    ## [117] lme4_1.1-38               sparseMatrixStats_1.23.0
    ## [119] lpsymphony_1.39.0         mashr_0.2.79
    ## [121] fastmap_1.2.0             rlang_1.1.7
    ## [123] GlobalOptions_0.1.3       DelayedMatrixStats_1.33.0
    ## [125] farver_2.1.2              jquerylib_0.1.4
    ## [127] IHW_1.39.0                zoo_1.8-15
    ## [129] jsonlite_2.0.0            BiocSingular_1.27.1
    ## [131] RCurl_1.98-1.17           magrittr_2.0.4
    ## [133] Rcpp_1.1.1                babelgene_22.9
    ## [135] viridis_0.6.5             EnrichmentBrowser_2.41.0
    ## [137] stringi_1.8.7             MASS_7.3-65
    ## [139] plyr_1.8.9                listenv_0.10.0
    ## [141] parallel_4.5.2            ggrepel_0.9.6
    ## [143] Biostrings_2.79.4         splines_4.5.2
    ## [145] hms_1.1.4                 circlize_0.4.17
    ## [147] locfit_1.5-9.12           buildtools_1.0.0
    ## [149] reshape2_1.4.5            ScaledMatrix_1.19.0
    ## [151] BiocVersion_3.23.1        XML_3.99-0.22
    ## [153] evaluate_1.0.5            RcppParallel_5.1.11-1
    ## [155] BiocManager_1.30.27       nloptr_2.2.1
    ## [157] foreach_1.5.2             tidyr_1.3.2
    ## [159] purrr_1.2.1               future_1.69.0
    ## [161] clue_0.3-67               scattermore_1.2
    ## [163] ashr_2.2-63               rsvd_1.0.5
    ## [165] broom_1.0.12              xtable_1.8-8
    ## [167] fANCOVA_0.6-1             viridisLite_0.4.3
    ## [169] truncnorm_1.0-9           tibble_3.3.1
    ## [171] lmerTest_3.2-0            glmmTMB_1.1.14
    ## [173] memoise_2.0.1             beeswarm_0.4.0
    ## [175] AnnotationDbi_1.73.0      cluster_2.1.8.2
    ## [177] globals_0.19.0            GSEABase_1.73.0