Tutorial 08 — Differential Abundance with miloR

Do cell-type proportions change between STIM and CTRL? — neighbourhood-level testing

Author

Single Cell RNA-seq Workshop

NoteRunning the code in this tutorial

Every code chunk here is tagged with #| eval: false, so the published page shows the code without running it. To run it yourself, open the downloaded .qmd in RStudio:

  • Run one chunk: click the green ▶ (Run Current Chunk) at the chunk’s top-right, or press Ctrl/Cmd + Enter. This runs the chunk regardless of its eval setting — the easiest way to work through the tutorial interactively.
  • Run a chunk on render: change that chunk’s #| eval: false to #| eval: true (or delete the line) so it executes when the document is rendered.
  • Render the whole tutorial: click the Render button in RStudio (or run quarto render in a terminal) to execute the chunks top-to-bottom and knit a finished HTML. To run everything on render, set eval: true once in the YAML header at the top.

About this tutorial

A short companion to Tutorial 06 — DESeq2: Bulk Primer + Pseudobulk DE. Pseudobulk DE asks “for cells of the same type, which genes change?” — but a perturbation can also change the cell-type composition itself (e.g. expansion of activated T cells, depletion of naive subsets). That’s the differential abundance question, and it’s a different test.

Note

Companion lecture: Lecture 08 — Differential Abundance: Do Cell-Type Proportions Change? · Companion book chapter: Chapter 8 — Differential Abundance — the long-form prose treatment of this tutorial’s material, with cross-references to the prerequisite appendices. · HPC version (Talapas): Talapas analysis pipeline — 08_differential_abundance.R

We use miloR (Dann et al. 2022, Nat Biotechnol), the modern method of choice. Instead of testing on hard cluster labels, miloR:

  1. Builds a k-nearest-neighbour graph on the integrated PCs
  2. Defines overlapping neighbourhoods of cells (~1000 of them)
  3. Counts how many cells from each sample fall into each neighbourhood
  4. Tests each neighbourhood for differential abundance with a negative binomial GLM, using a spatial FDR correction to handle the overlap

The output is a per-neighbourhood log-fold change in abundance and a p-value, plus a beautiful neighbourhood graph plot showing where the composition shifts live.

Note

Why neighbourhoods, not clusters? Clusters force a hard partition — every cell belongs to exactly one cluster, and cluster boundaries are arbitrary. A real composition shift along a continuum (naive → effector T cells, say) gets blurred by hard clustering. Neighbourhoods overlap, follow the underlying graph, and recover the resolution that clusters lose.

TipHow to use this page
  1. Download the .qmd source: Tutorial_08_DifferentialAbundance.qmd. If your browser saves the file as Tutorial_08_DifferentialAbundance.qmd.txt, drop the trailing .txt so the filename ends in .qmd, then open it in RStudio.
  2. Make sure you have completed Tutorial 05 — Integration — it writes ifnb_integrated.rds.
  3. Work through the chunks, flipping eval: true as you go.

Dataset — resuming from Tutorial 05

File Produced by What it is
../data/ifnb_integrated.rds Tutorial 05 Seurat object with the integrated harmony reduction, the original RNA assay, the stim covariate, and the seurat_annotations ground-truth labels.
WarningCommon errors / things that bite

could not find function "plotReducedDim" — that function is from scater, which the setup must load (library(scater)). plotNhoodSizeHist / plotNhoodGraphDA come from miloR (already loaded); only the UMAP panel needs scater.

reducedDimNames(milo) doesn’t include HARMONY — the SCE conversion sometimes drops reductions, especially on Seurat v4 → v5 transitions. Check reducedDimNames(milo); if HARMONY isn’t there, re-add it manually: reducedDim(milo, "HARMONY") <- Embeddings(ifnb, "harmony").

makeNhoods peaks the size histogram at < 30 cells — neighborhoods are too small. Increase k from 30 → 50 or 75 in buildGraph and re-run. Median neighbourhood size between 50 and 200 is the target.

testNhoods returns all NA p-values — your milo_design has only 1 level per condition (no replicates). miloR can’t fit the GLM with n = 1. Need biological replicates. The real ind donors give eight replicates (the synthetic fallback gives 4-vs-4); if you have only 2 samples, the test cannot work — this is a fundamental statistical constraint, not a code bug.

