bio-spatial-transcriptomics-spatial-statistics

name	bio-spatial-transcriptomics-spatial-statistics
description	Detects spatially variable genes, spatial autocorrelation, and cell-type colocalization for spatial transcriptomics using Squidpy with PySAL/esda for local statistics. Use when choosing an SVG method by its null and scaling (SpatialDE/SPARK GP variance-component vs SPARK-X/nnSVG linear vs Moran/Geary graph autocorrelation); separating genes that are spatially variable because of cell-type composition from genes regulated within a cell type; choosing the right autocorrelation statistic (global Moran/Geary vs Getis-Ord hot/cold spots vs local LISA and its FDR trap); and choosing a colocalization null strong enough to defeat the abundance/compartment confound (conditional or toroidal vs the weak Squidpy default permutation).
tool_type	python
primary_tool	squidpy

Version Compatibility

Reference examples tested with: squidpy 1.4+, scanpy 1.10+, anndata 0.10+, esda 2.5+, libpysal 4.9+ (SPARK, SPARK-X, and nnSVG are R/Bioconductor packages)

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Spatial Statistics

"Find spatially variable genes / run Moran's I / test which cell types colocalize" -> Quantify spatial structure in expression or in cell-type arrangement against an explicit null model.

• Python: squidpy.gr.spatial_autocorr (Moran/Geary), squidpy.gr.nhood_enrichment, squidpy.gr.co_occurrence, esda.Moran_Local/esda.getisord.G_Local (LISA, Getis-Ord)
• R: SPARK / SPARK-X / nnSVG (SVG), spdep (Moran/Geary/Getis-Ord/LISA)

Governing Principle

Two reframes decide whether a spatial-statistics result means anything. Both are silent failures: the code runs, the numbers look clean, and the interpretation is wrong.

A "spatially variable gene" is not necessarily spatially REGULATED. Moran's I, Geary's C, SPARK-X, and SpatialDE all detect that a gene's expression is spatially autocorrelated -- but a gene that is simply a marker of a spatially clustered cell type scores as "spatially variable" with zero cell-intrinsic spatial regulation, purely because the CELL TYPE is spatially organized. A hepatocyte gene in zonated liver, MAG in white matter, KRT17 in epithelium: all top SVGs, none of them regulated in space. This cell-type-driven signal swamps the genuinely interesting cell-type-INDEPENDENT signal (a gene graded across a niche WITHIN one cell type). Sample-wide SVG lists therefore largely re-derive marker genes and overlap heavily with HVGs -- if the SVG list is roughly the HVG list, the spatial test added almost nothing. The interesting question (within-type spatial regulation) needs cell-type-aware methods (C-SIDE, CTSV, CELINA/Celina), which are themselves unsettled and carry their own false positives. Decision: if the question is "where is tissue organized," sample-wide SVG is correct (the cell-type structure IS the answer); if the question is "which genes are regulated beyond cell identity," sample-wide SVG is the WRONG tool -- test within cell type or regress out composition first.

The null and the graph define the result. Every spatial statistic is computed on a neighbor graph (a weights matrix W) against a null distribution, and both are researcher choices, not properties of the tissue. Change kNN k from 6 to 30 and Moran's I, the enrichment z-scores, and the SVG ranking all move. Change the colocalization null from global label-shuffle to within-compartment and most "A is near B" claims evaporate. SVG methods disagree heavily across the literature precisely because each tests a DIFFERENT null (variance-component-zero vs covariance-independence vs graph-autocorrelation-zero) -- that disagreement is expected, not a bug. The honest workflow names the null, names the graph (cross-ref spatial-neighbors), and reports whether a hit survives a second graph or a stronger null.

Choosing an SVG Method

Goal: Pick a spatially-variable-gene test whose null hypothesis and computational scaling match the platform and the biological question.

Approach: Match GP variance-component methods to small Gaussian/count data, linear-scaling methods to single-cell-resolution data, and treat the SVG list as method-conditional -- cross-method intersection is more trustworthy than any single ranking, though it is small.

