Spatial Domain Detection

Spatial Domain Detection

Overview

Split a tissue section into spatially coherent regions (cortical layers, tumour core vs stroma, germinal centres) from expression plus location. Each method returns one domain label per spot or cell. Output: domain labels from one or more methods, their agreement, domain marker genes and domain maps.

Domains are clusters, not ground truth: run two methods and report where they agree.

When to Use

• Dividing a Visium or ST section into tissue regions such as layers, tumour core and stroma
• Using the H&E image together with expression to draw region boundaries (SpaGCN)
• Finding cellular niches in Xenium, MERFISH, CosMx or IMC data without labels (UTAG)
• Getting region labels to summarise deconvolved compositions or cell types by
• Checking whether two independent methods agree on the main boundaries
• Use Step 7 of spatial-transcriptomics-cell-type-annotation instead when cells are already annotated and niches should be defined by cell type mix
• For neighbourhood enrichment or spatially variable genes, use spatial-transcriptomics-analysis-squidpy

Prerequisites

• Input: raw counts in X (or layers["counts"]) and coordinates in obsm["spatial"]
• SpaGCN: pulls torch. Its Louvain init needs louvain, which has no Python 3.13 conda build or wheel and cannot be compiled here; without it use the k-means init (Step 3)
• GraphST: pulls torch; a CUDA GPU makes training ~10x faster; the dense spot-by-spot graph limits it to ~30k observations
• UTAG: not on PyPI; install from GitHub

bashCopy code
# Install only what `import` reports missing. `--pypi` keeps the env's scanpy/numba
# stack; a conda add can swap it and bare pip can shadow it.
pixi add --pypi scanpy SpaGCN GraphST
pixi add --pypi "utag @ git+https://github.com/ElementoLab/utag.git@main"

Workflow

1. Load Visium data (or your own)
2. Preprocess
3. SpaGCN domains (spots, with histology)
4. GraphST domains (spots or cells)
5. UTAG domains (single-cell resolution)
6. Compare methods and characterise domains

Step 1: Load Visium Data

For your own Space Ranger output set OUTS to the outs/ folder and skip the download. The demo is the 10x mouse brain (sagittal anterior) section.

pythonCopy code
import json
import tarfile
import urllib.request
from pathlib import Path

import numpy as np
import pandas as pd
import scanpy as sc

CACHE = Path("data_cache/V1_Mouse_Brain_Sagittal_Anterior")
BASE = ("https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Mouse_Brain_Sagittal_Anterior/"
        "V1_Mouse_Brain_Sagittal_Anterior_")
CACHE.mkdir(parents=True, exist_ok=True)
for name in ["filtered_feature_bc_matrix.h5", "spatial.tar.gz"]:
    if not (CACHE / name).exists():
        urllib.request.urlretrieve(BASE + name, CACHE / name)          # ~20 MB + ~9 MB
if not (CACHE / "spatial").exists():
    with tarfile.open(CACHE / "spatial.tar.gz") as tar:
        tar.extractall(CACHE, filter="data")
OUTS = CACHE                                    # your data: Path("path/to/outs")

SPATIAL_DIR = OUTS / "spatial"
adata = sc.read_10x_h5(OUTS / "filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()
cols = ["in_tissue", "array_row", "array_col", "pxl_row", "pxl_col"]
if (SPATIAL_DIR / "tissue_positions.csv").exists():                    # Space Ranger >= 2.0
    pos = pd.read_csv(SPATIAL_DIR / "tissue_positions.csv", index_col=0)
    pos.columns = cols
else:                                                                   # older, no header
    pos = pd.read_csv(SPATIAL_DIR / "tissue_positions_list.csv", header=None, index_col=0, names=cols)
adata.obs = adata.obs.join(pos, how="left")
adata.obsm["spatial"] = adata.obs[["pxl_col", "pxl_row"]].to_numpy(float)   # x, y in full-res pixels
scale = json.loads((SPATIAL_DIR / "scalefactors_json.json").read_text())
print(adata)

Step 2: Preprocess

pythonCopy code
adata.var["mt"] = adata.var_names.str.upper().str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)
sc.pp.filter_cells(adata, min_counts=500)       # Visium; use ~10 for imaging panels
sc.pp.filter_genes(adata, min_cells=10)
adata.layers["counts"] = adata.X.copy()         # SpaGCN and GraphST start from counts
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
print(f"{adata.n_obs} spots, {adata.n_vars} genes, median mt% {adata.obs['pct_counts_mt'].median():.1f}")

