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.

Loading

load_xenium

load_xenium(path: Union[str, Path], assay: str = 'Xenium', 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)

By default only Gene Expression features are kept in the assay (matching Seurat's LoadXenium, which routes Negative Control / Blank codewords to separate assays); set keep_controls=True to retain every feature row.

fov_column, if given and present in the cells table, splits the object into one image per FOV; otherwise a single image is created.

Source code in truecell/spatial/loaders.py

def load_xenium(
    path: Union[str, Path],
    assay: str = "Xenium",
    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)

    By default only ``Gene Expression`` features are kept in the assay (matching
    Seurat's ``LoadXenium``, which routes Negative Control / Blank codewords to
    separate assays); set ``keep_controls=True`` to retain every feature row.

    ``fov_column``, if given and present in the cells table, splits the object
    into one image per FOV; otherwise a single image is created.
    """
    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.")
    cdf["cell"] = cdf["cell"].astype(str)

    fov = fov_column if (fov_column and fov_column in cdf.columns) else None
    coords = cdf[["cell", "x", "y"] + ([fov] if fov 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=fov, meta_data=meta)

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_column: str = 'fov', project: str = 'CosMx')

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

expr_file is a cell×gene CSV (rows = cells, first columns cell/fov ids); meta_file carries CenterX_global_px / CenterY_global_px and a FOV column. If names are omitted, *exprMat_file.csv / *metadata_file.csv are auto-detected in path.

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_column: str = "fov",
    project: str = "CosMx",
):
    """Load NanoString CosMx output (exprMat + metadata CSVs) into a Truecell object.

    ``expr_file`` is a cell×gene CSV (rows = cells, first columns cell/fov ids);
    ``meta_file`` carries ``CenterX_global_px`` / ``CenterY_global_px`` and a FOV
    column. If names are omitted, ``*exprMat_file.csv`` / ``*metadata_file.csv``
    are auto-detected in ``path``.
    """
    path = Path(path)
    expr = Path(expr_file) if expr_file else next(iter(path.glob("*exprMat_file.csv")), None)
    meta = Path(meta_file) if meta_file else next(iter(path.glob("*metadata_file.csv")), None)
    if expr is None or meta is None:
        raise FileNotFoundError("Could not locate exprMat_file.csv / metadata_file.csv.")

    edf = pd.read_csv(expr)
    mdf = pd.read_csv(meta)
    id_cols = [c for c in ("fov", "cell_ID", "cell_id", "cell") if c in edf.columns]
    gene_cols = [c for c in edf.columns if c not in id_cols]

    def _cid(df):
        fovc = fov_column if fov_column in df.columns else id_cols[0]
        cidc = next((c for c in ("cell_ID", "cell_id", "cell") if c in df.columns), None)
        return (df[fovc].astype(str) + "_" + df[cidc].astype(str)).to_numpy()

    cell_ids = _cid(edf)
    counts = sp.csc_matrix(edf[gene_cols].to_numpy(dtype=float).T)   # genes × cells

    mdf = mdf.copy()
    mdf["cell"] = _cid(mdf)
    mcoord = mdf.rename(columns={"CenterX_global_px": "x", "CenterY_global_px": "y"})
    coords = mcoord[["cell", "x", "y"] + ([fov_column] if fov_column in mcoord else [])]
    return _build_spatial_object(counts, gene_cols, list(cell_ids), coords,
                                 assay, project,
                                 fov=fov_column if fov_column in coords else None,
                                 meta_data=mdf.set_index("cell"))

load_merscope

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

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

Mirrors Seurat's LoadVizgen. 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)

Blank/control barcodes (Blank-* columns) are dropped by default, matching LoadVizgen; set keep_controls=True to retain them.

fov_column, if present in the metadata, splits the object into one image per FOV; otherwise a single image is created.

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_column: str = "fov",
    project: str = "MERSCOPE",
    keep_controls: bool = False,
):
    """Load a Vizgen MERSCOPE output into a Truecell object with images.

    Mirrors Seurat's ``LoadVizgen``. 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``)

    Blank/control barcodes (``Blank-*`` columns) are dropped by default, matching
    ``LoadVizgen``; set ``keep_controls=True`` to retain them.

    ``fov_column``, if present in the metadata, splits the object into one image
    per FOV; otherwise a single image is created.
    """
    path = Path(path)
    expr = Path(expr_file) if expr_file else _first_existing(
        path, ["cell_by_gene.csv", "cell_by_gene.csv.gz"])
    meta = Path(meta_file) if meta_file else _first_existing(
        path, ["cell_metadata.csv", "cell_metadata.csv.gz"])
    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)

    ecid = _cell_id_column(edf)
    gene_cols = [c for c in edf.columns if c != ecid]
    if not keep_controls:
        gene_cols = [g for g in gene_cols if not str(g).lower().startswith("blank")]
    if not gene_cols:
        raise ValueError(f"No gene columns found in {expr}.")
    cell_ids = edf[ecid].astype(str).to_numpy()
    counts = sp.csc_matrix(edf[gene_cols].to_numpy(dtype=float).T)   # genes × cells

    mdf = mdf.copy()
    mdf["cell"] = mdf[_cell_id_column(mdf)].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.")
    fov = fov_column if fov_column in mcoord.columns else None
    coords = mcoord[["cell", "x", "y"] + ([fov] if fov else [])]
    return _build_spatial_object(counts, gene_cols, list(cell_ids), coords,
                                 assay, project, fov=fov,
                                 meta_data=mdf.set_index("cell"))

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. 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. 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.")
        labels = coords[fov].astype(str).to_numpy()
    else:
        labels = pd.Series(fov).astype(str).to_numpy()
        if len(labels) != len(coords):
            raise ValueError("fov label length must match number of rows in coords.")

    images: dict[str, FOV] = {}
    for name in pd.unique(labels):
        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())

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

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: Optional[str] = None, 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 (Optional[str], default: None ) –

  expression layer (default: the normalized data).

• 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

Run this on log-normalized data (the default). 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: Optional[str] = None,
    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 (default: the normalized ``data``).
    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
    -----
    Run this on log-normalized ``data`` (the default). 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_matrix

    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]
    data, feature_names = _get_expression_matrix(assay_obj, layer)

    pos = {c: i for i, c in enumerate(seurat.cell_names())}
    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) -> 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'].

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.

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,
) -> 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']``.

    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.

    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):
        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]

    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 = pd.DataFrame(rows)
    df["padj"] = _bh(df["p"].to_numpy())
    bins = [-np.inf, 0.001, 0.01, 0.05, np.inf]
    df["sig"] = pd.cut(df["padj"], bins, labels=["***", "**", "*", "ns"])
    df["enriched_in"] = np.where(df["log2_ratio"] > 0, test, ref)
    df = df.sort_values("log2_ratio", ascending=False).reset_index(drop=True)
    df.attrs["chisq_p"] = float(stats.chi2_contingency(tab.to_numpy())[1])
    df.attrs["reference"] = ref
    df.attrs["test"] = test
    return df