/trajectory-inference

name: "trajectory-inference"
description: "Trajectory inference and cell fate mapping for single-cell and spatial transcriptomics: PAGA graph + diffusion pseudotime (no velocity needed), scVelo RNA velocity (stochastic / dynamical, from spliced/unspliced counts), CellRank 2 terminal states, fate probabilities and lineage drivers, then pseudotime and fates mapped on tissue with a spatial autocorrelation check. Outputs pseudotime, velocity streams, fate probability tables. Imaging panels (Xenium, MERFISH) have no unspliced counts: use PAGA/DPT or the CellRank pseudotime kernel. For a basic PAGA recipe inside a standard scRNA-seq workflow see scanpy-scrna-seq."
license: "CC-BY-4.0"

Trajectory Inference

Overview

Order cells along differentiation or activation paths and estimate where each cell is heading. Output: pseudotime, velocity streams, terminal states and fate probabilities, and on spatial data their maps on the tissue with a check of whether the ordering follows space. Three layers, each optional on top of the previous one:

| Method | Needs | Gives | |---|---|---| | PAGA + DPT | Clusters, neighbour graph, a chosen root cluster | Cluster connectivity graph, pseudotime | | scVelo | `layers["spliced"]` and `layers["unspliced"]` | Direction of change per cell (velocity) | | CellRank 2 | A transition kernel (velocity or pseudotime) | Terminal states, fate probabilities, driver genes |

Spliced/unspliced counts come from velocyto, STARsolo (`--soloFeatures Gene Velocyto`) or kallisto|bustools run on the reads; for Visium and Slide-seq run them on the Space Ranger / pipeline BAM. Imaging platforms (Xenium, MERFISH, CosMx) count probes, not introns, so velocity is not available there.

When to Use

• Ordering cells along a differentiation or activation path from a chosen start
• Estimating where cells are heading from spliced/unspliced counts (RNA velocity)
• Finding terminal states, fate probabilities and the genes driving each fate
• Checking whether a developmental or activation ordering is laid out in tissue space
• Finding genes that change along the inferred path
• Use `scanpy-scrna-seq` instead when you only need a quick PAGA graph inside a standard scRNA-seq workflow
• For tissue regions rather than cell-state paths, use `spatial-domain-detection`

Prerequisites

• **Input**: raw counts plus a cluster or cell type column; for velocity, the two layers above
• **Spatial data**: coordinates in `obsm["spatial"]` for Step 6
• **Python packages**: `scvelo`, `cellrank` (CellRank 2.x needs Python >= 3.12; CPU is fine)

# Install only what `import` reports missing. `--pypi` keeps the env's scanpy/numba
# stack; a conda add can swap it and bare pip can shadow it.
pixi add --pypi scvelo cellrank scanpy

Workflow

1. Load data and check layers 2. Preprocess and build the neighbour graph 3. PAGA and diffusion pseudotime 4. RNA velocity (scVelo) 5. Terminal states and fate probabilities (CellRank) 6. Map pseudotime and fates on the tissue

Step 1: Load Data

import numpy as np
import pandas as pd
import scanpy as sc

adata = sc.read_h5ad("path/to/data.h5ad")        # your data
# Demo with spliced/unspliced layers (mouse pancreas endocrinogenesis, downloaded to ./data on first use):
# import scvelo as scv; adata = scv.datasets.pancreas()
CLUSTER_KEY = "clusters"                          # cell type / cluster column

has_velocity_layers = {"spliced", "unspliced"} <= set(adata.layers)
print(adata)
print(f"velocity layers: {has_velocity_layers}; spatial: {'spatial' in adata.obsm}")
print(adata.obs[CLUSTER_KEY].value_counts())

Step 2: Preprocess and Neighbour Graph

Uses the scVelo order (filter on spliced+unspliced, normalise all layers, then log `X`). Without velocity layers, the scanpy lines alone apply.

import scvelo as scv

if has_velocity_layers:
    scv.pp.filter_genes(adata, min_shared_counts=20)
    scv.pp.normalize_per_cell(adata)              # normalises X, spliced, unspliced
else:
    adata.layers["counts"] = adata.X.copy()
    sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)     # imaging panels: skip, use all genes
sc.pp.pca(adata, n_comps=30)
sc.pp.neighbors(adata, n_pcs=30, n_neighbors=30, random_state=0)
if "X_umap" not in adata.obsm:
    sc.tl.umap(adata)
adata.obs[CLUSTER_KEY] = adata.obs[CLUSTER_KEY].astype("category")

Step 3: PAGA and Diffusion Pseudotime

The root cluster is a biological choice (progenitor, naive state); ask the user when it is not obvious from markers.

ROOT_GROUP = "Ngn3 low EP"                        # root cluster name

sc.tl.paga(adata, groups=CLUSTER_KEY)
sc.pl.paga(adata, threshold=0.05, show=False, save="_graph.png")
groups = adata.obs[CLUSTER_KEY].cat.categories
paga_adj = pd.DataFrame(adata.uns["paga"]["connectivities"].toarray(), index=groups, columns=groups)
paga_adj.to_csv("paga_adjacency.csv")

# Root = the root-cluster cell at the extreme of diffusion component 1
# (the first cell of the cluster is arbitrary and shifts the whole pseudotime)
sc.tl.diffmap(adata)
dc1 = adata.obsm["X_diffmap"][:, 1]
in_root = (adata.obs[CLUSTER_KEY] == ROOT_GROUP).to_numpy()
assert in_root.any(), f"{ROOT_GROUP} not in {list(groups)}"
side = np.sign(dc1[in_root].mean() - dc1.mean()) or 1.0
adata.uns["iroot"] = int(np.flatnonzero(in_root)[np.argmax(side * dc1[in_root])])
sc.tl.dpt(adata)
print(adata.obs.groupby(CLUSTER_KEY, observed=True)["dpt_pseudotime"].median().sort_values())
sc.pl.umap(adata, color=[CLUSTER_KEY, "dpt_pseudotime"], show=False, save="_pseudotime.png")

If the median pseudotime per cluster does not follow the expected order, the root or the clustering is wrong; DPT is not a test of the order.

Step 4: RNA Velocity (scVelo)

MODE = "stochastic"                               # "dynamical" is slower but gives latent time

scv.pp.moments(adata, n_pcs=None, n_neighbors=None)      # reuses the Step