Tutorial 02 — scRNA-seq: Dimensionality Reduction & Clustering

PCA, how many PCs, neighbor graph, graph clustering at multiple resolutions, UMAP — and seeing the batch effect

Author

Single Cell RNA-seq Workshop

NoteRunning 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 two of the series. We start from the preprocessed object saved in Tutorial 01 (the ifnb dataset — CTRL + STIM PBMCs joined together) and walk through:

Note

Companion lecture: Lecture 02 — Dimensionality Reduction & Clustering · Companion book chapter: Chapter 2 — Dimensionality Reduction & Clustering — the long-form prose treatment of this tutorial’s material, with cross-references to the prerequisite appendices. · HPC version (Talapas): Talapas analysis pipeline — 02_dimreduction_clustering.R

  • Linear dimensionality reduction with PCA
  • Choosing how many PCs to keep
  • Building the shared nearest-neighbor graph
  • Clustering at multiple resolutions and comparing them
  • UMAP embedding — and an explicit look at the batch effect between samples
Important

We deliberately do NOT run integration in this tutorial. Running PCA + clustering + UMAP on the unintegrated joint matrix lets you see what a batch effect looks like — cell types splitting into per-sample islands instead of co-clustering. Tutorial 05 will fix it. Until then, every UMAP we produce is a teaching artifact, not a final figure.

TipHow to use this page

The rendered HTML shows the code but does not execute it. To run it:

  1. Download the .qmd source: Tutorial_02_DimReduction_Clustering.qmd. If your browser saves the file as Tutorial_02_DimReduction_Clustering.qmd.txt, drop the trailing .txt so the filename ends in .qmd, then open it in RStudio.
  2. Make sure you have completed Tutorial 01 — it writes the ifnb_preprocessed.rds that this tutorial loads.
  3. Work through the chunks, flipping eval: true as you go.

Dataset — resuming from Tutorial 01

File Produced by What it is
../data/ifnb_preprocessed.rds Tutorial 01 A Seurat object with raw + log-normalized counts for both samples (stim ∈ {"CTRL","STIM"}), 2,000 HVGs, scaled data, and the seurat_annotations ground-truth labels.

Data dictionary — the Seurat object at this stage

Slot Class Contents
seu@assays$RNA@counts dgCMatrix Raw UMI counts (cells that survived QC only).
seu@assays$RNA@data dgCMatrix log1p((counts / total_per_cell) * 10000) values from NormalizeData().
seu@assays$RNA@scale.data matrix z-scored expression from ScaleData(). Input to RunPCA().
seu@assays$RNA@var.features character ~2,000 highly variable gene symbols.
seu@meta.data$stim factor "CTRL" or "STIM". The covariate we expect to see drive PC1/PC2.
seu@meta.data$seurat_annotations factor Author-curated cell-type labels.

Learning goals

  • Run RunPCA() and read PCA diagnostic plots
  • Use ElbowPlot() to pick a defensible number of PCs
  • Run FindNeighbors() → FindClusters() at several resolutions and inspect the difference
  • Switch the active identity between clusterings with Idents()
  • Generate a UMAP and overlay cluster labels
  • Diagnose a batch effect by coloring the same UMAP by sample, by cluster, and by cell-type label
WarningCommon errors / things that bite

readRDS("../data/ifnb_preprocessed.rds") fails — Tutorial 01 writes it to your project’s data/ folder (../data/ relative to scripts/), and this tutorial reads it back from there. So (a) run Tutorial 01 first, and (b) keep all the tutorial .qmd files together in scripts/ so they share one project data/. Run getwd() to confirm R is in your scripts/ folder — the .. paths resolve one level up from there — and check that ../data/ifnb_preprocessed.rds exists.

FindClusters is slow — presto (installed in the Tutorial 00 setup) gives a much faster Wilcoxon implementation; Seurat detects it automatically.

DimPlot shows the wrong cluster column — Idents() defaults to seurat_clusters, which is the most recent FindClusters() resolution. To plot a different resolution, use DimPlot(seu, group.by = "RNA_snn_res.0.3") instead of relying on Idents().

TipSolutions / 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_02_DimReduction_Clustering.html. See Exercise_Folder/_quarto-solutions.yml for the build profile

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  Mod2_C<chunk>_<name>
# into ../output/Mod2/ so it can be cross-referenced from the rest of the site.
#   Mod2 = Module 2 (this tutorial);  C<n> = the nth code chunk.
# ---------------------------------------------------------------------------
out_dir <- "../output/Mod2"
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_preprocessed.rds")
ifnb

Step 1 — Linear dimensionality reduction (PCA)

ifnb <- RunPCA(
  ifnb,
  features = VariableFeatures(object = ifnb)
)

print(ifnb[["pca"]], dims = 1:5, nfeatures = 5)

