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.

Cell2location: Deconvolving Visium Spots to Single-Cell Resolution

By Lociven · SpatiaBio · August 2026

Visium spots are not single cells. Each spot captures RNA from multiple cells, and without knowing which cell types contributed what, your spatial analysis is blurry. Cell2location solves this by deconvolving Visium spots using a scRNA-seq reference. This post explains how it works and how to run it.

---

The problem: spots contain multiple cells

A standard Visium v2 spot is 55 µm in diameter. In most human tissues, that area contains 2–15 cells depending on cell size and packing density. When you cluster spots and call marker genes, you're clustering mixtures — and the "spatial domains" you find may reflect the composition of those mixtures rather than true biological boundaries.

Cell type deconvolution maps each spot to a probability distribution over cell types: "this spot is 60% cancer epithelial, 25% macrophage, 15% fibroblast." That's what cell2location produces.

Why cell2location over other methods? RCTD and NNLS are faster but assume cell type proportions sum to 1 (full tissue coverage). Cell2location uses a hierarchical Bayesian model that accounts for differences in RNA capture efficiency between scRNA-seq and Visium — which matters a lot in practice.

---

What you need

• Visium AnnData — raw counts (not normalized), spatial coordinates
• scRNA-seq reference — same tissue type, annotated cell types. Does not need to be the same sample.
• GPU strongly recommended (CPU works but is slow for large datasets)

The scRNA-seq reference is the critical input. If your reference has poorly annotated cell types, or doesn't cover all cell types present in the tissue, the deconvolution will be wrong. More cell types in the reference is generally better, as long as they're well-annotated.

---

Step 1: Prepare the scRNA-seq reference

import cell2location
import scanpy as sc
import numpy as np

# Load scRNA-seq reference
adata_ref = sc.read_h5ad("scrna_reference.h5ad")

# Keep raw counts
adata_ref.X = adata_ref.layers["counts"].copy()

# Filter to cell types with enough cells
cell_type_counts = adata_ref.obs["cell_type"].value_counts()
keep = cell_type_counts[cell_type_counts >= 5].index
adata_ref = adata_ref[adata_ref.obs["cell_type"].isin(keep)].copy()

# Filter genes
sc.pp.filter_genes(adata_ref, min_cells=1)

print(f"Reference: {adata_ref.n_obs} cells, {adata_ref.n_vars} genes")
print(f"Cell types: {adata_ref.obs['cell_type'].nunique()}")

---

Step 2: Train the reference model

from cell2location.models import RegressionModel

# Set up the regression model on the reference
RegressionModel.setup_anndata(
    adata=adata_ref,
    batch_key="sample",        # batch correction if multiple samples
    labels_key="cell_type",
    layer="counts"
)

mod = RegressionModel(adata_ref)
mod.train(max_epochs=250, use_gpu=True)

# Export the estimated expression signatures
adata_ref = mod.export_posterior(
    adata_ref,
    sample_kwargs={"num_samples": 1000, "batch_size": 2500, "use_gpu": True}
)

# Save the cell type signatures
inf_aver = adata_ref.varm["means_per_cluster_mu_fg"][
    [f"means_per_cluster_mu_fg_{i}" for i in adata_ref.uns["mod"]["factor_names"]]
].copy()
inf_aver.columns = adata_ref.uns["mod"]["factor_names"]
inf_aver.iloc[0:5, 0:5]

---

Step 3: Run deconvolution on Visium

import squidpy as sq

# Load Visium data
adata_vis = sq.read.visium("path/to/visium/")
adata_vis.var_names_make_unique()

# Keep only genes present in both datasets
intersect = np.intersect1d(adata_vis.var_names, inf_aver.index)
adata_vis = adata_vis[:, intersect].copy()
inf_aver = inf_aver.loc[intersect, :]

# Set up cell2location model on Visium
cell2location.models.Cell2location.setup_anndata(
    adata=adata_vis, layer="counts"
)

