{
  "skill": "alterlab-squidpy-spatial",
  "evals": [
    {
      "id": "xenium-neighborhood-enrichment",
      "prompt": "I have a 10x Xenium run with cells already clustered into cell types. I want to know which cell types tend to sit next to each other in the tissue — is there spatial co-localization between my tumor and immune clusters?",
      "expected_output": "Invokes alterlab-squidpy-spatial. Recognizes Xenium as single-cell resolution, builds the graph with sq.gr.spatial_neighbors(coord_type='generic', delaunay=True), then runs sq.gr.nhood_enrichment on the cluster_key and reports the permutation z-score matrix (positive z = spatial attraction/co-localization, negative = avoidance). Does NOT use coord_type='grid' (that is for Visium spots).",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "nhood_enrichment" },
        { "type": "behavior", "value": "Chooses coord_type='generic' (Delaunay) for single-cell Xenium data and reports neighborhood-enrichment z-scores with their permutation context, not raw adjacency counts." }
      ]
    },
    {
      "id": "visium-spatial-variable-genes",
      "prompt": "From my Visium slide I want the list of genes that vary across tissue space — spatially variable genes ranked by Moran's I. Walk me through it in squidpy.",
      "expected_output": "Invokes alterlab-squidpy-spatial. Builds a Visium spot graph with sq.gr.spatial_neighbors(coord_type='grid', n_neighs=6), runs sq.gr.spatial_autocorr(mode='moran'), and returns genes ranked by Moran's I from adata.uns['moranI']. Notes Visium spots are multi-cell and that cell-type-level claims require deconvolution via alterlab-scvi-tools.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "spatial_autocorr" },
        { "type": "behavior", "value": "Uses coord_type='grid' for Visium and ranks spatially variable genes by Moran's I via spatial_autocorr, flagging that Visium spots are multi-cell." }
      ]
    },
    {
      "id": "merfish-cooccurrence-distance",
      "prompt": "Using squidpy on my MERFISH data, how do I measure whether two cell types co-occur as a function of distance, and also build the neighbor graph correctly for single-molecule imaging data?",
      "expected_output": "Invokes alterlab-squidpy-spatial. Treats MERFISH as single-cell (coord_type='generic', delaunay=True or n_neighs), runs sq.gr.co_occurrence over the cluster_key to get the co-occurrence probability ratio across radial distance, and visualizes with sq.pl.co_occurrence.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "co_occurrence" }
      ]
    },
    {
      "id": "spatial-ligand-receptor",
      "prompt": "I want to run a ligand-receptor interaction analysis between my annotated clusters in spatial transcriptomics data using squidpy. Which function and how do I read the output?",
      "expected_output": "Invokes alterlab-squidpy-spatial. Uses sq.gr.ligrec with the cluster_key (Omnipath-backed permutation test of L-R co-expression between cluster pairs) and sq.pl.ligrec to visualize, explaining the means/p-values output per ligand-receptor pair per cluster pair.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "ligrec" }
      ]
    },
    {
      "id": "near-miss-alterlab-scanpy",
      "prompt": "I just loaded a Xenium dataset into AnnData. Run the standard QC, normalize, find highly variable genes, PCA, build the kNN graph, and Leiden-cluster the cells so I have cell types.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-scanpy. The request is the standard non-spatial single-cell pipeline (QC, normalization, HVGs, PCA, sc.pp.neighbors, Leiden) that produces the cluster labels; squidpy's spatial graph and neighborhood statistics come afterward. alterlab-squidpy-spatial assumes labels already exist.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-scanpy" }
      ]
    },
    {
      "id": "near-miss-alterlab-scvi-tools",
      "prompt": "My Visium spots are mixtures of cells. I have a matching scRNA-seq reference with cell-type labels and want to deconvolve each spot into cell-type proportions using destVI.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-scvi-tools. Spot deconvolution / mapping cell-type proportions onto Visium spots (destVI, Tangram) is a deep generative modeling task covered by alterlab-scvi-tools (references/models-spatial.md). alterlab-squidpy-spatial does neighborhood/graph statistics, not deconvolution.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-scvi-tools" }
      ]
    },
    {
      "id": "near-miss-alterlab-anndata",
      "prompt": "I have several .h5ad files from different Visium slides. Help me concatenate them into one AnnData, align the var index, and move the raw counts into a layer — no analysis yet, just the data structure.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-anndata. The request is pure AnnData object construction — concatenation, var-index alignment, layer management — with no spatial analysis. alterlab-anndata is the data-format skill; alterlab-squidpy-spatial only runs once the object is ready and labels exist.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-anndata" }
      ]
    }
  ]
}