plotNhoodGraphDA errors “UMAP_HARMONY not found” — the dimred name in the milo object isn’t UMAP_HARMONY. Try umap_harmony (lowercase) or pass layout = "UMAP". reducedDimNames(milo) shows what’s actually available.

TipSolutions / instructor copy

The Think about it prompts on this page have inline answers (click “Show answer” / “Click for answer” to expand). For the full runnable instructor copy with every chunk pre-evaluated and outputs shown:

  • For students: the same .qmd you downloaded already contains every solution inline — just expand the collapsed answers
  • For instructors: render with quarto render --profile solutions from the Exercise_Folder/ directory. The fully-evaluated HTML is written to docs/Exercise_Folder/_solutions/Tutorial_08_DifferentialAbundance.html. See Exercise_Folder/_quarto-solutions.yml for the build profile

Setup

library(Seurat)
library(SingleCellExperiment)
library(scater)        # plotReducedDim() for the UMAP panel
library(tidyverse)
library(patchwork)
library(miloR)

set.seed(2026)

# ---------------------------------------------------------------------------
# Output directory for this module's figures and tables.
# Every figure/table chunk below writes a file named  Mod8_C<chunk>_<name>
# into ../output/Mod8/ so it can be cross-referenced from the rest of the site.
#   Mod8 = Module 8 (this tutorial);  C<n> = the nth code chunk.
# ---------------------------------------------------------------------------
out_dir <- "../output/Mod8"
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
dir.create("../data",  showWarnings = FALSE)   # shared .rds hand-off folder (sibling of scripts/, created one level up)
message("This module writes its figures & tables to: ", normalizePath(out_dir))

# Prefer dplyr's verbs over identically-named functions from the Bioconductor
# packages loaded above: S4Vectors/IRanges/AnnotationDbi/matrixStats (pulled in by
# org.Hs.eg.db, DESeq2, miloR, SingleCellExperiment, ...) mask select, filter,
# count, desc, slice, rename, first. Re-bind the ones used below so the dplyr
# pipelines work regardless of package load order.
select <- dplyr::select; filter <- dplyr::filter; count <- dplyr::count
desc   <- dplyr::desc;   slice  <- dplyr::slice
rename <- dplyr::rename; first  <- dplyr::first

# --- Figure saving --------------------------------------------------------
# save_fig() writes every figure as a .png (for viewing) AND a .svg (vector,
# editable in Illustrator / Inkscape for a manuscript). Pass the .png path; the
# .svg is written alongside with the same basename. (.svg uses the 'svglite'
# package, installed in Tutorial 00.)
save_fig <- function(filename, plot, width, height, dpi = 300, ...) {
  ggplot2::ggsave(filename, plot, width = width, height = height, dpi = dpi, ...)
  svg_path <- paste0(tools::file_path_sans_ext(filename), ".svg")
  tryCatch(ggplot2::ggsave(svg_path, plot, width = width, height = height, ...),
           error = function(e) message("  (could not write ", basename(svg_path),
                                       " - install 'svglite'? ", conditionMessage(e), ")"))
}

ifnb <- readRDS("../data/ifnb_integrated.rds")
DefaultAssay(ifnb) <- "RNA"
ifnb

Step 1 — Define donor / sample IDs

Like Tutorial 06’s pseudobulk DE, miloR needs biological replicates within each condition. The muscData::Kang18_8vs8() loader from Tutorial 01 carries the real donor IDs in the ind column, so we use the eight real donors as replicates. A synthetic-donor fallback is kept only for objects that somehow arrive without donor metadata.

# Real donors (column `ind`, from the muscData loader) are the biological
# replicates. Fall back to synthetic pseudo-donors only if no donor column is
# present (e.g. an older SeuratData `ifnb`).
if (!"donor" %in% colnames(ifnb@meta.data)) {
  if ("ind" %in% colnames(ifnb@meta.data)) {
    ifnb$donor <- ifnb$ind
  } else {
    ifnb$donor <- paste0(ifnb$stim, "_d",
                         sample(1:4, size = ncol(ifnb), replace = TRUE))
  }
}

