Bio Single Cell Trajectory Inference

Version Compatibility

Reference examples tested with: scanpy 1.10+, scVelo 0.3+, CellRank 2.0+, Monocle3 1.3+, Slingshot 2.x

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Trajectory Inference

"Find the developmental trajectory in my data" -> Order cells along a continuous manifold and assign a pseudotime, fate probabilities, or velocity direction.

• Python: sc.tl.paga/sc.tl.dpt (scanpy), palantir, scvelo, cellrank (kernels + GPCCA)
• R: slingshot, monocle3, tradeSeq

Governing Principle

A snapshot is not a movie. Every method substitutes transcriptomic similarity for temporal adjacency: it takes a single static sample and imposes an ordering or a Markov process on the kNN graph. The output is meaningful only when the population is a genuine continuum of asynchronously-progressing cells, sampled densely enough to bridge intermediate states. Five rules follow and they drive every downstream decision.

1. Pseudotime is geometry, not a clock. The axis is monotone in progression at best; equal pseudotime intervals are NOT equal real-time intervals, and rates differ across lineages. Reading an interval as a duration is a category error.
2. The existence of a continuum is a BIOLOGICAL judgment made BEFORE ordering. Any method returns numbers on discrete cell types or on noise; no algorithm tests whether a continuum exists. Decide topology first with PAGA connectivity, then order.
3. Root-cell choice flips every gene trend. Pseudotime is defined only up to an origin; the sign of every trend and which cells are "early" invert with the root. Anchor the root with orthogonal evidence (known marker, real sampling time, velocity, or a stemness score), never by eye on a UMAP, which distorts global geometry. When no progenitor marker and no real timepoints are available and velocity is invalid (mature or non-cycling tissue), fall back to a model-free stemness score (CytoTRACE/CytoTRACE2) to nominate the root, and treat the resulting direction as a hypothesis rather than an established origin.
4. Near bifurcations use fate PROBABILITIES, not hard branch labels. A multipotent progenitor's fate is genuinely undetermined, so a hard assignment to one lineage is biologically false. Represent each cell as a distribution over terminal fates (Palantir, CellRank).
5. Transcriptomic state does not fully predict fate (Weinreb 2020). Sister cells in an indistinguishable state systematically diverge in fate, so a state-based branch call can be systematically wrong at exactly the decision point it most wants to resolve. Frame fate calls as predictions, not measurements, and validate with orthogonal evidence.

Snapshot dynamics are formally non-identifiable (Weinreb 2018, gauge freedom): a single snapshot constrains a family of dynamics, not one, and a unique answer requires extra assumptions (a potential field, known birth-death rates, or real timepoints). Treat any single-method pseudotime as a hypothesis to cross-validate, never a measurement.

Method Decision Table

Saelens 2019 benchmarked 45 methods across 110 real + 229 synthetic datasets: no single method wins across all topologies, so selection is topology-first.

Choosing by topology: linear -> Slingshot/DPT; bifurcating/multifurcating -> Slingshot/Palantir/CellRank 2; tree (>2 branches) -> Slingshot/Monocle3/PAGA-Tree; cyclic -> PAGA (Slingshot/Monocle2 cannot); disconnected/unknown -> PAGA first, then per-component pseudotime. Methodology evolves; verify current best practice against installed package docs before committing to one method, and report multi-method concordance.

Decide Topology First With PAGA

Goal: Test whether putative branches are truly connected before committing to a continuous model. Approach: Partition the kNN graph, compute PAGA connectivity, prune weak edges by threshold, then seed a global-faithful UMAP from PAGA.

import scanpy as sc, numpy as np
sc.pp.neighbors(adata, n_neighbors=15, use_rep='X_pca')
sc.tl.leiden(adata, resolution=1.0, flavor='igraph', n_iterations=2, directed=False)
sc.tl.paga(adata, groups='leiden')
sc.pl.paga(adata, threshold=0.03, color='leiden')   # prune low-connectivity (likely spurious) edges
sc.tl.umap(adata, init_pos='paga')                  # global topology preserved, local detail kept

The threshold in sc.pl.paga is the key judgment call: isolated clusters with no surviving edges are discrete cell types, not trajectory branches, and must not be forced into one ordering.

Diffusion Pseudotime From an Anchored Root

