---
title: "Tutorial 05 — scRNA-seq: Two-Sample Integration"
subtitle: "Align CTRL and STIM with Harmony and Seurat anchors; compare before/after; check for over-correction"
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

Tutorial 05 in the laptop-friendly `ifnb` series. So far we've:


::: callout-note
**Companion lecture:** [Lecture 05 — Multi-Sample Integration](../Lecture_Folder/Lecture_05_Integration.html) · **Companion book chapter:** [Chapter 5 — Multi-Sample Integration](../Resources_Folder/Chapter_05_Integration.html) — the long-form prose treatment of this tutorial's material, with cross-references to the prerequisite appendices. · **HPC version (Talapas):** [Talapas analysis pipeline — `05_integration.R`](Tutorial_10_Talapas_Pipeline.html)
:::


-   Loaded a two-sample dataset and run QC (Tutorial 01)
-   Run PCA + clustering on the unintegrated joint matrix and **seen the batch effect** — cell types splitting into per-sample islands (Tutorial 02)
-   Annotated those unintegrated clusters with canonical markers (Tutorial 03)

Now we **integrate**. Integration is the step that aligns shared cell-type structure across samples while preserving real biological differences. Done right, the same cell type in CTRL and STIM ends up in the *same* cluster instead of two parallel ones.

We'll run **two** common integration methods and compare:

1.  **Harmony** — fast, soft-clustering–based, runs in PC space, near-default in modern Seurat workflows
2.  **Seurat anchors (CCA)** — Seurat's original integration; slower, but useful for moderate-batch / few-sample datasets

We'll then re-cluster, re-UMAP, and **diagnose over-correction** by overlaying the cell-type ground truth and the IFN-β response signal.

::: {.callout-tip title="How to use this page"}
1.  Download the `.qmd` source: <a href="Tutorial_05_Integration.qmd.txt" download="Tutorial_05_Integration.qmd">Tutorial_05_Integration.qmd</a>. If your browser saves the file as `Tutorial_05_Integration.qmd.txt`, **drop the trailing `.txt`** so the filename ends in `.qmd`, then open it in RStudio.
2.  Make sure you have completed **[Tutorial 03](Tutorial_03_Markers_Annotation.html)** — it writes `ifnb_annotated.rds`.
3.  Work through the chunks, flipping `eval: true` as you go.
:::

::: callout-note
**Companion notebook.** [scNotebooks Module 05](https://integrativebioinformatics.github.io/scNotebooks/modules/en/module05/module05.html) runs the same Seurat-anchor + Harmony comparison on the same `ifnb` dataset, in a Jupyter notebook you can open directly in Colab. If you want a second walkthrough — or want the workflow in Spanish or Portuguese — that's a good place to look. They also compute a quantitative neighborhood-mixing metric we don't show here.
:::

## Dataset — resuming from Tutorial 03

| File | Produced by | What it is |
|---|---|---|
| `../data/ifnb_annotated.rds` | Tutorial 03 | Seurat object with normalized counts, scaled HVGs, unintegrated PCA / UMAP / clusters, manual labels, and ground-truth `seurat_annotations`. |

## Learning goals

-   Split a multi-sample Seurat object into per-sample objects and run normalization independently per sample
-   Run **Harmony** as a fast batch-correction in PC space
-   Run **Seurat CCA anchor-based** integration as an alternative
-   Re-cluster on the integrated reduction and compare to the unintegrated clustering
-   **Diagnose over- and under-correction** with side-by-side UMAPs colored by sample, by cell type, and by ISG score
-   Save an integrated, annotated object ready for pseudobulk DE in Tutorial 06


::: {.callout-warning title="Common errors / things that bite"}
**`harmony` is not installed** — it's installed in the [Tutorial 00 setup](Tutorial_00_Setup_RStudio_Packages.html); Method A here needs it (the Seurat-anchor path does not).

**`DefaultAssay()` is `"integrated"` when running `FindMarkers`** — the `IntegrateData` step switches the default assay. Always `DefaultAssay(seu) <- "RNA"` *before* any DE / marker work. The integrated assay holds *batch-corrected* values that are valid for *geometry* (PCA, UMAP, clustering) but NOT for differential expression.

**ISG score gap looks zero after integration** — you scored on the integrated assay instead of the `RNA` assay. `AddModuleScore` uses `DefaultAssay(seu)`. Switch back to `RNA` first.

**Harmony's `RunHarmony` errors with "object has no `pca` reduction"** — you skipped the `RunPCA` step. Harmony works on PC space, not on the raw scaled counts.

**`FindAllMarkers` on the Harmony object returns "No DE genes identified" + "data layers are not joined"** — `ifnb_h` is built with `merge()`, which in Seurat v5 leaves the RNA assay with per-sample split layers (`counts.CTRL`/`counts.STIM`, `data.*`). `FindAllMarkers` then silently skips every test and returns an empty table, so the downstream `group_by(cluster)` errors with *"Column `cluster` is not found"*. Fix: `ifnb_h <- JoinLayers(ifnb_h)` after the merge (the Harmony chunk now does this). `Layers(ifnb_h)` should show single `counts`/`data` layers.

**`FindIntegrationAnchors` seems to hang at "Running CCA"** — it isn't hung; CCA is an SVD-based step that grows with cell number, and this `ifnb` is ~23k cells, so it can take **several minutes** on a laptop (much slower than Harmony's seconds). As long as a CPU core is busy, let it finish — "Finding anchors" and "Finding integration vectors" follow. The `Different features in new layer data than already exists for scale.data` warnings during "Scaling features" are benign (Seurat v5 re-scales the integration features over the HVG `scale.data` from Tutorial 01).