# miloR needs a single per-cell sample ID. Combine donor + stim:
ifnb$sample <- paste(ifnb$donor, ifnb$stim, sep = "_")
table(ifnb$sample, ifnb$stim)

# --- Table out: cells per sample x condition -------------------------------
as.data.frame(table(sample = ifnb$sample, stim = ifnb$stim)) |>
  dplyr::rename(n_cells = Freq) |>
  readr::write_csv(file.path(out_dir, "Mod8_C2_cells_per_sample.csv"))
Warning

On the real ind donors — the default with the Tutorial 01 loader — the test is valid. The synthetic-donor fallback only triggers when no donor metadata is present; on synthetic donors the neighbourhood SpatialFDRs are not biologically meaningful, so you’ll see the workflow run correctly but should not interpret them as real findings.

Step 2 — Convert to a Milo object

miloR is built on SingleCellExperiment, so we cast first:

# Drop assay slots we don't need to keep memory low
sce  <- as.SingleCellExperiment(ifnb)
milo <- Milo(sce)
milo

# What reductions did we carry over?
reducedDimNames(milo)

Step 3 — Build the kNN graph and define neighbourhoods

The kNN graph is built on the integrated reduction (harmony), so neighbourhoods reflect biological similarity after batch correction.

# Build the kNN graph on the harmony reduction.
# k = 30, d = 20 are reasonable defaults for ~23k cells.
milo <- buildGraph(milo, k = 30, d = 20, reduced.dim = "HARMONY")

# Sample ~1000 representative neighbourhoods (refined sampling for stability)
milo <- makeNhoods(milo,
                   prop    = 0.1,    # sample 10% of cells as neighbourhood centres
                   k       = 30,
                   d       = 20,
                   refined = TRUE,
                   reduced_dims = "HARMONY")

# Diagnostic: distribution of neighbourhood sizes (want most >50, peak ~100)
p_nhood_size <- plotNhoodSizeHist(milo) +
  labs(
    title    = "Neighbourhood size distribution",
    subtitle = "Most neighbourhoods should hold 50-200 cells for a powered NB test",
    x        = "Neighbourhood size (cells)",
    y        = "Number of neighbourhoods",
    caption  = "Module 8 · Differential Abundance"
  )
p_nhood_size

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod8_C4_nhood_size_hist.png"), p_nhood_size,
       width = 6, height = 4, dpi = 300)
TipReading the output

The neighbourhood-size histogram (code) shows the distribution of how many cells belong to each sampled neighbourhood. The x-axis is neighbourhood size (cells) and the y-axis is the count of neighbourhoods at that size. You want the peak to sit between 50 and 200 cells: neighbourhoods smaller than ~30 have too few cells per sample for the negative-binomial GLM to estimate dispersion reliably, giving low power or NA p-values; very large neighbourhoods (>500) over-smooth the graph and reduce spatial resolution. For ifnb at ~23k cells, a peak around 100 with k = 30 is normal. If your histogram peaks at <30, increase k and re-run from buildGraph.

ImportantThink about it

What does the neighbourhood-size histogram tell you, and what should you do if the peak is at ~10–20?

Show answers

The histogram shows how many cells each sampled neighbourhood contains. If the peak is too small (~10–20), individual neighbourhoods don’t have enough cells per sample for the negative-binomial test to be powered, and you’ll get few or no significant hits. Increase k (try k = 50 or 75) and re-run buildGraph + makeNhoods until the peak sits between 50 and 200. For ifnb at ~23k cells, k = 30 is usually fine.

Step 4 — Count cells per neighbourhood per sample

milo <- countCells(milo,
                   meta.data = data.frame(colData(milo)),
                   sample    = "sample")

# Sanity check: rows = neighbourhoods, cols = samples
head(nhoodCounts(milo))

Step 5 — Build the design matrix and test

# One row per sample, with stim/donor metadata
milo_design <- data.frame(colData(milo)) |>
  distinct(sample, stim, donor)
rownames(milo_design) <- milo_design$sample
milo_design