Step 3: SpaGCN (Visium with H&E)

SpaGCN smooths expression over a graph whose edges combine spot distance and histology colour, then clusters with a graph convolutional network. Coordinates passed with the image must be (image row, image column) at the image's scale.

pythonCopy code
import random

import SpaGCN as spg
import torch
from PIL import Image

N_DOMAINS = 7
USE_HISTOLOGY = True

sp = adata.copy()
sp.X = sp.layers["counts"].toarray()            # dense: SpaGCN calls .A on sparse, removed in scipy>=1.14
hires = scale["tissue_hires_scalef"]
x_pix = (sp.obs["pxl_row"] * hires).round().astype(int).to_numpy()   # image row
y_pix = (sp.obs["pxl_col"] * hires).round().astype(int).to_numpy()   # image column
if USE_HISTOLOGY:
    img = np.array(Image.open(SPATIAL_DIR / "tissue_hires_image.png").convert("RGB"))
    beta = max(3, int(scale["spot_diameter_fullres"] * hires))       # spot size in hires pixels
    adj = spg.calculate_adj_matrix(x=x_pix, y=y_pix, x_pixel=x_pix, y_pixel=y_pix,
                                   image=img, beta=beta, alpha=1, histology=True)
else:
    adj = spg.calculate_adj_matrix(x=x_pix, y=y_pix, histology=False)

spg.prefilter_genes(sp, min_cells=3)
spg.prefilter_specialgenes(sp)                  # drops ERCC / MT- genes
sc.pp.normalize_total(sp)
sc.pp.log1p(sp)

l = spg.search_l(0.5, adj, start=0.01, end=1000, tol=0.01, max_run=100)   # p=0.5
r_seed = t_seed = n_seed = 100
res = spg.search_res(sp, adj, l, N_DOMAINS, start=0.7, step=0.1, tol=5e-3, lr=0.05,
                     max_epochs=20, r_seed=r_seed, t_seed=t_seed, n_seed=n_seed)
random.seed(r_seed); torch.manual_seed(t_seed); np.random.seed(n_seed)
clf = spg.SpaGCN()
clf.set_l(l)
clf.train(sp, adj, init_spa=True, init="louvain", res=res, tol=5e-3, lr=0.05, max_epochs=200)
y_pred, prob = clf.predict()
adata.obs["spagcn"] = pd.Categorical(y_pred.astype(str))

# Refinement: relabel a spot when most of its ring disagrees (spot grids only)
adj_2d = spg.calculate_adj_matrix(x=sp.obs["array_row"].to_numpy(), y=sp.obs["array_col"].to_numpy(),
                                  histology=False)
refined = spg.refine(sample_id=sp.obs_names.tolist(), pred=y_pred.tolist(), dis=adj_2d, shape="hexagon")
adata.obs["spagcn_refined"] = pd.Categorical(np.asarray(refined).astype(str))
print(adata.obs["spagcn_refined"].value_counts())

Without louvain, skip search_res and initialise with k-means: clf.train(sp, adj, init_spa=True, init="kmeans", n_clusters=N_DOMAINS, tol=5e-3, lr=0.05, max_epochs=200).

shape="hexagon" for Visium, "square" for legacy ST. Raise alpha to give histology more weight; set USE_HISTOLOGY = False when the image is missing or poorly registered.

Step 4: GraphST

GraphST learns a spatially smoothed embedding by contrasting the real graph with a corrupted one, then clusters it. Give it raw counts: it selects 3,000 HVGs (seurat_v3), normalises, log-transforms and scales internally.

pythonCopy code
import torch
from GraphST import GraphST
from GraphST.utils import clustering

gst = sc.AnnData(adata.layers["counts"].copy(), obs=adata.obs[[]].copy(), var=adata.var[[]].copy())
gst.obsm["spatial"] = adata.obsm["spatial"].copy()
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
# datatype: "10X" (spot grids, imaging cells); "Slide" or "Stereo" for Slide-seq / Stereo-seq (KNN graph)
model = GraphST.GraphST(gst, device=device, random_seed=42, datatype="10X")
gst = model.train()                             # 600 epochs; embedding in obsm["emb"]

