RNA velocity (scVelo) — directional flow from unspliced/spliced mRNA ratios. Note: scVelo is in maintenance mode (no new features); APIs remain stable.
CellRank v2 — unified kernel API combining velocity, pseudotime, real-time, or CytoTRACE signals into a Markov chain for macrostates, terminal states, and fate probabilities. CytoTRACEKernel enables velocity-free trajectory inference.
Inputs assume a preprocessed AnnData (normalized, log-transformed, HVGs, neighbors computed). Velocity additionally requires spliced/unspliced layers from velocyto or STARsolo.
Usage
Describe the analysis in natural language:
"Compute diffusion pseudotime rooted in the HSC cluster."
"Run PAGA on leiden clusters and re-embed with force-directed layout."
"Run dynamical scVelo and plot the velocity stream."
"Use CellRank with CytoTRACEKernel to find terminal states without velocity data."
"Combine VelocityKernel and ConnectivityKernel to compute fate probabilities."
The skill produces a runnable Python script or notebook cell.
Response Format
Method selection — state which approach fits the user's data (use decision tree below)
Gotchas — warn about common failure modes specific to the chosen method
Decision Tree — Choosing a Trajectory Method
Have spliced/unspliced layers?
├─ YES → RNA velocity (scVelo dynamical) → CellRank VelocityKernel
│ Best for: steady-state or dynamical kinetics inference
│
└─ NO
├─ Have experimental time points?
│ └─ YES → CellRank RealTimeKernel
│ Best for: time-course experiments (reprogramming, perturbation)
│
└─ NO time points
├─ Want directionality (root → terminal)?
│ ├─ CytoTRACEKernel (uses gene-count gradient as differentiation proxy)
│ │ Best for: differentiation without velocity; no extra input needed
│ └─ PseudotimeKernel (uses precomputed DPT as direction)
│ Best for: when you trust your DPT root choice
│
└─ Just want cluster connectivity?
└─ PAGA (no directionality, shows which clusters connect)
Core Concepts
1. Prerequisites — spliced/unspliced counts
Velocity needs spliced and unspliced layers. Generate with velocyto or STARsolo:
import scvelo as scv
adata = sc.read_h5ad("processed.h5ad")
ldata = scv.read("sample_id.loom", cache=True)
# scv.utils.merge inner-joins on barcodes — ensure obs_names match between objects
adata = scv.utils.merge(adata, ldata)
assert"spliced"in adata.layers and"unspliced"in adata.layers
2. Diffusion pseudotime (DPT)
import numpy as np
import scanpy as sc
sc.tl.diffmap(adata)
# Root cell: pick one cell in the known starting population
root_cluster = "HSC"
adata.uns["iroot"] = np.flatnonzero(adata.obs["celltype"] == root_cluster)[0]
sc.tl.dpt(adata) # writes adata.obs['dpt_pseudotime']
sc.pl.umap(adata, color=["celltype", "dpt_pseudotime"])
3. PAGA
sc.tl.paga(adata, groups="leiden")
# threshold controls edge pruning: higher = fewer edges shown# Default 0.01 shows almost everything; 0.03-0.1 is typical for clean plots
sc.pl.paga(adata, threshold=0.03, show=False)
# PAGA-initialized layout for cleaner trajectory visualization
sc.tl.draw_graph(adata, init_pos="paga")
sc.pl.draw_graph(adata, color=["leiden", "dpt_pseudotime"], legend_loc="on data")
# Gene trends along a chosen path
sc.pl.paga_path(
adata,
nodes=["HSC", "MPP", "GMP", "Mono"],
keys=["Elane", "Mpo", "Gata1"],
)
PAGA threshold note:sc.pl.paga(threshold=...) only affects visualization (which edges to draw), not the underlying graph. The connectivity weights live in adata.uns['paga']['connectivities']. Set threshold=0 to see all edges, increase to declutter.
4. RNA velocity (scVelo)
Status: scVelo is in maintenance mode. APIs are stable but no new features are planned. For new projects without velocity data, consider CellRank's CytoTRACEKernel as a velocity-free alternative.
import scvelo as scv
scv.settings.n_jobs = 8# Preprocessing — operates on spliced/unspliced layers
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
# Dynamical mode (recommended): fits full kinetic model per gene
scv.tl.recover_dynamics(adata, n_jobs=8)
scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="celltype")
# Latent time (dynamical mode only)
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color="latent_time", color_map="gnuplot")
5. CellRank v2 — unified kernel API
CellRank v2 provides a unified kernel interface. All kernels:
Accept an AnnData and call .compute_transition_matrix()
Can be combined with + and * (weighted sum)
Feed into cr.estimators.GPCCA for macrostate decomposition
5a. VelocityKernel (requires scVelo velocity)
import cellrank as cr
from cellrank.kernels import VelocityKernel, ConnectivityKernel
vk = VelocityKernel(adata).compute_transition_matrix()
ck = ConnectivityKernel(adata).compute_transition_matrix()
kernel = 0.8 * vk + 0.2 * ck # weighted combination
g = cr.estimators.GPCCA(kernel)
g.compute_schur(n_components=20)
g.compute_macrostates(n_states=5, cluster_key="celltype")
g.predict_terminal_states()
g.compute_fate_probabilities()
g.plot_fate_probabilities(same_plot=False)
# Lineage driver genes
drivers = g.compute_lineage_drivers(lineages=["Mono"], use_raw=False)
5b. CytoTRACEKernel (velocity-free, new in v2)
Uses gene-count gradient as a proxy for differentiation potential. No velocity data needed.
from cellrank.kernels import CytoTRACEKernel, ConnectivityKernel
ctk = CytoTRACEKernel(adata).compute_transition_matrix()
ck = ConnectivityKernel(adata).compute_transition_matrix()
kernel = 0.8 * ctk + 0.2 * ck
g = cr.estimators.GPCCA(kernel)
g.compute_schur(n_components=15)
g.compute_macrostates(n_states=4, cluster_key="celltype")
g.predict_terminal_states()
g.compute_fate_probabilities()
5c. RealTimeKernel (experimental time points)
For time-course experiments where cells are collected at discrete time points.
from cellrank.kernels import RealTimeKernel
# adata.obs["day"] must contain numeric time labels (e.g., 0, 2, 4, 7)
rtk = RealTimeKernel(adata, time_key="day").compute_transition_matrix()
g = cr.estimators.GPCCA(rtk)
g.compute_macrostates(n_states=6, cluster_key="celltype")
g.predict_terminal_states()
g.compute_fate_probabilities()
5d. PseudotimeKernel (DPT-directed)
Uses precomputed pseudotime to define directionality.
from cellrank.kernels import PseudotimeKernel, ConnectivityKernel
# Requires adata.obs["dpt_pseudotime"] from sc.tl.dpt
pk = PseudotimeKernel(adata, time_key="dpt_pseudotime").compute_transition_matrix()
ck = ConnectivityKernel(adata).compute_transition_matrix()
kernel = 0.8 * pk + 0.2 * ck
g = cr.estimators.GPCCA(kernel)
g.compute_macrostates(n_states=4, cluster_key="celltype")
g.predict_terminal_states()
g.compute_fate_probabilities()
5e. Kernel combination rules
# All kernels support arithmetic combination after .compute_transition_matrix()
combined = 0.6 * vk + 0.2 * ck + 0.2 * pk # weights must sum to 1# Inspect kernel quality
combined.plot_projection(basis="umap") # visualize transition arrows
6. Spatial neighbors (squidpy)
For spatial transcriptomics, replace kNN with a spatial graph:
import squidpy as sq
sq.gr.spatial_neighbors(adata, coord_type="generic", n_neighs=6)
# Feed spatial connectivities into CellRank ConnectivityKernel
Wrong: Calling sc.tl.dpt without setting adata.uns['iroot']Right: Set the root cell index before computing DPT: adata.uns['iroot'] = np.flatnonzero(adata.obs['celltype'] == root_cluster)[0]Why: Without a root, pseudotime has no direction and results are meaningless
Wrong: Running scVelo velocity without spliced/unspliced layers in the AnnData
Right: Check adata.layers.keys() after scv.utils.merge; ensure both layers exist with matching barcodes
Why: Velocity computation requires both layers; missing layers cause errors or silent failures
Wrong: Using mode='dynamical' in scv.tl.velocity without first running scv.tl.recover_dynamicsRight: Always run scv.tl.recover_dynamics(adata) before scv.tl.velocity(adata, mode='dynamical')Why: The dynamical model needs per-gene rate parameters that are only computed by recover_dynamics
Wrong: Computing scv.tl.latent_time after running velocity in stochastic mode
Right: Use mode='dynamical' for velocity if you need latent time
Why: Latent time is only defined for the dynamical model; stochastic mode does not produce the required parameters
Wrong: Calling scv.tl.velocity without running scv.pp.moments first
Right: Run scv.pp.moments(adata, n_pcs=30, n_neighbors=30) before velocity computation
Why: Moments (Ms, Mu layers) are required inputs for velocity estimation
Wrong: Running sc.tl.paga before computing cluster labels
Right: Run sc.tl.leiden (or equivalent clustering) first; pass the cluster key via groups=Why: PAGA computes connectivity between existing groups; without clusters it has nothing to connect
Wrong: Thinking the threshold parameter in sc.pl.paga() affects the underlying graph
Right: Understand that threshold is purely visual — it controls which edges are drawn, not computed
Why: Misinterpreting this leads to incorrect claims about cluster connectivity
Wrong: Re-running sc.pp.neighbors after computing velocity graph or CellRank transition matrices
Right: If neighbors must change, recompute velocity_graph and CellRank kernels afterward
Why: Velocity graph and transition matrices depend on the neighbor graph; changing it silently invalidates them
Wrong: Using n_cells parameter in CellRank v2's compute_macrostates()Right: Use n_states — the v2 API renamed this parameter
Why:n_cells was v1 syntax and will raise an error or be silently ignored in v2
Wrong: Calling predict_terminal_states() without first inspecting macrostates
Right: Always run g.plot_macrostates() and review coarse-grained transitions before predicting terminals
Why: Terminal state prediction depends on macrostate quality; garbage macrostates produce garbage fate maps
Wrong: Forcing velocity-based analysis when velocity arrows look random or incoherent
Right: Switch to CytoTRACEKernel or PseudotimeKernel when velocity signal is absent
Why: Random velocity arrows indicate the data lacks kinetic signal; forcing velocity produces meaningless trajectories
