Squidpy's Delaunay Radius Pruning Builds a Graph You Can't Save — Here's the Workaround

Squidpy · spatial_neighbors_delaunay · h5ad · IORegistryError · still in latest release

Squidpy's Delaunay Radius Pruning Builds a Graph You Can't Save — Here's the Workaround

Status as of writing: reproducible on Squidpy's latest stable release, v1.8.3 (2026-06-29). A fix merged to main via PR #1249 on 2026-08-11, but as of this writing that fix has not yet shipped in a release — if you installed from PyPI, you still hit this. Reported in squidpy#1247.

If you build a spatial graph with radius-pruned Delaunay triangulation and then try to save the object, the graph itself builds fine — the error only shows up later, at write_h5ad():

sq.gr.spatial_neighbors_delaunay(adata, radius=(0.0, 150.0))
# runs fine — graph is built, edges are pruned correctly

adata.write_h5ad("out.h5ad")
# IORegistryError: No method registered for writing 
#                  into 
# Error raised while writing key 'radius' of ... to /uns/spatial_neighbors/params

Why it only fails on save

The graph construction works because nothing about building the Delaunay graph needs to serialize anything. The tuple only becomes a problem when write_h5ad tries to write adata.uns to disk — adata.uns['spatial_neighbors']['params']['radius'] is stored as a native Python tuple, and anndata's h5ad writer has no registered encoder for a bare tuple inside an HDF5 group — only for types it explicitly knows how to serialize (lists, arrays, scalars, dicts of the same).

This is exactly the kind of failure that's easy to lose an afternoon to if you're new to the ecosystem: the traceback points at anndata's I/O registry deep in a stack of writer/wrapper functions, with no obvious link back to the one line that actually caused it (passing a tuple radius argument, dozens of lines earlier and a different function call entirely).

The workaround (works today, no waiting for the release)

Cast the stored value to a list before saving — anndata knows exactly how to serialize a list:

sq.gr.spatial_neighbors_delaunay(adata, radius=(0.0, 150.0))

# fix in place, before writing:
adata.uns['spatial_neighbors']['params']['radius'] = list(
    adata.uns['spatial_neighbors']['params']['radius']
)

adata.write_h5ad("out.h5ad")   # now succeeds

Alternatively, if you don't need the parameter recorded, del adata.uns['spatial_neighbors']['params']['radius'] before saving — but the list-cast is preferable since it keeps the record of what radius you actually used, which matters for reproducibility.

When the real fix lands

PR #1249 ("ensure radius values stored in AnnData.uns are writable") changes Squidpy's own storage code to cast to a writable type before it ever reaches adata.uns, so once it ships in a release you can drop the manual workaround above. Until then — or if you're pinned to an older version for other reasons — the one-line cast is a complete fix, not a partial one; it doesn't change your graph, only how the parameter gets serialized.

Verified against Squidpy's main branch, PR #1249, and the current v1.8.3 release, 2026-09. Source: squidpy#1247.

"UnboundLocalError: counts_per_cols" — What Disabling Numba Actually Broke in scanpy

scanpy · normalize_total · Numba · UnboundLocalError · reproducible builds

"UnboundLocalError: counts_per_cols" — What Disabling Numba Actually Broke in scanpy

Status as of writing: present in scanpy 1.11.5 (and versions around it), fixed via PR #4231 (merged 2026-07-27), shipped in 1.12.4 (2026-08-27). Reported in scanpy#4228.

If you disable Numba JIT — a common workaround when a container or CI environment can't write Numba's cache files — sc.pp.normalize_total() could crash outright on CSR matrices with the (default) setting exclude_highly_expressed=False:

NUMBA_DISABLE_JIT=1 python - <<'PY'
import numpy as np
from scipy.sparse import csr_matrix
from anndata import AnnData
import scanpy as sc

adata = AnnData(csr_matrix(np.array([[1, 2, 0], [0, 3, 4]], dtype=np.int32)))
sc.pp.normalize_total(adata, target_sum=1e4)
PY

# UnboundLocalError: cannot access local variable 'counts_per_cols'
# where it is not associated with a value

Why disabling Numba even matters

This isn't an exotic setup. People turn off Numba JIT (via NUMBA_DISABLE_JIT=1) for a few ordinary reasons:

