Trajectory Inference Methods

📋 Overview

Trajectory inference (TI) methods reconstruct the dynamic processes of cellular differentiation, development, and state transitions from single-cell data. These computational approaches allow researchers to understand how cells progress through different states over time, identify key transition points, and discover the genes that drive these changes.

🔬 Major Method Categories

1. RNA Velocity-Based Methods

These methods infer future cell states by modeling the relationship between unspliced and spliced mRNA:

2. Optimal Transport-Based Methods

These methods use optimal transport theory to match cells across conditions or timepoints:

3. Graph-Based & Pseudotime Methods

Traditional approaches that construct trajectories using dimensionality reduction and graph structures:

4. Hybrid & Multi-Method Approaches

Methods that combine multiple signals or integrate different trajectory inference approaches:

5. Specialized Methods

Purpose-built tools for specific biological questions or data types:

📊 Detailed Method Comparison

This table provides a comprehensive comparison of key features across different trajectory inference methods:

🗺️ Method Selection Guide

💡 Best Practices

General Recommendations

• Start simple: Begin with scVelo for standard scRNA-seq before trying more complex methods
• Validate predictions: Use experimental data or known biology to validate trajectory predictions
• Combine methods: Tools like CellRank 2 allow integration of multiple signals
• Consider scale: For atlas-scale data (>500K cells), use MOSCOT or Monocle 3
• Quality control: Ensure high-quality data preprocessing before trajectory inference
• Assess topology: Use Monocle 3 for complex branching or disconnected structures
• Pre-check with a potency score: Before trusting a velocity result, consider running CytoTRACE 2's no-velocity developmental-potential score as a sanity check. Low velocity_confidence does not mean the biological signal is absent — cells with low sequencing depth simply produce a weak velocity signal, but their developmental-potential score can still be meaningful. Used together: velocity tells you direction, CytoTRACE 2 tells you who is plausibly a progenitor.

When to Use RNA Velocity vs Optimal Transport vs Graph-Based Methods

RNA Velocity (scVelo, veloVI):

• Single timepoint with sufficient sequencing depth
• Interested in immediate transcriptional dynamics
• Need to identify driving genes
• Standard droplet-based scRNA-seq

Optimal Transport (MOSCOT, GENOT):

• Multiple experimental timepoints available
• Very large datasets (>200K cells)
• Spatial or multimodal data
• Need to model growth/death explicitly
• Low sequencing depth (insufficient for velocity)

Graph-Based (Monocle 3, PAGA, Slingshot):

• Complex developmental trajectories with multiple branches
• Discontinuous trajectories or convergent developmental paths (Monocle 3)
• Large-scale datasets (100K-2M+ cells with Monocle 3)
• Exploratory analysis of trajectory topology (PAGA)
• When data spans multiple timepoints with gaps
• Need for pseudotime ordering along complex paths

🛠️ Hands-On Practice

The walkthrough below takes a QC'd AnnData object end to end through a standard Python trajectory stack: neighbor graph and UMAP, Leiden clusters, PAGA for coarse topology, diffusion pseudotime for a continuous ordering, RNA velocity with scVelo, and finally CellRank to turn velocity plus connectivity into terminal states and fate probabilities. Every step is retrospective and computational — the point is to interrogate an existing dataset, not to generate new material.

Environment & packages

Scanpy, scVelo, and CellRank evolve quickly and pin overlapping dependencies, so install them into one dedicated environment rather than adding them to a general analysis env. CellRank 2 expects scvelo >= 0.3 and scanpy >= 1.9; installing all three in a single pip call lets the resolver pick a mutually consistent set.

# conda / mamba recommended
conda create -n traj python=3.10 -y
conda activate traj

# install together so the resolver reconciles shared pins
pip install "scanpy>=1.9" "scvelo>=0.3" "cellrank>=2.0" leidenalg igraph

# only needed if you must GENERATE spliced/unspliced counts yourself
# pip install velocyto scvelo[louvain] loompy
# (alternative: STARsolo --soloFeatures Gene Velocyto, or kb-python with --workflow lamanno)