**Harmless chatter during the Harmony chunk** — three things you can ignore: (1) `Different features in new layer data than already exists for scale.data` from `ScaleData` (it re-scales on the integration features, replacing Tutorial 01's HVG `scale.data` — exactly as intended); (2) lines like `4036exists, retrying for cluster 91 …` are **Harmony2**'s internal soft-clustering retrying a centroid collision (leaks through even with `verbose = FALSE`); (3) `The default method for RunUMAP has changed … UWOT … cosine` is a one-time informational note. The run worked if you see `Number of communities: …` and `Optimization finished`.
:::


::: {.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_05_Integration.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: M5-setup

library(Seurat)
library(tidyverse)
library(harmony)
library(patchwork)

set.seed(2026)

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

ifnb <- readRDS("../data/ifnb_annotated.rds")
ifnb
```

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

## Step 1 — Split by sample, then normalize each sample independently

Both Harmony and Seurat anchors work best when each sample's normalization, HVGs, and PCA are computed **per sample first**, then combined. The reasoning: HVGs are about within-sample variance, and per-sample HVGs avoid letting one dominant sample dictate the gene set.

```{r}
#| label: M5-split_norm

# Split object by sample of origin (the `stim` column from Tut 01)
ifnb.list <- SplitObject(ifnb, split.by = "stim")

# Per-sample LogNormalize + HVG selection
ifnb.list <- lapply(ifnb.list, function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
  x
})
ifnb.list
```

::: {.callout-important title="Think about it"}
1.  Why pick HVGs **per sample** instead of on the joint matrix?
2.  Suppose CTRL has 11,000 cells and STIM has 12,000. Will the per-sample HVG lists overlap?

<details>

<summary>Show answers</summary>

1.  HVGs measure variance *across cells*. If one sample is twice as big, joint-HVG selection over-weights its variance pattern. Per-sample HVGs let each sample contribute equally to the integration anchors.
2.  Strongly — many top HVGs (cell-type markers like `MS4A1`, `LYZ`, `GNLY`) appear in both samples' top HVGs because they vary within both. The interesting ones are the HVGs unique to each sample (e.g. ISGs in STIM but not CTRL) — `SelectIntegrationFeatures()` will keep features that are HVG in *most* samples.

</details>
:::

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

## Method A — Harmony

Harmony runs **after** PCA on the joint matrix. It iteratively soft-clusters cells in PC space and shifts each cluster centroid so that batches mix within clusters — producing a corrected `harmony` reduction that is a drop-in replacement for `pca` in `FindNeighbors` and `RunUMAP`.

```{r}
#| label: M5-harmony

# Re-merge the per-sample objects, run joint normalization / HVG / scale / PCA
features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 2000)

