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.
# Remove cells with no assigned cell type.
sce.seger <- sce.seger[,!is.na(sce.seger$`cell type`)]
table(sce.seger$`cell type`)##
## 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: RangedSummarizedExperiment
## dim: 26179 109
## 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 109 rows and 3 columns
## factor.individual factor.cell type counts
## <character> <character> <integer>
## 1 H1 MHC class II cell 1
## 2 H1 PSC cell 1
## 3 H1 acinar cell 4
## 4 H1 alpha cell 28
## 5 H1 beta cell 12
## ... ... ... ...
## 105 T2D4 delta cell 35
## 106 T2D4 ductal cell 47
## 107 T2D4 epsilon cell 1
## 108 T2D4 gamma cell 34
## 109 T2D4 mast cell 1## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## SGIP1 0 389 1 0 0 0 0 0 0 205
## AZIN2 0 0 0 125 0 0 0 0 0 19
## CLIC4 0 52 503 856 401 332 4 130 368 2
## AGBL4 0 0 0 399 0 0 0 0 0 0
## NECAP2 91 0 284 1835 491 28 3 354 0 0
## SLC45A1 0 0 0 383 36 0 0 0 0 1
## TGFBR3 0 0 341 1 232 262 541 0 0 1
## DBT 0 0 113 274 12 164 333 56 0 76
## RFWD2 0 0 162 340 91 9 6 0 0 32
## C1orf21 0 138 27 2125 0 149 123 207 0 100Once 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 fluorescence-activated cell sorting (FACS). 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
Speaking of FACS, we all know that immunologists love to create FACS plots showing some change in cell type percentages between treatments, e.g., Figure 1A of Richard et al. (2018). Now, we can do the same type of analysis with scRNA-seq data. Let’s use our pancreas dataset to create a count matrix of the number of cells assigned to each cell type in each sample. Our aim is to find significant differences in cell type abundance between conditions, i.e., differential abundance (DA).
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 0Testing for DA is bread-and-butter stuff in the microbiome field, so we’d recommend checking out some of their best practices. Alternatively, 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. 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)33.
y.ab.seger <- DGEList(ab.count.seger)
y.ab.seger$samples$disease <- sce.seger$disease[match(colnames(y.ab.seger), sce.seger$individual)]
# Omit the filtering, we have so few subpopulations that we don't want to
# remove any of them. We also don't call normLibSizes(), see explanation below.
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
## endothelial cell -1.9862141 13.98399 3.5324823 0.06277655 0.5687385
## unclassified endocrine cell -1.2247533 14.73732 2.0382351 0.15616987 0.5687385
## co-expression cell -1.1086693 14.70168 1.6843622 0.19701234 0.5687385
## ductal cell 0.8195745 17.40873 1.5525920 0.21535330 0.5687385
## beta cell -0.8209702 16.98834 1.3570687 0.24652134 0.5687385
## gamma cell 0.7887476 16.53608 1.2028361 0.27510809 0.5687385
## MHC class II cell 1.4139126 13.44115 1.1571497 0.28436926 0.5687385
## epsilon cell -1.0648313 13.55123 0.7851160 0.37747958 0.6605893
## acinar cell -0.4493706 16.45058 0.3720378 0.54313116 0.8448707
## PSC cell 0.3397519 15.00618 0.1822520 0.67026450 0.8697769The trickiest part of the DA analysis is the choice of normalization strategy.
We did not call normLibSizes() in the code above, which causes edgeR to default to normalization by the total count.
This assumes that any variation in the total number of cells between samples is not biologically interesting (e.g., differences in dissociation efficiency) and should be removed.
If this assumption is not valid, our normalization could introduce spurious differences between samples due to composition bias.
For example, if an experimental treatment induced a large expansion of a particular subpopulation in one condition, the total number of cells would be inflated in that condition.
Normalization by the total count would then decrease the abundance of all other subpopulations, even if nothing about those other subpopulations changed between conditions.
We might consider other methods like TMM normalization to explicitly remove composition biases,
though these have their own assumptions (i.e., that most subpopulations are not DA) that may not always be appropriate.
Unfortunately, there’s no universally good answer to this problem so some caution is required when interpreting the DA results.
It’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 pathological34.
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 annotation35. 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 replicates36, 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 environments37. 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.
9.5 Ambient contamination in DE analyses
Alright, back to some prosaic concerns. Ambient contamination is often present in scRNA-seq datasets generated by high-throughput protocols. During dissociation, cell lysis releases RNA molecules into the ambient solution. These molecules are captured and sequenced in each cell’s reaction chamber, contributing counts to genes that are not otherwise expressed in that cell. In most situations, ambient contamination is just a mild nuisance that can be ignored. It adds a constant baseline expression profile to all cells from the same sample, but our analysis focuses on the differences between cell types/states, so any addition more-or-less cancels out in the end.
For DE analyses of multi-condition datasets, ambient contamination is more problematic if the ambient expression profile varies across samples. After a major experimental perturbation, strong upregulation of a gene in one cell type can contaminate other cell types in the same sample. Subsequently, DE genes detected for a particular cell type may be driven by differences in the ambient profiles rather than any intrinsic change in gene regulation. To illustrate, let’s consider the Tal1-knockout (KO) chimera data from Pijuan-Sala et al. (2019). Tal1 is a transcription factor that has known roles in erythroid differentiation, and the aim of the experiment was to determine if blocking of the erythroid lineage in the early mouse embryo diverted cells to other developmental fates.
## class: SingleCellExperiment
## dim: 29453 56122
## metadata(0):
## assays(1): counts
## rownames(29453): ENSMUSG00000051951 ENSMUSG00000089699 ...
## ENSMUSG00000095742 tomato-td
## rowData names(2): ENSEMBL SYMBOL
## colnames(56122): cell_1 cell_2 ... cell_56121 cell_56122
## colData names(9): cell barcode ... pool sizeFactor
## reducedDimNames(1): pca.corrected
## mainExpName: NULL
## altExpNames(0):We perform a pseudo-bulk DE analysis between wild-type (WT) and KO cells labelled as “neural crest”. The strongest DE genes are the hemoglobins, which are downregulated in the KO cells. This is rather surprising as these cells do not belong to the erythroid lineage and should not express hemoglobins at all. The most sober explanation is that the WT samples contain more hemoglobin transcripts in the ambient solution due to leakage from erythrocytes or their precursors. (As an aside, it’s worth mentioning that the “replicates” in this study are more technical than biological. This results in some exaggeration of the significance of the effects, as evidenced by the very low \(p\)-values. For the purposes of this section, though, we’ll just pretend that these are proper replicates.)
pseudo.tal1 <- aggregateAcrossCells.se(
sce.tal1,
factor=colData(sce.tal1)[,c("sample", "celltype.mapped")]
)
# Add blocking factor (based on my second-hand knowledge of the study).
pseudo.tal1$block <- factor(pseudo.tal1$sample %% 2 == 0)
# Subset to our neural crest cells.
pseudo.neural <- pseudo.tal1[,pseudo.tal1$celltype.mapped == "Neural crest"]
# Standard voom analysis. 'tomato' is the factor of interest here, as
# tdTomato-marked cells are those in which Tal1 knockout was induced.
y.neural <- DGEList(
assay(pseudo.neural, "sums"),
samples=as.data.frame(colData(pseudo.neural)),
genes=as.data.frame(rowData(pseudo.neural))
)
keep.neural <- filterByExpr(y.neural, group=y.neural$samples$tomato)
y.neural <- y.neural[keep.neural,]
y.neural <- normLibSizes(y.neural)
# Using a log-fold change threshold of 0.5 here, to avoid prioritizing genes
# with very small log-fold changes when the variance is also very low.
design.neural <- model.matrix(~block + tomato, y.neural$samples)
v.neural <- voomWithQualityWeights(y.neural, design.neural)
fit.neural <- lmFit(v.neural)
fit.neural <- treat(fit.neural, robust=TRUE, lfc=0.5)
res.neural <- topTreat(fit.neural, sort.by="none", n=Inf, coef="tomatoTRUE")
head(res.neural[order(res.neural$P.Value),])## ENSEMBL SYMBOL logFC AveExpr t
## ENSMUSG00000052217 ENSMUSG00000052217 Hbb-bh1 -8.064177 6.028610 -31.77490
## ENSMUSG00000055609 ENSMUSG00000055609 Hba-x -7.691496 5.551962 -27.45492
## ENSMUSG00000096768 ENSMUSG00000096768 Erdr1 1.908131 7.399343 25.95790
## ENSMUSG00000086503 ENSMUSG00000086503 Xist -7.537245 5.211324 -24.55175
## ENSMUSG00000052187 ENSMUSG00000052187 Hbb-y -8.357119 5.030808 -23.26475
## ENSMUSG00000069919 ENSMUSG00000069919 Hba-a1 -8.539531 3.213202 -15.87795
## P.Value adj.P.Val
## ENSMUSG00000052217 4.523957e-212 4.818014e-208
## ENSMUSG00000055609 1.070635e-160 5.701134e-157
## ENSMUSG00000096768 2.113232e-144 7.501972e-141
## ENSMUSG00000086503 7.914819e-130 2.107321e-126
## ENSMUSG00000052187 4.004540e-117 8.529670e-114
## ENSMUSG00000069919 1.979477e-56 3.513572e-53To quantify the impact of ambient contamination on the interpretation of our DE results, we obtain an estimate of the ambient expression profile from the unfiltered count matrix for each sample. We consider all barcodes with total counts at or below 100 to represent empty droplets (Lun et al. 2019), and we sum the counts for each gene across these barcodes to obtain an expression vector representing the ambient profile for each sample.
# Extracting an ambient profile from the unfiltered count matrix,
# i.e., prior to any cell calling from Cellranger or emptyDrops.
ambient.tal1 <- list()
for (s in pseudo.neural$sample) {
raw.tal1 <- Tal1ChimeraData(type="raw", samples=s)[[1]]
mat <- counts(raw.tal1)
ambient.tal1[[s]] <- SummarizedExperiment(
list(sums=cbind(rowSums(mat[,colSums(mat) <= 100]))),
rowData=rowData(raw.tal1)
)
}
ambient.tal1 <- do.call(cbind, ambient.tal1)
ambient.tal1$tomato <- pseudo.neural$tomato
ambient.tal1$block <- pseudo.neural$blockOnce we have the ambient counts, we incorporate them into the DE analysis to mitigate the impact of contamination. The standard approach would be to use a “difference of differences” contrast, but this particular dataset is a bit too messy for that. Instead, we’ll just identify genes that are DE between the WT and KO ambient profiles and remove them from our pseudo-bulk DE list. The aim is to filter out potentially problematic genes for which DE might be caused by ambient contamination. This is a simple but effective strategy that gets rid of the hemoglobins from our top DE genes.
y.ambient <- DGEList(
assay(ambient.tal1),
genes=rowData(ambient.tal1),
samples=colData(ambient.tal1)
)
# Re-using the same filter that we used for the pseudo-bulk analysis,
# just to make sure we're using the same set of genes.
y.ambient <- y.ambient[keep.neural,]
y.ambient <- normLibSizes(y.ambient)
design.ambient <- model.matrix(~block + tomato, y.ambient$samples)
v.ambient <- voomWithQualityWeights(y.ambient, design.ambient)
fit.ambient <- lmFit(v.ambient)
fit.ambient <- treat(fit.ambient, robust=TRUE, lfc=0.5)
res.ambient <- topTreat(fit.ambient, sort.by="none", n=Inf, coef="tomatoTRUE")
# Decorate our pseudo-bulk result table with some extra ambient statistics.
res.neural.extra <- res.neural
res.neural.extra$ambient.logFC <- res.ambient$logFC
res.neural.extra$ambient.P.Value <- res.ambient$P.Value
# Only keeping genes that are either (i) not DE in the ambient analysis,
# or (ii) DE in the wrong direction compared to the pseudo-bulk analysis.
filtered.neural <- sign(res.ambient$logFC) != sign(res.neural$logFC) | res.ambient$adj.P.Val >= 0.1
res.neural.extra <- res.neural.extra[filtered.neural,]
head(res.neural.extra[order(res.neural.extra$P.Value),])## ENSEMBL SYMBOL logFC AveExpr
## ENSMUSG00000086429 ENSMUSG00000086429 Gt(ROSA)26Sor 1.4808518 5.681105
## ENSMUSG00000031604 ENSMUSG00000031604 Msmo1 1.4930488 5.399914
## ENSMUSG00000042814 ENSMUSG00000042814 Mcts2 1.1375719 6.419303
## ENSMUSG00000093930 ENSMUSG00000093930 Hmgcs1 1.2526567 5.687960
## ENSMUSG00000059743 ENSMUSG00000059743 Fdps 0.9816048 7.207519
## ENSMUSG00000058258 ENSMUSG00000058258 Idi1 1.1743859 5.376187
## t P.Value adj.P.Val ambient.logFC
## ENSMUSG00000086429 11.289199 1.089441e-29 1.289172e-26 0.20815585
## ENSMUSG00000031604 10.481173 7.019974e-26 6.230227e-23 0.45862857
## ENSMUSG00000042814 9.288706 9.328503e-21 5.844033e-18 0.50419317
## ENSMUSG00000093930 8.742397 1.312587e-18 7.766141e-16 0.71946376
## ENSMUSG00000059743 8.652106 2.891564e-18 1.620797e-15 0.12809431
## ENSMUSG00000058258 7.118483 5.807591e-13 2.290772e-10 -0.05481356
## ambient.P.Value
## ENSMUSG00000086429 0.999912246
## ENSMUSG00000031604 0.687997586
## ENSMUSG00000042814 0.477107058
## ENSMUSG00000093930 0.005795168
## ENSMUSG00000059743 0.999999996
## ENSMUSG00000058258 0.999999852Alternatively, we could attempt to remove the ambient contamination from the count matrix itself, e.g., using tools like cellbender or SoupX.
I don’t have much to say on this topic other than to point out that
(i) estimating the contaminating proportion in each single cell is hard and
(ii) removing it while preserving the mean-variance relationship is hard38.
The difficulty of the problem motivates some rather complicated models and I prefer a lighter intervention.
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] MouseGastrulationData_1.27.0 SpatialExperiment_1.23.0
## [3] edgeR_4.11.6 limma_3.69.4
## [5] scrapper_1.7.8 scRNAseq_2.27.0
## [7] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
## [9] Biobase_2.73.2 GenomicRanges_1.65.3
## [11] Seqinfo_1.3.2 IRanges_2.47.5
## [13] S4Vectors_0.51.9 BiocGenerics_0.59.12
## [15] generics_0.1.4 MatrixGenerics_1.25.0
## [17] 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 magick_2.9.1 dbplyr_2.6.0
## [19] XVector_0.53.0 Rsamtools_2.29.0 rmarkdown_2.31
## [22] UCSC.utils_1.9.0 purrr_1.2.2 bit_4.6.0
## [25] xfun_0.60 beachmat_2.29.2 cachem_1.1.0
## [28] cigarillo_1.3.1 GenomeInfoDb_1.49.1 jsonlite_2.0.0
## [31] blob_1.3.0 rhdf5filters_1.25.4 DelayedArray_0.39.6
## [34] Rhdf5lib_2.1.0 BiocParallel_1.47.0 parallel_4.6.1
## [37] R6_2.6.1 bslib_0.12.0 rtracklayer_1.73.0
## [40] jquerylib_0.1.4 Rcpp_1.1.2 bookdown_0.47
## [43] knitr_1.51 BiocBaseUtils_1.15.1 Matrix_1.7-6
## [46] tidyselect_1.2.1 abind_1.4-8 yaml_2.3.12
## [49] codetools_0.2-20 curl_8.0.0 lattice_0.23-1
## [52] alabaster.sce_1.13.0 tibble_3.3.1 withr_3.0.3
## [55] BumpyMatrix_1.21.0 KEGGREST_1.53.6 evaluate_1.0.5
## [58] BiocFileCache_3.3.0 alabaster.schemas_1.13.0 ExperimentHub_3.3.2
## [61] Biostrings_2.81.7 pillar_1.11.1 BiocManager_1.30.27
## [64] filelock_1.0.3 RCurl_1.98-1.20 BiocVersion_3.24.0
## [67] ensembldb_2.37.3 alabaster.base_1.13.2 glue_1.8.1
## [70] alabaster.ranges_1.13.0 alabaster.matrix_1.13.1 lazyeval_0.2.3
## [73] tools_4.6.1 AnnotationHub_4.3.2 BiocIO_1.23.3
## [76] BiocNeighbors_2.7.3 locfit_1.5-9.12 GenomicAlignments_1.49.1
## [79] XML_3.99-0.24 rhdf5_2.57.12 grid_4.6.1
## [82] AnnotationDbi_1.75.2 HDF5Array_1.41.3 restfulr_0.0.17
## [85] cli_3.6.6 rappdirs_0.3.4 S4Arrays_1.13.0
## [88] dplyr_1.2.1 AnnotationFilter_1.37.0 alabaster.se_1.13.0
## [91] sass_0.4.10 digest_0.6.39 SparseArray_1.13.2
## [94] rjson_0.2.23 memoise_2.0.1 htmltools_0.5.9
## [97] lifecycle_1.0.5 h5mread_1.5.2 httr_1.4.8
## [100] statmod_1.5.2 bit64_4.8.4References
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., S. Riesenfeld, T. Andrews, T. P. Dao, T. Gomes, participants in the 1st Human Cell Atlas Jamboree, and J. Marioni. 2019. “EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data.” Genome Biol. 20 (1): 63.
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.
Pijuan-Sala, B., J. A. Griffiths, C. Guibentif, T. W. Hiscock, W. Jawaid, F. J. Calero-Nieto, C. Mulas, et al. 2019. “A Single-Cell Molecular Map of Mouse Gastrulation and Early Organogenesis.” Nature 566 (7745): 490–95.
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.↩︎
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.↩︎
This reminds me of the old problem of input subtraction in ChIP-seq data, which I wrestled with as part of my PhD. In the end, we just decided to filter out regions with high input coverage rather than attempting subtraction. It was a lose-lose situation; either the input counts were small and subtraction had no effect, or the input counts were large and subtraction broke the expected mean-variance relationship in our count data.↩︎