---
title: "Tutorial 03 — scRNA-seq: Markers & Cell Type Annotation"
subtitle: "Find cluster markers, visualize, label by hand and check against ground truth"
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

Part three of the series. We start from the **unintegrated clustered object saved in Tutorial 02** and walk through:


::: callout-note
**Companion lecture:** [Lecture 03 — Markers & Manual Cell-Type Annotation](../Lecture_Folder/Lecture_03_Markers_Annotation.html) · **Companion book chapter:** [Chapter 3 — Markers & Manual Annotation](../Resources_Folder/Chapter_03_Markers_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 — `03_markers_annotation.R`](Tutorial_10_Talapas_Pipeline.html)
:::


-   Finding marker genes per cluster
-   The canonical marker visualizations (Vln / Feature / Heatmap / Ridge / Dot)
-   Manual cell-type annotation with `RenameIdents`
-   Cross-checking your manual labels against the **ground-truth `seurat_annotations`** that ship with `ifnb`
-   Saving the annotated object so Tutorial 04 (Reference Annotation) can pick it up

::: callout-note
We're annotating the **unintegrated** object. Why? Because **markers are cell-type-driven, not batch-driven** — even when an embedding is dominated by a batch effect, the genes that distinguish T cells from B cells are still T-cell genes. Working markers up before integration is a useful sanity check: if a marker gene cleanly identifies a cluster on the unintegrated data, it should still cleanly identify the *biologically equivalent* cluster after integration in Tutorial 05.
:::

