# Bio Applied Trajectory Analysis

> Order scRNA-seq cells with diffusion pseudotime (scanpy sc.tl.dpt), PAGA graphs, and RNA velocity (scVelo) on spliced/unspliced counts. Use for pseudotime, root cell selection, or RNA velocity streamline plots.

- Skill: `pavel-kravchenko/bio-applied-trajectory-analysis` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-trajectory-analysis`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-trajectory-analysis/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: AI & ML
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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).

```python
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)

```python
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.

```python
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).

```python
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.

```python
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.