Spatial transcriptomics · trajectory analysis · pseudotime · Squidpy · Spateo
The Problem: Clustering Loses Spatial Continuity
Standard spatial analysis assigns each spot to a discrete cluster. But tissue biology is continuous — there is no hard boundary between stroma and tumor, between hypoxic core and normoxic periphery, between vascularized and avascular zones. Trajectory analysis reconstructs this continuum by ordering spots along a path through gene expression space, anchored to their physical location in tissue.
The result is pseudotime: a one-dimensional coordinate that captures how far each spot is along a biological transition. Paired with spatial coordinates, it reveals not just that gene expression changes, but exactly where in the tissue those changes happen.
Figure 1. Left: spatial trajectory colored by pseudotime from stroma (early) to tumor core (late). Right: gene expression dynamics along the trajectory — ECM genes decline, proliferation and immune markers peak at different positions.
Two Approaches: PAGA + Diffusion Pseudotime vs. Spateo
You have two practical options for spatial trajectory analysis:
PAGA + DPT (Scanpy/Squidpy) — Works on the spatial graph. Mature, well-documented. Best when you already have a spatial graph and want trajectories between defined regions. Requires a root cell/spot (the biological starting point).
Spateo — Purpose-built for spatial trajectories. Handles 3D spatial data and multiple trajectory branches. More flexible root selection. Heavier dependency chain.
For most Visium workflows, start with PAGA + DPT — it integrates cleanly with your existing Squidpy AnnData.
Running PAGA + DPT on Visium Data
import scanpy as sc
import squidpy as sq
import numpy as np
# Start from annotated Visium AnnData
adata = sc.read_h5ad("visium_annotated.h5ad")
# 1. Build neighborhood graph on HVGs (not spatial — for PAGA topology)
sc.pp.neighbors(adata, n_neighbors=15, use_rep="X_pca")
# 2. Run PAGA to get cluster connectivity
sc.tl.paga(adata, groups="spatial_domain")
sc.pl.paga(adata, threshold=0.1, node_size_scale=2)
# 3. Initialize UMAP from PAGA for better trajectory layout
sc.tl.umap(adata, init_pos="paga")
# 4. Set root: spot with highest stroma marker expression
stroma_score = adata[:, "COL1A1"].X.toarray().flatten()
adata.uns["iroot"] = int(np.argmax(stroma_score))
# 5. Run diffusion pseudotime
sc.tl.dpt(adata)
# Result: adata.obs["dpt_pseudotime"]Spatial Visualization of Pseudotime
Once pseudotime is computed, project it back onto tissue coordinates:
sq.pl.spatial_scatter(
adata,
color="dpt_pseudotime",
cmap="plasma",
size=1.5,
title="Pseudotime — stroma to tumor core",
)
# Overlay with cell type to validate the trajectory makes biological sense
sq.pl.spatial_scatter(
adata,
color=["dpt_pseudotime", "spatial_domain"],
ncols=2,
size=1.5,
)If pseudotime increases smoothly from stroma toward tumor core, the trajectory is capturing real biology. If it looks noisy or inverted, try a different root spot — pseudotime direction is set by your root choice, not by the algorithm.
Finding Genes That Change Along the Trajectory
The most useful output is the list of genes whose expression correlates with pseudotime — these are candidates for signals driving the stroma-to-tumor transition:
from scipy.stats import spearmanr
import pandas as pd
pt = adata.obs["dpt_pseudotime"].values
gene_matrix = adata.X.toarray() # or adata.layers["normalized"]
correlations = []
for i, gene in enumerate(adata.var_names):
rho, pval = spearmanr(pt, gene_matrix[:, i])
correlations.append({"gene": gene, "rho": rho, "pval": pval})
corr_df = pd.DataFrame(correlations)
corr_df["qval"] = corr_df["pval"] * len(corr_df) # Bonferroni
corr_df = corr_df[corr_df["qval"] < 0.05].sort_values("rho", ascending=False)
print("Top increasing along trajectory:")
print(corr_df.head(10)[["gene", "rho"]])
print("
Top decreasing:")
print(corr_df.tail(10)[["gene", "rho"]])Interpreting the Results
From the figure, four biologically coherent patterns emerge along the stroma-to-tumor trajectory:
- COL1A1 (ECM marker): Highest in stroma, decays toward tumor core — validates that the trajectory direction is correct
- CD8A (cytotoxic T cells): Peaks at the interface zone, consistent with T cell exclusion from the tumor core by immunosuppressive signals
- CXCL10 (interferon response): Transient peak at the interface — the zone of active immune engagement
- MKI67 (proliferation): Rises sharply in the tumor core, confirming that high-pseudotime spots correspond to actively dividing cancer cells
Combining with Ligand-Receptor Analysis
Trajectory analysis and LR analysis are complementary. Once you know which spatial zones correspond to early vs. late pseudotime, you can ask: which ligand-receptor interactions are active specifically at the interface? Filter spots by pseudotime bin before running gr.ligrec:
# Subset to interface zone (pseudotime 0.4–0.6)
interface_mask = (adata.obs["dpt_pseudotime"] > 0.4) & (adata.obs["dpt_pseudotime"] < 0.6)
adata_interface = adata[interface_mask].copy()
sq.gr.spatial_neighbors(adata_interface, coord_type="generic", n_neighs=6)
sq.gr.ligrec(adata_interface, n_perms=500, cluster_key="cell_type")This reveals whether the CD8A peak and CXCL10 burst at the interface are accompanied by specific T cell - cancer cell interactions that might be candidates for checkpoint blockade potentiation.
For the immunotherapy biology behind the interface dynamics — why T cells stall at the tumor boundary and what checkpoint inhibitors do about it — see NeoantigenLab.
From the NeoantigenLab sister blog
Working with tumor transcriptomics for neoantigen research?
NeoantigenLab covers neoantigen biology, WES pipelines, pVACseq, and HLA typing for experimental researchers — the biology behind what you're sequencing.
Visit NeoantigenLab →