Goal: Assign a scalar pseudotime once topology is fixed. Approach: Run a diffusion map, set the root as a positional index on adata.uns['iroot'], then run DPT.

sc.tl.diffmap(adata, n_comps=15)
adata.uns['iroot'] = np.flatnonzero(adata.obs['cell_type'] == 'HSC')[0]   # root anchored by a known marker, not by eye
sc.tl.dpt(adata, n_dcs=10, n_branchings=0)          # n_branchings=0 -> pure pseudotime; branch mode is fragile

iroot is a positional integer into adata.obs_names, set on adata.uns BEFORE dpt. The entire ordering and the sign of every gene trend flip with this choice.

Fate Probabilities With Palantir

Goal: Represent each cell as a distribution over terminal fates near bifurcations. Approach: Build a diffusion-map multiscale space, run Palantir from an early cell, and read pseudotime, entropy, and branch probabilities.

import palantir
dm_res = palantir.utils.run_diffusion_maps(adata, n_components=5)
ms_data = palantir.utils.determine_multiscale_space(dm_res)
pr_res = palantir.core.run_palantir(ms_data, early_cell='HSC_cell_id', terminal_states=None, num_waypoints=1200)
# pr_res.pseudotime, pr_res.entropy (differentiation potential), pr_res.branch_probs

Entropy of the fate-probability vector is the differentiation-potential proxy: high near multipotent cells, falling toward 0 as cells commit. Auto terminal-state detection can miss real fates or invent spurious ones, so verify terminals against markers.

Directed Fate Mapping With CellRank 2

Goal: Infer initial states, terminal states, and uncertainty-aware fate probabilities from any directional evidence source. Approach: Build a directed transition matrix from one or more kernels, combine with a connectivity kernel for smoothing, then coarse-grain into macrostates with GPCCA.

import cellrank as cr
pk = cr.kernels.PseudotimeKernel(adata, time_key='dpt_pseudotime').compute_transition_matrix()
ck = cr.kernels.ConnectivityKernel(adata).compute_transition_matrix()
combined = 0.8 * pk + 0.2 * ck                      # weights are a researcher choice; sweep them

g = cr.estimators.GPCCA(combined)
g.compute_macrostates(n_states=10, cluster_key='leiden')   # n_states from the Schur/eigenvalue spectral gap
g.predict_terminal_states(method='stability')
g.predict_initial_states(n_states=1)
g.compute_fate_probabilities()
g.compute_lineage_drivers()

Kernels decouple WHERE direction comes from (RealTime when timepoints exist, Pseudotime/CytoTRACE otherwise, Velocity only when trustworthy, Connectivity for smoothing) from WHAT is computed (GPCCA macrostates + fate probabilities). Prefer the RealTimeKernel for time courses. Fate probabilities are a deterministic function of the transition matrix, so a wrong kernel yields confidently wrong, well-formed probabilities with no internal warning; check that conclusions survive dropping the velocity kernel.

Slingshot and Monocle3 (R)

Goal: Fit smooth lineage curves (Slingshot) or a principal graph (Monocle3) and order cells. Approach: Slingshot needs user-supplied dimred + cluster labels + a start cluster; Monocle3 learns its own graph and roots by node.

library(slingshot)
sce <- slingshot(sce, clusterLabels='seurat_clusters', reducedDim='UMAP', start.clus='HSC')
pt  <- slingPseudotime(sce)     # cells x lineages; NA off-lineage; a trunk cell scores in every descendant lineage

library(monocle3)
cds <- cluster_cells(cds)                           # produces clusters AND partitions
cds <- learn_graph(cds, use_partition = TRUE)       # TRUE allows disconnected trajectories
cds <- order_cells(cds, root_pr_nodes = root_node)  # root via graph node name, anchored by biology
graph_test_res <- graph_test(cds, neighbor_graph = 'principal_graph', cores = 4)   # Moran's I trajectory DE

start.clus is mandatory in practice for Slingshot; downstream DE goes through tradeSeq (fitGAM then associationTest for any-variation-along-pseudotime or startVsEndTest for endpoint contrasts), not Slingshot itself. Monocle3's own trajectory DE is graph_test above. Monocle3's principal graph is learned in UMAP space, so loops and branches can be embedding artifacts.