# A single re-merged object is fine for Harmony — it just needs PCA + a batch column.
ifnb_h <- merge(ifnb.list[[1]], y = ifnb.list[[2]],
                merge.data = TRUE,
                project    = "ifnb_harmony")
# In Seurat v5, merge() leaves the RNA assay with per-sample split layers
# (counts.CTRL/counts.STIM, data.*). Join them back into single counts/data layers
# so later RNA-assay steps (ISG module score, FindAllMarkers) work. Harmony aligns
# in PCA space, so joining the expression layers doesn't affect the alignment.
ifnb_h <- JoinLayers(ifnb_h)
VariableFeatures(ifnb_h) <- features
ifnb_h <- ScaleData(ifnb_h, features = features, verbose = FALSE)
ifnb_h <- RunPCA(ifnb_h, features = features, npcs = 30, verbose = FALSE)

# Harmony alignment in PC space — `group.by.vars` is the batch column
ifnb_h <- RunHarmony(ifnb_h, group.by.vars = "stim", reduction = "pca",
                     reduction.save = "harmony", verbose = FALSE)

# Re-cluster + UMAP on the harmony reduction
ifnb_h <- FindNeighbors(ifnb_h, reduction = "harmony", dims = 1:20)
ifnb_h <- FindClusters(ifnb_h, resolution = 0.5)
ifnb_h <- RunUMAP(ifnb_h, reduction = "harmony", dims = 1:20,
                  reduction.name = "umap_harmony")
```

::: {.callout-important title="Think about it"}
1.  We pass `reduction = "harmony"` to `FindNeighbors` and `RunUMAP`. Why not `reduction = "pca"`?
2.  What's `group.by.vars = "stim"` doing in `RunHarmony`?

<details>

<summary>Show answers</summary>

1.  Because the whole point is that the **harmony reduction** is the batch-corrected version of PCA. Passing `pca` here would re-do the analysis on the uncorrected coordinates and recreate the per-sample islands.
2.  Telling Harmony which metadata column encodes the batch / sample axis. Harmony will iteratively shift centroids so that *each Harmony cluster* contains a roughly equal proportion of CTRL and STIM cells — i.e. it's correcting along the `stim` axis.

</details>
:::

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

## Method B — Seurat CCA anchor-based integration

The classical Seurat approach: find **anchors** (mutual nearest neighbors after CCA) between sample pairs, then learn a transformation that maps each sample into a shared space.

```{r}
#| label: M5-seurat_anchors

ifnb.anchors <- FindIntegrationAnchors(object.list = ifnb.list,
                                       anchor.features = features,
                                       reduction = "cca",
                                       dims = 1:20)
ifnb_s <- IntegrateData(anchorset = ifnb.anchors, dims = 1:20)

# After IntegrateData the default assay becomes "integrated".
# Run the rest of the workflow on it.
DefaultAssay(ifnb_s) <- "integrated"
ifnb_s <- ScaleData(ifnb_s, verbose = FALSE)
ifnb_s <- RunPCA(ifnb_s, npcs = 30, verbose = FALSE)
ifnb_s <- FindNeighbors(ifnb_s, reduction = "pca", dims = 1:20)
ifnb_s <- FindClusters(ifnb_s, resolution = 0.5)
ifnb_s <- RunUMAP(ifnb_s, reduction = "pca", dims = 1:20,
                  reduction.name = "umap_seurat")