# Scans Leiden resolution from `end` down to `start` until exactly N_DOMAINS clusters appear
clustering(gst, n_clusters=N_DOMAINS, radius=50, method="leiden",
           start=0.1, end=2.0, increment=0.01, refinement=False)
adata.obs["graphst"] = pd.Categorical(gst.obs["domain"].astype(str).to_numpy())
adata.obsm["graphst_emb"] = gst.obsm["emb"]
print(adata.obs["graphst"].value_counts())

The GraphST authors recommend method="mclust" (needs R with the mclust package and rpy2); refinement=True relabels each observation by its radius nearest neighbours.

Step 5: UTAG (Single-Cell Resolution)

UTAG averages each cell's expression over neighbours within max_dist, then clusters. Input: one .h5ad of cells with log-normalised X and obsm["spatial"] in microns (e.g. the output of spatial-transcriptomics-cell-type-annotation).

pythonCopy code
from scipy.spatial import cKDTree
from sklearn.neighbors import NearestNeighbors
from utag import utag

cells = sc.read_h5ad("celltype_annotated.h5ad")
SLIDE_KEY = None                     # obs column with sample/section IDs, if several
MIN_NICHES, MIN_CLUSTER_SIZE = 5, 100

# max_dist: radius where a typical cell has ~3 neighbours (estimated per slide)
slides = cells.obs[SLIDE_KEY].unique() if SLIDE_KEY else [None]
d3 = []
for s in slides:
    xy = cells.obsm["spatial"] if s is None else cells[cells.obs[SLIDE_KEY] == s].obsm["spatial"]
    d3.append(cKDTree(xy).query(xy[:20000], k=4)[0][:, 3])
MAX_DIST = float(np.ceil(np.median(np.concatenate(d3))))
print(f"max_dist = {MAX_DIST}")

for resolution in [0.05, 0.1, 0.3, 0.5]:          # coarsest resolution giving enough niches
    out = utag(cells, slide_key=SLIDE_KEY, max_dist=MAX_DIST, normalization_mode="l1_norm",
               apply_umap=False, apply_clustering=True, clustering_method=["leiden"],
               resolutions=[resolution])
    key = f"UTAG Label_leiden_{resolution}"
    print(resolution, out.obs[key].nunique())
    if out.obs[key].nunique() >= MIN_NICHES:
        break

# Merge clusters smaller than MIN_CLUSTER_SIZE (within each slide) into the nearest larger one
lab = out.obs[key].astype(str).copy()
for s in slides:
    m = np.ones(out.n_obs, bool) if s is None else (out.obs[SLIDE_KEY] == s).to_numpy()
    counts = lab[m].value_counts()
    small = lab[m].isin(counts[counts < MIN_CLUSTER_SIZE].index).to_numpy()
    if small.any() and (~small).any():
        xy = out.obsm["spatial"][m]
        nn = NearestNeighbors(n_neighbors=1).fit(xy[~small])
        idx = nn.kneighbors(xy[small], return_distance=False).ravel()
        lab.iloc[np.flatnonzero(m)[small]] = lab[m].to_numpy()[~small][idx]
cells.obs["utag"] = pd.Categorical(lab.reindex(cells.obs_names))   # utag may reorder cells by slide
print(cells.obs["utag"].value_counts())

Step 6: Compare Methods and Characterise Domains

pythonCopy code
from sklearn.metrics import adjusted_rand_score

ari = adjusted_rand_score(adata.obs["spagcn_refined"], adata.obs["graphst"])
print(f"SpaGCN vs GraphST ARI: {ari:.2f}")      # > ~0.5: the main boundaries agree
print(pd.crosstab(adata.obs["spagcn_refined"], adata.obs["graphst"]))

import matplotlib.pyplot as plt

axes = sc.pl.embedding(adata, basis="spatial", color=["spagcn_refined", "graphst"], size=1.5, show=False)
for ax in axes:
    ax.invert_yaxis()                            # pixel y grows downwards
plt.savefig("domains_spatial.png", dpi=150, bbox_inches="tight")

DOMAIN = "spagcn_refined"
sc.tl.rank_genes_groups(adata, DOMAIN, method="wilcoxon")
markers = sc.get.rank_genes_groups_df(adata, group=None)
markers.to_csv("domain_markers.csv", index=False)
top = (markers[markers["logfoldchanges"] > 0].groupby("group").head(8)
       .groupby("group")["names"].apply(", ".join))