Hardware. A laptop with 16 GB RAM comfortably handles ~50k cells through PAGA, DPT, and the stochastic velocity model. The scVelo dynamical model is the bottleneck: it fits a full splicing ODE per gene by EM and can take hours on 100k cells, so run it on a compute node with n_jobs set to the available cores. CellRank's kernels and GPCCA decomposition are sparse-matrix operations and stay cheap by comparison. No GPU is required anywhere in this pipeline.

Data structures & formats

• AnnData — the container everything reads and writes: adata.X (expression), adata.obs (per-cell annotation, where pseudotime and cluster labels land), adata.obsm (embeddings such as X_pca, X_umap, X_diffmap), adata.uns (graph and PAGA results)
• adata.layers["spliced"] / adata.layers["unspliced"] — hard prerequisite for RNA velocity. These are separate count matrices produced by a velocity-aware quantifier (velocyto on the CellRanger BAM, STARsolo Velocyto mode, or kb-python's La Manno workflow). They cannot be derived from a standard filtered count matrix: the spliced/unspliced split requires read-to-intron alignment, and that information is gone once you have only a gene-by-cell matrix.
• loom — velocyto's native output (.loom), loaded with scv.read() and merged onto your QC'd object with scv.utils.merge(); barcode suffixes rarely match out of the box and usually need cleaning first
• Neighbor graph — adata.obsp["connectivities"] and adata.obsp["distances"]; PAGA, DPT, and CellRank's connectivity kernel are all functions of this single graph, so its parameters propagate everywhere
• Velocity outputs — adata.layers["velocity"] (per-gene rates), adata.uns["velocity_graph"] (cell-cell transition scores); the dynamical model additionally writes adata.obs["latent_time"]
• CellRank outputs — a transition matrix on the kernel, then adata.obsm["lineages_fwd"] holding per-cell absorption probabilities toward each terminal state

Minimal code walkthrough

Start from the QC-filtered object, build the graph, run PAGA for topology, anchor a root cell using biology and compute diffusion pseudotime, then bring in the spliced/unspliced layers for velocity and hand both signals to CellRank.

import scanpy as sc
import scvelo as scv
import cellrank as cr
import numpy as np

# 1. Start from a QC'd object (doublets removed, ambient corrected).
#    Trajectory methods are extremely sensitive to junk cells: a stray
#    doublet sits "between" two clusters and looks exactly like a
#    transitional state to any graph-based method.
adata = sc.read_h5ad("adata_qc.h5ad")
adata.layers["counts"] = adata.X.copy()          # keep raw counts for later

# 2. Standard normalisation -> HVGs -> PCA
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata                                 # log-norm values for plots/markers
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata = adata[:, adata.var.highly_variable].copy()
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50, svd_solver="arpack")

# 3. Neighbour graph -> UMAP -> Leiden.
#    n_neighbors is the single most consequential knob downstream: PAGA,
#    DPT and CellRank's connectivity kernel all read this same graph.
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=1.0, key_added="leiden")

# 4. PAGA: coarse-grained topology between clusters. Read this BEFORE
#    committing to a pseudotime -- it tells you whether the data is one
#    connected continuum or several disconnected islands.
sc.tl.paga(adata, groups="leiden")
sc.pl.paga(adata, threshold=0.05, show=False)     # prune weak edges
sc.tl.umap(adata, init_pos="paga")                # PAGA-initialised embedding

# 5. Diffusion pseudotime. The root must come from BIOLOGY, not convenience:
#    pick the cluster expressing your known progenitor/stem markers, then
#    take an extreme cell within it along the first diffusion component.
sc.tl.diffmap(adata)
root_cluster = "3"                                 # e.g. HSC/progenitor cluster
mask = (adata.obs["leiden"] == root_cluster).values
dc1 = adata.obsm["X_diffmap"][:, 1]
adata.uns["iroot"] = np.flatnonzero(mask)[np.argmin(dc1[mask])]
sc.tl.dpt(adata, n_dcs=15)                         # writes adata.obs["dpt_pseudotime"]

# Sanity check: pseudotime should INCREASE away from the progenitors.
# If your known mature markers sit at low pseudotime, the root is wrong
# and the whole trajectory is silently reversed.
print(adata.obs.groupby("leiden")["dpt_pseudotime"].median())

