Bio applied trajectory analysis

scRNA-seq: Trajectory Analysis and RNA Velocity

When to Use

• Ordering cells along a continuous differentiation or activation process (pseudotime).
• Summarizing cluster-level connectivity/topology of a trajectory that is too tangled to read cell-by-cell in UMAP (PAGA).
• Estimating the direction of transcriptional change per cell from spliced/unspliced count ratios (RNA velocity).
• Deciding whether a dataset even supports trajectory inference vs. discrete cell-type calling.
• Debugging a pseudotime result that looks flipped, noisy, or dominated by batch effects.

Version Compatibility

• scanpy ≥ 1.10, anndata ≥ 0.10, Python ≥ 3.10
• scvelo ≥ 0.3 (dynamical model requires numba)
• scikit-learn ≥ 1.3, scipy ≥ 1.11

Prerequisites

• pip install scanpy scvelo anndata scikit-learn scipy
• A normalized, log-transformed expression matrix (adata.X) with PCA (adata.obsm['X_pca']) and a k-NN graph already computed (sc.pp.neighbors) — see bio-applied-scrna-preprocessing and bio-applied-single-cell-scanpy.
• Cluster labels (e.g., Leiden) if using PAGA.
• Spliced/unspliced count layers if using RNA velocity (from STARsolo --soloFeatures Gene Velocyto, Velocyto CLI, or alevin-fry).

Method Overview

Method	Tool	Input	Output
Diffusion pseudotime (DPT)	scanpy sc.tl.dpt	k-NN graph + root cell	Pseudotime per cell
PAGA	scanpy sc.tl.paga	Cluster labels + k-NN graph	Cluster connectivity graph
RNA velocity (steady-state)	scVelo mode='stochastic'	Spliced + unspliced counts	Velocity arrows
RNA velocity (dynamical)	scVelo mode='dynamical'	Spliced + unspliced counts	Velocity + kinetics

When Pseudotime Is Valid / Invalid

Valid	Invalid
Continuous developmental transitions sampled	Only start + end states, no intermediates
Trajectory visible as smooth continuum in UMAP	Discrete blobs in UMAP
Differentiation is the dominant variation source	Batch effects / cell cycle dominate PCs
≥ 50 cells in intermediate stages	Too few intermediate cells

Goal: Order cells along a trajectory and identify the root without relying on the black-box scanpy call. Approach: Build a k-NN affinity graph, symmetrize it, convert affinities to distances, and take shortest-path distance from a root cell as pseudotime (this is what sc.tl.dpt approximates internally).

import numpy as np
import scipy.sparse as sp
from sklearn.neighbors import NearestNeighbors
from scipy.sparse.csgraph import dijkstra


def compute_dpt(X_pca, root_idx, n_neighbors=20):
    """Diffusion-pseudotime approximation via Gaussian-kernel affinity + Dijkstra shortest paths.

    Parameters
    ----------
    X_pca : (n_cells, n_pcs) array of PCA coordinates.
    root_idx : int, index of the root cell (pseudotime = 0 there).
    n_neighbors : int, k-NN graph size.

    Returns
    -------
    dpt : (n_cells,) array of pseudotime values scaled to [0, 1].
    """
    knn = NearestNeighbors(n_neighbors=n_neighbors, metric='euclidean')
    knn.fit(X_pca)
    distances, indices = knn.kneighbors(X_pca)

    sigma = np.median(distances[:, 1:])  # bandwidth = median NN distance
    affinities = np.exp(-distances**2 / (2 * sigma**2))

    n = X_pca.shape[0]
    A = sp.lil_matrix((n, n))
    for i in range(n):
        for j_local, j in enumerate(indices[i]):
            A[i, j] = affinities[i, j_local]
    A = sp.csr_matrix(A)
    A = (A + A.T) / 2  # symmetrize

    # Convert affinities -> distances for shortest-path search
    dist_graph = sp.csr_matrix(1.0 / (A.toarray() + 1e-10))
    dist_graph = sp.csr_matrix(np.where(A.toarray() > 0, dist_graph.toarray(), 0))

    graph_dist = dijkstra(dist_graph, directed=False, indices=root_idx)
    graph_dist[np.isinf(graph_dist)] = graph_dist[~np.isinf(graph_dist)].max()
    dpt = (graph_dist - graph_dist.min()) / (graph_dist.max() - graph_dist.min())
    return dpt