# Squidpy Analysis Recipes

Copy-paste recipes for the core `sq.gr` graph statistics and their `sq.pl`
visualizations, with the parameters that matter and how to read each output. These
assume an AnnData (or SpatialData `table`) with cluster/cell-type labels in
`adata.obs[cluster_key]` and a spatial graph already built (see
`platform_routing.md`). Function names are verified against the squidpy 1.8 API.

## 1. Spatial neighbor graph (prerequisite)

```python
import squidpy as sq

# Single-cell (Xenium / MERFISH / CosMx)
sq.gr.spatial_neighbors(adata, coord_type="generic", delaunay=True)

# Visium spots
sq.gr.spatial_neighbors(adata, coord_type="grid", n_neighs=6, n_rings=1)
```

Writes `adata.obsp["spatial_connectivities"]` and `adata.obsp["spatial_distances"]`,
plus `adata.uns["spatial_neighbors"]`.

## 2. Neighborhood enrichment

```python
sq.gr.nhood_enrichment(adata, cluster_key="leiden")
sq.pl.nhood_enrichment(adata, cluster_key="leiden")
z = adata.uns["leiden_nhood_enrichment"]["zscore"]   # cluster x cluster z-scores
```

Permutation z-score per cluster pair: positive = adjacent more than chance
(attraction), negative = avoidance. Report with the permutation context.

## 3. Co-occurrence across distance

```python
sq.gr.co_occurrence(adata, cluster_key="leiden")
sq.pl.co_occurrence(adata, cluster_key="leiden", clusters="0")
```

Computes the co-occurrence probability ratio of cluster pairs as a function of
radial distance. Use it when you care about *how adjacency decays with distance*,
not just binary neighbor counts. Stored in `adata.uns["leiden_co_occurrence"]`.

## 4. Spatially variable genes — Moran's I

```python
sq.gr.spatial_autocorr(adata, mode="moran")
top = adata.uns["moranI"].sort_values("I", ascending=False).head(20)
```

`mode="moran"` writes the `adata.uns["moranI"]` table (Moran's I, p-value, FDR per
gene); `mode="geary"` writes `gearyC`. High positive Moran's I = strong spatial
autocorrelation = a spatially variable gene. Restrict to highly variable genes
first to cut runtime.

### Alternative: Sepal

```python
sq.gr.sepal(adata, max_neighs=6)   # max_neighs=6 hex (Visium) / 4 square grid
```

`sq.gr.sepal` identifies spatially variable genes via a diffusion-based score; it is
an alternative to Moran's I, useful as a cross-check.

## 5. Ligand-receptor analysis

```python
sq.gr.ligrec(adata, cluster_key="leiden")
sq.pl.ligrec(adata, cluster_key="leiden")
```

Tests ligand-receptor co-expression between cluster pairs using an Omnipath-backed
interaction database (a permutation test in the spirit of CellPhoneDB). Returns
means and p-values per L-R pair per cluster pair.

## 6. Other graph statistics

```python
sq.gr.interaction_matrix(adata, cluster_key="leiden")   # raw inter-cluster edge counts
sq.gr.centrality_scores(adata, cluster_key="leiden")    # degree / closeness / clustering per cluster
sq.gr.ripley(adata, cluster_key="leiden", mode="L")     # Ripley's L vs. complete spatial randomness
```

`ripley` (modes `"F"`, `"G"`, `"L"`) tests whether a cluster is more clustered or
dispersed than complete spatial randomness. `centrality_scores` summarizes each
cluster's position in the spatial graph.

## 7. Visualizing tissue

```python
sq.pl.spatial_scatter(adata, color="leiden")    # spot/point overlay (Visium, points)
sq.pl.spatial_segment(adata, color="leiden", seg_cell_id="cell_id")  # segmented cells (Xenium)
```

## 8. Image features (optional, H&E / IF)

```python
from squidpy.im import ImageContainer
img = ImageContainer("tissue_image.tif")
sq.im.process(img, layer="image", method="smooth")
sq.im.segment(img, layer="image", method="watershed")
sq.im.calculate_image_features(adata, img, features="summary")
```

The `sq.im` module operates on an `ImageContainer` and writes per-observation image
features back into `adata.obs` for joint expression+morphology analysis. Skip this
entirely for pure transcriptomic neighborhood work.

# Platform → coord_type Routing

