diff --git a/CHANGELOG.md b/CHANGELOG.md index c9bb2cf..c0d1810 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,11 +8,28 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] Substrate / microenvironment access for PhysiCell simulations — the main -functional gap called out in the v0.2.0 defect report. Additive; nothing -existing changes behavior. +functional gap called out in the v0.2.0 defect report — plus a new +`microenvironment` analysis module. Additive; nothing existing changes +behavior. ### Added +- **`microenvironment` analysis module.** The previously-empty stub is now + implemented with three composed, domain-framed analyses: + - `identify_niches(data, n_niches, radius=...)` — spatial niches (recurring + local cell-type compositions) by k-means over neighborhood composition + vectors, returning a `NicheResult` with per-cell labels and per-niche + composition profiles. Complementary to the topic-model view in + `spatialtissuepy.lda`. + - `detect_boundaries(data, radius, labels=None)` — cells at the interface + between groups, scored by neighborhood "foreignness". Generalizes the + pairwise `spatial.interface_cells` to any grouping; pass niche labels to + find niche interfaces. Returns a `BoundaryResult`. + - `spatial_gradient` / `substrate_gradient` / `density_gradient` — spatial + gradients of a scalar field (substrate concentration or local cell density) + by local linear regression, returning a `GradientField` with `magnitude` + and `direction`. Validated against a finite-difference reference on the + regular voxel grid (correlation > 0.97 on interior points). - **Environmental substrate sampling on `PhysiCellTimeStep`.** The microenvironment `.mat` (oxygen and other diffusible fields) is now reachable through the high-level API, having previously been parseable only by calling diff --git a/spatialtissuepy/microenvironment/__init__.py b/spatialtissuepy/microenvironment/__init__.py index 44a0450..3576d9a 100644 --- a/spatialtissuepy/microenvironment/__init__.py +++ b/spatialtissuepy/microenvironment/__init__.py @@ -1,15 +1,37 @@ """ -Tumor microenvironment analysis module. +Tumor microenvironment analysis. -.. note:: +Higher-level, domain-framed analyses of tissue microenvironment structure, +composed from the spatial, statistics, and PhysiCell-substrate machinery: - This module is a planned feature for the v0.3.0 release and is not yet - implemented. It will provide niche identification, boundary detection, and - spatial gradient analysis. In the meantime, much of this functionality can - be assembled from the :mod:`spatialtissuepy.spatial` and - :mod:`spatialtissuepy.statistics` modules (e.g. neighborhood composition, - Getis-Ord Gi* hotspots, and boundary cell detection). +- **Niche identification** (:func:`identify_niches`): recurring local cell-type + compositions, found by clustering neighborhood composition vectors. +- **Boundary detection** (:func:`detect_boundaries`): cells at the interface + between cell types or niches. +- **Gradient analysis** (:func:`spatial_gradient`, :func:`substrate_gradient`, + :func:`density_gradient`): spatial gradients of substrate concentration or + cell density. """ -# Planned for v0.3.0 — see module docstring. -__all__ = [] +from .boundaries import BoundaryResult, detect_boundaries +from .gradients import ( + GradientField, + density_gradient, + spatial_gradient, + substrate_gradient, +) +from .niches import NicheResult, identify_niches + +__all__ = [ + # Niches + 'identify_niches', + 'NicheResult', + # Boundaries + 'detect_boundaries', + 'BoundaryResult', + # Gradients + 'spatial_gradient', + 'substrate_gradient', + 'density_gradient', + 'GradientField', +] diff --git a/spatialtissuepy/microenvironment/boundaries.py b/spatialtissuepy/microenvironment/boundaries.py new file mode 100644 index 0000000..1f5ee5a --- /dev/null +++ b/spatialtissuepy/microenvironment/boundaries.py @@ -0,0 +1,132 @@ +""" +Boundary and interface detection. + +A boundary cell sits at the interface between regions of different identity -- +different cell types, or different niches. Each cell's *foreignness* is the +fraction of its spatial neighbors whose label differs from its own; cells above +a threshold are boundary cells. Passing niche labels (from +:func:`~spatialtissuepy.microenvironment.niches.identify_niches`) detects niche +interfaces; the default uses cell type, generalizing the pairwise +:func:`~spatialtissuepy.spatial.interface_cells` to any grouping. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +from ..core import SpatialTissueData +from ..spatial import compute_neighborhoods + + +@dataclass +class BoundaryResult: + """ + Result of :func:`detect_boundaries`. + + Attributes + ---------- + is_boundary : np.ndarray + ``(n_cells,)`` boolean mask; ``True`` for boundary cells. + foreignness : np.ndarray + ``(n_cells,)`` fraction of each cell's neighbors whose label differs + from its own. ``0`` for cells with too few neighbors. + labels : np.ndarray + ``(n_cells,)`` grouping labels used for the comparison. + radius : float + Neighborhood radius used. + """ + + is_boundary: np.ndarray + foreignness: np.ndarray + labels: np.ndarray + radius: float + + @property + def boundary_indices(self) -> np.ndarray: + """Indices of the boundary cells.""" + return np.where(self.is_boundary)[0] + + def boundary_fraction(self) -> float: + """Fraction of all cells that are boundary cells.""" + if self.is_boundary.size == 0: + return 0.0 + return float(np.mean(self.is_boundary)) + + +def detect_boundaries( + data: SpatialTissueData, + radius: float, + labels: Optional[np.ndarray] = None, + threshold: float = 0.0, + min_neighbors: int = 1, +) -> BoundaryResult: + """ + Detect cells at the interface between groups. + + For each cell, ``foreignness`` is the fraction of neighbors within + ``radius`` whose group label differs from the cell's own. A cell is a + boundary cell when its foreignness exceeds ``threshold`` and it has at least + ``min_neighbors`` neighbors. + + Parameters + ---------- + data : SpatialTissueData + Spatial tissue data. + radius : float + Neighborhood radius. + labels : np.ndarray, optional + ``(n_cells,)`` group label per cell. Defaults to the cell types. Pass + niche labels to detect niche interfaces. + threshold : float, default 0.0 + Minimum foreignness (exclusive) to call a cell a boundary. The default + of ``0.0`` flags any cell with at least one differing neighbor; raise it + to require a more mixed neighborhood. + min_neighbors : int, default 1 + Minimum neighbor count; cells with fewer are never boundary cells and + get foreignness ``0``. + + Returns + ------- + BoundaryResult + + Raises + ------ + ValueError + If ``labels`` is given but its length does not match ``n_cells``. + """ + if labels is None: + labels = np.asarray(data.cell_types) + else: + labels = np.asarray(labels) + if labels.shape[0] != data.n_cells: + raise ValueError( + f"labels length {labels.shape[0]} != n_cells {data.n_cells}" + ) + + neighborhoods = compute_neighborhoods( + data, method='radius', radius=radius, include_self=False + ) + + n = data.n_cells + foreignness = np.zeros(n) + is_boundary = np.zeros(n, dtype=bool) + + for i, neighbors in enumerate(neighborhoods): + neighbors = np.asarray(neighbors) + if neighbors.size < min_neighbors: + continue + differing = np.count_nonzero(labels[neighbors] != labels[i]) + frac = differing / neighbors.size + foreignness[i] = frac + if frac > threshold: + is_boundary[i] = True + + return BoundaryResult( + is_boundary=is_boundary, + foreignness=foreignness, + labels=labels, + radius=radius, + ) diff --git a/spatialtissuepy/microenvironment/gradients.py b/spatialtissuepy/microenvironment/gradients.py new file mode 100644 index 0000000..b50ded5 --- /dev/null +++ b/spatialtissuepy/microenvironment/gradients.py @@ -0,0 +1,254 @@ +""" +Spatial gradient analysis. + +Estimates the spatial gradient of a scalar field -- a substrate concentration, +or local cell density -- by fitting a local linear model to the values at each +point's nearest neighbors. The gradient vector points in the direction of +steepest increase, and its magnitude is the rate of change per unit distance. +This is grid-agnostic: it works on the regular voxel mesh of a microenvironment +field and on scattered cell positions alike. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Optional + +import numpy as np + +from ..core import SpatialTissueData +from ..spatial import compute_neighborhoods + +if TYPE_CHECKING: + from ..synthetic.physicell import PhysiCellTimeStep + + +@dataclass +class GradientField: + """ + Result of a gradient estimation. + + Attributes + ---------- + points : np.ndarray + ``(n_points, n_dims)`` coordinates where the gradient was evaluated. + gradients : np.ndarray + ``(n_points, n_dims)`` gradient vector at each point. + values : np.ndarray + ``(n_points,)`` field value sampled at each point. + """ + + points: np.ndarray + gradients: np.ndarray + values: np.ndarray + + @property + def magnitude(self) -> np.ndarray: + """``(n_points,)`` gradient magnitude (steepness).""" + return np.linalg.norm(self.gradients, axis=1) + + @property + def direction(self) -> np.ndarray: + """ + ``(n_points, n_dims)`` unit vectors along each gradient. + + Points whose gradient magnitude is negligible relative to the steepest + gradient in the field -- flat regions, and the round-off noise a local + fit produces on a constant field -- yield zero vectors rather than + normalized noise or NaNs. + """ + mag = self.magnitude + scale = mag.max() if mag.size else 0.0 + if scale <= 1e-12: + # The whole field is flat; no meaningful directions. + return np.zeros_like(self.gradients) + tol = scale * 1e-6 + safe = np.where(mag > tol, mag, 1.0) + unit = self.gradients / safe[:, None] + unit[mag <= tol] = 0.0 + return unit + + +def spatial_gradient( + positions: np.ndarray, + values: np.ndarray, + query_points: Optional[np.ndarray] = None, + k: int = 12, +) -> GradientField: + """ + Estimate the gradient of a scalar field by local linear regression. + + Around each query point, the ``k`` nearest source points are found and a + linear model ``value ~ a + g . (x - q)`` is fit by least squares; the slope + ``g`` is the gradient estimate. Points are centered on the query for + numerical stability. + + Parameters + ---------- + positions : np.ndarray + ``(n, n_dims)`` coordinates where the field is known. + values : np.ndarray + ``(n,)`` field values at ``positions``. + query_points : np.ndarray, optional + ``(m, n_dims)`` points at which to evaluate the gradient. Defaults to + ``positions``. The reported ``values`` are then nearest-source samples. + k : int, default 12 + Number of neighbors for each local fit. Raised to at least ``n_dims + 1`` + (a fit needs that many points to be well-posed) and capped at the number + of source points. + + Returns + ------- + GradientField + Points whose local fit is under-determined -- too few points, or + collinear/degenerate neighbors -- get a ``NaN`` gradient rather than a + misleading minimum-norm estimate. + + Raises + ------ + ValueError + If shapes are inconsistent or there are no source points. + """ + from scipy.spatial import cKDTree + + positions = np.asarray(positions, dtype=float) + values = np.asarray(values, dtype=float) + if positions.ndim != 2: + raise ValueError("positions must be (n, n_dims)") + n, ndim = positions.shape + if n == 0: + raise ValueError("no source points") + if values.shape != (n,): + raise ValueError(f"values must be ({n},), got {values.shape}") + + if query_points is None: + query = positions + sampled = values + else: + query = np.asarray(query_points, dtype=float) + if query.ndim != 2 or query.shape[1] != ndim: + raise ValueError(f"query_points must be (m, {ndim})") + + # A well-posed local linear fit needs at least ndim + 1 points; ask for that + # many even if the caller passed a smaller k. Still capped at n. + k_eff = min(max(k, ndim + 1), n) + tree = cKDTree(positions) + _, idx = tree.query(query, k=k_eff) + idx = np.atleast_2d(idx.T).T # ensure (m, k_eff) even when k_eff == 1 + + if query_points is not None: + # Nearest-source value at each query point (consistent with + # substrate_at's nearest-voxel convention). + sampled = values[idx[:, 0]] + + gradients = np.zeros((query.shape[0], ndim)) + for i in range(query.shape[0]): + neigh = idx[i] + # Local design matrix [1, dx, dy, ...] centered on the query point. + offsets = positions[neigh] - query[i] + design = np.column_stack([np.ones(len(neigh)), offsets]) + coeffs, _, rank, _ = np.linalg.lstsq(design, values[neigh], rcond=None) + if rank < ndim + 1: + # Under-determined (too few points, or collinear/degenerate + # neighbors): the gradient is not identifiable. Flag rather than + # return the misleading minimum-norm solution. + gradients[i] = np.nan + else: + gradients[i] = coeffs[1:] + + return GradientField(points=query, gradients=gradients, values=sampled) + + +def substrate_gradient( + timestep: PhysiCellTimeStep, + name: str, + query_points: Optional[np.ndarray] = None, + k: int = 12, +) -> GradientField: + """ + Gradient of a substrate concentration field. + + Evaluated at the voxel centers by default. Constant spatial axes (e.g. ``z`` + in a 2-D simulation) are dropped, so a 2-D field yields 2-D gradients. + + Parameters + ---------- + timestep : PhysiCellTimeStep + A time step exposing :attr:`voxel_positions` and :attr:`substrates`. + name : str + Substrate name. + query_points : np.ndarray, optional + Points at which to evaluate, with one column per non-constant axis. + Defaults to the voxel centers. + k : int, default 12 + Neighbors per local fit. + + Returns + ------- + GradientField + + Raises + ------ + ValueError + If the substrate is unknown or there is no microenvironment field. + """ + voxels = timestep.voxel_positions + substrates = timestep.substrates + if name not in substrates: + available = ', '.join(substrates) or '(none)' + raise ValueError(f"Unknown substrate {name!r}. Available: {available}") + if voxels.shape[0] == 0: + raise ValueError("No microenvironment voxels available") + + # Keep only axes that vary, so a planar (2-D) field is treated as 2-D. + active = np.where(np.ptp(voxels, axis=0) > 0)[0] + if active.size == 0: + active = np.array([0]) + positions = voxels[:, active] + + return spatial_gradient(positions, substrates[name], query_points, k=k) + + +def density_gradient( + data: SpatialTissueData, + radius: float, + query_points: Optional[np.ndarray] = None, + k: int = 12, +) -> GradientField: + """ + Gradient of local cell density. + + Local density is the number of neighbors within ``radius`` divided by the + neighborhood area (2-D) or volume (3-D); its gradient points toward denser + regions. + + Parameters + ---------- + data : SpatialTissueData + Spatial tissue data. + radius : float + Radius defining local density. + query_points : np.ndarray, optional + Points at which to evaluate the gradient. Defaults to the cell + positions. + k : int, default 12 + Neighbors per local fit. + + Returns + ------- + GradientField + """ + coords = data.coordinates + ndim = coords.shape[1] + if ndim == 2: + volume = np.pi * radius ** 2 + else: + volume = (4.0 / 3.0) * np.pi * radius ** 3 + + neighborhoods = compute_neighborhoods( + data, method='radius', radius=radius, include_self=False + ) + counts = np.array([len(np.asarray(n)) for n in neighborhoods], dtype=float) + density = counts / volume + + return spatial_gradient(coords, density, query_points, k=k) diff --git a/spatialtissuepy/microenvironment/niches.py b/spatialtissuepy/microenvironment/niches.py new file mode 100644 index 0000000..1a0e261 --- /dev/null +++ b/spatialtissuepy/microenvironment/niches.py @@ -0,0 +1,175 @@ +""" +Spatial niche identification. + +A *niche* (also called a cellular neighborhood) is a recurring local +composition of cell types. Cells are described by the mix of types in their +spatial neighborhood, and those composition vectors are clustered; each cluster +is a niche. This is the composition-clustering approach popularized for +multiplexed imaging (e.g. Schürch/Nolan cellular neighborhoods), and is +complementary to the topic-model view in :mod:`spatialtissuepy.lda`. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Dict, List, Optional, Tuple + +import numpy as np +import pandas as pd +from sklearn.cluster import KMeans + +from ..core import SpatialTissueData +from ..spatial import neighborhood_composition + + +@dataclass +class NicheResult: + """ + Result of :func:`identify_niches`. + + Attributes + ---------- + labels : np.ndarray + ``(n_cells,)`` integer niche id for each cell. + n_niches : int + Number of niches. + cell_types : np.ndarray + ``(n_types,)`` cell-type names, giving the column order of + :attr:`profiles` and of the composition vectors. + profiles : np.ndarray + ``(n_niches, n_types)`` mean neighborhood composition of each niche. + radius : float or None + Radius used, if the radius method was chosen. + k : int or None + Neighbor count used, if the knn method was chosen. + """ + + labels: np.ndarray + n_niches: int + cell_types: np.ndarray + profiles: np.ndarray + radius: Optional[float] = None + k: Optional[int] = None + + def niche_sizes(self) -> Dict[int, int]: + """Number of cells assigned to each niche.""" + ids, counts = np.unique(self.labels, return_counts=True) + return {int(i): int(c) for i, c in zip(ids, counts)} + + def dominant_types(self, top: int = 3) -> Dict[int, List[Tuple[str, float]]]: + """ + The most abundant cell types in each niche's mean composition. + + Parameters + ---------- + top : int, default 3 + How many types to report per niche. + + Returns + ------- + dict + ``{niche_id: [(type_name, mean_fraction), ...]}``, largest first. + """ + result: Dict[int, List[Tuple[str, float]]] = {} + for niche in range(self.n_niches): + order = np.argsort(self.profiles[niche])[::-1][:top] + result[niche] = [ + (str(self.cell_types[j]), float(self.profiles[niche, j])) + for j in order + ] + return result + + def profiles_dataframe(self) -> pd.DataFrame: + """Niche composition profiles as a DataFrame (niches x types).""" + return pd.DataFrame( + self.profiles, + index=[f'niche_{i}' for i in range(self.n_niches)], + columns=list(self.cell_types), + ) + + +def identify_niches( + data: SpatialTissueData, + n_niches: int, + radius: Optional[float] = None, + k: Optional[int] = None, + method: str = 'radius', + include_self: bool = True, + random_state: Optional[int] = None, +) -> NicheResult: + """ + Identify spatial niches by clustering neighborhood composition. + + Each cell's local neighborhood composition (the fraction of each cell type + within its neighborhood) is computed, and cells are grouped into ``n_niches`` + clusters by k-means. A niche is therefore a characteristic local cell-type + mixture, wherever in the tissue it recurs. + + Parameters + ---------- + data : SpatialTissueData + Spatial tissue data. + n_niches : int + Number of niches (k-means clusters). Must be at least 1 and at most the + number of cells. + radius : float, optional + Neighborhood radius (required when ``method='radius'``). + k : int, optional + Number of neighbors (required when ``method='knn'``). + method : str, default 'radius' + Neighborhood definition, ``'radius'`` or ``'knn'``. + include_self : bool, default True + Whether a cell's own type counts toward its neighborhood composition. + random_state : int, optional + Seed for k-means, for reproducibility. + + Returns + ------- + NicheResult + Per-cell niche labels and per-niche composition profiles. + + Raises + ------ + ValueError + If ``n_niches`` is out of range, or the neighborhood parameter for the + chosen method is missing. + """ + if method == 'radius' and radius is None: + raise ValueError("radius is required when method='radius'") + if method == 'knn' and k is None: + raise ValueError("k is required when method='knn'") + if not 1 <= n_niches <= data.n_cells: + raise ValueError( + f"n_niches must be in [1, n_cells={data.n_cells}], got {n_niches}" + ) + + composition = neighborhood_composition( + data, + method=method, + radius=radius, + k=k, + include_self=include_self, + normalize=True, + ) + + kmeans = KMeans(n_clusters=n_niches, random_state=random_state, n_init=10) + labels = kmeans.fit_predict(composition) + + # Per-niche mean composition (from the data, not the k-means centroids, so + # the profile is a genuine average of member cells even if a niche is empty + # of some type). + cell_types = data.cell_types_unique + profiles = np.zeros((n_niches, len(cell_types))) + for niche in range(n_niches): + members = labels == niche + if np.any(members): + profiles[niche] = composition[members].mean(axis=0) + + return NicheResult( + labels=labels.astype(int), + n_niches=n_niches, + cell_types=cell_types, + profiles=profiles, + radius=radius, + k=k, + ) diff --git a/spatialtissuepy/synthetic/physicell/reader.py b/spatialtissuepy/synthetic/physicell/reader.py index 839b3b4..e58d8ae 100644 --- a/spatialtissuepy/synthetic/physicell/reader.py +++ b/spatialtissuepy/synthetic/physicell/reader.py @@ -231,15 +231,17 @@ def substrate_at( name : str Substrate name; must be one of :attr:`substrate_names`. x, y : float or array-like - Query coordinates. Scalars or equal-length arrays. + Query coordinates. Scalars, arrays, or any broadcastable + combination (e.g. scalar ``x`` with an array ``y``). z : float or array-like, optional Query z. If given, nearest-voxel search is 3-D; otherwise it uses the x-y plane (appropriate for 2-D simulations). Returns ------- - np.ndarray - Concentration at each query point, matching the input shape. + float or np.ndarray + Concentration at each query point. A scalar in yields a scalar out; + otherwise the result has the broadcast shape of the inputs. """ me = self._load_microenvironment() concentrations = me['concentrations'] @@ -258,18 +260,26 @@ def substrate_at( x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) if z is None: - query = np.column_stack([np.ravel(x), np.ravel(y)]) + x, y = np.broadcast_arrays(x, y) + out_shape = x.shape + query = np.column_stack([x.ravel(), y.ravel()]) ndim = 2 else: z = np.asarray(z, dtype=float) - query = np.column_stack([np.ravel(x), np.ravel(y), np.ravel(z)]) + x, y, z = np.broadcast_arrays(x, y, z) + out_shape = x.shape + query = np.column_stack([x.ravel(), y.ravel(), z.ravel()]) ndim = 3 tree = self._voxel_kdtree(ndim) _, idx = tree.query(query) result = concentrations[name][idx] - return result if x.ndim else result[0] + # Return shape mirrors the query: scalar in -> scalar out, and an + # (m,) or (h, w) grid of points in -> the same shape out. + if out_shape == (): + return result[0] + return result.reshape(out_shape) def _voxel_kdtree(self, ndim: int): """Build (and cache) a KD-tree over voxel centers for `ndim` axes.""" @@ -320,16 +330,41 @@ def internalized_substrates( (the model did not write the ``internalized_total_substrates`` field, or the XML labels were unavailable). """ + from .parser import expand_cell_labels + data = self._load_cell_data() columns = data.get('columns') or {} names = self.substrate_names - field_keys = [f'internalized_total_substrates_{i}' for i in range(len(names))] - if not names or not all(k in columns for k in field_keys): + # Resolve the field's column names the same way the parser named them. + # A size-1 vector is stored under the bare name, a size-3 vector under + # _x/_y/_z, others under _0.._{n-1}; hardcoding one convention breaks + # for single- and triple-substrate models. + field_label = None + try: + custom_labels = self._load_metadata().extra.get('custom_labels') or {} + except (OSError, ET.ParseError): + custom_labels = {} + for start, (label_name, size) in custom_labels.items(): + if label_name == 'internalized_total_substrates': + field_label = (start, size) + break + + field_keys = None + if names and field_label is not None: + start, size = field_label + if size == len(names): + expanded = expand_cell_labels( + {start: ('internalized_total_substrates', size)} + ) + field_keys = [expanded[start + off] for off in range(size)] + + if field_keys is None or not all(k in columns for k in field_keys): raise ValueError( "This frame does not record internalized substrates. It " "requires the PhysiCell 'internalized_total_substrates' field " - "and readable XML ; neither is optional here." + "(sized to the substrate count) and readable XML ; " + "neither is optional here." ) if include_dead_cells is None: diff --git a/tests/test_microenvironment.py b/tests/test_microenvironment.py new file mode 100644 index 0000000..ce2b84a --- /dev/null +++ b/tests/test_microenvironment.py @@ -0,0 +1,260 @@ +""" +Tests for the microenvironment analysis module (v0.4.0). + +Covers niche identification, boundary detection, and gradient estimation. +Where possible, correctness is checked against a construction with a known +answer (two separated composition regions; a linear field with a known +gradient) rather than a golden value. +""" + +from pathlib import Path + +import numpy as np +import pytest + +from spatialtissuepy.core import SpatialTissueData +from spatialtissuepy.microenvironment import ( + BoundaryResult, + GradientField, + NicheResult, + density_gradient, + detect_boundaries, + identify_niches, + spatial_gradient, + substrate_gradient, +) + + +@pytest.fixture +def two_region_data(): + """Left half tumor-dominant, right half immune-dominant.""" + rng = np.random.default_rng(0) + left = np.column_stack([rng.uniform(0, 50, 100), rng.uniform(0, 100, 100)]) + right = np.column_stack([rng.uniform(50, 100, 100), rng.uniform(0, 100, 100)]) + coords = np.vstack([left, right]) + types = np.array( + ['Tumor'] * 85 + ['Immune'] * 15 + ['Immune'] * 85 + ['Tumor'] * 15 + ) + return SpatialTissueData(coordinates=coords, cell_types=types) + + +@pytest.fixture +def example_timestep(): + path = ( + Path(__file__).parent.parent + / 'examples' / 'sample_data' / 'example_physicell_sim' + ) + if not path.exists(): + pytest.skip("Example PhysiCell data not found") + from spatialtissuepy.synthetic.physicell import read_physicell_timestep + xmls = sorted(path.glob('output*.xml')) + return read_physicell_timestep(xmls[len(xmls) // 2]) + + +class TestNiches: + def test_labels_shape_and_range(self, two_region_data): + res = identify_niches(two_region_data, n_niches=3, radius=20.0, + random_state=0) + assert isinstance(res, NicheResult) + assert res.labels.shape == (two_region_data.n_cells,) + assert set(np.unique(res.labels)).issubset(set(range(3))) + + def test_separates_two_regions(self, two_region_data): + """With two composition regions, 2 niches should split them.""" + res = identify_niches(two_region_data, n_niches=2, radius=20.0, + random_state=0) + # Each niche's dominant type should differ between the two niches. + dom = res.dominant_types(top=1) + assert dom[0][0][0] != dom[1][0][0] + + def test_profiles_are_distributions(self, two_region_data): + res = identify_niches(two_region_data, n_niches=3, radius=20.0, + random_state=0) + # Non-empty niche profiles sum to ~1 (composition proportions). + for niche, size in res.niche_sizes().items(): + if size > 0: + assert np.isclose(res.profiles[niche].sum(), 1.0, atol=1e-6) + + def test_profiles_dataframe(self, two_region_data): + res = identify_niches(two_region_data, n_niches=2, radius=20.0, + random_state=0) + df = res.profiles_dataframe() + assert list(df.columns) == list(two_region_data.cell_types_unique) + assert len(df) == 2 + + def test_reproducible(self, two_region_data): + a = identify_niches(two_region_data, n_niches=3, radius=20.0, + random_state=42) + b = identify_niches(two_region_data, n_niches=3, radius=20.0, + random_state=42) + np.testing.assert_array_equal(a.labels, b.labels) + + def test_requires_radius(self, two_region_data): + with pytest.raises(ValueError, match="radius is required"): + identify_niches(two_region_data, n_niches=2, method='radius') + + def test_rejects_bad_n_niches(self, two_region_data): + with pytest.raises(ValueError, match="n_niches"): + identify_niches(two_region_data, n_niches=0, radius=20.0) + with pytest.raises(ValueError, match="n_niches"): + identify_niches(two_region_data, n_niches=10 ** 6, radius=20.0) + + +class TestBoundaries: + def test_interface_has_higher_foreignness(self, two_region_data): + res = detect_boundaries(two_region_data, radius=15.0) + assert isinstance(res, BoundaryResult) + # Cells near the x=50 interface should be more "foreign" on average + # than cells deep in a region. + coords = two_region_data.coordinates + near = np.abs(coords[:, 0] - 50) < 10 + deep = np.abs(coords[:, 0] - 50) > 35 + assert res.foreignness[near].mean() > res.foreignness[deep].mean() + + def test_foreignness_bounded(self, two_region_data): + res = detect_boundaries(two_region_data, radius=15.0) + assert np.all(res.foreignness >= 0) and np.all(res.foreignness <= 1) + + def test_boundary_indices_and_fraction(self, two_region_data): + res = detect_boundaries(two_region_data, radius=15.0) + assert res.boundary_indices.shape[0] == res.is_boundary.sum() + assert 0.0 <= res.boundary_fraction() <= 1.0 + + def test_accepts_custom_labels(self, two_region_data): + niche = identify_niches(two_region_data, n_niches=2, radius=20.0, + random_state=0) + res = detect_boundaries(two_region_data, radius=15.0, + labels=niche.labels) + np.testing.assert_array_equal(res.labels, niche.labels) + + def test_rejects_mismatched_labels(self, two_region_data): + with pytest.raises(ValueError, match="labels length"): + detect_boundaries(two_region_data, radius=15.0, + labels=np.array([0, 1, 2])) + + def test_threshold_reduces_boundaries(self, two_region_data): + low = detect_boundaries(two_region_data, radius=15.0, threshold=0.0) + high = detect_boundaries(two_region_data, radius=15.0, threshold=0.5) + assert high.is_boundary.sum() <= low.is_boundary.sum() + + def test_min_neighbors_excludes_isolated(self): + # Two far-apart cells: neither has neighbors within radius. + data = SpatialTissueData( + coordinates=np.array([[0.0, 0.0], [1000.0, 1000.0]]), + cell_types=np.array(['A', 'B']), + ) + res = detect_boundaries(data, radius=5.0, min_neighbors=1) + assert not res.is_boundary.any() + assert np.all(res.foreignness == 0) + + +class TestSpatialGradient: + def test_linear_field_recovers_exact_gradient(self): + """A local linear fit is exact for a globally linear field.""" + rng = np.random.default_rng(0) + pos = rng.uniform(0, 100, (300, 2)) + values = 3.0 * pos[:, 0] + 2.0 * pos[:, 1] + 5.0 # grad = (3, 2) + gf = spatial_gradient(pos, values, k=8) + # Interior points recover (3, 2) closely; use median for robustness. + np.testing.assert_allclose(np.median(gf.gradients, axis=0), [3, 2], + atol=1e-6) + + def test_magnitude_and_direction(self): + pos = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) + values = pos[:, 0] # gradient (1, 0) + gf = spatial_gradient(pos, values, k=4) + assert np.allclose(gf.magnitude, 1.0) + # direction unit vectors + norms = np.linalg.norm(gf.direction, axis=1) + assert np.allclose(norms, 1.0) + + def test_zero_field_direction_is_zero(self): + pos = np.random.default_rng(0).uniform(0, 10, (20, 2)) + gf = spatial_gradient(pos, np.ones(20), k=5) + assert np.allclose(gf.gradients, 0.0) + assert np.allclose(gf.direction, 0.0) # no NaNs + + def test_query_points(self): + pos = np.random.default_rng(1).uniform(0, 10, (50, 2)) + values = pos[:, 0] * 2 + q = np.array([[5.0, 5.0], [2.0, 8.0]]) + gf = spatial_gradient(pos, values, query_points=q, k=8) + assert gf.points.shape == (2, 2) + assert gf.gradients.shape == (2, 2) + np.testing.assert_allclose(gf.gradients, [[2, 0], [2, 0]], atol=1e-6) + + def test_shape_validation(self): + with pytest.raises(ValueError): + spatial_gradient(np.zeros((0, 2)), np.zeros(0)) + with pytest.raises(ValueError): + spatial_gradient(np.zeros((5, 2)), np.zeros(4)) + + def test_too_few_points_is_nan_not_wrong(self): + """2 points in 2-D can't determine a 2-D gradient -> NaN, not (g, 0).""" + pos = np.array([[0.0, 0.0], [1.0, 0.0]]) + values = np.array([0.0, 3.0]) # true gradient underdetermined in y + gf = spatial_gradient(pos, values, k=2) + assert np.isnan(gf.gradients).all() + + def test_collinear_neighbors_are_nan(self): + """Collinear points give a rank-deficient fit -> NaN, not a silent + wrong minimum-norm gradient.""" + pos = np.column_stack([np.arange(5.0), np.zeros(5)]) # all on y=0 + values = 3.0 * pos[:, 0] + gf = spatial_gradient(pos, values, k=5) + assert np.isnan(gf.gradients).all() + + +class TestSubstrateGradient: + def test_planar_field_gives_2d_gradient(self, example_timestep): + gf = substrate_gradient(example_timestep, 'oxygen') + assert gf.gradients.shape[1] == 2 # z is constant, dropped + assert gf.points.shape[0] == example_timestep.voxel_positions.shape[0] + + def test_agrees_with_finite_difference(self, example_timestep): + """Local-fit gradient correlates strongly with a grid finite diff.""" + vox = example_timestep.voxel_positions[:, :2] + ox = example_timestep.substrates['oxygen'] + xs, ys = np.unique(vox[:, 0]), np.unique(vox[:, 1]) + if len(xs) * len(ys) != len(vox): + pytest.skip("Non-grid microenvironment") + ix = np.searchsorted(xs, vox[:, 0]) + iy = np.searchsorted(ys, vox[:, 1]) + grid = np.full((len(xs), len(ys)), np.nan) + grid[ix, iy] = ox + gx_true, gy_true = np.gradient(grid, xs, ys, edge_order=2) + + gf = substrate_gradient(example_timestep, 'oxygen') + gxm = np.full_like(grid, np.nan) + gxm[ix, iy] = gf.gradients[:, 0] + interior = np.zeros_like(grid, dtype=bool) + interior[2:-2, 2:-2] = True + corr = np.corrcoef(gxm[interior], gx_true[interior])[0, 1] + assert corr > 0.9 + + def test_unknown_substrate_raises(self, example_timestep): + with pytest.raises(ValueError, match="Unknown substrate"): + substrate_gradient(example_timestep, 'nope') + + +class TestDensityGradient: + def test_shape(self, two_region_data): + gf = density_gradient(two_region_data, radius=20.0) + assert isinstance(gf, GradientField) + assert gf.gradients.shape == (two_region_data.n_cells, 2) + + def test_points_toward_denser_region(self): + """A cluster next to sparse space: density gradient points inward.""" + rng = np.random.default_rng(0) + dense = rng.uniform(0, 20, (200, 2)) + sparse = rng.uniform(60, 100, (10, 2)) + coords = np.vstack([dense, sparse]) + data = SpatialTissueData( + coordinates=coords, + cell_types=np.array(['A'] * len(coords)), + ) + gf = density_gradient(data, radius=15.0) + # At a dense-cluster edge cell (near x=20), gradient x-component should + # be negative (pointing back into the dense cluster at lower x). + edge = np.argmin(np.abs(coords[:200, 0] - 20)) + assert gf.gradients[edge, 0] <= 0 diff --git a/tests/test_physicell_microenvironment.py b/tests/test_physicell_microenvironment.py index a60abc1..7ce84b6 100644 --- a/tests/test_physicell_microenvironment.py +++ b/tests/test_physicell_microenvironment.py @@ -60,6 +60,44 @@ def _injected_timestep(voxels, concentrations): return ts +def _build_substrate_frame(dir_path, n_subs, start=6, n_cells=8): + """ + Write a minimal PhysiCell-like frame with `n_subs` substrates and an + internalized_total_substrates field of size n_subs. Value in substrate i, + cell j is (i+1)*100 + j, so the mapping can be checked exactly. + + Returns the substrate names. + """ + subs = [f'sub{i}' for i in range(n_subs)] + n_rows = start + n_subs # declared variable count == matrix rows + matrix = np.zeros((n_rows, n_cells)) + matrix[0] = np.arange(n_cells) + matrix[1] = np.linspace(0, 70, n_cells) + matrix[2] = np.linspace(0, 70, n_cells) + for i in range(n_subs): + matrix[start + i] = (i + 1) * 100 + np.arange(n_cells) + savemat(str(dir_path / 'output00000000_cells.mat'), {'cells': matrix}) + + var_xml = ''.join( + f'' for i, s in enumerate(subs) + ) + labels = [ + ('ID', 0, 1), ('position', 1, 3), ('cell_type', 5, 1), + ('internalized_total_substrates', start, n_subs), + ] + label_xml = ''.join( + f'' for nm, idx, sz in labels + ) + (dir_path / 'output00000000.xml').write_text( + '' + f'{var_xml}' + '' + f'{label_xml}' + '' + ) + return subs + + class TestSubstrateFields: def test_substrate_names(self, frame_xml): ts = read_physicell_timestep(frame_xml) @@ -131,6 +169,24 @@ def test_3d_query(self): assert ts.substrate_at('o2', 0.0, 0.0, 1.0) == 5.0 assert ts.substrate_at('o2', 0.0, 0.0, 9.0) == 9.0 + def test_output_shape_mirrors_input(self): + # 4 voxels on a line; query with a 2-D grid of points. + voxels = [[0, 0, 0], [10, 0, 0], [20, 0, 0], [30, 0, 0]] + ts = _injected_timestep(voxels, {'o2': [1.0, 2.0, 3.0, 4.0]}) + + grid = ts.substrate_at('o2', np.array([[0.0, 10.0], [20.0, 30.0]]), + np.array([[0.0, 0.0], [0.0, 0.0]])) + assert grid.shape == (2, 2) + np.testing.assert_array_equal(grid, [[1.0, 2.0], [3.0, 4.0]]) + + def test_broadcasts_scalar_and_array(self): + voxels = [[0, 0, 0], [0, 10, 0]] + ts = _injected_timestep(voxels, {'o2': [1.0, 2.0]}) + # scalar x, array y -> broadcast to array + out = ts.substrate_at('o2', 0.0, np.array([0.0, 10.0])) + assert out.shape == (2,) + np.testing.assert_array_equal(out, [1.0, 2.0]) + class TestInternalizedSubstrates: def test_columns_match_raw_matrix(self, frame_xml, example_physicell_dir): @@ -173,6 +229,23 @@ def test_raises_without_field(self, tmp_path): with pytest.raises(ValueError, match="internalized"): ts.internalized_substrates() + @pytest.mark.parametrize('n_subs', [1, 3, 5]) + def test_maps_field_for_any_substrate_count(self, tmp_path, n_subs): + """The field's columns are named by the label's own size convention: + bare name for size 1, _x/_y/_z for size 3, _0.._n otherwise. Hardcoding + one convention breaks the single- and triple-substrate cases (the two + most common PhysiCell setups).""" + subs = _build_substrate_frame(tmp_path, n_subs) + ts = read_physicell_timestep( + tmp_path / 'output00000000.xml', include_dead_cells=True + ) + inte = ts.internalized_substrates(include_dead_cells=True) + + assert list(inte.columns) == subs + for i, name in enumerate(subs): + expected = (i + 1) * 100 + np.arange(len(inte)) + np.testing.assert_array_equal(inte[name].to_numpy(), expected) + class TestFileResolution: def test_resolves_from_xml_filename(self, frame_xml, example_physicell_dir):