RNA Velocity & Fate Mapping in Single-Cell Data: scVelo, CellRank and PAGA

Go beyond pseudotime — add direction with RNA velocity, ask where each cell is heading with CellRank, and map the coarse topology with PAGA, plus the honest R interoperability path

Pseudotime orders cells but does not tell you which way they are moving. RNA velocity does. This lesson explains RNA velocity (scVelo), fate mapping (CellRank 2) and topology (PAGA) — what they add beyond Slingshot and Monocle3, why they need spliced and unspliced counts, and how you would run them. These are the Python state of the art, so the code here is illustrative; you also get the R-to-Python bridge and the managed-compute handoff for running velocity for real.

TipKey takeaways

• Pseudotime gives order; velocity gives direction. Trajectory methods like Slingshot and Monocle3 order cells along a continuum but cannot say which way the arrow of change points. RNA velocity — the time derivative of the gene-expression state — recovers that direction from the balance of unspliced and spliced mRNA (La Manno et al., 2018).
• Velocity needs a special kind of data. It requires spliced and unspliced count matrices, which come from velocyto or kallisto | bustools, not from an ordinary gene-by-cell count matrix. If your data has no unspliced layer, you cannot compute velocity — this is the first thing to check.
• The state of the art is a Python stack. scVelo estimates velocity, CellRank 2 turns velocity into fate probabilities (“where is each cell heading?”), and PAGA draws a coarse-grained topology map. All three are Python tools built on AnnData (the annotated-data container).
• This lesson is illustrative — nothing here executes. The velocity workflow runs in Python on velocity-ready data, which this reproducible R lesson does not have. Every snippet below is shown for reading, not run; copy them into a real Python environment with spliced/unspliced counts to reproduce the results.
• R users are not locked out. You can move a Seurat or SingleCellExperiment object into AnnData with zellkonverter / anndata / sceasy, and the velociraptor bridge calls scVelo from R. But you still need velocity-ready data and a real Python environment — so the practical answer is to run velocity in a managed-compute environment, not in the browser.

Introduction

You have ordered your cells. Using Slingshot or Monocle3 you fitted a trajectory through a differentiating population and assigned each cell a pseudotime — its position along the continuum. But pseudotime has a blind spot: it is an ordering, not a direction. Nothing in the geometry of a fitted curve tells you whether cells flow from the stem end to the differentiated end or the other way around. You supply the direction yourself, by rooting the trajectory at a cluster you already believe is the origin. When you do not know the origin — regeneration, reprogramming, a disease process running in reverse — pseudotime alone cannot help.

RNA velocity removes that blind spot. It reads the direction of change directly from the data, cell by cell, and points an arrow into the near future. Layered on top, fate mapping asks the next question — given its velocity, where is each cell most likely to end up? — and topology methods draw a simplified map of how the major cell groups connect. Together they are the level beyond the ordering you already have.

There is one honest complication, and it shapes this entire lesson. The mature, well-validated tools for velocity, fate mapping and topology are Python packages, and they need a kind of input — spliced and unspliced counts — that a standard single-cell RNA sequencing (scRNA-seq) count matrix does not contain. Because of that, this lesson teaches the methods illustratively: no code here is executed. We explain what each tool does and show exactly how you would call it, then hand execution off to a real Python environment. The section why nothing here runs makes the reasoning explicit — it is part of the lesson, not an apology for it.

If you have not yet met the ordering methods this builds on, read the Slingshot and Monocle3 deep dives and the Slingshot vs Monocle3 chooser first. This lesson assumes you already have a trajectory and want to add direction to it.

What RNA velocity actually measures

The idea is elegant and rests on a fact of molecular biology. When a gene is transcribed, the cell first makes unspliced pre-messenger RNA (pre-mRNA, still carrying introns), then the spliceosome removes the introns to produce mature spliced mRNA, which is eventually degraded. Standard scRNA-seq protocols capture both species — you can count reads that map across intron–exon boundaries (unspliced) separately from reads that map to mature transcripts (spliced).

That separation is the whole trick. Consider one gene in one cell:

• If the gene has just been switched on, unspliced pre-mRNA accumulates faster than it can be spliced and degraded — so the unspliced-to-spliced ratio is above its steady-state value. The gene is being induced: its future spliced level will rise.
• If the gene has just been switched off, unspliced production drops, the existing spliced pool keeps degrading, and the ratio falls below steady state. The gene is being repressed: its future spliced level will fall.

