Imported from tony-zhelonkin/scio (skills/rna-velocity-trajectory/SKILL.md). Install upstream with npx skills add tony-zhelonkin/scio --skill rna-velocity-trajectory. Copyright stays with the author (MIT).

RNA Velocity & Trajectory Inference: Cellular Dynamics from Splicing Kinetics

Foundation: RNA velocity infers cellular differentiation trajectories and future cell states by modeling the ratio of unspliced (nascent) to spliced (mature) RNA. The change in spliced RNA over time—termed RNA velocity—predicts where cells are heading in gene expression space.

When to Use RNA Velocity

Infer developmental trajectories from single snapshot scRNA-seq data
Predict terminal cell fates and differentiation endpoints
Identify driver genes for cell state transitions
Validate lineage relationships in differentiation systems
Analyze cell cycle dynamics and proliferation states
Study transient cell populations in development or disease

Limitations:

Requires alignment to capture intronic reads (spliced/unspliced quantification)
Time scale of process must match RNA half-life (~hours)
Fails in steady-state systems (e.g., mature PBMCs with no transitions)
Phase portraits must show expected almond shape for reliable inference
Low-dimensional projections can be misleading; use CellRank for rigorous analysis

Quick Start: scVelo

import scvelo as scv
import scanpy as sc

# Load data with spliced/unspliced layers
adata = scv.read_loom("velocyto_output.loom")
# Or: adata.layers["spliced"], adata.layers["unspliced"] from kallisto|bustools

# Preprocessing
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# Compute velocity (dynamical model recommended)
scv.tl.recover_dynamics(adata, n_jobs=8)  # Infer kinetic rates
scv.tl.velocity(adata, mode="dynamical")

# Build velocity graph and visualize
scv.tl.velocity_graph(adata, n_jobs=8)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="clusters")

Velocity Models Comparison

Model	Command	Assumptions	Use When
Steady-state	`mode="deterministic"`	Equilibrium reached, common splicing rate	Quick exploration, simple systems
EM/Dynamical	`mode="dynamical"` after `recover_dynamics()`	Gene-specific rates, full kinetics	Publication-quality, complex trajectories
Stochastic	`mode="stochastic"`	Adds noise estimation	Uncertainty quantification

Dynamical model parameters:

scv.tl.recover_dynamics(
    adata,
    n_jobs=8,                    # Parallel processing
    fit_scaling=True,            # Gene-wise scaling
    var_names='velocity_genes',  # Subset to velocity genes
    max_iter=100                 # EM iterations
)
# Outputs: adata.var["fit_alpha", "fit_beta", "fit_gamma", "fit_t_", "fit_likelihood"]

scVelo AnnData Output Fields

After running the scVelo workflow, the following fields are added to the AnnData object:

Location	Key	Description
`adata.layers`	`velocity`	RNA velocity per gene per cell
`adata.layers`	`fit_t`	Fitted latent time per gene per cell (dynamical only)
`adata.obsm`	`velocity_umap`	2D velocity vectors projected onto UMAP
`adata.obs`	`velocity_pseudotime`	Pseudotime inferred from velocity graph
`adata.obs`	`latent_time`	Shared latent time from dynamical model
`adata.obs`	`velocity_length`	Speed of each cell (magnitude)
`adata.obs`	`velocity_confidence`	Coherence of velocity direction
`adata.var`	`fit_likelihood`	Gene-level model fit quality (dynamical only)
`adata.var`	`fit_alpha`	Transcription rate per gene
`adata.var`	`fit_beta`	Splicing rate per gene
`adata.var`	`fit_gamma`	Degradation rate per gene
`adata.uns`	`velocity_graph`	Cell-cell transition probability matrix

Quick inspection:

# After running velocity workflow
print([k for k in adata.obs.columns if 'velocity' in k or 'latent' in k])
print([k for k in adata.layers.keys()])
print([k for k in adata.var.columns if 'fit_' in k])

VeloVI: Variational Inference for Velocity

VeloVI provides uncertainty quantification via a variational autoencoder framework.

from velovi import VELOVI, preprocess_data
import scvelo as scv