The single most consequential choice in a squidpy analysis is how the spatial
neighbor graph is built, because every downstream statistic (`nhood_enrichment`,
`co_occurrence`, `spatial_autocorr`, `ripley`, `interaction_matrix`,
`centrality_scores`) consumes that graph. The right builder depends on the
measurement geometry of the platform. All parameter semantics below are from the
squidpy 1.8 `sq.gr.spatial_neighbors` API.

## Decision table

| Platform | Spatial unit | Geometry | `coord_type` | Key args | Notes |
|----------|--------------|----------|--------------|----------|-------|
| Visium | spot (10–100s of cells) | hexagonal grid | `"grid"` | `n_neighs=6`, `n_rings=1` (or `2`) | 6 immediate neighbors per spot; raise `n_rings` for wider neighborhoods |
| Visium HD | 2/8/16 µm bin | square grid | `"grid"` | `n_neighs=4` or `6` | square lattice; pick the bin size before graphing |
| Xenium | single cell | irregular points | `"generic"` | `delaunay=True` | Delaunay triangulation gives adjacency without a k cutoff |
| MERFISH / MERSCOPE | single cell | irregular points | `"generic"` | `delaunay=True` or `n_neighs=k` | k-NN (`delaunay=False`) when you want a fixed degree |
| CosMx (Nanostring) | single cell | irregular points | `"generic"` | `delaunay=True` | same as Xenium/MERFISH |

## Parameter semantics (squidpy 1.8 `sq.gr.spatial_neighbors`)

- **`coord_type`** — `"grid"`, `"generic"`, or `None`.
  - `None` auto-selects `"grid"` **only** when `spatial` is present in `adata.uns`
    with `n_neighs == 6` (the Visium signature); otherwise it uses `"generic"`.
  - Do not rely on auto-detection — pass `coord_type` explicitly so the geometry is
    unambiguous and reproducible.
- **`n_neighs`** (default `6`) — for `"grid"`, the number of neighboring tiles; for
  `"generic"`, the number of nearest neighbors, applied only when `delaunay=False`.
- **`n_rings`** (default `1`) — number of rings of neighbors; **only used for
  `coord_type="grid"`**. `n_rings=2` includes second-shell spots.
- **`delaunay`** (default `False`) — build the graph from a Delaunay triangulation;
  **only used for `coord_type="generic"`**. Preferred for single-cell platforms
  because it adapts to local density without a fixed `k`.

## Per-platform gotchas

- **Visium spots are multi-cell.** Any cell-type-level claim from raw Visium needs
  **deconvolution first** (destVI / Tangram via `alterlab-scvi-tools`,
  `references/models-spatial.md`). Squidpy graph statistics on undeconvolved spots
  describe *spot-cluster* adjacency, not single-cell adjacency.
- **Single-cell platforms (Xenium/MERFISH/CosMx)** carry true per-cell coordinates,
  so `nhood_enrichment` and `co_occurrence` are directly interpretable at cell-type
  resolution — no deconvolution step.
- **Visium HD** must have its bin size chosen before analysis; 2 µm bins are near
  single-cell but sparse, 8/16 µm bins trade resolution for counts.
- **Clusters/labels must already exist.** `nhood_enrichment`, `co_occurrence`,
  `interaction_matrix`, `centrality_scores`, and `ligrec` all take a `cluster_key`
  pointing at a categorical `adata.obs` column produced upstream (Leiden via
  `alterlab-scanpy`, or imported cell-type annotations).

## Reading `nhood_enrichment` output

`sq.gr.nhood_enrichment` runs a permutation test and stores a **z-score** matrix in
`adata.uns[f"{cluster_key}_nhood_enrichment"]["zscore"]`. Positive z = the two
clusters are adjacent more often than expected under the permuted null (attraction);
negative z = avoidance. Always report z-scores with that permutation context, never
raw adjacency counts.

# Reading Modern Platforms with spatialdata-io

Squidpy 1.8 requires `spatialdata>=0.7.2` and accepts `SpatialData` objects, which
is the recommended representation for single-cell-resolution platforms (Xenium,
Visium HD) that carry images, shapes, and points alongside the expression table.
Squidpy's own `sq.read` module covers `visium`, `vizgen`, and `nanostring` only —
**there is no `sq.read.xenium`** — so Xenium and Visium HD come in through
`spatialdata-io`. Reader names below are verified against the spatialdata-io stable
API.

## Readers

```python
from spatialdata_io import xenium, visium, visium_hd, merscope

sdata = xenium("path/to/xenium_outs/")        # 10x Genomics Xenium
sdata = visium("path/to/visium_outs/")        # 10x Genomics Visium (spatialdata form)
sdata = visium_hd("path/to/visium_hd_outs/")  # 10x Genomics Visium HD
sdata = merscope("path/to/merscope/")         # Vizgen MERSCOPE
```

