Spatial Trajectory Analysis: PAGA and Diffusion Pseudotime vs. Spateo

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 →

 

SPATIAL TRANSCRIPTOMICS · TUTORIAL

BANKSY: Using Neighborhood Context to Find Spatial Domains

Why adding just 20% neighbor signal outperforms pure expression clustering for tissue layer detection

July 2026  ·  9 min read  ·  SpatiaBio

Standard Leiden clustering on spatial data ignores one obvious fact: cells that sit next to each other tend to be in the same tissue layer. BANKSY (Building Aggregates with a Neighborhood Kernel for Yeoman Segmentation) fixes this by mixing each cell's own expression with a smoothed average of its neighbors — then clustering on the combined signal. The result is cleaner, more spatially coherent domains without requiring any graph-based post-processing.

Figure 1. BANKSY spatial domain map (left) and performance comparison vs. other methods by ARI and NMI (right).

The Core Idea: Lambda Controls the Mix

BANKSY introduces one key parameter: lambda (0 to 1), which controls how much neighbor context to blend in:

lambda = 0

Pure cell-intrinsic expression. Identical to standard Leiden.

lambda = 0.2

Sweet spot. Adds spatial coherence while preserving cell identity.

lambda = 0.8

Strong spatial smoothing. Best for laminar structures (brain cortex).

The embedding is constructed as: z = (1 - lambda) * x_self + lambda * x_neighbors, where x_neighbors is a weighted mean of expression in the spatial neighborhood (typically k=15-30 nearest neighbors in physical space).

Running BANKSY with Squidpy

import squidpy as sq
import scanpy as sc
import numpy as np
from banksy import Banksy

# Load and preprocess
adata = sc.read_visium('path/to/visium/')
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata = adata[:, adata.var.highly_variable]

# Build spatial graph (required for BANKSY)
sq.gr.spatial_neighbors(adata, coord_type='generic', n_neighs=15)

# Run BANKSY
banksy = Banksy(
    adata,
    banksy_lambda=0.2,   # neighbor mix ratio
    resolution=0.8,       # Leiden resolution
    pca_dims=20,
    n_neighbors=15,
)
banksy.fit()
banksy.plot_domains(adata, color='banksy_labels')

Installation: pip install banksy-py — works with AnnData/Squidpy natively.

BANKSY vs. Other Spatial Domain Methods

	BANKSY	BayesSpace	STAGATE	Leiden only
Approach	Neighbor smoothing + clustering	Bayesian HMRF	Graph autoencoder	Expression only
Speed	Fast	Slow	Slow (GPU)	Fastest
Scalability	Visium HD / Xenium OK	Visium only	Large datasets OK	Any
Interpretability	High (lambda is intuitive)	Low	Low (latent space)	High
Best for	General use, layered tissues	Small datasets, fine structure	Complex, non-laminar	Quick baseline

Practical Tips

• Start with lambda=0.2 for most tissues. Go higher (0.5-0.8) only for clearly laminar structures like brain cortex or retina.
• Run both BANKSY and Leiden and compare — BANKSY domains should look spatially cleaner. If they don't, your spatial graph might be poorly constructed.
• Use sq.gr.spatial_neighbors with coord_type='generic' for Visium. The default 'grid' mode can misrepresent actual tissue distances.
• BANKSY domains are not cell types — they are tissue microenvironments. Combine with cell type deconvolution (cell2location or Tangram) to get the full picture.

Key Takeaways

• BANKSY adds spatial coherence with a single parameter (lambda) on top of standard clustering
• lambda=0.2 is a reliable default; increase for laminar tissues
• Faster and more interpretable than deep learning alternatives for most use cases
• Best used alongside cell type deconvolution, not as a replacement

Interested in the immunology behind those spatial domains?

See NeoantigenLab for checkpoint inhibitors, CAR-T, and neoantigen biology.

Explore 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 →

SPATIAL TRANSCRIPTOMICS · TUTORIAL

Tangram: Mapping Single-Cell Annotations onto Spatial Coordinates

Transfer cell type labels, gene programs, and custom scores from scRNA-seq directly onto your Visium or Xenium data

July 2026  ·  10 min read  ·  SpatiaBio

Cell2location (covered in Post 12) infers cell type abundances per spot using a Bayesian model. Tangram takes a complementary approach: it learns a mapping from single-cell space to spatial space by aligning gene expression profiles, then transfers any annotation — cell types, transcription factor activity, pseudotime, custom gene scores — onto spatial coordinates. If you have a well-annotated scRNA-seq atlas and want to ask "where do these cells sit in tissue?", Tangram is your tool.

Figure 1. Tangram-mapped dominant cell type per spatial location (left) and deconvolution method comparison by accuracy and compute speed (right).