Do this for every gene and you get, for each cell, a high-dimensional velocity vector — the predicted rate of change of its expression state over the next few hours (La Manno et al., 2018). Project those vectors onto your uniform manifold approximation and projection (UMAP, the 2-D map) and you get the famous velocity streamlines: arrows flowing across the embedding, showing the direction of differentiation without you having to root anything.

The relationship between unspliced and spliced mRNA for a single gene is often drawn as a phase portrait — a scatter of cells with spliced level on the x-axis and unspliced level on the y-axis. Cells above the diagonal steady-state line are inducing the gene; cells below it are repressing it. scVelo fits this phase portrait per gene and reads off the velocity.

The original method estimated velocity from a steady-state (deterministic) model. scVelo generalized it: its stochastic model uses the second moments of the unspliced/spliced distribution, and its dynamical model solves the full splicing kinetics with an expectation-maximization (EM, an iterative likelihood-fitting procedure) algorithm — which makes velocity valid even for transient cell states that are not at steady state (Bergen et al., 2020).

The velocity-to-fate workflow at a glance

Before any code, here is the shape of the whole pipeline — what feeds what, from raw spliced/unspliced counts through to fate probabilities and a topology map.

flowchart TD
  raw["Spliced + unspliced counts<br/>(velocyto / kallisto|bustools)"] --> adata["AnnData object<br/>(spliced & unspliced layers)"]
  adata --> moments["scVelo: moments"]
  moments --> velocity["scVelo: velocity<br/>(stochastic / dynamical)"]
  velocity --> vgraph["scVelo: velocity graph"]
  vgraph --> stream["Velocity streamlines<br/>on the UMAP"]
  vgraph --> vk["CellRank: VelocityKernel<br/>+ ConnectivityKernel"]
  vk --> gpcca["CellRank: GPCCA estimator<br/>macrostates -> terminal states"]
  gpcca --> fate["Per-cell fate probabilities<br/>(where is each cell heading?)"]
  adata --> paga["scanpy: PAGA"]
  paga --> topo["Coarse-grained topology graph<br/>of the clusters"]

Figure 1

Read it left to right in three stages: scVelo turns spliced/unspliced counts into velocity and streamlines; CellRank turns velocity into fate probabilities; PAGA gives an independent, coarse topology map. The next three sections show the illustrative code for each stage.

Estimate velocity with scVelo

scVelo (Bergen et al., 2020) is the standard velocity estimator. It works on an AnnData object (adata) that carries spliced and unspliced layers — the two matrices produced upstream by velocyto or kallisto | bustools. The workflow is four calls plus a plot: preprocess and normalize, compute the moments used by the model, estimate velocity, build the velocity graph, then draw the streamlines.

# ILLUSTRATIVE — Python, not executed here. Run in a Python environment
# with an AnnData object that has `spliced` and `unspliced` layers.
import scvelo as scv

# `adata` already has spliced/unspliced layers (from velocyto or kallisto|bustools)
# and a UMAP in adata.obsm["X_umap"], plus a cluster label in adata.obs["clusters"].

# 1. Filter genes and normalize spliced + unspliced counts
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)

# 2. First/second-order moments among nearest neighbours (needed by the model)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# 3. Estimate velocity. mode="stochastic" uses second moments; "dynamical"
#    solves the full splicing kinetics with EM (run scv.tl.recover_dynamics first).
scv.tl.velocity(adata, mode="stochastic")

# 4. Build the velocity graph (cosine similarities of velocity vs cell transitions)
scv.tl.velocity_graph(adata)

# 5. Project velocity onto the UMAP as streamlines
scv.pl.velocity_embedding_stream(adata, basis="umap", color="clusters")

Each step maps onto the biology above. filter_and_normalize keeps genes with enough spliced and unspliced signal and normalizes both layers. moments smooths the counts across each cell’s neighbours — the stochastic model needs the second moments (the variances), the deterministic model only the first. velocity fits the phase portrait per gene and reads off the per-cell velocity vector. velocity_graph computes, for every pair of neighbouring cells, whether the velocity vector points toward that neighbour — this is what lets the arrows be projected into 2-D. Finally velocity_embedding_stream draws the streamlines on your UMAP.

