Showing posts with label R. Show all posts
Showing posts with label R. Show all posts

Tuesday, June 30, 2026

CellChat on Spatial Data: A Step-by-Step Tutorial (With Real Errors)

If you've already got a clustered Visium object from Seurat (like the one we built in our Squidpy vs Seurat comparison), the natural next question is: which cell types are actually talking to each other, and does that change with physical distance? That's exactly what CellChat's spatial mode answers — and unlike our earlier posts, this one comes with three real errors we hit along the way, not just the happy path.

What CellChat Does

CellChat infers, analyzes, and visualizes cell-cell communication networks from single-cell or spatial transcriptomics data, using a curated ligand-receptor interaction database. The spatial mode adds physical distance as a constraint — two cell types might both express a matching ligand-receptor pair, but if they're never actually near each other in the tissue, CellChat won't call that a likely real interaction.

We ran this on the exact same dataset from our Seurat post — the official stxBrain mouse brain Visium demo (2,696 spots, 15 clusters) — so this is a direct continuation, not a new dataset.

The Real Pipeline (With Real Errors)

library(Seurat)
library(SeuratData)
library(CellChat)

brain <- LoadData("stxBrain", type = "anterior1")
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)
brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30, verbose = FALSE)
brain <- FindClusters(brain, verbose = FALSE)

data.input <- GetAssayData(brain, layer = "data", assay = "SCT")
meta <- data.frame(labels = factor(paste0("C", Idents(brain))), row.names = names(Idents(brain)))
spatial.locs <- as.matrix(GetTissueCoordinates(brain)[, c("x","y")])

cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels",
                            datatype = "spatial", coordinates = spatial.locs,
                            spatial.factors = data.frame(ratio = 1, tol = 5))
cellchat@DB <- CellChatDB.mouse
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

cellchat <- computeCommunProb(cellchat, type = "truncatedMean", trim = 0.1,
                               distance.use = TRUE, interaction.range = 250,
                               scale.distance = 0.01, contact.dependent = TRUE, contact.range = 100)
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- aggregateNet(cellchat)

This looks clean here because we already fixed it. The first time through, we hit three separate real errors:

Error 1: GetAssayData() with slot is defunct

Using slot = "data" (which is what most older tutorials still show) throws a hard error in current SeuratObject — it was deprecated in 5.0.0 and is now fully removed. Fix: use layer = "data" instead. If you're following any Seurat tutorial written before 2024, check for this.

Error 2: CellChat rejects cluster label "0"

Error in setIdent(object, ident.use = group.by) :
  Cell labels cannot contain `0`!

Seurat's default cluster identities are numeric starting at 0 — completely normal and fine for Seurat itself. But CellChat's internal identity handling breaks if any cluster is literally named "0". Fix: prefix cluster labels with a letter before handing them to CellChat, e.g. factor(paste0("C", Idents(brain))) turns cluster 0 into "C0". Easy to miss because the error message doesn't make the cause obvious.

Error 3: missing presto dependency

CellChat's identifyOverExpressedGenes() wants the presto package for a faster Wilcoxon test implementation and throws an error (not just a warning) if it's missing — unlike many R packages that silently fall back. Fix: devtools::install_github('immunogenomics/presto') before running CellChat, or pass do.fast = FALSE if you don't want the extra dependency.

Real Results

On this dataset, the full pipeline (preprocessing + spatial-distance-constrained communication probability calculation) took ~360 seconds (6 minutes) on a standard machine — almost all of it in computeCommunProb(), which is the step that actually factors in spatial distance between every pair of spots.

  • 13,110 significant ligand-receptor interactions identified across the 15 spatial clusters
  • Top signaling pathways by interaction count: Glutamate (4,787), GABA-A (1,361), LAMININ (874), WNT (538), COLLAGEN (489), GABA-B (396)
  • This makes biological sense without any cherry-picking: glutamate and GABA are the two dominant excitatory/inhibitory neurotransmitter systems in brain tissue, so seeing them as the top communication signals in a mouse brain section is exactly what you'd expect from a sane result, not a red flag
CellChat circle plot showing cell-cell communication strength between all 15 clusters in mouse brain Visium tissue

