GPTomics/bio-spatial-transcriptomics-spatial-statistics

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| by chance | Report cell-type abundances; discount enrichment involving rare types |

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