Tutorial 05 — scRNA-seq: Two-Sample Integration

Align CTRL and STIM with Harmony and Seurat anchors; compare before/after; check for over-correction

Author

Single Cell RNA-seq Workshop

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:

Note

Companion lecture: Lecture 05 — Multi-Sample Integration · Companion book chapter: Chapter 5 — Multi-Sample Integration — 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

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:

Harmony — fast, soft-clustering–based, runs in PC space, near-default in modern Seurat workflows
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.

How to use this page

Download the .qmd source: Tutorial_05_Integration.qmd. 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.
Make sure you have completed Tutorial 03 — it writes ifnb_annotated.rds.
Work through the chunks, flipping eval: true as you go.

Note

Companion notebook. scNotebooks Module 05 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

Common errors / things that bite

harmony is not installed — it’s installed in the Tutorial 00 setup; 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.

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 for the build profile

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), ")"))
}

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.

# 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

Think about it

Why pick HVGs per sample instead of on the joint matrix?
Suppose CTRL has 11,000 cells and STIM has 12,000. Will the per-sample HVG lists overlap?

Show answers

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

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.

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

Think about it

We pass reduction = "harmony" to FindNeighbors and RunUMAP. Why not reduction = "pca"?
What’s group.by.vars = "stim" doing in RunHarmony?

Show answers

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

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.

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

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.

Think about it

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?
Why does DefaultAssay(ifnb_s) <- "integrated" matter for clustering but not for marker testing?

Show answers

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

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

The six-panel comparison (code) shows each integration method colored by sample (top row) and by cell type (bottom row).

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

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.

Think about it

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?
In the by-cell-type row, each label should occupy a single region. After integration, has the doubling-up from Tutorial 02 collapsed?
Do Harmony and Seurat anchors produce visually similar UMAPs? Where do they differ?

Show answers

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

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 scores each cell on a nine-gene IFN-β module and checks that the STIM–CTRL gap survives integration.

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

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.

Think about it

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

Show answers

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

Step 4 — Re-find markers on integrated clusters

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

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

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.

Think about it

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?
Are the cell-type markers (CD14, MS4A1, CD8A, GNLY, …) cleaner / more specific to single clusters now?

Show answers

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

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

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

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

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

How to load and QC a multi-sample scRNA-seq dataset (Tut 01)
How to recognize a batch effect from PCs / clusters / UMAP (Tut 02)
How to identify cell types from canonical markers (Tut 03)
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 and Tutorial 10 — the analysis pipeline on Talapas — 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 (Korsunsky et al. 2019, Nat Methods); Seurat anchors (Stuart, Butler et al. 2019, Cell)
Reference vignette: Seurat integration tutorial
Inspiration: scNotebooks