Method	Null it tests	Scaling	Best when	Fails when
SpatialDE (GP)	spatial variance component = 0 at the tested length scale	O(n^3); infeasible past ~1e4 locations	Small Visium-scale, Gaussian on log-normalized	Sparse/low counts violate Gaussian; fixed length-scale grid misses other scales
SPARK (count GLSM)	no pattern matching any of 10 fixed kernels	O(n^3)-ish (PQL); slow at large n	Small count data; want Poisson model, not Gaussian	Pattern unlike its 10 kernels; still cell-type-confounded
SPARK-X	expression covariance independent of location covariance	LINEAR in n and genes	1e4-1e6 cells (MERFISH/Xenium/CosMx) needing scalability	Fixed location kernels miss unusual length scales; low power on small focal hotspots
nnSVG (NNGP)	spatial variance = 0, with a per-gene length scale	LINEAR in n	Length scales genuinely differ across genes; large single-cell data	Still cell-type-confounded; more compute per gene than SPARK-X; needs adequate counts
Moran's I / Geary's C	no autocorrelation on graph W	Fast (sparse W)	Quick screen on an existing neighbor graph	Single fixed scale (the graph); misses multi-focal/small hotspots

There is no uniformly best SVG method; power is pattern-specific (SPARK-X, nnSVG, and Moran's I all have LOW power for genes high in small focal areas). Methods evolve fast -- verify the current benchmark before committing. Threshold on EFFECT SIZE (fraction of spatial variance), not p alone: with thousands of cells, trivial autocorrelation reaches tiny p-values.

Computing Spatial Autocorrelation with Squidpy

Goal: Rank genes by graph-based spatial autocorrelation as a fast SVG screen.

Approach: Build a neighbor graph, run Moran's I (or Geary's C) per gene with a permutation/analytic p-value and FDR, then read effect size before significance.

import squidpy as sq
import scanpy as sc

# The graph IS the model: k, coord_type, and units all change the result (see spatial-neighbors)
sq.gr.spatial_neighbors(adata, coord_type='generic', n_neighs=6)   # 6 mimics the Visium hex lattice

# genes=None defaults to highly_variable if present; n_perms adds a permutation null, corr_method applies FDR
sq.gr.spatial_autocorr(adata, mode='moran', n_perms=100, corr_method='fdr_bh')   # statsmodels name, not 'benjamini-hochberg'
moran = adata.uns['moranI']                  # columns: I, pval_norm, pval_norm_fdr_bh, ...

# Threshold on effect size (I) AND FDR, not p alone -- large n makes trivial autocorrelation 'significant'
svg = moran[(moran['I'] > 0.1) & (moran['pval_norm_fdr_bh'] < 0.05)].sort_values('I', ascending=False)

A top-ranked SVG here is a hypothesis about spatial structure, NOT evidence of spatial regulation. Before interpreting, compare the SVG list to the HVG list: the overlap is cell-type marker genes; the SVG-not-HVG subset (modest-amplitude gradients) is where spatial information actually lives.

Choosing an Autocorrelation Statistic

Goal: Match the statistic to the spatial question -- "is this gene clustered" is a different question from "where is it HIGH" and from "is THIS region a cluster."

Approach: Use a global statistic for one tissue-wide number, Getis-Ord when the sign (hot vs cold) matters, and local LISA for non-stationary tissue -- but pay the local multiple-testing tax correctly.

Statistic	Global / local	Hot vs cold?	Use when
Moran's I	global	NO (clustering of like values only)	One number: is this gene spatially structured across the whole section
Geary's C	global	NO	Same as Moran but more sensitive to LOCAL/short-range differences; disagreement with Moran is informative
Getis-Ord Gi*	local	YES -- separates high-clusters from low-clusters	"Where is this gene HIGH" -- hot/cold spot mapping
Local Moran / LISA	local	partial (HH/LL/HL/LH quadrants)	Non-stationary tissue: per-location clusters and spatial outliers

Moran's I and Geary's C cannot tell a hot spot from a cold spot -- both flag "similar values cluster" regardless of high or low. Choosing Moran when the question is "where is this gene HIGH" is a category error; use Getis-Ord Gi*. A non-significant GLOBAL Moran's I does NOT mean "no spatial structure": over heterogeneous tissue, positive autocorrelation in one region cancels negative in another, so use local statistics for non-stationary sections.

from esda.getisord import G_Local
from esda.moran import Moran_Local
from libpysal.weights import KNN

coords = adata.obsm['spatial']
w = KNN.from_array(coords, k=6)
w.transform = 'r'                            # row-standardized; changes the value AND its variance vs binary W