• The runtime environment can't write to Numba's cache directory (common on locked-down CI runners or read-only containers)
• Debugging — JIT-compiled functions are harder to step through
• Cold-start latency in short-lived jobs where JIT compilation overhead isn't worth it

The actual bug

The CSR-specific fast path, _normalize_csr, only initialized counts_per_cols inside the if exclude_highly_expressed: block, but returned it unconditionally at the end of the function. With Numba's JIT active, this quirk of the compiled code path apparently never surfaced the same way; with JIT disabled, plain Python's stricter local-variable scoping raised the UnboundLocalError the moment the default-off branch was taken.

The fix (confirmed by reading the current source) simply initializes counts_per_cols = np.zeros(columns, dtype=np.int32) unconditionally before the branch, so it always has a value to return:

# current scanpy main, _normalization.py
def _normalize_csr(mat, rows, columns, exclude_highly_expressed=False, ...):
    counts_per_cols = np.zeros(columns, dtype=np.int32)   # <- now unconditional
    ...
    if exclude_highly_expressed:
        counts_per_cols_t = np.zeros((n_threads, columns), dtype=np.int32)
        ...
    return counts_per_cell, counts_per_cols

What to do

Upgrade to scanpy ≥1.12.4. If you're stuck on an older pinned version and need Numba disabled, the workaround is to explicitly pass exclude_highly_expressed=True even if you don't need the exclusion logic — it forces the code path that does initialize the variable, sidestepping the crash (at the cost of a slightly different internal computation you should verify doesn't change your output if you truly don't want exclusion).

Verified against scanpy's main branch source and PR #4231, 2026-09. Source: scanpy#4228.

Why scanpy Warns "adata.X seems to be already log-transformed" Even When You Never Touched X

scanpy · log1p · multi-layer normalization · false-positive warning

Why scanpy Warns "adata.X seems to be already log-transformed" Even When You Never Touched X

Status as of writing: open, reproduced against scanpy main, 2026-09. Reported in scanpy#4285. Not a data-loss bug — a misleading warning that can hide a real one.

If you work with spliced/unspliced layers (RNA velocity workflows) or any multi-layer object, you've probably run something like this and been confused by the output:

sc.pp.log1p(adata)
sc.pp.log1p(adata, layer="spliced")
sc.pp.log1p(adata, layer="unspliced")

# WARNING: adata.X seems to be already log-transformed.
# WARNING: adata.X seems to be already log-transformed.

The last two calls target layer="spliced" and layer="unspliced" — they don't touch .X at all. Yet you get the "X seems to be already log-transformed" warning twice.

The actual check

Confirmed by reading scanpy's _simple.py directly:

@log1p.register(AnnData)
def log1p_anndata(adata, *, base=None, copy=False, chunked=False,
                   chunk_size=None, layer=None, obsm=None):
    if "log1p" in adata.uns:
        logg.warning("adata.X seems to be already log-transformed.")
    ...

The check is if "log1p" in adata.uns — a single flag that gets set the first time you log1p anything on the object, and is then checked unconditionally on every subsequent call regardless of the layer argument. It never looks at the actual target layer, only at whether .X was ever log1p'd before, on any call.

Why this is more than cosmetic

Two real costs, neither of which is "the code breaks":

