Chapter 9 Multi-condition analysis
9.1 Motivation
Some of our most interesting scRNA-seq datasets consist of multiple samples collected across different experimental conditions. The idea is to partition cells into putative cell types or states as previously described (Chapters 6, 8), and then identify changes in gene expression or abundance between conditions for each cell type/state. This can yield some useful insights on the differences between the conditions, which is a nice change from the descriptive nature of most single-cell analyses. In fact, a multi-condition analysis represents one of the rare cases where we will actually do some formal hypothesis testing. So, this chapter is our chance to claw our way back towards some semblance of statistical rigor.
9.2 Differential expression
Differential expression (DE) is an obvious low-hanging fruit when it comes to detecting differences between conditions. Specifically, our goal is to test for DE between conditions within each cell type/state identified from the single-cell data. This resolves changes in expression to specific subpopulations, which is more informative than the corresponding bulk RNA-seq analysis32. To illustrate, let’s pull out some pancreas data generated from normal donors and patients with type II diabetes (Segerstolpe et al. 2016):
library(scRNAseq)
sce.seger <- SegerstolpePancreasData()
table(sce.seger$individual, sce.seger$disease)##
## normal type II diabetes mellitus
## H1 96 0
## H2 352 0
## H3 383 0
## H4 383 0
## H5 383 0
## H6 383 0
## T2D1 0 383
## T2D2 0 383
## T2D3 0 384
## T2D4 0 384Typically, we would cluster the cells (possibly after batch correction) and then assign some biological interpretation to each cluster. Happily enough, the authors provided cell type labels so we’ll use those directly instead of defining clusters ourselves.
##
## MHC class II cell PSC cell
## 5 54
## acinar cell alpha cell
## 185 886
## beta cell co-expression cell
## 270 39
## delta cell ductal cell
## 114 386
## endothelial cell epsilon cell
## 16 7
## gamma cell mast cell
## 197 7
## unclassified cell unclassified endocrine cell
## 2 41We compute “pseudo-bulk” expression profiles (Tung et al. 2017) by summing counts together for all cells with the same combination of cell type and sample.
As their name suggests, these pseudo-bulk profiles are intended to mimic bulk RNA-seq data so that they can be analyzed with existing DE workflows, e.g., edgeR, voom().
We use the sum of counts for several reasons:
- Larger counts are more amenable to analysis workflows designed for bulk RNA-seq data. Normalization is more straightforward and certain statistical approximations are more accurate, e.g., the saddlepoint approximation for quasi-likelihood methods or normality for linear models.
- Collapsing cells into samples reflects the fact that our biological replication occurs at the sample level (Lun and Marioni 2017). Each sample is represented no more than once for each condition, avoiding problems from unmodelled correlations between samples. Supplying the per-cell counts directly to a bulk RNA-seq workflow would imply that each cell is an independent biological replicate, which is not true from an experimental perspective. (A mixed effects model can handle this variance structure but involves extra complexity, typically for little benefit - see Crowell et al. (2020).)
- Variance across cells within each sample is ignored, provided it does not affect the variance across (replicate) samples. This is generally desirable for multi-condition analyses where the primary goal is to find consistent differences between conditions. Consider a gene with a strong and consistent change in expression upon treatment but only in a subset of the cells. We would happily consider this as a DE gene between conditions, despite the fact that the variance across cells in the treated samples would be inflated. By comparison, marker genes should be consistently up- or down-regulated in all cells between subpopulations, hence the use of the per-cell expression values in Chapter 7.
library(scrapper)
pseudo.bulk.seger <- aggregateAcrossCells.se(sce.seger, colData(sce.seger)[,c("individual","cell type")])
pseudo.bulk.seger## class: SummarizedExperiment
## dim: 26179 119
## metadata(1): aggregated
## assays(2): sums detected
## rownames(26179): SGIP1 AZIN2 ... BIVM-ERCC5 eGFP
## rowData names(2): refseq symbol
## colnames: NULL
## colData names(12): factor.individual factor.cell type ... submitted
## single cell quality cell type## DataFrame with 119 rows and 3 columns
## factor.individual factor.cell type counts
## <character> <character> <integer>
## 1 H1 NA 23
## 2 H1 MHC class II cell 1
## 3 H1 PSC cell 1
## 4 H1 acinar cell 4
## 5 H1 alpha cell 28
## ... ... ... ...
## 115 T2D4 delta cell 35
## 116 T2D4 ductal cell 47
## 117 T2D4 epsilon cell 1
## 118 T2D4 gamma cell 34
## 119 T2D4 mast cell 1## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## SGIP1 0 0 389 1 0 0 0 0 0 0
## AZIN2 0 0 0 0 125 0 0 0 0 0
## CLIC4 1 0 52 503 856 401 332 4 130 368
## AGBL4 0 0 0 0 399 0 0 0 0 0
## NECAP2 0 91 0 284 1835 491 28 3 354 0
## SLC45A1 0 0 0 0 383 36 0 0 0 0
## TGFBR3 0 0 0 341 1 232 262 541 0 0
## DBT 0 0 0 113 274 12 164 333 56 0
## RFWD2 0 0 0 162 340 91 9 6 0 0
## C1orf21 0 0 138 27 2125 0 149 123 207 0Once we’ve generated the pseudo-bulk count matrix, we test for differences between conditions - in this case, disease status - within each cell type.
Any DE analysis method that works with bulk RNA-seq data can be used, provided that we have replicates within each condition.
Here, we’ll be using voom() from the limma package (Law et al. 2014).
We won’t go into too much detail as there is plentiful documentation elsewhere, e.g., in the limmaUsersGuide(),
though there are a few pieces of advice that are specific to pseudo-bulk samples:
- Consider removing unreliable pseudo-bulk profiles with very few cells. The exact threshold depends on the dataset, the rarity of the cell type, the variance of the assay technology (e.g., UMIs versus reads), and whether the DE analysis supports downweighting of low-quality profiles. A good rule of thumb seems to be 10 cells (Crowell et al. 2020).
- Perform a separate analysis for each cell type instead of cramming all cell types into the same design matrix. This protects against differences in the mean-variance relationship across cell types. It also ensures that any odd behavior for one cell type does not affect results for the other cell types.
- Get used to higher variances and fewer DE genes compared to actual bulk RNA-seq data. The number of cells contributing to each pseudo-bulk profile is often orders of magnitude less than that used in bulk RNA-seq, so the latter will be a more precise assay of the population transcriptome.
We test for disease-associated DE genes in beta cells using voom() with additional weighting for sample quality.
Perhaps unsurprisingly, INS is one of the top DE genes.
We could then repeat this analysis on all available cell types, though for brevity’s sake, we won’t show that here.
Check out other Bioconductor packages like muscat, which implement some convenient functions for iterative DE analyses within cell types.
pseudo.beta.seger <- pseudo.bulk.seger[,which(pseudo.bulk.seger$`factor.cell type` == "beta cell")]
# We can have a look at the number of cells contributing to each profile, in
# case we want to remove low-abundance profiles.
pseudo.beta.seger$counts## [1] 12 48 32 34 10 35 10 14 11 64library(edgeR)
y.beta.seger <- DGEList(assay(pseudo.beta.seger, "sums"), samples=as.data.frame(colData(pseudo.beta.seger)))
keep.beta.seger <- filterByExpr(y.beta.seger, group=y.beta.seger$samples$disease)
y.beta.seger <- y.beta.seger[keep.beta.seger,]
y.beta.seger <- normLibSizes(y.beta.seger)
design.beta.seger <- model.matrix(~disease, y.beta.seger$samples)
v.beta.seger <- voomWithQualityWeights(y.beta.seger, design.beta.seger)
fit.beta.seger <- lmFit(v.beta.seger)
fit.beta.seger <- eBayes(fit.beta.seger, robust=TRUE)
res.beta.seger <- topTable(fit.beta.seger, sort.by="p", n=Inf, coef=2)
head(res.beta.seger)## ID logFC AveExpr t P.Value adj.P.Val B
## 7287 INS -2.761129 16.655021 -7.500553 4.296911e-06 0.05020081 4.616178
## 7689 FXYD2 -3.519464 5.676002 -6.864834 1.072013e-05 0.05020081 2.602945
## 7688 FXYD2 -2.589501 7.284994 -6.795743 1.195542e-05 0.05020081 3.388534
## 8349 ARL6IP4 -1.726072 7.811593 -5.660748 7.428443e-05 0.11330881 1.868303
## 11413 HPN -1.797720 6.132052 -5.654938 7.501799e-05 0.11330881 1.664214
## 11187 TRAPPC5 -2.121947 7.037605 -5.645027 7.628700e-05 0.11330881 1.768896In effect, we treat our scRNA-seq data analysis as a kind of super-powered in silico FACS33. Here, experimental isolation of cell types based on a few surface markers is replaced by computational assignment of cell types based on their transcriptomic profiles. Our pseudo-bulk DE analysis is analogous to a FACS experiment followed by bulk RNA-seq on the isolated populations.
9.3 Differential abundance
We can also test for differences in cell type abundance between conditions, i.e., differential abundance (DA). As we all know, immunologists love to create FACS plots showing some change in the percentages between treatments, e.g., Figure 1A of Richard et al. (2018). Now, we can do the same with scRNA-seq data. For our pancreas dataset, we create a count matrix of the number of cells assigned to each cell type in each sample. (If we didn’t already have annotated cell types, we could instead consider using a tool like miloR, which performs a DA analysis without requiring explicit assignment of each cells to clusters.)
ab.count.seger <- countGroupsByBlock(colData(sce.seger)[,"cell type"], colData(sce.seger)$individual)
ab.count.seger <- unclass(ab.count.seger) # get rid of the weird table class.
ab.count.seger## block
## groups H1 H2 H3 H4 H5 H6 T2D1 T2D2 T2D3 T2D4
## MHC class II cell 1 0 0 0 0 0 0 2 1 1
## PSC cell 1 1 2 6 3 10 2 12 13 4
## acinar cell 4 20 80 3 2 3 8 28 24 13
## alpha cell 28 117 26 136 44 92 141 119 87 96
## beta cell 12 48 32 34 10 35 10 14 11 64
## co-expression cell 3 3 5 6 3 6 1 5 1 6
## delta cell 7 21 2 7 10 12 9 6 5 35
## ductal cell 4 19 67 8 23 14 3 76 125 47
## endothelial cell 1 1 0 1 2 8 1 1 1 0
## epsilon cell 0 1 1 0 0 3 1 0 0 1
## gamma cell 7 19 15 2 1 31 70 8 10 34
## mast cell 0 4 0 0 0 0 0 2 0 1
## unclassified cell 0 0 0 0 0 1 1 0 0 0
## unclassified endocrine cell 5 15 4 0 0 5 3 3 6 0We then apply standard DA pipelines to see which cell types are affected by disease. In particular, testing for DA is bread-and-butter stuff in the microbiome field, so we’d recommend checking out some of their best practices. Right now, though, this book is hard enough to compile without adding extra dependencies, so we’ll just re-use edgeR’s statistical machinery to test for differences in the cell abundance matrix (Robinson, McCarthy, and Smyth 2010)34.
y.ab.seger <- DGEList(ab.count.seger)
y.ab.seger$samples$disease <- sce.seger$disease[match(colnames(y.ab.seger), sce.seger$individual)]
keep.ab.seger <- filterByExpr(y.ab.seger, group=y.ab.seger$samples$disease)
y.ab.seger <- y.ab.seger[keep.ab.seger,]
# If we don't normalize, our results will be affected by composition bias. But
# if we use TMM normalization, that would assume that most cell types do not
# have any change in their abundance. Hard to tell which one's worse here.
design.ab.seger <- model.matrix(~disease, y.ab.seger$samples)
fit.ab.seger <- glmQLFit(y.ab.seger, design.ab.seger)
res.ab.seger <- glmQLFTest(fit.ab.seger, coef=2)
topTags(res.ab.seger)## Coefficient: diseasetype II diabetes mellitus
## logFC logCPM F PValue FDR
## ductal cell 0.82000345 17.40873 1.38453794 0.2451320 0.6327002
## beta cell -0.82063001 16.98834 1.08166503 0.3035359 0.6327002
## gamma cell 0.78555501 16.53608 1.02526643 0.3163501 0.6327002
## acinar cell -0.45764795 16.45058 0.33776409 0.5638419 0.8457628
## delta cell -0.29945991 15.88049 0.13536851 0.7145470 0.8574564
## alpha cell -0.02224822 18.64197 0.00231903 0.9617915 0.9617915It’s worth noting that DA and DE are two sides of the same coin as they are both inferred from the per-cell expression profiles. Consider a scRNA-seq experiment involving two biological conditions with several shared cell types. We focus on a cell type \(X\) that is present in both conditions but contains some DE genes between conditions. This leads to two possible outcomes:
- The DE between conditions is strong enough to split \(X\) into two separate clusters (say, \(X_1\) and \(X_2\)) in expression space. This manifests as DA where \(X_1\) is enriched in one condition and \(X_2\) is enriched in the other condition.
- The DE between conditions is not sufficient to split \(X\) into two separate clusters, e.g., because our batch correction algorithm identifies them as corresponding cell types and merges them together. Thus, the differences between conditions manifest as DE within the single cluster corresponding to \(X\).
It is difficult to predict whether a difference between conditions will manifest as DE or DA. For example, we might see DE for coarser clusters but DA for finer clusters. We’d recommend performing both DE and DA analyses to ensure that we can detect either possibility.
9.4 More thoughts on statistical rigor
In the sections above, we were fortunate enough to use pre-existing cell type annotations from Segerstolpe et al. (2016). For actual multi-condition analyses, we would first have to assign our own biological identities to our cell subpopulations. Typically, this involves some kind of correction to merge shared cell types across samples (Chapter 8), clustering on the corrected data (Chapters 6), and finally examination of each cluster’s marker genes (Chapter 7). Each cluster is used as a proxy for a cell state/type identity that is common across samples, serving as the basis for cell type-specific differential expression or abundance across conditions.
The most obvious concern here is that the hypothesis testing is performed on the same data used to define the subpopulations. This represents a form of data snooping that complicates the interpretation of the DE/DA \(p\)-values. For example, each cluster will consist of cells with similar expression profiles, which may (i) artificially deflate the variance across replicates but also (ii) understate the differences between conditions in a DE analysis. In most cases, though, it’s probably fine as the process of defining the cell types/states (clustering or otherwise) is blind to the condition label of each cell. Any arbitrary placement of cell type/state boundaries in high-dimensional expression space should be more-or-less independent of any differences between conditions. Of course, we can easily imagine exceptions to this rule but these seem slightly pathological35.
In our opinion, the use of a common clustering is the real Achilles heel of this strategy in terms of statistical rigor. We fail to capture the uncertainty in the clustering and its biological interpretation, which reduces confidence in the reproducibility of the results. Say we discover significant DE/DA for a cell type in our dataset. If an independent party were to repeat our experiment and analysis, would they be able to reach the same conclusion? More specifically, would they be able to partition an equivalent cluster and assign the same cell type identity? Weakly separated cell subtypes might not manifest as separate clusters in a new dataset, or the ranking of markers might change in a manner that causes the analyst to assign a different biological identity. We wouldn’t know - we can’t evaluate the reproducibility of our cell type annotations because we only did the clustering and interpretation once. The same criticism applies to the interpretation of any common manifold, even if no explicit clustering is performed.
To perform multi-condition analyses “more correctly”, we need to process each sample independently to capture the variation in interpretation. Consider a dataset that has multiple replicate samples for each of multiple conditions. Our analysis strategy would look something like this:
- Analyze each sample independently, from quality control to identification of cell types/states from the clusters. Specifically, biological meaning should be assigned to cell subpopulations without any information from other samples. Indeed, if we were being very careful, we would blind and randomize samples across multiple analysts so that variances in human bias are also modelled during manual annotation36. Alternatively, we could use automated cell type annotation tools like SingleR; these do not require any clustering and can be applied to each sample independently, but assume that our cell types of interest exist in the reference annotation.
- Match corresponding cell types or states across samples. This is pretty straightforward if the per-sample cell type/state assignments use a controlled vocabulary, e.g, from the Cell Ontology, where the biology of interest is explicitly defined for all samples (and analysts, if more than one person is involved). We avoid the use of batch correction to merge cells across samples, which means that we aren’t affected by the assumptions and errors of the correction algorithm. It’s usually at this point that people often start complaining about their favorite cell types/states not showing up consistently across samples. But frankly, if the biology is real, it had better be reproducible across your replicates37, otherwise it’s just wishful thinking.
- Create a pseudo-bulk or cell abundance count matrix based on the annotated cell types/states from all samples. Any variability in the per-sample analysis will manifest as greater variance across replicates in these count matrices. For example, if a cell subtype is weakly defined, we may not be able to identify it consistently across replicates, increasing the variance in the cell type abundances. Similarly, if a subtype is poorly separated from its relatives, its cluster may occasionally include cells from neighboring subtypes, increasing the variance of the pseudo-bulk profiles. The increased variance is important as it properly reflects our uncertainty about the existence of the cell subtype itself.
In practice, this kind of analysis is pretty exhausting, especially for larger studies. We’ve only seen this approach used on a handful of occasions over the years because it’s just too inconvenient. Besides, the incentives for reproducibility don’t exist in most scientific environments38. We typically settle on a compromise between convenience and rigor, where we still use a common clustering from corrected PCs but invest the extra time and resources into independent validation experiments (see also suggestions in Section 7.7). As long as our conclusions can be validated, we can say that our preceding analyses were “exploratory” and give ourselves a pass for any statistical impropriety.
Session information
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB 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: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] edgeR_4.11.4 limma_3.69.2
## [3] scrapper_1.7.3 scRNAseq_2.27.0
## [5] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
## [7] Biobase_2.73.2 GenomicRanges_1.65.1
## [9] Seqinfo_1.3.0 IRanges_2.47.2
## [11] S4Vectors_0.51.6 BiocGenerics_0.59.10
## [13] generics_0.1.4 MatrixGenerics_1.25.0
## [15] matrixStats_1.5.0 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.3.0 bitops_1.1-0 httr2_1.3.0
## [4] rlang_1.3.0 magrittr_2.0.5 otel_0.2.0
## [7] gypsum_1.9.0 compiler_4.6.1 RSQLite_3.53.3
## [10] GenomicFeatures_1.65.0 png_0.1-9 vctrs_0.7.3
## [13] ProtGenerics_1.45.0 pkgconfig_2.0.3 crayon_1.5.3
## [16] fastmap_1.2.0 dbplyr_2.6.0 XVector_0.53.0
## [19] Rsamtools_2.29.0 rmarkdown_2.31 UCSC.utils_1.9.0
## [22] bit_4.6.0 xfun_0.60 beachmat_2.29.0
## [25] cachem_1.1.0 cigarillo_1.3.1 GenomeInfoDb_1.49.1
## [28] jsonlite_2.0.0 blob_1.3.0 rhdf5filters_1.25.2
## [31] DelayedArray_0.39.3 Rhdf5lib_2.1.0 BiocParallel_1.47.0
## [34] parallel_4.6.1 R6_2.6.1 bslib_0.11.0
## [37] rtracklayer_1.73.0 jquerylib_0.1.4 Rcpp_1.1.2
## [40] bookdown_0.47 knitr_1.51 BiocBaseUtils_1.15.1
## [43] Matrix_1.7-6 tidyselect_1.2.1 abind_1.4-8
## [46] yaml_2.3.12 codetools_0.2-20 curl_7.1.0
## [49] lattice_0.22-9 alabaster.sce_1.13.0 tibble_3.3.1
## [52] KEGGREST_1.53.6 evaluate_1.0.5 BiocFileCache_3.3.0
## [55] alabaster.schemas_1.13.0 ExperimentHub_3.3.1 Biostrings_2.81.6
## [58] pillar_1.11.1 BiocManager_1.30.27 filelock_1.0.3
## [61] RCurl_1.98-1.19 BiocVersion_3.24.0 ensembldb_2.37.3
## [64] alabaster.base_1.13.1 glue_1.8.1 alabaster.ranges_1.13.0
## [67] alabaster.matrix_1.13.0 lazyeval_0.2.3 tools_4.6.1
## [70] AnnotationHub_4.3.2 BiocIO_1.23.3 BiocNeighbors_2.7.2
## [73] locfit_1.5-9.12 GenomicAlignments_1.49.1 XML_3.99-0.23
## [76] rhdf5_2.57.3 grid_4.6.1 AnnotationDbi_1.75.2
## [79] HDF5Array_1.41.0 restfulr_0.0.17 cli_3.6.6
## [82] rappdirs_0.3.4 S4Arrays_1.13.0 dplyr_1.2.1
## [85] AnnotationFilter_1.37.0 alabaster.se_1.13.0 sass_0.4.10
## [88] digest_0.6.39 SparseArray_1.13.2 rjson_0.2.23
## [91] memoise_2.0.1 htmltools_0.5.9 lifecycle_1.0.5
## [94] h5mread_1.5.0 httr_1.4.8 statmod_1.5.2
## [97] bit64_4.8.2References
Crowell, H. L., C. Soneson, P.-L. Germain, D. Calini, L. Collin, C. Raposo, D. Malhotra, and M. D. Robinson. 2020. “muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data.” Nat. Commun. 11: 6077.
Law, C. W., Y. Chen, W. Shi, and G. K. Smyth. 2014. “voom: Precision weights unlock linear model analysis tools for RNA-seq read counts.” Genome Biol. 15 (2): R29.
Lun, A. T. L., and J. C. Marioni. 2017. “Overcoming confounding plate effects in differential expression analyses of single-cell RNA-seq data.” Biostatistics 18 (3): 451–64.
Richard, A. C., A. T. L. Lun, W. W. Y. Lau, B. Gottgens, J. C. Marioni, and G. M. Griffiths. 2018. “T cell cytolytic capacity is independent of initial stimulation strength.” Nat. Immunol. 19 (8): 849–58.
Robinson, M. D., D. J. McCarthy, and G. K. Smyth. 2010. “edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.” Bioinformatics 26 (1): 139–40.
Segerstolpe, A., A. Palasantza, P. Eliasson, E. M. Andersson, A. C. Andreasson, X. Sun, S. Picelli, et al. 2016. “Single-cell transcriptome profiling of human pancreatic islets in health and type 2 diabetes.” Cell Metab. 24 (4): 593–607.
Tung, P. Y., J. D. Blischak, C. J. Hsiao, D. A. Knowles, J. E. Burnett, J. K. Pritchard, and Y. Gilad. 2017. “Batch effects and the effective design of single-cell gene expression studies.” Sci. Rep. 7 (January): 39921.
Though also more expensive.↩︎
Fluorescence-activated cell sorting, duh.↩︎
With a hammer like edgeR, everything kind of looks like a nail.↩︎
For example, if one condition has many more cells than the other, it would dictate the definition of each subpopulation, which could bias the DA analysis. Or, at very fine clusterings, the distributional assumptions of some DE tools are violated due to underdispersion.↩︎
Though this is so exceptionally laborious, it’s probably not worth doing anything less serious than a clinical trial.↩︎
If you cant’t even find your cell type in your own replicates, what chance is there of reproducing it in an independent study?↩︎
Why should we do more work to introduce more variance and reduce the number of significant hits? This is antithetical to the raison d’être of single-cell genomics, which is to create publishable results.↩︎