# 6. RNA velocity. This step REQUIRES spliced/unspliced layers -- they come
#    from velocyto/STARsolo/kb-python, and cannot be recovered from a plain
#    filtered count matrix. Merge the loom onto the QC'd object by barcode.
ldata = scv.read("sample.loom", cache=True)
ldata.var_names_make_unique()
adata = scv.utils.merge(adata, ldata)              # intersects on barcodes
assert "spliced" in adata.layers and "unspliced" in adata.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)    # kNN-smoothed first/second moments

# stochastic: seconds-to-minutes, assumes a common steady state
scv.tl.velocity(adata, mode="stochastic")

# dynamical: fits the full splicing ODE per gene by EM. Much slower
# (hours at 100k cells) but relaxes the steady-state assumption and is the
# only mode that yields latent time.
# scv.tl.recover_dynamics(adata, n_jobs=8)
# scv.tl.velocity(adata, mode="dynamical")
# scv.tl.latent_time(adata)                        # adata.obs["latent_time"]

scv.tl.velocity_graph(adata, n_jobs=8)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="leiden")

# Confidence diagnostics -- low values mean the arrows are not trustworthy.
scv.tl.velocity_confidence(adata)
scv.pl.scatter(adata, c=["velocity_length", "velocity_confidence"], cmap="coolwarm")

# 7. CellRank: combine a directional signal (velocity) with a similarity
#    signal (connectivity). The connectivity kernel regularises velocity,
#    which is noisy at the single-cell level.
vk = cr.kernels.VelocityKernel(adata).compute_transition_matrix()
ck = cr.kernels.ConnectivityKernel(adata).compute_transition_matrix()
combined = 0.8 * vk + 0.2 * ck                     # weights are a real choice

g = cr.estimators.GPCCA(combined)
g.compute_schur(n_components=20)
g.compute_macrostates(n_states=6, cluster_key="leiden")
g.predict_terminal_states()                        # inspect against known endpoints
g.compute_fate_probabilities()                     # adata.obsm["lineages_fwd"]
g.plot_fate_probabilities(same_plot=False)

# 8. Genes correlated with commitment to one fate
drivers = g.compute_lineage_drivers(lineages=["Erythroid"])
print(drivers.head(20))

adata.write_h5ad("adata_trajectory.h5ad")

Common pitfalls & tips

• Pseudotime is a rank ordering, not real time. DPT and latent time produce a monotone coordinate along the graph, not hours or days. Two cells one pseudotime unit apart may be minutes or weeks apart biologically, and pseudotime is not comparable across datasets. Never report it in physical units or use differences in it as rates.
• A wrong root silently reverses the trajectory. sc.tl.dpt will happily run from any iroot and return a smooth, plausible-looking ordering. Choose the root from progenitor markers, then verify that known mature populations land at high pseudotime — nothing in the output warns you when the direction is flipped.
• Ask whether a trajectory is warranted at all. These methods fit a continuum to whatever they are given, including discrete, terminally differentiated populations with no transitions between them. Check PAGA connectivity first: sparse, near-zero edges between clusters mean you are looking at distinct cell types, and any pseudotime you compute across them is an artifact.
• Velocity fails when kinetics vary or genes are off steady state. The steady-state (stochastic) model assumes one global splicing/degradation ratio per gene; if a gene has different kinetics across cell types, or the system was never at equilibrium, the arrows can point confidently backwards. Mature neurons and other slow-turnover populations are classic failure cases. Always inspect velocity_confidence and check a few phase portraits with scv.pl.velocity() instead of trusting the streamline plot alone.
• Smoothing and kNN parameters reshape the topology. n_neighbors in sc.pp.neighbors and scv.pp.moments, the number of PCs, and the PAGA edge threshold all change which branches appear. Small neighborhoods fragment a real continuum; large ones bridge unrelated populations. Re-run the pipeline across a parameter grid and only report branches that survive.
• Batch effects manufacture false branches. Two samples that separate in PCA become two "lineages" in any graph-based method. Integrate (Harmony, scVI, BBKNN) before building the neighbor graph, and colour the final trajectory by batch or orig.ident — if a branch is one sample, it is technical.
• Validate branch assignments with known markers. Terminal states from GPCCA and fate probabilities are model outputs, not ground truth. Score each predicted endpoint with an established lineage marker set before naming it; if a terminal state has no marker support, treat it as a numerical artifact of the macrostate count you chose.