How to Make and Read a UMAP of Your Own RNA-seq Data

By the end of this you will have a two-dimensional embedding of RNA-seq samples or cells, produced by a pipeline you can rerun with different parameters, plus a small set of checks that tell you which features of the picture are real and which are artifacts of the algorithm. You need a count matrix (genes × samples for bulk, genes × cells for single-cell), Python 3.10+ with scanpy, umap-learn ≥ 0.5, leidenalg, and scikit-learn, and, if you are starting from bulk FASTQ files, salmon and R with tximport and DESeq2. The one prerequisite people skip is sample count. UMAP on twelve longitudinal blood draws from one person will produce a picture, and that picture will be noise. If you have fewer than a few hundred points, either embed your samples into a public reference cohort (recount3, GTEx) or skip UMAP entirely and look at PCA.

1. Build a count matrix you trust

For bulk RNA-seq, quantify against a transcriptome index rather than aligning to the genome unless you specifically need splice-junction evidence. Salmon’s selective alignment with bias correction is what we use:

salmon index -t gencode.v44.transcripts.fa.gz -d decoys.txt \
  -i salmon_idx_v44 -k 31 --gencode -p 16

salmon quant -i salmon_idx_v44 -l A \
  -1 s01_R1.fastq.gz -2 s01_R2.fastq.gz \
  --validateMappings --gcBias --seqBias --posBias \
  --numGibbsSamples 20 -p 16 -o quants/s01

Then collapse transcripts to genes in R, keeping the offsets that correct for effective length differences between samples:

library(tximport); library(DESeq2)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")
dds <- DESeqDataSetFromTximport(txi, colData, ~ 1)
dds <- dds[rowSums(counts(dds) >= 10) >= 3, ]
vsd <- vst(dds, blind = TRUE)
write.csv(assay(vsd), "vst_matrix.csv")

The variance-stabilizing transform matters more than most people expect. Raw counts have variance proportional to the mean, so a plain Euclidean distance on counts is dominated by a few highly expressed genes (in whole blood, hemoglobin and ribosomal transcripts will decide your entire embedding). VST, or rlog for small sample counts, flattens that relationship so that a distance computed across genes reflects coordinated expression changes rather than sequencing depth. blind = TRUE is the right choice when the embedding is exploratory and you do not want the design to inform the dispersion estimates.

For single-cell data starting from a CellRanger or STARsolo output directory, load the filtered matrix and do quality control before anything else:

import scanpy as sc, numpy as np
adata = sc.read_10x_mtx("filtered_feature_bc_matrix/", var_names="gene_symbols")
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True, log1p=False)
adata = adata[(adata.obs.n_genes_by_counts > 200) &
              (adata.obs.n_genes_by_counts < 6000) &
              (adata.obs.pct_counts_mt < 15)].copy()
adata.layers["counts"] = adata.X.copy()

Those thresholds are starting points, not rules. Plot the distributions first: a 15% mitochondrial cutoff discards most cardiomyocytes and hepatocytes, and a 6000-gene ceiling that is meant to catch doublets will also remove legitimately large cells. Run a dedicated doublet caller (scrublet or DoubletFinder) rather than relying on gene-count ceilings, because doublets are the single most common source of fake intermediate populations in a UMAP.

2. Select features and reduce with PCA first

UMAP should almost never see the raw gene space. Run principal component analysis first, then run UMAP on the top PCs. This is the step that answers the “PCA versus UMAP” question: they are not competitors in a pipeline, they are sequential. PCA is a linear projection that finds orthogonal directions of maximum variance, so its axes have loadings you can read gene by gene, distances in PC space are meaningful, and it removes most of the technical noise cheaply. UMAP is a nonlinear graph embedding that gives up interpretable axes and meaningful global distances in exchange for showing local neighborhood structure in two dimensions.

sc.pp.highly_variable_genes(adata, flavor="seurat_v3",
                            n_top_genes=2000, layer="counts")
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata
adata = adata[:, adata.var.highly_variable].copy()
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])  # optional, see below
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50, svd_solver="arpack")
sc.pl.pca_variance_ratio(adata, n_pcs=50, log=True)

Note the order: seurat_v3 highly variable gene selection expects raw counts and models the mean-variance relationship directly, so it runs on the counts layer before normalization. The flavor="seurat" default expects log-transformed data instead. Mixing these up is a quiet way to select the wrong 2000 genes.

We usually skip regress_out. It is slow, it densifies sparse matrices, and regressing out total counts after normalization frequently removes real biological signal from cell types that genuinely differ in RNA content. Use it only when a QC covariate visibly structures your PCA.

For bulk VST data, the equivalent is to keep the top 1000 to 3000 genes by variance and run PCA on the centered matrix. With 50 samples, 10 to 20 PCs is plenty. For single-cell data, read the variance ratio plot and take PCs up to where the curve flattens, typically 20 to 40. The choice is less consequential than it feels: going from 30 to 40 PCs usually shifts cluster boundaries slightly and rarely changes which populations exist.

3. Build the neighbor graph