For transient systems (development, perturbation response), prefer the dynamical model: call scv.tl.recover_dynamics(adata) first, then scv.tl.velocity(adata, mode="dynamical"). It is slower but does not assume steady state, which is exactly the assumption that breaks during rapid differentiation.

Map fate probabilities with CellRank

Velocity streamlines are a visual summary. CellRank (Lange et al., 2022; Weiler et al., 2024) turns them into quantitative fate mapping: for each cell it computes a fate probability — the chance that the cell, following its velocity, ends up in each terminal state. That directly answers “where is each cell heading?” in a way a streamline plot cannot.

CellRank 2 is built around a clean two-part design — kernels compute a cell-to-cell transition matrix, and estimators analyze it:

# ILLUSTRATIVE — Python, not executed here.
import cellrank as cr

# A) A VelocityKernel builds a transition matrix from scVelo's velocity...
vk = cr.kernels.VelocityKernel(adata).compute_transition_matrix()

# B) ...combined with a ConnectivityKernel (expression similarity) for robustness.
ck = cr.kernels.ConnectivityKernel(adata).compute_transition_matrix()
combined_kernel = 0.8 * vk + 0.2 * ck

# C) A GPCCA estimator coarse-grains that matrix into macrostates
#    (Generalized Perron Cluster Cluster Analysis).
g = cr.estimators.GPCCA(combined_kernel)
g.compute_macrostates(n_states=6, cluster_key="clusters")

# D) Classify some macrostates as terminal (the fates), then compute,
#    for every cell, its probability of reaching each terminal state.
g.predict_terminal_states()
g.compute_fate_probabilities()
g.plot_fate_probabilities(same_plot=False)

The mental model: the VelocityKernel says “cells flow in the direction their velocity points”; the ConnectivityKernel adds “cells also move between transcriptionally similar neighbours.” Weighting them (here 0.8 / 0.2) makes the transition matrix robust to noisy velocity vectors. The GPCCA estimator then finds macrostates (metastable groups), flags the ones with no outflow as terminal states, and computes the fate probabilities by simulating random walks along the transition matrix. The result is a per-cell, per-fate probability you can plot, correlate with genes, or use to find lineage drivers. CellRank 2 generalizes this beyond velocity — the same estimator runs on a pseudotime kernel or a similarity kernel — but the velocity kernel is the classic directed-fate use case.

Draw the topology with PAGA

RNA velocity and CellRank operate at single-cell resolution. Sometimes you want the opposite: a coarse map of how the major cell groups connect — which clusters are adjacent on the manifold, and which are islands. PAGA (partition-based graph abstraction) (Wolf et al., 2019) provides exactly that. It treats each cluster as a node and draws a weighted edge between two clusters proportional to how connected their cells are, giving a simplified graph that preserves the global topology a UMAP can distort.

PAGA lives in scanpy, the general single-cell toolkit:

# ILLUSTRATIVE — Python, not executed here.
import scanpy as sc

# `adata` has a neighbourhood graph and a cluster label (here Leiden clusters).
sc.tl.paga(adata, groups="leiden")

# The abstracted graph: nodes = clusters, edge weight = connectivity
sc.pl.paga(adata, color="leiden")

# Optionally lay out a PAGA-initialized force-directed graph for a topology-faithful embedding
sc.tl.draw_graph(adata, init_pos="paga")
sc.pl.draw_graph(adata, color="leiden")

PAGA is complementary to velocity, not a competitor. It gives no direction — it is undirected connectivity — but it answers a different question cleanly: is this dataset one connected manifold or several islands, and which clusters bridge to which? A common pattern is to read the PAGA graph first to understand the coarse topology, then use velocity and CellRank to put direction and fate on top of it. Using PAGA to initialize a force-directed layout (init_pos="paga") — or a UMAP via sc.tl.umap(adata, init_pos="paga") — also produces embeddings whose global structure is more trustworthy for trajectory work.

Why this lesson is illustrative

Every other trajectory lesson on this site runs its code and commits the exact output. This one does not, and the reason is worth understanding — it is a real property of velocity analysis, not a limitation of the site.

