Dimensionality reduction — the intuition

Dimensionality reduction projects high-dimensional data into 2-3D

• Cells live in a very high-dimensional space (thousands of genes)
• Dimensionality reduction projects them into 2–50 dimensions while preserving the main structure
• Different projections emphasize different axes — the full signal isn’t visible in any single 2D view

Why PCA before everything else?

• The HVG-scaled matrix has ~2,000 genes × tens of thousands of cells. Most pairs of cells are nearly equidistant in that space — the curse of dimensionality at work.
• PCA does three jobs simultaneously:
  • Compresses the data: 2,000 genes → ~50 PCs that explain most of the structured variance
  • Denoises by discarding low-variance directions dominated by technical noise
  • Decorrelates the axes — PCs are orthogonal, so distances in PC space behave well for the neighbour-graph and clustering steps that follow
• Everything downstream — Harmony, FindNeighbors, clustering, UMAP — operates on the PC matrix, not the gene matrix

PCA finds axes of maximum variance

PCA finds axes of maximum variance

• The PCs are ordered latent variables
  • PC1 = direction in gene space along which the data spread the most.
  • PC2 = next-most variance, orthogonal to PC1.
  • PC3 perpendicular to both, etc.
• Same algebraic object two ways: eigendecomposition of the covariance matrix, or SVD of the centered data matrix \(\tilde X = U \Sigma V^T\)
  • Right singular vectors \(V\) → loadings (PC directions in gene space)
  • \(U \Sigma\) → scores (each cell’s coordinates in PC space)
  • \(\sigma_i^2 / \sum \sigma_j^2\) → proportion of variance explained by PC \(i\)
• RunPCA(), sc.tl.pca(), and prcomp()
  • all call SVD under the hood
  • it scales better than eigendecomposition when the gene count exceeds the cell count

In scRNA-seq what does each PC capture?

• PC1–PC10 typically separate major cell-type axes (T vs B vs myeloid; epithelial vs stromal vs immune)
• PC10–PC30 resolve subtypes and states (CD4 vs CD8, naïve vs memory, cycling vs resting)
• PC30–PC50 capture finer biology + a growing share of noise — keep them only if your downstream clustering still benefits
• Beyond PC50 is overwhelmingly noise

Center, then scale — and why this matters

• In Seurat, ScaleData() z-scores each gene across cells (mean 0, SD 1) before RunPCA(). In scanpy, sc.pp.scale() does the same.
• Without scaling, a few high-mean / high-variance genes (MALAT1, ribosomal proteins, some lncRNAs) dominate every PC simply because they have larger raw variance — not because they are biologically informative.
• This is why the canonical order HVG → ScaleData → RunPCA is fixed: HVGs select genes that carry biology, ScaleData puts them on equal footing, then PCA finds the structure.

Reading the scree plot

Reading the scree plot

• The elbow is where each additional PC stops adding meaningful variance — usually between PC15 and PC30 for 10x runs
• Cumulative variance is a sanity check, not a goal: a “90% rule” is misleading because the noise floor is high in scRNA-seq
• ElbowPlot(seu, ndims = 50) is the visual; JackStraw() is a permutation-based statistical alternative when you want a defensible cutoff

Loadings — which genes drive each PC?

Loadings — which genes drive each PC?

• Loadings are the weights that mix each gene into a PC. Big positive (or negative) loadings = genes moving cells along that axis.
• In the previous cartoon:
  • PC1 separates T-cell markers (CD3D, CD3E, TRAC) from B-cell markers (CD79A, MS4A1);
  • PC2 separates myeloid markers (LYZ, S100A8, S100A9) from lymphoid.
• Workflow tools:
  • VizDimLoadings(seu, dims = 1:4) for bar plots,
  • DimHeatmap(seu, dims = 1:9, balanced = TRUE) for a per-cell heatmap of the top loaders.
• This is the step where PCA earns its reputation as interpretable: t-SNE and UMAP cannot tell you which genes drive a separation — PCA can.

SNN graph → clusters (Louvain / Leiden)

• Build a shared nearest-neighbor (SNN) graph from the integrated PCs — each cell connects to its \(k\) nearest neighbours, weighted by how many of those neighbours are shared
• Cluster the graph at a chosen resolution — with Seurat’s default Louvain, or the modern Leiden (preferred, but needs a Python dependency)
• The graph (not the gene matrix nor UMAP) is the object the clustering actually runs on

KNN — the building block

• Each cell’s identity is set by the cells it is closest to in PC space
• The KNN graph is the substrate the SNN refinement and Leiden/Louvain modularity all run on top of

Why graph-based clustering instead of K-means?

• K-means assumes clusters are roughly spherical and equally sized — a poor match for scRNA-seq, where rare cell types coexist with abundant ones and clusters have irregular shapes in PC space
• K-means also requires you to know k in advance; for scRNA-seq the whole point is often discovering how many populations there are
• Graph-based methods (Louvain, Leiden) discover clusters as densely-connected communities, naturally handle uneven cluster sizes, and let resolution control granularity instead of \(k\)