```

::: callout-note
After `IntegrateData()` the **default assay** becomes `"integrated"`. That assay holds the corrected expression values used for clustering / UMAP — but for **marker-gene tests** you should switch back to the `"RNA"` assay, because integration's corrected values aren't the right input to `FindMarkers()`. This is a common gotcha.
:::

::: {.callout-important title="Think about it"}
1.  Anchor-based integration is widely considered the most rigorous of the alignment methods, but it's also slow on large datasets. For `ifnb` (\~23k cells, 2 samples), how does runtime compare to Harmony?
2.  Why does `DefaultAssay(ifnb_s) <- "integrated"` matter for clustering but not for marker testing?

<details>

<summary>Show answers</summary>

1.  Harmony finishes in **well under a minute**; Seurat CCA anchors are much slower — on this ~23k-cell `ifnb` the `FindIntegrationAnchors` "Running CCA" step alone can take **several minutes** on a laptop (it's an SVD-based step that grows quickly with cell number), so don't be alarmed if it sits there. For atlas-scale data (1M cells, 50 samples) the difference balloons to hours-vs-minutes, and Harmony or scVI become the practical choice.
2.  The integrated assay holds **batch-corrected** expression values that are valid for *geometry* (PCA, neighbor graph, UMAP). They're **not** valid for differential expression because the correction can dampen or distort real biological variation. Always run `FindMarkers()` on the original `RNA` assay.

</details>
:::

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

## Step 2 — Compare unintegrated, Harmony, and Seurat-anchor UMAPs

The six-panel comparison ([code](#lst-M5-compare)) shows each integration method colored by sample (top row) and by cell type (bottom row).

[]{#lst-M5-compare}

```{r}
#| label: M5-compare
#| fig-cap: "Unintegrated vs Harmony vs Seurat-anchor UMAPs, coloured by sample (mixing check) and by cell type (structure check)."

# Unintegrated UMAP from Tut 02 (still inside the `ifnb` object)
p_un_sample   <- DimPlot(ifnb, reduction = "umap", group.by = "stim") +
  ggtitle("Unintegrated — by sample") + labs(x = "UMAP 1", y = "UMAP 2", colour = "Sample")
p_un_cell     <- DimPlot(ifnb, reduction = "umap", group.by = "seurat_annotations",
                         label = TRUE, repel = TRUE) +
  NoLegend() + ggtitle("Unintegrated — by cell type") + labs(x = "UMAP 1", y = "UMAP 2")

# Harmony
p_h_sample    <- DimPlot(ifnb_h, reduction = "umap_harmony", group.by = "stim") +
  ggtitle("Harmony — by sample") + labs(x = "UMAP 1", y = "UMAP 2", colour = "Sample")
p_h_cell      <- DimPlot(ifnb_h, reduction = "umap_harmony", group.by = "seurat_annotations",
                         label = TRUE, repel = TRUE) +
  NoLegend() + ggtitle("Harmony — by cell type") + labs(x = "UMAP 1", y = "UMAP 2")

# Seurat anchors
p_s_sample    <- DimPlot(ifnb_s, reduction = "umap_seurat", group.by = "stim") +
  ggtitle("Seurat anchors — by sample") + labs(x = "UMAP 1", y = "UMAP 2", colour = "Sample")
p_s_cell      <- DimPlot(ifnb_s, reduction = "umap_seurat", group.by = "seurat_annotations",
                         label = TRUE, repel = TRUE) +
  NoLegend() + ggtitle("Seurat anchors — by cell type") + labs(x = "UMAP 1", y = "UMAP 2")

# Render whichever pair is most useful — or use patchwork to show all six
p_compare_sample <- (p_un_sample + p_h_sample + p_s_sample) +
  patchwork::plot_annotation(
    title    = "Sample mixing across integration methods",
    subtitle = "CTRL vs STIM should interleave within clusters after integration",
    caption  = "Module 5 · Two-Sample Integration"
  )
p_compare_sample

p_compare_cell <- (p_un_cell + p_h_cell + p_s_cell) +
  patchwork::plot_annotation(
    title    = "Cell-type structure across integration methods",
    subtitle = "Each cell type should collapse from two per-sample islands into one region",
    caption  = "Module 5 · Two-Sample Integration"
  )
p_compare_cell

# --- Figures out -----------------------------------------------------------
save_fig(file.path(out_dir, "Mod5_C6_compare_by_sample.png"), p_compare_sample,
       width = 15, height = 5, dpi = 300)
