---
title: "Tutorial 08 — Differential Abundance with miloR"
subtitle: "Do cell-type proportions change between STIM and CTRL? — neighbourhood-level testing"
author: "Single Cell RNA-seq Workshop"
format:
  html:
    toc: true
    toc-depth: 3
    code-fold: false
    code-overflow: wrap
    highlight-style: github
    embed-resources: true
execute:
  eval: false
  echo: true
  warning: false
  message: false
editor: visual
---

::: {.callout-note title="Running 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](Tutorial_06_DESeq2_DE.html)**. 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.


::: callout-note
**Companion lecture:** [Lecture 08 — Differential Abundance: Do Cell-Type Proportions Change?](../Lecture_Folder/Lecture_08_DifferentialAbundance.html) · **Companion book chapter:** [Chapter 8 — Differential Abundance](../Resources_Folder/Chapter_08_DifferentialAbundance.html) — 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`](Tutorial_10_Talapas_Pipeline.html)
:::


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.

::: callout-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.
:::

::: {.callout-tip title="How to use this page"}
1.  Download the `.qmd` source: <a href="Tutorial_08_DifferentialAbundance.qmd.txt" download="Tutorial_08_DifferentialAbundance.qmd">Tutorial_08_DifferentialAbundance.qmd</a>. 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](Tutorial_05_Integration.html)** — 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. |


::: {.callout-warning title="Common 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.
:::


::: {.callout-tip title="Solutions / 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](https://github.com/wcresko/scRNAseq_tutorial/blob/main/Exercise_Folder/_quarto-solutions.yml) for the build profile
:::

## Setup

```{r}
#| label: M8-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](Tutorial_01_QC_Preprocessing.html) 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.

```{r}
#| label: M8-donor

# 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"))
```

::: callout-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:

```{r}
#| label: M8-milo_build

# 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.

[]{#lst-M8-build_graph}

```{r}
#| label: M8-build_graph
#| fig-cap: "Distribution of neighbourhood sizes after refined sampling (target peak 50-200 cells)."

# 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)
```

::: {.callout-tip title="Reading the output"}
The neighbourhood-size histogram ([code](#lst-M8-build_graph)) 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`.
:::

::: {.callout-important title="Think about it"}
What does the neighbourhood-size histogram tell you, and what should you do if the peak is at \~10–20?

<details>

<summary>Show answers</summary>

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.

</details>
:::

## Step 4 — Count cells per neighbourhood per sample

```{r}
#| label: M8-count_cells

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

```{r}
#| label: M8-design

# 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

```{r}
#| label: M8-annotate

# 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"))
```

::: {.callout-important title="Think 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?

<details>

<summary>Show answers</summary>

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.

</details>
:::

## Step 7 — Visualize: neighbourhood graph

The signature miloR plot ([code](#lst-M8-plot_graph)). 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):

[]{#lst-M8-plot_graph}

```{r}
#| label: M8-plot_graph
#| fig-cap: "Cell-type UMAP (left) alongside the DA neighbourhood graph (right); 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)
```

::: callout-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)`.
:::

::: {.callout-tip title="Reading the output"}
The two-panel plot ([code](#lst-M8-plot_graph)) 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](#lst-M8-beeswarm)) that strips out the spatial graph and asks "which annotated cell type is enriched / depleted on average?":

[]{#lst-M8-beeswarm}

```{r}
#| label: M8-beeswarm
#| fig-cap: "Beeswarm of per-neighbourhood log-fold change grouped by annotated cell type (STIM vs CTRL)."

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)
```

::: {.callout-tip title="Reading the output"}
The beeswarm ([code](#lst-M8-beeswarm)) 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.
:::

::: {.callout-important title="Think 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?

<details>

<summary>Show answers</summary>

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.

</details>
:::

## Step 9 — Save the DA result

```{r}
#| label: M8-save

# 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.

::: {.callout-tip title="Going 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

-   Dann *et al.* 2022 — miloR paper, [doi:10.1038/s41587-021-01033-z](https://doi.org/10.1038/s41587-021-01033-z)
-   [scNotebooks Module 04](https://integrativebioinformatics.github.io/scNotebooks/modules/en/module04/module04.html) — covers miloR on COVID-19 mild vs severe BAL data
-   [Tutorial 06 — DESeq2: Bulk Primer + Pseudobulk DE](Tutorial_06_DESeq2_DE.html) — the per-cell-type DE counterpart
-   [Tutorial 07 — Functional Analysis](Tutorial_07_FunctionalAnalysis.html) — what to do with the DE table