Aggregated communication network across all signaling pathways — edge thickness reflects interaction strength between cluster pairs. Real output from the pipeline above.

CellChat circle plot showing Glutamate signaling specifically between clusters

Glutamate signaling specifically — the single largest signaling category in this dataset by interaction count.

Why the Spatial Constraint Matters

Without distance.use = TRUE, CellChat would just look at co-expression across the whole dataset, the same way it would for dissociated single-cell data — any cluster pair with matching ligand-receptor expression gets flagged, regardless of whether they're ever physically close in the tissue. With the spatial constraint on, interactions only count if the cell types are within the specified interaction.range (250 units here) or contact.range (100 units, for direct contact-dependent signaling). For a tissue like brain where spatial organization is functionally meaningful, this is the difference between a biologically plausible result and a list of statistically-matched-but-physically-impossible interactions.

Bottom Line

CellChat's spatial mode works, and the three errors above are exactly the kind of thing that wastes an afternoon if you don't know they're coming. If you're already running Seurat on Visium data, adding CellChat on top is a relatively small additional step — the real cost is the ~6 minute runtime for computeCommunProb() on a dataset this size, which will scale up on larger tissue sections.


Last updated: July 2026. Tested with CellChat (jinworks/CellChat, GitHub HEAD), Seurat 5.5.1, R 4.6.1, SeuratData stxBrain 0.1.2.

Squidpy vs Seurat for Spatial Transcriptomics: Which Tool Should You Use in 2026?

Introduction

If you're starting a spatial transcriptomics project in 2026, you've probably run into the same fork in the road every single-cell person eventually hits: Seurat or Squidpy? Both are mature, well-maintained, and backed by active communities — but they come from different ecosystems (R vs Python) and take different philosophical approaches to spatial data.

This isn't a "which one is objectively better" post. It's a practical breakdown of where each tool actually shines, based on what each is built for — so you can pick the right one for your pipeline instead of just going with whatever your lab used last.

To make this concrete instead of just theoretical, I actually ran Squidpy's core spatial workflow end-to-end on a real public dataset (squidpy's mouse brain Visium demo) and reported the real numbers below — not just textbook claims.

Quick Comparison

Seurat Squidpy
Language R Python
Core strength scRNA-seq + spatial integration Spatial-first, neighborhood/graph analysis
Best for Labs already in R/Bioconductor ecosystem Pipelines needing scanpy/AnnData compatibility
Spatial stats Basic (SpatialFeaturePlot, cluster overlay) Advanced (Moran's I, neighborhood enrichment, ligand-receptor via squidpy.gr)
Learning curve Lower if you know R Lower if you know scanpy
GitHub stars (as of 2026) ~2.8k ~580
Integration Seurat objects, SeuratDisk AnnData, works natively with scvi-tools

When Seurat Makes Sense

  • Your lab's existing scRNA-seq pipeline is already in Seurat — switching ecosystems mid-project has a real cost
  • You need quick, publication-ready spatial overlay plots without much custom code
  • You're working with 10x Visium data and want the well-documented Load10X_Spatial() workflow
library(Seurat)
library(SeuratData)

brain <- LoadData("stxBrain", type = "anterior1")
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)
brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30, verbose = FALSE)
brain <- FindClusters(brain, verbose = FALSE)
SpatialDimPlot(brain, label = TRUE)