gene = adata[:, 'GENE1'].X.toarray().ravel()
gi = G_Local(gene, w, transform='B', star=True, permutations=999)   # star=True -> Gi* (includes self); binary weights for Getis-Ord
lisa = Moran_Local(gene, w, permutations=999)            # conditional-permutation local null
adata.obs['GENE1_hotspot'] = gi.Zs                       # positive Z = hot spot, negative = cold spot
adata.obs['GENE1_lisa_q'] = lisa.q                       # 1=HH, 2=LH, 3=LL, 4=HL

Local statistics carry a DOUBLE trap. There are n tests (one per location), so uncorrected LISA/Gi* maps are mostly false positives -- FDR is mandatory, and Anselin recommends stricter base cutoffs (0.01/0.005/0.001), not 0.05. Worse, the local statistics are themselves spatially autocorrelated (adjacent locations share neighbors, so adjacent I_i values are correlated), which violates the independence assumption of standard BH-FDR; the effective number of tests is far below n. Conditional permutation gives the correct local null but does not fix the cross-location dependence. Treat the cluster map as exploratory, not a set of independent discoveries.

Testing Cell-Type Colocalization

Goal: Decide whether two cell types are SPECIFICALLY associated in space, not merely both abundant or both in the same compartment.

Approach: Choose a permutation null strong enough to defeat the abundance/compartment confound; the Squidpy default answers only the weak question, and a co-occurrence distance profile is more informative than a single z-score.

Null model	What it permutes	Controls for	Misses
Global label permutation (Squidpy nhood_enrichment default)	all labels over all positions	graph topology, marginal counts	tissue compartmentalization -- two abundant co-compartment types pass trivially
Conditional / within-compartment permutation	labels within a region only	shared-compartment forcing	cross-compartment questions; the region choice is itself a degree of freedom
Toroidal shift	whole label field translated (wrapped)	each type's first-order density pattern	anisotropy; boundary realism (wrapping a bounded tissue is artificial)
Grid-tile shuffle across scales (CRAWDAD)	labels within tiles of size s	structure above scale s -- isolates colocalization AT scale s	within-tile structure below s

The single most common error in spatial-omics papers is reading a positive nhood_enrichment z-score as a specific A-B interaction. Under the weak global-permutation null it usually reflects co-compartmentalization plus abundance: two stromal populations, or tumor plus tumor-associated macrophages both in the tumor bed, pass trivially. A specific-affinity claim must SURVIVE a stronger null (conditional/within-compartment, toroidal shift preserving each type's density, or CRAWDAD's scale-explicit tiles). Rare-type enrichment is the least trustworthy: few edges give high-variance, often spuriously large |z| -- be most skeptical exactly where the biology is most exciting (a rare type near the tumor).

For clustering as a function of distance rather than a single graph z, Ripley's K/L (squidpy.gr.ripley, mode 'L') counts within-cluster neighbors within radius r against a complete-spatial-randomness expectation, so it reads as clustering vs dispersion ACROSS scale per cell type (squidpy computes the univariate L per cluster; a true bivariate cross-K answering "are A and B closer than chance, at what radius" needs a dedicated point-process tool such as spatstat). Its assumptions are the geostatistics ones tissue violates: a bounded, holey, non-stationary window. EDGE CORRECTION is mandatory and usually omitted -- without it, counts near the tissue boundary or a necrotic hole are biased DOWN and read as false depletion, so an ROI with gaps needs an edge-corrected estimator (or restrict analysis to the interior).

sq.gr.spatial_neighbors(adata, coord_type='generic', delaunay=True)   # nhood_enrichment reuses this stored graph

# Global label-permutation null -- the WEAK question: more adjacent than complete spatial randomness?
sq.gr.nhood_enrichment(adata, cluster_key='cell_type', n_perms=1000)
z = adata.uns['cell_type_nhood_enrichment']['zscore']

# co_occurrence gives a DISTANCE PROFILE (at what scale colocalization appears/vanishes), unlike a single z
sq.gr.co_occurrence(adata, cluster_key='cell_type')

Cellular Neighborhoods (Schurch/Nolan: cluster per-cell windows of neighbor composition) have NO inferential null at all -- they are descriptive k-means clusters whose number and window size are user knobs. They are useful summaries but routinely over-read as tested findings; the number of neighborhoods is chosen, not discovered, and identity shifts with window size. Test them downstream (neighborhood composition vs outcome), do not report them as significant in themselves.

Common Errors