# Compute neighbourhood-graph distances (needed for SpatialFDR)
milo <- calcNhoodDistance(milo, d = 20, reduced.dim = "HARMONY")

# Test each neighbourhood for differential abundance.
# The contrast is STIM vs CTRL (alphabetical reference: CTRL).
da_results <- testNhoods(milo,
                         design     = ~ stim,
                         design.df  = milo_design,
                         reduced.dim = "HARMONY")

# Top hits
da_top10 <- da_results |>
  arrange(SpatialFDR) |>
  head(10)
da_top10

# --- Table out: top-10 neighbourhoods by SpatialFDR ------------------------
readr::write_csv(da_top10, file.path(out_dir, "Mod8_C6_da_top10.csv"))

The SpatialFDR column is miloR’s analogue of padj — it controls the false discovery rate while accounting for the fact that overlapping neighbourhoods have correlated tests.

Step 6 — Annotate neighbourhoods with cell-type labels

# For each neighbourhood, take the most common cell type among its cells
da_results <- annotateNhoods(milo,
                              da_results,
                              coldata_col = "seurat_annotations")

# Per-cell-type fraction of significantly DA neighbourhoods
da_by_celltype <- da_results |>
  filter(SpatialFDR < 0.1) |>
  count(seurat_annotations, sign = sign(logFC)) |>
  pivot_wider(names_from = sign, values_from = n,
              names_prefix = "logFC_") |>
  rename(up_in_STIM = logFC_1, up_in_CTRL = `logFC_-1`)
da_by_celltype

# --- Table out: DA neighbourhood counts per cell type, by direction --------
readr::write_csv(da_by_celltype,
                 file.path(out_dir, "Mod8_C7_da_by_celltype.csv"))
ImportantThink about it

Pseudobulk DE (Tutorial 06) and differential abundance (this tutorial) sometimes give opposite answers for the same cell type. How can both be right?

Show answers

Pseudobulk DE answers “within cells annotated as cell type X, which genes shift?” Differential abundance answers “do proportions of cell type X change?” These are independent questions. A cell type can be unchanged in proportion but every cell in it shifts its expression (DE-positive, DA-negative — common for IFN-β–stimulated mature cell types) or it can shift dramatically in proportion without changing per-cell expression (DA-positive, DE-negative — common for cell-cycle expansion). Most often you see a mix: e.g. STIM expands monocytes (DA-positive) and every monocyte expresses ISGs (DE-positive). They’re complementary, not redundant.

Step 7 — Visualize: neighbourhood graph

The signature miloR plot (code). Each node is a neighbourhood, sized by the number of cells in it, coloured by log-fold-change (red = enriched in STIM, blue = enriched in CTRL):

milo <- buildNhoodGraph(milo)

# Plot the UMAP and the DA-neighbourhood graph side by side
um  <- plotReducedDim(milo, dimred = "UMAP_HARMONY",
                       colour_by = "seurat_annotations",
                       text_by   = "seurat_annotations") +
        ggtitle("Cell types") +
        guides(colour = "none")

# Note: the DA plot uses the same UMAP layout
da  <- plotNhoodGraphDA(milo, da_results, alpha = 0.1,
                         layout = "UMAP_HARMONY") +
        ggtitle("DA neighbourhoods (red = ↑ STIM, blue = ↑ CTRL)")

p_graph <- (um + da) +
  patchwork::plot_annotation(
    title    = "Differential abundance across the harmony UMAP — STIM vs CTRL",
    subtitle = "Each node is a neighbourhood, sized by cell count, coloured by log-fold change",
    caption  = "Module 8 · Differential Abundance"
  )
p_graph

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod8_C8_da_nhood_graph.png"), p_graph,
       width = 12, height = 5, dpi = 300)
Tip

If plotReducedDim complains that UMAP_HARMONY doesn’t exist, the harmony UMAP didn’t make it into the SCE conversion. Run RunUMAP(ifnb, reduction = "harmony", dims = 1:20, reduction.name = "umap_harmony") before as.SingleCellExperiment(), or pass any other 2-D embedding present in reducedDimNames(milo).

TipReading the output