# --- Table out: top +/- gene loadings for the first five PCs ---------------
pca_loadings <- Loadings(ifnb[["pca"]])[, 1:5]
top_pc_loadings <- purrr::map_dfr(1:5, function(d) {
  ord <- order(pca_loadings[, d])
  tibble::tibble(
    PC        = paste0("PC", d),
    direction = c(rep("negative", 5), rep("positive", 5)),
    gene      = c(rownames(pca_loadings)[head(ord, 5)],
                  rownames(pca_loadings)[tail(ord, 5)]),
    loading   = c(pca_loadings[head(ord, 5), d],
                  pca_loadings[tail(ord, 5), d])
  )
})
readr::write_csv(top_pc_loadings, file.path(out_dir, "Mod2_C3_pca_loadings.csv"))

# PC1 loadings heatmap (DimHeatmap draws to the active device); add a title and
# capture it as a patchwork element so it can be annotated and saved.
p_pca_heatmap <- DimHeatmap(ifnb, dims = 1, cells = 500, balanced = TRUE,
                            fast = FALSE) +
  ggtitle("PC1 loadings — top genes across 500 cells")
p_pca_heatmap <- patchwork::wrap_elements(p_pca_heatmap) +
  patchwork::plot_annotation(
    title    = "PCA diagnostic — PC1 gene-loading heatmap",
    subtitle = "Cells ordered by PC1 score; expect interferon-stimulated genes to dominate",
    caption  = "Module 2 · Dimensionality Reduction & Clustering"
  )
p_pca_heatmap

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod2_C3_pca_heatmap.png"), p_pca_heatmap,
       width = 7, height = 6, dpi = 300)
TipReading the output

The print() call shows the top 5 positive and negative gene loadings for PCs 1–5: these are the genes that push cells to the high or low end of each PC. In the ifnb data, PC1’s positive loadings are dominated by interferon-stimulated genes (ISG15, IFI6, IFIT1, etc.), telling you that the IFN-β condition effect — not cell-type identity — explains the most variance in the unintegrated data. The heatmap rows are genes and columns are the 500 cells with the most extreme PC1 scores (ordered by score); a clear split into two expression blocks on PC1 is the visual signature of the batch/condition axis. Later PCs should start to show cell-type-specific patterns (e.g. monocyte genes on one PC, T-cell genes on another).

ImportantThink about it
  1. Look at the top genes loading PC1. Do you recognize any interferon-stimulated genes (ISG15, IFI6, IFIT1, IFIT3, ISG20, MX1, OAS1)? What does that suggest about what PC1 is encoding?
  2. PCA here runs on HVGs, not all genes. What would change if we included all genes?
Show answers
  1. PC1 in the unintegrated ifnb object is heavily loaded with ISGs because the IFN-β response is the largest source of variance in this dataset. PC1 is essentially stim (CTRL vs STIM) — i.e. the batch / condition axis, not a cell-type axis. That’s the batch effect we’ll fix in Tutorial 05.
  2. PCs would be dominated by highly-expressed housekeeping genes whose variance simply reflects sequencing depth. HVG-restricted PCA isolates biologically-variable signal.

Step 2 — Choose how many PCs to keep

p_elbow <- ElbowPlot(ifnb, ndims = 40) +
  labs(
    title    = "PCA elbow plot — variance captured per principal component",
    subtitle = "The elbow (~PC15–20) marks where additional PCs stop adding signal",
    x        = "Principal component",
    y        = "Standard deviation"
  )
p_elbow

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod2_C4_elbow.png"), p_elbow,
       width = 7, height = 5, dpi = 300)
TipReading the output

The x-axis is the principal component number (1 to 40) and the y-axis is its standard deviation (how much variance it explains). A healthy elbow plot starts high on the left, drops steeply for the first few PCs, then flattens into a near-horizontal plateau. The “elbow” — the kink where the curve stops falling noticeably — is your cutoff: PCs before it carry signal, PCs after it carry mostly noise. For ifnb, the elbow sits around PC 15–20; choosing anywhere from 15 to 30 gives similar downstream results, so err on the generous side.

ImportantThink about it
  1. Where is “the elbow” on an ElbowPlot, and why does it suggest a cutoff?
  2. Is it risky to choose too few PCs? Too many?
Show answers
  1. The elbow is where the variance explained per PC flattens out — additional PCs stop adding much signal. Picking the elbow captures most of the biological structure without the noise-heavy tail.
  2. Too few loses cell-type resolution (rare populations blur together). Too many pulls in noisy dimensions and inflates distances in the neighbor graph. Being slightly generous (e.g. 15–30) is usually safe; for ifnb the elbow sits around PC15–PC20.
n_pcs <- 20

Step 3 — Build the neighbor graph

ifnb <- FindNeighbors(ifnb, dims = 1:n_pcs)
ImportantThink about it
  1. What graph does FindNeighbors() build?
  2. Why operate on PC space rather than raw expression?