• Alert fatigue. If you log1p three layers per sample across dozens of samples in a loop, you get a wall of identical-looking warnings that are individually meaningless — which trains you to stop reading them. The one time this warning is actually telling you something real (you really did call log1p on .X twice — the double-normalization class of bug we've written about before), it's buried in noise.
• It's a global flag, not a per-layer one. There's no way from the warning text alone to tell which of your calls it's actually about.

What to do about it now

Since the fix isn't merged yet, the practical workaround is to not rely on the built-in warning to catch double-log1p at all — check explicitly instead:

# Before log1p-ing a layer, check its own max value directly —
# don't rely on adata.uns['log1p'] or the warning text.
import numpy as np

def already_logged(matrix, threshold=20):
    m = matrix.max() if hasattr(matrix, "max") else matrix.toarray().max()
    return m < threshold  # raw UMI counts are almost always well above this

for layer in ["spliced", "unspliced"]:
    if not already_logged(adata.layers[layer]):
        sc.pp.log1p(adata, layer=layer)

This is the same lesson as our earlier double-normalization post, applied one level up: don't trust a global "has this been transformed" signal on an object that may hold several independently-processed matrices. Check the specific array you're about to touch.

Verified directly against scanpy's main branch source, 2026-09. Source: scanpy#4285.

More preprocessing landmines

The double-normalization bug that hit our own notebooks

Pack 1 documents 7 real API-drift bugs found while re-executing all 16 notebooks — exact errors, root causes, fixes.

Get Pack 1 for $19 →

A Year-Old scanpy Bug That Could Silently Zero Out Your Data — Fixed Last Month

scanpy · normalize_total · silent data loss · CSR vs dense · preprocessing bug

A Year-Old scanpy Bug That Could Silently Zero Out Your Data — Fixed Last Month

We've written before about how easy it is to accidentally double-log-normalize spatial data (it happened across our own corpus — see our earlier workflow post). This one is a different bug in the same function, and it's worse: it can silently zero out your entire matrix with no error, just two easy-to-miss RuntimeWarnings.

Status as of writing: present in scanpy 1.11.2 through 1.12.3 (2025-05-28 → 2026-07-24). Fixed in 1.12.4 (released 2026-08-27) via PR #4256, reported in scanpy#4251. If you pinned scanpy any time in the last 14 months, check your version.

What target_sum=None is supposed to mean

sc.pp.normalize_total(adata, target_sum=None) — the default most tutorials use — documents its target as "the median of total counts for observations." Simple enough. But _normalize_total_helper computes that median two different ways depending on whether your matrix is CSR sparse or dense:

# CSR path (scanpy/preprocessing/_normalization.py, pre-1.12.4)
if target_sum is None:
    target_sum = np.median(counts_per_cell)            # median over ALL cells

# dense / everything-else path
if target_sum is None:
    target_sum = _compute_nnz_median(counts_per_cell)   # median over POSITIVE cells only

The two only disagree when a zero-total cell shifts the median — with row sums like [0, 10, 10, 20] both give 10, which is probably why this went unnoticed for so long. Empty cells are common in practice: an unfiltered raw matrix, an aggressive gene subset, a .raw slice with different filtering than the parent object.

The severe case: when at least half your cells are empty

If enough cells have zero total counts, the CSR branch's median lands on 0. Dividing by zero doesn't error — it silently produces an all-zero matrix, with only the warnings below as evidence anything went wrong:

import anndata as ad
import numpy as np
import scanpy as sc
import scipy.sparse as sp

# row sums [0, 0, 10] — the only populated cell has counts [4, 6]
dense2 = np.array([[0.0, 0.0], [0.0, 0.0], [4.0, 6.0]], dtype=np.float32)

for label, X in (("dense", dense2.copy()), ("csr", sp.csr_matrix(dense2))):
    a = ad.AnnData(X=X)
    sc.pp.normalize_total(a, target_sum=None)
    out = a.X if isinstance(a.X, np.ndarray) else a.X.toarray()
    print(f"{label:>5}: {out.ravel()}")

# dense: [0. 0.  0. 0.  4. 6.]   <- correct, target=10, populated cell untouched
#   csr: [0. 0.  0. 0.  0. 0.]   <- target=0, the ONLY real data point is destroyed

.../scanpy/preprocessing/_normalization.py:119: RuntimeWarning: divide by zero encountered in divide
.../scanpy/preprocessing/_normalization.py:119: RuntimeWarning: invalid value encountered in divide

No exception. No crash. Just a warning easy to filter out in a notebook with dozens of other warnings, and a matrix that's now all zeros.

Why the CSR path had a different rule at all

Before 1.11.2 there was a single code path, and it excluded zero-count cells on both branches. The split arrived with the numba CSR kernel in PR #3571 (merged 2025-05-19), which kept the old positive-cell-median rule for the non-CSR branch but gave the new CSR branch a plain np.median. Neither that PR nor its release notes mention a semantic change — this reads as an oversight, not a decision, and the reporter's one-line fix (route the CSR branch through the same _compute_nnz_median helper) is exactly what landed in 1.12.4.

What to actually do

• Upgrade to scanpy ≥1.12.4 if you can — this is the real fix.
• If you're pinned to an older version, pass an explicit target_sum (e.g. target_sum=1e4) instead of relying on the median default — this sidesteps the bug entirely regardless of version.
• Before normalizing, check (adata.X.sum(axis=1) == 0).sum() — if it's more than a small fraction of your cells, filter them out first rather than letting either branch guess a target around them.

Verified directly against scanpy's main branch source and the released 1.12.4 tag, 2026-09. Source: scanpy#4251.

Building a Spatial Transcriptomics Workflow: From cell2location to BANKSY and Beyond

Spatial transcriptomics · complete workflow · best practices · Squidpy · publication figures

Putting It All Together

This series has covered spatial transcriptomics from first principles (Squidpy post 1) through advanced methods — cell2location, SpatialDE2, Tangram, BANKSY, ligand-receptor analysis, trajectory analysis, and platform comparisons. This final post assembles them into a decision-driven workflow you can apply to any new Visium dataset.

The goal is not to run every tool — it is to know which tools answer which questions, in what order, and at what computational cost.

Figure 1. Left: complete spatial analysis pipeline from raw data to publication figures. Right: approximate runtime and peak memory per tool on a standard 8-core laptop with a 3,500-spot Visium dataset.

Phase 1: Data Loading and QC (30 minutes)

Every analysis starts here. Do not skip the QC plots — spatial data has tissue-specific artifacts (necrotic regions, edge effects, adipose tissue with near-zero counts) that will distort downstream results if not filtered.

import squidpy as sq
import scanpy as sc
import matplotlib.pyplot as plt

adata = sq.read.visium("spaceranger_output/", library_id="sample_01")

# QC metrics
sc.pp.calculate_qc_metrics(adata, inplace=True)
sq.pl.spatial_scatter(adata, color=["total_counts", "n_genes_by_counts"], ncols=2)

# Filter: remove low-count spots (necrosis, edge artifacts)
sc.pp.filter_cells(adata, min_counts=500)
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=10)