The two-panel plot (code) shows the same harmony UMAP layout twice. Left: cells coloured by seurat_annotations — use this as the reference to orient yourself (where are the monocytes, T cells, B cells?). Right: the DA neighbourhood graph where each node is a neighbourhood (not an individual cell), node size reflects how many cells are in the neighbourhood, and node colour shows the log-fold change (red = enriched in STIM, blue = enriched in CTRL). Grey nodes are non-significant at alpha = 0.1 (SpatialFDR). For ifnb, you should see clusters of red nodes in the monocyte region (expanded in STIM) and mostly grey nodes in the lymphocyte regions. A patchwork of adjacent red/blue nodes in the same cell-type region indicates within-type heterogeneity that the beeswarm plot (Step 8) will summarise more cleanly.

Step 8 — Beeswarm: DA per cell type

A clean summary (code) that strips out the spatial graph and asks “which annotated cell type is enriched / depleted on average?”:

p_beeswarm <- plotDAbeeswarm(da_results, group.by = "seurat_annotations") +
  labs(
    title    = "Differential abundance per cell type — STIM vs CTRL",
    subtitle = "Each dot is a neighbourhood; above zero = enriched in STIM, below = enriched in CTRL",
    x        = "Cell type",
    y        = "Log-fold change (STIM vs CTRL)",
    colour   = "SpatialFDR"
  )
p_beeswarm

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod8_C9_da_beeswarm.png"), p_beeswarm,
       width = 8, height = 6, dpi = 300)
TipReading the output

The beeswarm (code) places each neighbourhood as a dot at its log-fold change (y-axis), grouped by the most common cell type in that neighbourhood (x-axis). Dot colour encodes SpatialFDR (darker = more significant). Cell types whose dots cluster tightly above zero are consistently enriched in STIM across the entire region; those below zero are depleted. An even spread of all dots across the zero line means no significant abundance shift for that type. Compare these results to your pseudobulk DE (Tutorial 06): monocytes typically show both strong DA (expanded in STIM) and strong DE (ISG induction), while other types may show DE without DA or vice versa. propeller (Phipson et al. 2022) is a lightweight alternative that tests cluster proportions with limma and can cross-validate your top miloR hits.

ImportantThink about it

The beeswarm shows each neighbourhood as a dot at its log-fold-change. A cell type with most dots above zero is enriched in STIM; below zero, depleted. What would a cell type with dots scattered widely on both sides of zero suggest?

Show answers

Within-cell-type heterogeneity: a subset of that cell type is expanding while a different subset is contracting. This is exactly the kind of structure you’d miss with hard cluster-level proportions tests but recover with neighbourhood-level testing. A common example in ifnb is the T-cell compartment: naive T cells often shrink slightly while activated T cells expand — both labelled T at coarse resolution, but moving in opposite directions.

Step 9 — Save the DA result

# Pipeline hand-off: keep the data/ copies so downstream tutorials and the
# Talapas pipeline can pick the DA result and Milo object up unchanged.
write_csv(da_results, "../data/ifnb_milo_da.csv")
saveRDS(milo,         "../data/ifnb_milo.rds")

# --- Table out: full DA result table for this module's outputs --------------
readr::write_csv(da_results, file.path(out_dir, "Mod8_C10_milo_da_results.csv"))

Wrap-up

You’ve now run the two complementary condition-level tests on the same dataset:

Question Tool Output
Within a cell type, which genes change? DESeq2 on pseudobulk (Tut 06) per-cell-type, per-gene log2FC + padj
Do cell-type / state proportions themselves change? miloR (this tutorial) per-neighbourhood log2FC + SpatialFDR

Both should be reported for any condition contrast. They tell different parts of the story.

TipGoing further
  • More than two conditions: miloR supports any DESeq2-style design. Replace ~ stim with ~ donor + stim (paired design), ~ batch + stim, etc.
  • Continuous covariates: if your phenotype is a gradient (severity score, dose, time), pass it as a numeric column in design.df.
  • Other DA tools: scCODA (Python) tests on cluster proportions with a Bayesian framework, more conservative than miloR; propeller (CRAN) is a fast classical method using limma on cluster proportions. Cross-check your top miloR hits with propeller to make sure they’re not artefacts of the neighbourhood definition.

See also