print(top)
adata.obs[["spagcn", "spagcn_refined", "graphst"]].to_csv("spatial_domains.csv")
adata.write_h5ad("spatial_domains.h5ad", compression="gzip")

Name domains from their markers (e.g. Mbp / Plp1 = white matter, layer markers in cortex). If cell type labels or cell2location abundances exist, pd.crosstab(domain, cell_type, normalize="index") shows what each domain is made of. With an expert annotation, report ARI against it instead.

Key Parameters

Common Recipes

Recipe: Choose the Number of Domains

N_DOMAINS is a guess unless the anatomy fixes it. Scan k on the GraphST embedding and look for a plateau or peak in silhouette:

pythonCopy code
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

emb = adata.obsm["graphst_emb"]
rng = np.random.default_rng(0)
sub = rng.choice(emb.shape[0], min(emb.shape[0], 5000), replace=False)
scores = {}
for k in range(4, 16):
    labels = KMeans(k, n_init=10, random_state=0).fit_predict(emb)
    scores[k] = silhouette_score(emb[sub], labels[sub])
print(pd.Series(scores).round(3))

Prefer a k that matches known anatomy over the highest score; report the scan.

Recipe: What Each Domain Is Made Of

With cell type labels (single cells) or cell2location abundances (spots):

pythonCopy code
DOMAIN = "spagcn_refined"
if "cell_type" in adata.obs:
    comp = pd.crosstab(adata.obs[DOMAIN], adata.obs["cell_type"], normalize="index")
else:
    abund = pd.read_csv("cell2location_abundance.csv", index_col=0).reindex(adata.obs_names)
    prop = abund.div(abund.sum(axis=1), axis=0)
    comp = prop.groupby(adata.obs[DOMAIN], observed=True).mean()
comp.round(3).to_csv("domain_composition.csv")
print(comp.apply(lambda r: ", ".join(f"{k} {v:.0%}" for k, v in r.nlargest(3).items()), axis=1))

Recipe: Domain Boundaries and Interfaces

Which domains touch, and which spots sit on a boundary (interfaces such as the tumour edge are often where the biology is):

pythonCopy code
import squidpy as sq

sq.gr.spatial_neighbors(adata, coord_type="grid", n_rings=1)   # imaging cells: coord_type="generic", n_neighs=10
A = adata.obsp["spatial_connectivities"].tocoo()
lab = adata.obs[DOMAIN].astype(str).to_numpy()
edges = pd.DataFrame({"a": lab[A.row], "b": lab[A.col]})
touch = pd.crosstab(edges["a"], edges["b"])
print(touch)                                       # off-diagonal counts = shared boundary length

diff = lab[A.row] != lab[A.col]
boundary = np.zeros(adata.n_obs, bool)
boundary[A.row[diff]] = True
adata.obs["domain_boundary"] = boundary
print(f"boundary spots: {boundary.mean():.0%}")

Expected Outputs

• spatial_domains.csv - spagcn, spagcn_refined, graphst labels per spot
• spatial_domains.h5ad - data with domain labels and obsm["graphst_emb"]
• domain_markers.csv - Wilcoxon markers per domain
• domains_spatial.png - domain maps on tissue coordinates
• domain_composition.csv - cell type mix per domain (recipe)

Troubleshooting

Related Skills

• spatial-transcriptomics-cell-type-annotation - cell labels; its Step 7 builds composition niches
• spatial-transcriptomics-deconvolution-cell2location - cell type abundance per spot to interpret domains
• spatial-transcriptomics-analysis-squidpy - neighbourhood enrichment, co-occurrence, spatially variable genes
• spatial-transcriptomics-cell-cell-communication - signalling between cell types inside or across domains
• trajectory-inference - pseudotime within a domain

References

• Hu et al. (2021) SpaGCN, Nature Methods: https://doi.org/10.1038/s41592-021-01255-8
• Long et al. (2023) GraphST, Nature Communications: https://doi.org/10.1038/s41467-023-36796-3
• Kim et al. (2022) UTAG, Nature Methods: https://doi.org/10.1038/s41592-022-01657-2
• SpaGCN source: https://github.com/jianhuupenn/SpaGCN
• GraphST source: https://github.com/JinmiaoChenLab/GraphST
• UTAG source: https://github.com/ElementoLab/utag