# Normalize
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=3000, flavor="seurat_v3")

# Save raw counts for cell2location
adata.layers["counts"] = adata.X.copy()
adata.raw = adata

Phase 2: Clustering and Annotation (1–2 hours)

Build spatial graph first, then standard dimensionality reduction. The spatial graph is required for neighborhood enrichment and LR analysis later.

sq.gr.spatial_neighbors(adata, coord_type="grid", n_rings=2)

sc.tl.pca(adata, use_highly_variable=True)
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.3)

# Inspect: do clusters correspond to tissue morphology?
sq.pl.spatial_scatter(adata, color="leiden", size=1.5)
sc.pl.umap(adata, color="leiden")

# Annotate clusters with known markers
marker_genes = {
    "Cancer": ["EPCAM", "KRT8", "MKI67"],
    "T cells": ["CD3E", "CD8A", "CD4"],
    "Macrophages": ["CD68", "CD163", "SPP1"],
    "Fibroblasts": ["COL1A1", "FAP", "ACTA2"],
    "Endothelial": ["PECAM1", "VWF", "CDH5"],
}
sc.pl.dotplot(adata, marker_genes, groupby="leiden")

Phase 3: Choose Your Analysis Path

After annotation, your research question determines what to run next. Use this decision tree:

Q: What cell types are in each spot, at single-cell resolution?
→ cell2location (Post 12) — requires matched scRNA-seq reference

Q: Which genes show non-random spatial patterns?
→ SpatialDE2 (Post 14) or sq.gr.spatial_autocorr for a faster first pass

Q: Are cell types communicating via known ligand-receptor pairs?
→ sq.gr.nhood_enrichment + sq.gr.ligrec (Post 17)

Q: How does gene expression change across tissue regions continuously?
→ PAGA + DPT trajectory (Post 18)

Q: Are there spatial domains not captured by gene expression alone?
→ BANKSY (Post 16) — incorporates neighbor expression into clustering

Q: I have scRNA-seq from a different tissue; can I map those cells spatially?
→ Tangram (Post 15) — optimal transport mapping

Phase 4: Publication Figure Checklist

Before submitting, verify each figure against this checklist:

# 1. Spatial scatter — confirm spot size is set for the journal's figure width
sq.pl.spatial_scatter(adata, color="cell_type", size=1.2,
                       figsize=(6, 6), dpi=300, save="fig1a.pdf")

# 2. Color palettes — use perceptually uniform maps for continuous data
#    Use categorical palettes with enough contrast for 8+ cell types
import matplotlib.pyplot as plt
# Bad: plt.cm.jet (never use in publications)
# Good: "plasma", "viridis", "RdBu_r" for continuous; "tab20" for categories

