Scrna trajectory inference

Single-Cell Trajectory Inference

When to Use This Skill

Use when you have preprocessed scRNA-seq data and want to:

• ✅ Order cells along a differentiation or disease trajectory (pseudotime)
• ✅ Identify branching points and terminal cell fates
• ✅ Discover genes driving cell state transitions
• ✅ Visualize RNA velocity (direction of cell state change)
• ✅ Compute cell fate probabilities with CellRank
• ✅ Chain from scrnaseq-scanpy-core-analysis output

Do NOT use when:

• ❌ Data is not yet preprocessed (use scrnaseq-scanpy-core-analysis first)
• ❌ You have bulk RNA-seq (use disease-progression-longitudinal instead)
• ❌ Cells are terminally differentiated with no trajectory (e.g., resting PBMCs)
• ❌ Fewer than 200 cells

Installation

pip install scanpy anndata scvelo cellrank numpy pandas matplotlib seaborn scipy statsmodels reportlab

Package	Version	License	Commercial Use	Notes
scanpy	≥1.9	BSD-3	✅ Permitted	Core trajectory (PAGA, DPT)
anndata	≥0.8	BSD-3	✅ Permitted	Data container
scvelo	≥0.2.5	BSD-3	✅ Permitted	RNA velocity (optional but recommended)
cellrank	≥2.0	BSD-3	✅ Permitted	Fate mapping (optional)
matplotlib	≥3.4	PSF	✅ Permitted	Plotting
seaborn	≥0.11	BSD-3	✅ Permitted	Statistical plotting, heatmaps
scipy	≥1.7	BSD-3	✅ Permitted	Statistics
statsmodels	≥0.13	BSD-3	✅ Permitted	FDR correction
reportlab	≥3.6	BSD	✅ Permitted	PDF report (optional)

Graceful degradation: Core analysis (PAGA + pseudotime) requires only scanpy. scVelo and CellRank are optional — scripts detect availability and skip gracefully.

Inputs

Required:

• Preprocessed AnnData (.h5ad) with PCA, UMAP, and cluster annotations
  • Output from scrnaseq-scanpy-core-analysis (adata_processed.h5ad) works directly
  • Must have ≥200 cells, ≥100 genes, cluster labels in .obs

For RNA velocity (optional):

• Spliced/unspliced count layers (adata.layers['spliced'], adata.layers['unspliced'])
• Generated by STARsolo, Cell Ranger, or velocyto

Outputs

Analysis objects (for downstream skills):

• adata_trajectory.h5ad — AnnData with pseudotime, PAGA, diffusion map embedded
  • Load with: adata = sc.read_h5ad('adata_trajectory.h5ad')
  • Required for: downstream enrichment, regulatory network analysis
• trajectory_results.pkl — Full results dict (pseudotime, gene lists, model objects)
  • Load with: results = pickle.loads(Path('trajectory_results.pkl').read_bytes())

Primary results (CSV):

• pseudotime_assignments.csv — Cell barcode, pseudotime, cell type
• trajectory_genes.csv — Genes correlated with pseudotime (gene, correlation, FDR, direction)
• velocity_genes.csv — Top RNA velocity genes (if scVelo ran)
• fate_probabilities.csv — Cell fate probabilities (if CellRank ran)
• driver_genes_*.csv — Driver genes per terminal fate (if CellRank ran)

Visualizations (PNG + SVG at 300 DPI):

• paga_graph.png/.svg — PAGA cluster connectivity + UMAP
• pseudotime_umap.png/.svg — UMAP colored by pseudotime
• pseudotime_violin.png/.svg — Pseudotime distribution per cell type
• diffusion_components.png/.svg — Diffusion map components
• gene_heatmap.png/.svg — Top trajectory genes heatmap
• gene_trends.png/.svg — Gene expression trends along pseudotime
• paga_connectivity.png/.svg — PAGA connectivity heatmap
• velocity_stream.png/.svg — RNA velocity stream plot (if scVelo)
• velocity_confidence.png/.svg — Velocity confidence (if scVelo)
• latent_time.png/.svg — scVelo latent time (if dynamical model)
• velocity_top_genes.png/.svg — Phase portraits for top velocity genes (if scVelo)
• fate_probabilities.png/.svg — Cell fate UMAP (if CellRank)
• fate_heatmap.png/.svg — Fate probability heatmap (if CellRank)
• driver_genes.png/.svg — Top driver genes per terminal fate (if CellRank)

Reports:

• trajectory_analysis_report.pdf — Publication-quality PDF (requires reportlab)
• trajectory_analysis_report.md — Markdown fallback report
• analysis_metadata.json — Parameters and quality metrics

Clarification Questions

🚨 ALWAYS ask Question 1 FIRST. Do not ask about analysis parameters before the user has answered Question 1.

1. Input Files (ASK THIS FIRST)

• Do you have a preprocessed scRNA-seq object (.h5ad)?
  • If uploaded: Is this your processed AnnData file?
  • Expected: .h5ad with PCA, UMAP, and cluster annotations
• Or use example/demo data?
  • Pancreatic endocrinogenesis (3,696 cells, branching trajectory into alpha/beta/delta/epsilon cells)

🚨 IF EXAMPLE DATA SELECTED: All parameters are pre-defined. Only ask Question 2. Then proceed to Step 1.

2. Analysis Scope (structured — works for both demo and user data)

• a) Core trajectory only — PAGA + pseudotime (~3 min)
• b) Core + RNA velocity — adds scVelo streams (~5 min) (recommended)
• c) Full analysis — adds CellRank fate mapping (~8 min)

Questions 3-5 are ONLY for users providing their own data:

3. Root Cell Type

• Which cell type represents the starting point of the trajectory?
• (For demo data: Ductal cells are the progenitors — pre-selected)

