Modeling continuous cell-level covariates
Collapse using mean value for pseudobulk data
Developed by Gabriel Hoffman
Run on 2024-11-20 14:13:13
Source:vignettes/cell_covs.Rmd
cell_covs.RmdIntroduction
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$stimIn 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:
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.127 0.0241
## 2 B cells ctrl1015 4.00 -0.0539 -0.00302
## 3 B cells ctrl1016 4 0.00571 0.0366
## 4 B cells ctrl1039 4.04 -0.102 -0.0970
## 5 B cells ctrl107 4 0.0438 0.0163
## 6 B cells ctrl1244 4 -0.0928 0.138
## 7 B cells ctrl1256 4.01 0.0344 -0.0432
## 8 B cells ctrl1488 4.02 -0.0599 -0.00684
## 9 B cells stim101 4.09 0.128 0.0350
## 10 B cells stim1015 4.06 0.0267 0.0681
## # ℹ 118 more rowsAnalysis
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 value2Details
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.
- [constant within each donor/cell type] it is stored in
Session Info
## R version 4.3.0 (2023-04-21)
## Platform: x86_64-apple-darwin22.4.0 (64-bit)
## Running under: macOS Ventura 13.5
##
## Matrix products: default
## BLAS: /Users/gabrielhoffman/prog/R-4.3.0/lib/libRblas.dylib
## LAPACK: /usr/local/Cellar/r/4.3.0_1/lib/R/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] muscData_1.14.0 scater_1.28.0
## [3] scuttle_1.10.1 ExperimentHub_2.8.1
## [5] AnnotationHub_3.8.0 BiocFileCache_2.8.0
## [7] dbplyr_2.3.2 muscat_1.14.0
## [9] dreamlet_1.1.5 SingleCellExperiment_1.22.0
## [11] SummarizedExperiment_1.30.1 Biobase_2.60.0
## [13] GenomicRanges_1.52.0 GenomeInfoDb_1.36.1
## [15] IRanges_2.34.1 S4Vectors_0.38.1
## [17] BiocGenerics_0.46.0 MatrixGenerics_1.12.0
## [19] matrixStats_1.0.0 variancePartition_1.33.2
## [21] BiocParallel_1.34.2 limma_3.56.2
## [23] ggplot2_3.4.4 BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.2 bitops_1.0-7
## [3] httr_1.4.6 RColorBrewer_1.1-3
## [5] doParallel_1.0.17 Rgraphviz_2.44.0
## [7] numDeriv_2016.8-1.1 tools_4.3.0
## [9] sctransform_0.3.5 backports_1.4.1
## [11] utf8_1.2.3 R6_2.5.1
## [13] GetoptLong_1.0.5 withr_2.5.0
## [15] prettyunits_1.1.1 gridExtra_2.3
## [17] cli_3.6.1 textshaping_0.3.6
## [19] sandwich_3.0-2 labeling_0.4.2
## [21] sass_0.4.6 KEGGgraph_1.60.0
## [23] SQUAREM_2021.1 mvtnorm_1.2-2
## [25] blme_1.0-5 pkgdown_2.0.7
## [27] mixsqp_0.3-48 systemfonts_1.0.4
## [29] zenith_1.4.1 parallelly_1.36.0
## [31] invgamma_1.1 RSQLite_2.3.1
## [33] generics_0.1.3 shape_1.4.6
## [35] gtools_3.9.4 dplyr_1.1.2
## [37] Matrix_1.5-4.1 ggbeeswarm_0.7.2
## [39] fansi_1.0.4 abind_1.4-5
## [41] lifecycle_1.0.3 multcomp_1.4-23
## [43] yaml_2.3.7 edgeR_3.42.4
## [45] gplots_3.1.3 grid_4.3.0
## [47] blob_1.2.4 promises_1.2.0.1
## [49] crayon_1.5.2 lattice_0.21-8
## [51] beachmat_2.16.0 msigdbr_7.5.1
## [53] annotate_1.78.0 KEGGREST_1.40.0
## [55] pillar_1.9.0 knitr_1.43
## [57] ComplexHeatmap_2.16.0 rjson_0.2.21
## [59] boot_1.3-28.1 estimability_1.4.1
## [61] corpcor_1.6.10 future.apply_1.11.0
## [63] codetools_0.2-19 glue_1.6.2
## [65] data.table_1.14.8 vctrs_0.6.3
## [67] png_0.1-8 Rdpack_2.4
## [69] gtable_0.3.3 assertthat_0.2.1
## [71] cachem_1.0.8 xfun_0.39
## [73] mime_0.12 rbibutils_2.2.13
## [75] S4Arrays_1.0.4 Rfast_2.0.7
## [77] coda_0.19-4 survival_3.5-5
## [79] iterators_1.0.14 ellipsis_0.3.2
## [81] interactiveDisplayBase_1.38.0 TH.data_1.1-2
## [83] nlme_3.1-162 pbkrtest_0.5.2
## [85] bit64_4.0.5 filelock_1.0.2
## [87] progress_1.2.2 EnvStats_2.7.0
## [89] rprojroot_2.0.3 bslib_0.4.2
## [91] TMB_1.9.4 irlba_2.3.5.1
## [93] vipor_0.4.5 KernSmooth_2.23-21
## [95] colorspace_2.1-0 rmeta_3.0
## [97] DBI_1.1.3 DESeq2_1.40.1
## [99] tidyselect_1.2.0 emmeans_1.8.7
## [101] curl_5.0.0 bit_4.0.5
## [103] compiler_4.3.0 graph_1.78.0
## [105] BiocNeighbors_1.18.0 desc_1.4.2
## [107] DelayedArray_0.26.3 bookdown_0.34
## [109] scales_1.2.1 caTools_1.18.2
## [111] remaCor_0.0.17 rappdirs_0.3.3
## [113] stringr_1.5.0 digest_0.6.33
## [115] minqa_1.2.5 rmarkdown_2.22
## [117] aod_1.3.2 XVector_0.40.0
## [119] RhpcBLASctl_0.23-42 htmltools_0.5.5
## [121] pkgconfig_2.0.3 lme4_1.1-34
## [123] sparseMatrixStats_1.12.0 highr_0.10
## [125] mashr_0.2.69 fastmap_1.1.1
## [127] rlang_1.1.1 GlobalOptions_0.1.2
## [129] shiny_1.7.4 DelayedMatrixStats_1.22.0
## [131] farver_2.1.1 jquerylib_0.1.4
## [133] zoo_1.8-12 jsonlite_1.8.5
## [135] BiocSingular_1.16.0 RCurl_1.98-1.12
## [137] magrittr_2.0.3 GenomeInfoDbData_1.2.10
## [139] munsell_0.5.0 Rcpp_1.0.11
## [141] babelgene_22.9 viridis_0.6.3
## [143] EnrichmentBrowser_2.30.1 RcppZiggurat_0.1.6
## [145] stringi_1.7.12 zlibbioc_1.46.0
## [147] MASS_7.3-60 plyr_1.8.8
## [149] listenv_0.9.0 parallel_4.3.0
## [151] ggrepel_0.9.3 Biostrings_2.68.1
## [153] splines_4.3.0 hms_1.1.3
## [155] circlize_0.4.15 locfit_1.5-9.7
## [157] reshape2_1.4.4 ScaledMatrix_1.8.1
## [159] BiocVersion_3.17.1 XML_3.99-0.14
## [161] evaluate_0.21 BiocManager_1.30.20
## [163] httpuv_1.6.11 nloptr_2.0.3
## [165] foreach_1.5.2 tidyr_1.3.0
## [167] purrr_1.0.2 future_1.32.0
## [169] clue_0.3-64 scattermore_1.1
## [171] ashr_2.2-54 rsvd_1.0.5
## [173] broom_1.0.5 xtable_1.8-4
## [175] fANCOVA_0.6-1 later_1.3.1
## [177] viridisLite_0.4.2 ragg_1.2.5
## [179] truncnorm_1.0-9 tibble_3.2.1
## [181] lmerTest_3.1-3 glmmTMB_1.1.7
## [183] memoise_2.0.1 beeswarm_0.4.0
## [185] AnnotationDbi_1.62.1 cluster_2.1.4
## [187] globals_0.16.2 GSEABase_1.62.0