1. The data. Velocity needs spliced and unspliced count matrices, produced by a special quantification step — velocyto or kallisto | bustools run against the reads, separating intronic (unspliced) from exonic (spliced) signal. An ordinary gene-by-cell count matrix — the kind every other lesson in this series uses — does not contain the unspliced layer and simply cannot produce velocity. Velocity-ready datasets are distributed as Python AnnData (.h5ad) objects with those two layers already computed.

2. The environment. The sanctioned R bridge, velociraptor (Bioconductor), runs scVelo from R — but it does so by provisioning a conda environment at render time through basilisk. That means a network download and a Python stack materialized during the build, which is neither offline-deterministic nor reproducible in the byte-identical way our committed render freeze requires.

3. The stochasticity. Velocity estimation, the velocity graph, and CellRank’s random-walk fate probabilities are stochastic. Even pinned to a seed, the embedded arrows and probabilities can shift across versions and platforms, so an “executed” number here would not reproduce stably.

Put together, the honest and reproducible choice is to teach the Python standard illustratively and hand execution to an environment built for it. That is not a workaround — it is how this analysis is actually done in practice: on a real Python stack, with real velocity-ready data.

The R interoperability and managed-compute path

None of this locks R users out; it just means the execution happens in Python. Two things bridge the gap.

Move your object into AnnData. If your analysis lives in a Seurat or SingleCellExperiment (SCE) object, you can convert it to the Python AnnData format that scVelo, CellRank and scanpy expect:

# ILLUSTRATIVE — shows the bridge; run where the Python stack is available.

# Option 1: zellkonverter writes a SingleCellExperiment to an .h5ad file.
library(zellkonverter)
writeH5AD(sce, file = "for_scvelo.h5ad")   # then read it in Python: adata = sc.read_h5ad(...)

# Option 2: the anndata R package + reticulate build/inspect AnnData from R.
library(anndata)
ad <- read_h5ad("for_scvelo.h5ad")

# Option 3: sceasy::convertFormat() converts Seurat <-> AnnData directly.
# sceasy::convertFormat(seurat_obj, from = "seurat", to = "anndata", out = "for_scvelo.h5ad")

Crucially, conversion moves the object, not the missing biology. If your SCE has no unspliced assay, the resulting AnnData has no unspliced layer, and velocity is still impossible — you must go back and re-quantify with velocyto / kallisto | bustools to get spliced and unspliced counts in the first place.

Run velocity in a managed-compute environment. Because velocity needs a real Python environment (scVelo, CellRank, scanpy — or velociraptor + basilisk from R) and velocity-ready data, the practical home for it is a managed-compute environment: a hosted Python session or a container with the stack preinstalled and the .h5ad data mounted. That is where you copy the snippets above, run them for real, and get the streamlines and fate probabilities. This lesson gives you the map; the execution belongs there.

Common issues

“My count matrix has no unspliced layer.” This is the single most common blocker, and it is not fixable downstream. RNA velocity is computed from the ratio of unspliced pre-mRNA to spliced mRNA, so both must be quantified from the reads. A standard gene-by-cell matrix (from Cell Ranger’s filtered output, for example) carries only the total counts. You must re-run quantification with velocyto (which adds a spliced/unspliced/ambiguous split) or kallisto | bustools in its lamanno/nucleus workflow. No unspliced layer, no velocity — check this first, before anything else.

“The velocity arrows point the wrong way / confidence is low.” Velocity can genuinely fail when its assumptions break: very low unspliced capture (common in some 3′ chemistries), a system far from steady state under the deterministic model, or genes with unusual kinetics. Diagnose with scv.tl.velocity_confidence() and inspect per-gene phase portraits (scv.pl.velocity(adata, ["GeneA", "GeneB"])) — if the portraits do not show the expected curved trajectory, the estimate is unreliable. Switching to the dynamical model (scv.tl.recover_dynamics then mode="dynamical") helps for transient states; if arrows still contradict known biology, treat velocity as inconclusive for that dataset rather than forcing an interpretation.

“I work only in R — how do I even run this?” You bridge, then execute in Python. Convert your object to AnnData (zellkonverter, anndata, or sceasy), or drive scVelo from R with velociraptor — but either way you need velocity-ready data (spliced + unspliced) and a real Python environment, so plan to run it in a managed-compute session rather than inside a lightweight R render. There is no pure-R, no-Python route to RNA velocity today.

Frequently asked questions

NoteCan I compute RNA velocity in R?