4. Cluster Key

• Which column in .obs contains your cell type annotations?
• Common: clusters, cell_type, leiden, louvain

5. Expected Trajectory Structure

• a) Linear differentiation (one lineage)
• b) Branching (multiple fates from one progenitor)
• c) Not sure — let the data decide

Standard Workflow

Note: Run from the OmicsClaw root directory and add the workflow scripts to sys.path:

import sys; import os; sys.path.insert(0, os.path.abspath('knowledge_base/scripts/scrna-trajectory-inference'))

🚨 MANDATORY: USE SCRIPTS EXACTLY AS SHOWN — DO NOT WRITE INLINE CODE 🚨

Step 1 — Load data:

from load_example_data import load_example_data
adata = load_example_data()

DO NOT write inline data loading code. Just use the script.

✅ VERIFICATION: You MUST see: "✓ Data loaded successfully!"

Step 2 — Run trajectory analysis:

from run_trajectory_analysis import run_trajectory
results = run_trajectory(adata, root_cell_type="Ductal", cluster_key="clusters")

DO NOT write inline trajectory code. Just use the script.

✅ VERIFICATION: You MUST see: "✓ Trajectory analysis completed successfully!"

Step 3 — Generate visualizations:

from generate_all_plots import generate_all_plots
generate_all_plots(adata, results, output_dir="trajectory_results", cluster_key="clusters")

🚨 DO NOT write inline plotting code. Just use the script. 🚨

✅ VERIFICATION: You MUST see: "✓ All plots generated successfully!"

Step 4 — Export results:

from export_results import export_all
export_all(adata, results, output_dir="trajectory_results")

DO NOT write custom export code. Use export_all().

✅ VERIFICATION: You MUST see: "=== Export Complete ==="

---

⚠️ CRITICAL — DO NOT:

• ❌ Write inline data loading code → STOP: Use load_example_data() or load_user_data()
• ❌ Write inline PAGA/DPT/scVelo code → STOP: Use run_trajectory()
• ❌ Write inline plotting code (plt.savefig, sc.pl, etc.) → STOP: Use generate_all_plots()
• ❌ Write custom export code → STOP: Use export_all()

⚠️ IF SCRIPTS FAIL — Script Failure Hierarchy:

1. Fix and Retry (90%) — Install missing package, re-run script
2. Modify Script (5%) — Edit the script file itself, document changes
3. Use as Reference (4%) — Read script, adapt approach, cite source
4. Write from Scratch (1%) — Only if genuinely impossible, explain why

NEVER skip directly to writing inline code without trying the script first.

Common Issues

Error	Cause	Solution
"Root cell type not found"	Cluster name mismatch	Check adata.obs['clusters'].unique() for exact names
"scvelo not installed"	Missing optional dependency	pip install scvelo — or skip velocity (core trajectory still works)
"cellrank not installed"	Missing optional dependency	pip install cellrank — or skip fate mapping
scVelo dynamical model fails	Insufficient spliced/unspliced counts	Script auto-falls back to stochastic model
SVG export failed	Missing system library	Normal — PNG still generated. Script handles fallback.
"Too few cells"	Dataset too small	Need ≥200 cells for meaningful trajectory
CellRank terminal states incorrect	Auto-detection picked wrong states	Check PAGA graph to verify expected terminal fates

Suggested Next Steps

1. Functional enrichment → Use functional-enrichment-from-degs skill

   • Input: trajectory_genes.csv (top up/down genes along pseudotime)
   • Find pathways driving differentiation

2. Gene regulatory networks → Use grn-pyscenic skill

   • Input: adata_trajectory.h5ad
   • Identify transcription factors controlling cell fate decisions

3. Upstream regulator analysis → Use upstream-regulator-analysis skill

   • Input: trajectory_genes.csv
   • Predict upstream regulators of trajectory dynamics

4. Differential expression at branch points → Use scrnaseq-scanpy-core-analysis (marker finding)

   • Compare cells at branch points to identify fate-determining genes

Related Skills

Upstream (data generation):

• scrnaseq-scanpy-core-analysis — Preprocessing, clustering, UMAP → feeds .h5ad into this skill
• scrnaseq-seurat-core-analysis — R alternative (convert with SeuratDisk)

Downstream (interpretation):

• functional-enrichment-from-degs — Pathway analysis of trajectory genes
• grn-pyscenic — Gene regulatory networks
• upstream-regulator-analysis — Upstream regulators of trajectory genes

Alternative trajectory methods:

• disease-progression-longitudinal — Bulk/multi-omics longitudinal trajectories (TimeAx)

References

Primary Citations

1. PAGA: Wolf FA, Hamey FK, Plass M, et al. PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol. 2019;20:59.

2. Diffusion Pseudotime: Haghverdi L, Büttner M, Wolf FA, et al. Diffusion pseudotime robustly reconstructs lineage branching. Nat Methods. 2016;13:845-848.

3. scVelo: Bergen V, Lange M, Peidli S, et al. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat Biotechnol. 2020;38:1408-1414.

4. CellRank: Lange M, Bergen V, Klein M, et al. CellRank for directed single-cell fate mapping. Nat Methods. 2022;19:159-170.

5. Example dataset: Bastidas-Ponce A, Tritschler S, Dony L, et al. Comprehensive single cell mRNA profiling reveals a detailed roadmap for pancreatic endocrinogenesis. Development. 2019;146:dev173849.

Software

Software	Version	License	Commercial Use
scanpy	≥1.9	BSD-3	✅ Permitted
scVelo	≥0.2.5	BSD-3	✅ Permitted
CellRank	≥2.0	BSD-3	✅ Permitted
seaborn	≥0.11	BSD-3	✅ Permitted
reportlab	≥3.6	BSD	✅ Permitted