UMAP’s input is a k-nearest-neighbor graph, and n_neighbors is the parameter that controls everything downstream. It sets the size of the local neighborhood the algorithm tries to preserve: small values (5 to 10) emphasize fine local structure and fragment the embedding, large values (50 to 100) emphasize broader relationships and merge subpopulations.

sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30, metric="cosine")

We prefer cosine distance over Euclidean on log-normalized scRNA-seq PCs, because it is insensitive to residual depth differences that survive normalization. For bulk VST data, Euclidean is fine and more interpretable.

If your samples come from different sequencing runs, library prep batches, or draw dates months apart, integrate before this step. Run scanpy.external.pp.harmony_integrate or scvi-tools, then pass use_rep="X_pca_harmony" to sc.pp.neighbors. Do not judge whether you have a batch effect by looking at the UMAP alone. Color the PCA by batch and compute a metric such as kBET or the average silhouette width across batch labels, because UMAP will happily produce a picture in which batches look mixed when the underlying graph says otherwise.

4. Run UMAP with a deterministic initialization

The default in most pipelines is a spectral initialization, which is stochastic in sign and scale and is the main reason two runs of the same data produce embeddings that look globally different. Initializing from the first two principal components instead makes the global arrangement of clusters reproducible and much more faithful to the data, which Kobak and Linderman demonstrated across both t-SNE and UMAP: informative initialization is what preserves global structure, and the nonlinear optimization itself mostly handles local neighborhoods.¹

import numpy as np
X_pca = adata.obsm["X_pca"][:, :2].copy()
X_pca = X_pca / np.std(X_pca[:, 0]) * 1e-4   # small scale, as recommended
adata.obsm["X_pca_init"] = X_pca

sc.tl.umap(adata, min_dist=0.3, spread=1.0,
           init_pos="X_pca_init", random_state=0, n_components=2)
sc.pl.umap(adata, color=["leiden", "total_counts", "pct_counts_mt"], ncols=3)

min_dist controls how tightly points are allowed to pack within a cluster. It is purely cosmetic with respect to the graph: 0.1 gives dense, well-separated blobs that are easier to see clusters in, and 0.5 gives a diffuse cloud that is better for seeing gradients such as differentiation trajectories. Neither is more correct. Pick one, state it in your figure caption, and do not read anything into the resulting density.

Then cluster on the same graph UMAP used, so that the colors and the layout come from the same object:

sc.tl.leiden(adata, resolution=1.0, flavor="igraph",
             n_iterations=2, directed=False, key_added="leiden")
sc.tl.rank_genes_groups(adata, "leiden", method="wilcoxon", pts=True)
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5)

Cluster identity comes from marker genes, not from position on the plot. Sweep resolution from 0.2 to 1.5 and watch which clusters are stable across the sweep. Clusters that appear only at resolution 1.4 and have no distinguishing markers with adjusted p < 0.01 and a meaningful log fold change are over-partitioning.

5. Check the embedding before you believe it

Three checks catch most of the ways a UMAP misleads. First, quantify neighborhood preservation rather than eyeballing it:

from sklearn.manifold import trustworthiness
tw = trustworthiness(adata.obsm["X_pca"][:, :30],
                     adata.obsm["X_umap"], n_neighbors=15)

Values above roughly 0.9 indicate that most k-nearest neighbors in PC space remain neighbors in the embedding. A low value means the 2D picture is discarding structure you care about.

Second, rerun with n_neighbors in {5, 15, 50} and random_state in {0, 1, 2}, and overlay the cluster labels from the original run. Features that survive all six runs are real. A gap that closes when you move from 15 to 50 neighbors is a property of the algorithm.

Third, compare against a second method. Benchmarks of dimensionality reduction on transcriptome data show that no single method dominates across tasks, and that which representation best captures a biological signal depends on the question being asked.² PaCMAP and TriMap make different local-versus-global tradeoffs than UMAP and are worth running as a sanity check; PaCMAP embeddings have been used as input representations for downstream multi-omics models for exactly this reason.³ If a structure appears in UMAP, PaCMAP, and the first three PCs, it is in the data.

6. Read the plot correctly

Here is what a UMAP legitimately tells you. Points that sit together were near-neighbors in the high-dimensional space, so cluster membership is meaningful. The presence of distinct, well-populated groups indicates discrete states. Continuous ribbons and bridges suggest gradual transitions, though you should confirm with a trajectory method rather than trusting the shape.

Here is what it does not tell you. The distance between two clusters is not proportional to their transcriptional dissimilarity, because UMAP’s loss function only constrains the attractive and repulsive forces on the neighbor graph. With PCA initialization the coarse arrangement carries some signal,¹ but you still cannot say cluster A is twice as far from C as from B. The area a cluster occupies has nothing to do with how many cells are in it or how variable it is. The axes have no units and no meaning, so never report “UMAP1” as a variable. Above all, never run a statistical test on UMAP coordinates. Differential expression, correlation, and regression belong in the normalized expression space or in PC space, where distances are metric.

