Complete Squidpy Spatial Transcriptomics Workflow: From SVGs to Cell Communication

This post wraps up the SpatiaBio Squidpy series. Over ten posts we've built a complete spatial transcriptomics analysis pipeline from scratch — from loading raw Visium data to identifying cell-cell communication patterns. Here's the full workflow in one place.

The complete pipeline

import squidpy as sq
import scanpy as sc

# ── 1. Load data ──────────────────────────────────────────────────────────────
adata = sq.datasets.visium_hne_adata()

# ── 2. Preprocessing ──────────────────────────────────────────────────────────
# NOTE: do NOT normalize again. visium_hne_adata() ships .X ALREADY log-normalized
# (max ~8.4, non-integer); the raw integer counts live in adata.raw.X.
# normalize_total + log1p here would double-normalize — silently, and it changes
# every result below, including which genes come out as top SVGs.
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="cell_ranger")
sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.leiden(adata, key_added="cluster")   # pip install leidenalg

# ── 3. Spatially variable genes ───────────────────────────────────────────────
sq.gr.spatial_neighbors(adata)
sq.gr.spatial_autocorr(adata, mode="moran")          # post 5: Moran's I
sq.pl.spatial_scatter(adata, color="Prkcd", use_raw=False)   # post 6: SVG visualization
#   use_raw=False matters: spatial_scatter defaults to use_raw=True and would plot
#   raw counts, which is not what Moran's I was computed on.

# ── 4. Spatial neighborhood analysis ─────────────────────────────────────────
sq.gr.nhood_enrichment(adata, cluster_key="cluster")  # post 7
sq.pl.nhood_enrichment(adata, cluster_key="cluster")

sq.gr.co_occurrence(adata, cluster_key="cluster")     # post 8
sq.pl.co_occurrence(adata, cluster_key="cluster", clusters=["0","1","2"])

# ── 5. Cell-cell communication ────────────────────────────────────────────────
sq.gr.ligrec(adata, cluster_key="cluster", n_perms=1000)  # post 9
sq.pl.ligrec(adata, cluster_key="cluster", source_groups=["0","1","2"])

What the pipeline above actually produces, run end-to-end on the built-in mouse-brain section. 1 the 15 annotated regions that define the units of analysis; 2 the top spatially variable gene among 500 HVGs — Prkcd (Moran's I = 0.824), localising cleanly to thalamus; 3 neighborhood enrichment with the self-enrichment diagonal masked, which is what makes the real adjacencies visible (Cortex_4–Cortex_5, Hippocampus–pyramidal layers); 4 co-occurrence, showing how Cortex_1's partners fall off with distance.

What each step tells you

Step	Function	Question answered
SVG detection	sq.gr.spatial_autocorr	Which genes vary across space?
SVG visualization	sq.pl.spatial_scatter	Where exactly do SVGs express?
Neighborhood enrichment	sq.gr.nhood_enrichment	Which clusters physically border each other?
Co-occurrence	sq.gr.co_occurrence	At what distance do clusters start co-occurring?
Ligand-receptor	sq.gr.ligrec	Which cells are actively signaling to each other?

How the analyses connect

These aren't independent analyses — they form a natural progression:

1. Clustering defines the units of analysis (Leiden clusters = putative cell types or regions)
2. SVGs tell you which genes drive spatial structure
3. Neighborhood enrichment tells you which clusters are physically adjacent at the spot scale (~55 µm)
4. Co-occurrence extends that to a range of distances, revealing spatial gradients
5. Ligand-receptor asks: given that two clusters are adjacent, are they actually communicating?

A common workflow is to identify spatially adjacent cluster pairs from neighborhood enrichment, then specifically test those pairs for ligand-receptor interactions — focusing statistical power where anatomy suggests contact.

Common pitfalls

• Normalizing data that is already normalized — this one throws no error at all. visium_hne_adata() ships .X log-normalized, so a reflexive normalize_total + log1p logs it twice and quietly changes your results — including which genes rank as top SVGs. Check adata.X.max(): single digits and non-integer means it is already logged.
• Skipping sq.gr.spatial_neighbors before neighborhood/co-occurrence analyses — both require the spatial graph to be built first.
• Running sc.tl.leiden without leidenalg — install it separately: pip install leidenalg.
• Co-occurrence on too many clusters — it scales as O(n²) per distance bin. Subset to clusters of interest if slow.
• Interpreting ligrec p-values without a mean threshold — low p-value + near-zero mean is noise. Always filter by means_range.
• Treating Visium clusters as single cell types — each 55 µm spot contains ~5–20 cells. Cluster labels are regional, not single-cell.

What's next

The Squidpy ecosystem continues to expand. Upcoming areas worth watching:

• Image features (sq.im) — extract morphological features from H&E images alongside transcriptomics
• Deconvolution (cell2location, RCTD) — resolve cell type mixtures within spots
• Multi-sample integration — comparing spatial patterns across conditions or replicates
• Sub-spot resolution (10x Xenium, MERFISH) — transcript-level spatial data where each cell is individually resolved

Full series notebooks on GitHub

Every analysis in this series — with runnable code and real outputs — is available at the link below.

github.com/Lociven/spatiabio-tutorials

This walkthrough, plus 15 more

Every step above as a standalone notebook

If this end-to-end version was useful, Pack 1 has each stage broken out separately with more depth — plus setup, memory handling, and error debugging.

Get Pack 1 for $19 →

---

Data: squidpy built-in visium_hne_adata (10x Genomics Visium, mouse brain H&E)