# Bio Applied Dimensionality Reduction

> Compute PCA/UMAP embeddings and Leiden clusters for scRNA-seq with scanpy/Seurat; tune n_pcs/n_neighbors/resolution, find markers via rank_genes_groups. Use for UMAP plots, clustering single-cell data, or picking PCs/resolution.

- Skill: `pavel-kravchenko/bio-applied-dimensionality-reduction` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-dimensionality-reduction`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-dimensionality-reduction/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

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

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

```python
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:**

| | UMAP | t-SNE |
|--|------|-------|
| Speed | Fast (minutes) | Slow (hours for >50k cells) |
| Global structure | Partially preserved | Not preserved |
| Between-cluster distances | Approximate (directionally meaningful) | Meaningless |
| Reproducibility | Seed-controlled | Less stable |

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

| Parameter | Default | Increase when | Decrease when |
|-----------|---------|--------------|---------------|
| `n_pcs` | 15-20 | Complex tissue, many cell types | Simple dataset, elbow is early |
| `n_neighbors` | 15 | Need global structure | Need fine local detail |
| `min_dist` | 0.3-0.5 | Clusters too compressed | Want tighter visual clusters |
| Leiden `resolution` | 0.5 | Known subtypes being merged | Too many spurious clusters |

## Seurat (R) Equivalent

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