Show answers
  1. It builds a k-nearest-neighbor graph (default k = 20) and then a shared-nearest-neighbor graph: two cells are connected if they share many nearest neighbors. The Jaccard-weighted SNN graph is what clustering runs on.
  2. PC space is denoised, low-dimensional, and has well-behaved Euclidean distance. Raw expression is sparse, high-dimensional, and distance-distorted by the curse of dimensionality.

Step 4 — Cluster at several resolutions

ifnb <- FindClusters(
  ifnb,
  resolution = c(0.1, 0.3, 0.5, 0.7, 1.0)
)

head(ifnb@meta.data)
TipReading the output

The head() call shows the first few rows of meta.data; the key columns to notice are RNA_snn_res.0.1, RNA_snn_res.0.3, RNA_snn_res.0.5, RNA_snn_res.0.7, and RNA_snn_res.1 — one integer cluster assignment per cell for each resolution you ran. Lower-resolution columns will show fewer unique values (coarser clusters); higher-resolution columns will show more. The seurat_clusters column always tracks the most recent resolution called, so it currently equals RNA_snn_res.1. Use Idents(ifnb) <- "RNA_snn_res.0.5" to switch the active identity before plotting.

Note

FindClusters() here uses Seurat’s default Louvain algorithm (algorithm = 1). Lecture 02 introduces Leiden (algorithm = 4) as the modern preferred method; we use the Louvain default for portability because Leiden needs the leidenalg Python package via reticulate. To switch, add algorithm = 4 to the call above once that dependency is installed.

ImportantThink about it
  1. What does resolution control intuitively?
  2. Each resolution creates its own RNA_snn_res.XYZ column in the metadata. How do you make one of them the active identity?
Show answers
  1. The granularity knob for the graph-based community detection. Higher resolution ⇒ smaller, more numerous clusters.
  2. Idents(ifnb) <- "RNA_snn_res.0.5" — or any of the resolution columns. DimPlot(..., group.by = "RNA_snn_res.0.5") also works without changing active idents.

Step 5 — UMAP

ifnb <- RunUMAP(ifnb, dims = 1:n_pcs)

Idents(ifnb) <- "RNA_snn_res.0.5"
p_umap_clusters <- DimPlot(ifnb, reduction = "umap", label = TRUE) +
  labs(
    title    = "Unintegrated UMAP — graph clusters (res = 0.5)",
    subtitle = "Clusters computed on the joint CTRL + STIM matrix without integration",
    x        = "UMAP 1",
    y        = "UMAP 2",
    colour   = "Cluster"
  )
p_umap_clusters

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod2_C8_umap_clusters.png"), p_umap_clusters,
       width = 7, height = 6, dpi = 300)
TipReading the output

Each point is a cell projected into 2-D by UMAP; nearby points share similar transcriptomes in PC space. Colors represent the Louvain clusters at resolution 0.5 — each island or peninsula should correspond to a cell type or subpopulation. Well-separated, compact blobs indicate a clean cluster; ragged or poorly-separated shapes suggest either that the resolution is too high or that those populations are transcriptionally similar. Because this is the unintegrated data, you may already notice that some clusters appear in two nearby sub-clouds — a preview of the CTRL/STIM split you will quantify in Step 6.

Step 6 — Diagnose the batch effect

This is the key exercise. The UMAP in code is colored by cluster. Now color the same UMAP by sample and by author cell-type label and compare.

p_by_sample <- DimPlot(ifnb, reduction = "umap", group.by = "stim") +
  labs(
    title    = "Unintegrated UMAP — coloured by sample (stim)",
    subtitle = "CTRL and STIM cells form largely parallel islands",
    x        = "UMAP 1",
    y        = "UMAP 2",
    colour   = "Sample"
  )

p_by_celltype <- DimPlot(ifnb, reduction = "umap", group.by = "seurat_annotations",
                         label = TRUE, repel = TRUE) +
  NoLegend() +
  labs(
    title    = "Unintegrated UMAP — coloured by author cell-type",
    subtitle = "Most cell types split into two sister clusters, tracking stim",
    x        = "UMAP 1",
    y        = "UMAP 2"
  )

p_by_sample
p_by_celltype

# Combined two-panel batch-effect diagnostic, shared annotation
p_batch <- (p_by_sample | p_by_celltype) +
  patchwork::plot_annotation(
    title    = "Batch-effect diagnostic — sample vs cell-type on the same embedding",
    subtitle = "If cell types co-clustered, the two panels would overlay; here they don't",
    caption  = "Module 2 · Dimensionality Reduction & Clustering"
  )

# --- Figures out -----------------------------------------------------------
save_fig(file.path(out_dir, "Mod2_C9_umap_by_sample.png"),   p_by_sample,
       width = 7, height = 6, dpi = 300)