Not natively — the estimation happens in Python. The Bioconductor package velociraptor lets you call scVelo from R via its scvelo() function, but under the hood it provisions a Python environment with basilisk and converts your object with zellkonverter. So you can drive velocity from an R session, but you are still running Python, and you still need spliced and unspliced counts plus a real environment for it. There is no pure-R implementation of RNA velocity.

NoteDo I need spliced and unspliced counts for RNA velocity?

Yes — they are mandatory. RNA velocity is defined from the balance between unspliced pre-messenger RNA and mature spliced mRNA, so both must be quantified. You generate them with velocyto or kallisto | bustools, which separate intronic (unspliced) from exonic (spliced) reads. A standard gene-by-cell count matrix does not contain the unspliced layer, so it cannot produce velocity no matter what tool you point at it.

NoteRNA velocity vs pseudotime — what is the difference?

Pseudotime orders cells along a trajectory but has no inherent direction — you tell it where the start is by rooting it. RNA velocity infers the direction of change per cell from the unspliced/spliced ratio, so it points the arrow for you without a root, even when you do not know the origin. In practice they are complementary: pseudotime (from Slingshot or Monocle3) gives the ordering; velocity gives the direction; and CellRank can combine both.

NoteWhat does CellRank do that scVelo does not?

scVelo estimates velocity and draws streamlines — a visual sense of flow. CellRank turns that flow into quantitative fate mapping: it builds a transition matrix from velocity (its VelocityKernel), coarse-grains it into macrostates and terminal states (the GPCCA estimator), and computes a per-cell fate probability of reaching each terminal state. Where scVelo shows the direction, CellRank tells you, cell by cell, where each cell is most likely to end up.

NoteWhat is PAGA used for?

PAGA (partition-based graph abstraction) builds a coarse-grained map of your data: each cluster becomes a node, and edge weights encode how connected the clusters are. It is the tool for reading global topology — whether the dataset is one connected manifold or several islands, and which clusters bridge to which. It gives no direction (it is undirected), so it complements RNA velocity rather than replacing it; a common workflow reads the PAGA topology first, then adds velocity and fate on top.

Test your understanding

ImportantYour turn: plan a velocity analysis from an R object

You have a Seurat object of differentiating cells with a UMAP and Leiden clusters, produced by the standard scRNA-seq pipeline. A collaborator asks you to “add RNA velocity arrows and tell us which progenitors are heading to which fate.” Write down the plan: what you must check about the data first, how you would move the object to the Python side, which tool computes the arrows, which computes the fates, and where you would actually run it. (This is a planning exercise — there is no code to execute here.)

TipHint

Start with the data question, not the tools: does the object carry an unspliced layer? A standard Seurat object built from a Cell Ranger filtered matrix does not. Then think about the bridge to AnnData, the division of labour between scVelo (velocity/arrows) and CellRank (fate probabilities), and why execution has to leave the lightweight R render.

TipSolution (conceptual — illustrative, not executed)

1. Check the data first. A standard Seurat object has only total counts — no unspliced layer — so velocity is impossible as-is. You must re-quantify the reads with velocyto or kallisto | bustools to get spliced and unspliced matrices. This is the gating step; if it cannot be done, velocity is off the table.
2. Bridge to AnnData. Convert the (re-quantified) object to .h5ad — zellkonverter::writeH5AD() from an SCE, or sceasy::convertFormat() from Seurat — so the Python stack can read it.
3. Compute the arrows with scVelo: filter_and_normalize → moments → velocity(mode="stochastic") (or dynamical for transient states) → velocity_graph → velocity_embedding_stream.
4. Compute the fates with CellRank: a VelocityKernel (combined with a ConnectivityKernel) → GPCCA estimator → compute_macrostates → predict_terminal_states → compute_fate_probabilities. The fate probabilities answer “which progenitors head to which fate.”
5. Run it in a managed-compute environment — a hosted Python session or a container with scVelo/CellRank and the .h5ad mounted — not in the browser and not in a lightweight R render, because velocity needs a real Python environment plus velocity-ready data.

NoteQuick check: your Seurat object was built from a standard Cell Ranger filtered count matrix. Why can’t you compute RNA velocity from it directly?

