---
title: "Tutorial 02 — scRNA-seq: Dimensionality Reduction & Clustering"
subtitle: "PCA, how many PCs, neighbor graph, graph clustering at multiple resolutions, UMAP — and seeing the batch effect"
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 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:


::: callout-note
**Companion lecture:** [Lecture 02 — Dimensionality Reduction & Clustering](../Lecture_Folder/Lecture_02_DimReduction_Clustering.html) · **Companion book chapter:** [Chapter 2 — Dimensionality Reduction & Clustering](../Resources_Folder/Chapter_02_DimReduction.html) — 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`](Tutorial_10_Talapas_Pipeline.html)
:::


-   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

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

::: {.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_02_DimReduction_Clustering.qmd.txt" download="Tutorial_02_DimReduction_Clustering.qmd">Tutorial_02_DimReduction_Clustering.qmd</a>. 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](Tutorial_01_QC_Preprocessing.html)** — 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


::: {.callout-warning title="Common 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](Tutorial_00_Setup_RStudio_Packages.html)) 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()`.
:::


::: {.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_02_DimReduction_Clustering.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: M2-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), ")"))
}
```

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

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

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

## Step 1 — Linear dimensionality reduction (PCA)

```{r}
#| label: M2-pca
#| fig-cap: "PC1 gene-loading heatmap on 500 cells — the axis is dominated by interferon-stimulated genes (the batch/condition signal)."

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

::: {.callout-tip title="Reading 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).
:::

::: {.callout-important title="Think 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?

<details>

<summary>Show answers</summary>

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.

</details>
:::

## Step 2 — Choose how many PCs to keep

```{r}
#| label: M2-elbow
#| fig-cap: "Elbow plot of per-PC standard deviation — used to pick a defensible number of PCs (the elbow sits around PC15–20)."

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

::: {.callout-tip title="Reading 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.
:::

::: {.callout-important title="Think 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*?

<details>

<summary>Show answers</summary>

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.

</details>
:::

```{r}
#| label: M2-n_pcs

n_pcs <- 20
```

## Step 3 — Build the neighbor graph

```{r}
#| label: M2-neighbors

ifnb <- FindNeighbors(ifnb, dims = 1:n_pcs)
```

::: {.callout-important title="Think about it"}
1.  What graph does `FindNeighbors()` build?
2.  Why operate on PC space rather than raw expression?

<details>

<summary>Show answers</summary>

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

</details>
:::

## Step 4 — Cluster at several resolutions

```{r}
#| label: M2-cluster

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

head(ifnb@meta.data)
```

::: {.callout-tip title="Reading 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.
:::

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

::: {.callout-important title="Think 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?

<details>

<summary>Show answers</summary>

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.

</details>
:::

## Step 5 — UMAP

[]{#lst-M2-umap}

```{r}
#| label: M2-umap
#| fig-cap: "UMAP of the unintegrated ifnb cells coloured by graph cluster (resolution 0.5)."

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

::: {.callout-tip title="Reading 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](#lst-M2-umap) is colored by cluster. Now color the *same UMAP* by **sample** and by **author cell-type label** and compare.

```{r}
#| label: M2-batch_diagnostic
#| fig-cap: "The same unintegrated UMAP coloured by sample (left/top) and by author cell-type label — revealing per-sample islands driven by the batch effect."

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

::: {.callout-tip title="Reading 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.
:::

::: {.callout-important title="Think 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?

<details>

<summary>Show answers</summary>

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.

</details>
:::

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

```{r}
#| label: M2-crosstab

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

::: {.callout-tip title="Reading 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.
:::

::: {.callout-important title="Think 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?

<details>

<summary>Show answers</summary>

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.

</details>
:::

## Save and continue to Tutorial 03

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

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

::: callout-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

-   Dataset: `ifnb` via Bioconductor's `muscData::Kang18_8vs8()` — Kang *et al.* (2017) *Nat Biotechnol* 36: 89–94
-   Reference vignettes: [Seurat integration tutorial](https://satijalab.org/seurat/articles/integration_introduction.html) · [clustree](https://cran.r-project.org/package=clustree)
-   Inspiration: [scNotebooks](https://integrativebioinformatics.github.io/scNotebooks/)