Symptom	Cause	Fix
Top SVGs are all known cell-type markers; SVG list ~= HVG list	Sample-wide SVG re-derives markers of spatially clustered cell types (cell-type-driven, not regulated)	Ask within-type: regress out cell-type composition or use ctSVG (C-SIDE/CTSV/CELINA); report the SVG-not-HVG subset
Moran's I finds "clustering" but cannot locate where the gene is HIGH	Moran/Geary detect clustering of like values, blind to high vs low	Use Getis-Ord Gi* for hot/cold spots
LISA/Gi* map is mostly "significant"; thousands of hits	n tests, AND local statistics are spatially autocorrelated so naive BH-FDR is invalid	FDR with stricter cutoffs (0.001); conditional permutation; treat map as exploratory, not independent discoveries
Confident "cell type A interacts with B" that vanishes on a second look	Default nhood_enrichment global-permutation null only beats complete randomness; abundant co-compartment types pass trivially	Demand survival under a conditional/within-compartment or toroidal null; report abundances
Global Moran's I near zero but tissue is clearly structured	Non-stationarity: opposite-sign local regions cancel in one global number	Use local statistics (LISA/Gi*); stratify by region
SVG ranking changes completely between two runs/tools	Different graph (k, Delaunay vs kNN) or different null -- methods test different hypotheses	Name the graph and null; report graph-robust hits; expect cross-method intersection to be small
Radius/length-scale statistic gives nonsense	Coordinates in pixels/array units, parameter in microns; Visium array coords are not microns	Convert to microns via scale factors before any distance parameter (see spatial-neighbors)
Rare cell type shows a huge enrichment z-score	Few edges -> high-variance estimate -> large	z

Related Skills

• spatial-neighbors - builds the graph W that every statistic here inherits; the choice propagates
• spatial-domains - region-level structure; a domain is not a colocalization result
• spatial-communication - ligand-receptor in space; the spillover/abundance confounds recur there
• single-cell/markers-annotation - cell-type labels feeding colocalization, and the marker overlap that confounds SVG

References

• Svensson V, Teichmann SA, Stegle O (2018) SpatialDE: identification of spatially variable genes. Nature Methods 15(5):343-346. DOI 10.1038/nmeth.4636
• Sun S, Zhu J, Zhou X (2020) Statistical analysis of spatial expression patterns for spatially resolved transcriptomic studies (SPARK). Nature Methods 17(2):193-200. DOI 10.1038/s41592-019-0701-7
• Zhu J, Sun S, Zhou X (2021) SPARK-X: non-parametric modeling enables scalable and robust detection of spatial expression patterns for large spatial transcriptomic studies. Genome Biology 22:184. DOI 10.1186/s13059-021-02404-0
• Weber LM, Saha A, Datta A, Hansen KD, Hicks SC (2023) nnSVG for the scalable identification of spatially variable genes using nearest-neighbor Gaussian processes. Nature Communications 14:4059. DOI 10.1038/s41467-023-39748-z
• Palla G, Spitzer H, Klein M, et al. (2022) Squidpy: a scalable framework for spatial omics analysis. Nature Methods 19(2):171-178. DOI 10.1038/s41592-021-01358-2
• Schurch CM, Bhate SS, Barlow GL, et al. (2020) Coordinated cellular neighborhoods orchestrate antitumoral immunity at the colorectal cancer invasive front. Cell 182(5):1341-1359. DOI 10.1016/j.cell.2020.07.005
• Dos Santos Peixoto R, Miller BF, Brusko MA, et al. (2025) Characterizing cell-type spatial relationships across length scales in spatially resolved omics data (CRAWDAD). Nature Communications 16:350. DOI 10.1038/s41467-024-55700-1
• Moran PAP (1950) Notes on continuous stochastic phenomena. Biometrika 37(1/2):17-23. DOI 10.2307/2332142
• Geary RC (1954) The contiguity ratio and statistical mapping. The Incorporated Statistician 5(3):115-145. DOI 10.2307/2986645
• Getis A, Ord JK (1992) The analysis of spatial association by use of distance statistics. Geographical Analysis 24(3):189-206. DOI 10.1111/j.1538-4632.1992.tb00261.x
• Anselin L (1995) Local indicators of spatial association -- LISA. Geographical Analysis 27(2):93-115. DOI 10.1111/j.1538-4632.1995.tb00338.x