Bio Applied Dimensionality Reduction

scRNA-seq: Dimensionality Reduction and Clustering

When to Use

• Reducing a highly variable gene (HVG) matrix to a PCA embedding before clustering or trajectory analysis
• Building a UMAP plot to visualize cell populations, batches, or conditions
• Clustering cells with Leiden/Louvain and deciding on a resolution
• Debugging "too many/too few clusters", "clusters don't match known cell types", or UMAP looking like a blob
• Finding marker genes per cluster (rank_genes_groups, FindAllMarkers) to annotate cell types

Version Compatibility

scanpy ≥1.10, anndata ≥0.10, umap-learn ≥0.5, leidenalg ≥0.10, scikit-learn ≥1.3, Python ≥3.10. Seurat ≥5.0, R ≥4.3.

Prerequisites

• pip install scanpy anndata umap-learn leidenalg python-igraph scikit-learn statsmodels
• Input: a QC'd, normalized (log1p, CPM/median-scaled), HVG-subset expression matrix (cells × genes). See bio-single-cell-preprocessing for how to get there.
• Concepts: what a k-NN graph is, why PCA needs scaled input, what "resolution" means in community detection.

Pipeline Overview

HVG matrix (cells × ~2000 genes)
  -> Scale (z-score, clip +/-10)
  -> PCA (top 10-50 PCs)
  -> k-NN graph (n_neighbors=15-20)
  -> UMAP (visualization only)
  -> Leiden clustering (on the graph, not the UMAP coords)

Goal: turn a cells x genes matrix into a low-dimensional embedding and discrete clusters. Approach: scale -> PCA for a robust distance metric -> k-NN graph in PCA space -> UMAP for 2D plotting -> Leiden for clustering on the graph -> rank_genes_groups for marker genes.

import numpy as np
import scanpy as sc

def embed_and_cluster(adata, n_pcs=20, n_neighbors=15, resolution=0.5, seed=42):
    """Standard PCA -> neighbors -> UMAP -> Leiden pipeline on an HVG-subset AnnData.

    Assumes adata.X is normalized/log1p'd and already subset to HVGs.
    Mutates and returns adata with obsm['X_pca'], obsm['X_umap'], obs['leiden'].
    """
    sc.pp.scale(adata, max_value=10)          # z-score + clip at +/-10 SD
    sc.tl.pca(adata, n_comps=n_pcs, svd_solver='arpack', random_state=seed)
    sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs, random_state=seed)
    sc.tl.umap(adata, min_dist=0.3, random_state=seed)
    sc.tl.leiden(adata, resolution=resolution, random_state=seed, flavor='igraph', n_iterations=2)
    return adata


def top_markers(adata, groupby='leiden', method='wilcoxon', n_genes=5):
    """Rank marker genes per cluster with FDR-corrected significance."""
    sc.tl.rank_genes_groups(adata, groupby=groupby, method=method)
    names = adata.uns['rank_genes_groups']['names']
    return {group: list(names[group][:n_genes]) for group in names.dtype.names}

Equivalent manual PCA (useful outside scanpy, e.g. plain AnnData/NumPy):

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.neighbors import NearestNeighbors

def pca_knn(X, n_pcs=20, n_neighbors=15, seed=42):
    """Scale -> PCA -> k-NN graph, returning PCA coords and neighbor indices/distances."""
    X_dense = X.toarray() if hasattr(X, 'toarray') else X
    X_scaled = np.clip(StandardScaler().fit_transform(X_dense), -10, 10)
    pca = PCA(n_components=n_pcs, random_state=seed)
    X_pca = pca.fit_transform(X_scaled)
    knn = NearestNeighbors(n_neighbors=n_neighbors, metric='euclidean', n_jobs=-1).fit(X_pca)
    distances, indices = knn.kneighbors(X_pca)
    return X_pca, distances, indices, pca.explained_variance_ratio_

Choosing n_pcs: plot cumulative variance explained; pick ~2x the "elbow" PC. Under-including merges cell types; over-including adds noise but rarely destroys structure. Typical: 10-20 PCs for PBMC, 30-50 for complex tissue.

n_neighbors effect:

• Small (5-10): emphasizes local structure, tight clusters, may fragment continuous populations
• Large (30-50): emphasizes global structure, smoother but may merge distinct types
• Default 15: good starting point for 5k-50k cells

UMAP vs t-SNE:

Resolution tuning: start at 0.5; increase to 0.8-1.0 if known cell types are merged; decrease to 0.3 if over-fragmented. Validate against marker genes, not just cluster count.

Parameter Decision Table

Seurat (R) Equivalent

library(Seurat)

seu <- ScaleData(seu, features = VariableFeatures(seu), vars.to.regress = NULL)
seu <- RunPCA(seu, npcs = 30, verbose = FALSE)
seu <- FindNeighbors(seu, dims = 1:20, k.param = 15)
seu <- RunUMAP(seu, dims = 1:20, min.dist = 0.3)
seu <- FindClusters(seu, resolution = 0.5, algorithm = 4)  # algorithm 4 = Leiden

markers <- FindAllMarkers(seu, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top5 <- markers |> dplyr::group_by(cluster) |> dplyr::slice_max(avg_log2FC, n = 5)

Pitfalls

• UMAP distances are not transcriptional distances: two clusters far apart in UMAP are not necessarily more different. Use PCA space or expression values for quantitative comparisons.
• Scaling required before PCA: high-mean housekeeping genes dominate PCs without scaling, burying rare-marker variation.
• Clip at +/-10 SD: doublets that survived QC can have extreme z-scores (50+) that distort PCs. Clipping at 10 is standard (sc.pp.scale(max_value=10)).
• Leiden requires leidenalg/python-igraph: pip install leidenalg python-igraph. Louvain can produce internally disconnected communities — prefer Leiden.
• Cluster on the graph, not the UMAP coordinates: sc.tl.leiden operates on the neighbor graph in adata.uns['neighbors']; clustering directly on 2D UMAP coordinates throws away information and is not reproducible across UMAP re-runs.
• Resolution is not portable: resolution=0.5 on one dataset does not give the same number of clusters on another; always re-tune per dataset.
• Batch effects appear as clusters: run batch correction (Harmony, scVI, BBKNN) before clustering if samples come from different batches — see bio-single-cell-batch-integration.
• Multiple testing: apply FDR correction (Benjamini-Hochberg) for marker gene testing; sc.tl.rank_genes_groups reports pvals_adj — use that column, not raw pvals.

See Also

• bio-single-cell-preprocessing — QC, normalization, and HVG selection feeding into this pipeline
• bio-single-cell-clustering — deeper Leiden/Louvain tuning and cluster validation
• bio-single-cell-batch-integration — Harmony/scVI correction before embedding multi-batch data
• bio-single-cell-markers-annotation — annotating clusters into cell types from marker genes