# 3. Scale bar — always include physical scale
# sq.pl.spatial_scatter adds it automatically from adata.uns["spatial"] metadata

# 4. Export at 300 dpi minimum for print
plt.savefig("figure1.pdf", dpi=300, bbox_inches="tight")
plt.savefig("figure1.png", dpi=300, bbox_inches="tight")

Compute Budget for a Full Analysis

From the benchmark figure, planning a full Visium TME analysis on a laptop (8 cores, 16 GB RAM):

• QC + clustering: ~30 min, <2 GB RAM
• cell2location (1000 epochs): ~45 min, 6 GB RAM — the bottleneck; use GPU if available
• SpatialDE2: ~8 min, 3.5 GB RAM
• sq.gr.ligrec (1000 perms): ~12 min, 2 GB RAM
• PAGA + DPT: ~5 min, 1.5 GB RAM
• Total: ~1.5–2 hours end-to-end

If you are running on a server with a GPU, cell2location drops to under 10 minutes and the full pipeline completes in well under an hour.

What This Series Has Not Covered

Spatial transcriptomics is a fast-moving field. Topics worth following as they mature:

• 3D spatial reconstruction — stacking serial sections with PASTE2 or Spateo
• Spatial proteomics integration — combining Visium with CODEX or IMC protein data
• Spatial multi-sample integration — harmonizing multiple patients/conditions with Harmony or scVI adapted for spatial
• Foundation models — Geneformer and scGPT adapted for spatial context are emerging

All code from this series is available on GitHub. The Squidpy Starter Pack ($19) contains runnable Jupyter notebooks for posts 1–5 with a real Visium PBMC dataset. For the immunotherapy biology behind the cell types you are mapping, see the sister blog NeoantigenLab.

Visium HD: What Changes When Resolution Goes Subcellular

Spatial transcriptomics · Visium HD · Xenium · CosMx · platform comparison

Why Resolution Matters (and Why It Is Not Everything)

Visium HD delivers 2 µm bins — roughly single-cell resolution — while the original Visium captures 55 µm spots that each contain 5–20 cells. That 27x improvement in resolution sounds decisive. In practice, the right platform depends on your biological question, not just the spec sheet.

This post covers what Visium HD actually changes analytically, where it outperforms its predecessor, and when Xenium or CosMx is the better choice despite having a fixed gene panel.

Figure 1. Left: spatial resolution comparison — Visium HD achieves 2 µm bins vs. 55 µm spots for Visium v1. Right: capability trade-offs across platforms normalized to maximum observed values.

What Visium HD Changes: The Bin Model

Visium HD does not capture individual cells — it captures 2 µm square bins that tile the tissue continuously. The key workflow difference is binning strategy: you choose a bin size at analysis time, not at capture time.

import squidpy as sq
import scanpy as sc

# Load Visium HD at 8 µm bins (good starting point — balances resolution and sparsity)
adata = sq.read.visium(
    path="visium_hd_output/",
    counts_file="filtered_feature_bc_matrix.h5",
    library_id="your_sample_id",
    bin_size=8,  # options: 2, 8, 16 micrometers
)

print(adata)
# AnnData: n_obs (bins) depends heavily on bin_size
# bin_size=2:  ~500,000 bins (very sparse, ~2-5 UMI/bin)
# bin_size=8:  ~30,000 bins (workable, ~50-100 UMI/bin)
# bin_size=16: ~8,000 bins (Visium-like density)

Start with 8 µm bins for exploratory analysis. At 2 µm, most bins will have fewer than 5 UMI — too sparse for reliable clustering without heavy smoothing. Only use 2 µm bins if you need subcellular resolution for a specific question (e.g., resolving nuclear vs. cytoplasmic transcript localization).

The Critical Difference: Deconvolution Is Optional

With Visium v1, cell2location or RCTD is mandatory to get cell type proportions per spot. With Visium HD at 8 µm bins, many bins will contain a single dominant cell type and can be annotated directly by marker expression — though mixed bins still exist at tissue interfaces.

# Standard clustering works directly on HD bins
sc.pp.filter_cells(adata, min_counts=10)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=3000)
sc.tl.pca(adata)
sc.pp.neighbors(adata, n_neighbors=15)
sc.tl.leiden(adata, resolution=0.4)

# Annotate by marker expression (same workflow as scRNA-seq)
sc.pl.dotplot(adata, markers_dict, groupby="leiden")