save_fig(file.path(out_dir, "Mod5_C6_compare_by_celltype.png"), p_compare_cell,
       width = 15, height = 5, dpi = 300)
```

::: {.callout-tip title="Reading the output"}
Two three-panel figures appear side-by-side. The **by-sample** row (x = UMAP 1, y = UMAP 2, colour = CTRL/STIM) is your **mixing diagnostic**: in the unintegrated panel the two colours form parallel islands per cell type; after successful integration they should interleave within each cluster. The **by-cell-type** row uses labelled cluster colours — in the unintegrated panel each cell type is split into two sample-specific blobs; after integration each type should collapse into one. Comparing Harmony and Seurat-anchor columns lets you see whether the two methods differ in how aggressively they mix samples or how cleanly they separate rare populations.
:::

::: {.callout-important title="Think about it"}
1.  In the **by-sample** row, after integration the two colors should be **interleaved** within each cluster. Are they? Are CTRL and STIM mixed at the same ratio in every cluster, or do some clusters still skew?
2.  In the **by-cell-type** row, each label should occupy a **single** region. After integration, has the doubling-up from Tutorial 02 collapsed?
3.  Do Harmony and Seurat anchors produce visually similar UMAPs? Where do they differ?

<details>

<summary>Show answers</summary>

1.  After good integration the colors interleave nicely. A cluster that's still \>80% one sample after integration is suspicious — either (a) integration didn't fully correct, or (b) the cells in that cluster genuinely only exist in one condition (e.g. a STIM-only activated state).
2.  Yes — most cell types collapse from two CTRL/STIM islands into one. The exception is `T activated` and `B activated`, which exist primarily in STIM and may stay STIM-skewed even after integration. **That's not over-correction — that's real biology.**
3.  Harmony tends to mix samples slightly more aggressively (it's tuning toward batch-balance); Seurat anchors are slightly more conservative on rare populations. Pick based on your priorities: Harmony when sample mixing matters most, anchors when preserving rare cell types matters most.

</details>
:::

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

## Step 3 — Diagnose over-correction

The risk of integration is **erasing real biology**. The ISG response is real — STIM cells *should* differ from CTRL cells in every cell type — and integration should **mix the cells in space** without zeroing out their expression. The code in [code](#lst-M5-isg_score) scores each cell on a nine-gene IFN-β module and checks that the STIM–CTRL gap survives integration.

[]{#lst-M5-isg_score}

```{r}
#| label: M5-isg_score
#| fig-cap: "Per-cell interferon-stimulated-gene (ISG) module score on the Harmony and Seurat-anchor UMAPs — integration should mix cells in space without flattening the ISG response."

isg_set <- c("ISG15","IFI6","IFIT1","IFIT3","ISG20","MX1","OAS1","RSAD2","IFIT2")

# Score each cell on its IFN-β response — using the original RNA assay
DefaultAssay(ifnb_h) <- "RNA"
ifnb_h <- AddModuleScore(ifnb_h, features = list(isg_set), name = "ISG_score")

DefaultAssay(ifnb_s) <- "RNA"
ifnb_s <- AddModuleScore(ifnb_s, features = list(isg_set), name = "ISG_score")

# Score on the unintegrated object too, for a reference point
DefaultAssay(ifnb)   <- "RNA"
ifnb   <- AddModuleScore(ifnb,   features = list(isg_set), name = "ISG_score")

p_isg_h <- FeaturePlot(ifnb_h, features = "ISG_score1", reduction = "umap_harmony") +
  ggtitle("Harmony UMAP, ISG score (per-cell)") +
  labs(x = "UMAP 1", y = "UMAP 2", colour = "ISG score")
p_isg_s <- FeaturePlot(ifnb_s, features = "ISG_score1", reduction = "umap_seurat") +
  ggtitle("Seurat-anchors UMAP, ISG score (per-cell)") +
  labs(x = "UMAP 1", y = "UMAP 2", colour = "ISG score")

