UMAP/t-SNE clustering of thousands of cells by gene expression — mapping the transcriptome to find rare, therapeutically relevant cell populations
Every scRNA-seq experiment begins by converting a piece of tissue or a blood sample into a suspension of individual, viable cells, then physically isolating thousands of them one-by-one so their RNA can be labeled separately. Droplet microfluidics — commercialized by 10x Genomics' Chromium platform — made this step routine, fast, and cheap enough to profile tens of thousands of cells in a single afternoon.
Solid tissue is enzymatically and mechanically dissociated (collagenase/trypsin digestion, gentleMACS homogenization) into single cells, then filtered through a 40 µm strainer to remove clumps and debris. Blood samples are simpler: density-gradient centrifugation (Ficoll) isolates peripheral blood mononuclear cells (PBMCs) directly. Viability is checked by trypan blue or AO/PI staining — a healthy input suspension should be >80–90% viable, since dying cells leak RNA and generate ambient background contamination that confounds downstream analysis.
Cell concentration is precisely counted (hemocytometer or automated counter) and diluted to a target loading concentration, because everything downstream depends on getting the cell:droplet ratio right.
The 10x Genomics Chromium controller pushes a suspension of cells, a suspension of barcoded gel beads, and oil through a microfluidic chip with a T-junction. At the junction, aqueous phase pinches off into nanoliter droplets suspended in oil — Gel Beads-in-Emulsion, or GEMs.
Each gel bead carries ~750,000 copies of an oligonucleotide primer containing: a 16 bp cell barcode (identical across all oligos on one bead, unique per bead), a 10–12 bp unique molecular identifier (UMI, different on every oligo copy), and a poly-dT sequence that primes reverse transcription of polyadenylated mRNA. Because bead barcode sequences are drawn from a curated whitelist of ~3.5 million sequences, the probability that two beads share the same barcode in one run is negligible — every captured cell's transcriptome can be uniquely traced back through sequencing.
Cell encapsulation into droplets is a random, independent process, so the number of cells per droplet follows a Poisson distribution: P(k cells) = e^(−λ) λᵏ / k!, where λ is the average number of cells loaded per droplet.
To keep multiplet (doublet/triplet) rates low, operators deliberately under-load: λ is kept around 0.05–0.1, which wastes most droplets (only ~5–10% capture a cell at all) but keeps the doublet rate to roughly 0.4–1.6% per 1,000 cells targeted. Loading more cells to increase yield rapidly increases the doublet rate — at λ≈0.3 (targeting ~10,000 cells) doublet rates climb to several percent, so computational doublet detection (Scrublet, DoubletFinder) is run downstream as a safety net regardless.
Alternative capture platforms trade throughput for other advantages: plate-based Smart-seq2/3 captures full-length transcripts at much higher sensitivity but only hundreds of cells per plate; combinatorial indexing (sci-RNA-seq, SPLiT-seq) avoids microfluidics entirely and scales to millions of cells via split-pool barcoding; spatial platforms (10x Visium, Xenium) sacrifice single-cell resolution or throughput to preserve tissue location.
Inside each droplet, captured mRNA is reverse-transcribed, the droplets are broken, and all barcoded cDNA is pooled into one sequencing library. After Illumina sequencing, computational pipelines convert raw reads into a clean, cell-by-gene expression matrix — but only after aggressively filtering out empty droplets, doublets, and dying cells that would otherwise corrupt every downstream analysis.
CellRanger count (or the open-source STARsolo / Salmon alevin) processes raw sequencing reads through a fixed pipeline: extract and error-correct the 16 bp cell barcode against the whitelist (up to 1 mismatch tolerated); extract the UMI; splice-aware align the cDNA read (typically 91 bp, Read 2) to a reference genome (e.g. GRCh38) using STAR; assign each read to a gene by exonic overlap.
The critical deduplication step: multiple reads sharing the same cell barcode + UMI + gene originated from PCR amplification of the same original mRNA molecule, so they are collapsed and counted once. This UMI counting removes PCR amplification bias entirely — unlike bulk RNA-seq, the final count directly reflects molecule number, not amplification efficiency. Output is a sparse matrix (~20,000–36,000 genes × thousands of cells, >90% zeros) stored in .h5 or MatrixMarket (.mtx) format.
Three metrics computed per barcode drive the standard QC gate (Seurat / Scanpy defaults, tissue-dependent):
• nFeature (genes detected): filter out barcodes with too few genes (<200–500, likely empty droplets or ambient RNA) and unusually high genes (>2.5x median, likely doublets) • nCount (total UMIs): extreme outliers in either direction flagged • percent.mt (fraction of UMIs from mitochondrial genes): cells with ruptured membranes leak cytoplasmic mRNA but retain mitochondria, so a high mitochondrial fraction (commonly >10–20%, tissue-dependent) flags a dying or stressed cell
Ambient RNA — free-floating transcripts from lysed cells in the original suspension that get incorporated into unrelated droplets — is corrected computationally with tools like SoupX or CellBender before or alongside this filtering. From a typical ~14,000 loaded cells, doublet removal + QC filtering + ambient correction commonly leaves 60–75% (roughly 8,000–10,000) as high-confidence single-cell profiles.
Raw UMI counts are not directly comparable between cells because sequencing depth varies. Standard log-normalization scales each cell's counts to a common total (e.g. 10,000), then applies log1p; the newer SCTransform (Hafemeister & Satija, Genome Biology 2019) instead models counts with a regularized negative binomial regression, better handling the mean-variance relationship inherent to count data.
Of the ~20,000 detected genes, most are uninformative for distinguishing cell types (housekeeping genes, low-count noise). A highly variable gene (HVG) selection step — using a variance-stabilizing transform to rank genes by variance after accounting for mean expression — reduces the matrix to roughly 2,000 informative genes, which become the input to PCA in the next stage.
A single cell's expression profile lives in a space with as many dimensions as there are genes measured — around 2,000 after variable-gene selection. No human can look at a 2,000-dimensional point cloud. Dimensionality reduction algorithms compress that space down to 2 (or 3) dimensions for visualization while trying to preserve the property that matters most: cells with similar expression profiles should end up close together.
Principal Component Analysis is almost always run first, on the ~2,000-gene normalized, scaled expression matrix. PCA finds orthogonal linear combinations of genes ("components") ordered by how much expression variance they explain. The first PC alone often separates the broadest biological distinction in the dataset (e.g. immune vs. non-immune, or cell-cycle phase); subsequent PCs capture progressively finer structure.
An "elbow plot" of variance explained per component (or a statistical test like JackStraw) is used to choose how many PCs to keep — typically 30–50. Every cell is now represented by a vector of 30–50 numbers instead of 2,000 gene counts: a large denoising step, since much of single-gene variation is technical noise averaged away by combining many correlated genes.
Working in the reduced PCA space, each cell is connected to its k nearest neighbors by Euclidean distance (k typically 15–30). This kNN graph is the shared computational backbone for two seemingly different downstream tasks: it is the direct input to UMAP's and t-SNE's low-dimensional embedding optimization in this stage, and it is also the graph that Leiden clustering will partition into communities in the next stage. The choice of k (n_neighbors) directly trades off local detail against global structure — small k emphasizes fine local neighborhoods (more fragmented, textured layouts); large k smooths over local detail to better preserve broad relationships between distant clusters.
UMAP (Uniform Manifold Approximation and Projection; McInnes, Healy & Melville, 2018) constructs a fuzzy topological representation of the high-dimensional kNN graph (a "fuzzy simplicial set"), then uses stochastic gradient descent to find a low-dimensional embedding whose fuzzy simplicial set is as similar as possible (minimizing cross-entropy) to the high-dimensional one. Its n_neighbors parameter directly controls this graph's locality, and min_dist controls how tightly points are allowed to pack together in the embedding.
t-SNE (t-distributed Stochastic Neighbor Embedding; van der Maaten & Hinton, 2008) instead converts pairwise distances into conditional probabilities in high-dimensional space, and optimizes a low-dimensional layout so that a Student t-distribution over pairwise distances matches those probabilities as closely as possible, minimizing KL divergence via gradient descent. Its "perplexity" parameter (typically 30–50) plays a role similar to n_neighbors. The classic O(N²) computation is intractable beyond a few thousand cells, but the Barnes-Hut approximation (and later FIt-SNE) reduced this to O(N log N), making t-SNE practical at single-cell atlas scale.
In practice: UMAP tends to better preserve some global/inter-cluster distance relationships and runs faster at large N; t-SNE often produces more tightly separated, visually crisp local clusters but can distort relative distances between distant clusters. Most modern pipelines (Seurat, Scanpy) default to UMAP for the final published figure while still computing t-SNE for comparison.
A UMAP plot only shows that groups of cells exist — it does not, by itself, define discrete clusters or tell you what those groups are biologically. That requires a graph-partitioning algorithm to formally assign cluster labels, followed by differential expression testing against known marker genes to translate anonymous cluster numbers into recognized cell types like "CD8+ cytotoxic T cell" or "classical monocyte."
Clustering operates on the same kNN graph built in the previous stage, treating cluster discovery as a graph community-detection problem: find groups of nodes (cells) that are more densely connected to each other than to the rest of the graph, by optimizing a quality function called modularity.
The Louvain algorithm (Blondel et al., 2008) greedily moves nodes between communities to locally improve modularity, then aggregates communities into super-nodes and repeats — fast and widely used, but with a known flaw: it can produce clusters that are internally disconnected or badly connected subgraphs, an artifact of its aggregation step. The Leiden algorithm (Traag, Waltman & van Eck, Scientific Reports 2019) refines this by explicitly guaranteeing every output community is well-connected, while also converging faster and finding higher-quality (higher modularity) partitions. A user-set "resolution" parameter directly controls granularity: low resolution merges cells into a few broad clusters (e.g. 6 major immune lineages), high resolution splits them into many fine-grained subclusters (e.g. 20+, distinguishing naive vs. memory vs. effector T-cell subsets).
Once clusters exist, each is annotated by testing which genes are significantly over-expressed in that cluster versus all other cells (or versus each other cluster pairwise) — most commonly with a Wilcoxon rank-sum test (Seurat's FindMarkers, Scanpy's rank_genes_groups), thresholded at log2 fold-change >0.25 and adjusted p-value <0.05 by default. The resulting ranked marker gene lists are cross-referenced against decades of immunology and cell-biology literature: a cluster enriched for CD3D/CD3E/IL7R is confidently called T cells; MS4A1/CD19/CD79A marks B cells; CD14/LYZ/FCGR3A marks monocytes; NKG7/GNLY/KLRD1 marks NK cells; CD1C/FCER1A marks dendritic cells. This manual, literature-grounded annotation step is still considered the gold standard, though it does not scale trivially to atlas-sized datasets with dozens of clusters.
To speed up and standardize annotation, automated classifiers map a query dataset's cells onto a pre-labeled reference atlas rather than requiring manual marker inspection each time: SingleR correlates each cell's expression profile against reference cell-type signatures; CellTypist uses a logistic-regression model trained on millions of immune cells; Seurat's Azimuth performs reference-based label transfer using an integrated reference UMAP. These tools are increasingly trained on massive community reference atlases — the Human Cell Atlas and Tabula Sapiens (~500,000 cells spanning 24 human tissues) — turning cell-type annotation from a bespoke manual exercise into a standardized, reproducible lookup against a shared reference map.
| Product | Indication | Trial Design | Key Result |
|---|---|---|---|
| T cells | CD3D, CD3E, IL7R | ~45% of PBMCs · adaptive cellular immunity | CD4 helper / CD8 cytotoxic subsets resolvable at higher resolution |
| B cells | MS4A1, CD19, CD79A | ~14% of PBMCs · antibody production | Naive vs. memory B cells separable via CD27 |
| Monocytes | CD14, LYZ, FCGR3A | ~24% of PBMCs · phagocytosis, antigen presentation | Classical (CD14+) vs. non-classical (FCGR3A+) split |
| NK cells | NKG7, GNLY, KLRD1 | ~12% of PBMCs · innate cytotoxicity | Viral / tumor surveillance without prior sensitization |
| Dendritic cells | CD1C, FCER1A, CLEC9A | ~2–5% of PBMCs · antigen presentation, T-cell priming | Bridges innate detection and adaptive activation |
| Rare / stem-like | PROM1, CD34 (context-dependent) | <2% of cells · progenitor / persister / CTC-like state | Frequently the actual therapeutic target population |
The most clinically valuable cell populations are often the rarest ones on the plot: a handful of drug-persister cells surviving chemotherapy, a stem-like progenitor subclone driving relapse, or circulating tumor cells at a frequency of one in a billion blood cells. Sub-clustering, careful artifact control, and trajectory inference let these populations be resolved, characterized, and tracked — turning a descriptive atlas into an actionable target-discovery tool for precision dosing.
Detecting a rare subpopulation reliably requires enough total cells sequenced to expect several representatives of it: a rough statistical rule of thumb is to profile at least ~10× the reciprocal of the expected population frequency, so resolving a 0.5% population comfortably wants several thousand total cells, and resolving a 0.05% population wants tens of thousands. Coarse clustering often buries small subpopulations inside a larger, transcriptionally similar cluster (e.g. a stem-like subset hidden within the broader monocyte cluster); iterative sub-clustering — reclustering only the cells within one coarse cluster at higher resolution — frequently reveals structure invisible in the first pass.
A critical caveat: apparent rare clusters are disproportionately likely to be technical artifacts rather than real biology — residual doublets, ambient RNA contamination, or stressed/dying cells with distorted expression profiles can all masquerade as a distinct rare population. Rigorous QC (Stage 2) and doublet scoring must be revisited specifically for any candidate rare cluster before it is trusted.
Many biological processes — differentiation, activation, disease progression — are continuous rather than discrete, so cells captured at a single time point actually sample many different points along an underlying continuum. Trajectory inference algorithms exploit this: Monocle3 (Trapnell lab) builds a principal graph through the expression manifold and orders cells along it into "pseudotime," an inferred progression coordinate that substitutes for real time-series sampling. PAGA (Partition-based Graph Abstraction; Wolf et al., 2019) instead abstracts the whole dataset into a simplified graph of clusters and their connectivity, useful for mapping coarse lineage relationships in very large, complex atlases.
RNA velocity (La Manno et al., 2018; scVelo, Bergen et al., 2020) adds directionality that pseudotime alone cannot provide: by comparing the ratio of unspliced (intronic) to spliced (mature) mRNA per gene per cell, it estimates each cell's instantaneous transcriptional "velocity" — predicting which cell state it is moving toward next, producing an arrow field over the UMAP embedding rather than just an ordering.
Once a rare population is confidently resolved and characterized by its marker genes, it becomes an actionable target rather than a curiosity. In oncology, serial scRNA-seq biopsies can track the fractional expansion of a rare drug-resistant clone during treatment, directly informing dose escalation, combination therapy timing, or early switching decisions — a core input to AI-driven precision dosing pipelines. In cell therapy manufacturing, scRNA-seq profiling of a CAR-T cell product characterizes the proportion of stem-like memory cells (associated with durable clinical response) versus terminally exhausted cells, feeding potency-based dosing decisions rather than cell count alone. Population-scale reference atlases — the Human Cell Atlas, Tabula Sapiens, and disease-specific atlases from the Human Tumor Atlas Network — provide the baseline maps against which a rare population in a new patient sample is recognized and interpreted at all.