Trajectory Inference & Pseudotime Analysis

When to Use

• Analyzing developmental processes, differentiation, or cell state transitions
• Ordering cells along a temporal trajectory
• Identifying branch points and cell fate decisions
• Finding genes that change along pseudotime

Do NOT use when you need to:

• 数据没有连续分化过程 → use scrna-clustering-umap for discrete populations
• 需要 RNA velocity → use scVelo/dynamo (not covered by this skill)
• 数据未聚类 → (use scrna-clustering-umap first)
• 需要伪时间基因动力学 → (use scrna-gene-regulatory-network for regulon activity)

Prerequisites

• Python >= 3.10
• Required packages: scanpy, anndata, numpy, pandas, matplotlib, scipy; optional: palantir
• Input: h5ad file with PCA embeddings in obsm["X_pca"] and neighbor graph in uns["neighbors"]

Core Rules

1. Always verify the data contains a continuous biological process before running trajectory inference — discrete cell types without transitions do not need pseudotime.
2. Root cell selection is critical: use known stem/progenitor markers or the cell with highest diffusion component 1 (DC1) as root.
3. Check PAGA connectivity first to confirm clusters form a connected trajectory before computing pseudotime.
4. Use n_branchings=0 for linear trajectories; increase to 1+ only when branching is biologically expected.

Common Mistakes

• Running on non-continuous data: If cell types are discrete with no transitions, trajectory inference produces meaningless pseudotime.
• Wrong root cell: A wrong root inverts the trajectory — validate with known early-stage marker expression.
• Ignoring disconnected components: If PAGA shows disconnected clusters, some cells will get arbitrary pseudotime values.

Troubleshooting

Output

写入 --output 指定的 JSON 文件：

{
  "skill": "scrna-trajectory-inference",
  "status": "success",
  "steps": ["load_data", "paga", "diffusion_pseudotime", "palantir", "visualize", "pseudotime_genes", "heatmap", "save"],
  "result": {
    "input_file": "clustered.h5ad",
    "output_file": "trajectory.h5ad",
    "n_cells": 12000,
    "pseudotime_range": [0.0, 1.0],
    "root_cell_index": 42,
    "n_branches": 0,
    "top_pseudotime_genes": ["GeneA", "GeneB", "GeneC"],
    "pseudotime_genes_csv": "pseudotime_genes.csv",
    "plots": ["trajectory_analysis.png", "pseudotime_heatmap.png"]
  }
}

Workflow

Step 1: Load Data

import scanpy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

adata = sc.read_h5ad("clustered.h5ad")
assert "X_pca" in adata.obsm, "PCA embeddings required"
assert "neighbors" in adata.uns, "Neighbor graph required"

print(f"Input: {adata.shape}")
print(f"Clusters: {adata.obs['leiden'].nunique() if 'leiden' in adata.obs.columns else 'N/A'}")

Step 2: PAGA (Partition-based Graph Abstraction)

# Run PAGA to understand cluster connectivity
sc.tl.paga(adata, groups="leiden")
sc.pl.paga(
    adata,
    threshold=0.1,
    node_size_scale=2,
    edge_width_scale=0.5,
    title="PAGA Graph",
    save="_paga.png",
)

# Use PAGA-initialized UMAP for better trajectory visualization
sc.tl.umap(adata, init_pos="paga")

Step 3: Diffusion Pseudotime (DPT)

# Compute diffusion map
# ARPACK may fail to converge on some datasets; try with fewer components
use_dpt = False
for n_comps in [15, 10, 8]:
    try:
        sc.tl.diffmap(adata, n_comps=n_comps)
        root_idx = adata.obsm["X_diffmap"][:, 1].argmax()
        adata.uns["iroot"] = root_idx
        sc.tl.dpt(adata, n_branchings=0, n_dcs=n_comps)
        use_dpt = True
        break
    except Exception:
        continue

if not use_dpt:
    # Fallback: PCA-distance pseudotime
    from scipy.spatial.distance import cdist
    print("DPT failed, using PCA-distance pseudotime fallback")
    pca_coords = adata.obsm["X_pca"][:, :10]
    dists = cdist(pca_coords[0:1], pca_coords).flatten()
    adata.obs["dpt_pseudotime"] = dists / (dists.max() + 1e-10)

# Identify root cell (highest expression of stem/progenitor markers)
# Option 1: Manual selection based on known markers
# Option 2: Automatic - cell with highest DC1
root_cell = adata.obs.index[adata.obsm["X_diffmap"][:, 1].argmax()]
adata.uns["iroot"] = np.where(adata.obs.index == root_cell)[0][0]

# Compute pseudotime
sc.tl.dpt(adata, n_branchings=0, n_dcs=10)

print(f"Pseudotime range: {adata.obs['dpt_pseudotime'].min():.3f} - {adata.obs['dpt_pseudotime'].max():.3f}")

Step 4: Palantir Trajectory (Optional, more advanced)