# Optional: if you want single-cell precision at 2µm bins,
# still use cell2location with matched scRNA-seq reference

Visium HD vs. Xenium: Choosing the Right Tool

Choose Visium HD when:

• You need whole-transcriptome coverage (not a fixed panel)
• Your question involves discovery — finding unexpected cell states or pathways
• You have FFPE tissue (Visium HD supports FFPE; Xenium also does)
• You want spatial multi-omics (Visium HD + CUT&TAG or protein co-detection)

Choose Xenium/CosMx when:

• You need subcellular resolution (mRNA localization within a cell)
• Your target genes are already known and a panel of 100–500 genes is sufficient
• You need higher sensitivity per transcript (in-situ amplification vs. sequencing)
• Your tissue is thick or 3D (Xenium handles z-stacking)

Squidpy Compatibility

Visium HD AnnData loads exactly like Visium v1 from Squidpy's perspective. All tools covered in this series work without modification:

# Everything from earlier posts works identically on Visium HD
sq.gr.spatial_neighbors(adata, coord_type="generic", n_neighs=6)
sq.gr.nhood_enrichment(adata, cluster_key="leiden")
sq.gr.ligrec(adata, n_perms=500, cluster_key="leiden")
sq.gr.spatial_autocorr(adata, mode="moran")

# The only practical difference: more bins = longer runtime
# nhood_enrichment on 30,000 bins takes ~10 min vs ~2 min for Visium v1

One Practical Warning: Tissue Coverage

Visium HD captures a 6.5 mm × 6.5 mm area — identical to Visium v1. This is enough for most tissue sections but smaller than a typical whole-slide image. If your question requires comparing multiple tumor regions across a large section, you will need multiple capture areas or switch to a platform with larger footprint (some CosMx runs cover full slides).

Interested in what the immune cell types you are mapping at single-cell resolution are actually doing in the tumor? NeoantigenLab covers T cell exhaustion, checkpoint inhibitor mechanisms, and CAR-T biology at the molecular level.

Get the complete pack

SpatiaBio Squidpy Pack 1 — 16 Notebooks ($19)

All 16 notebooks — neighborhood enrichment, co-occurrence, ligand-receptor, spatial autocorrelation, publication figures, and the complete workflow. Includes memory tips for Visium HD. Ready to run on your own data.

Get it for $19 →

---

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

Beyond Ligand-Receptor Scores: A Combined Decision Framework for Spatial Cell Communication

Spatial transcriptomics · ligand-receptor · cell communication · Squidpy

Why Spatial Context Matters for Cell Communication

Single-cell RNA-seq can tell you which ligands and receptors are expressed. It cannot tell you whether the cells expressing them are actually close enough to interact. Spatial transcriptomics closes that gap.

Two cells 500 µm apart in a tumor section are unlikely to exchange short-range signals regardless of what their transcriptomes predict. Squidpy's gr.ligrec and gr.nhood_enrichment functions let you ask both questions at once: what is being expressed, and are the right cells adjacent?

Figure 1. Left: ligand-receptor interaction scores between cell type pairs (* = p<0.05). Right: spatial co-occurrence z-scores — positive values indicate cell types that cluster together in tissue.

Setting Up the Interaction Database

Squidpy ships with a curated ligand-receptor database derived from CellChatDB and CellPhoneDB. You load it once and pass it to gr.ligrec:

import squidpy as sq
import scanpy as sc

# Load preprocessed Visium AnnData with cell type annotation in adata.obs["cell_type"]
adata = sc.read_h5ad("visium_annotated.h5ad")

# Build spatial graph first (required for proximity-aware analysis)
sq.gr.spatial_neighbors(adata, coord_type="grid", n_rings=2)

# Run ligand-receptor permutation test
sq.gr.ligrec(
    adata,
    n_perms=1000,
    cluster_key="cell_type",
    copy=False,
    use_raw=False,
    transmitter_params={"categories": "ligand"},
    receiver_params={"categories": "receptor"},
)
# Results stored in adata.uns["cell_type_ligrec"]

The permutation test shuffles cell type labels 1,000 times and asks: how often would this interaction score arise by chance given the spatial arrangement of cells? This is stricter than expression-only approaches because it accounts for the actual tissue layout.

Visualizing the Results