p_isg <- (p_isg_h + p_isg_s) +
  patchwork::plot_annotation(
    title    = "Per-cell IFN-β response (ISG module score) after integration",
    subtitle = "Scored on the RNA assay; high-ISG STIM cells should co-cluster with low-ISG CTRL cells",
    caption  = "Module 5 · Two-Sample Integration"
  )
p_isg

# Check that mean ISG score per cell type is still ~0 in CTRL and ~positive in STIM
# (i.e. the biology survived integration)
isg_by_celltype <- ifnb_h@meta.data |>
  group_by(seurat_annotations, stim) |>
  summarise(ISG_mean = mean(ISG_score1), .groups = "drop") |>
  pivot_wider(names_from = stim, values_from = ISG_mean) |>
  arrange(desc(STIM - CTRL))
isg_by_celltype

# --- Outputs: ISG FeaturePlots + per-cell-type STIM-vs-CTRL ISG means -------
save_fig(file.path(out_dir, "Mod5_C7_isg_featureplots.png"), p_isg,
       width = 12, height = 5, dpi = 300)

readr::write_csv(isg_by_celltype,
                 file.path(out_dir, "Mod5_C7_isg_mean_by_celltype.csv"))
```

::: {.callout-tip title="Reading the output"}
Two **FeaturePlots** appear side-by-side (Harmony left, Seurat-anchors right); each cell is coloured by its ISG module score (low = dark/cool, high = bright/warm). A successful integration shows cells of the same type co-localised on UMAP regardless of score, but STIM cells still carrying higher scores than their CTRL neighbours. The **`isg_by_celltype` table** below prints the mean ISG score per cell-type × condition pair, sorted by STIM–CTRL gap; monocytes typically show the largest gap (> 0.5 in module-score units) while T cells show a smaller but still positive gap. A gap of zero in any row is the key over-correction warning sign.
:::

::: {.callout-important title="Think about it"}
1.  In the FeaturePlot, do high-ISG cells cluster with low-ISG cells of the same type? They should — that's exactly what integration is supposed to do (put them in the same cluster but at different positions in *gene-expression space*, not erase the difference).
2.  Look at the per-cell-type CTRL-vs-STIM mean ISG score table. The STIM column should be **noticeably higher** than CTRL for every cell type. If a row's STIM-CTRL gap is near zero, integration has flattened that cell type's IFN response — over-correction.

<details>

<summary>Show answers</summary>

1.  In a healthy integration, low-ISG (CTRL) and high-ISG (STIM) cells of the same type **co-cluster on UMAP** (good — sample mixing) but the STIM cells are still **higher on the ISG score** (good — biology preserved). It's the spatial mixing that integration changes; the per-cell expression vectors remain real.
2.  CD14 Mono and CD16 Mono are the strongest IFN responders and should show the largest STIM-CTRL gap (often \> 0.5 in module-score units). T/B/NK gaps will be smaller but still positive. If you see a gap of 0 anywhere, check that you scored on the `RNA` assay rather than the integrated assay.

</details>
:::

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

## Step 4 — Re-find markers on integrated clusters

The top markers per integrated cluster ([code](#lst-M5-integrated_markers)) let you check whether cell-type signals are now cleaner than in the unintegrated Tutorial 03 results.

[]{#lst-M5-integrated_markers}

```{r}
#| label: M5-integrated_markers

# Switch back to the RNA assay for marker testing
DefaultAssay(ifnb_h) <- "RNA"

# ifnb_h was built with merge() (in the Harmony step), and Seurat v5 merge leaves
# the RNA assay with per-sample split layers (counts.CTRL/counts.STIM, data.*).
# FindAllMarkers needs a single joined layer, so rejoin before testing — otherwise
# every test fails with "data layers are not joined" and returns no markers.
if (length(Layers(ifnb_h, search = "data")) > 1) ifnb_h <- JoinLayers(ifnb_h)

Idents(ifnb_h) <- "seurat_clusters"

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

ifnb_h_top_markers <- ifnb_h_markers |>
  group_by(cluster) |>
  slice_max(avg_log2FC, n = 5) |>
  ungroup()