Each returns a `SpatialData` object — a container of elements: `images`, `labels`
(segmentation masks), `shapes` (cell/nucleus boundaries), `points` (transcript
locations), and `tables` (the AnnData expression table(s)).

## Getting the AnnData table squidpy operates on

Squidpy graph functions run on the expression `table` (an AnnData) inside the
`SpatialData` object. Access it via the `tables` mapping:

```python
adata = sdata.tables["table"]            # the AnnData squidpy analyses use
# ... run scanpy QC/clustering (alterlab-scanpy) to populate adata.obs["leiden"] ...
import squidpy as sq
sq.gr.spatial_neighbors(adata, coord_type="generic", delaunay=True)  # Xenium = single cell
sq.gr.nhood_enrichment(adata, cluster_key="leiden")
```

The exact key under `tables` depends on the reader/dataset; inspect `sdata.tables`
(a dict-like mapping) to find it. Spatial coordinates live in `adata.obsm["spatial"]`,
which `spatial_neighbors` reads by default (`spatial_key="spatial"`).

## When to stay in AnnData vs. SpatialData

- **Visium (legacy spot data):** `sq.read.visium(...)` returns a plain AnnData that
  already has `adata.uns["spatial"]` and `adata.obsm["spatial"]`; you can skip
  SpatialData entirely. Use `coord_type="grid"`.
- **Xenium / Visium HD / large MERSCOPE:** prefer `spatialdata-io` + `SpatialData`
  so images, segmentation shapes, and transcript points stay aligned; extract the
  `table` for squidpy graph statistics. Use `coord_type="generic"` (Xenium/MERSCOPE)
  or `"grid"` (Visium HD bins) per `platform_routing.md`.

## Hand-offs

- Non-spatial QC, normalization, PCA/UMAP, Leiden clustering, marker genes on the
  extracted `table` → `alterlab-scanpy`.
- Spot deconvolution for Visium (destVI / Tangram) → `alterlab-scvi-tools`
  (`references/models-spatial.md`).
- Pure `.h5ad` AnnData construction/slicing with no spatial analysis →
  `alterlab-anndata`.

#!/usr/bin/env python3
"""
Spatial neighborhood analysis for squidpy (1.8.x).

Given a clustered AnnData (.h5ad) that already carries a categorical cluster /
cell-type column in ``adata.obs`` (produce it upstream with alterlab-scanpy), this
script:

  1. builds the spatial neighbor graph with the *correct* ``coord_type`` for the
     platform (the one decision that invalidates everything downstream if wrong),
  2. runs ``sq.gr.nhood_enrichment`` (permutation z-scores per cluster pair),
  3. runs ``sq.gr.co_occurrence`` (adjacency vs. radial distance),
  4. runs ``sq.gr.spatial_autocorr`` (Moran's I -> spatially variable genes),

then writes a compact JSON summary and (optionally) the updated .h5ad.

Platform -> coord_type (squidpy 1.8 sq.gr.spatial_neighbors semantics):
  visium / visium_hd        -> coord_type="grid"     (n_neighs / n_rings)
  xenium / merfish / cosmx  -> coord_type="generic"  (delaunay=True)

This runs locally and sends no data anywhere. squidpy + scanpy must be installed
(``uv run --with squidpy python ...`` or a project env). For large single-cell
panels this is a good candidate to run on local compute rather than burning API
calls.

Usage:
    uv run python spatial_neighborhood.py clustered.h5ad \
        --platform xenium --cluster-key leiden --out spatial_report.json
"""
from __future__ import annotations

import argparse
import json
import sys

# Platform -> (coord_type, kwargs for sq.gr.spatial_neighbors)
GRID_PLATFORMS = {"visium", "visium_hd"}
GENERIC_PLATFORMS = {"xenium", "merfish", "merscope", "cosmx"}


def graph_kwargs(platform: str, n_neighs: int, n_rings: int):
    """Return the spatial_neighbors kwargs for the platform."""
    platform = platform.lower()
    if platform in GRID_PLATFORMS:
        return {"coord_type": "grid", "n_neighs": n_neighs, "n_rings": n_rings}
    if platform in GENERIC_PLATFORMS:
        # Delaunay adapts to local density without a fixed k -> preferred single-cell.
        return {"coord_type": "generic", "delaunay": True}
    raise ValueError(
        f"unknown platform {platform!r}; "
        f"expected one of {sorted(GRID_PLATFORMS | GENERIC_PLATFORMS)}"
    )