::: {.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_03_Markers_Annotation.qmd.txt" download="Tutorial_03_Markers_Annotation.qmd">Tutorial_03_Markers_Annotation.qmd</a>. If your browser saves the file as `Tutorial_03_Markers_Annotation.qmd.txt`, **drop the trailing `.txt`** so the filename ends in `.qmd`, then open it in RStudio.
2.  Make sure you have completed **[Tutorial 02](Tutorial_02_DimReduction_Clustering.html)** — it writes the `ifnb_clustered.rds` that this tutorial loads.
3.  Work through the chunks, flipping `eval: true` as you go.
:::

## Dataset — resuming from Tutorial 02

| File | Produced by | What it is |
|---|---|---|
| `../data/ifnb_clustered.rds` | Tutorial 02 | Seurat object with PCA, SNN graph, clusters at five resolutions, UMAP, plus the `stim` and `seurat_annotations` metadata. |

## Learning goals

-   Use `FindMarkers` / `FindAllMarkers` appropriately
-   Read each of the five standard marker-visualization plots and know when to pick which
-   Assign cell-type labels manually using PBMC canonical markers
-   Validate against the `seurat_annotations` ground truth and quantify agreement


::: {.callout-warning title="Common errors / things that bite"}
**`FindAllMarkers` takes 5+ minutes** — `presto` (installed in the [Tutorial 00 setup](Tutorial_00_Setup_RStudio_Packages.html)) makes it \~10× faster; Seurat detects it automatically.

**ISGs (`ISG15`, `IFI6`, `IFIT3`, `MX1`) top *some* clusters' marker lists** — expected on the *unintegrated* `ifnb` data: the cell types with the strongest interferon response (monocytes especially) split into CTRL vs STIM clusters, so those clusters' top markers are condition genes rather than cell-type genes. Most clusters are still topped by cell-type markers. Tutorial 05 (Integration) collapses the condition split.

**`names(new.cluster.ids) <- levels(ifnb)` errors with "'names' attribute [N] must be the same length as the vector [M]"** — your clustering produced more clusters than the template's label list. **This is real biology, not a bug.** You're clustering the *unintegrated* data, where the interferon-β response is a strong, genome-wide transcriptional signal — so the cell types that respond most (monocytes above all) separate into *two* clusters, a CTRL one and a STIM one, instead of one. That's the very same condition/batch split you visualized on the UMAP in Tutorial 02, now showing up as extra clusters: you commonly end up with ~15–20 clusters for ~13 cell types, often with two clusters of the same type. It's expected, and Tutorial 05's integration aligns CTRL and STIM so each cell type collapses back to a single cluster. Practically: the Step 4 template no longer uses that brittle line — it maps labels by cluster id and applies only the ids that exist, so it works for any cluster count. Run `levels(ifnb)` to see how many you have, and label each (giving the CTRL and STIM copies of a type the same label is fine).

**`Idents()` is set to the wrong column** — clusters at multiple resolutions all live in `meta.data` (`RNA_snn_res.0.1`, `_0.3`, etc.) but only one is the "active" identity. `Idents(seu) <- "RNA_snn_res.0.5"` switches; `levels(Idents(seu))` confirms.

**`FindMarkers` errors with "data layers are not joined. Please run JoinLayers"** — Seurat v5 keeps per-sample split layers (`counts.CTRL`, `counts.STIM`, `data.CTRL`, `data.STIM`) after `merge()`, and `FindMarkers` refuses to run until they're collapsed. Fix with `ifnb <- JoinLayers(ifnb)` before calling `FindMarkers` / `FindAllMarkers`. The resume chunk above does this automatically. `Layers(ifnb)` shows the current layer list.
:::


::: {.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_03_Markers_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: M3-setup

library(Seurat)
library(tidyverse)
library(patchwork)   # combine + annotate multi-panel Seurat plots

set.seed(2026)

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

```{r}
#| label: M3-resume

ifnb <- readRDS("../data/ifnb_clustered.rds")
Idents(ifnb) <- "RNA_snn_res.0.5"

# In Seurat v5, merging per-sample objects (as Tutorial 01 does after
# scDblFinder) leaves the RNA assay with split layers. FindMarkers requires
# a single joined layer, so rejoin here if it hasn't been done already.
if (length(Layers(ifnb, search = "counts")) > 1) {
  ifnb <- JoinLayers(ifnb)
}

ifnb
```

------------------------------------------------------------------------

## Step 1 — Find markers for one cluster

```{r}
#| label: M3-markers_one

cluster1.markers <- FindMarkers(ifnb, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, 10)

# --- Table out: top markers for cluster 1 (gene moved from rownames) -------
cluster1.markers |>
  tibble::rownames_to_column("gene") |>
  readr::write_csv(file.path(out_dir, "Mod3_C3_cluster1_markers.csv"))
```

::: {.callout-tip title="Reading the output"}
Each row is a gene; the five columns are: `p_val` (raw Wilcoxon p-value), `avg_log2FC` (log₂ fold-change, cluster 1 vs all other cells), `pct.1` (fraction of cluster-1 cells expressing the gene), `pct.2` (fraction of all other cells expressing it), and `p_val_adj` (Bonferroni-corrected p-value). The best *specific* markers have high `avg_log2FC`, high `pct.1` (expressed in most of your cluster), and low `pct.2` (rare outside it) — genes with high `pct.2` may be significant but are not cluster-specific. Rows are sorted by p-value by default, not fold-change, so scroll or re-sort by `avg_log2FC` to find the strongest effect-size markers.
:::

::: {.callout-important title="Think about it"}
1.  What does `min.pct = 0.25` do, and why is that a useful default?
2.  The output has `avg_log2FC`, `pct.1`, `pct.2`. Which gene would you pick as a *specific* marker for this cluster?

<details>

<summary>Show answers</summary>

1.  Only tests genes detected in ≥ 25% of cells in at least one of the two groups. Prunes low-expression noise and speeds the test up significantly.
2.  A specific marker has **high `pct.1` and low `pct.2`** (expressed in almost all cluster-1 cells, almost none elsewhere) *and* high positive `avg_log2FC`. Markers that are broadly expressed (high `pct.2`) may be real but aren't *specific*.

</details>
:::

## Step 2 — Markers for every cluster

```{r}
#| label: M3-markers_all

ifnb.markers <- FindAllMarkers(
  ifnb,
  only.pos        = TRUE,
  min.pct         = 0.25,
  logfc.threshold = 0.25
)

top_markers <- ifnb.markers |>
  group_by(cluster) |>
  slice_max(avg_log2FC, n = 5)
top_markers

# --- Tables out: all positive markers + the per-cluster top 5 --------------
readr::write_csv(ifnb.markers, file.path(out_dir, "Mod3_C4_all_markers.csv"))
readr::write_csv(top_markers,  file.path(out_dir, "Mod3_C4_top5_markers_per_cluster.csv"))
```

::: {.callout-tip title="Reading the output"}
`top_markers` is a tibble grouped by `cluster`, showing the five genes with the highest `avg_log2FC` in each cluster. Scan the `gene` column cluster-by-cluster: most clusters will be topped by recognizable PBMC markers (`LYZ`/`S100A8` for monocytes, `MS4A1`/`CD79A` for B cells, `GNLY`/`NKG7` for NK, `IL7R`/`CCR7` for naive T). If a cluster's top genes are interferon-stimulated genes (`ISG15`, `IFI6`, `IFIT3`, `MX1`) rather than cell-type markers, that cluster is likely a condition-split copy of a cell type rather than a biologically distinct population — expected on unintegrated data and resolved in Tutorial 05.
:::

::: {.callout-important title="Think about it"}
1.  Most of the top per-cluster markers are **cell-type** genes (e.g. `LYZ`/`S100A8` for monocytes, `MS4A1`/`CD79A` for B cells, `GNLY`/`NKG7` for NK). But for *some* clusters — especially the **monocytes** — you'll also see **interferon-stimulated genes** (`ISG15`, `IFI6`, `IFIT3`, `MX1`) near the top. If those ISGs aren't cell-type markers, what are they marking?
2.  For condition-level DE (STIM vs CTRL), is `FindAllMarkers` the right tool?

<details>

<summary>Show answers</summary>

1.  Those ISGs are **condition markers** — they distinguish IFN-β–stimulated cells from control cells, *not* cell types. They surface because, on the **unintegrated** data, the cell types with the strongest interferon response (monocytes most of all) split into separate CTRL and STIM clusters — so a "STIM monocyte" cluster's top markers are ISGs. Most clusters are still driven by cell-type identity, which is why cell-type markers dominate the lists overall. Tutorial 05 integrates CTRL and STIM so clusters track biology, not condition, and the ISGs drop out of the per-cluster marker lists.
2.  No. `FindAllMarkers` / `FindMarkers` treat each cell as an independent observation, which is **anti-conservative** for condition-level comparisons across biological samples. For STIM vs CTRL within a cell type, use **pseudobulk** + DESeq2/edgeR/limma — that's exactly what Tutorial 06 covers, on this dataset.

</details>
:::

## Step 3 — Visualizing markers

PBMC canonical markers — useful to keep handy:

| Cell type | Canonical markers |
|---|---|
| CD14+ Monocyte | `CD14`, `LYZ`, `S100A8`, `S100A9` |
| CD16+ Monocyte | `FCGR3A`, `MS4A7`, `CDKN1C` |
| Dendritic cell (DC) | `FCER1A`, `CST3`, `CLEC10A` |
| pDC | `LILRA4`, `IL3RA`, `CLEC4C` |
| CD4 Naive T | `IL7R`, `CCR7`, `LEF1` |
| CD4 Memory T | `IL7R`, `S100A4` |
| CD8 T | `CD8A`, `CD8B`, `GZMK` |
| NK | `GNLY`, `NKG7`, `KLRD1` |
| B | `MS4A1`, `CD79A`, `CD79B` |
| T-activated | `CD69`, `IL2RA` |
| Mk (megakaryocyte) | `PPBP`, `PF4` |
| Eryth | `HBA1`, `HBB` |

### 3a. Violin plots — expression distribution per cluster

```{r}
#| label: M3-viz_vln
#| fig-cap: "Violin plots of six canonical PBMC marker genes across clusters (resolution 0.5)."

p_vln <- VlnPlot(ifnb,
        features = c("CD14", "LYZ", "MS4A1", "CD79A", "CD8A", "GNLY"),
        pt.size = 0, ncol = 3) &
  xlab("Cluster (RNA_snn_res.0.5)") &
  ylab("Log-normalized expression") &
  theme(plot.title = element_text(size = 11))
p_vln <- p_vln +
  patchwork::plot_annotation(
    title    = "Canonical PBMC marker expression by cluster",
    subtitle = "Monocyte (CD14/LYZ), B (MS4A1/CD79A), CD8 T (CD8A) and NK (GNLY) markers",
    caption  = "Module 3 · Markers & Annotation"
  )
p_vln

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod3_C5_marker_violins.png"), p_vln,
       width = 10, height = 6, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Each panel is one marker gene; x-axis = cluster number, y-axis = log-normalized expression. The width of a violin at a given height represents how many cells in that cluster have that expression level, so a tall, narrow spike means most cells in the cluster are consistently expressing the gene. Look for clusters with a violin that rises visibly higher than all others for a given marker — for example, one or two clusters dominating `CD14` identify the monocyte clusters, and one or two dominating `MS4A1`/`CD79A` identify B cells. Clusters with flat (near-zero) violins for all six markers are likely T cells or NK cells not covered by this panel — check `CD8A` and `GNLY` for those.
:::

### 3b. Feature plots — expression on UMAP

```{r}
#| label: M3-viz_feature
#| fig-cap: "UMAP feature plots showing where four canonical markers (CD14, MS4A1, CD8A, GNLY) are expressed."

p_feature <- FeaturePlot(ifnb,
            features = c("CD14", "MS4A1", "CD8A", "GNLY"),
            ncol = 2) &
  xlab("UMAP 1") &
  ylab("UMAP 2") &
  theme(plot.title = element_text(size = 11))
p_feature <- p_feature +
  patchwork::plot_annotation(
    title    = "Marker expression overlaid on the UMAP",
    subtitle = "CD14 (monocytes), MS4A1 (B), CD8A (CD8 T), GNLY (NK)",
    caption  = "Module 3 · Markers & Annotation"
  )
p_feature

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod3_C6_marker_featureplots.png"), p_feature,
       width = 9, height = 8, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Each of the four panels places every cell at its UMAP coordinates, coloring it from light (low expression) to dark blue/purple (high expression). Healthy marker geography looks like a bright patch localized to one or a few islands — `CD14` should light up the monocyte region, `MS4A1` the B-cell island, `CD8A` the cytotoxic T region, and `GNLY` the NK region. Cells that appear bright in two non-adjacent islands may indicate doublets or a cross-reactive probe; diffuse high expression across the whole UMAP usually means the gene is broadly expressed and not a good discriminating marker for annotation.
:::

### 3c. Heatmap — top markers across clusters

```{r}
#| label: M3-viz_heatmap
#| fig-cap: "Heatmap of the top-10 markers per cluster (z-scored expression), cells grouped by cluster."

top10 <- ifnb.markers |> group_by(cluster) |> slice_max(avg_log2FC, n = 10)

# Tutorial 01 only scaled the HVGs; re-scale the top-marker genes so any
# non-HVG markers also appear in the heatmap rather than being dropped.
ifnb <- ScaleData(ifnb, features = unique(top10$gene))

p_heatmap <- DoHeatmap(ifnb, features = top10$gene) +
  NoLegend() +
  labs(
    title    = "Top-10 markers per cluster — scaled expression heatmap",
    subtitle = "Each column is a cell, grouped by cluster; rows are marker genes",
    caption  = "Module 3 · Markers & Annotation"
  )
p_heatmap

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod3_C7_marker_heatmap.png"), p_heatmap,
       width = 12, height = 10, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Columns are individual cells grouped by cluster (cluster boundaries are marked by thin colored bars on the x-axis); rows are the top-10 marker genes for each cluster. Color encodes z-scored expression: red = high relative expression, blue = low. A well-behaved marker block shows a bright red diagonal band — each set of marker rows lights up in exactly its own cluster and is blue elsewhere. Marker rows that stay red across many clusters are broadly expressed genes (less informative for annotation); rows that are consistently blue everywhere except one cluster are the most specific markers. Note that because columns represent cells rather than aggregated clusters, large clusters dominate the x-axis width, which can make small populations hard to see — the dot plot ([code](#lst-M3-viz_dot)) normalizes for cluster size.
:::

### 3d. Dot plot — % expressing × mean expression

[]{#lst-M3-viz_dot}

```{r}
#| label: M3-viz_dot
#| fig-cap: "Dot plot of a 16-gene canonical marker panel across clusters — dot size is percent expressing, colour is scaled mean expression."

canonical_markers <- c("CD14", "LYZ", "FCGR3A", "MS4A7", "FCER1A", "CST3",
                       "IL7R", "CCR7", "S100A4", "CD8A", "GZMK", "GNLY",
                       "NKG7", "MS4A1", "CD79A", "PPBP")
p_dot <- DotPlot(ifnb, features = canonical_markers) +
  RotatedAxis() +
  labs(
    title    = "Canonical PBMC marker panel across clusters",
    subtitle = "Dot size = % of cells expressing; colour = scaled average expression",
    x        = "Marker gene",
    y        = "Cluster (RNA_snn_res.0.5)",
    colour   = "Avg. expression",
    size     = "Percent expressed"
  )
p_dot

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod3_C8_marker_dotplot.png"), p_dot,
       width = 10, height = 6, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
The grid has clusters on the y-axis and marker genes on the x-axis. Each dot encodes two things simultaneously: **size** = the fraction of cells in that cluster expressing the gene (larger = more cells expressing it), and **color** = the scaled average expression level (darker = higher). To identify a cluster, look for the column(s) where it shows large, dark dots — those are its marker genes. A cluster with a big dark dot on `CD14` and `LYZ` is a monocyte; one showing `MS4A1` and `CD79A` is a B cell. Clusters with no large dots across the whole panel need additional markers (check `IL7R`, `CD3D`, or `CD3E` for T cells, which aren't in this panel). This is the primary visualization to use when deciding labels in [code](#lst-M3-rename) below.
:::

::: {.callout-important title="Think about it"}
1.  Which plot is best for **showing a 16-gene marker panel across 12 clusters on one page**?
2.  Which is best for **checking whether a single gene has spatial structure on the UMAP**?
3.  Which is most misleading if your cluster sizes are very unequal, and why?

<details>

<summary>Show answers</summary>

1.  **Dot plot** — compact, shows both percent-expressing (dot size) and mean expression (color) in a genes × clusters grid.
2.  **Feature plot** — overlays expression on the 2-D embedding so you can see geography.
3.  **Heatmap** can be misleading: it shows cells, not clusters, so a dominant cluster (say 40% of cells) visually monopolizes the plot. Dot plot normalizes per cluster and avoids this.

</details>
:::

## Step 4 — Manual annotation

Use the canonical marker table above and the dot plot ([code](#lst-M3-viz_dot)) to assign labels. Your cluster **count and numbering** depend on the random seed, Seurat version, and the fact that this is *unintegrated* data — so you'll often have **more clusters than cell types** (e.g. 15–20, with two clusters of the same type that split CTRL vs STIM). **Do not copy these labels blindly** — run `levels(ifnb)` to see your clusters, then label each one from the table that matches *your* dot plot (it's fine to give two clusters the same label).

[]{#lst-M3-rename}

```{r}
#| label: M3-rename
#| eval: false

# First, see how many clusters you actually have, and their ids:
levels(ifnb)
# On the UNINTEGRATED data, cell types often split by CONDITION (e.g. CTRL vs
# STIM monocytes), so you'll usually get MORE clusters than cell types — often
# ~15-20, sometimes two clusters of the same type. That's expected here;
# Tutorial 05 integrates the samples and collapses these.

# TEMPLATE — give each cluster id a label. Edit this named vector so the names on
# the LEFT (cluster ids, in quotes) cover YOUR clusters from the dot plot. The
# left side is the cluster id, the right side is the label you assign. You can
# reuse a label (e.g. two monocyte clusters both "CD14 Mono"). It's fine to label
# only some clusters while you iterate — unlisted ones keep their number.
new.cluster.ids <- c(
  "0"  = "CD14 Mono",
  "1"  = "CD4 Naive T",
  "2"  = "CD4 Memory T",
  "3"  = "B",
  "4"  = "CD8 T",
  "5"  = "NK",
  "6"  = "CD16 Mono",
  "7"  = "T activated",
  "8"  = "DC",
  "9"  = "Mk",
  "10" = "B activated",
  "11" = "pDC",
  "12" = "Eryth"
  # ... add entries for any remaining cluster ids your run produced (13, 14, ...)
)

# Apply only the labels whose cluster id actually exists. This keys on the vector
# NAMES (cluster ids), so it works for any number of clusters and never throws the
# "'names' attribute [N] must be the same length as the vector [M]" error.
new.cluster.ids <- new.cluster.ids[names(new.cluster.ids) %in% levels(ifnb)]
ifnb <- RenameIdents(ifnb, new.cluster.ids)
ifnb$celltype_manual <- Idents(ifnb)

DimPlot(ifnb, reduction = "umap", label = TRUE, repel = TRUE) + NoLegend()
```

::: {.callout-tip title="Reading the output"}
The UMAP now shows your cell-type labels instead of cluster numbers — each island should carry a readable name. A well-annotated map has compact, non-overlapping label positions (aided by `repel = TRUE`) and each island dominated by a single biological category. If two neighboring islands share the same label, that is expected on unintegrated data (the same cell type split by condition); if a large island has a label that doesn't feel right (e.g. "CD14 Mono" sitting next to the T-cell area), revisit the dot plot ([code](#lst-M3-viz_dot)) and check whether a different marker combination fits better.
:::

::: {.callout-important title="Think about it"}
1.  Two clusters share `CD3D` / `CD3E` as top markers. Are they the same cell type?
2.  What's the risk of assigning a label from a single marker gene?

<details>

<summary>Show answers</summary>

1.  Not necessarily — `CD3D/E` are pan-T markers. Both clusters are likely T cells, but could be CD4 Naive vs CD4 Memory vs CD8 vs T-activated. You need *discriminating* markers (`IL7R`+`CCR7` for naive, `S100A4` for memory, `CD8A` for cytotoxic, `IL2RA` for activated) on top of the shared T-cell markers.
2.  Many markers are non-specific, and biological "marker genes" are often graded rather than binary. Label from **combinations** (a marker *plus* absence of another marker) and always cross-check.

</details>
:::

## Step 5 — Validate against the ground-truth labels

The `ifnb` dataset ships with author-curated labels in `seurat_annotations`. This is rare in real life — but it makes a beautiful teaching moment to compare your manual labels to a published ground truth.

```{r}
#| label: M3-validate
#| fig-cap: "Side-by-side UMAPs comparing your manual cluster annotation against the author-curated ground-truth labels that ship with ifnb."

# Side-by-side UMAPs
p_manual <- DimPlot(ifnb, reduction = "umap", label = TRUE, repel = TRUE) +
  NoLegend() +
  labs(
    title  = "Manual annotation (yours)",
    x      = "UMAP 1",
    y      = "UMAP 2"
  )
p_truth  <- DimPlot(ifnb, reduction = "umap", group.by = "seurat_annotations",
                    label = TRUE, repel = TRUE) +
  NoLegend() +
  labs(
    title  = "Author labels (ground truth)",
    x      = "UMAP 1",
    y      = "UMAP 2"
  )

p_manual
p_truth

# Combined two-panel manual-vs-truth comparison, shared annotation
p_validate <- (p_manual | p_truth) +
  patchwork::plot_annotation(
    title    = "Manual annotation vs author ground truth on the same embedding",
    subtitle = "Cells should land in matching positions if your manual labels agree with the authors'",
    caption  = "Module 3 · Markers & Annotation"
  )

# Cross-tabulate. Each row should ideally concentrate on a single column.
manual_vs_truth <- table(manual = Idents(ifnb), truth = ifnb$seurat_annotations)
manual_vs_truth

# --- Figures out -----------------------------------------------------------
save_fig(file.path(out_dir, "Mod3_C10_umap_manual.png"),         p_manual,
       width = 7,  height = 6, dpi = 300)
save_fig(file.path(out_dir, "Mod3_C10_umap_ground_truth.png"),   p_truth,
       width = 7,  height = 6, dpi = 300)
save_fig(file.path(out_dir, "Mod3_C10_umap_manual_vs_truth.png"), p_validate,
       width = 13, height = 6, dpi = 300)

# --- Table out: cross-tab of manual labels vs ground-truth labels ----------
as.data.frame(manual_vs_truth) |>
  readr::write_csv(file.path(out_dir, "Mod3_C10_manual_vs_truth.csv"))
```

::: {.callout-tip title="Reading the output"}
The two UMAP panels side-by-side show your labels (left) and the author ground-truth (right) on the same embedding: identical islands in matching positions mean agreement; an island that carries a different name in each panel is a disagreement. The `manual_vs_truth` cross-tab printed below has your labels as rows and the author labels as columns — each cell is a count of cells assigned that label pair. A perfect annotator would see counts concentrated on a diagonal (each of your labels maps cleanly to one author label); off-diagonal counts identify the specific cell types where your calls and the authors' disagree.
:::

::: {.callout-important title="Think about it"}
1.  Where do you and the authors **agree**? Where do you **disagree**?
2.  Which cell type was easiest to call from markers alone? Hardest?
3.  In Tutorial 05, after integration, the same cluster numbers will mean different things. Why is `seurat_annotations` (a per-cell label) more durable than your `RNA_snn_res.0.5` manual labels (a per-cluster mapping)?

<details>

<summary>Show answers</summary>

1.  Big myeloid (CD14 / CD16 / DC) and lymphoid (B / NK) clusters usually agree perfectly. T-cell subdivisions (CD4 Naive vs CD4 Memory vs T activated) and rare populations (pDC, Eryth, Mk) are the common disagreements — they require multi-gene logic and have soft boundaries.
2.  **Easiest:** B cells (`MS4A1`+, `CD79A`+, otherwise quiet) and NK (`GNLY`+, `NKG7`+, no `CD3`). **Hardest:** the various CD4 T subsets, because they share most markers and differ in subtle, graded ways.
3.  Per-cell labels travel with the cell — if a cell is `CD14 Mono` in `seurat_annotations` here, it stays `CD14 Mono` after integration too. Cluster numbers are recomputed every time you re-cluster (and integration recomputes), so *cluster-level* labels are tied to a specific clustering and have to be re-applied.

</details>
:::

## Step 6 — Save the annotated object

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

saveRDS(ifnb, file = "../data/ifnb_annotated.rds")
```

## Continue to Tutorial 04 — Reference Annotation

You've named the clusters by hand. The immediate next step, **[Tutorial 04 — Reference-based Annotation](Tutorial_04_Reference_Annotation.html)**, cross-checks those manual labels against an automated reference (Azimuth) on this same annotated object — confirming the confident calls and flagging the cells the reference can't place.

After that, **[Tutorial 05 — Two-Sample Integration](Tutorial_05_Integration.html)** re-does the PCA + neighbor-graph + clustering steps on an **integrated** representation that aligns CTRL and STIM, and you'll watch what happens to:

-   the UMAP (CTRL/STIM islands collapse onto each other)
-   the clusters (cell types stop splitting in two)
-   the markers (ISGs stop being top markers; cell-type markers reclaim the top)

That integrated `.rds` is then the input to **pseudobulk DE** in Tutorial 06.

## Credits

-   Dataset: `ifnb` via Bioconductor's `muscData::Kang18_8vs8()` — Kang *et al.* (2017) *Nat Biotechnol*
-   PBMC canonical markers: [Seurat PBMC 3k tutorial](https://satijalab.org/seurat/articles/pbmc3k_tutorial.html), [Azimuth PBMC reference](https://azimuth.hubmapconsortium.org/references/human_pbmc/), [PanglaoDB](https://panglaodb.se/)
-   Inspiration: [scNotebooks](https://integrativebioinformatics.github.io/scNotebooks/)