Using crumblr in practice
Introduction
Changes in cell type composition play an important role in health and
disease. Recent advances in single cell technology have enabled
measurement of cell type composition at increasing cell lineage
resolution across large cohorts of individuals. Yet this raises new
challenges for statistical analysis of these compositional data to
identify changes associated with a phenotype. We introduce
crumblr, a scalable statistical method for analyzing count
ratio data using precision-weighted linear models incorporating random
effects for complex study designs. Uniquely, crumblr
performs tests of association at multiple levels of the cell lineage
hierarchy using multivariate regression to increase power over tests of
a single cell component. In simulations, crumblr increases
power compared to existing methods, while controlling the false positive
rate.
The crumblr package integrates with the variancePartition
and dreamlet
packages in the Bioconductor ecosystem.
Here we consider counts for 8 cell types from quantified using single cell RNA-seq data from unstimulated and interferon β stimulated PBMCs from 8 subjects (Kang, et al., 2018).
The functions here incorporate the precision weights:
variancePartition::fitExtractVarPartModel()variancePartition::dream()limma::lmFit()
Installation
To install this package, start R and enter:
Analysis workflow
Process data
Here we evaluate whether the observed cell proportions change in response to interferon β. Given the results here, we cannot reject the null hypothesis that interferon β does not affect the cell type proportions.
library(crumblr)
# Load cell counts, clustering and metadata
# from Kang, et al. (2018) https://doi.org/10.1038/nbt.4042
data(IFNCellCounts)
# Apply crumblr transformation
# cobj is an EList object compatable with limma workflow
# cobj$E stores transformed values
# cobj$weights stores precision weights
# corresponding to the regularized inverse variance
cobj <- crumblr(df_cellCounts)Variance partitioning
Decomposing the variance illustrates that more variation is explained by subject than stimulation status.
library(variancePartition)
# Partition variance into components for Subject (i.e. ind)
# and stimulation status, and residual variation
form <- ~ (1 | ind) + (1 | StimStatus)
vp <- fitExtractVarPartModel(cobj, form, info)
# Plot variance fractions
fig.vp <- plotPercentBars(vp)
fig.vp
PCA
Performing PCA on the transformed cell counts indicates that the
samples cluster based on subject rather than stimulation status. Here,
we use standardize() so that each observation has
approximately equal variance (i.e. homoscedasticity) by dividing the CLR
transformed frequencies by their estimated sampling standard deviation.
Transforming the data to be approximately homoscedastic has been shown
to improve performance of PCA.
library(ggplot2)
# Perform PCA
# use crumblr::standardize() to get values with
# approximately equal sampling variance,
# which is a key property for downstream PCA and clustering analysis.
pca <- prcomp(t(standardize(cobj)))
# merge with metadata
df_pca <- merge(pca$x, info, by = "row.names")
# Plot PCA
# color by Subject
# shape by Stimulated vs unstimulated
ggplot(df_pca, aes(PC1, PC2, color = as.character(ind), shape = StimStatus)) +
geom_point(size = 3) +
theme_classic() +
theme(aspect.ratio = 1) +
scale_color_discrete(name = "Subject") +
xlab("PC1") +
ylab("PC2")
Differential testing
# Use variancePartition workflow to analyze each cell type
# Perform regression on each cell type separately
# then use eBayes to shrink residual variance
# Also compatible with limma::lmFit()
fit <- dream(cobj, ~ StimStatus + ind, info)
fit <- eBayes(fit)
# Extract results for each cell type
topTable(fit, coef = "StimStatusstim", number = Inf)## logFC AveExpr t P.Value adj.P.Val B
## CD8 T cells -0.25085170 0.0857175 -4.0787416 0.002436375 0.01949100 -1.279815
## Dendritic cells 0.37386979 -2.1849234 3.1619195 0.010692544 0.02738587 -2.638507
## CD14+ Monocytes -0.10525402 1.2698117 -3.1226341 0.011413912 0.02738587 -2.709377
## B cells -0.10478652 0.5516882 -3.0134349 0.013692935 0.02738587 -2.940542
## CD4 T cells -0.07840101 2.0201947 -2.2318104 0.050869691 0.08139151 -4.128069
## FCGR3A+ Monocytes 0.07425165 -0.2567492 1.6647681 0.128337022 0.17111603 -4.935304
## NK cells 0.10270672 0.3797777 1.5181860 0.161321761 0.18436773 -5.247806
## Megakaryocytes 0.01377768 -1.8655172 0.1555131 0.879651456 0.87965146 -6.198336Multivariate testing along a tree
We can gain power by jointly testing multiple cell types using a
multivariate statistical model, instead of testing one cell type at a
time. Here we construct a hierarchical clustering between cell types
based on gene expression from pseudobulk, and perform a multivariate
test for each internal node of the tree based on its leaf nodes. The
results for the leaves are the same as from topTable()
above. At each internal node treeTest() performs a fixed
effects meta-analysis of the coefficients of the leaves while modeling
the covariance between coefficient estimates. In the backend, this is
implemented using variancePartition::mvTest() and remaCor
package.
Here the hierarchical clustering, hcl, is precomputed
from pseudobulk gene expression using
buildClusterTreeFromPB().
# Perform multivariate test across the hierarchy
res <- treeTest(fit, cobj, hcl, coef = "StimStatusstim")
# Plot hierarchy and testing results
plotTreeTest(res)
Hierarchical clustering
The hierarchical clustering used by treeTest() can be
computed a number of ways, depending on the available data and
biological question. For example, see Article for details about how hierarchical
clustering was run in this dataset.
In general, hierarchical clustering can be computed from
- pseudobulked single cell gene expression
- cell type frequencies
# correlation matrix between all cell types
C <- cor(t(standardize(cobj)))
# convert to distance
dm <- as.dist(1 - abs(C))
# eval hierarchical clustering
hcl <- hclust(dm)- Newick formated tree computed from external data
# Make sure packages are installed
# BiocManager::install(c("ctc", "ape", "phylogram"))
library(ape)
library(ctc)
library(phylogram)
library(tidyverse)
# Write tree to file,
# edit manually
# then read back into R
#
# Specify tree as text in Newick format
txt = "((CD14+ Monocytes,(B cells,(Dendritic cells,Megakaryocytes))),(CD8 T cells,(NK cells,(CD4 T cells,FCGR3A+ Monocytes))));"
# read from text
hcl_from_txt <- read.tree(text = txt) %>%
as.dendrogram.phylo %>%
as.hclust
# Alternatively,
# write existing tree to file
# and edit manaully
write(hc2Newick(hcl),file='hcl.nwk')
hcl_from_txt2 <- read.tree(file = 'hcl.nwk') %>%
as.dendrogram.phylo %>%
as.hclustConsiderations
Computational scaling
The crumblr() workflow is very fast, especially compared
to more demanding differential expression analysies. Applying
crumblr() requires <1 sec, even for very large datsets.
Differential testing with dream() takes < 10 seconds for
fixed effects models and < 30 seconds for mixed models with typical
sample sizes and number of cell types. Running treeTest()
can be a little more demanding, but should finish in < 30 seconds
with 20 cell types.
Complex study designs
The crumblr() workflow can handle complex study designs
with repeated measures or multiple random effects. dream()
uses lme4::lmer() in the backend to fit weighted linear
mixed models. Considerations about defining the regression model are
described in documentation to variancePartition
or this book
by the author of lme4.
Tuning parameters
crumblr() uses two tuning parameters accessable to the
user. These are fixed to default values in all simulations and data
analysis. The user is strongly recommended to keep these dfault
values.
In order to deal with zero counts,
crumblr()uses a pseudocount (default: 0.5) added to all observed counts.For real data, the asymptotic variance formula can give weights that vary substantially across samples and give very high weights for a subset of samples. In order to address this, we regularize the weights to reduce the variation in the weights to have a maximum ratio (default of 5) between the maximum and specified quantile value (default of 5%).
Session Info
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04 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 LC_TIME=en_US.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] variancePartition_1.43.1 BiocParallel_1.47.0 limma_3.69.2
## [4] lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0
## [7] dplyr_1.2.1 purrr_1.2.2 readr_2.2.0
## [10] tidyr_1.3.2 tibble_3.3.1 tidyverse_2.0.0
## [13] glue_1.8.1 crumblr_1.5.3 ggplot2_4.0.3
## [16] BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 sys_3.4.3 jsonlite_2.0.0
## [4] magrittr_2.0.5 farver_2.1.2 nloptr_2.2.1
## [7] rmarkdown_2.31 fs_2.1.0 vctrs_0.7.3
## [10] minqa_1.2.8 ggtree_4.3.0 htmltools_0.5.9
## [13] S4Arrays_1.13.0 broom_1.0.13 SparseArray_1.13.2
## [16] gridGraphics_0.5-1 sass_0.4.10 KernSmooth_2.23-26
## [19] bslib_0.11.0 htmlwidgets_1.6.4 pbkrtest_0.5.5
## [22] plyr_1.8.9 cachem_1.1.0 buildtools_1.0.0
## [25] lifecycle_1.0.5 iterators_1.0.14 pkgconfig_2.0.3
## [28] Matrix_1.7-5 R6_2.6.1 fastmap_1.2.0
## [31] rbibutils_2.4.1 MatrixGenerics_1.25.0 digest_0.6.39
## [34] numDeriv_2016.8-1.1 aplot_0.3.1 patchwork_1.3.2
## [37] S4Vectors_0.51.5 GenomicRanges_1.65.1 labeling_0.4.3
## [40] timechange_0.4.0 abind_1.4-8 compiler_4.6.1
## [43] fontquiver_0.2.1 aod_1.3.3 withr_3.0.3
## [46] S7_0.2.2 backports_1.5.1 viridis_0.6.5
## [49] gplots_3.3.0 MASS_7.3-65 rappdirs_0.3.4
## [52] DelayedArray_0.39.3 corpcor_1.6.10 gtools_3.9.5
## [55] caTools_1.18.3 tools_4.6.1 otel_0.2.0
## [58] ape_5.8-1 remaCor_0.0.20 nlme_3.1-169
## [61] grid_4.6.1 reshape2_1.4.5 generics_0.1.4
## [64] gtable_0.3.6 tzdb_0.5.0 hms_1.1.4
## [67] XVector_0.53.0 BiocGenerics_0.59.10 pillar_1.11.1
## [70] yulab.utils_0.2.4 splines_4.6.1 treeio_1.37.0
## [73] lattice_0.22-9 dirmult_0.1.3-5 tidyselect_1.2.1
## [76] fontLiberation_0.1.0 SingleCellExperiment_1.35.1 maketools_1.3.2
## [79] knitr_1.51 fontBitstreamVera_0.1.1 reformulas_0.4.4
## [82] gridExtra_2.3.1 IRanges_2.47.2 Seqinfo_1.3.0
## [85] SummarizedExperiment_1.43.0 RhpcBLASctl_0.23-42 stats4_4.6.1
## [88] xfun_0.60 Biobase_2.73.1 statmod_1.5.2
## [91] matrixStats_1.5.0 stringi_1.8.7 lazyeval_0.2.3
## [94] ggfun_0.2.1 yaml_2.3.12 boot_1.3-32
## [97] evaluate_1.0.5 codetools_0.2-20 gdtools_0.5.1
## [100] BiocManager_1.30.27 ggplotify_0.1.3 cli_3.6.6
## [103] RcppParallel_5.1.11-2 systemfonts_1.3.2 Rdpack_2.6.6
## [106] jquerylib_0.1.4 Rcpp_1.1.2 zigg_0.0.2
## [109] EnvStats_3.1.0 Rfast_2.1.5.2 bitops_1.0-9
## [112] lme4_2.0-1 viridisLite_0.4.3 mvtnorm_1.4-1
## [115] tidytree_0.4.8 ggiraph_0.9.6 lmerTest_3.2-1
## [118] scales_1.4.0 fANCOVA_0.6-1 rlang_1.3.0