Spatially Variable Genes with Squidpy: Moran's I on Mouse Brain Visium Data

One of the first questions you ask after clustering a Visium dataset is: which genes actually vary across space? Not just between clusters — but continuously, across the tissue. That's what spatially variable gene (SVG) analysis answers.

In this post we run Moran's I via Squidpy on the 10x Genomics mouse brain Visium demo dataset and look at what comes out.

What is Moran's I?

Moran's I is a spatial autocorrelation statistic. For a given gene, it asks: are spots with high expression near other spots with high expression? A score near 1 means strong spatial clustering. Near 0 means random. Near -1 means a checkerboard pattern (rare in transcriptomics).

The formula is:

I = (N / W) * (Σᵢ Σⱼ wᵢⱼ(xᵢ - x̄)(xⱼ - x̄)) / Σᵢ(xᵢ - x̄)²

Where N is the number of spots, wᵢⱼ is the spatial weight between spots i and j (1 if neighbors, 0 otherwise), and x is normalized expression.

Squidpy wraps this with permutation-based p-values so you get significance alongside the score.

The setup

Dataset: visium_hne_adata() from Squidpy's built-in demos — 2,688 spots, 18,078 genes, mouse brain coronal section.

We filtered to the top 500 highly variable genes (HVGs) before running Moran's I. Running on all 18k genes is possible but slow; HVG filtering keeps the analysis focused on genes that vary meaningfully across the dataset.

import squidpy as sq
import scanpy as sc

adata = sc.read_h5ad("visium_hne_adata.h5ad")

# Filter to top 500 HVGs
sc.pp.highly_variable_genes(adata, n_top_genes=500, flavor="cell_ranger")
adata = adata[:, adata.var.highly_variable].copy()

# NOTE: do NOT normalize again. This dataset ships .X already log-normalized
# (max ~8.9, non-integer). Raw counts live in adata.raw.X.
# Running normalize_total + log1p here double-normalizes and changes the results.

# Build spatial graph
sq.gr.spatial_neighbors(adata)

# Run Moran's I with permutation testing
if __name__ == '__main__':
    sq.gr.spatial_autocorr(adata, mode="moran", n_perms=200, n_jobs=1)

svg_results = adata.uns["moranI"]
top20 = svg_results.sort_values("I", ascending=False).head(20)

One note on if __name__ == '__main__': — this guard is required on Windows because Squidpy's permutation step uses multiprocessing. Without it you get a RuntimeError about the bootstrapping phase. On Linux/Mac it's not needed.

Results: top 20 SVGs

Of the 500 genes tested, 486 have an FDR-corrected p-value indistinguishable from zero under the normal approximation (pval_norm_fdr_bh); the largest is 0.027. Note that the permutation-based p-values (pval_sim) can never be exactly zero — they bottom out at roughly 1/(n_perms+1), so with 200 permutations the floor is ~0.005. Either way the ranking by I score is what carries the information here.

Rank	Gene	Moran's I	FDR p (norm)
1	Prkcd	0.824	0.0
2	Pmch	0.818	0.0
3	Mobp	0.811	0.0
4	Agrp	0.799	0.0
5	Tcf7l2	0.794	0.0
6	Agt	0.759	0.0
7	Tnnt1	0.753	0.0
8	Hcrt	0.747	0.0
9	Baiap3	0.744	0.0
10	Gal	0.737	0.0
11	Itpka	0.727	0.0
12	Gpx3	0.712	0.0
13	Gm5741	0.710	0.0
14	Trf	0.705	0.0
15	Fezf1	0.698	0.0
16	Stx1a	0.689	0.0
17	Ctxn3	0.680	0.0
18	Mal	0.680	0.0
19	Cldn11	0.674	0.0
20	Calb2	0.673	0.0

Figure 1. What a high vs. low Moran's I actually looks like on tissue. Left: Prkcd (I = 0.824), the top-ranked gene — expression forms one compact, contiguous territory (the thalamus). Right: Gzma (I = 0.022) — expressed at similar sparsity but scattered at random, with no spatial structure. Same dataset, same spatial graph; only the spatial arrangement differs.

What the results tell us

The striking thing about the corrected ranking is that it isn't a grab-bag — the top hits fall into three anatomically coherent compartments, which is exactly what you'd hope a spatial statistic would find in a coronal section:

Thalamus — Prkcd (rank 1), Tcf7l2 (5). Both are well-established thalamic markers, and the thalamus is a large, spatially contiguous block in this section. That's the ideal shape for a high Moran's I: a compact, on/off territory.

Hypothalamic neuropeptides — Pmch (2), Agrp (4), Hcrt (8), Gal (10). These are classic hypothalamic neuropeptides, expressed in tight, well-defined nuclei. A gene that is essentially silent everywhere except one small dense region scores extremely well here.

Myelin / white matter — Mobp (3), Trf (14), Mal (18), Cldn11 (19). Oligodendrocyte and myelin genes trace the fiber tracts, which are spatially continuous structures. Notably, this is a non-neuronal signal — Moran's I doesn't care whether the pattern comes from neurons; it only cares that the pattern is spatially coherent.

The through-line: Moran's I rewards compact, contiguous territories — not "importance." A gene restricted to one anatomical nucleus outranks a broadly-expressed cortical layer marker, because the layer marker's signal is spread across a thinner, more interdigitated band. That's a property of the statistic, not a statement about biology.

(The anatomical grouping above is our read of the gene list, not an output of the analysis. The I values and p-values in the table are computed; the compartment labels are interpretation.)

A note on HVG pre-filtering

We filtered to 500 HVGs before running Moran's I. This is a practical choice, not a methodological requirement. Running on all genes would likely surface similar top hits but would take proportionally longer. The tradeoff: genes that are spatially variable but not highly variable (e.g., lowly expressed region-specific markers) could be missed. For a complete SVG screen, running on all expressed genes with a minimum count filter is more thorough.

The trap: this dataset is already normalized

Here is the mistake that is very easy to make, and that this post originally made: visium_hne_adata() ships with .X already log-normalized. You can check in one line — the values top out around 8.9 and are not integers, while the untouched counts sit in adata.raw.X (max ~23,700):

a = sq.datasets.visium_hne_adata()
a.X.max()        # ~8.86, non-integer  -> already log-normalized
a.raw.X.max()    # ~23703              -> these are the raw counts

If you follow the usual muscle memory and run normalize_total + log1p anyway, you are normalizing normalized data and logging it twice. It does not error. It just quietly gives you different answers:

Rank	Double-normalized (wrong)	Correct
1	Itpka 0.674	Prkcd 0.824
2	Fezf1 0.634	Pmch 0.818
3	Baiap3 0.622	Mobp 0.811

Different genes at the top, and every I value is depressed — the extra normalization squashes the dynamic range the statistic depends on. Itpka itself doesn't vanish; it moves from rank 1 (0.674) to rank 11 (0.727). Nothing warns you.

The general lesson is worth more than the specific fix: before you normalize, check whether your object already is. Convenience datasets and processed .h5ad files from collaborators frequently arrive pre-normalized, and a "runs fine, results look plausible" outcome is the worst kind of failure.

What's next

The natural follow-up is to visualize where these top SVGs are expressed on the tissue — Squidpy's sq.pl.spatial_scatter does this in a few lines. That'll be the next post.

The full code for this analysis is at github.com/Lociven/spatiabio-tutorials.

Next step from here

Moran's I is one piece of the pipeline

Pack 1 has the SVG visualization notebook that follows this one, plus 14 others — clustering through publication figures, all outputs included.

See what's in Pack 1 →