def main(argv: list[str] | None = None) -> int:
    p = argparse.ArgumentParser(
        description="Build a squidpy spatial graph and run core neighborhood statistics."
    )
    p.add_argument("h5ad", help="Path to a clustered AnnData .h5ad (labels in obs).")
    p.add_argument(
        "--platform",
        required=True,
        help="visium | visium_hd | xenium | merfish | merscope | cosmx",
    )
    p.add_argument(
        "--cluster-key",
        default="leiden",
        help="Categorical obs column with cluster / cell-type labels (default: leiden).",
    )
    p.add_argument("--n-neighs", type=int, default=6, help="grid n_neighs (default: 6).")
    p.add_argument("--n-rings", type=int, default=1, help="grid n_rings (default: 1).")
    p.add_argument(
        "--moran-top", type=int, default=20, help="How many top SVGs to report (default: 20)."
    )
    p.add_argument("--out", default=None, help="Write JSON summary here (default: stdout).")
    p.add_argument(
        "--write-h5ad", default=None, help="Optionally write the updated AnnData here."
    )
    args = p.parse_args(argv)

    try:
        import scanpy as sc
        import squidpy as sq
    except ImportError as exc:  # pragma: no cover - dependency guard
        print(
            f"ERROR: this script needs scanpy + squidpy installed ({exc}). "
            "Try: uv run --with squidpy python " + __file__,
            file=sys.stderr,
        )
        return 2

    try:
        kwargs = graph_kwargs(args.platform, args.n_neighs, args.n_rings)
    except ValueError as exc:
        print(f"ERROR: {exc}", file=sys.stderr)
        return 2

    adata = sc.read_h5ad(args.h5ad)

    if args.cluster_key not in adata.obs:
        print(
            f"ERROR: cluster key {args.cluster_key!r} not found in adata.obs "
            f"(have: {list(adata.obs.columns)}). Cluster upstream with alterlab-scanpy.",
            file=sys.stderr,
        )
        return 2

    # 1. spatial graph (the consequential choice)
    sq.gr.spatial_neighbors(adata, **kwargs)

    # 2. neighborhood enrichment (permutation z-scores)
    sq.gr.nhood_enrichment(adata, cluster_key=args.cluster_key)
    nhood = adata.uns.get(f"{args.cluster_key}_nhood_enrichment", {})
    zscore = nhood.get("zscore")
    categories = list(adata.obs[args.cluster_key].cat.categories) if hasattr(
        adata.obs[args.cluster_key], "cat"
    ) else sorted(set(map(str, adata.obs[args.cluster_key])))

    # 3. co-occurrence across distance
    sq.gr.co_occurrence(adata, cluster_key=args.cluster_key)

    # 4. spatially variable genes via Moran's I
    sq.gr.spatial_autocorr(adata, mode="moran")
    moran = adata.uns.get("moranI")
    top_svgs = []
    if moran is not None:
        moran_sorted = moran.sort_values("I", ascending=False).head(args.moran_top)
        for gene, row in moran_sorted.iterrows():
            entry = {"gene": str(gene), "morans_I": float(row["I"])}
            if "pval_norm" in moran_sorted.columns:
                entry["pval_norm"] = float(row["pval_norm"])
            if "pval_norm_fdr_bh" in moran_sorted.columns:
                entry["fdr_bh"] = float(row["pval_norm_fdr_bh"])
            top_svgs.append(entry)

    report = {
        "tool": "alterlab-squidpy-spatial/spatial_neighborhood.py",
        "version": "1.0.0",
        "input": args.h5ad,
        "platform": args.platform.lower(),
        "spatial_neighbors_kwargs": kwargs,
        "cluster_key": args.cluster_key,
        "n_obs": int(adata.n_obs),
        "n_vars": int(adata.n_vars),
        "clusters": categories,
        "nhood_enrichment_zscore": (
            zscore.tolist() if zscore is not None else None
        ),
        "top_spatially_variable_genes": top_svgs,
        "notes": (
            "z-scores are permutation-test based: positive = spatial attraction, "
            "negative = avoidance. For Visium spot data, deconvolve "
            "(alterlab-scvi-tools) before cell-type-level claims."
        ),
    }

    if args.write_h5ad:
        adata.write(args.write_h5ad)
        report["written_h5ad"] = args.write_h5ad

    payload = json.dumps(report, indent=2)
    if args.out:
        with open(args.out, "w") as fh:
            fh.write(payload)
        print(f"Wrote {args.out}", file=sys.stderr)
    else:
        print(payload)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())