try:
    import palantir

    # Run Palantir
    dm_res = palantir.utils.run_diffusion_maps(
        sc.utils.view_to_actual(adata),
        n_components=10,
    )
    ms_data = palantir.utils.determine_multiscale_space(adata)

    # Infer trajectory
    start_cell = root_cell  # use same root as DPT
    pr_res = palantir.core.run_palantir(
        adata,
        start_cell,
        num_waypoints=500,
    )

    # Palantir results
    adata.obs["palantir_pseudotime"] = pr_res.pseudotime
    adata.obs["palantir_entropy"] = pr_res.entropy

    print("Palantir trajectory computed")
except ImportError:
    print("Palantir not available, using DPT only")

Step 5: Visualize Trajectory

fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# Pseudotime on UMAP
sc.pl.umap(adata, color="dpt_pseudotime", ax=axes[0, 0], show=False,
           title="Diffusion Pseudotime", color_map="viridis")

# Clusters on UMAP
if "leiden" in adata.obs.columns:
    sc.pl.umap(adata, color="leiden", ax=axes[0, 1], show=False, title="Clusters")

# Pseudotime distribution
axes[1, 0].hist(adata.obs["dpt_pseudotime"], bins=50, alpha=0.7, edgecolor="black")
axes[1, 0].set_xlabel("Pseudotime")
axes[1, 0].set_ylabel("Cell Count")
axes[1, 0].set_title("Pseudotime Distribution")

# Pseudotime by cluster
if "leiden" in adata.obs.columns:
    cluster_pseudotime = adata.obs.groupby("leiden")["dpt_pseudotime"].mean().sort_values()
    axes[1, 1].barh(cluster_pseudotime.index.astype(str), cluster_pseudotime.values)
    axes[1, 1].set_xlabel("Mean Pseudotime")
    axes[1, 1].set_title("Pseudotime by Cluster")

plt.tight_layout()
plt.savefig("trajectory_analysis.png", dpi=150, bbox_inches="tight")
plt.close()

Step 6: Genes Correlated with Pseudotime

# Find genes that change along pseudotime
from scipy.stats import spearmanr

# Use raw data for gene expression
if adata.raw is not None:
    gene_data = adata.raw.to_adata()
else:
    gene_data = adata

# Compute correlation with pseudotime
pseudotime = adata.obs["dpt_pseudotime"]
correlations = {}

for gene in gene_data.var_names[:5000]:  # limit for speed
    if sparse.issparse(gene_data.X):
        expr = gene_data[:, gene].X.toarray().flatten()
    else:
        expr = gene_data[:, gene].X.flatten()
    corr, pval = spearmanr(pseudotime, expr)
    if not np.isnan(corr):
        correlations[gene] = {"correlation": corr, "pvalue": pval}

corr_df = pd.DataFrame(correlations).T
corr_df["abs_corr"] = corr_df["correlation"].abs()
corr_df = corr_df.sort_values("abs_corr", ascending=False)

# Save top genes
corr_df.head(100).to_csv("pseudotime_genes.csv")
print("\nTop 10 genes correlated with pseudotime:")
print(corr_df.head(10)[["correlation", "pvalue"]].to_string())

Step 7: Heatmap of Genes Along Pseudotime

# Plot top genes along pseudotime
top_genes = corr_df.head(20).index.tolist()

# Sort cells by pseudotime
cell_order = pseudotime.argsort()
adata_sorted = adata[cell_order]

fig, ax = plt.subplots(figsize=(12, 8))
if adata.raw is not None:
    plot_data = adata.raw[cell_order, :][:, top_genes]
else:
    plot_data = adata_sorted[:, top_genes]

if sparse.issparse(plot_data.X):
    plot_matrix = plot_data.X.toarray()
else:
    plot_matrix = plot_data.X

# Smooth along pseudotime
from scipy.ndimage import uniform_filter1d
smoothed = uniform_filter1d(plot_matrix, size=50, axis=0)

im = ax.imshow(smoothed.T, aspect="auto", cmap="RdBu_r",
               extent=[0, 1, 0, len(top_genes)])
ax.set_yticks(range(len(top_genes)))
ax.set_yticklabels(top_genes, fontsize=8)
ax.set_xlabel("Pseudotime (normalized)")
ax.set_title("Gene Expression Along Pseudotime")
plt.colorbar(im, ax=ax, label="Expression")
plt.savefig("pseudotime_heatmap.png", dpi=150, bbox_inches="tight")
plt.close()

Step 8: Save Results

adata.write_h5ad("trajectory.h5ad")
corr_df.to_csv("pseudotime_genes.csv")
print(f"Saved: {adata.shape}")

Key Parameters

Choosing the Right Method

Common Issues

• Wrong root cell: Check known stem/progenitor markers; manually set root
• Disconnected trajectory: Clusters may be unrelated; check PAGA connectivity
• No branching detected: Increase n_branchings or use Palantir
• Pseudotime reversed: Root cell may be wrong; flip with 1 - pseudotime