---
title: "Tutorial 04 — Reference-based Annotation (Azimuth)"
subtitle: "Map cells against a labeled reference for fast, reproducible cell-type calls"
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 hands-on companion to **[Lecture 04 — Reference-based Annotation](../Lecture_Folder/Lecture_04_Reference_Annotation.html)**. We start from the **annotated object saved at the end of Tutorial 03** (`ifnb_annotated.rds`) — reference-based methods like Azimuth run per cell and don't need your own integration, so this step slots in right after manual annotation and before integration. We:


::: callout-note
**Companion book chapter:** [Chapter 4 — Reference Annotation](../Resources_Folder/Chapter_04_Reference_Annotation.html) — the long-form prose treatment of this tutorial's material, with cross-references to the prerequisite appendices. · **HPC version (Talapas):** [Talapas analysis pipeline — `04_reference_annotation.R`](Tutorial_10_Talapas_Pipeline.html)
:::


-   Run `Azimuth` against the PBMC reference (a perfect fit for this dataset)
-   Compare the automatic calls against the manual labels you assigned in Tutorial 03 *and* against the author-curated `seurat_annotations` ground truth
-   Reconcile disagreements and save a final, reference-validated annotation

::: callout-note
**Why this works particularly cleanly for `ifnb`.** The dataset is human PBMCs — exactly what the Azimuth `pbmcref` reference is built for. You should expect **high agreement** between the automatic labels and your manual ones. That's a feature, not a coincidence: it's the easy regime, where reference-based annotation shines. In the **bonus raw-data track** — the [Full Talapas run from raw FASTQs](../Resources_Folder/Talapas_FullRun_FromFASTQs.html) on the NSCLC tumor dataset — you'll see what happens when the data includes tumor / non-immune cells outside the reference's support: the harder regime.
:::

::: {.callout-tip title="How to use this page"}
The rendered HTML shows the code but does **not** execute it. To run it:

1.  Download the `.qmd` source: <a href="Tutorial_04_Reference_Annotation.qmd.txt" download="Tutorial_04_Reference_Annotation.qmd">Tutorial_04_Reference_Annotation.qmd</a>. If your browser saves the file as `Tutorial_04_Reference_Annotation.qmd.txt`, **drop the trailing `.txt`** so the filename ends in `.qmd`, then open it in RStudio.
2.  Make sure you've completed **[Tutorial 03 — Markers & Manual Annotation](Tutorial_03_Markers_Annotation.html)** first — it writes `ifnb_annotated.rds`.
3.  Work through the chunks, flipping `eval: true` as you go.
:::


::: {.callout-warning title="Common errors / things that bite"}
**`RunAzimuth` errors during reference download** — on first call it installs the `pbmcref.SeuratData` package (~75 MB) from the SeuratData repo. Confirm internet and disk. If it fails midway, retry; you can also install it directly with `SeuratData::InstallData("pbmcref")`.

**Lots of warnings while `RunAzimuth` runs (all expected)** — a successful run is *noisy*. You'll see `944 features … not present … Continuing with remaining 4056` (normal — your data has 4056 of the reference's model genes, which is plenty), `Feature names cannot have underscores … replacing with dashes` (from Azimuth's **own** internal assays, not your data), plus `refUMAP`/`integrated_dr_` key and `RunUMAP method changed` notices. None are errors. The tell that it worked: `detected … HUMAN with id type Gene.name` for **both** query and reference, `Found … anchors`, and new `predicted.celltype.l1/l2/l3` + `mapping.score` columns in `ifnb@meta.data`.

**STIM cells get lower mapping scores than CTRL** — that's expected, not a bug: the interferon-β response shifts cells off the resting-PBMC reference manifold even though their *cell-type* identity is unchanged. Read it as a confidence flag, not a wrong label.

**Many cells labeled with a coarse type** — `predicted.celltype.l1` is broad by design. Use `predicted.celltype.l2` (or `l3`) for finer resolution.
:::


::: {.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_04_Reference_Annotation.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: M4-setup
#| fig-cap: "Manual cell-type annotations (from Tutorial 03) on the unintegrated UMAP — the starting point for reference-based annotation."
library(Seurat)
library(tidyverse)
library(patchwork)

set.seed(2026)

# ---------------------------------------------------------------------------
# Output directory for this module's figures and tables.
# Every figure/table chunk below writes a file named  Mod4_C<chunk>_<name>
# into ../output/Mod4/ so it can be cross-referenced from the rest of the site.
#   Mod4 = Module 4 (this tutorial);  C<n> = the nth code chunk.
# ---------------------------------------------------------------------------
out_dir <- "../output/Mod4"
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))

# --- 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_annotated.rds")
DefaultAssay(ifnb) <- "RNA"

