--- name: bio-single-cell-clustering description: Dimensionality reduction and graph-based clustering for single-cell RNA-seq with Scanpy (Python) and Seurat (R). Resolves which algorithm to use (Leiden vs Louvain), how many PCs and neighbors to set, how to sweep and validate resolution, when a split is over-clustering, and why post-clustering marker p-values are not valid inference. Use when clustering cells, choosing a clustering resolution, deciding whether two clusters are one population, building a UMAP/tSNE, or judging whether clusters are real. tool_type: mixed primary_tool: Seurat --- ## Version Compatibility Reference examples tested with: scanpy 1.10+, Seurat 5.0+, anndata 0.10+ Before using code patterns, verify installed versions match. If versions differ: - Python: `pip show ` then `help(module.function)` to check signatures - R: `packageVersion('')` then `?function_name` to verify parameters If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying. # Single-Cell Clustering **"Cluster my cells"** -> Build a k-nearest-neighbor graph in PCA space, partition it into communities, and embed in 2D for display. - Python: `sc.pp.neighbors` + `sc.tl.leiden` + `sc.tl.umap` (Scanpy) - R: `FindNeighbors` + `FindClusters` + `RunUMAP` (Seurat) ## Governing Principle Resolution is not a truth knob, and clusters are not discoveries. Graph-based clustering partitions a kNN graph built in PCA space; there is no ground truth, no "correct" number of clusters, and resolution selects a scale of description rather than revealing one. Clusters are hypotheses that must be validated by markers, stability, and (where claims are made) a significance test before any partition is named a cell population. Over-clustering is the default failure mode: any homogeneous blob can be bisected, and a higher resolution will always split it further. UMAP/tSNE distances, cluster sizes, and apparent gaps are artifacts of a non-linear neighbor-preserving objective and are not metric data (Chari & Pachter 2023) - cluster on the graph, never on embedding coordinates. The deepest trap: clustering chooses labels to maximize between-group separation, so running a marker test on those same clusters tests a hypothesis built from the data used to test it (double-dipping), and the resulting p-values are not merely inflated, they are invalid inference. ## Leiden vs Louvain Leiden is the current default for graph community detection because Louvain can return internally disconnected communities (Traag 2019). Seurat still ships Louvain as its default algorithm; Scanpy uses Leiden but is mid-migration between backends. Methodology evolves - verify the current default and backend against the installed package docs before pinning a pipeline. | Method | Model / assumption | Use when | Fails when | |--------|-------------------|----------|------------| | Leiden | Modularity/CPM optimization with a refinement phase guaranteeing connected communities | Default for scRNA-seq; reproducibility matters; large graphs (faster) | Backend/iteration count left unpinned -> silently different labels across Scanpy versions | | Louvain | Modularity optimization without refinement | Legacy pipelines; Seurat default (`algorithm=1`) | Can yield internally disconnected communities; superseded by Leiden (Traag 2019) | | SLM | Smart Local Moving refinement | Seurat option (`algorithm=3`) for tighter modularity optima | Slower; rarely needed over Leiden | Scanpy 1.10 backend migration (pin for reproducibility): `sc.tl.leiden` still defaults to the `leidenalg` backend through 1.10-1.12 and emits a FutureWarning that the default will switch to igraph. The exact flip version is unconfirmed, so pin the backend explicitly: `sc.tl.leiden(adata, flavor='igraph', n_iterations=2, directed=False)`. Switching backend or `n_iterations` changes the labels - a pipeline that pins neither is non-reproducible across versions. Seurat's Leiden (`algorithm=4`) requires the `leidenalg` Python module via reticulate, which is why most Seurat pipelines still run Louvain. ## Parameter Reference `n_pcs` dominates the result far more than `n_neighbors` and is the largest under-tuned lever - too few collapses real structure, too many reintroduces technical noise. | Parameter | Typical range | Rationale | Validation | |-----------|--------------|-----------|------------| | n_pcs | 30-50 (check elbow) | Captures biological variance while denoising; effect dwarfs n_neighbors | Elbow plot; cluster stability across nearby n_pcs values | | n_neighbors | 10-30 (15 default) | Higher = smoother, fewer fine clusters; lower = more local, fragmented | Secondary lever; vary only after n_pcs is set | | resolution | 0.2-2.0 (sweep, do not fix) | Higher = more, smaller clusters; has no biological meaning | clustree across the sweep; marker check; significance test | | min_dist (UMAP) | 0.1-0.5 | Visualization only; lower = tighter visual clusters | Affects display, never the partition | Resolution is an unidentifiable nuisance parameter: it cannot be validated internally (no ground truth), so tuning it until clusters "match known cell types" is confirmation bias laundered as analysis. Sweep a range, visualize cell flow with clustree, and pick the coarsest level whose populations are defensible by orthogonal evidence - label finer splits as hypotheses. ## Cluster Cells with Scanpy **Goal:** Reduce dimensions, build the neighbor graph, partition with Leiden, and embed for display. **Approach:** PCA -> kNN graph on a chosen `n_pcs` -> Leiden with a pinned backend -> UMAP. ```python import scanpy as sc adata = sc.read_h5ad('preprocessed.h5ad') sc.tl.pca(adata, n_comps=50, svd_solver='arpack') sc.pl.pca_variance_ratio(adata, n_pcs=50, log=True) # elbow to choose n_pcs sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30) sc.tl.leiden(adata, resolution=0.5, flavor='igraph', n_iterations=2, directed=False) adata.obs['leiden'].value_counts() sc.tl.umap(adata, min_dist=0.3) sc.pl.umap(adata, color=['leiden', 'CD3D', 'MS4A1', 'CD14']) ``` ## Sweep and Validate Resolution **Goal:** Choose a defensible granularity instead of a single tuned-to-taste resolution. **Approach:** Cluster across a resolution range, inspect cell flow (clustree), and confirm each cluster carries distinct markers. ```python import scanpy as sc for res in [0.2, 0.4, 0.6, 0.8, 1.0]: sc.tl.leiden(adata, resolution=res, key_added=f'leiden_r{res}', flavor='igraph', n_iterations=2, directed=False) print(res, adata.obs[f'leiden_r{res}'].nunique(), 'clusters') sc.pl.umap(adata, color=['leiden_r0.2', 'leiden_r0.6', 'leiden_r1.0'], ncols=3) # clustree (R) or sc.tl.dendrogram for flow across resolutions; merge clusters # whose top markers are indistinguishable -> they are one population over-split ``` ## Cluster Cells with Seurat **Goal:** Run PCA, build the SNN graph, partition, and embed in Seurat. **Approach:** RunPCA -> FindNeighbors on chosen dims -> FindClusters (sweep resolutions) -> RunUMAP. ```r library(Seurat) seurat_obj <- readRDS('preprocessed.rds') seurat_obj <- RunPCA(seurat_obj, npcs = 50, verbose = FALSE) ElbowPlot(seurat_obj, ndims = 50) # choose dims seurat_obj <- FindNeighbors(seurat_obj, dims = 1:30) seurat_obj <- FindClusters(seurat_obj, resolution = c(0.2, 0.4, 0.6, 0.8, 1.0)) seurat_obj <- RunUMAP(seurat_obj, dims = 1:30) library(clustree) clustree(seurat_obj, prefix = 'RNA_snn_res.') # cell flow across the sweep DimPlot(seurat_obj, reduction = 'umap', label = TRUE) ``` `FindClusters` defaults to Louvain (`algorithm=1`); pass `algorithm=4` for Leiden (requires the leidenalg Python module). Resolutions stored as `RNA_snn_res.` columns feed clustree directly. ## Subclustering **Goal:** Resolve fine states inside a coarse cluster without importing global axes. **Approach:** Subset the cluster, then recompute HVGs, PCA, and the kNN graph on the subset. ```python sub = adata[adata.obs['leiden'] == '3'].copy() sc.pp.highly_variable_genes(sub, n_top_genes=2000) sub = sub[:, sub.var.highly_variable] sc.pp.scale(sub, max_value=10) sc.tl.pca(sub, n_comps=30) sc.pp.neighbors(sub, n_neighbors=15, n_pcs=20) sc.tl.leiden(sub, resolution=0.4, flavor='igraph', n_iterations=2, directed=False) ``` Reusing the global PCA imports axes uninformative within a homogeneous subset and manufactures artifactual sub-splits. Subclustering compounds double-dipping (cells selected twice), so stop when splits lose distinct markers or fail a significance test - not when resolution can technically still split (it always can). ## Validating That Clusters Are Real Stability and significance are separate questions, and both differ from biological reality. - Stability (necessary, not sufficient): bootstrap cells, re-cluster, measure per-cluster label agreement (Jaccard >= ~0.6-0.7 = stable). A perfectly reproducible split can still be technical - driven by cell-cycle phase, dissociation stress (FOS/JUN/HSPA1A), mitochondrial fraction, ambient RNA, or batch. Stable does not mean real. - Significance (whether a split is two populations or one): scSHC (Grabski 2023) and CHOIR (2025) test each split under a null with error control; a failed split is over-clustering and the clusters should be merged. - Biological-vs-technical adjudication: check that a split survives regressing out cycle/mito and carries non-stress markers before claiming a population. Double-dipping (post-clustering inference is invalid, not just inflated): `rank_genes_groups`/`FindAllMarkers` p-values are conditioned on a clustering chosen to maximize separation, so under one homogeneous population they do not follow their nominal null - type-I error approaches 1 as resolution rises, and BH correction does nothing because the p-values are invalid before correction. Use these marker tests for ranking and labeling only. To make a defensible claim that a cluster is a distinct population, pair a cluster significance test (scSHC/CHOIR) with a double-dipping-robust DE method (ClusterDE's synthetic null, or count splitting where the noise model holds - Poisson thinning on overdispersed counts silently reinstates the bias). See markers-annotation for marker testing and the pseudobulk path for cross-condition DE. ## UMAP and tSNE Are Visualization Only Cluster on the graph or PCA, never on embedding coordinates. Inter-cluster distances, relative sizes, and apparent gaps in UMAP/tSNE are artifacts of the embedding objective and are not metric (Chari & Pachter 2023) - do not read them as lineage, evolutionary distance, or population separation. "Cells far apart in UMAP are more different" and "UMAP preserves global structure" are folklore; PCA initialization plus high perplexity makes tSNE less misleading but does not make distances trustworthy (Kobak & Berens 2019). Embeddings display graph-derived labels; back every structural claim with the graph, validly-tested markers, or quantitative analysis in PCA space. ## Common Errors | Symptom | Cause | Fix | |---------|-------|-----| | One giant blob, no structure | Too few HVGs or wrong n_pcs (too few), or the sample is genuinely one cell type | Increase HVGs (~2000), raise n_pcs, check the elbow plot; if markers stay uniform across a resolution sweep, the blob may be a single real population, not a parameter bug | | Far too many clusters | Resolution too high; n_pcs too high (noise) | Lower resolution; sweep with clustree; reduce n_pcs to the elbow | | Adjacent clusters share all top markers | Over-clustering one population | Merge them; lower resolution; significance-test the split (scSHC/CHOIR) | | Labels change between runs/versions | Leiden backend or n_iterations unpinned | Pin `flavor='igraph', n_iterations=2, directed=False`; set random_state | | A cluster maps to one sample/lane only | Batch effect, not biology | Integrate batches first (batch-integration); inspect QC covariates | | A "stable" cluster of stress/cycle genes | Technical split (dissociation, cycle, mito) | Regress out cycle/mito or score and exclude; require non-stress markers | | Cluster expresses two lineages' markers | Doublets clustering together | Run doublet detection before clustering (doublet-detection) | | Marker p-values quoted as proof clusters are real | Double-dipping (selective inference) | Use markers for ranking only; validate with scSHC/CHOIR + ClusterDE | ## Related Skills - preprocessing - QC, normalization, and HVG selection that must precede clustering - doublet-detection - Remove doublets before clustering so they do not form fake intermediate clusters - batch-integration - Integrate batches before clustering when a cluster tracks a single sample - markers-annotation - Find and test markers per cluster (with the double-dipping caveat) - cell-annotation - Assign cell-type identities to validated clusters - single-cell/differential-abundance - Test whether cluster proportions shift across conditions - data-visualization/dimensionality-reduction-plots - Publication-quality UMAP/tSNE/PCA figures - pathway-analysis/go-enrichment - Interpret per-cluster marker sets ## References - Traag, Waltman & van Eck (2019). From Louvain to Leiden: guaranteeing well-connected communities. Sci Rep 9:5233. - Chari & Pachter (2023). The specious art of single-cell genomics. PLoS Comput Biol 19(8):e1011288. - Kobak & Berens (2019). The art of using t-SNE for single-cell transcriptomics. Nat Commun 10:5416. - Grabski, Street & Irizarry (2023). Significance analysis for clustering with single-cell RNA-sequencing data. Nat Methods 20:1196-1202. - Neufeld, Gao, Popp, Battle & Witten (2024). Inference after latent variable estimation for single-cell RNA-seq (count splitting). Biostatistics 25(1):270-287. - Zappia & Oshlack (2018). Clustering trees: a visualization for evaluating clusterings at multiple resolutions. GigaScience 7(7):giy083.