Scrna Trajectory Inference

---
id: scrna-trajectory-inference
name: Single-Cell Trajectory Inference
category: transcriptomics
short-description: "Infer differentiation trajectories, pseudotime ordering, RNA velocity, and cell fate probabilities from scRNA-seq data."
detailed-description: "Reconstruct developmental or differentiation trajectories from single-cell RNA-seq data using PAGA, diffusion pseudotime, scVelo RNA velocity, and CellRank fate mapping. Discovers cell ordering along pseudotime, identifies branching points and terminal fates, and reveals gene expression dynamics along trajectories. Chains from scrnaseq-scanpy-core-analysis or any preprocessed AnnData (.h5ad). Produces publication-ready trajectory visualizations and structured PDF reports."
starting-prompt: Infer differentiation trajectories from my single-cell RNA-seq data using PAGA, pseudotime, and RNA velocity
---

# 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

```bash
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`:
> ```python
> 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:**
```python
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:**
```python
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:**
```python
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:**
```python
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 |