p_setup_umap <- DimPlot(ifnb, group.by = "seurat_annotations", reduction = "umap",
        label = TRUE, repel = TRUE) +
  NoLegend() +
  labs(
    title = "Manual cell-type annotations on the unintegrated UMAP",
    subtitle = "Author-curated seurat_annotations carried from Tutorial 03",
    x = "UMAP 1",
    y = "UMAP 2"
  )
p_setup_umap

save_fig(file.path(out_dir, "Mod4_C1_manual_annotations_umap.png"), p_setup_umap,
       width = 7, height = 6, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Each point is a cell, positioned by the unintegrated UMAP from Tutorial 02–03, and coloured/labelled by the author-curated `seurat_annotations`. You should see well-separated clusters for the major PBMC types — monocytes, T cells, B cells, NK cells, dendritic cells — each occupying its own region of the embedding. This is the baseline you'll compare Azimuth's automatic labels against in the next steps.
:::

## Azimuth — Seurat reference mapping

Azimuth wraps Seurat's anchor-based reference mapping with curated atlases, and returns hierarchical labels (`predicted.celltype.l1` / `l2` / `l3`) plus a per-cell **mapping score**. The PBMC reference (`pbmcref`) is the most widely tested and is a perfect fit for `ifnb`.

```{r}
#| label: M4-azimuth
#| fig-cap: "Azimuth (pbmcref) level-2 predicted cell types projected onto the unintegrated UMAP."
library(Azimuth)

# Azimuth (loaded here, *after* tidyverse) and its dependencies mask several dplyr
# verbs. Rebind the ones used below so they resolve to dplyr — otherwise
# `ann |> count(...)` calls the wrong `count` and errors with
# "Argument 'x' is not a vector: list". (Same pattern as Tutorials 07/08/13.)
count  <- dplyr::count
select <- dplyr::select
filter <- dplyr::filter
arrange <- dplyr::arrange

# Installs the pbmcref.SeuratData reference package on first use (~75 MB download)
ifnb <- RunAzimuth(ifnb, reference = "pbmcref")

# New columns: predicted.celltype.l1 / l2 / l3 + scores
head(ifnb@meta.data[, grep("predicted", colnames(ifnb@meta.data))])

# --- Table out: per-cell Azimuth predictions (l1/l2/l3 + scores) -----------
ifnb@meta.data[, grep("predicted", colnames(ifnb@meta.data)), drop = FALSE] |>
  tibble::rownames_to_column("cell") |>
  readr::write_csv(file.path(out_dir, "Mod4_C2_azimuth_predictions.csv"))

p_azimuth_umap <- DimPlot(ifnb, group.by = "predicted.celltype.l2", reduction = "umap",
        label = TRUE, repel = TRUE) +
  NoLegend() +
  labs(
    title = "Azimuth predicted cell types (level 2) on the unintegrated UMAP",
    subtitle = "Anchor-based mapping against the Azimuth pbmcref reference",
    x = "UMAP 1",
    y = "UMAP 2"
  )
p_azimuth_umap

save_fig(file.path(out_dir, "Mod4_C2_azimuth_labels_umap.png"), p_azimuth_umap,
       width = 7, height = 6, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
The `head()` table shows each cell's Azimuth predictions at three granularity levels (`predicted.celltype.l1/l2/l3`) and the per-cell `mapping.score` (0–1): higher means the cell's neighborhood projects cleanly into the reference manifold. In the UMAP, each point is coloured by its `predicted.celltype.l2` label; for a clean PBMC dataset like `ifnb` you should see the major types (`CD14 Mono`, `CD4 T`, `CD8 T`, `B naive`, `NK`, `DC`) mapping onto the same clusters you labelled manually in Tutorial 03. Large discrepancies between the UMAP positions of Azimuth labels and your manual labels indicate either an off-manifold cell type or a granularity mismatch between `l1`/`l2`/`l3`.
:::

### Reading the mapping score — a built-in confidence check

```{r}
#| label: M4-azimuth_confidence
#| fig-cap: "Azimuth mapping scores on the UMAP. Low-scoring cells project poorly into the reference manifold — treat their labels with caution."

p_mapping <- FeaturePlot(ifnb, features = "mapping.score", reduction = "umap") +
  labs(
    title = "Azimuth mapping score per cell",
    subtitle = "High (~0.7–1.0) = projects cleanly into the reference; low = off-manifold",
    x = "UMAP 1",
    y = "UMAP 2"
  )
p_mapping

save_fig(file.path(out_dir, "Mod4_C3_azimuth_mapping_score.png"), p_mapping,
       width = 7, height = 6, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
The colour scale runs from low (dark/cool) to high (bright/warm) mapping score, overlaid on the same UMAP positions as before. Most cells in a clean PBMC dataset should score ≥ 0.7; cells that fall below ~0.5 project poorly into the `pbmcref` manifold and their type labels deserve extra scrutiny. In `ifnb`, watch for IFN-β-stimulated STIM cells showing systematically lower scores than their CTRL counterparts — this is expected because the resting-PBMC reference doesn't model the activated state, even though the cell-type identity is unchanged.
:::

::: {.callout-important title="Think about it"}
Azimuth's PBMC reference is built for **resting** PBMCs. The `ifnb` STIM cells are **interferon-β stimulated**. Will Azimuth confidently place STIM cells, or might it struggle? How would you spot a problem?

<details>

<summary>Click for answer</summary>

Azimuth still assigns the best-matching label, but the **mapping score** (`mapping.score`) reflects how well the cell's neighborhood projects into the reference space. In the plot above, if STIM cells have **systematically lower mapping scores** than CTRL cells, that's evidence that the IFN response is shifting STIM cells off the reference manifold even though their *cell-type* identity is unchanged.

In practice the PBMC reference is robust enough that this effect is mild — but it's a great teaching example of how a "good" reference can still under-represent a perturbation state. Compare the means with `tapply(ifnb$mapping.score, ifnb$stim, mean)`.

</details>
:::

## Compare Azimuth against the manual / ground-truth labels

[]{#lst-M4-compare}

```{r}
#| label: M4-compare

ann <- ifnb@meta.data |>
  select(seurat_annotations,                  # ground truth (Kang authors)
         predicted.celltype.l2)               # Azimuth

table(ann$seurat_annotations, ann$predicted.celltype.l2, useNA = "ifany")

# --- Table out: ground-truth vs Azimuth cross-tabulation -------------------
as.data.frame(
  table(seurat_annotations = ann$seurat_annotations,
        predicted_celltype_l2 = ann$predicted.celltype.l2, useNA = "ifany"),
  responseName = "n_cells"
) |>
  readr::write_csv(file.path(out_dir, "Mod4_C4_truth_vs_azimuth.csv"))
```

::: {.callout-tip title="Reading the output"}
The `table()` printout ([code](#lst-M4-compare)) is a contingency matrix with ground-truth cell types on rows and Azimuth `predicted.celltype.l2` labels on columns; each cell value is the number of cells with that combination. For a clean PBMC dataset like `ifnb`, most cells should land in one off-diagonal cell per row — you want large numbers concentrated near the "correct" column. Zero rows or columns indicate a label present in one system but not the other. The heatmap below ([code](#lst-M4-compare_heatmap)) makes patterns easier to read at a glance.
:::

A cleaner cross-tabulation as a heatmap ([code](#lst-M4-compare_heatmap)):

[]{#lst-M4-compare_heatmap}

```{r}
#| label: M4-compare_heatmap
#| fig-cap: "Confusion heatmap of author ground-truth labels (rows) against Azimuth level-2 predictions (columns)."

p_compare_heatmap <- ann |>
  count(seurat_annotations, predicted.celltype.l2) |>
  ggplot(aes(predicted.celltype.l2, seurat_annotations, fill = n)) +
  geom_tile() +
  scale_fill_viridis_c() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Ground truth (rows) vs. Azimuth (cols)",
       subtitle = "Tile colour = number of cells with that (truth, Azimuth) label pair",
       x = "Azimuth predicted.celltype.l2",
       y = "Author seurat_annotations",
       fill = "Cells")
p_compare_heatmap

save_fig(file.path(out_dir, "Mod4_C5_truth_vs_azimuth_heatmap.png"), p_compare_heatmap,
       width = 8, height = 6, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
The rows are the author ground-truth labels (`seurat_annotations`) and the columns are Azimuth's `predicted.celltype.l2` calls; tile colour encodes the number of cells with that combination. A well-behaved mapping shows **bright tiles on or near the diagonal**, meaning Azimuth agrees with the ground truth for most cells. Off-diagonal bright tiles mark genuine disagreements — common cases are coarse-vs-fine granularity mismatches (e.g. `B cells` vs `B naive`/`B memory`) rather than biological misclassification. Cells in rows with no single bright tile dominating are the ones most worth reviewing.
:::

## Reconcile disagreements

A simple "agree-or-flag" rule: keep the ground-truth label where Azimuth agrees; otherwise mark for review ([code](#lst-M4-reconcile)).

::: callout-important
**You can't compare the labels as raw strings.** The author ground truth (`seurat_annotations`, from `muscData`) and Azimuth use **different vocabularies and granularities** — e.g. `CD14+ Monocytes` (truth) vs `CD14 Mono` (Azimuth l2) vs `Mono` (Azimuth l1); `B cells` vs `B naive`/`B memory`. A direct `==` would flag *everything* as a disagreement even when the calls actually agree. So we **harmonize both to a common coarse vocabulary** and reconcile at Azimuth **level 1** (`predicted.celltype.l1`), the granularity that matches the coarse ground truth.
:::

[]{#lst-M4-reconcile}

```{r}
#| label: M4-reconcile

# Map the author ground-truth labels onto Azimuth's coarse level-1 vocabulary
# (B, CD4 T, CD8 T, NK, DC, Mono). Anything without a clean level-1 counterpart
# (e.g. Megakaryocytes — platelets aren't a level-1 class) maps to NA and is left
# for review. Adjust this crosswalk if your reference/labels differ.
to_coarse <- function(x) dplyr::recode(as.character(x),
  "B cells"           = "B",
  "CD4 T cells"       = "CD4 T",
  "CD8 T cells"       = "CD8 T",
  "NK cells"          = "NK",
  "Dendritic cells"   = "DC",
  "CD14+ Monocytes"   = "Mono",
  "FCGR3A+ Monocytes" = "Mono",
  "Megakaryocytes"    = "Mk",
  .default            = NA_character_)

truth_coarse <- to_coarse(ifnb$seurat_annotations)
agree <- !is.na(truth_coarse) &
         truth_coarse == as.character(ifnb$predicted.celltype.l1)

ifnb$celltype_final <- ifelse(
  agree,
  as.character(ifnb$seurat_annotations),
  paste0("REVIEW:", as.character(ifnb$seurat_annotations))
)

# Fraction of cells where the coarse calls agree (expect high — this is a clean
# PBMC dataset against a PBMC reference), then the final call counts.
round(mean(agree), 3)
table(ifnb$celltype_final)

# --- Table out: final reconciled cell-type call counts ---------------------
tibble::enframe(table(ifnb$celltype_final),
                name = "celltype_final", value = "n_cells") |>
  dplyr::mutate(n_cells = as.integer(n_cells)) |>
  dplyr::arrange(dplyr::desc(n_cells)) |>
  readr::write_csv(file.path(out_dir, "Mod4_C6_celltype_final_counts.csv"))
```

::: {.callout-tip title="Reading the output"}
The first printed value is the **agreement fraction** — the proportion of cells where the coarsened ground-truth label matches Azimuth's level-1 call. For a clean PBMC dataset against the PBMC reference you should see ≥ 0.90; values below \~0.80 suggest vocabulary mismatches or genuine off-reference populations. The `table(ifnb$celltype_final)` below it shows how many cells kept their ground-truth label vs how many were prefixed with `REVIEW:` — the `REVIEW:` entries are your work queue for the next annotation check.
:::

::: {.callout-important title="Think about it"}
A reviewer asks you to "use Azimuth labels instead of your manual ones because they're more reproducible". How do you respond?

<details>

<summary>Click for answer</summary>

Reproducibility is necessary but not sufficient. Azimuth is reproducible *and* biased by its reference: any cell type absent from the reference will be silently re-labeled. The right answer is to use both, as the tutorial does — manual marker-based labels and the automatic reference call — and document where they disagree, not to pick one because it's mechanical.

</details>
:::

## Save the final annotated object

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

saveRDS(ifnb, "../data/ifnb_annotated_final.rds")
```

## Wrap-up

You now have:

-   An Azimuth label per cell (with a mapping score)
-   A reconciled `celltype_final` you can carry into Tutorial 06 (DESeq2 + pseudobulk DE) or downstream analyses

::: {.callout-tip title="Going further"}
-   Try a tissue-appropriate reference: another Azimuth atlas (`lungref`, `bonemarrowref`, …), the [Tabula Sapiens](https://tabula-sapiens-portal.ds.czbiohub.org/) atlas, or the [HCA Lung Cell Atlas](https://data.humancellatlas.org/).
-   For a second opinion from a different annotation paradigm, try a **correlation-based** method (`SingleR` with a `celldex` reference) or a **classifier-based** one (**CellTypist**, Python). These aren't installed for the workshop — Azimuth alone is enough here — but they're worth knowing; see [Chapter 4](../Resources_Folder/Chapter_04_Reference_Annotation.html) for how they compare.
-   For tumor data, look at **TumorDecon** or **ProjecTILs** to handle the tumor/immune mix more carefully — see also the parallel reference-annotation example in the **Talapas Advanced series**.
:::

## See also

-   [Lecture 04 — Reference-based Annotation](../Lecture_Folder/Lecture_04_Reference_Annotation.html)
-   [Lecture 03 — Markers & Manual Annotation](../Lecture_Folder/Lecture_03_Markers_Annotation.html)
-   [Glossary — Cell type annotation](../Resources_Folder/Glossary.html#c)
-   Azimuth web tool: [azimuth.hubmapconsortium.org](https://azimuth.hubmapconsortium.org/)