ifnb_h_top_markers

# --- Tables out: all integrated-cluster markers + top-5 per cluster ---------
readr::write_csv(ifnb_h_markers,
                 file.path(out_dir, "Mod5_C8_integrated_markers.csv"))
readr::write_csv(ifnb_h_top_markers,
                 file.path(out_dir, "Mod5_C8_integrated_top5_markers.csv"))
```

::: {.callout-tip title="Reading the output"}
`ifnb_h_top_markers` prints as a tibble with columns `cluster`, `gene`, `avg_log2FC`, `pct.1`, `pct.2`, and `p_val_adj` — five rows per integrated cluster, sorted by highest average log2 fold change. A clean integration result shows canonical PBMC markers (`LYZ`/`CD14` for monocytes, `MS4A1` for B cells, `GNLY`/`NKG7` for NK, `CD8A` for cytotoxic T) at the top of their respective clusters, with ISGs (`ISG15`, `IFIT1`, …) absent or much lower ranked than they were in the unintegrated Tutorial 03 markers.
:::

::: {.callout-important title="Think about it"}
1.  Compare the top markers per cluster here to the top markers from Tutorial 03 (the unintegrated version). Did ISGs disappear from the top of the lists?
2.  Are the cell-type markers (`CD14`, `MS4A1`, `CD8A`, `GNLY`, …) cleaner / more specific to single clusters now?

<details>

<summary>Show answers</summary>

1.  Yes — ISGs no longer dominate top markers because the *clusters* now correspond to cell types, not to STIM-vs-CTRL. ISGs still exist in the data (and we can find them in DE between STIM and CTRL within each cell type — that's Tutorial 06's job) but they don't drive cluster identity anymore.
2.  Yes. `CD14` should top exactly one cluster (the merged CD14 Mono cluster) instead of being split between two unintegrated CTRL-CD14 / STIM-CD14 clusters. Same for B / T / NK markers.

</details>
:::

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

## Step 5 — Save the integrated, annotated object

For Tutorial 06's pseudobulk DE we want an object with:

-   The original `RNA` assay as the default (for any DE work)
-   A clean cluster column from the integrated clustering
-   `seurat_annotations` carried through

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

# Make sure RNA is the default assay before saving
DefaultAssay(ifnb_h) <- "RNA"

saveRDS(ifnb_h, file = "../data/ifnb_integrated.rds")
```

::: callout-tip
Tutorial 06 (Pseudobulk DE) loads `ifnb_integrated.rds` and tests `STIM` vs `CTRL` **within each cell type** using DESeq2 on cell-type × sample × donor pseudobulk counts.
:::

## What you've now learned end-to-end

1.  How to load and QC a multi-sample scRNA-seq dataset (Tut 01)
2.  How to recognize a batch effect from PCs / clusters / UMAP (Tut 02)
3.  How to identify cell types from canonical markers (Tut 03)
4.  How to integrate two samples with two different methods, compare them, and check for over-correction (Tut 05 — this one)

The same workflow scales to atlas-scale data — but at \~23k cells, a laptop is fine. When you graduate to bigger datasets, see the **Talapas HPC modules** — [Tutorial 09 — SLURM Basics](Tutorial_09_SLURM_Basics.html) and [Tutorial 10 — the analysis pipeline on Talapas](Tutorial_10_Talapas_Pipeline.html) — for running this same pipeline on the cluster.

## Credits

-   Dataset: `ifnb` via Bioconductor's `muscData::Kang18_8vs8()` — Kang *et al.* (2017) *Nat Biotechnol*
-   Methods: [Harmony](https://portals.broadinstitute.org/harmony/) (Korsunsky *et al.* 2019, *Nat Methods*); Seurat anchors (Stuart, Butler *et al.* 2019, *Cell*)
-   Reference vignette: [Seurat integration tutorial](https://satijalab.org/seurat/articles/integration_introduction.html)
-   Inspiration: [scNotebooks](https://integrativebioinformatics.github.io/scNotebooks/)