How Tangram Works

Tangram frames cell-to-space mapping as an optimal transport problem. Given a scRNA-seq dataset with C cells and a spatial dataset with S spots, it learns a mapping matrix M (C × S) that minimizes the difference between:

• The spatially reconstructed gene expression (M^T × scRNA matrix)
• The observed spatial gene expression

The key insight: you only need a shared set of marker genes to align the two modalities. Tangram then uses the learned mapping to project everything else — including genes not measured in the spatial assay.

import tangram as tg
import scanpy as sc
import squidpy as sq

# Load data
sc_adata  = sc.read_h5ad('sc_annotated.h5ad')   # scRNA-seq with cell_type labels
sp_adata  = sc.read_h5ad('visium_filtered.h5ad') # Visium spatial data

# Identify shared marker genes (top 100 per cell type works well)
sc.tl.rank_genes_groups(sc_adata, groupby='cell_type', method='wilcoxon')
markers = tg.pp.select_genes(sc_adata, n_genes=100)

tg.pp.filter_gene_list(sc_adata, sp_adata, gene_list=markers)

# Train mapping (GPU recommended for large datasets)
tg.mapping.map_cells_to_space(
    sc_adata,
    sp_adata,
    mode='cells',        # or 'clusters' for memory efficiency
    target_count=sp_adata.obs.cell_count.sum(),
    density_prior='rna_count_based',
    num_epochs=500,
    device='cpu',
)

Three Mapping Modes

cells mode

Maps each individual cell to a spot. Best for datasets <50k cells. Gives highest resolution but memory-intensive.

clusters mode

Maps cluster centroids instead of individual cells. Scalable to large atlases. Recommended starting point.

constrained mode

Uses prior cell density estimates (e.g. from DAPI) to constrain the mapping. More accurate when cell counts per spot are known.

Projecting Cell Types and Custom Scores

After mapping, Tangram stores results in sp_adata.obsm['tangram_ct_pred']. You can then project back to the spatial grid:

# Project cell type probabilities onto spatial coords
tg.pl.plot_cell_annotation_sc(sp_adata, annotation='cell_type', perc=0.02)

# Project any custom score from scRNA-seq (e.g. exhaustion score)
sc_adata.obs['exhaustion_score'] = (
    sc_adata[:, ['PDCD1','HAVCR2','LAG3','TIGIT','TOX']].X.mean(axis=1)
)
tg.pp.project_genes(adata_map=sp_adata, adata_sc=sc_adata)

# Now sp_adata has exhaustion_score — plot spatially
import squidpy as sq
sq.pl.spatial_scatter(sp_adata, color='exhaustion_score', cmap='Reds')

Pro tip: Project T cell exhaustion scores, EMT signatures, or pVACseq neoantigen burden per cell — then map them spatially. This links neoantigen-rich tumor zones to immune infiltration patterns.

Tangram vs cell2location: When to Use Which

	Tangram	cell2location
Approach	Optimal transport mapping	Bayesian deconvolution
Output	Cell-to-spot probability matrix	Cell type abundance per spot
Custom score transfer	✓ Any obs/obsm column	✗ Cell types only
Gene imputation	✓ Projects unmeasured genes	✗ No
Scalability	Moderate (<200k cells)	Large atlases OK
Best for	Transferring rich annotations	Quantifying cell composition

Common Pitfalls

• Too few marker genes — use at least 50-100 per cell type. With <20, mapping is unreliable.
• Batch effects between sc and spatial — if both datasets are from different labs or platforms, normalize carefully before running Tangram. Harmony or scVI integration first helps.
• Overinterpretation of imputed genes — genes projected via Tangram are predicted, not measured. Don't use imputed values for differential expression.
• Using cells mode on >100k cells — switch to clusters mode. The mapping matrix becomes O(C × S) and will OOM on most machines.

Key Takeaways

• Tangram maps any scRNA-seq annotation onto spatial coordinates via optimal transport
• Use clusters mode for large atlases, cells mode for maximum resolution
• Combine with Squidpy for spatial statistics on projected scores
• Complements cell2location — use both for a complete picture

Working on neoantigen biology or cancer immunotherapy?

The spatial exhaustion score workflow above connects directly to neoantigen research — see NeoantigenLab for the immunology side.

Explore 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 →

SpatialDE2: Finding Spatially Variable Genes Without Clustering

By Lociven · SpatiaBio · August 2026

Finding genes that vary spatially — not just between clusters, but continuously across tissue — is one of the core tasks in spatial transcriptomics. SpatialDE2 is the current standard for this. It's faster than its predecessor, handles Visium HD and single-cell resolution data, and gives you interpretable spatial patterns. This post covers how it works and how to run it.

---

What is a spatially variable gene?