mod_vis = cell2location.models.Cell2location(
    adata_vis,
    cell_state_df=inf_aver,
    N_cells_per_location=10,    # expected cells per spot
    detection_alpha=20
)

mod_vis.train(
    max_epochs=30000,
    batch_size=None,
    train_size=1,
    use_gpu=True
)

---

Step 4: Extract and visualize cell type abundances

# Export posterior
adata_vis = mod_vis.export_posterior(
    adata_vis,
    sample_kwargs={"num_samples": 1000, "batch_size": 10, "use_gpu": True}
)

# Add 5% quantile cell type abundances to obsm
adata_vis.obs[adata_vis.uns["mod"]["factor_names"]] = adata_vis.obsm["q05_cell_abundance_w_sf"]

# Plot spatial distribution of each cell type
import matplotlib.pyplot as plt

cell_types_to_plot = ["Cancer_epithelial", "CD8_T_cell", "Macrophage", "Fibroblast"]

fig, axes = plt.subplots(1, len(cell_types_to_plot), figsize=(16, 4))
for ax, ct in zip(axes, cell_types_to_plot):
    sq.pl.spatial_scatter(
        adata_vis, color=ct, ax=ax,
        cmap="magma", size=1.5, title=ct
    )
plt.tight_layout()
plt.savefig("cell2location_output.png", dpi=150, bbox_inches="tight")

---

How deconvolution changes what a spot means. Left: a 55 µm Visium spot physically contains several cells of different types, and sequencing returns one averaged profile for all of them. Right: cell2location returns a composition instead — the worked example from this post. Read the q05 column, the conservative 5th-percentile estimate, rather than the mean.

Interpreting the output

The key output column is q05_cell_abundance_w_sf — the 5th percentile of the posterior distribution of cell type abundance per spot. Using the 5th percentile (rather than the mean) gives conservative, high-confidence estimates.

Common patterns to look for:

• Cancer cells concentrated in the tumor core, stromal cells at the margin
• CD8+ T cells enriched at the invasive margin vs excluded from the tumor core (a key immunotherapy biomarker)
• Macrophage subtypes (M1/M2) in different spatial zones
• Spots with very low total abundance — these may be necrotic areas or technical artifacts

Common mistake: Setting N_cells_per_location too high or too low drastically affects results. Use histology or DAPI staining to estimate the actual cell density in your tissue before running the model.

---

Using cell2location output with Squidpy

Once you have cell type abundances per spot, you can feed them directly into Squidpy's neighborhood enrichment and co-occurrence functions — treating cell types as categorical annotations rather than gene expression clusters.

import squidpy as sq

# Assign dominant cell type per spot
adata_vis.obs["dominant_cell_type"] = (
    adata_vis.obs[adata_vis.uns["mod"]["factor_names"]]
    .idxmax(axis=1)
)

# Now run neighborhood enrichment on deconvolved cell types
sq.gr.spatial_neighbors(adata_vis, coord_type="grid", n_rings=2)
sq.gr.nhood_enrichment(adata_vis, cluster_key="dominant_cell_type")
sq.pl.nhood_enrichment(adata_vis, cluster_key="dominant_cell_type", method="ward")

---

Get the complete pack

SpatiaBio Pack 2 — Deconvolution & Spatial Pattern Detection

The full cell2location deconvolution notebook — real scRNA-seq reference mapping, tissue zone clustering, and spatial neighborhood enrichment between zones. Plus a SpatialDE2 SVG detection notebook.

Get it for $19 →

From the NeoantigenLab sister blog

Applying spatial transcriptomics to tumor immunology?

NeoantigenLab covers the biology behind what you're profiling — neoantigen pipelines, TIL therapy, checkpoint inhibitors, and cancer immunotherapy for experimental researchers.

Visit NeoantigenLab →

Tags: cell2location, cell type deconvolution, Visium, spatial transcriptomics, scRNA-seq, RCTD, Bayesian deconvolution, tumor microenvironment, Squidpy