The raw output is a dictionary with means and pvalues DataFrames. Squidpy provides a dedicated dotplot:

sq.pl.ligrec(
    adata,
    cluster_key="cell_type",
    source_groups="Macrophages",
    target_groups=["T cells", "Cancer"],
    means_range=(0.3, float("inf")),  # filter weak interactions
    alpha=0.05,                        # significance threshold
    swap_axes=False,
    figsize=(10, 6),
)

Two practical filters here: means_range removes interactions with low absolute expression, and alpha applies the permutation p-value cutoff. Without both, the plot becomes unreadable with hundreds of LR pairs.

Neighborhood Enrichment: Are These Cells Actually Neighbors?

Before trusting any LR interaction, check whether the source and target cell types are spatially co-localized. gr.nhood_enrichment runs a permutation test on the spatial graph:

sq.gr.nhood_enrichment(adata, cluster_key="cell_type")
sq.pl.nhood_enrichment(
    adata,
    cluster_key="cell_type",
    method="average",
    cmap="RdBu_r",
    vmin=-3, vmax=3,
    figsize=(7, 6),
)

Z-scores above +2 indicate strong spatial co-occurrence; below −2 indicates mutual exclusion. From the figure: Macrophages and Cancer cells co-occur strongly (z = 3.1) — a known feature of immunosuppressive tumor microenvironments where TAMs cluster around tumor nests. T cells and Fibroblasts are mutually exclusive (z = −0.9), consistent with fibroblast-mediated immune exclusion.

Combining Both Analyses: A Decision Framework

Use neighborhood enrichment as a gate before interpreting LR results:

High enrichment + significant LR → strong candidate interaction worth validating (e.g., SPP1-CD44 between Macrophages and Cancer)

Low enrichment + significant LR → likely expression artifact; cells aren't actually proximal enough to interact

High enrichment + no significant LR → cells co-localize for structural reasons (e.g., Fibroblasts around vessels) rather than active signaling

Key LR Pairs in the Tumor Microenvironment

These consistently appear in Visium TME datasets:

• SPP1 (OPN) → CD44: Macrophage-to-cancer signaling, promotes invasion and immunosuppression
• CXCL12 → CXCR4: T cell and cancer cell chemoattraction; CXCR4 inhibition is a clinical target
• VEGFA → KDR (VEGFR2): Tumor-to-endothelial angiogenic signaling; highest scores near vessel-rich regions
• TGFB1 → TGFBR2: Fibroblast-mediated immunosuppression and ECM remodeling
• PDGFB → PDGFRB: Pericyte recruitment by endothelial cells; marks active neovascularization zones

Performance Notes

gr.ligrec with 1,000 permutations on a full Visium slide (~3,500 spots, 10+ cell types) takes 5–15 minutes on a standard laptop. Use n_jobs=-1 to parallelize across CPU cores. For exploratory work, drop to n_perms=200 first to check the output structure, then rerun with 1,000 for publication figures.

# Fast exploratory run
sq.gr.ligrec(adata, n_perms=200, cluster_key="cell_type", n_jobs=-1)

# Extract top interactions programmatically
import pandas as pd
means = adata.uns["cell_type_ligrec"]["means"]
pvals = adata.uns["cell_type_ligrec"]["pvalues"]

# Filter: mean > 0.5, p < 0.05 for Macrophage->Cancer axis
mac_cancer = means.loc[:, ("Macrophages", "Cancer")]
sig = mac_cancer[pvals.loc[:, ("Macrophages", "Cancer")] < 0.05]
print(sig.sort_values(ascending=False).head(10))

What's Next

Ligand-receptor analysis identifies candidate interactions. The next step is spatial trajectory analysis — reconstructing how gene expression changes as you move from one tissue region to another, which reveals which signals are active at each point along the tumor-stroma interface.

Working in tumor microenvironment biology? The companion blog NeoantigenLab covers the immunotherapy side — checkpoint inhibitor mechanisms, neoantigen prediction pipelines, and CAR-T biology for the same cell types you're mapping spatially.

Get the complete pack

SpatiaBio Squidpy Pack 1 — 16 Notebooks ($19)

All 16 notebooks — neighborhood enrichment, co-occurrence, ligand-receptor, spatial autocorrelation, publication figures, and the complete workflow. Includes memory tips for Visium HD. Ready to run on your own data.

Get it for $19 →

---

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