# Root selection
# Option 1: use a known early marker (e.g., a stemness gene)
root_idx = int(np.argmax(adata.X[:, adata.var_names.get_loc('SOX2')]))
# Option 2: most extreme cell on diffusion component 1
root_idx = int(np.argmin(adata.obsm['X_diffmap'][:, 0]))

scanpy DPT Workflow (production path)

import scanpy as sc

sc.pp.neighbors(adata, n_neighbors=20, n_pcs=30)
sc.tl.diffmap(adata)
adata.uns['iroot'] = root_idx
sc.tl.dpt(adata)
# result: adata.obs['dpt_pseudotime']

# Validate against ground truth or known ordering (e.g., experimental time points)
from scipy.stats import spearmanr
corr, p = spearmanr(adata.obs['dpt_pseudotime'], adata.obs['known_stage_order'])
print(f"Spearman r = {corr:.3f}, p = {p:.2e}")

PAGA

Goal: Summarize how clusters connect when >20 clusters or per-cell UMAP is too cluttered to read a trajectory from. Approach: Compute connectivity between Leiden clusters, then optionally use it to seed a more reproducible UMAP layout.

sc.tl.paga(adata, groups='leiden')
sc.pl.paga(adata, threshold=0.03, layout='fr')

# PAGA-initialized UMAP: preserves global topology, more reproducible across seeds
sc.tl.umap(adata, init_pos='paga')

Use PAGA when: >20 clusters, trajectory topology is unclear, or per-cell UMAP is too cluttered to interpret.

RNA Velocity (scVelo)

Goal: Estimate the direction and speed of a cell's transcriptional change from the ratio of unspliced (pre-mRNA) to spliced (mature mRNA) counts. Approach: du/dt = alpha - beta*u, ds/dt = beta*u - gamma*s; a gene with more unspliced than its splicing/degradation steady-state predicts is being up-regulated (positive velocity).

import scvelo as scv

scv.pp.filter_and_normalize(adata)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# Dynamical model: slower but recovers per-gene kinetics and transient states
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis='umap')

Data requirements: spliced + unspliced count layers, ≥ 50,000 reads/cell. Obtain via STARsolo (--soloFeatures Gene Velocyto), the Velocyto CLI, or alevin-fry.

import numpy as np


def toy_velocity(unspliced, spliced, beta=0.5, gamma=0.3):
    """Steady-state RNA velocity: ds/dt ~= beta*u - gamma*s.

    Positive velocity => gene is being up-regulated (more unspliced than
    the spliced/degradation balance predicts); negative => down-regulated.
    """
    return beta * np.asarray(unspliced) - gamma * np.asarray(spliced)

Pitfalls

• Root cell choice dominates results: wrong root flips pseudotime direction. Validate using known markers or experimental time points.
• Pseudotime != real time: ordering reflects transcriptional similarity, not clock time. Do not compare pseudotime values across datasets.
• RNA velocity on fully differentiated cells is noise: velocity requires cells in transition. Stable mature types show circular/short arrows — this is correct, not a bug.
• UMAP distances are not biological distances: cells far apart in UMAP may not be transcriptionally distant. Use PC space for quantitative comparisons.
• Batch effects corrupt trajectories: correct batch effects before trajectory inference; a batch effect can look like a developmental branch.
• sc.tl.dpt requires sc.tl.diffmap first: DPT reads adata.obsm['X_diffmap']; running it before diffmap raises a KeyError.

See Also

• bio-applied-scrna-preprocessing — normalization/QC that must precede trajectory inference.
• bio-applied-single-cell-scanpy — clustering (Leiden) and neighbor graph construction used as PAGA/DPT input.
• bio-applied-cell-type-annotation — labeling clusters before/after trajectory analysis.
• bio-applied-sc-integration — batch correction, required before trusting a trajectory across samples.