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.

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.

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.