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.