A spatially variable gene (SVG) is a gene whose expression level is not randomly distributed across tissue — it shows spatial structure. That structure might be a sharp boundary (tumor vs stroma), a gradient (hypoxia from center to edge), or a periodic pattern (cortical layers in brain).

Standard differential expression tests (comparing cluster A vs cluster B) can find SVGs indirectly, but they depend on your clustering choices. SpatialDE2 tests for spatial variation directly, without needing predefined clusters, using a Gaussian process model.

SpatialDE2 vs Moran's I: Moran's I (available in Squidpy via sq.gr.spatial_autocorr) is faster and good for a quick screen. SpatialDE2 is slower but models the spatial covariance structure more accurately and provides pattern classification (linear gradient, periodic, etc.).

---

Installation

pip install SpatialDE
# SpatialDE2 requires Python 3.9+ and JAX
pip install jax jaxlib
pip install git+https://github.com/PMBio/SpatialDE.git

---

Running SpatialDE2

import squidpy as sq
import scanpy as sc
import SpatialDE
import numpy as np
import pandas as pd

# Load and prepare Visium data
adata = sq.read.visium("path/to/visium/")
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# Use highly variable genes to limit runtime
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata_hvg = adata[:, adata.var["highly_variable"]].copy()

# Get spatial coordinates
coords = pd.DataFrame(
    adata_hvg.obsm["spatial"],
    columns=["x", "y"],
    index=adata_hvg.obs_names
)

# Expression matrix (spots x genes)
expr = pd.DataFrame(
    adata_hvg.X.toarray() if hasattr(adata_hvg.X, "toarray") else adata_hvg.X,
    index=adata_hvg.obs_names,
    columns=adata_hvg.var_names
)

# Run SpatialDE2
results = SpatialDE.run(coords, expr)
results = results.sort_values("qvalue")

print(results.head(10)[["g", "l", "qvalue", "FSV"]])

Key output columns:

• g — gene name
• l — estimated spatial length scale (how broad the pattern is)
• qvalue — FDR-adjusted p-value for spatial variation
• FSV — fraction of variance explained by spatial structure (0–1; higher = more spatially structured)

---

The difference in one view. Left: Moran's I returns a ranking — these are the real top genes from the mouse-brain dataset (top 500 HVGs, n_perms=200). Right: a single score cannot distinguish a compact hotspot from a smooth gradient; all three shapes can score high. Classifying the pattern is what SpatialDE2's Gaussian-process test adds, and what the extra runtime pays for.

Pattern classification with SPARK-X

SpatialDE2 classifies SVGs into spatial patterns: linear gradients, radial patterns, and periodic patterns. For most tissue sections, linear and radial patterns dominate. Periodic patterns are most relevant in brain (cortical layers) and intestine (crypt-villus axis).

# Get top SVGs
sig_results = results[results["qvalue"] < 0.05].sort_values("FSV", ascending=False)
top_svgs = sig_results["g"].head(20).tolist()

print(f"Significant SVGs: {len(sig_results)}")
print(f"Top 5 by FSV:\n{sig_results[['g','FSV','l']].head()}")

# Visualize top SVGs on tissue
sq.pl.spatial_scatter(
    adata,
    color=top_svgs[:4],
    ncols=2,
    cmap="magma",
    size=1.5
)

---

Quick alternative: Moran's I in Squidpy

If you need a fast screen before committing to SpatialDE2, Squidpy's sq.gr.spatial_autocorr with Moran's I runs in minutes on a full Visium dataset:

sq.gr.spatial_neighbors(adata, coord_type="grid", n_rings=2)
sq.gr.spatial_autocorr(adata, mode="moran", genes=adata.var_names[adata.var.highly_variable])

# Top spatially autocorrelated genes
top_morans = (
    adata.uns["moranI"]
    .sort_values("I", ascending=False)
    .head(20)
)
print(top_morans[["I", "pval_norm_fdr_bh"]])

Use Moran's I to screen, then run SpatialDE2 on the top candidates for pattern classification and accurate FSV estimates.

---

Get the complete pack

SpatiaBio Pack 2 — Deconvolution & Spatial Pattern Detection

The full SpatialDE2 spatial pattern classification notebook — Moran's I screening, hotspot/gradient pattern detection, and spatial co-expression modules. Plus a cell2location deconvolution notebook.

Get it for $19 →

From the NeoantigenLab sister blog

Using spatial data for cancer immunology research?

NeoantigenLab covers neoantigen biology, pVACseq, HLA typing, and checkpoint inhibitors — the biology behind your spatial findings.

Visit NeoantigenLab →

Tags: SpatialDE2, spatially variable genes, SVG, Moran's I, spatial autocorrelation, Visium, Squidpy, Gaussian process, spatial transcriptomics