Graph-based clustering — step by step

• Start with the KNN graph over cells in PCA space
• Find communities → initial partition; refine merges/splits based on modularity
• Aggregate the network, refine again → final partition

Choosing a resolution

• A lower resolution means fewer clusters, higher means more
• No single “right” resolution exists — the question is what cell types you want to resolve
• Typical starting range: 0.4 – 1.2 for Louvain/Leiden on a 10–30k cell dataset
• Use clustree to visualize cluster stability across a sweep of resolutions
• Cross-check by inspecting markers of merged vs. split clusters

From PCs to a 2D picture for human consumption

• PCA gave us a denoised ~30-D space; clustering already happened on the SNN graph in that space. What’s left is a picture — a 2D embedding so we can show clusters, gradients, and trajectories on a slide.

• Two non-linear methods dominate:

  • t-SNE (van der Maaten & Hinton 2008) and
  • UMAP (McInnes, Healy & Melville 2018)

• Both share the same philosophy:

  • In high-D: “for each cell, who are its neighbours, and how similar?”
  • In low-D: “place points in 2D so those neighbour relationships are preserved”

• Crucially, t-SNE/UMAP are visualizations, not analyses. They don’t replace PCA or clustering — they sit on top of them.

t-SNE: model neighbours as probabilities

t-SNE: model neighbours as probabilities

• For each cell \(i\), t-SNE turns distances to all other cells into a probability distribution \(p_{j|i}\) — Gaussian in high-D, with bandwidth \(\sigma_i\) chosen per cell so the effective neighbourhood size matches the user’s perplexity
• In the 2D embedding, the same idea gives \(q_{j|i}\) from the current 2D positions
• t-SNE then moves the 2D points to minimize the Kullback–Leibler divergence

\[D_{KL}(P \| Q) = \sum_{i \neq j} p_{ij} \log \frac{p_{ij}}{q_{ij}}\]

• KL divergence is asymmetric — and that’s the point

• Heavily penalised: \(p_{ij}\) large but \(q_{ij}\) small → “you put true neighbours far apart” → strong attractive gradient pulls them together

• Weakly penalised: \(p_{ij}\) small but \(q_{ij}\) large → “you put unrelated cells nearby” → near-zero contribution to the loss

UMAP

UMAP

• For each cell
  • UMAP builds a local weighted graph out to its \(k\) nearest neighbours, but normalises distances per cell
  • so points in sparse and dense regions of the data are still connected to their own local neighbours
• The local graphs are merged into one global fuzzy weighted graph (\(w_{ij} = a + b - a \cdot b\))
• The 2D embedding starts from a spectral initialization (or random) and is optimised by stochastic gradient descent

UMAP’s cross-entropy preserves more global structure

UMAP’s cross-entropy preserves more global structure

\[CE = -\sum_{i,j} \left[ w_{ij} \log v_{ij} + (1 - w_{ij}) \log(1 - v_{ij}) \right]\]

• Attraction term \(w_{ij} \log v_{ij}\): connected cells should stay close
• Repulsion term \((1 - w_{ij}) \log(1 - v_{ij})\): unconnected cells should stay apart
• Unlike KL, both errors are penalised → cluster positions and approximate inter-cluster distances are more faithful than t-SNE’s

n_neighbors and min_dist, in practice

(Seurat’s RunUMAP defaults; the underlying Python umap-learn uses different defaults of 15 and 0.1, so don’t be surprised if other tutorials quote those.)

• Small n_neighbors (5–10): very local detail; clusters fragment more. Large (30–50): broader topology; subtypes may merge.
• Small min_dist (≤ 0.1): tight clusters with clear gaps — good for discrete cell types. Large (0.4–0.8): spread-out — good for continuous gradients (development, pseudotime).
• UMAP is noticeably more stable across n_neighbors than t-SNE is across perplexity — usually one well-chosen UMAP is enough; with t-SNE you really do want to sweep.

Reading a UMAP responsibly

• What a UMAP is good for:
  • showing the cluster solution
  • mapping a known marker onto the cells (FeaturePlot)
  • spotting batch effects (color by sample after integration)
• What a UMAP is not good for
  • differential testing
  • distance-based statistics
  • conclusions about lineage relationships

What to remember from Lecture 02

• PCA on HVGs is the substrate every downstream step builds on — read the scree, read the loadings before you trust the embedding
• Build the SNN graph on the PCs, then cluster with Louvain/Leiden — graph-based methods beat K-means for scRNA-seq
• t-SNE preserves local neighbourhoods via KL divergence + the t-distribution; cluster distances in t-SNE are not meaningful
• UMAP preserves more global topology via cross-entropy on a fuzzy graph; today it’s the default
• Both are visualizations, not analyses — they sit on top of PCA + clustering, never replace them