Chapter 8 — Differential Abundance: Do Cell-Type Proportions Change?
Where this chapter sits. Companion to Lecture 08 and Tutorial 08. Prerequisites: Chapter 5 — Integration (you need an integrated reduction), Chapter 6 — DESeq2: Bulk Fundamentals + Pseudobulk DE (parallel statistical machinery).
8.1 The complementary question
Differential expression (Chapters 5–6) asks: for cells of the same type, which genes change between conditions? Differential abundance asks: do cell-type proportions themselves change between conditions?
These are independent questions. A perturbation can change cell-type composition, gene expression within types, both, or neither. A complete condition-level analysis reports both. The four logical combinations are worth keeping in mind as you design and interpret:
- DE positive, DA negative: every immune cell shifts toward an antiviral state, but proportions are unchanged — the biology lives entirely in the transcriptional response.
- DA positive, DE negative: a memory T-cell population expands in infected samples, but expression within each cell type is unchanged — the biology lives entirely in cellular composition.
- Both positive: populations expand and the surviving or expanding cells express differently — common in strong perturbations, tumors, and disease progression.
- Both negative: no detectable change in either — either the biology is subtle, the experiment is underpowered, or the comparison is inappropriate.
Reviewers increasingly expect both analyses in a single-cell manuscript. Running only DE while ignoring composition, or only DA while ignoring per-cell transcriptional state, leaves half the biology unreported.
Differential expression (DE) and differential abundance (DA) address orthogonal aspects of a perturbation’s effect on a single-cell dataset. DE asks: for cells of the same annotated type, which genes are expressed at different levels between conditions? — a question answered by pseudobulk DESeq2 on per-(donor × condition × cell-type) aggregated counts. DA asks: do the relative proportions of cell types shift between conditions? — a question answered by testing whether the fraction of cells belonging to each neighborhood or cluster changes. These analyses measure different biological signals: a drug might expand a regulatory T cell population (DA-positive) while simultaneously shutting down cytokine production in every T cell (DE-positive), or might transcriptionally activate monocytes (DE-positive) without changing how many monocytes there are (DA-negative). Both analyses are routinely expected in high-quality single-cell manuscripts, and their combined interpretation — which populations expanded and what changed within them — provides a far more complete picture of the biology than either alone.
8.2 Why cluster-level proportion tests fail
The naive approach — count cells per (cluster, sample), compute each cluster’s fraction of the total, and run a Fisher’s exact or chi-squared test per cluster — is tempting because it requires no new tools. It is also systematically unreliable for four distinct reasons.
Arbitrary cluster boundaries. Cluster resolution is a tuning parameter, not a biological truth. A cell type that spans two clusters at one resolution sits in one cluster at another; the proportion test result changes accordingly. Resolution-dependent inferences are analytically fragile, and reviewers can ask for a different resolution to challenge any reported result.
Within-cluster heterogeneity is invisible. A “T cell” cluster containing both an expanding activated subpopulation and a shrinking naive subpopulation will show no net proportion change at the cluster level, even though biologically important DA is happening at finer resolution. The aggregation hides the signal.
Pseudoreplication. Cluster-level Fisher tests treat every cell as an independent sample. With 5,000 T cells from eight donors, the effective N becomes huge and the variance collapses to near zero, yielding extremely small p-values for compositional differences that would not replicate in a new cohort. The correct unit is the donor, not the cell, for exactly the same reason as in pseudobulk DE1.
The compositional constraint. Cell-type fractions within a sample must sum to 1. They are therefore not independent: when one population expands, the others must mathematically shrink, even if their absolute counts are unchanged. Naive proportion tests treat the fractions as independent observations, which produces spurious negative results for non-changing populations that appear to shrink purely as a consequence of another population’s expansion. This is the compositional problem, well known from microbiome analysis and increasingly recognized in single-cell DA. Methods like scCODA2 and sccomp3 model this explicitly using Bayesian compositional frameworks.
A dedicated benchmark of six DA methods evaluated these failure modes across real and simulated datasets and is the reference for method selection4. The benchmark’s main finding is that neighborhood-based methods and compositional-aware methods outperform naive cluster proportion tests, especially when biological gradients are present.
When cluster-level proportion tests are acceptable. The four failures above are most severe for continuous biology (gradients, trajectories) and small cohorts. For well-separated, discrete cell types in a large, well-powered cohort, cluster-level proportion tests using propeller’s logit-transform approach produce results that agree well with neighborhood-based methods. If your clusters are biologically well-defined (e.g. separate B cells, T cells, and NK cells), cell numbers per donor are large, and you have > 8 donors per condition, a cluster-level screen with propeller as the first pass is defensible — with miloR for follow-up on any uncertain result. The benchmark4 confirms this: cluster-level methods perform acceptably when assumptions are met, but neighborhood methods are consistently more robust when any of the four assumptions is violated.
8.3 The compositional problem in depth
The compositional constraint deserves a more precise treatment because it underlies the design choices for every DA method. In any sample, if there are \(C\) cell types with proportions \(p_1, p_2, \ldots, p_C\), then \(\sum_c p_c = 1\) by definition. This means the proportion vector lives on a \(C-1\) dimensional simplex, not in \(C\)-dimensional Euclidean space. Standard statistical models (linear regression, t-tests, logistic regression) assume an unconstrained observation space and are technically misspecified when applied to proportions.
In practice, the consequences are:
- Spurious negative DA. If CD14 monocytes expand from 10% to 30% of a PBMC sample, the remaining 90%→70% total fraction must distribute across all other cell types. Even if B cells are absolutely unchanged in number, their fraction drops because the denominator grew. A naive proportion test reports B cells as depleted — a false positive.
- Sign ambiguity. Methods that test proportions relative to the total (rather than modeling absolute counts) cannot distinguish “population A grew” from “everything else shrank relative to A.”
- Reference composition sensitivity. Removing one cell type from the analysis changes the proportions of all remaining types and can reverse test conclusions.
The practical response is to either (a) use absolute counts where possible, (b) apply a compositional transformation (additive log-ratio, isometric log-ratio) before testing, or (c) use a method that models the simplex geometry explicitly. Neighborhood-based methods like miloR sidestep the constraint by testing local abundance against a negative-binomial model of cell counts rather than proportions, which is one reason they perform robustly.
8.4 Available methods and when to use each
Table 1 summarizes the four main DA tools used in single-cell analysis, their statistical frameworks, and recommended use cases. No single method is best for all datasets; the benchmark4 provides guidance, but cross-validating two complementary approaches on your own data is the most defensible practice.
Understanding how each method addresses the compositional constraint and the pseudoreplication problem is the key to choosing among them. miloR sidesteps both issues by working on absolute cell counts per neighborhood (not proportions) and requiring a per-sample design analogous to pseudobulk DESeq2. scCODA addresses the compositional constraint directly through its Dirichlet-multinomial prior but works at cluster level, requiring you to trust your cluster annotations. propeller makes a partial correction via logit transformation but is otherwise a classical linear model — simple and fast but with the weakest theoretical treatment of compositionality. sccomp goes further by jointly modeling composition and within-type expression, capturing dependencies between the two that the other methods ignore entirely.
| Method | Statistical framework | Resolution | Handles composition? | Best suited to |
|---|---|---|---|---|
| miloR5 | NB GLM on kNN neighborhoods + SpatialFDR | Sub-cluster (continuous) | Indirectly (models counts) | Continuous trajectories; any dataset where cluster boundaries are uncertain |
| scCODA2 | Bayesian Dirichlet-multinomial compositional model | Cluster level | Yes — explicitly | Discrete cell types; small cohort with strong prior; need posterior DA probability |
| sccomp3 | Joint Bayesian model of composition and variability | Cluster level | Yes — explicitly | Datasets where both composition and within-cluster expression vary; large cohorts |
| propeller6 | Linear model on logit-transformed proportions | Cluster level | Partially (log-ratio transform) | Fast cluster-level screening; large well-powered cohorts; validation of miloR hits |
Design consideration: a practical two-step workflow is to run miloR for the neighborhood-resolution analysis and then cross-validate your top cell-type-level hits with propeller (which is faster and produces results directly on cluster labels). Concordance between the two methods strengthens the claim; discordance highlights cases where the neighborhood boundary definition matters and merits closer inspection.
8.5 The miloR approach — neighborhoods, not clusters
miloR5 replaces hard clusters with overlapping neighborhoods on the kNN graph. Each neighborhood is a representative cell plus its \(k\) nearest neighbors; roughly 1,000 neighborhoods are sampled to cover the full embedding. The key insight is that neighborhood membership is a continuous function of position in the reduced-dimension space — two cells sharing a neighborhood are, by definition, transcriptionally similar — so DA at this resolution is independent of any arbitrary clustering decision.
Table 2 summarizes the complete miloR workflow. Each step is explained in the sections that follow.
| Step | Action | Key function | Pitfall to watch |
|---|---|---|---|
| 1. Convert | Convert Seurat object to SingleCellExperiment, add integrated reduction |
as.SingleCellExperiment() + reducedDim() |
Harmony/PCA reduction must be explicitly transferred to the SCE |
| 2. Build graph | Construct kNN graph on integrated reduction | buildGraph(milo, k=30, d=20) |
Use the integrated (corrected) reduction, not raw PCA |
| 3. Sample nhoods | Sample representative neighborhoods | makeNhoods(milo, prop=0.1, refined=TRUE) |
refined=TRUE improves neighborhood coverage; don’t omit it |
| 4. Check sizes | Inspect neighborhood size distribution | plotNhoodSizeHist(milo) |
Target median 50–200; increase k if too small |
| 5. Count cells | Count cells per (sample × neighborhood) | countCells(milo, meta.data=..., sample=...) |
Sample column must match colData exactly |
| 6. Build design | Create per-sample metadata data frame | manual data.frame(...) |
One row per sample (donor × condition), not per cell |
| 7. Test | Fit NB GLM + SpatialFDR per neighborhood | testNhoods(milo, design=~donor+stim) |
Check resultsNames for correct contrast |
| 8. Annotate | Assign modal cell-type label to each neighborhood | annotateNhoods(milo, da_results, coldata_col=...) |
mixed_cell_types column flags heterogeneous neighborhoods |
| 9. Visualize | Plot DA beeswarm and neighborhood graph | plotDAbeeswarm(), plotNhoodGraphDA() |
Neighborhood graph layout requires a 2D reduction stored in reducedDim |
8.5.1 The kNN graph and neighborhood sampling
The kNN graph is the same structure used for Leiden/Louvain clustering in Chapter 2, built on the integrated PC reduction. Each cell is connected to its \(k\) nearest neighbors by Euclidean distance in that space. A neighborhood is a cell (the “index cell”) plus all its \(k\) neighbors — a set of \(k+1\) cells that are collectively close in transcriptional space.
Because the graph is built on the integrated reduction (Harmony, Seurat integration, or scVI), the neighborhoods represent biological proximity with batch effects removed. This is the critical dependency on Chapter 5: if integration was poor or over-corrected, the neighborhoods reflect artifacts rather than biology, and the DA test will detect those artifacts.
The makeNhoods step samples a subset of cells as index cells (typically 10% of the total, controlled by prop). With refined = TRUE, the sampling uses a second-round refinement step that ensures better spatial coverage across the graph. For a 25,000-cell dataset with prop = 0.1, this produces roughly 2,500 neighborhoods — enough to cover the embedding densely without testing every cell individually.
makeNhoods(milo, prop = 0.1, k = 30, d = 20, refined = TRUE) controls the neighborhood sampling. Default: prop = 0.1, k = 30, d = 20, refined = TRUE. prop sets what fraction of cells are sampled as neighborhood index cells — 10% of 25,000 cells produces ~2,500 neighborhoods, which is sufficient for dense coverage while keeping computation tractable. k is the number of nearest neighbors in the kNN graph (same parameter as for clustering); increasing it produces larger neighborhoods and smoothes the DA signal, decreasing it produces smaller, more granular neighborhoods. d is the number of integrated PCs used for graph construction; this should match the d used in the buildGraph call. refined = TRUE runs a second-pass refinement that improves the spatial distribution of sampled neighborhoods — omitting it produces neighborhoods that are clustered in dense regions of the embedding and sparse elsewhere. Always inspect plotNhoodSizeHist(milo) after makeNhoods to verify that the median neighborhood size falls between 50 and 200 cells before proceeding to testNhoods.
8.5.2 The negative-binomial GLM
For each neighborhood, miloR fits a generalized linear model analogous to the pseudobulk DESeq2 model of Chapter 6:
\[\log(\mu_{sn}) = X_s^\top \beta_n + \log(N_s)\]
where \(\mu_{sn}\) is the expected number of cells from sample \(s\) in neighborhood \(n\), \(X_s\) is the design matrix row for sample \(s\), \(\beta_n\) is the per-neighborhood coefficient vector (the DA effect), and \(N_s\) is the total cell count in sample \(s\) (an offset for sequencing depth / sample size). The likelihood is negative binomial, absorbing over-dispersion in counts just as in DESeq2.
The output coefficient for the condition contrast is the log-fold-change in neighborhood abundance: a positive LFC means the neighborhood is more populated in the perturbed condition, a negative LFC means it is depleted.
8.5.3 SpatialFDR — FDR correction for overlapping tests
Standard BH correction assumes independent tests. Neighborhoods overlap by construction — two adjacent neighborhoods share most of their constituent cells — so their test statistics are correlated. Applying BH naively would under-correct the FDR.
SpatialFDR5 addresses this by weighting each neighborhood’s p-value by the inverse of its graph connectivity: a neighborhood with many overlapping neighbors is downweighted because its result is not truly independent. The corrected p-values are stored as SpatialFDR in the results; the standard significance threshold is SpatialFDR < 0.1 (10% FDR), slightly more lenient than the conventional 5% because the correction is conservative.
To understand why 10% rather than 5%: the spatial weighting in SpatialFDR is inherently conservative — it assigns each neighborhood a weight proportional to the inverse of how many neighborhoods it overlaps with, which tends to over-penalize tests in dense regions of the graph. Simulation studies in the original miloR paper found that SpatialFDR < 0.1 corresponds approximately to a true FDR of 5–7% in practice, while SpatialFDR < 0.05 is too stringent and misses many genuine hits. In practice, use SpatialFDR < 0.1 as the primary filter and treat neighborhoods between 0.05 and 0.1 as borderline results warranting visual inspection of the neighborhood graph.
Filtering on da_results[SpatialFDR < 0.1,] controls which neighborhoods are called significant after correction for overlapping tests. De-facto default threshold: 0.1 (10% FDR). This is more lenient than the conventional 5% BH threshold because SpatialFDR’s spatial weighting is inherently conservative and a threshold of 0.1 corresponds empirically to a true FDR of approximately 5–7% in benchmarked simulations. Tightening to SpatialFDR < 0.05 is appropriate when you want the most conservative call set (e.g. for reporting in a clinical manuscript where false positives are costly). Loosening to SpatialFDR < 0.2 is occasionally justified for exploratory analyses in small cohorts where power is limited, but should be clearly stated as exploratory.
testNhoods(milo, design = ~ stim, design.df = milo_design) specifies the model to fit per neighborhood. Default: required argument with no default. The design formula follows exactly the same rules as in DESeq2 (see §6.5): ~ condition for a simple two-group contrast, ~ donor + condition for a paired design that blocks on donor-level variation (the recommended upgrade when donors are balanced across conditions), ~ batch + condition to block on batch effects, and continuous covariates (e.g. ~ age + condition) for quantitative phenotypes. The design.df must have one row per sample (donor × condition combination), not one row per cell. Adding donor to the design substantially increases power in paired experiments: donor-level differences in baseline cell-type composition are large (often 2–5-fold between donors), and failing to account for them makes the condition variance estimate too wide, reducing sensitivity. The design.df must contain all variables referenced in the design formula as columns.
8.6 Neighborhood size — the only parameter you’ll tune
The defaults (k = 30, d = 20, prop = 0.1) work well for most 10–50k-cell datasets. The diagnostic is plotNhoodSizeHist(milo): the distribution of cells per neighborhood should be approximately unimodal with a median between 50 and 200 cells.
Below ~30 cells per neighborhood the test is underpowered — a neighborhood that spans only a handful of cells is dominated by stochastic variation and the NB GLM cannot distinguish biological signal from noise. Above ~500, you are spatially smoothing so broadly that neighborhoods span multiple biologically distinct cell states, and DA detected at that scale reflects a blunt average rather than the fine-grained shift you want to capture.
If the histogram peak is below 30, increase k (try 50 or 75) and rebuild the graph before re-running makeNhoods. If the peak is above 500, decrease k to 20 or 15. For very large datasets (> 100,000 cells), reducing prop to 0.05 or 0.02 prevents testing an unwieldy number of neighborhoods without sacrificing coverage.
Why neighborhood size matters for the GLM. The NB GLM that testNhoods fits per neighborhood estimates a dispersion parameter for the cell counts across samples. With fewer than ~20–30 cells per sample in a neighborhood, that cell count is dominated by Poisson sampling noise from a handful of cells, and the dispersion estimate becomes unreliable. A neighborhood might show zero cells from one donor in CTRL simply by chance, which the GLM incorrectly interprets as a condition effect. With 50–100 cells per sample, the cell counts average out the stochastic variation of individual cells and reflect genuine compositional biology. This is the DA analogy of the pseudobulk minimum-cells-per-group threshold in Chapter 6: both exist to ensure that the input to the statistical model represents a real biological signal rather than sampling noise.
Always run plotNhoodSizeHist before testNhoods. Checking the size distribution takes seconds and catches mis-parameterized runs that would otherwise produce unreliable or uninterpretable results.
8.7 The design and replication requirements
The design formula passed to testNhoods follows exactly the same logic as the DESeq2 design in Chapter 6. Common designs are:
~ stim— two-condition contrast, no blocking covariates~ donor + stim— paired design; preferred when donors are balanced across conditions; blocks donor-level variation and increases power~ batch + stim— block on batch effects when samples were processed across runs~ severity— continuous covariate for dose-response or disease-progression analyses
The replication requirements are the same as for pseudobulk DE: at minimum 3–4 biological replicates per condition for any meaningful power, and at least 2 per condition for the test to be mathematically possible. With 1 vs 1 (one sample per condition), there are no degrees of freedom, and testNhoods returns NA p-values for every neighborhood.
Don’t pseudo-replicate by splitting cells into “fake donors”. The workshop loads ifnb via muscData::Kang18_8vs8(), which carries the real per-cell donor IDs from Kang et al. (2017) in the ind column — eight lupus donors, each measured in both CTRL and STIM. Tutorial 08 uses these real donors automatically and produces valid inferences. A synthetic-donor fallback is retained in the code only for edge cases where an object arrives without donor metadata (e.g. an older SeuratData-sourced copy); in that fallback, pseudo-replication makes the test anti-conservative — its p-values have no biological meaning. For your own data with real donor IDs, the workflow gives valid inferences without modification.
8.8 Why neighborhoods recover what clusters miss
To see concretely why neighborhoods outperform cluster-level tests for continuous biology, consider a naive-to-effector T cell activation gradient. Hard clustering will draw a single boundary somewhere along the gradient and label one side “naive T” and the other “effector T.” A cluster-level proportion test reports the proportions of those two bins and nothing more.
If a perturbation shifts the position of cells along the gradient — pushing more cells toward the effector end — neighborhoods near the effector end of the embedding show up-DA (more cells from the perturbed condition) and neighborhoods near the naive end show down-DA. The cluster-level test, depending on where the boundary fell, might report the change correctly, or might miss it entirely if the arbitrary boundary happened to sit in a region that did not shift. miloR’s neighborhood graph will show a gradient of log-fold-changes transitioning from negative (depleted naive end) to positive (enriched effector end), a much more informative result.
For very discrete biology — B cells vs T cells, where the populations are separated by thousands of gene-expression distances — neighborhoods give essentially the same answer as cluster-level tests. So neighborhoods are at worst equivalent to cluster tests and often substantially better. A key diagnostic is the third pattern in the DA beeswarm (§8.10.2): a cell-type cluster with neighborhoods scattered both above and below zero is the signature of within-cluster heterogeneity, a case where cluster-level testing would have returned a null result.
8.9 The integration prerequisite
miloR builds neighborhoods on whichever low-dimensional reduction you provide. Using the unintegrated PCA is a common error. In a multi-sample experiment where batch effects are present, the unintegrated kNN graph will partly cluster by sample rather than by cell type: cells from sample A appear as neighbors of other cells from sample A even if those cells are not transcriptionally similar. Neighborhoods will be enriched for sample-of-origin rather than cell-state, and testNhoods will detect composition shifts that are entirely batch artifacts.
The miloR neighborhood graph defines “similar cells” as those close together in the low-dimensional reduction you provide. If that reduction retains batch effects — as raw PCA often does in multi-sample experiments — then two cells from the same donor will appear as nearest neighbors even if they belong to different cell types, simply because they share a donor-specific transcriptional background. testNhoods would then detect abundance differences between conditions that are actually driven by donor-to-donor composition differences, not the condition effect you intended to test. Building the graph on the integrated reduction (Harmony, Seurat CCA, or scVI latent space from Chapter 5) ensures that nearest-neighbor relationships reflect cell-type biology rather than sample identity. This is not optional — it is the fundamental prerequisite for the neighborhood-level DA test to be interpretable.
The fix is to build the graph on the integrated reduction — the Harmony embeddings, the Seurat-integration corrected PCA, or the scVI latent space from Chapter 5. In Seurat:
milo <- buildGraph(milo, k = 30, d = 20, reduced.dim = "HARMONY")If you used the Seurat-anchor path instead (RunPCA on the integrated assay, as in Tutorial 05), the integrated reduction is "PCA", so pass reduced.dim = "PCA". The key requirement is that the reduction used for graph construction is the batch-corrected one — with Harmony that is "HARMONY", not "PCA".
Two further consequences of the integration dependency: (a) over-integration (Chapter 5) will cause neighborhoods to be artificially mixed across biological groups, reducing DA sensitivity; (b) under-integration will leave batch structure in the graph and inflate DA false positives. Running the Chapter 5 integration diagnostics before DA is strongly recommended.
8.10 Reading the output
8.10.1 The neighborhood graph
plotNhoodGraphDA(milo, da_results, layout = "UMAP_HARMONY", alpha = 0.1) overlays DA results on a neighborhood-resolution graph (not the cell-level UMAP). Each dot is a neighborhood, sized by the number of cells it contains and colored by the neighborhood’s log-fold-change in abundance (red = enriched in the perturbed condition; blue = depleted). Neighborhoods with SpatialFDR > 0.1 are rendered faint or transparent.
Look for:
- Compact bright clusters — a coherent region of cell-type space that consistently shifts between conditions. This is a robust finding; the signal is spatially contiguous.
- Mixed bright/faint regions — neighborhoods near each other with opposite or non-significant effects. This signals within-cell-type heterogeneity: some sub-states expand while others contract. Flag for deeper annotation.
- Isolated bright dots — single significant neighborhoods surrounded by non-significant ones. Treat with caution; these are more likely false positives after SpatialFDR correction.
The neighborhood graph is the DA analog of a volcano plot: it maps statistical significance onto biological space rather than onto a gene-level axis.
Troubleshooting the neighborhood graph plot. The layout argument must match one of the names in reducedDimNames(milo) — use reducedDimNames(milo) to check what is available before calling the plot function. If the Harmony-based UMAP is stored under "UMAP_HARMONY" in the Seurat object, it should transfer to the SCE under the same name; if it was stored under "umap.harmony" (lowercase), that is what to pass. A common error is the reduction being named "UMAP" (the default Seurat UMAP slot) without the Harmony suffix, which produces a plot on the uncorrected UMAP — the layout will look wrong because samples will cluster by condition rather than by cell type.
8.10.2 The DA beeswarm
plotDAbeeswarm(da_results, group.by = "seurat_annotations") plots each neighborhood’s log-fold-change as a dot, grouped on the x-axis by the modal cell-type label of the neighborhood’s constituent cells (assigned by annotateNhoods; see §8.11). The y-axis is the log-fold-change in abundance. Dots are colored by significance (typically SpatialFDR < 0.1).
Reading rules:
- A cell type with most dots above zero → the majority of this cell type’s transcriptional space is enriched in the perturbed condition.
- A cell type with most dots below zero → the majority of this cell type’s space is depleted.
- A cell type with dots scattered above and below zero → within-cell-type heterogeneity: some sub-states expand while others contract. This is the case where neighborhood testing recovers signal that cluster-level testing misses, and it is the most biologically informative pattern — flag and investigate.
For the ifnb STIM-vs-CTRL contrast, the expected result is most CD14 monocyte, CD16 monocyte, and NK cell neighborhoods enriched in STIM (these cell types respond strongly to IFN-β); T-cell and B-cell DA tends to be more modest and variable across donors.
8.10.3 Reading the results table
The da_results data frame from testNhoods has columns described in Table 3. The primary filter is SpatialFDR < 0.1.
testNhoods output table. The key columns for downstream filtering and reporting are logFC and SpatialFDR. SpatialFDR < 0.1 (10% FDR) is the standard threshold because SpatialFDR is inherently conservative due to neighborhood overlap.
| Column | Meaning | How to use it |
|---|---|---|
logFC |
Log-fold-change in neighborhood abundance (perturbed vs reference) | Effect size; positive = enriched in perturbation |
logCPM |
Log counts per million across all samples | Low values flag neighborhoods with very few cells |
F or PValue |
Quasi-likelihood F-statistic or p-value from the NB GLM | Do not filter directly — use SpatialFDR |
SpatialFDR |
FDR-corrected p-value accounting for neighborhood overlap | Primary significance filter: SpatialFDR < 0.1 |
Nhood |
Neighborhood index | Join key to retrieve neighborhood members |
NhoodGroup |
Cluster of spatially adjacent neighborhoods (after groupNhoods) |
Use for summarizing spatially coherent DA regions |
8.11 Annotating neighborhoods
After the test, each neighborhood has a log-FC and SpatialFDR but no human-readable label. annotateNhoods(milo, da_results, coldata_col = "seurat_annotations") assigns each neighborhood the modal cell-type label of its constituent cells — the label held by the plurality of cells in the neighborhood.
The function also adds a mixed_cell_types column: the fraction of cells in each neighborhood that carry a label different from the modal one. Values above ~0.2 indicate a heterogeneous neighborhood straddling a cluster boundary, and these should be interpreted cautiously — the log-FC is a mixture of signals from multiple cell types rather than a pure single-cell-type signal. Filtering to mixed_cell_types < 0.2 before summarizing per-cell-type DA is a reasonable quality filter.
Use the annotation to:
- Generate the beeswarm plot (§8.10.2)
- Produce a per-cell-type count of significant neighborhoods
- Cross-tabulate against your manual annotations from Chapter 3
8.12 propeller — a fast cluster-level alternative
For experiments with high-quality, well-separated cluster annotations and a need for fast, publication-ready cluster-level DA results, propeller6 provides a linear-model test on logit-transformed cell-type proportions. The transformation partially addresses the compositional constraint by moving proportions from the [0,1] simplex to an unconstrained real-valued scale before fitting the linear model.
The workflow is minimal. The high-level propeller() wrapper logit-transforms the per-sample proportions and runs the appropriate test automatically — a moderated t-test for two groups, ANOVA for more:
library(speckle)
results <- propeller(clusters = seu$seurat_annotations,
sample = seu$donor,
group = seu$stim)If you need to adjust for additional covariates (e.g. sex or batch), call the lower-level pair directly: getTransformedProps() to build the transformed proportions, then propeller.ttest() with a design matrix (non-intercept format) and a contrasts vector identifying the two groups — see the speckle vignette for the design-matrix construction. propeller.ttest() has no covariates argument; covariates enter as extra columns of design.
The output is a table of cluster-level p-values, log-fold-changes in transformed proportion, and BH-adjusted p-values. propeller’s limitation is that it is cluster-level: it inherits all the boundary-sensitivity problems described in §8.2. Its advantage is speed and interpretability: reviewers who find neighborhood graphs confusing may find a cluster-level table straightforward. Use it for validation and for producing the “cell-type composition summary” that many journals require in a methods table.
8.13 scCODA — Bayesian compositional analysis
scCODA2 takes a different philosophical approach. Rather than testing each cluster or neighborhood independently, it models the full composition vector of each sample as a draw from a Dirichlet-multinomial distribution. The key features are:
- Full compositional modeling. All cell-type proportions are modeled jointly, so the constraint that they sum to one is respected exactly.
- A reference category. The model requires one cell type to serve as the reference against which others are compared. Typically, choose a cell type expected to be stable between conditions.
- Bayesian posterior inference. Rather than p-values, the output is the posterior inclusion probability (PIP) for each cell type’s condition coefficient. A PIP > 0.95 is the standard evidence threshold for DA.
- Particularly powerful for small cohorts. With only 3–4 donors per condition, frequentist tests have very low power; scCODA’s prior distributions provide regularization that prevents overfit.
scCODA is available as a Python package. For R workflows, the TreeSummarizedExperiment interface supports import/export to Python via reticulate, or results can be obtained from the original Python implementation.
Choosing a reference cell type in scCODA. The reference cell type is the one whose coefficient is fixed at zero; all other cell types’ DA effects are expressed relative to it. The choice matters: if the reference cell type itself changes in abundance, every other type’s estimate will be biased. A practical rule is to choose the cell type with the most stable proportion across all samples and conditions — check with a simple boxplot of proportions per condition. If no cell type is clearly stable, scCODA’s tree-structured prior (TreeSCCODA) can use a cell-type hierarchy to assign the reference more robustly. Alternatively, if you have a priori biological knowledge that a specific cell type is unaffected by your perturbation (e.g. mature erythrocytes in a monocyte activation experiment), that is the natural reference.
When cell-type proportions in a sample must sum to 1 (as they always do), a change in one cell type mathematically forces a compensating change in all others — even if those other cell types are biologically unchanged in absolute cell number. This compositional constraint means that naive proportion tests can report spurious negative DA results for cell types that are merely relatively smaller because another type expanded. Methods that test on raw cell counts per neighborhood (miloR) or model the full composition vector jointly (scCODA) avoid this artifact. Methods that test individual proportions independently (a simple t-test on fractions) are most affected. The compositional problem is best known from microbiome analysis, where it motivated the development of log-ratio transformations — the same mathematical framework applies directly to cell-type proportions in single-cell DA.
8.14 Common errors
Table 4 collects the errors most frequently encountered in the DA workflow.
| You see | What’s wrong | What to do |
|---|---|---|
reducedDimNames(milo) doesn’t include HARMONY |
SCE conversion dropped reductions | Re-add: reducedDim(milo, "HARMONY") <- Embeddings(seu, "harmony") |
| Neighborhood size histogram peaks at <30 | Neighborhoods too small | Increase k to 50 or 75 |
testNhoods returns all NA p-values |
No replication in design | Need ≥ 2 samples per condition; with 1 vs 1, no test possible |
plotNhoodGraphDA errors “UMAP_HARMONY not found” |
Dimred name mismatch | Check reducedDimNames(milo); pass the correct name to layout = |
| All neighborhoods show similar log-FC | Possibly real (uniform shift), or possibly under-correction in integration | Check Chapter 5 over-correction diagnostics; if integration is fine, trust the result |
countCells errors about missing sample column |
Sample column name mismatch between SCE colData and argument | Check colnames(colData(milo)) and pass the exact column name |
| Very few significant neighborhoods despite expected biology | k too large (over-smoothing) or poor integration |
Decrease k; re-check integration diagnostics |
8.15 Connecting DA to DE
The differential abundance and differential expression analyses are complementary lenses on the same dataset. After running both, three integration strategies produce a coherent narrative:
Overlay cell types. Which cell types show both significant DA and significant DE? These are the most strongly affected populations and deserve the most detailed follow-up. Which cell types show DA but no DE? Those populations change in number but not in transcriptional state — possibly a developmental or trafficking effect rather than a transcriptional perturbation. Which cell types show DE but no DA? Those populations are transcriptionally activated but their absolute numbers are unchanged — a common pattern for mature, long-lived cell types responding to a stimulus without dividing or being recruited.
Test whether DA-enriched neighborhoods have distinct expression. findNhoodMarkers(milo, da_results, assay = "logcounts") performs a differential expression test comparing DA-significant neighborhoods against background, directly within the miloR framework. This recovers the transcriptional signature of the expanding or shrinking cell states without relying on hard cluster labels. The advantage over pseudobulk DE per annotated cell type is that findNhoodMarkers works at neighborhood resolution and can identify transcriptional differences within a broadly annotated cell type — for example, the expanding monocyte subpopulation that is DA-positive might show higher expression of a specific activation program compared to the stable monocyte background.
Directional consistency. If a cell type is DA-enriched in the STIM condition and also shows up-regulated type I interferon response genes by DE, the two analyses are mutually reinforcing. If a cell type is DA-depleted but shows up-regulated survival/proliferation genes, that is a more complex picture — the fewer surviving cells may be transcriptionally altered, suggesting a selection effect. If a cell type is DA-enriched in one direction and DE-enriched in the opposite direction (i.e. the expanding population has lower average expression of a pathway), the most common explanation is cell maturation state: the newly recruited or proliferating cells are at an earlier developmental stage with a different transcriptional profile from the resident population.
A practical figure to include. A useful two-panel composite figure pairs the miloR beeswarm (DA effect per cell type) with a heatmap of pseudobulk LFCs for the top DE genes per cell type. The beeswarm provides the composition story; the heatmap provides the transcriptional story; together they constitute the complete condition-level characterization that high-quality single-cell manuscripts now require.
Reporting both analyses together — either in a single combined figure or as adjacent supplementary figures — is now expected in high-quality single-cell manuscripts.
8.16 Where this material is also discussed
- Lecture: Lecture 08
- Tutorial: Tutorial 08
- DE complement: Chapter 6 — DESeq2: Bulk Fundamentals + Pseudobulk DE, Lecture 06
- Integration prerequisite: Chapter 5
- Annotation prerequisite: Chapter 3 — Markers & Annotation
8.17 Going further
The tutorial’s method is Milo5; its main alternatives are scCODA2, sccomp3, and propeller6, and the benchmark of six DA methods4 is the reference for choosing among them. The pseudoreplication caution carries over directly from differential expression1. For the compositional problem in the microbiome context (where it is most thoroughly developed), the statistical foundations are directly applicable to single-cell DA. The full curated list is on Key Papers, Reviews & Benchmarks.