# Load and preprocess
adata = scv.datasets.pancreas()
scv.pp.filter_and_normalize(adata, min_shared_counts=30, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
adata = preprocess_data(adata)  # Runs scVelo velocity internally

# Setup and train
VELOVI.setup_anndata(adata, spliced_layer="Ms", unspliced_layer="Mu")
vae = VELOVI(adata)
vae.train()  # ~5-10 min, early stopping enabled

# Extract outputs
latent_time = vae.get_latent_time(n_samples=25)
velocities = vae.get_velocity(n_samples=25, velo_statistic="mean")
adata.layers["velocity"] = velocities

# Uncertainty quantification
uncertainty_df, _ = vae.get_directional_uncertainty(n_samples=100)
adata.obs["velocity_uncertainty"] = uncertainty_df["directional_cosine_sim_variance"].values

# Gene-level confidence via permutation
perm_df, _ = vae.get_permutation_scores(labels_key="clusters")
adata.var["permutation_score"] = perm_df.max(1).values

VeloVI advantages:

Probabilistic framework with uncertainty estimates
Intrinsic uncertainty (per-cell) and extrinsic uncertainty (model stability)
Permutation scores identify genes with reliable velocity signals
Better handling of noisy data

Chronocell: Biophysical Trajectory Modeling

Chronocell infers trajectories using explicit kinetic equations and topology specification. It provides pseudotime, state assignments, and kinetic parameter estimates.

Installation:

import sys
sys.path.insert(0, '/path/to/FGP_2024')  # Local copy
# Or: pip install git+https://github.com/pachterlab/Chronocell.git

Basic usage:

from Chronocell.inference import Trajectory
from Chronocell.utils import select_genes_by_G
from Chronocell.plotting import plot_t, plot_y
import numpy as np

# Prepare data (n_cells x n_genes x 2 array)
X = np.stack([
    adata.layers["unspliced"].toarray(),
    adata.layers["spliced"].toarray()
], axis=-1)

# Select informative genes via G-test
gene_mask = select_genes_by_G(adata, model="Poisson", n_genes=50)
X_select = X[:, gene_mask, :]

# Define topology: linear trajectory with 3 states
topo = np.array([[0, 1, 2]])           # State 0 → 1 → 2
tau = np.array([0.0, 1.0, 2.0, 3.0])   # Transition timepoints

# Fit trajectory
traj = Trajectory(topo, tau, model="two_species_ss", verbose=1)
traj = traj.fit(
    X_select,
    n_init=20,              # Multiple random restarts (critical)
    epoch=100,              # EM iterations
    params={'r': adata.obs["total_counts"].values},  # Read depth
    parallel=True,
    n_threads=8
)

# Extract results
cell_pseudotime = (traj.Q.sum(1) @ traj.t)  # Weighted average time
cell_states = traj.Q.sum(2).argmax(1)       # Most likely state
adata.obs["chronocell_time"] = cell_pseudotime

# Model comparison
print(f"ELBO: {traj.compute_lower_bound(X_select):.2f}")
print(f"AIC: {traj.compute_AIC(X_select):.2f}")
print(f"BIC: {traj.compute_BIC(X_select):.2f}")

# Visualize
plot_t(traj)                    # Cell-time-state posterior heatmap
plot_y(traj, idx=np.arange(5))  # Gene expression trajectories

Chronocell models:

Model	Description	Use Case
`one_species`	Single RNA species (spliced only)	Simplified analysis
`two_species_ss`	Two species, steady-state splicing	Standard velocity data
`two_species_ss_tau`	Gene-wise time shifts	Non-steady-state genes

Topology examples:

# Linear: 0 → 1 → 2
topo = np.array([[0, 1, 2]])

# Bifurcating: 0 → 1 → {2, 3}
topo = np.array([[0, 1, 2],
                 [0, 1, 3]])

# Merging: {0, 1} → 2
topo = np.array([[0, 2],
                 [1, 2]])

CellRank Integration

CellRank uses velocity to compute fate probabilities and terminal states, avoiding misleading 2D projections.

import cellrank as cr

# Initialize kernel from velocity
vk = cr.kernels.VelocityKernel(adata)
vk.compute_transition_matrix()

# Combine with connectivity kernel for robustness
ck = cr.kernels.ConnectivityKernel(adata)
ck.compute_transition_matrix()
combined_kernel = 0.8 * vk + 0.2 * ck

# Compute fate probabilities
g = cr.estimators.GPCCA(combined_kernel)
g.compute_macrostates(n_states=5)
g.set_terminal_states_from_macrostates()
g.compute_fate_probabilities()

# Visualize
cr.pl.fate_probabilities(adata)
g.compute_lineage_drivers()

Data Preparation

Spliced/Unspliced Quantification

Option 1: velocyto (BAM-based)

velocyto run10x /path/to/cellranger_output gtf_file.gtf
# Output: velocyto_output.loom with spliced/unspliced layers

Option 2: kallisto|bustools (Pseudoalignment, faster)

kb count -i index.idx -g t2g.txt -x 10xv3 \
  --workflow lamanno --loom \
  R1.fastq.gz R2.fastq.gz

Option 3: alevin-fry (Memory efficient)

alevin-fry quant --usa  # Unspliced/spliced/ambiguous mode

Loading into AnnData

# From loom file
adata = scv.read_loom("output.loom")

# From separate matrices
adata = sc.read_h5ad("processed.h5ad")
adata.layers["spliced"] = spliced_matrix
adata.layers["unspliced"] = unspliced_matrix

# Verify layers exist
assert "spliced" in adata.layers and "unspliced" in adata.layers

Downstream Analysis

Latent Time

scv.tl.latent_time(adata)
scv.pl.scatter(adata, color="latent_time", cmap="gnuplot")

Velocity Confidence & Genes

# Confidence per cell
scv.tl.velocity_confidence(adata)
scv.pl.scatter(adata, color="velocity_confidence")

# Top velocity genes
scv.tl.rank_velocity_genes(adata, groupby="clusters", min_corr=0.3)
top_genes = adata.uns["rank_velocity_genes"]["names"]["cluster_1"][:10]
scv.pl.scatter(adata, basis=top_genes[:3], color="clusters")

Phase Portraits

# Check model fit quality
scv.pl.velocity(adata, var_names=["Gene1", "Gene2", "Gene3"], color="clusters")

Critical Gotchas

Issue	Solution
No almond-shaped phase portrait	Gene violates kinetic assumptions; exclude from analysis
Backflow in velocity stream	Check if steady-state assumption violated; use dynamical model
All velocities point same direction	Insufficient dynamic range; check if system is at steady state
Low spliced/unspliced ratio	Poor intronic capture; check alignment to introns
Conflicting velocity directions	Use CellRank for robust fate inference instead of 2D projection
Missing `layers["spliced"]`	Need velocity-aware quantification (velocyto, kb count --lamanno)
Chronocell poor convergence	Increase `n_init` (20+), check topology matches biology
VeloVI training unstable	Reduce learning rate, increase n_epochs for warmup

Model Selection Guidelines

When to use each tool:

Scenario	Recommended Tool
Quick exploration	scVelo steady-state
Publication-quality velocity	scVelo dynamical
Uncertainty quantification	VeloVI
Explicit topology/pseudotime	Chronocell
Fate probability, terminal states	CellRank
Cross-method validation	Run multiple, compare

Velocity gene selection:

# scVelo: automatic during recover_dynamics
# Manual: filter by fit likelihood
good_genes = adata.var["fit_likelihood"] > 0.1

# VeloVI: permutation score
good_genes = adata.var["permutation_score"] > 0.05

Visualization Best Practices

# Stream plot (most common, but can mislead)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="clusters")

# Arrow plot (more precise direction)
scv.pl.velocity_embedding(adata, basis="umap", arrow_length=3, arrow_size=2)

# Grid plot (summarized flow)
scv.pl.velocity_embedding_grid(adata, basis="umap", color="clusters")

# Phase portraits for top genes (model validation)
top_genes = adata.var["fit_likelihood"].sort_values(ascending=False).index[:6]
scv.pl.scatter(adata, basis=top_genes, color="clusters", frameon=False)

# Latent time heatmap
scv.pl.heatmap(adata, var_names=top_genes, sortby="latent_time", col_color="clusters")

Resources

scVelo:

Documentation: https://scvelo.readthedocs.io
Paper: Bergen et al. 2020, Nature Biotechnology
Tutorial: https://scvelo.readthedocs.io/en/stable/DynamicalModeling/

VeloVI:

Documentation: https://docs.scvi-tools.org/en/stable/user_guide/models/velovi.html
Paper: Gayoso et al. 2022, Nature Biotechnology
Tutorial: https://docs.scvi-tools.org/en/stable/tutorials/notebooks/scrna/velovi.html

Chronocell:

GitHub: https://github.com/pachterlab/Chronocell
Paper: Fang, Gorin & Pachter 2025, PLoS Computational Biology
Local: /01_Scripts/FGP_2024/Chronocell/

CellRank:

Documentation: https://cellrank.readthedocs.io
Paper: Lange et al. 2022, Nature Methods

RNA velocity theory:

Original paper: La Manno et al. 2018, Nature
Review: Bergen et al. 2021, Molecular Systems Biology

When not to use

Do not use on steady-state systems (e.g. mature PBMCs) — velocity assumptions fail.
Do not use without spliced/unspliced layers — requires intron-aware quantification upstream.

rna-velocity-trajectory

About