<!--

COPYRIGHT NOTICE

This file is part of the "Universal Biomedical Skills" project.

Copyright (c) 2026 MD BABU MIA, PhD <md.babu.mia@mssm.edu>

All Rights Reserved.

#

This code is proprietary and confidential.

Unauthorized copying of this file, via any medium is strictly prohibited.

#

Provenance: Authenticated by MD BABU MIA

-->

Spatial Statistics

Compute spatial statistics and identify spatially variable features.

Required Imports

python
import squidpy as sq
import scanpy as sc
import pandas as pd
import numpy as np

Compute Spatial Autocorrelation (Moran's I)

python
# Requires spatial neighbors
sq.gr.spatial_neighbors(adata, coord_type='generic', n_neighs=6)

# Compute Moran's I for all genes (can be slow)
sq.gr.spatial_autocorr(adata, mode='moran')

# Or for specific genes
sq.gr.spatial_autocorr(adata, mode='moran', genes=['GENE1', 'GENE2', 'GENE3'])

# Results stored in adata.uns['moranI']
moran_results = adata.uns['moranI']
print(moran_results.head(20))

Interpret Moran's I

python
# Moran's I ranges from -1 to 1
# I > 0: positive spatial autocorrelation (similar values cluster)
# I = 0: random spatial distribution
# I < 0: negative spatial autocorrelation (dissimilar values cluster)

# Get significantly spatially variable genes
svg = moran_results[moran_results['pval_norm'] < 0.05].sort_values('I', ascending=False)
print(f'Found {len(svg)} spatially variable genes (p < 0.05)')
print('\nTop 10 spatially variable genes:')
print(svg.head(10)[['I', 'pval_norm']])

Compute Geary's C

python
# Alternative spatial autocorrelation measure
sq.gr.spatial_autocorr(adata, mode='geary')

# Results in adata.uns['gearyC']
geary_results = adata.uns['gearyC']
# C < 1: positive spatial autocorrelation
# C = 1: random
# C > 1: negative spatial autocorrelation

Co-occurrence Analysis

python
# Analyze co-localization of cell types/clusters
# First, ensure you have cluster labels
sc.pp.neighbors(adata)
sc.tl.leiden(adata)

# Compute co-occurrence
sq.gr.co_occurrence(adata, cluster_key='leiden')

# Results in adata.uns['leiden_co_occurrence']
# Visualize co-occurrence
sq.pl.co_occurrence(adata, cluster_key='leiden')

Interpret Co-occurrence

python
co_occ = adata.uns['leiden_co_occurrence']
occ_matrix = co_occ['occ']  # Occurrence matrix
interval = co_occ['interval']  # Distance intervals

# occ_matrix[i, j, k] = occurrence of cluster j around cluster i at distance interval k
print(f'Occurrence matrix shape: {occ_matrix.shape}')
print(f'Distance intervals: {interval}')

Neighborhood Enrichment

python
# Test if clusters are enriched in each other's neighborhoods
sq.gr.nhood_enrichment(adata, cluster_key='leiden')

# Results in adata.uns['leiden_nhood_enrichment']
# zscore > 0: clusters co-localize more than expected
# zscore < 0: clusters avoid each other

# Visualize
sq.pl.nhood_enrichment(adata, cluster_key='leiden')

Extract Enrichment Z-scores

python
enrichment = adata.uns['leiden_nhood_enrichment']
zscore = enrichment['zscore']
clusters = adata.obs['leiden'].cat.categories

# Convert to DataFrame
zscore_df = pd.DataFrame(zscore, index=clusters, columns=clusters)
print('Neighborhood enrichment z-scores:')
print(zscore_df)

Ripley's Statistics

python
# Ripley's K/L function for point pattern analysis (single-cell resolution data)
sq.gr.ripley(adata, cluster_key='leiden', mode='L')

# Results in adata.uns['leiden_ripley']
sq.pl.ripley(adata, cluster_key='leiden')

Centrality Scores

python
# Compute centrality of each cell type
sq.gr.centrality_scores(adata, cluster_key='leiden')

# Results in adata.uns['leiden_centrality_scores']
centrality = adata.uns['leiden_centrality_scores']
print(centrality)

Interaction Matrix

python
# Build interaction matrix between clusters
sq.gr.interaction_matrix(adata, cluster_key='leiden')

# Results in adata.uns['leiden_interactions']
interactions = adata.uns['leiden_interactions']
print(interactions)

Custom Spatial Statistic

python
from scipy.stats import pearsonr

def spatial_correlation(adata, gene1, gene2):
    '''Compute spatial correlation between two genes'''
    expr1 = adata[:, gene1].X.toarray().flatten()
    expr2 = adata[:, gene2].X.toarray().flatten()
    r, p = pearsonr(expr1, expr2)
    return r, p

r, p = spatial_correlation(adata, 'GENE1', 'GENE2')
print(f'Spatial correlation: r={r:.3f}, p={p:.2e}')

Local Moran's I (LISA)

python
from esda.moran import Moran_Local
from libpysal.weights import KNN

# Build weights matrix
coords = adata.obsm['spatial']
w = KNN.from_array(coords, k=6)
w.transform = 'r'

# Compute local Moran's I for a gene
gene_expr = adata[:, 'GENE1'].X.toarray().flatten()
lisa = Moran_Local(gene_expr, w)

# Add to adata
adata.obs['GENE1_lisa'] = lisa.Is
adata.obs['GENE1_lisa_q'] = lisa.q  # Quadrant (HH, HL, LH, LL)

Batch Spatial Statistics

python
# Compute Moran's I for top variable genes only
hvg = adata.var_names[adata.var['highly_variable']][:500]
sq.gr.spatial_autocorr(adata, mode='moran', genes=hvg)

results = adata.uns['moranI']
significant = results[results['pval_norm'] < 0.01]
print(f'{len(significant)} genes with significant spatial autocorrelation')

Related Skills

• spatial-neighbors - Build spatial graphs (prerequisite)
• spatial-domains - Identify spatial domains
• spatial-visualization - Visualize spatial statistics

<!-- AUTHOR_SIGNATURE: 9a7f3c2e-MD-BABU-MIA-2026-MSSM-SECURE -->