A. Seurat objects are incompatible with scVelo, so no velocity is possible from any Seurat object. B. The filtered count matrix has only total (spliced) counts and no unspliced pre-mRNA layer, which RNA velocity requires — you must re-quantify with velocyto or kallisto | bustools first. C. RNA velocity needs a pseudotime, and you have not run Slingshot yet.

TipShow answer

B. RNA velocity is computed from the ratio of unspliced to spliced mRNA, and a standard filtered count matrix carries only total counts — there is no unspliced layer to work with. The fix is upstream: re-run quantification with velocyto or kallisto | bustools so both species are counted. (A is wrong — a Seurat object can be bridged to AnnData once it has the layers; C is wrong — velocity infers direction on its own and does not require a precomputed pseudotime.)

Conclusion

Pseudotime told you the order of a differentiating population; RNA velocity tells you the direction. scVelo reads that direction from the unspliced/spliced balance and draws it as streamlines, CellRank turns it into per-cell fate probabilities that answer “where is each cell heading?”, and PAGA maps the coarse topology the whole thing sits on. These are Python tools that need spliced and unspliced counts, so this lesson taught them honestly — illustratively, with nothing executed — and pointed you to the R-to-AnnData bridge and the managed-compute environment where velocity actually runs. When your data has an unspliced layer and you have a real Python stack, the snippets here are the workflow; until then, the direction they add is the reason to reach for velocity once ordering alone is not enough.

Related lessons

• Trajectory & pseudotime with Slingshot — the R method that orders cells along branching lineage curves; velocity adds the direction Slingshot cannot infer on its own. · Monocle3 trajectory analysis — learns a principal graph and pseudotime in R; pair it with velocity for directed fate. · Slingshot vs Monocle3 — choose between the two R ordering methods before you add velocity on top. · Clustering & UMAP — builds the embedding and clusters that velocity streamlines and PAGA are drawn over. · Bioinformatics — the pillar.
• Where this fits: clustering & UMAP → marker genes & annotation → pseudotime ordering (with Slingshot or Monocle3) → RNA velocity, fate mapping & topology — this lesson closes the trajectory block by adding direction and fate to the ordering.

🟢 With an AI agent

Ask Prova “I have a Seurat object of differentiating cells and I want RNA velocity — do I have the right data, and what is the exact plan to run scVelo and CellRank from an R starting point?” — it checks whether your object carries an unspliced layer, walks you through the zellkonverter / sceasy bridge to AnnData, and lays out the scVelo and CellRank 2 calls you would run in a Python environment. The runtime is the judge. Ask Prova →

Was this page helpful?

Prove you can do it. Master the whole Single-Cell RNA-Seq Analysis in R series — track your path, build projects, and earn a certificate.

Start free →

Go Pro — unlimited Prova on your own data and a verifiable certificate that proves the skill.

from $15/mo billed yearly

Go Pro →

✓ You're Pro — keep going. The runtime is the judge.

Get new R & Python lessons by email

Practical, reproducible, no spam. Unsubscribe anytime.

Double opt-in. We never share your email.

Share this pageXLinkedInRedditHN

Note

This lesson is illustrative, by design. RNA velocity, CellRank fate mapping and PAGA run in Python on spliced/unspliced count matrices, which standard scRNA-seq pipelines — and this lesson’s data — do not produce. So the code here is shown for reading, not executed: no numbers or plots are reproduced on this page. Copy the snippets into a Python (or velociraptor + basilisk) environment with velocity-ready data to run them for real. The reasoning is spelled out in why this lesson is illustrative.

References

• La Manno, G., et al. (2018). RNA velocity of single cells. Nature, 560, 494–498. https://doi.org/10.1038/s41586-018-0414-6
• Bergen, V., et al. (2020). Generalizing RNA velocity to transient cell states through dynamical modeling (scVelo). Nature Biotechnology, 38, 1408–1414. https://doi.org/10.1038/s41587-020-0591-3
• Lange, M., et al. (2022). CellRank for directed single-cell fate mapping. Nature Methods, 19, 159–170. https://doi.org/10.1038/s41592-021-01346-6
• Weiler, P., et al. (2024). CellRank 2: unified fate mapping in multiview single-cell data. Nature Methods, 21, 1196–1205. https://doi.org/10.1038/s41592-024-02303-9
• Wolf, F. A., et al. (2019). PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biology, 20, 59. https://doi.org/10.1186/s13059-019-1663-x