Real output from this run (Seurat 5.5.1, R 4.6.1, SeuratData's stxBrain mouse brain Visium dataset — the official Seurat demo, conceptually the same tissue type as the squidpy demo below):

  • Dataset loaded: 31,053 genes × 2,696 spots
  • SCTransform took ~273 seconds (4.5 minutes) on this machine
  • Clustering (FindNeighbors + FindClusters) took only ~3.4 seconds after that — the bottleneck is entirely in SCTransform, not clustering
  • 15 clusters were identified, matching well with visually distinct anatomical regions in the tissue image
  • One real gotcha we hit: SCTransform printed a warning that it couldn't find the glmGamPoi package and fell back to a slower native implementation. Installing glmGamPoi via Bioconductor first is a meaningful speedup for anyone running this regularly — we didn't have it installed for this run, which is reflected in the ~273s figure above.
Seurat SpatialDimPlot showing 15 clusters identified in mouse brain Visium tissue

Seurat's spatial cluster plot on the official stxBrain demo dataset — real output from the code above, not a mockup.

(This workflow follows Seurat's official Visium tutorial structure — Seurat workflows for Visium data are extensively documented and battle-tested across thousands of published studies, which is itself part of Seurat's value: you're rarely the first person debugging a given step.)

When Squidpy Makes Sense

  • You need actual spatial statistics — neighborhood enrichment, co-occurrence, Moran's I autocorrelation — not just visualization
  • Your downstream pipeline is Python-based (e.g., feeding into scvi-tools or a custom ML model)
  • You're working with multiple spatial modalities (Visium, MERFISH, Xenium) and want one consistent API via AnnData
import squidpy as sq
import scanpy as sc

# Built-in public demo dataset (mouse brain, Visium)
adata = sq.datasets.visium_hne_adata()

sq.gr.spatial_neighbors(adata)
sq.gr.nhood_enrichment(adata, cluster_key="cluster")
sq.gr.co_occurrence(adata, cluster_key="cluster")
Spatial cluster distribution in mouse brain Visium dataset, generated with sq.pl.spatial_scatter

Spatial distribution of the 15 annotated brain regions, plotted directly from the squidpy run above.

Real output from this run (squidpy 1.x, mouse brain Visium demo dataset, 2,688 spots × 18,078 genes, 15 annotated regions including Cortex_1–5, Hippocampus, Thalamus_1–2, Striatum, etc.):

  • spatial_neighbors + nhood_enrichment completed in ~9.5 seconds on the full 2,688-spot dataset
  • co_occurrence completed in ~6.7 seconds
  • Neighborhood enrichment z-scores ranged up to ~78, meaning some region pairs (e.g. cortical layers next to each other) show very strong spatial co-localization — exactly the kind of signal Seurat's basic spatial plots won't quantify for you directly
  • One practical gotcha I hit: sq.gr.spatial_neighbors throws a FutureWarning as of squidpy 1.x — it's being split into spatial_neighbors_knn / _radius / _delaunay / _grid in v1.9. If you're starting a new project now, it's worth using the specific function instead of the generic one to avoid a breaking change later.
Neighborhood enrichment z-score heatmap between brain regions, generated with sq.pl.nhood_enrichment

Neighborhood enrichment z-scores between region pairs — darker/brighter cells indicate stronger spatial co-localization.

(Tested on squidpy's public demo dataset, not proprietary research data — same workflow applies to your own Visium/MERFISH/Xenium data via sq.read.visium() or the equivalent reader.)

The Honest Answer

If you're already deep in the R/Bioconductor world and just need solid visualization + basic spatial clustering, Seurat won't slow you down. If your work leans toward quantitative spatial statistics or you're building anything that touches a Python ML stack, Squidpy's purpose-built spatial functions save you from reinventing the wheel.

Based on actually running the Squidpy workflow above, the value isn't in the plots — it's in the quantified output (z-scores, co-occurrence stats) you can directly use in a paper or downstream model. That's the gap Seurat's spatial functions don't close on their own.

FAQ

Can I use both in the same project? Yes — many labs preprocess in Seurat, export via SeuratDisk to AnnData, then run Squidpy's spatial stats.

Does Squidpy work with non-10x platforms like MERFISH or Xenium? Yes, Squidpy has built-in readers for multiple spatial platforms, which is one of its main advantages over Seurat's more Visium-centric workflow.

Is one faster than the other for large datasets? We ran both on a similar-sized dataset (~2,700 spots) on the same machine. Squidpy's spatial graph + neighborhood enrichment finished in under 10 seconds. Seurat's SCTransform alone took ~273 seconds, though clustering after that was fast (~3.4s). Important caveat: these aren't doing identical work — SCTransform is a normalization step with no direct Squidpy equivalent in our comparison, so this isn't a clean apples-to-apples benchmark, just real numbers from real runs of each tool's typical workflow. Installing glmGamPoi would likely cut Seurat's SCTransform time significantly; we ran without it.


Last updated: June 2026. Tool versions referenced: Seurat 5.x, Squidpy 1.x — always check current documentation as APIs evolve.

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 G...