Two practical answers to the recurring comparisons. UMAP is typically faster than Barnes-Hut t-SNE on datasets above about 50,000 points, largely because of approximate nearest-neighbor search, though FIt-SNE closes most of that gap. t-SNE remains preferable when you want strict local-neighborhood fidelity and are willing to accept that global arrangement is meaningless. PCA is what you should use whenever you need interpretable loadings, when your sample count is small, or when the output feeds a downstream statistical model.

7. Embed new samples into an existing reference

For longitudinal personal data, the useful operation is projecting each new draw into a fixed reference embedding so that successive timepoints are comparable. Fit once, then transform:

import umap, joblib
reducer = umap.UMAP(n_neighbors=15, min_dist=0.3, init="pca",
                    metric="euclidean", random_state=0).fit(ref_pcs)
joblib.dump(reducer, "umap_reference.pkl")

new_pcs = pca_model.transform(new_vst_matrix)   # same PCA, same genes, same scaler
new_emb = reducer.transform(new_pcs)

The projection is only valid if the new sample passed through an identical preprocessing path: same gene set, same VST parameters, same PCA rotation. Refitting the PCA on the combined data invalidates the saved reducer. Also be aware that transform places new points using the existing graph, so a genuinely novel cell state will be forced into the nearest existing region rather than appearing as a new island. This approach works well for tracking sample-level heterogeneity in bulk transcriptomes, where UMAP has been shown to resolve subgroup structure that linear projections blur together.⁴ Applied to tissue single-cell data, the same workflow is how groups identify and then follow specific populations such as inflammatory and prehypertrophic chondrocytes across conditions.⁵

If your data include more than one modality, do not embed each separately and compare pictures. Joint generalizations of t-SNE and UMAP learn a single embedding with per-modality weights, which is a better-posed problem than eyeballing two plots side by side.⁶

Common problems

The embedding is dominated by one axis of technical variation. Color by total_counts, pct_counts_mt, and batch. If any of them paints a gradient across the whole plot, fix it upstream through filtering or integration rather than by rotating the picture.

Every rerun looks different. You did not set random_state, or you are using spectral initialization. Set both random_state and init_pos to a PCA-derived embedding.¹ Also pin your versions: umap-learn and numba releases have changed default behavior, and embeddings are not comparable across versions.

Clusters appear that have no markers. Check for doublets first, then drop the Leiden resolution. UMAP will separate groups that differ only by ambient RNA contamination, so run SoupX or CellBender on droplet data if you have not.

A cell type you expect is missing. Look at your QC thresholds before concluding anything. Also check whether n_neighbors is large enough to have merged it into a neighbor, and search your marker genes directly in the full expression matrix rather than in the plot.

Bulk samples separate perfectly by collection date. That is a batch effect until proven otherwise. With longitudinal personal data, collection date and biology are confounded by design, and no algorithm can separate them. Randomize library prep order across timepoints if you control the sample handling, and keep an aliquot of a fixed reference RNA in every batch.

You see a structure and want to know what it means for you. A UMAP is a visualization of variance, not a clinical finding. Marker genes, cluster proportions, and shifts across timepoints are hypotheses to investigate, and any interpretation touching your health belongs with a clinician who can see the full picture. Published multi-omics studies that use these embeddings pair them with orthogonal validation before drawing biological conclusions, and so should you.⁷

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Dmitry Kobak, George C. Linderman. Initialization is critical for preserving global data structure in both t-SNE and UMAP. Nature Biotechnology, 2021. https://doi.org/10.1038/s41587-020-00809-z ↩ ↩² ↩³

2. Yuseong Kwon, Sojeong Park, Soyoung Park, et al. Benchmarking of dimensionality reduction methods to capture drug response in transcriptome data. Scientific Reports, 2025. https://doi.org/10.1038/s41598-025-12021-7 ↩

3. Hazem Qattous, Mohammad Azzeh, Rahmeh Ibrahim, et al. PaCMAP-embedded convolutional neural network for multi-omics data integration. Heliyon, 2024. https://doi.org/10.1016/j.heliyon.2023.e23195 ↩

4. Yang Yang, Hongjian Sun, Yu Zhang, et al. Dimensionality reduction by UMAP reinforces sample heterogeneity analysis in bulk transcriptomic data. Cell Reports, 2021. https://doi.org/10.1016/j.celrep.2021.109442 ↩

5. Yue Fan, Xuzhao Bian, Xiaogao Meng, et al. Unveiling inflammatory and prehypertrophic cell populations as key contributors to knee cartilage degeneration in osteoarthritis using multi-omics data integration. Annals of the Rheumatic Diseases, 2024. https://doi.org/10.1136/ard-2023-224420 ↩

6. Van Hoan Do, Stefan Canzar. A generalization of t-SNE and UMAP to single-cell multimodal omics. Genome Biology, 2021. https://doi.org/10.1186/s13059-021-02356-5 ↩

7. Jiatao Zhang, Wenhao Zhou, Na Li, et al. Multi-omics analysis unveils immunosuppressive microenvironment in the occurrence and development of multiple pulmonary lung cancers. npj Precision Oncology, 2024. https://doi.org/10.1038/s41698-024-00651-5 ↩