Spatial

Imaging-based and spot-based assays: Xenium, Visium, CosMx and MERSCOPE. The container mirrors Seurat v5's — an FOV holding Centroids, Segmentation and Molecules boundaries, or a VisiumV2 holding the H&E image and its ScaleFactors.

radius on a Visium FOV

Seurat stores spot_diameter_fullres in the FOV's radius slot — a diameter where a radius is named — and Radius() on its own VisiumV2 returns NULL. load_visium stores a radius. The slide's fixed 100 µm spot pitch is what settles which reading is right; the working is in the Visium vignette.

find_spatially_variable_features computes Moran's I on R's inverse-square distance weights, not on a kNN graph. Those give different answers, and the kNN version was the bug. Like FindSpatiallyVariableFeatures on an object, it reads scale.data unless given layer="data".

Loading

What each loader has been checked against

Each loader builds the object Seurat's reader builds from the same files: the same cells under the same names, the same features, and one image under Seurat's name for it.

Loader	Seurat	Checked on
load_visium	Load10X_Spatial	10x's Visium mouse-brain slide (vignette)
load_xenium	LoadXenium	10x's Xenium mouse brain, Onboard Analysis 1.0.1 (vignette); synthetic bundles in the 1.x and 2.0+ layouts, cells.parquet included
load_cosmx	LoadNanostring	a synthetic bundle in CosMx's layout only
load_merscope	LoadVizgen	a synthetic bundle in MERSCOPE's layout only

The synthetic comparisons are tests/test_loaders_vs_seurat.py. They pin how each layout is read, but cannot show that a platform's current software still writes it, so load_cosmx, load_merscope and Xenium's 2.0+ layout have not been compared with R on a real run. The imaging loaders read centroids only: Seurat's also load cell polygons and transcripts into the image.

load_xenium

load_xenium(path: Union[str, Path], assay: str = 'Xenium', fov: str = 'fov', fov_column: Optional[str] = None, project: str = 'Xenium', keep_controls: bool = False)

Load a 10x Xenium output bundle into a Truecell object with images.

Expects (from the Xenium output folder): * cell_feature_matrix/ — 10x MTX triplet (barcodes/features/matrix) * cells.parquet or cells.csv[.gz] — with cell_id, x_centroid, y_centroid (and optionally fov / transcript QC)

Builds what Seurat's LoadXenium(molecule.coordinates = FALSE) builds from the same files:

• the assay holds the Gene Expression features. Seurat routes the control and codeword features to separate assays; load_xenium drops them, and keep_controls=True keeps every row in the one assay instead;
• cells.parquet is read when pandas has a parquet engine, as ReadXenium reads it when arrow is installed, and cells.csv.gz otherwise. A cell_id stored as bytes is decoded, as ReadXenium decodes it;
• one image, named fov (LoadXenium's argument, with its default), holds the centroids;
• segmentation_method comes from the cells table, or is "cell" when the table has none, as LoadXenium fills it.

truecell also carries the cells table's other columns into meta_data. fov_column, if given and present in the cells table, splits the cells into one image per value instead, named by the value.

Checked against LoadXenium on 10x's Xenium Onboard Analysis 1.0.1 mouse brain, and on synthetic bundles in the XOA 1.x and 2.0+ layouts. LoadXenium can also read cell_feature_matrix.h5; load_xenium needs the MTX directory.

Source code in truecell/spatial/loaders.py

def load_xenium(
    path: Union[str, Path],
    assay: str = "Xenium",
    fov: str = "fov",
    fov_column: Optional[str] = None,
    project: str = "Xenium",
    keep_controls: bool = False,
):
    """Load a 10x Xenium output bundle into a Truecell object with images.

    Expects (from the Xenium output folder):
      * ``cell_feature_matrix/`` — 10x MTX triplet (barcodes/features/matrix)
      * ``cells.parquet`` or ``cells.csv[.gz]`` — with ``cell_id``,
        ``x_centroid``, ``y_centroid`` (and optionally ``fov`` / transcript QC)

    Builds what Seurat's ``LoadXenium(molecule.coordinates = FALSE)`` builds from
    the same files:

    * the assay holds the ``Gene Expression`` features. Seurat routes the control
      and codeword features to separate assays; ``load_xenium`` drops them, and
      ``keep_controls=True`` keeps every row in the one assay instead;
    * ``cells.parquet`` is read when pandas has a parquet engine, as
      ``ReadXenium`` reads it when ``arrow`` is installed, and ``cells.csv.gz``
      otherwise. A ``cell_id`` stored as bytes is decoded, as ``ReadXenium``
      decodes it;
    * one image, named ``fov`` (``LoadXenium``'s argument, with its default),
      holds the centroids;
    * ``segmentation_method`` comes from the cells table, or is ``"cell"`` when
      the table has none, as ``LoadXenium`` fills it.

    truecell also carries the cells table's other columns into ``meta_data``.
    ``fov_column``, if given and present in the cells table, splits the cells into
    one image per value instead, named by the value.

    Checked against ``LoadXenium`` on 10x's Xenium Onboard Analysis 1.0.1 mouse
    brain, and on synthetic bundles in the XOA 1.x and 2.0+ layouts.
    ``LoadXenium`` can also read ``cell_feature_matrix.h5``; ``load_xenium``
    needs the MTX directory.
    """
    path = Path(path)
    mtx_dir = path / "cell_feature_matrix"
    if not mtx_dir.exists():
        raise FileNotFoundError(
            f"{mtx_dir} not found. load_xenium expects the unpacked "
            "cell_feature_matrix/ MTX directory."
        )
    counts, feats, cells = read_10x(mtx_dir)

    if not keep_controls:
        types = _feature_types(mtx_dir)
        if types is not None and len(types) == len(feats) and "Gene Expression" in types:
            mask = np.array([t == "Gene Expression" for t in types])
            counts = counts[mask, :]
            feats = [f for f, m in zip(feats, mask) if m]

    cell_candidates = [path / n for n in ("cells.parquet", "cells.csv.gz", "cells.csv")
                       if (path / n).exists()]
    if not cell_candidates:
        raise FileNotFoundError(f"No cells.parquet/csv found in {path}.")
    cdf = None
    for cf in cell_candidates:                      # prefer parquet, fall back to csv
        try:
            cdf = _read_table(cf)
            break
        except ImportError:                         # no parquet engine → try next
            continue
    if cdf is None:
        raise ImportError(
            "cells.parquet found but no parquet engine is installed. Install "
            "pyarrow, or provide cells.csv[.gz] alongside it."
        )
    rename = {"cell_id": "cell", "x_centroid": "x", "y_centroid": "y"}
    cdf = cdf.rename(columns={k: v for k, v in rename.items() if k in cdf.columns})
    if not {"cell", "x", "y"} <= set(cdf.columns):
        raise ValueError("cells table must contain cell_id, x_centroid, y_centroid.")
    # A cell_id column stored as bytes comes back from pandas as `bytes`, which
    # str() would turn into "b'aaabinlp-1'" and match no barcode.
    cdf["cell"] = [c.decode() if isinstance(c, bytes) else str(c) for c in cdf["cell"]]
    if "segmentation_method" not in cdf.columns:
        cdf["segmentation_method"] = "cell"

    split = fov_column if (fov_column and fov_column in cdf.columns) else None
    coords = cdf[["cell", "x", "y"] + ([split] if split else [])]
    meta = cdf.set_index("cell").drop(columns=["x", "y"], errors="ignore")
    return _build_spatial_object(counts, feats, [str(c) for c in cells], coords,
                                 assay, project, fov=split, meta_data=meta,
                                 image_name=fov, radius=_auto_radius(coords))

load_visium

load_visium(path: Union[str, Path], assay: str = 'Spatial', project: str = 'Visium', image: bool = True, image_resolution: str = 'lowres', filter_by_tissue: bool = True, slice_name: str = 'slice1')

Load a 10x Visium output into a Truecell object with spot coordinates.

Expects: * filtered_feature_bc_matrix/ — 10x MTX triplet * spatial/tissue_positions.csv (or tissue_positions_list.csv) with barcode, in_tissue, array row/col and pixel row/col columns

Optionally (image=True, the default) also reads spatial/tissue_{hires,lowres}_image.png and spatial/scalefactors_json.json, producing a VisiumV2 image that carries the H&E tissue photo — what spatial_dim_plot / spatial_feature_plot draw on. A bundle with no PNG still loads; you just get a plain FOV, as before.

Parameters:

• image (bool, default: True ) –

  read the tissue image + scale factors (default True).

• image_resolution (str, default: 'lowres' ) –

  'lowres' (default) or 'hires'; falls back to whichever is present. Matches Read10X_Image's image.name = "tissue_lowres_image.png".

• filter_by_tissue (bool, default: True ) –

  keep only spots with in_tissue == 1 (default True), matching Read10X_Image's filter.matrix = TRUE.

• slice_name (str, default: 'slice1' ) –

  key for the FOV in obj.images (default "slice1", the name Load10X_Spatial uses).

Notes

Spot coordinates stay in full-resolution pixels, matching tissue_positions.csv. The scale factors convert them to image pixels — see VisiumV2.scale_coordinates.

The three defaults above changed to match Seurat. Reading a bundle the way earlier versions did is load_visium(path, image_resolution="hires", filter_by_tissue=False, slice_name="spatial").

Source code in truecell/spatial/loaders.py

def load_visium(
    path: Union[str, Path],
    assay: str = "Spatial",
    project: str = "Visium",
    image: bool = True,
    image_resolution: str = "lowres",
    filter_by_tissue: bool = True,
    slice_name: str = "slice1",
):
    """Load a 10x Visium output into a Truecell object with spot coordinates.

    Expects:
      * ``filtered_feature_bc_matrix/`` — 10x MTX triplet
      * ``spatial/tissue_positions.csv`` (or ``tissue_positions_list.csv``) with
        barcode, in_tissue, array row/col and pixel row/col columns

    Optionally (``image=True``, the default) also reads
    ``spatial/tissue_{hires,lowres}_image.png`` and ``spatial/scalefactors_json.json``,
    producing a :class:`~truecell.spatial.visium.VisiumV2` image that carries the H&E
    tissue photo — what ``spatial_dim_plot`` / ``spatial_feature_plot`` draw on. A
    bundle with no PNG still loads; you just get a plain FOV, as before.

    Parameters
    ----------
    image            : read the tissue image + scale factors (default True).
    image_resolution : 'lowres' (default) or 'hires'; falls back to whichever is
                       present. Matches ``Read10X_Image``'s
                       ``image.name = "tissue_lowres_image.png"``.
    filter_by_tissue : keep only spots with ``in_tissue == 1`` (default True),
                       matching ``Read10X_Image``'s ``filter.matrix = TRUE``.
    slice_name       : key for the FOV in ``obj.images`` (default ``"slice1"``,
                       the name ``Load10X_Spatial`` uses).

    Notes
    -----
    Spot coordinates stay in **full-resolution pixels**, matching
    ``tissue_positions.csv``. The scale factors convert them to image pixels —
    see :meth:`VisiumV2.scale_coordinates`.

    The three defaults above changed to match Seurat. Reading a bundle the way
    earlier versions did is ``load_visium(path, image_resolution="hires",
    filter_by_tissue=False, slice_name="spatial")``.
    """
    from .visium import VisiumV2, read_scale_factors, read_tissue_image

    path = Path(path)
    mtx_dir = _first_existing(path, ["filtered_feature_bc_matrix", "raw_feature_bc_matrix"])
    if mtx_dir is None:
        raise FileNotFoundError(f"No filtered_feature_bc_matrix/ in {path}.")
    counts, feats, cells = read_10x(mtx_dir)

    pos_file = _first_existing(path / "spatial",
                               ["tissue_positions.csv", "tissue_positions_list.csv"])
    if pos_file is None:
        raise FileNotFoundError(f"No spatial/tissue_positions.csv in {path}.")
    header = 0 if pos_file.name == "tissue_positions.csv" else None
    pos = pd.read_csv(pos_file, header=header)
    if header is None:
        pos.columns = ["barcode", "in_tissue", "array_row", "array_col",
                       "pxl_row_in_fullres", "pxl_col_in_fullres"]
    pos = pos.rename(columns={"barcode": "cell", "pxl_col_in_fullres": "x",
                              "pxl_row_in_fullres": "y"})
    pos["cell"] = pos["cell"].astype(str)
    cells = [str(c) for c in cells]

    if filter_by_tissue and "in_tissue" in pos.columns:
        pos = pos[pos["in_tissue"].astype(int) == 1]
        keep = set(pos["cell"])
        idx = [i for i, c in enumerate(cells) if c in keep]
        counts = counts[:, idx]              # drop off-tissue spots from the matrix too
        cells = [cells[i] for i in idx]

    coords = pos[["cell", "x", "y"]]
    obj = _build_spatial_object(counts, feats, cells, coords, assay, project,
                                image_name=slice_name)

    if not image:
        return obj

    sf_file = path / "spatial" / "scalefactors_json.json"
    sf = read_scale_factors(sf_file) if sf_file.exists() else None
    read = read_tissue_image(path / "spatial", resolution=image_resolution)
    if read is None and sf is None:
        return obj                      # nothing image-ish on disk; plain FOV is right
    img, res = read if read is not None else (None, image_resolution)
    obj.images = {
        name: VisiumV2.from_fov(fov, image=img, scale_factors=sf, image_resolution=res)
        for name, fov in obj.images.items()
    }
    return obj

load_cosmx

load_cosmx(path: Union[str, Path], expr_file: Optional[str] = None, meta_file: Optional[str] = None, assay: str = 'Nanostring', fov: str = 'fov', fov_column: Optional[str] = None, project: str = 'CosMx')

Load NanoString CosMx output (exprMat + metadata CSVs) into a Truecell object.

Builds what Seurat's LoadNanostring builds from the same files, for the expression matrix and the cell centroids:

• cells are named <cell_ID>_<fov>, as ReadNanostring names them;
• the cell_ID 0 row each FOV carries, for transcripts outside any cell, is dropped, and so is every cell with no counts;
• every column but fov and cell_ID is a feature, the negative probes included;
• one image, named fov (LoadNanostring's argument), holds each cell's CenterX_global_px / CenterY_global_px, for the cells in the object.

LoadNanostring also loads the cell polygons (*-polygons.csv) and the transcript coordinates (*_tx_file.csv) into that image; load_cosmx reads neither. truecell carries the metadata file's columns into meta_data, which Seurat does not. fov_column splits the cells into one image per value of that metadata column instead, named by the value.

expr_file / meta_file default to the files ReadNanostring finds in path: the last name, in sorted order, matching *_exprMat_file.csv and *_metadata_file.csv.

Checked against LoadNanostring (Seurat 5.5.1) on a synthetic bundle in CosMx's file layout; not compared with R on a real run.

Source code in truecell/spatial/loaders.py

def load_cosmx(
    path: Union[str, Path],
    expr_file: Optional[str] = None,
    meta_file: Optional[str] = None,
    assay: str = "Nanostring",
    fov: str = "fov",
    fov_column: Optional[str] = None,
    project: str = "CosMx",
):
    """Load NanoString CosMx output (exprMat + metadata CSVs) into a Truecell object.

    Builds what Seurat's ``LoadNanostring`` builds from the same files, for the
    expression matrix and the cell centroids:

    * cells are named ``<cell_ID>_<fov>``, as ``ReadNanostring`` names them;
    * the ``cell_ID`` 0 row each FOV carries, for transcripts outside any cell, is
      dropped, and so is every cell with no counts;
    * every column but ``fov`` and ``cell_ID`` is a feature, the negative probes
      included;
    * one image, named ``fov`` (``LoadNanostring``'s argument), holds each cell's
      ``CenterX_global_px`` / ``CenterY_global_px``, for the cells in the object.

    ``LoadNanostring`` also loads the cell polygons (``*-polygons.csv``) and the
    transcript coordinates (``*_tx_file.csv``) into that image; ``load_cosmx``
    reads neither. truecell carries the metadata file's columns into
    ``meta_data``, which Seurat does not. ``fov_column`` splits the cells into one
    image per value of that metadata column instead, named by the value.

    ``expr_file`` / ``meta_file`` default to the files ``ReadNanostring`` finds in
    ``path``: the last name, in sorted order, matching ``*_exprMat_file.csv`` and
    ``*_metadata_file.csv``.

    Checked against ``LoadNanostring`` (Seurat 5.5.1) on a synthetic bundle in
    CosMx's file layout; not compared with R on a real run.
    """
    path = Path(path)
    expr = (Path(expr_file) if expr_file
            else _file_like_r(path, r"[_a-zA-Z0-9]*_exprMat_file.csv"))
    meta = (Path(meta_file) if meta_file
            else _file_like_r(path, r"[_a-zA-Z0-9]*_metadata_file.csv"))
    if expr is None or meta is None:
        raise FileNotFoundError(
            f"Could not locate *_exprMat_file.csv / *_metadata_file.csv in {path}.")

    edf = pd.read_csv(expr)
    mdf = pd.read_csv(meta)
    for label, df in (("exprMat", edf), ("metadata", mdf)):
        missing = [c for c in ("fov", "cell_ID") if c not in df.columns]
        if missing:
            raise ValueError(f"The CosMx {label} file has no {' or '.join(missing)} column.")

    edf = edf[edf["cell_ID"] != 0]
    gene_cols = [c for c in edf.columns if c not in ("fov", "cell_ID")]
    counts = sp.csc_matrix(edf[gene_cols].to_numpy(dtype=float).T)   # genes × cells
    cell_ids = _cosmx_cell_names(edf)
    has_counts = np.asarray(counts.sum(axis=0)).ravel() != 0
    counts, cell_ids = counts[:, has_counts], cell_ids[has_counts]

    mdf = mdf.copy()
    mdf["cell"] = _cosmx_cell_names(mdf)
    mcoord = mdf.rename(columns={"CenterX_global_px": "x", "CenterY_global_px": "y"})
    if not {"x", "y"} <= set(mcoord.columns):
        raise ValueError("The CosMx metadata file has no CenterX_global_px / "
                         "CenterY_global_px columns.")
    split = fov_column if (fov_column and fov_column in mcoord.columns) else None
    coords = mcoord[["cell", "x", "y"] + ([split] if split else [])]
    return _build_spatial_object(counts, gene_cols, list(cell_ids), coords,
                                 assay, project, fov=split,
                                 meta_data=mdf.set_index("cell"),
                                 image_name=fov, radius=_auto_radius(coords))

load_merscope

load_merscope(path: Union[str, Path], expr_file: Optional[str] = None, meta_file: Optional[str] = None, assay: str = 'Vizgen', fov: str = 'fov', fov_column: Optional[str] = None, project: str = 'MERSCOPE', keep_controls: bool = False)

Load a Vizgen MERSCOPE output into a Truecell object with images.

Expects, in path: * cell_by_gene.csv — cell × gene counts (leading column = cell id) * cell_metadata.csv — with center_x / center_y (and usually fov, volume)

Builds what Seurat's LoadVizgen builds from the same files, for the expression matrix and the cell centroids:

• cell ids are the leading column of both files, as ReadVizgen reads them;
• features matching ^Blank-, the blank barcodes, are dropped as LoadVizgen drops them. The match is case-sensitive, so blank-3 stays; keep_controls=True keeps them all;
• cells with no counts stay, as they do in LoadVizgen;
• one image, named fov (LoadVizgen's argument), holds the center_x / center_y centroids of the cells in the object.

LoadVizgen also loads the cell boundaries (cell_boundaries/*.hdf5) and the transcripts of one z-plane into that image, and drops from the image any cell without a boundary; load_merscope reads neither and keeps those cells. truecell carries the metadata file's columns into meta_data, which Seurat does not. fov_column splits the cells into one image per value of that metadata column instead, named by the value.

expr_file / meta_file default to the files ReadVizgen finds in path: the last name, in sorted order, matching cell_by_gene*.csv and cell_metadata*.csv.

Checked against LoadVizgen (Seurat 5.5.1) on a synthetic bundle in MERSCOPE's file layout; not compared with R on a real run.

Source code in truecell/spatial/loaders.py

def load_merscope(
    path: Union[str, Path],
    expr_file: Optional[str] = None,
    meta_file: Optional[str] = None,
    assay: str = "Vizgen",
    fov: str = "fov",
    fov_column: Optional[str] = None,
    project: str = "MERSCOPE",
    keep_controls: bool = False,
):
    """Load a Vizgen MERSCOPE output into a Truecell object with images.

    Expects, in ``path``:
      * ``cell_by_gene.csv`` — cell × gene counts (leading column = cell id)
      * ``cell_metadata.csv`` — with ``center_x`` / ``center_y`` (and usually
        ``fov``, ``volume``)

    Builds what Seurat's ``LoadVizgen`` builds from the same files, for the
    expression matrix and the cell centroids:

    * cell ids are the leading column of both files, as ``ReadVizgen`` reads them;
    * features matching ``^Blank-``, the blank barcodes, are dropped as
      ``LoadVizgen`` drops them. The match is case-sensitive, so ``blank-3``
      stays; ``keep_controls=True`` keeps them all;
    * cells with no counts stay, as they do in ``LoadVizgen``;
    * one image, named ``fov`` (``LoadVizgen``'s argument), holds the
      ``center_x`` / ``center_y`` centroids of the cells in the object.

    ``LoadVizgen`` also loads the cell boundaries (``cell_boundaries/*.hdf5``) and
    the transcripts of one z-plane into that image, and drops from the image any
    cell without a boundary; ``load_merscope`` reads neither and keeps those
    cells. truecell carries the metadata file's columns into ``meta_data``, which
    Seurat does not. ``fov_column`` splits the cells into one image per value of
    that metadata column instead, named by the value.

    ``expr_file`` / ``meta_file`` default to the files ``ReadVizgen`` finds in
    ``path``: the last name, in sorted order, matching ``cell_by_gene*.csv`` and
    ``cell_metadata*.csv``.

    Checked against ``LoadVizgen`` (Seurat 5.5.1) on a synthetic bundle in
    MERSCOPE's file layout; not compared with R on a real run.
    """
    path = Path(path)
    expr = (Path(expr_file) if expr_file
            else _file_like_r(path, r"cell_by_gene[_a-zA-Z0-9]*.csv"))
    meta = (Path(meta_file) if meta_file
            else _file_like_r(path, r"cell_metadata[_a-zA-Z0-9]*.csv"))
    if expr is None or meta is None:
        raise FileNotFoundError(
            f"Could not locate cell_by_gene.csv / cell_metadata.csv in {path}."
        )

    edf = pd.read_csv(expr)
    mdf = pd.read_csv(meta)

    gene_cols = [str(c) for c in edf.columns[1:]]
    if not keep_controls:
        gene_cols = [g for g in gene_cols if not re.match(r"Blank-", g)]
    if not gene_cols:
        raise ValueError(f"No gene columns found in {expr}.")
    cell_ids = edf.iloc[:, 0].astype(str).to_numpy()
    counts = sp.csc_matrix(edf[gene_cols].to_numpy(dtype=float).T)   # genes × cells

    mdf = mdf.copy()
    id_column = mdf.columns[0]
    mdf["cell"] = mdf[id_column].astype(str)
    mcoord = mdf.rename(columns={"center_x": "x", "center_y": "y"})
    if not {"x", "y"} <= set(mcoord.columns):
        raise ValueError("cell_metadata must contain center_x / center_y columns.")
    split = fov_column if (fov_column and fov_column in mcoord.columns) else None
    coords = mcoord[["cell", "x", "y"] + ([split] if split else [])]
    meta_data = mdf.set_index("cell").drop(columns=[id_column], errors="ignore")
    return _build_spatial_object(counts, gene_cols, list(cell_ids), coords,
                                 assay, project, fov=split, meta_data=meta_data,
                                 image_name=fov, radius=_auto_radius(coords))

Containers

FOV

FOV(boundaries: Optional[dict[str, BoundaryType]] = None, molecules: Optional[dict[str, Molecules]] = None, coords_x_orientation: str = 'horizontal', assay: str = '', key: str = 'fov_', misc: Optional[dict] = None)

Bases: SpatialImage

Field of view — modern container for spatially-resolved single-cell coordinates.

Mirrors R's FOV class from fov.R. Can hold multiple segmentation boundaries (Centroids, Segmentation) and molecule-level FISH data (Molecules).

Slots

• molecules (dict[str, Molecules])
• boundaries (dict[str, BoundaryType])
• coords_x_orientation (str) — which axis x maps to in visualisation

Source code in truecell/spatial/fov.py

def __init__(
    self,
    boundaries: Optional[dict[str, BoundaryType]] = None,
    molecules: Optional[dict[str, Molecules]] = None,
    coords_x_orientation: str = "horizontal",
    assay: str = "",
    key: str = "fov_",
    misc: Optional[dict] = None,
) -> None:
    super().__init__(assay=assay, key=key, misc=misc)
    self.boundaries: dict[str, BoundaryType] = boundaries or {}
    self.molecules: dict[str, Molecules] = molecules or {}
    self.coords_x_orientation = coords_x_orientation
    # Default boundary is first entry by insertion order
    self._default_boundary: Optional[str] = (
        next(iter(self.boundaries)) if self.boundaries else None
    )

overlay

overlay(query: 'FOV') -> list[str]

Return cells from self whose centroids fall within any boundary of query.

Source code in truecell/spatial/fov.py

def overlay(self, query: "FOV") -> list[str]:
    """Return cells from self whose centroids fall within any boundary of query."""
    query_coords = query.get_tissue_coordinates()
    if query_coords.empty:
        return []
    x_min, x_max = query_coords["x"].min(), query_coords["x"].max()
    y_min, y_max = query_coords["y"].min(), query_coords["y"].max()
    self_coords = self.get_tissue_coordinates()
    mask = (
        (self_coords["x"] >= x_min)
        & (self_coords["x"] <= x_max)
        & (self_coords["y"] >= y_min)
        & (self_coords["y"] <= y_max)
    )
    return list(self_coords.index[mask])

create_fov

create_fov(coords: DataFrame, type_: str = 'centroids', nsides: int = 0, radius: Optional[float] = None, theta: Optional[float] = None, assay: str = '', key: str = 'fov_') -> FOV

Create an FOV from a coordinate DataFrame.

type_ : 'centroids' | 'segmentation' | 'molecules'

Source code in truecell/spatial/fov.py

def create_fov(
    coords: pd.DataFrame,
    type_: str = "centroids",
    nsides: int = 0,
    radius: Optional[float] = None,
    theta: Optional[float] = None,
    assay: str = "",
    key: str = "fov_",
) -> FOV:
    """Create an FOV from a coordinate DataFrame.

    type_ : 'centroids' | 'segmentation' | 'molecules'
    """
    boundary: BoundaryType
    if type_ == "centroids":
        boundary = create_centroids(coords, nsides=nsides, radius=radius, theta=theta, assay=assay)
        return FOV(boundaries={"centroids": boundary}, assay=assay, key=key)
    elif type_ == "segmentation":
        boundary = create_segmentation(coords, assay=assay)
        return FOV(boundaries={"segmentation": boundary}, assay=assay, key=key)
    elif type_ == "molecules":
        mol = create_molecules(coords, assay=assay)
        return FOV(molecules={"molecules": mol}, assay=assay, key=key)
    else:
        raise ValueError(f"Unknown type_ '{type_}'. Choose centroids, segmentation, or molecules.")

create_fovs

create_fovs(coords: DataFrame, fov: Optional[Union[str, 'pd.Series', Sequence]] = None, assay: str = '', default_name: str = 'fov') -> dict[str, FOV]

Build a {name: FOV} dict of centroid FOVs from a coordinate frame.

coords must have columns x, y, cell. When fov is given (a column name in coords or a per-row array of labels) the cells are split into one FOV per distinct label — matching a multi-FOV Xenium/CosMx run. The FOVs come in order of first appearance, or in category order when the labels are categorical, which is how as_anndata records the order of the images. Otherwise a single FOV named default_name is returned.

Shared by the spatial loaders and from_anndata so both build identical, accessor-ready seurat.images structures.

Source code in truecell/spatial/fov.py

def create_fovs(
    coords: pd.DataFrame,
    fov: Optional[Union[str, "pd.Series", Sequence]] = None,
    assay: str = "",
    default_name: str = "fov",
) -> dict[str, FOV]:
    """Build a ``{name: FOV}`` dict of centroid FOVs from a coordinate frame.

    ``coords`` must have columns ``x, y, cell``. When ``fov`` is given (a column
    name in ``coords`` or a per-row array of labels) the cells are split into one
    FOV per distinct label — matching a multi-FOV Xenium/CosMx run. The FOVs come
    in order of first appearance, or in category order when the labels are
    categorical, which is how ``as_anndata`` records the order of the images.
    Otherwise a single FOV named ``default_name`` is returned.

    Shared by the spatial loaders and ``from_anndata`` so both build identical,
    accessor-ready ``seurat.images`` structures.
    """
    coords = coords.copy()
    for col in ("x", "y", "cell"):
        if col not in coords.columns:
            raise ValueError(f"coords must have a '{col}' column.")

    if fov is None:
        return {default_name: create_fov(coords, type_="centroids", assay=assay,
                                         key=f"{default_name}_")}

    if isinstance(fov, str):
        if fov not in coords.columns:
            raise ValueError(f"fov column '{fov}' not in coords.")
        series = coords[fov]
    else:
        series = pd.Series(fov)
        if len(series) != len(coords):
            raise ValueError("fov label length must match number of rows in coords.")
    labels = series.astype(str).to_numpy()
    names = list(pd.unique(labels))
    if isinstance(series.dtype, pd.CategoricalDtype):
        present = set(names)
        ordered = [c for c in series.cat.categories.astype(str) if c in present]
        names = ordered + [n for n in names if n not in set(ordered)]

    images: dict[str, FOV] = {}
    for name in names:
        sub = coords.loc[labels == name, ["x", "y", "cell"]]
        safe = str(name).replace(" ", "_")
        images[safe] = create_fov(sub, type_="centroids", assay=assay, key=f"{safe}_")
    return images

Centroids

Centroids(coords: DataFrame, nsides: int = 0, radius: Optional[float] = None, theta: Optional[float] = None, assay: str = '', key: str = 'centroids_', misc: Optional[dict] = None)

Bases: SpatialImage

Cell centroid coordinates.

Mirrors R's Centroids class from centroids.R. Stores one (x, y) coordinate per cell representing the center of that cell.

Slots

• _coords (pd.DataFrame) — columns: x, y, cell
• nsides (int) — number of polygon sides (0 = circle)
• radius_ (Optional[float]) — spot radius (for spot-based technologies)
• theta_ (Optional[float]) — angle offset

Source code in truecell/spatial/centroids.py

def __init__(
    self,
    coords: pd.DataFrame,
    nsides: int = 0,
    radius: Optional[float] = None,  # None → SeuratObject's .AutoRadius
    theta: Optional[float] = None,
    assay: str = "",
    key: str = "centroids_",
    misc: Optional[dict] = None,
) -> None:
    super().__init__(assay=assay, key=key, misc=misc)
    coords = coords.copy()
    for col in ("x", "y", "cell"):
        if col not in coords.columns:
            raise ValueError(f"coords must have a '{col}' column.")
    self._coords = coords[["x", "y", "cell"]].copy()
    self.nsides = nsides
    self.radius_ = radius if radius is not None else _auto_radius(self._coords)
    self.theta_ = theta

create_centroids

create_centroids(coords: DataFrame, nsides: int = 0, radius: Optional[float] = None, theta: Optional[float] = None, assay: str = '', key: str = 'centroids_') -> Centroids

Source code in truecell/spatial/centroids.py

def create_centroids(
    coords: pd.DataFrame,
    nsides: int = 0,
    radius: Optional[float] = None,
    theta: Optional[float] = None,
    assay: str = "",
    key: str = "centroids_",
) -> Centroids:
    return Centroids(
        coords=coords,
        nsides=nsides,
        radius=radius,
        theta=theta,
        assay=assay,
        key=key,
    )

Segmentation

Segmentation(coords: DataFrame, assay: str = '', key: str = 'segmentation_', misc: Optional[dict] = None)

Bases: SpatialImage

Cell boundary polygon coordinates.

Mirrors R's Segmentation class from segmentation.R. Each cell may have multiple (x, y) polygon vertices.

Slots

• _coords (pd.DataFrame) — columns: x, y, cell (multiple rows per cell)

Source code in truecell/spatial/segmentation.py

def __init__(
    self,
    coords: pd.DataFrame,
    assay: str = "",
    key: str = "segmentation_",
    misc: Optional[dict] = None,
) -> None:
    super().__init__(assay=assay, key=key, misc=misc)
    for col in ("x", "y", "cell"):
        if col not in coords.columns:
            raise ValueError(f"coords must have a '{col}' column.")
    self._coords = _close_rings(coords[["x", "y", "cell"]].copy())

as_centroids

as_centroids() -> Centroids

One point per cell, where SeuratObject puts it.

Mirrors as(segmentation, "Centroids"): the area centroid of each cell's ring, or its first vertex when the ring has no area. The radius is left to .AutoRadius, as CreateCentroids leaves it.

Checked against SeuratObject 5.4.0 on 40 random polygons, where it agrees to 5.2e-16 relative, and exactly on rings with no area. The last bits differ because CRAN's arm64 build fuses multiply-adds, which truecell's plain IEEE arithmetic does not.

Source code in truecell/spatial/segmentation.py

def as_centroids(self) -> Centroids:
    """One point per cell, where SeuratObject puts it.

    Mirrors ``as(segmentation, "Centroids")``: the area centroid of each cell's
    ring, or its first vertex when the ring has no area. The radius is left to
    ``.AutoRadius``, as ``CreateCentroids`` leaves it.

    Checked against SeuratObject 5.4.0 on 40 random polygons, where it agrees to
    5.2e-16 relative, and exactly on rings with no area. The last bits differ
    because CRAN's arm64 build fuses multiply-adds, which truecell's plain IEEE
    arithmetic does not.
    """
    return Centroids(coords=_ring_centroids(self._coords), assay=self.assay,
                     misc=dict(self.misc))

simplify

simplify(tol: float = 0.5) -> 'Segmentation'

Reduce polygon vertex count by removing vertices closer than tol.

Simple Douglas–Peucker–style approximation per cell.

Source code in truecell/spatial/segmentation.py

def simplify(self, tol: float = 0.5) -> "Segmentation":
    """Reduce polygon vertex count by removing vertices closer than tol.

    Simple Douglas–Peucker–style approximation per cell.
    """
    groups = []
    for cell_id, group in self._coords.groupby("cell", sort=False):
        pts = group[["x", "y"]].values
        if len(pts) <= 3:
            groups.append(group)
            continue
        # keep[i] = True means retain pts[i]; always keep first and last
        keep = np.zeros(len(pts), dtype=bool)
        keep[0] = True
        keep[-1] = True
        diffs = np.linalg.norm(np.diff(pts, axis=0), axis=1)  # len(pts)-1
        keep[1:-1] = diffs[:-1] > tol
        groups.append(group.iloc[np.where(keep)[0]])

    new_coords = pd.concat(groups, ignore_index=True)
    return Segmentation(
        coords=new_coords,
        assay=self.assay,
        key=self._key,
        misc=dict(self.misc),
    )

create_segmentation

create_segmentation(coords: DataFrame, assay: str = '', key: str = 'segmentation_') -> Segmentation

Source code in truecell/spatial/segmentation.py

def create_segmentation(
    coords: pd.DataFrame,
    assay: str = "",
    key: str = "segmentation_",
) -> Segmentation:
    return Segmentation(coords=coords, assay=assay, key=key)

Molecules

Molecules(coords: DataFrame, assay: str = '', key: str = 'molecules_', misc: Optional[dict] = None)

Bases: SpatialImage

Spatially-resolved molecule (FISH) data.

Mirrors R's Molecules class from molecules.R. Each row represents one molecule detection event.

Slots

• _coords (pd.DataFrame) — columns: x, y, gene, [cell] 'cell' is optional (not all FISH protocols assign cells)

Source code in truecell/spatial/molecules.py

def __init__(
    self,
    coords: pd.DataFrame,
    assay: str = "",
    key: str = "molecules_",
    misc: Optional[dict] = None,
) -> None:
    super().__init__(assay=assay, key=key, misc=misc)
    for col in ("x", "y", "gene"):
        if col not in coords.columns:
            raise ValueError(f"coords must have a '{col}' column.")
    keep = ["x", "y", "gene"]
    if "cell" in coords.columns:
        keep.append("cell")
    self._coords = coords[keep].copy()

create_molecules

create_molecules(coords: DataFrame, assay: str = '', key: str = 'molecules_') -> Molecules

Source code in truecell/spatial/molecules.py

def create_molecules(
    coords: pd.DataFrame,
    assay: str = "",
    key: str = "molecules_",
) -> Molecules:
    return Molecules(coords=coords, assay=assay, key=key)

SpatialImage

SpatialImage(assay: str = '', key: str = 'image_', misc: Optional[dict] = None)

Bases: KeyMixin, ABC

Abstract base class for all spatial image objects.

Mirrors R's SpatialImage virtual class from spatial.R. Subclasses must implement: cells, dim, get_tissue_coordinates, rename_cells, subset.

Slots

• assay (str) — associated assay name
• misc (dict) — miscellaneous storage
• _key (str) — inherited from KeyMixin

Source code in truecell/spatial/base.py

def __init__(
    self,
    assay: str = "",
    key: str = "image_",
    misc: Optional[dict] = None,
) -> None:
    self.assay = assay
    self._key = key
    self.misc = misc or {}

VisiumV2

VisiumV2(boundaries: Optional[dict] = None, molecules: Optional[dict[str, Molecules]] = None, image: Optional[ndarray] = None, scale_factors: Optional[ScaleFactors] = None, image_resolution: str = 'lowres', coords_x_orientation: str = 'horizontal', assay: str = '', key: str = 'slice1_', misc: Optional[dict] = None)

Bases: FOV

An FOV that also carries the Visium tissue image and its scale factors.

Mirrors Seurat v5's VisiumV2. Behaves as a normal FOV everywhere (its coordinates stay in fullres pixels); the extra slots are what spatial_dim_plot / spatial_feature_plot need to draw spots on tissue.

Slots

• image (np.ndarray) — (H, W[, C]), or None
• scale_factors (ScaleFactors) — or None
• image_resolution (str) — 'hires' or 'lowres', which image is stored

Source code in truecell/spatial/visium.py

def __init__(
    self,
    boundaries: Optional[dict] = None,
    molecules: Optional[dict[str, Molecules]] = None,
    image: Optional[np.ndarray] = None,
    scale_factors: Optional[ScaleFactors] = None,
    image_resolution: str = "lowres",
    coords_x_orientation: str = "horizontal",
    assay: str = "",
    key: str = "slice1_",
    misc: Optional[dict] = None,
) -> None:
    super().__init__(
        boundaries=boundaries,
        molecules=molecules,
        coords_x_orientation=coords_x_orientation,
        assay=assay,
        key=key,
        misc=misc,
    )
    self.image = image
    self.scale_factors = scale_factors
    self.image_resolution = image_resolution

from_fov classmethod

from_fov(fov: FOV, image: Optional[ndarray] = None, scale_factors: Optional[ScaleFactors] = None, image_resolution: str = 'lowres') -> 'VisiumV2'

Upgrade a plain FOV in place-ish: same boundaries, plus the image.

Source code in truecell/spatial/visium.py

@classmethod
def from_fov(
    cls,
    fov: FOV,
    image: Optional[np.ndarray] = None,
    scale_factors: Optional[ScaleFactors] = None,
    image_resolution: str = "lowres",
) -> "VisiumV2":
    """Upgrade a plain FOV in place-ish: same boundaries, plus the image."""
    obj = cls(
        boundaries=dict(fov.boundaries),
        molecules=dict(fov.molecules),
        image=image,
        scale_factors=scale_factors,
        image_resolution=image_resolution,
        coords_x_orientation=fov.coords_x_orientation,
        assay=fov.assay,
        key=fov._key,
        misc=dict(fov.misc),
    )
    # Spot boundaries know their own radius, in the same fullres pixel space.
    r = obj.radius()
    if r is not None:
        for b in obj.boundaries.values():
            if isinstance(b, Centroids):
                b.radius_ = r
    return obj

get_image

get_image() -> Optional[ndarray]

The stored tissue image, or None when the bundle had none.

Source code in truecell/spatial/visium.py

def get_image(self) -> Optional[np.ndarray]:
    """The stored tissue image, or None when the bundle had none."""
    return self.image

radius

radius() -> Optional[float]

Spot radius in full-resolution pixels (half the spot diameter).

Source code in truecell/spatial/visium.py

def radius(self) -> Optional[float]:
    """Spot radius in full-resolution pixels (half the spot diameter)."""
    if self.scale_factors is None or not np.isfinite(self.scale_factors.spot):
        return None
    return float(self.scale_factors.spot) / 2.0

scale_factor

scale_factor(resolution: Optional[str] = None) -> float

Fullres → image-pixel multiplier for the stored (or given) resolution.

Source code in truecell/spatial/visium.py

def scale_factor(self, resolution: Optional[str] = None) -> float:
    """Fullres → image-pixel multiplier for the stored (or given) resolution."""
    if self.scale_factors is None:
        return 1.0
    return self.scale_factors.scale_factor(resolution or self.image_resolution)

scale_coordinates

scale_coordinates(cells: Optional[list[str]] = None, resolution: Optional[str] = None) -> DataFrame

Tissue coordinates rescaled into the pixel space of the stored image.

Same frame as get_tissue_coordinates, with x/y multiplied by the scale factor — i.e. ready to overlay on get_image.

Source code in truecell/spatial/visium.py

def scale_coordinates(
    self,
    cells: Optional[list[str]] = None,
    resolution: Optional[str] = None,
) -> pd.DataFrame:
    """Tissue coordinates rescaled into the pixel space of the stored image.

    Same frame as :meth:`get_tissue_coordinates`, with ``x``/``y`` multiplied
    by the scale factor — i.e. ready to overlay on :meth:`get_image`.
    """
    coords = self.get_tissue_coordinates(cells=cells).copy()
    f = self.scale_factor(resolution)
    for col in ("x", "y"):
        if col in coords.columns:
            coords[col] = coords[col].astype(float) * f
    return coords

spot_radius

spot_radius(resolution: Optional[str] = None) -> Optional[float]

Spot radius in the pixel space of the stored image (None if unknown).

Source code in truecell/spatial/visium.py

def spot_radius(self, resolution: Optional[str] = None) -> Optional[float]:
    """Spot radius in the pixel space of the stored image (None if unknown)."""
    r = self.radius()
    return None if r is None else r * self.scale_factor(resolution)

ScaleFactors dataclass

ScaleFactors(spot: float, fiducial: float, hires: float, lowres: float)

The four values in a Visium spatial/scalefactors_json.json.

spot : spot diameter in full-resolution pixels fiducial : fiducial-marker diameter in full-resolution pixels hires : multiply a fullres coordinate by this to land on tissue_hires_image.png lowres : ditto for tissue_lowres_image.png

from_dict classmethod

from_dict(d) -> ScaleFactors

From scalefactors_json.json's keys. A missing key is NaN.

Source code in truecell/spatial/visium.py

@classmethod
def from_dict(cls, d) -> ScaleFactors:
    """From ``scalefactors_json.json``'s keys. A missing key is NaN."""
    return cls(
        spot=float(d.get("spot_diameter_fullres", np.nan)),
        fiducial=float(d.get("fiducial_diameter_fullres", np.nan)),
        hires=float(d.get("tissue_hires_scalef", np.nan)),
        lowres=float(d.get("tissue_lowres_scalef", np.nan)),
    )

to_dict

to_dict() -> dict[str, float]

Under scalefactors_json.json's keys, as Scanpy keeps them in uns.

Source code in truecell/spatial/visium.py

def to_dict(self) -> dict[str, float]:
    """Under ``scalefactors_json.json``'s keys, as Scanpy keeps them in ``uns``."""
    return {
        "spot_diameter_fullres": self.spot,
        "fiducial_diameter_fullres": self.fiducial,
        "tissue_hires_scalef": self.hires,
        "tissue_lowres_scalef": self.lowres,
    }

Spatial analysis

get_tissue_coordinates

get_tissue_coordinates(seurat, image: Optional[Union[str, Sequence[str]]] = None) -> DataFrame

Return centroid coordinates for the object, concatenated across images.

Mirrors object-level GetTissueCoordinates. Returns a DataFrame with columns x, y, cell, image (one row per cell per image).

Source code in truecell/spatial/analysis.py

def get_tissue_coordinates(
    seurat,
    image: Optional[Union[str, Sequence[str]]] = None,
) -> pd.DataFrame:
    """Return centroid coordinates for the object, concatenated across images.

    Mirrors object-level ``GetTissueCoordinates``. Returns a DataFrame with
    columns ``x, y, cell, image`` (one row per cell per image).
    """
    frames = []
    for nm in _image_names(seurat, image):
        coords = seurat.images[nm].get_tissue_coordinates()
        c = coords.copy()
        c["cell"] = list(coords.index)
        c["image"] = nm
        frames.append(c.reset_index(drop=True))
    if not frames:
        return pd.DataFrame(columns=["x", "y", "cell", "image"])
    return pd.concat(frames, ignore_index=True)[["x", "y", "cell", "image"]]

spatial_knn

spatial_knn(coords: ndarray, k: int = 10, query: Optional[ndarray] = None) -> tuple[ndarray, ndarray]

k-nearest-neighbour distances and indices on a set of coordinates.

Mirrors FNN::get.knn (query=None, self excluded) and FNN::get.knnx (query given, searched against coords).

Returns (distances, indices) each shape (n_query, k); indices point into coords.

Source code in truecell/spatial/analysis.py

def spatial_knn(
    coords: np.ndarray,
    k: int = 10,
    query: Optional[np.ndarray] = None,
) -> tuple[np.ndarray, np.ndarray]:
    """k-nearest-neighbour distances and indices on a set of coordinates.

    Mirrors ``FNN::get.knn`` (``query=None``, self excluded) and
    ``FNN::get.knnx`` (``query`` given, searched against ``coords``).

    Returns ``(distances, indices)`` each shape ``(n_query, k)``; indices point
    into ``coords``.
    """
    coords = np.asarray(coords, dtype=float)
    tree = cKDTree(coords)
    if query is None:
        d, i = tree.query(coords, k=min(k + 1, len(coords)))
        d = np.atleast_2d(d)
        i = np.atleast_2d(i)
        return d[:, 1:], i[:, 1:]           # drop self
    query = np.asarray(query, dtype=float)
    d, i = tree.query(query, k=k)
    if d.ndim == 1:                          # k == 1 → make 2-D (n_query, 1)
        d = d[:, None]
        i = i[:, None]
    return d, i

nearest_neighbor_distance

nearest_neighbor_distance(seurat, group_by: str, reference, target=None, image: Optional[Union[str, Sequence[str]]] = None) -> DataFrame

Distance from each reference cell to the nearest target cell.

Computed per image (so distances never cross FOVs), then concatenated. If target is None it defaults to reference (nearest same-type cell, self excluded) — this is the mast-to-nearest-mast idiom.

Returns a DataFrame: cell, image, reference, target, distance.

Source code in truecell/spatial/analysis.py

def nearest_neighbor_distance(
    seurat,
    group_by: str,
    reference,
    target=None,
    image: Optional[Union[str, Sequence[str]]] = None,
) -> pd.DataFrame:
    """Distance from each ``reference`` cell to the nearest ``target`` cell.

    Computed per image (so distances never cross FOVs), then concatenated.
    If ``target`` is None it defaults to ``reference`` (nearest same-type cell,
    self excluded) — this is the mast-to-nearest-mast idiom.

    Returns a DataFrame: ``cell, image, reference, target, distance``.
    """
    target = reference if target is None else target
    labels = seurat.meta_data[group_by].astype(str)
    same = str(target) == str(reference)
    out = []
    for nm in _image_names(seurat, image):
        coords = seurat.images[nm].get_tissue_coordinates()
        cells = list(coords.index)
        lab = labels.reindex(cells).astype(str).to_numpy()
        xy = coords[["x", "y"]].to_numpy(dtype=float)
        ref_mask = lab == str(reference)
        tgt_mask = lab == str(target)
        if ref_mask.sum() == 0 or tgt_mask.sum() < (2 if same else 1):
            continue
        tree = cKDTree(xy[tgt_mask])
        if same:
            dist = tree.query(xy[ref_mask], k=2)[0][:, 1]
        else:
            dist = tree.query(xy[ref_mask], k=1)[0]
            dist = np.atleast_1d(dist)
        out.append(pd.DataFrame({
            "cell": [c for c, m in zip(cells, ref_mask) if m],
            "image": nm, "reference": str(reference),
            "target": str(target), "distance": dist,
        }))
    if not out:
        return pd.DataFrame(columns=["cell", "image", "reference", "target", "distance"])
    return pd.concat(out, ignore_index=True)

local_neighborhood

local_neighborhood(seurat, group_by: str, reference=None, k: int = 10, image: Optional[Union[str, Sequence[str]]] = None) -> DataFrame

Composition of each cell's k nearest neighbours (self excluded).

For every reference cell (all cells if reference is None) returns the count and proportion of each group_by category among its k spatial neighbours. Columns: cell, image, n_<group>…, prop_<group>….

prop_<reference> is the "local density" of that type; the n_<group> columns are the raw neighbourhood composition.

Source code in truecell/spatial/analysis.py

def local_neighborhood(
    seurat,
    group_by: str,
    reference=None,
    k: int = 10,
    image: Optional[Union[str, Sequence[str]]] = None,
) -> pd.DataFrame:
    """Composition of each cell's ``k`` nearest neighbours (self excluded).

    For every ``reference`` cell (all cells if ``reference`` is None) returns
    the count and proportion of each ``group_by`` category among its k spatial
    neighbours. Columns: ``cell, image, n_<group>…, prop_<group>…``.

    ``prop_<reference>`` is the "local density" of that type; the ``n_<group>``
    columns are the raw neighbourhood composition.
    """
    labels = seurat.meta_data[group_by].astype(str)
    groups = sorted(labels.dropna().unique())
    out = []
    for nm in _image_names(seurat, image):
        coords = seurat.images[nm].get_tissue_coordinates()
        cells = list(coords.index)
        if len(cells) < k + 1:
            continue
        lab = labels.reindex(cells).astype(str).to_numpy()
        xy = coords[["x", "y"]].to_numpy(dtype=float)
        ref_idx = (np.arange(len(cells)) if reference is None
                   else np.where(lab == str(reference))[0])
        if ref_idx.size == 0:
            continue
        _, idx = cKDTree(xy).query(xy[ref_idx], k=k + 1)
        neigh = lab[idx[:, 1:]]                       # (n_ref, k), self dropped
        df = pd.DataFrame({"cell": [cells[j] for j in ref_idx], "image": nm})
        for g in groups:
            df[f"n_{g}"] = (neigh == g).sum(axis=1)
        tot = df[[f"n_{g}" for g in groups]].sum(axis=1).replace(0, np.nan)
        for g in groups:
            df[f"prop_{g}"] = df[f"n_{g}"] / tot
        out.append(df)
    if not out:
        return pd.DataFrame(columns=["cell", "image"])
    return pd.concat(out, ignore_index=True)

build_niche_assay

build_niche_assay(seurat, group_by: str, image: Optional[Union[str, Sequence[str]]] = None, k: int = 20, niches: int = 4, assay_name: str = 'niche', cluster: bool = True, seed: int = 0) -> 'object'

Build a neighbourhood-composition assay and cluster cells into niches.

Mirrors Seurat v5's BuildNicheAssay: each cell's feature vector is the count of every group_by category among its k spatial neighbours. The composition matrix is stored as a new assay (assay_name); when cluster is True the proportions are k-means clustered into niches groups and written to meta_data['niches'].

Source code in truecell/spatial/analysis.py

def build_niche_assay(
    seurat,
    group_by: str,
    image: Optional[Union[str, Sequence[str]]] = None,
    k: int = 20,
    niches: int = 4,
    assay_name: str = "niche",
    cluster: bool = True,
    seed: int = 0,
) -> "object":
    """Build a neighbourhood-composition assay and cluster cells into niches.

    Mirrors Seurat v5's ``BuildNicheAssay``: each cell's feature vector is the
    count of every ``group_by`` category among its ``k`` spatial neighbours.
    The composition matrix is stored as a new assay (``assay_name``); when
    ``cluster`` is True the proportions are k-means clustered into ``niches``
    groups and written to ``meta_data['niches']``.
    """
    from ..assay5 import create_assay5_object

    comp = local_neighborhood(seurat, group_by, reference=None, k=k, image=image)
    if comp.empty:
        raise ValueError("No cells with enough neighbours to build a niche assay.")
    groups = [c[2:] for c in comp.columns if c.startswith("n_")]
    comp = comp.set_index("cell")
    mat = comp[[f"n_{g}" for g in groups]].to_numpy(dtype=float).T   # groups × cells
    cells = list(comp.index)

    assay = create_assay5_object(
        counts=sp.csc_matrix(mat), feature_names=groups, cell_names=cells,
        key=f"{assay_name.lower()}_",
    )
    seurat.assays[assay_name] = assay

    if cluster:
        try:
            from sklearn.cluster import KMeans
        except ImportError as e:  # pragma: no cover
            raise ImportError("scikit-learn is required for niche clustering "
                              "(pip install 'truecell[analysis]').") from e
        props = comp[[f"prop_{g}" for g in groups]].fillna(0).to_numpy()
        km = KMeans(n_clusters=niches, random_state=seed, n_init=10).fit(props)
        niche_lab = pd.Series([f"niche_{c + 1}" for c in km.labels_], index=cells)
        seurat.meta_data["niches"] = niche_lab.reindex(seurat.meta_data.index).values
    return seurat

find_spatially_variable_features

find_spatially_variable_features(seurat, features: Optional[list[str]] = None, method: str = 'moransi', k: int = 10, weights: str = 'inverse_square', assay: Optional[str] = None, layer: str = 'scale.data', image: Optional[Union[str, Sequence[str]]] = None, r_metric: float = 5.0, bandwidth: float = 1.0) -> DataFrame

Rank features by how spatially structured their expression is.

Mirrors R's FindSpatiallyVariableFeatures(obj, method = ...).

Parameters:

• features (Optional[list[str]], default: None ) –

  restrict to these genes (default: all in the assay). Passing the variable features keeps this fast on large panels.

• method (str, default: 'moransi' ) –

  "moransi" (default) or "markvariogram".

• weights (str, default: 'inverse_square' ) –

  moransi only — "inverse_square" (default) reproduces R exactly: 1/d² between every pair of cells, row-standardised. "knn" uses a k-nearest-neighbour graph instead, an approximation that trades R's O(n²) cost for O(nk). Prefer it only when the slide is too large to wait on; see the note below on how far the two answers drift apart.

• k (int, default: 10 ) –

  moransi only, weights="knn" only — neighbours per cell.

• r_metric (float, default: 5.0 ) –

  markvariogram only — the distance at which to read the variogram, in units of the median nearest-neighbour spacing. The default of 5 therefore means "five cells apart".

• bandwidth (float, default: 1.0 ) –

  markvariogram only — half-width of the distance band around r_metric, in the same units. Widen it if the slide is sparse and the band catches too few pairs to average over.

• layer (str, default: 'scale.data' ) –

  expression layer. "scale.data" by default, as in Seurat's object method, so only the features that were scaled are ranked; "data" ranks every gene on its log-normalized values.

• image (Optional[Union[str, Sequence[str]]], default: None ) –

  image name(s) to draw coordinates from (default: all).

Returns:

• DataFrame indexed by gene, sorted so that rank 1 is the most spatially –
• variable, with the columns for the chosen method: –
• ``moransi`` –

  moransi (I; +ve = spatially clustered), moransi_pval (two-sided, normality assumption), moransi_padj (Benjamini-Hochberg) and moransi_rank.

• ``markvariogram`` –

  markvariogram (γ at r_metric; lower = more spatially structured, ≈ 1 = none) and markvariogram_rank. There is no p-value — the variogram has no closed-form null, and R does not offer one either.

• The same columns are also written into the assay's feature-level metadata, –
• as ``find_variable_features`` does. –

Notes

Seurat ranks scale.data by default, and so does this. Moran's I does not change when a gene is shifted and rescaled, so data gives the same statistic except where ScaleData clipped a value at 10. Be aware that when a few strongly spatial genes dominate a cell's library size, log-normalization divides every gene by a spatially-structured total and leaks that structure into otherwise-flat genes — inflating their score. That is a property of compositional normalization, not of either statistic; it is negligible on real panels but can bite on small synthetic ones.

The moransi statistic matches R to ~1e-15 with the default weights. The p-value deliberately does not. R runs a 999-permutation test through Rfast2::moranI, which on a 248-gene panel returns just 14 distinct values and ties 233 genes at its 1/1025 floor — it cannot separate the most spatially variable gene from the two-hundredth. The normal-approximation p-value here is continuous, deterministic, and standard for Moran's I, so it is kept; matching R would cost information and gain nothing. weights="knn" changes the statistic itself, not just its cost. It is a decent approximation — on 2,000 Xenium cells it tracks R at Pearson 0.986 and recovers 46 of R's top 50 genes — but it runs a median 1.23× high, differs by up to 0.14 in absolute I, and agrees on only 7 of R's top 10. Close enough to mislead, not close enough to call parity, which is why it is no longer the default.

markvariogram differs from R's in two deliberate ways. R passes r.metric straight through to spatstat in raw coordinate units, so the same script gives different answers on pixel and micron coordinates; here r is measured in nearest-neighbour spacings and is scale-free. And γ is a kernel-weighted ratio estimator rather than spatstat's translation-corrected one, so the absolute γ values are close to but not identical with R's. The gene ranking — the thing the function is for — is what carries over.

Source code in truecell/spatial/variable_features.py

def find_spatially_variable_features(
    seurat,
    features: Optional[list[str]] = None,
    method: str = "moransi",
    k: int = 10,
    weights: str = "inverse_square",
    assay: Optional[str] = None,
    layer: str = "scale.data",
    image: Optional[Union[str, Sequence[str]]] = None,
    r_metric: float = 5.0,
    bandwidth: float = 1.0,
) -> pd.DataFrame:
    """Rank features by how spatially structured their expression is.

    Mirrors R's ``FindSpatiallyVariableFeatures(obj, method = ...)``.

    Parameters
    ----------
    features : restrict to these genes (default: all in the assay). Passing the
               variable features keeps this fast on large panels.
    method   : ``"moransi"`` (default) or ``"markvariogram"``.
    weights  : *moransi only* — ``"inverse_square"`` (default) reproduces R
               exactly: 1/d² between every pair of cells, row-standardised.
               ``"knn"`` uses a k-nearest-neighbour graph instead, an
               approximation that trades R's O(n²) cost for O(nk). Prefer it only
               when the slide is too large to wait on; see the note below on how
               far the two answers drift apart.
    k        : *moransi only, weights="knn" only* — neighbours per cell.
    r_metric : *markvariogram only* — the distance at which to read the
               variogram, **in units of the median nearest-neighbour spacing**.
               The default of 5 therefore means "five cells apart".
    bandwidth: *markvariogram only* — half-width of the distance band around
               ``r_metric``, in the same units. Widen it if the slide is sparse
               and the band catches too few pairs to average over.
    layer    : expression layer. ``"scale.data"`` by default, as in Seurat's
               object method, so only the features that were scaled are
               ranked; ``"data"`` ranks every gene on its log-normalized values.
    image    : image name(s) to draw coordinates from (default: all).

    Returns
    -------
    DataFrame indexed by gene, sorted so that rank 1 is the most spatially
    variable, with the columns for the chosen method:

    ``moransi``
        ``moransi`` (I; +ve = spatially clustered), ``moransi_pval`` (two-sided,
        normality assumption), ``moransi_padj`` (Benjamini-Hochberg) and
        ``moransi_rank``.
    ``markvariogram``
        ``markvariogram`` (γ at ``r_metric``; **lower** = more spatially
        structured, ≈ 1 = none) and ``markvariogram_rank``. There is no p-value —
        the variogram has no closed-form null, and R does not offer one either.

    The same columns are also written into the assay's feature-level metadata,
    as ``find_variable_features`` does.

    Notes
    -----
    Seurat ranks ``scale.data`` by default, and so does this. Moran's I does not
    change when a gene is shifted and rescaled, so ``data`` gives the same
    statistic except where ScaleData clipped a value at 10. Be aware that when a few
    strongly spatial genes dominate a cell's library size, log-normalization
    divides every gene by a spatially-structured total and leaks that structure
    into otherwise-flat genes — inflating their score. That is a property of
    compositional normalization, not of either statistic; it is negligible on
    real panels but can bite on small synthetic ones.

    The ``moransi`` **statistic** matches R to ~1e-15 with the default weights.
    The **p-value** deliberately does not. R runs a 999-permutation test through
    ``Rfast2::moranI``, which on a 248-gene panel returns just 14 distinct values
    and ties 233 genes at its 1/1025 floor — it cannot separate the most spatially
    variable gene from the two-hundredth. The normal-approximation p-value here is
    continuous, deterministic, and standard for Moran's I, so it is kept; matching
    R would cost information and gain nothing. ``weights="knn"`` changes the
    statistic itself, not just its cost. It is a decent approximation — on 2,000
    Xenium cells it tracks R at Pearson 0.986 and recovers 46 of R's top 50 genes
    — but it runs a median 1.23× high, differs by up to 0.14 in absolute I, and
    agrees on only 7 of R's top 10. Close enough to mislead, not close enough to
    call parity, which is why it is no longer the default.

    ``markvariogram`` differs from R's in two deliberate ways. R passes
    ``r.metric`` straight through to ``spatstat`` in raw coordinate units, so the
    same script gives different answers on pixel and micron coordinates; here r
    is measured in nearest-neighbour spacings and is scale-free. And γ is a
    kernel-weighted ratio estimator rather than ``spatstat``'s
    translation-corrected one, so the absolute γ values are close to but not
    identical with R's. The gene *ranking* — the thing the function is for — is
    what carries over.
    """
    from ..markers import _get_expression_layer

    if method not in METHODS:
        raise NotImplementedError(
            f"method={method!r} is not implemented; use one of {list(METHODS)}."
        )
    if weights not in WEIGHTS:
        raise ValueError(
            f"weights={weights!r} is not recognised; use one of {list(WEIGHTS)}."
        )

    coords = get_tissue_coordinates(seurat, image)
    if coords.empty:
        raise ValueError("Object has no spatial coordinates.")
    # A cell can appear once per image; keep its first placement.
    coords = coords.drop_duplicates(subset="cell")

    assay_name = assay or seurat.active_assay
    assay_obj = seurat.assays[assay_name]
    if not _has_layer(assay_obj, layer):
        raise ValueError(
            f"Layer {layer!r} not found in assay {assay_name!r}. Like Seurat, this "
            "reads scale.data by default: run scale_data() first, or pass "
            "layer='data'."
        )
    data, feature_names, layer_cells = _get_expression_layer(assay_obj, layer)

    # Columns are found by the layer's *own* cell names. `seurat.cell_names()`
    # names the metadata rows, not the matrix columns; the two agree on a
    # well-formed object, and when a subset let them drift apart this paired
    # every cell's coordinates with another cell's expression.
    pos = {c: i for i, c in enumerate(layer_cells)}
    keep = [c for c in coords["cell"] if c in pos]
    if len(keep) < 3:
        raise ValueError("Fewer than 3 cells have both coordinates and expression.")
    col_idx = [pos[c] for c in keep]
    xy = coords.set_index("cell").loc[keep, ["x", "y"]].to_numpy(dtype=float)

    if features is not None:
        want = set(features)
        rows = [i for i, f in enumerate(feature_names) if f in want]
        if not rows:
            raise ValueError("None of the requested features are in the assay.")
    else:
        rows = list(range(len(feature_names)))
    genes = [feature_names[i] for i in rows]

    sub = data[rows, :][:, col_idx]
    X = np.asarray(sub.toarray() if sp.issparse(sub) else sub, dtype=float)
    Z = X - X.mean(axis=1, keepdims=True)

    if method == "moransi":
        res = _moransi_table(Z, xy, genes, k=k, weights=weights)
    else:
        res = _markvariogram_table(Z, xy, genes, r_metric=r_metric, bandwidth=bandwidth)

    _write_feature_meta(assay_obj, res)
    return res

composition_test

composition_test(seurat, group_by: str, split_by: str, reference: Optional[str] = None, sample_col: str | None = None) -> DataFrame

Directional abundance test of group_by categories across split_by.

Mirrors the enrichment table an analyst builds by hand: for a two-level split_by (e.g. condition), each group_by category (e.g. cluster) gets a log2(prop_test / prop_reference) and a Fisher exact p (that category vs all others). p-values are BH-adjusted. The overall chi-square p is stored in df.attrs['chisq_p'].

Cells are not replicates. Without sample_col the Fisher test counts every cell as an independent observation. With thousands of cells nearly every difference is then significant, including one that is only donor-to-donor variation, because the unit that replicates is the donor, not the cell. In a simulation with no condition effect at all (12 donors, 6 against 6, each with its own cell-type proportions) at least one cell type came out significant in 35 of 40 runs. Read the p-values as a ranking, or as a test for one sample per condition. For several donors per condition pass sample_col.

Parameters:

• group_by (str) –

  categorical metadata column tested for enrichment (rows).

• split_by (str) –

  metadata column with exactly two levels (the conditions).

• reference (Optional[str], default: None ) –

  which split_by level is the denominator (default: the first sorted level). log2 > 0 ⇒ enriched in the other level.

• sample_col (str | None, default: None ) –

  metadata column naming the sample (donor) each cell came from. The test then has one observation per sample: each group's proportion in each sample is compared between the two levels with a two-sided Mann-Whitney U test, BH-adjusted across groups. Every sample must sit in one level of split_by, and a cell with no sample is left out. prop_<level> is then the mean of the per-sample proportions, odds_ratio and chisq_p are not reported (they are cell-level), and df.attrs['n_samples'] counts the samples per level. With few samples the test cannot reach significance whatever the data: three against three cannot give a p-value below 0.1, and a warning says so.

Returns:

• A DataFrame ordered by ``log2_ratio``, with the columns ``group``, –
• ``n_<ref>``, ``n_<test>``, ``prop_<ref>``, ``prop_<test>``, ``log2_ratio``, –
• ``odds_ratio``, ``p``, ``padj``, ``sig`` and ``enriched_in``. –

Source code in truecell/composition.py

def composition_test(
    seurat,
    group_by: str,
    split_by: str,
    reference: Optional[str] = None,
    sample_col: str | None = None,
) -> pd.DataFrame:
    """Directional abundance test of ``group_by`` categories across ``split_by``.

    Mirrors the enrichment table an analyst builds by hand: for a two-level
    ``split_by`` (e.g. condition), each ``group_by`` category (e.g. cluster) gets
    a ``log2(prop_test / prop_reference)`` and a Fisher exact p (that category vs
    all others). p-values are BH-adjusted. The overall chi-square p is stored in
    ``df.attrs['chisq_p']``.

    **Cells are not replicates.** Without ``sample_col`` the Fisher test counts every
    cell as an independent observation. With thousands of cells nearly every
    difference is then significant, including one that is only donor-to-donor
    variation, because the unit that replicates is the donor, not the cell. In a
    simulation with no condition effect at all (12 donors, 6 against 6, each with its
    own cell-type proportions) at least one cell type came out significant in 35 of 40
    runs. Read the p-values as a ranking, or as a test for one sample per condition.
    For several donors per condition pass ``sample_col``.

    Parameters
    ----------
    group_by   : categorical metadata column tested for enrichment (rows).
    split_by   : metadata column with exactly two levels (the conditions).
    reference  : which ``split_by`` level is the denominator (default: the first
                 sorted level). log2 > 0 ⇒ enriched in the *other* level.
    sample_col : metadata column naming the sample (donor) each cell came from. The
                 test then has one observation per sample: each group's proportion
                 in each sample is compared between the two levels with a two-sided
                 Mann-Whitney U test, BH-adjusted across groups. Every sample must sit
                 in one level of ``split_by``, and a cell with no sample is left out.
                 ``prop_<level>`` is then the mean of the per-sample proportions,
                 ``odds_ratio`` and ``chisq_p`` are not reported (they are cell-level),
                 and ``df.attrs['n_samples']`` counts the samples per level. With few
                 samples the test cannot reach significance whatever the data: three
                 against three cannot give a p-value below 0.1, and a warning says so.

    Returns
    -------
    A DataFrame ordered by ``log2_ratio``, with the columns ``group``,
    ``n_<ref>``, ``n_<test>``, ``prop_<ref>``, ``prop_<test>``, ``log2_ratio``,
    ``odds_ratio``, ``p``, ``padj``, ``sig`` and ``enriched_in``.
    """
    md = seurat.meta_data
    for col in (group_by, split_by) + ((sample_col,) if sample_col is not None else ()):
        if col not in md.columns:
            raise KeyError(f"'{col}' not in meta_data.")
    tab = pd.crosstab(md[group_by].astype(str), md[split_by].astype(str))
    conds = list(tab.columns)
    if len(conds) != 2:
        raise ValueError(
            f"split_by='{split_by}' must have exactly 2 levels, found {conds}."
        )
    ref = reference if reference is not None else conds[0]
    if ref not in conds:
        raise ValueError(f"reference '{ref}' not a level of {split_by}: {conds}.")
    test = [c for c in conds if c != ref][0]
    if sample_col is not None:
        return _across_samples(md, group_by, split_by, sample_col, ref, test)

    n_ref, n_test = tab[ref].sum(), tab[test].sum()
    rows = []
    for grp in tab.index:
        a, b = tab.loc[grp, test], tab.loc[grp, ref]          # this group
        c, d = n_test - a, n_ref - b                          # all other groups
        odds, p = stats.fisher_exact([[a, b], [c, d]])
        prop_test = a / n_test if n_test else np.nan
        prop_ref = b / n_ref if n_ref else np.nan
        with np.errstate(divide="ignore"):
            log2 = np.log2(prop_test / prop_ref) if prop_ref else np.nan
        rows.append({
            "group": grp, f"n_{ref}": b, f"n_{test}": a,
            f"prop_{ref}": prop_ref, f"prop_{test}": prop_test,
            "log2_ratio": log2, "odds_ratio": odds, "p": p,
        })
    df = _finish(pd.DataFrame(rows), ref, test)
    df.attrs["chisq_p"] = float(stats.chi2_contingency(tab.to_numpy())[1])
    df.attrs["reference"] = ref
    df.attrs["test"] = test
    df.attrs["method"] = "fisher_exact"
    return df