save_fig(file.path(out_dir, "Mod2_C9_umap_by_celltype.png"), p_by_celltype,
       width = 7, height = 6, dpi = 300)
save_fig(file.path(out_dir, "Mod2_C9_umap_batch_diagnostic.png"), p_batch,
       width = 13, height = 6, dpi = 300)
TipReading the output

The left panel colors each cell by sample (CTRL vs STIM): if the two colors are well-mixed across every island the data is batch-free; if they form largely separate regions or parallel sub-clouds, a condition/batch effect is dominating the embedding. The right panel uses the author cell-type labels: each label should ideally form one coherent blob, but on unintegrated data you’ll typically see each cell type split into two sister clouds — one CTRL-colored, one STIM-colored — because PC1 separates by condition rather than cell type. Seeing that split here is the intended outcome: it motivates integration in Tutorial 05.

ImportantThink about it
  1. In the by-sample plot, do the CTRL and STIM cells overlap, or do they form parallel islands?
  2. In the by-cell-type plot, does each cell type form one cluster or several? If a cell type splits, how does its split correlate with sample?
  3. Cross-tabulate: table(ifnb$RNA_snn_res.0.5, ifnb$stim). Are any clusters mostly CTRL? Mostly STIM?
  4. Based on these three diagnostics, do you need to integrate?
Show answers
  1. CTRL and STIM cells form largely parallel structures. Many cell-type clusters appear twice — once for CTRL, once for STIM — because the IFN-β response shifts every cell along PC1.
  2. Most cell types split into two sister clusters, with the split tracking stim. T cells appear once on the CTRL side and once on the STIM side; same for monocytes, B cells, NK cells, etc. This is the classic “per-sample islands” pattern — the batch / condition axis dominates the embedding.
  3. Most clusters at res = 0.5 will be >80% one sample. Some are >95% one sample.
  4. Yes. The batch axis is dominating PC1 and is splitting biologically identical cell types across two clusters. Tutorial 05 will integrate the data with Harmony and Seurat anchors, and you’ll see the same cells re-arrange into shared, sample-mixed clusters.

Step 7 — Cross-tab clusters vs cell types vs samples

# Clusters × samples — how mixed is each cluster across CTRL / STIM?
cluster_by_sample <- prop.table(
  table(cluster = ifnb$RNA_snn_res.0.5, sample = ifnb$stim), margin = 1) |>
  round(2)
cluster_by_sample

# Clusters × cell types — how clean is each cluster's biological identity?
cluster_by_celltype <- table(cluster = ifnb$RNA_snn_res.0.5,
                             celltype = ifnb$seurat_annotations)
cluster_by_celltype

# --- Tables out: cluster composition by sample and by cell type ------------
as.data.frame(cluster_by_sample) |>
  readr::write_csv(file.path(out_dir, "Mod2_C10_cluster_by_sample.csv"))

as.data.frame(cluster_by_celltype) |>
  readr::write_csv(file.path(out_dir, "Mod2_C10_cluster_by_celltype.csv"))
TipReading the output

cluster_by_sample is a clusters × samples proportion table (rows sum to 1): each number is the fraction of that cluster’s cells coming from CTRL or STIM. A value near 1.0 in one column means that cluster is almost entirely one sample — the hallmark of an unintegrated, condition-split cluster. cluster_by_celltype is a raw count table: ideally each row (cluster) would have counts concentrated in a single column (cell type), indicating a pure cluster; spread across multiple columns flags a mixed or ambiguous cluster.

ImportantThink about it
  1. A “good” cluster on integrated data would be ~50/50 CTRL/STIM (because both samples have all cell types) and dominated by a single cell-type label. Are any of your unintegrated clusters like that already? Why might some be more sample-mixed than others?
  2. What’s the lowest-abundance cell-type label in seurat_annotations? Does it split across two clusters or stay together?
Show answers
  1. A few clusters — typically very small populations like dendritic cells (DC) or megakaryocytes (Mk) — may already mix CTRL/STIM, because their absolute interferon response is small enough that other biology dominates their distance to PC1. Most clusters won’t.
  2. Often DC, Mk, or eryth (<2% each). Whether they split tells you something about how strong the IFN response is in those cells; rare cell types with weak ISG responses sometimes survive the unintegrated UMAP intact, while big subsets like CD14 monocytes split cleanly along the CTRL/STIM axis.

Save and continue to Tutorial 03

saveRDS(ifnb, file = "../data/ifnb_clustered.rds")
Tip

Continue with Tutorial 03 — Markers & Cell Type Annotation. We’ll do annotation on the unintegrated clusters first (it’s still possible — markers are cell-type-driven even when the embedding is sample-driven), then in Tutorial 05 we’ll redo PCA + clustering on the integrated representation and watch the embedding reorganize.

Credits