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.

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)