diff --git a/CHANGELOG.md b/CHANGELOG.md index d298d1ad..543e616e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,12 +22,32 @@ Release notes for `v0.7.1` and earlier are available on the [Releases][] page. and the raster cell/nucleus labels are reconstructed by rasterizing the boundary polygons. The cell feature matrix falls back to a `cell_feature_matrix/` MatrixMarket directory or `cell_feature_matrix.tar.gz` when `cell_feature_matrix.h5` is missing. +- `phenocycler()` reader and CLI command for Akoya PhenoCycler-Fusion 1.x/2.x `.qptiff` / `.er.qptiff` scans + (incl. Opal mode) and other PerkinElmer/Akoya QPTIFFs (Vectra Polaris, PhenoImager HT, inForm component data): + lazy images with channel names from the page XML, scan metadata in `sdata.attrs["akoya"]`, optional `.xpd` / + `*channels.csv` sidecars and segmentation masks (`labels`); `read_all_bit_depths=True` also reads the 8-bit scan. ### Changed - Adopted the current `cookiecutter-scverse` template: `hatch`-managed test environments, `uv`-based CI and documentation builds, `mypy` type checking of `src` and `tests`, `biome`/`pyproject-fmt`/`zizmor` pre-commit hooks, and Dependabot updates. +- Requires `tifffile>=2026.5.2` and `zarr>=3.2`, the first versions whose zarr store the lazy Akoya TIFF reads work with. +- Requires `flowio>=1.4`, the first version with `FlowData.pnn_labels`, with which `codex()` detects CODEX FCS tables. +- `codex()` was rewritten to read CODEX Processor, CODEX Uploader and CODEX Instrument Manager (`.qptiff`) + outputs, sharing the QPTIFF core with `phenocycler()`. Breaking changes: + - Element names are `{region}_image`, `{region}_cells`, `{region}_labels` and `table` (were `images`, `{region}`). + - All regions are read (was the first file found); `regions=[...]` subsets. Unrelated regions get one coordinate + system each; regions of one slide are placed at their stage positions in `"global"`. + - `table.var_names` and image channels carry marker names (FCS `$PnS`, CSV `cycCCC_chHHH:`, + `channelNames.txt`), not `cycCCC_chHHH`. + - New keyword-only parameters `regions`, `channel_names`, `cells_table`, `labels`, `*_models_kwargs`; + `imread_kwargs` is ignored with a warning. + +### Fixed + +- `codex()` read the CSV table and the image relative to the current working directory instead of `path`, and an + operator-precedence bug in the table check made it ignore `fcs`. ### Removed diff --git a/README.md b/README.md index adb20e76..01b3c57d 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,8 @@ This package contains reader functions to load common spatial omics formats into - 10x Genomics Visium® - 10x Genomics Visium HD® - 10x Genomics Xenium® -- Akoya PhenoCycler® (formerly CODEX®) +- Akoya PhenoCycler® (formerly CODEX®): CODEX Processor/Uploader/Instrument Manager outputs and PhenoCycler-Fusion 1.x/2.x, + Vectra® Polaris™ / PhenoImager® HT QPTIFFs and inForm® component TIFFs - Curio Seeker® - DBiT-seq - MCMICRO (output data) diff --git a/docs/api.md b/docs/api.md index 95ff1650..71dd3bc5 100644 --- a/docs/api.md +++ b/docs/api.md @@ -25,6 +25,7 @@ I/O for the `spatialdata` project. experimental.iss mcmicro merscope + phenocycler seqfish steinbock stereoseq diff --git a/pyproject.toml b/pyproject.toml index 9c7135e1..2c666895 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,7 @@ dependencies = [ "anndata", "click", "dask-image", + "flowio>=1.4", "h5py", "imagecodecs", "joblib", @@ -40,8 +41,9 @@ dependencies = [ # for debug logging (referenced from the issue template) "session-info2", "spatialdata>=0.7.3a0", - "tifffile>=2023.8.12", + "tifffile>=2026.5.2", "xmltodict", + "zarr>=3.2", ] # https://docs.pypi.org/project_metadata/#project-urls urls.Documentation = "https://spatialdata-io.readthedocs.io/" @@ -156,6 +158,7 @@ exclude = "^tests/data/" [[tool.mypy.overrides]] module = [ "dask_image.*", + "flowio.*", "h5py.*", "multiscale_spatial_image.*", "pyarrow.*", diff --git a/src/spatialdata_io/__init__.py b/src/spatialdata_io/__init__.py index 1c10cc53..8c165b5b 100644 --- a/src/spatialdata_io/__init__.py +++ b/src/spatialdata_io/__init__.py @@ -13,6 +13,7 @@ "macsima": "spatialdata_io.readers.macsima", "mcmicro": "spatialdata_io.readers.mcmicro", "merscope": "spatialdata_io.readers.merscope", + "phenocycler": "spatialdata_io.readers.phenocycler", "seqfish": "spatialdata_io.readers.seqfish", "steinbock": "spatialdata_io.readers.steinbock", "stereoseq": "spatialdata_io.readers.stereoseq", @@ -38,6 +39,7 @@ "macsima", "mcmicro", "merscope", + "phenocycler", "seqfish", "steinbock", "stereoseq", @@ -88,6 +90,7 @@ def __dir__() -> list[str]: from spatialdata_io.readers.macsima import macsima from spatialdata_io.readers.mcmicro import mcmicro from spatialdata_io.readers.merscope import merscope + from spatialdata_io.readers.phenocycler import phenocycler from spatialdata_io.readers.seqfish import seqfish from spatialdata_io.readers.steinbock import steinbock from spatialdata_io.readers.stereoseq import stereoseq diff --git a/src/spatialdata_io/__main__.py b/src/spatialdata_io/__main__.py index f0a5f9f7..c871cb25 100644 --- a/src/spatialdata_io/__main__.py +++ b/src/spatialdata_io/__main__.py @@ -40,36 +40,112 @@ def _input_output_click_options(func: Callable[..., None]) -> Callable[..., None def _parse_json_param(value: str, param_name: str) -> dict[str, Any]: try: - result: dict[str, Any] = json.loads(value) - return result + result = json.loads(value) except json.JSONDecodeError as e: raise click.BadParameter(f"Invalid JSON for {param_name!r}: {e}") from e + if not isinstance(result, dict): + raise click.BadParameter(f"{param_name!r} must be a JSON object, got {value!r}.") + return result + + +def _parse_path_or_mapping_param(value: str | None, param_name: str) -> str | dict[str, Any] | None: + """A path, or a JSON mapping from region to path if ``value`` starts with ``{``.""" + if value is None or not value.lstrip().startswith("{"): + return value + mapping = _parse_json_param(value, param_name) + if not all(isinstance(path, str) for path in mapping.values()): + raise click.BadParameter(f"{param_name!r} must be a path or a JSON object mapping regions to paths.") + return mapping + + +def _akoya_click_options(func: Callable[..., None]) -> Callable[..., None]: + """Decorator adding the options shared by the Akoya readers (``codex``, ``phenocycler``).""" + options = [ + click.option("--output", "-o", type=click.Path(), help="Path to the output.zarr file.", required=True), + click.option("--regions", multiple=True, help="A region to read (repeat for several); all if omitted."), + click.option( + "--channel-names", + multiple=True, + help="A channel name overriding the discovered ones (repeat for every channel).", + ), + click.option( + "--labels", + type=str, + default=None, + help="Path to a segmentation mask, or a JSON mapping from region to path. [default: None]", + ), + click.option( + "--image-models-kwargs", + type=str, + default="{}", + help="JSON string of keyword arguments passed to Image2DModel. [default: {}]", + ), + click.option( + "--labels-models-kwargs", + type=str, + default="{}", + help="JSON string of keyword arguments passed to Labels2DModel. [default: {}]", + ), + ] + for option in reversed(options): + func = option(func) + return func @cli.command(name="codex") -@_input_output_click_options +@click.option( + "--input", + "-i", + type=click.Path(exists=True), + help="Path to a CODEX output directory or to a single .qptiff/.tif image.", + required=True, +) +@_akoya_click_options +@click.option( + "--cells-table", + type=str, + default=None, + help='Path to a CODEX cell table, or a JSON mapping from region to path, e.g. \'{"reg001": "a.csv"}\'. ' + "[default: None]", +) @click.option( "--fcs", type=bool, default=True, - help="Whether the .fcs file is provided if False a .csv file is expected. [default: True]", + help="Whether to prefer .fcs over .csv cell tables when both exist. [default: True]", ) @click.option( "--imread-kwargs", type=str, default="{}", - help="JSON string of keyword arguments passed to imread. [default: {}]", + help="Ignored, with a warning; kept for backward compatibility. [default: {}]", ) def codex_wrapper( input: str, output: str, fcs: bool = True, + regions: tuple[str, ...] = (), + channel_names: tuple[str, ...] = (), + cells_table: str | None = None, + labels: str | None = None, imread_kwargs: str = "{}", + image_models_kwargs: str = "{}", + labels_models_kwargs: str = "{}", ) -> None: """Codex conversion to SpatialData.""" from spatialdata_io.readers.codex import codex - sdata = codex(input, fcs=fcs, imread_kwargs=_parse_json_param(imread_kwargs, "imread_kwargs")) + sdata = codex( + input, + fcs=fcs, + regions=list(regions) or None, + channel_names=list(channel_names) or None, + cells_table=_parse_path_or_mapping_param(cells_table, "cells_table"), + labels=_parse_path_or_mapping_param(labels, "labels"), + imread_kwargs=_parse_json_param(imread_kwargs, "imread_kwargs"), + image_models_kwargs=_parse_json_param(image_models_kwargs, "image_models_kwargs"), + labels_models_kwargs=_parse_json_param(labels_models_kwargs, "labels_models_kwargs"), + ) sdata.write(output) @@ -357,6 +433,46 @@ def merscope_wrapper( sdata.write(output) +@cli.command(name="phenocycler") +@click.option( + "--input", + "-i", + type=click.Path(exists=True), + help="Path to a .qptiff file or to a directory containing .qptiff files.", + required=True, +) +@_akoya_click_options +@click.option( + "--read-all-bit-depths/--no-read-all-bit-depths", + default=False, + help="Read both the 8-bit .qptiff and the extended-range .er.qptiff of a scan (the 8-bit one as region " + "'_8bit'). [default: no-read-all-bit-depths]", +) +def phenocycler_wrapper( + input: str, + output: str, + regions: tuple[str, ...] = (), + channel_names: tuple[str, ...] = (), + labels: str | None = None, + read_all_bit_depths: bool = False, + image_models_kwargs: str = "{}", + labels_models_kwargs: str = "{}", +) -> None: + """PhenoCycler-Fusion conversion to SpatialData.""" + from spatialdata_io.readers.phenocycler import phenocycler + + sdata = phenocycler( + input, + regions=list(regions) or None, + channel_names=list(channel_names) or None, + labels=_parse_path_or_mapping_param(labels, "labels"), + read_all_bit_depths=read_all_bit_depths, + image_models_kwargs=_parse_json_param(image_models_kwargs, "image_models_kwargs"), + labels_models_kwargs=_parse_json_param(labels_models_kwargs, "labels_models_kwargs"), + ) + sdata.write(output) + + @cli.command(name="seqfish") @_input_output_click_options @click.option("--load-images", type=bool, default=True, help="Whether to load images. [default: True]") diff --git a/src/spatialdata_io/_constants/_constants.py b/src/spatialdata_io/_constants/_constants.py index 36e2a88f..cfbbd08a 100644 --- a/src/spatialdata_io/_constants/_constants.py +++ b/src/spatialdata_io/_constants/_constants.py @@ -10,6 +10,9 @@ class CodexKeys(ModeEnum): # files and directories FCS_FILE = ".fcs" FCS_FILE_CSV = ".csv" + CHANNEL_NAMES_FILE = "channelNames.txt" + EXPERIMENT_FILE = "experiment.json" + QPTIFF = ".qptiff" # metadata REGION_KEY = "region" INSTANCE_KEY = "cell_id" @@ -18,6 +21,26 @@ class CodexKeys(ModeEnum): IMAGE_TIF = ".tif" +@unique +class PhenocyclerKeys(ModeEnum): + """Keys for *Akoya PhenoCycler-Fusion* (and CIM QPTIFF) formatted dataset.""" + + # files and directories + QPTIFF_SUFFIX = ".qptiff" + ER_QPTIFF_SUFFIX = ".er.qptiff" + XPD_SUFFIX = ".xpd" + CHANNELS_CSV_SUFFIX = "channels.csv" + # elements + IMAGE_SUFFIX = "_image" + LABELS_SUFFIX = "_labels" + CELLS_SUFFIX = "_cells" + TABLE_KEY = "table" + # metadata + REGION_KEY = "region" + INSTANCE_KEY = "cell_id" + ATTRS_KEY = "akoya" + + class CurioKeys(ModeEnum): """Keys for *Curio* formatted dataset.""" diff --git a/src/spatialdata_io/readers/_utils/_akoya.py b/src/spatialdata_io/readers/_utils/_akoya.py new file mode 100644 index 00000000..760861f7 --- /dev/null +++ b/src/spatialdata_io/readers/_utils/_akoya.py @@ -0,0 +1,420 @@ +"""Channel, transformation, image and metadata helpers shared by the Akoya readers.""" + +from __future__ import annotations + +import io +import json +import math +import re +import threading +import warnings +from collections import Counter +from collections.abc import Iterable, Mapping, Sequence, Set +from dataclasses import asdict +from typing import TYPE_CHECKING, Any + +import dask.array as da +import numpy as np +import tifffile +import zarr +from dask.base import tokenize +from spatialdata import sanitize_name +from spatialdata._logging import logger +from spatialdata.models import Image2DModel +from spatialdata.transformations import BaseTransformation, Identity, Translation + +from spatialdata_io.readers._utils._utils import _initialize_raster_models_kwargs + +if TYPE_CHECKING: + from pathlib import Path + + from xarray import DataArray, DataTree + + from spatialdata_io.readers._utils._qptiff import AkoyaScan + +__all__ = [ + "IMAGE_EXTENSION", + "apply_channel_names", + "build_image_element", + "dedupe_channel_names", + "lazy_tiff_array", + "normalize_regions", + "read_json_object", + "read_tiff", + "region_id", + "region_transformations", + "scans_to_attrs", + "select_regions", + "spatial_chunks", +] + +# the extensions of image files, removed to name a region after its file +IMAGE_EXTENSION = re.compile(r"(?:\.er)?\.(?:qptiff|ome\.tiff?|tiff?)$", re.IGNORECASE) +# `reg001`, `exp_reg001` or `1_reg001`, but not `Treg2` or `coreg1` +_REGION = re.compile(r"(? str | None: + """The CODEX region id ``regRRR`` in a file name: ``reg`` and digits, not preceded by a letter; else `None`.""" + match = _REGION.search(name) + return f"reg{int(match.group(1)):03d}" if match else None + + +def dedupe_channel_names(names: Sequence[str]) -> list[str]: + """Make channel names unique as ``name``, ``name (2)``, ``name (3)``; the first occurrence is unchanged. + + A suffix is skipped if it would repeat a name of the list, e.g. ``["A", "A", "A (2)"]`` -> ``A``, ``A (3)``, + ``A (2)``. + """ + taken = set(names) + seen: Counter[str] = Counter() + deduped = [] + for name in names: + seen[name] += 1 + if seen[name] == 1: + deduped.append(name) + continue + while (unique := f"{name} ({seen[name]})") in taken: + seen[name] += 1 + taken.add(unique) + deduped.append(unique) + return deduped + + +def apply_channel_names(scan: AkoyaScan, channel_names: Sequence[str] | None) -> list[str]: + """Return the de-duplicated channel names of ``scan``: the user override if given, else the metadata names.""" + if channel_names is None: + return dedupe_channel_names([channel.name for channel in scan.channels]) + if isinstance(channel_names, Mapping | Set): + raise TypeError(f"`channel_names` must be a sequence of names, got {type(channel_names).__name__}.") + if len(channel_names) != len(scan.channels): + raise ValueError( + f"`channel_names` has {len(channel_names)} entries, but region '{scan.region}' ({scan.path}) has " + f"{len(scan.channels)} channels: expected {len(scan.channels)}." + ) + return dedupe_channel_names(list(channel_names)) + + +def region_transformations( + scans: Sequence[AkoyaScan], regions_without_image: Iterable[str] = () +) -> dict[str, dict[str, BaseTransformation]]: + """Compute the coordinate systems and transformations of every region, by region id. + + Akoya scans are already stitched, so scans of the same sample space are only placed: if there are several + scans and all have a stage position and a pixel size, the same slide and pixel sizes equal within 0.1%, each is + translated by its stage position (in pixels, relative to the top-left one) into the ``"global"`` coordinate + system. Otherwise the regions are unrelated and each gets its own coordinate system, named after the region, + with ``Identity``; a read of a single region uses ``"global"``. + + Parameters + ---------- + scans + The scans read. + regions_without_image + Regions read without an image (e.g. only their cells); they are unrelated to the scans. + + Returns + ------- + ``{region: {coordinate_system: transformation}}``, to be passed to every element of the region. + """ + imaged = [scan.region for scan in scans] + others = sorted(set(regions_without_image) - set(imaged)) + if len(imaged) + len(others) == 1: + return {region: {"global": Identity()} for region in [*imaged, *others]} + + result: dict[str, dict[str, BaseTransformation]] = {} + if _same_sample_space(scans): + min_x = min(scan.stage_position_um[0] for scan in scans) # type: ignore[index] + min_y = min(scan.stage_position_um[1] for scan in scans) # type: ignore[index] + for scan in scans: + x, y = scan.stage_position_um # type: ignore[misc] + pixel_size = scan.pixel_size_um + translation = [(x - min_x) / pixel_size, (y - min_y) / pixel_size] # type: ignore[operator] + result[scan.region] = {"global": Translation(translation, axes=("x", "y"))} + else: + if len(imaged) > 1: + logger.info( + "Regions %s do not share a sample space (missing stage position or pixel size, different slides " + "or pixel sizes), so each is stored in its own coordinate system, named after the region.", + imaged, + ) + result = {region: {region: Identity()} for region in imaged} + return result | {region: {region: Identity()} for region in others} + + +def _same_sample_space(scans: Sequence[AkoyaScan]) -> bool: + """Whether several scans can be placed by their stage positions in one pixel coordinate system.""" + if len(scans) < 2 or not all(scan.stage_position_um is not None and scan.pixel_size_um for scan in scans): + return False + slides = {scan.slide_id for scan in scans} + pixel_size = scans[0].pixel_size_um + return ( + len(slides) == 1 + and None not in slides + and all(math.isclose(scan.pixel_size_um, pixel_size, rel_tol=1e-3) for scan in scans) # type: ignore[arg-type] + ) + + +def build_image_element( + image: da.Array, + scan: AkoyaScan, + c_coords: list[str], + transformations: Mapping[str, BaseTransformation], + image_models_kwargs: Mapping[str, Any], +) -> DataArray | DataTree: + """Parse a lazy ``(c, y, x)`` array, or ``(y, x, 3)`` for brightfield RGB scans, into an ``Image2DModel``. + + ``transformations`` are the region's, see :func:`region_transformations`. + """ + image_models_kwargs, _ = _initialize_raster_models_kwargs(image_models_kwargs, {}) + if scan.rgb: + image = image.transpose(2, 0, 1) + c_coords = ["r", "g", "b"] + return Image2DModel.parse( + image, + dims=("c", "y", "x"), + c_coords=c_coords, + rgb=scan.rgb, + transformations=dict(transformations), + **image_models_kwargs, + ) + + +def scans_to_attrs(scans: Sequence[AkoyaScan], channel_names: Mapping[str, list[str]]) -> dict[str, Any]: + """Build the JSON-serialisable per-region metadata stored in ``sdata.attrs["akoya"]``.""" + regions: dict[str, Any] = {} + for scan in scans: + names = channel_names[scan.region] + regions[scan.region] = { + "generation": scan.generation, + "software": scan.software, + "slide_id": scan.slide_id, + "pixel_size_um": scan.pixel_size_um, + "nominal_pixel_size_um": scan.nominal_pixel_size_um, + "stage_position_um": list(scan.stage_position_um) if scan.stage_position_um is not None else None, + "path": str(scan.path), + "channels": [ + {**asdict(channel), "name_deduped": name} for channel, name in zip(scan.channels, names, strict=True) + ], + "raw": scan.raw, + } + return {"regions": regions} + + +def normalize_regions(regions: str | Sequence[str] | None) -> list[str] | None: + """Validate the ``regions`` argument of a reader: a region id or a sequence of them, without duplicates. + + The region ids are sanitised with :func:`spatialdata.sanitize_name`, like those of the datasets, so raw names + (e.g. ``"Slide 1"`` for the region ``Slide_1``) select their region. Raises a `TypeError` for a mapping and a + `ValueError` for an empty sequence. + """ + if regions is None: + return None + if isinstance(regions, str): + return [sanitize_name(regions)] + if isinstance(regions, Mapping): + raise TypeError(f"`regions` must be a sequence of region ids, got {type(regions).__name__}.") + sanitized = list(dict.fromkeys(sanitize_name(str(region)) for region in regions)) + if not sanitized: + raise ValueError("`regions` is empty; pass `None` to read all regions.") + return sanitized + + +def select_regions(available: Sequence[str], regions: str | Sequence[str] | None) -> list[str]: + """Return the requested regions (all if ``None``), raising if one is not available; see :func:`normalize_regions`.""" + regions = normalize_regions(regions) + if regions is None: + return list(available) + unknown = [region for region in regions if region not in available] + if unknown: + raise ValueError(f"Unknown region(s) {unknown}. Available regions: {list(available)}.") + return regions + + +def read_json_object(path: Path) -> dict[str, Any] | None: + """Read a JSON object from ``path``; warn and return ``None`` if the file is not valid UTF-8 JSON or not an object. + + A UTF-8 byte order mark, as written by some Windows tools, is ignored. + """ + try: + content = json.loads(path.read_text(encoding="utf-8-sig")) + except (json.JSONDecodeError, UnicodeDecodeError) as e: + warnings.warn(f"Ignoring {path}: it is not valid JSON ({e}).", UserWarning, stacklevel=2) + return None + if not isinstance(content, dict): + warnings.warn(f"Ignoring {path}: it is not a JSON object.", UserWarning, stacklevel=2) + return None + return content + + +def read_tiff(path: Path, chunks: tuple[int, int]) -> tuple[da.Array, tuple[int, ...], list[str] | None]: + """Lazily read a tiff as ``(c, y, x)``, all non-spatial axes flattened into ``c``. + + Uncompressed contiguous data (e.g. ImageJ hyperstacks and per-marker images, often a single strip per plane) + is memory-mapped and split into ``chunks``, so a chunk never needs a whole plane in memory; other layouts keep + their native tiles via zarr. Returns the image, the shape of the flattened axes and the ImageJ ``Labels`` (if + every plane has a non-blank one). The array opens the file only when its chunks are read and can be pickled, + e.g. for ``dask.distributed``. + """ + with tifffile.TiffFile(path) as tf: + series = tf.series[0] + axes, offset, shape = series.axes, series.dataoffset, series.shape + # tifffile reports the dtype in native byte order; a memory map needs the file's (big-endian for ImageJ) + dtype = series.dtype.newbyteorder(tf.byteorder) + metadata = tf.imagej_metadata or {} + labels = list(metadata["Labels"]) if isinstance(metadata.get("Labels"), list | tuple) else None + # ImageJ writes an empty label for every slice without one, so partial labels do not name the channels + if labels is not None and not all(isinstance(label, str) and label.strip() for label in labels): + labels = None + if offset is not None: + source = _ContiguousTiff(path, offset, shape, dtype) + tiles = dict.fromkeys(range(source.ndim), 1) | {axes.index("Y"): chunks[0], axes.index("X"): chunks[1]} + tile_chunks = tuple(tiles[i] for i in range(source.ndim)) + image = da.from_array(source, chunks=tile_chunks, name=f"{source.name}-{tokenize(tile_chunks)}") + else: + image = lazy_tiff_array(path, series=0, chunks=chunks) + image = da.moveaxis(image, [axes.index("Y"), axes.index("X")], [-2, -1]) + leading_shape = tuple(image.shape[:-2]) + return image.reshape((math.prod(leading_shape), *image.shape[-2:])), leading_shape, labels + + +class _ContiguousTiff: + """Array-like view of uncompressed contiguous tiff data, memory-mapped anew for every chunk. + + Handing an ``np.memmap`` to :func:`dask.array.from_array` would read the whole file (dask hashes and slices + numpy arrays eagerly), so dask only gets this lightweight object and a name derived from the file and chunks. + """ + + def __init__(self, path: Path, offset: int, shape: tuple[int, ...], dtype: np.dtype[Any]) -> None: + self.path = path + self.offset = offset + self.shape = shape + self.ndim = len(shape) + self._file_dtype = dtype + self.dtype = dtype.newbyteorder("=") + self.name = f"akoya-tiff-{tokenize(str(path), path.stat().st_mtime_ns, offset, shape, str(dtype))}" + + def __getitem__(self, key: Any) -> np.ndarray: + data = np.memmap(self.path, dtype=self._file_dtype, mode="r", offset=self.offset, shape=self.shape) + # a copy, as a view would keep the memory map and its file descriptor open as long as the chunk lives + return np.array(data[key], dtype=self.dtype) + + +def lazy_tiff_array( + path: Path, *, series: int = 0, page: int | None = None, chunks: tuple[int, int] = (4096, 4096) +) -> da.Array: + """Lazily read the level 0 of a tiff series, or one page, in chunks of whole native tiles or strips. + + The array can be pickled (e.g. for ``dask.distributed`` or process-based schedulers): its chunks are read by a + :class:`_TiffZarrSource`, which opens the file in each process that reads them. + + Parameters + ---------- + path + The tiff file. + series + The index of the series to read, if ``page`` is `None`. + page + The index of the page to read instead of a series. + chunks + The ``(y, x)`` chunk size to aim for: the native tiles or strips are grouped into chunks of about as many + pixels (the default is that of the raster models), so the graph does not get one task per strip. + + Returns + ------- + The array, with the shape of the series or page. + """ + source = _TiffZarrSource(path, series=series, page=page) + y, x = source.axes.index("Y"), source.axes.index("X") + grouped = list(source.chunks) + native_y, native_x = grouped[y], grouped[x] + grouped[x] = min(source.shape[x], max(native_x, chunks[1] // native_x * native_x)) + grouped[y] = min(source.shape[y], max(native_y, chunks[0] * chunks[1] // grouped[x] // native_y * native_y)) + return da.from_array( + source, + chunks=tuple(grouped), + name=f"{source.name}-{tokenize(grouped)}", + meta=np.empty((0,) * source.ndim, dtype=source.dtype), + ) + + +class _TiffZarrSource: + """Picklable array-like view of a tiff series (level 0) or page, read chunk-wise through tifffile's zarr store. + + Neither the ``TiffFile`` nor its zarr store are pickled, but opened again in each process that reads chunks, so + dask graphs built on it work across processes. The file is read through a :class:`_TransientFile`, so the parsed + tiff stays cached without holding a file descriptor. + """ + + def __init__(self, path: Path, *, series: int = 0, page: int | None = None) -> None: + self.path = path + self.series = series + self.page = page + self._lock = threading.Lock() + self._array: zarr.Array | None = None + array = self._open() + self.shape = tuple(array.shape) + self.chunks = tuple(array.chunks) + self.dtype = array.dtype + self.ndim = len(self.shape) + token = tokenize(str(path), path.stat().st_mtime_ns, series, page, self.shape, str(self.dtype)) + self.name = f"akoya-tiff-zarr-{token}" + + def _open(self) -> zarr.Array: + with self._lock: + if self._array is None: + tiff = tifffile.TiffFile(io.BufferedReader(_TransientFile(self.path))) + data = tiff.pages[self.page] if self.page is not None else tiff.series[self.series] + self.axes = data.axes + self._array = zarr.open_array(data.aszarr(level=0), mode="r") + return self._array + + def __getitem__(self, key: Any) -> np.ndarray: + return np.asarray(self._open()[key]) + + def __getstate__(self) -> dict[str, Any]: + return {key: value for key, value in self.__dict__.items() if key not in ("_lock", "_array")} + + def __setstate__(self, state: dict[str, Any]) -> None: + self.__dict__.update(state, _lock=threading.Lock(), _array=None) + + +class _TransientFile(io.RawIOBase): + """Read-only binary file that opens ``path`` for every read, so it holds no file descriptor in between.""" + + def __init__(self, path: Path) -> None: + super().__init__() + self.name = str(path) + self._size = path.stat().st_size + self._position = 0 + + def readable(self) -> bool: + return True + + def seekable(self) -> bool: + return True + + def tell(self) -> int: + return self._position + + def seek(self, offset: int, whence: int = io.SEEK_SET) -> int: + self._position = (0, self._position, self._size)[whence] + offset + return self._position + + def readinto(self, buffer: Any) -> int: + with open(self.name, "rb") as f: + f.seek(self._position) + n = f.readinto(buffer) + self._position += n + return n + + +def spatial_chunks(chunks: Any, default: Sequence[int]) -> tuple[int, int]: + """The ``(y, x)`` sizes of a raster model ``chunks`` argument, ``default`` when it does not give plain sizes.""" + if isinstance(chunks, int): + return chunks, chunks + if isinstance(chunks, Mapping) and isinstance(chunks.get("y"), int) and isinstance(chunks.get("x"), int): + return chunks["y"], chunks["x"] + if isinstance(chunks, list | tuple) and len(chunks) >= 2 and all(isinstance(c, int) for c in chunks[-2:]): + return chunks[-2], chunks[-1] + return default[-2], default[-1] diff --git a/src/spatialdata_io/readers/_utils/_cell_tables.py b/src/spatialdata_io/readers/_utils/_cell_tables.py new file mode 100644 index 00000000..eee74994 --- /dev/null +++ b/src/spatialdata_io/readers/_utils/_cell_tables.py @@ -0,0 +1,726 @@ +"""CODEX Processor cell tables and segmentation masks shared by the Akoya readers, and the assembly of their output.""" + +from __future__ import annotations + +import csv +import re +import warnings +from collections.abc import Mapping +from dataclasses import dataclass +from pathlib import Path +from typing import TYPE_CHECKING, Any + +import anndata as ad +import flowio # readfcs' parser, for the TEXT segment only +import numpy as np +import pandas as pd +import readfcs +from spatialdata import SpatialData, sanitize_name, sanitize_table +from spatialdata._logging import logger +from spatialdata.models import Labels2DModel, ShapesModel, TableModel +from spatialdata.transformations import Identity, Scale, set_transformation +from spatialdata.transformations import Sequence as TransformationSequence +from xarray import DataTree + +from spatialdata_io._constants._constants import CodexKeys, PhenocyclerKeys +from spatialdata_io.readers._utils._akoya import ( + build_image_element, + dedupe_channel_names, + read_tiff, + region_transformations, + scans_to_attrs, + spatial_chunks, +) +from spatialdata_io.readers._utils._utils import _initialize_raster_models_kwargs + +if TYPE_CHECKING: + from collections.abc import Callable, Collection, Sequence + + import dask.array as da + from geopandas import GeoDataFrame + from spatialdata.transformations import BaseTransformation + from xarray import DataArray + + from spatialdata_io.readers._utils._qptiff import AkoyaScan + +__all__ = [ + "CellTable", + "build_cells", + "build_spatialdata", + "is_codex_table", + "mask_downsample_factor", + "raster_shape", + "read_cell_table", + "read_labels", + "read_region_labels", + "read_region_tables", + "resolve_per_region", +] + +# `cyc001_ch001` (FCS `$PnN`, plain CSV header) or the name part of the Processor CSV header `cyc001_ch001:DAPI1` +_MARKER_COLUMN = re.compile(r"cyc(\d+)_ch(\d+)", re.IGNORECASE) +# non-marker columns of CODEX Processor tables, for tables whose marker columns are named after the markers +_META_COLUMNS = frozenset( + {"cell_id", "region", "tile_num", "x", "y", "z", "x_tile", "y_tile", "size", "tsne_x", "tsne_y", "homogeneity"} +) + + +@dataclass +class CellTable: + """A CODEX cell table. + + Attributes + ---------- + obs + The non-marker columns, indexed by the string instance ids. + X + The marker intensities, ``(n_cells, n_markers)`` float32. + var + Per-marker metadata indexed by the de-duplicated marker names: ``source_column`` and, for + ``cycCCC_chHHH`` columns, ``cycle`` and ``channel``. + region + The region id of every cell. + instance_id + The instance id of every cell (int64). + centroids_px + The ``(x, y)`` centroids in level-0 pixels of the region image. + radius_px + The radius of a circle with the cell's area, or `None` if the table has no area. + row_instance_ids + Whether the instance ids are row numbers because the table has no instance id column. + """ + + obs: pd.DataFrame + X: np.ndarray + var: pd.DataFrame + region: pd.Series + instance_id: pd.Series + centroids_px: np.ndarray + radius_px: np.ndarray | None + row_instance_ids: bool = False + + def take(self, rows: np.ndarray) -> CellTable: + """The cells at the positional indices ``rows``.""" + return CellTable( + obs=self.obs.iloc[rows], + X=self.X[rows], + var=self.var, + region=self.region.iloc[rows], + instance_id=self.instance_id.iloc[rows], + centroids_px=self.centroids_px[rows], + radius_px=None if self.radius_px is None else self.radius_px[rows], + row_instance_ids=self.row_instance_ids, + ) + + +def is_codex_table(path: Path) -> bool: + """Whether ``path`` is a CODEX Processor ``.fcs`` or ``.csv`` cell table. + + Its columns (the ``$PnN`` channel names of an FCS file) must include ``x`` and ``y`` and ``cell_id`` or a + ``cycCCC_chHHH`` column. Other CSV files (e.g. antibody lists) and FCS files of other instruments, also + unreadable ones, are not. + """ + if path.suffix.lower() == CodexKeys.FCS_FILE.value: + try: + columns = list(flowio.read_multiple_data_sets(str(path), only_text=True)[0].pnn_labels) + # flowio raises these, not only its own exceptions, for malformed files (e.g. truncated, without `$DATATYPE`) + except (ValueError, KeyError, IndexError, EOFError, OSError, flowio.exceptions.FlowIOException): + return False + else: + columns = _read_header(path) + # the name part of the CODEX Processor header form `cyc001_ch001:DAPI1` + names = {column.partition(":")[0].casefold() for column in columns} + return {"x", "y"} <= names and ( + CodexKeys.INSTANCE_KEY.value in names or any(_MARKER_COLUMN.fullmatch(name) for name in names) + ) + + +def read_cell_table(path: Path, *, default_region: str | None, regions: Collection[str] | None = None) -> CellTable: + """Read a CODEX Processor ``.fcs`` or ``.csv`` cell table. + + The markers are the ``cycCCC_chHHH`` columns, named by the FCS ``$PnS`` keywords or by the display part of the + Processor CSV header ``cycCCC_chHHH:`` (in ``var`` with their ``cycle`` and ``channel``), else, with a + `UserWarning`, the numeric columns that are not CODEX Processor columns (``x``, ``size``, ...). Centroids + ``x``/``y`` in pixels, instance ids ``cell_id`` (else the row number in the file, from 1), the region from + ``region`` (a numeric region ``r`` is the region id ``regRRR``), the radius from ``size``. + + Parameters + ---------- + path + The cell table. + default_region + The region id of all cells if the table has no region column. + regions + The region ids of the dataset, to which the values of the region column are matched, also after + :func:`spatialdata.sanitize_name`; values that name none of them are kept. If the table has + ``default_region`` and its region column holds a single value that names no region, that value names + ``default_region``. If `None`, the values are the region ids. + + Returns + ------- + The cell table. + """ + df = _read_fcs(path) if path.suffix.lower() == CodexKeys.FCS_FILE.value else _read_text_table(path) + headers = [str(column) for column in df.columns] + split = [header.partition(":") for header in headers] + # the source of an FCS column is its `$PnN` + fcs = path.suffix.lower() == CodexKeys.FCS_FILE.value + sources = [name if fcs else header for header, (name, _, _) in zip(headers, split, strict=True)] + names = dedupe_channel_names([name for name, _, _ in split]) + displays = [display.strip() or name for name, (_, _, display) in zip(names, split, strict=True)] + df.columns = pd.Index(names) + columns = {column.casefold(): column for column in names} + x, y = columns.get("x"), columns.get("y") + if x is None or y is None: + raise ValueError(f"{path} has no 'x' and 'y' centroid columns; found {names}.") + instance_column = columns.get(CodexKeys.INSTANCE_KEY.value) + region_column = columns.get(CodexKeys.REGION_KEY.value) + area_column = columns.get("size") + + matches = [_MARKER_COLUMN.fullmatch(name) for name in names] + if any(matches): + is_marker = [match is not None for match in matches] + else: + is_marker = [ + name.casefold() not in _META_COLUMNS and pd.api.types.is_numeric_dtype(df.iloc[:, i]) + for i, name in enumerate(names) + ] + taken = [name for name, marker in zip(names, is_marker, strict=True) if marker] + warnings.warn( + f"{path} has no 'cycCCC_chHHH' columns; reading the {len(taken)} numeric columns {taken} that are not " + "CODEX Processor columns as the markers.", + UserWarning, + stacklevel=4, + ) + positions = [i for i, marker in enumerate(is_marker) if marker] + markers = [names[i] for i in positions] + var = pd.DataFrame( + {"source_column": [sources[i] for i in positions]}, + index=pd.Index(dedupe_channel_names([displays[i] for i in positions])), + ) + if positions and all(matches[i] for i in positions): + var["cycle"] = [int(matches[i].group(1)) for i in positions] # type: ignore[union-attr] + var["channel"] = [int(matches[i].group(2)) for i in positions] # type: ignore[union-attr] + + if instance_column is not None: + instance_id = pd.Series(_instance_ids(path, df[instance_column])) + else: + # the row numbers in the file; label images use 0 for the background, so instance ids start at 1 + instance_id = pd.Series(np.arange(1, len(df) + 1, dtype=np.int64)) + index = pd.Index(instance_id.astype(str)) + instance_id.index = index + instance_id.name = CodexKeys.INSTANCE_KEY.value + + region = _region_values(df, region_column, default_region, regions) + if region is None: + raise ValueError( + f"{path} has no region column, so its cells cannot be assigned to regions: a table for several regions " + "needs a 'region' column. Alternatively, pass one table per region as a mapping {region: path}." + ) + region = pd.Series(region.to_numpy(), index=index, name=CodexKeys.REGION_KEY.value) + _check_unique_instance_ids(path, region, instance_id) + + # the instance ids and regions are stored under the keys of the table element + obs = df.drop(columns=[*markers, *(c for c in (instance_column, region_column) if c is not None)]) + obs.index = index + radius = None if area_column is None else np.sqrt(df[area_column].to_numpy(dtype=np.float64) / np.pi) + return CellTable( + obs=obs, + X=df[markers].to_numpy(dtype=np.float32), + var=var, + region=region, + instance_id=instance_id, + centroids_px=df[[x, y]].to_numpy(dtype=np.float64), + radius_px=radius, + row_instance_ids=instance_column is None, + ) + + +def _region_values( + df: pd.DataFrame, region_column: str | None, default_region: str | None, regions: Collection[str] | None +) -> pd.Series | None: + """The region of every row of ``df``, see :func:`read_cell_table`; `None` if it is unknown.""" + if region_column is None: + return None if default_region is None else pd.Series(default_region, index=df.index) + values = df[region_column] + if pd.api.types.is_numeric_dtype(values): + values = values.map({value: f"reg{int(value):03d}" for value in values.dropna().unique()}) + values = values.astype(str) + if regions is None: + return values + matched = {value: _match_region(value, regions) for value in values.unique()} + if default_region is not None and len(matched) == 1 and None in matched.values(): + # the region column of a table for one region names another region (e.g. `1` of an image `Kidney.tif`) + return pd.Series(default_region, index=df.index) + return values.map({value: value if region is None else region for value, region in matched.items()}) + + +def _match_region(value: str, regions: Collection[str]) -> str | None: + """The region named by ``value``, also after :func:`spatialdata.sanitize_name` (``A [1]`` names ``A__1_``).""" + return next((candidate for candidate in (value, sanitize_name(value)) if candidate in regions), None) + + +def _read_fcs(path: Path) -> pd.DataFrame: + """The events of an FCS file, with the columns named ``$PnN:$PnS`` like the CODEX Processor CSV header.""" + fcs = readfcs.ReadFCS(path) + df = fcs.data # positional columns, named by `$PnN` + labels = fcs.channels["PnS"] if "PnS" in fcs.channels.columns else [""] * df.shape[1] + names = fcs.channels["PnN"] + df.columns = pd.Index( + [f"{name}:{label.strip()}" if label.strip() else name for name, label in zip(names, labels, strict=True)] + ) + return df + + +def _instance_ids(path: Path, column: pd.Series) -> np.ndarray: + """The integer instance ids in ``column``; raises for missing or non-integer ids, which no label can have.""" + if pd.api.types.is_integer_dtype(column): + return column.to_numpy(dtype=np.int64) + ids = pd.to_numeric(column, errors="coerce").to_numpy(dtype=np.float64) + invalid = ~(np.isfinite(ids) & (ids == np.round(ids))) + if invalid.any(): + raise ValueError( + f"{path} has {invalid.sum()} cell(s) whose instance id in '{column.name}' is missing or not an integer " + f"(e.g. {column[invalid].tolist()[:5]}); every cell needs an integer instance id, the value of its label." + ) + return ids.astype(np.int64) + + +def _check_unique_instance_ids(path: Path, region: pd.Series, instance_id: pd.Series) -> None: + """Raise if two cells of the same region have the same instance id: they could not be told apart.""" + regions, ids = region.to_numpy(), instance_id.to_numpy() + duplicated = pd.DataFrame({"region": regions, "id": ids}).duplicated().to_numpy() + if duplicated.any(): + first = regions[duplicated][0] + found = sorted(set(ids[duplicated & (regions == first)].tolist())) + raise ValueError( + f"{path} has {len(found)} duplicate instance id(s) {found[:5]} in region '{first}'; every cell of a " + "region needs its own instance id." + ) + + +def _read_header(path: Path) -> list[str]: + """The column names of a comma-separated table.""" + with path.open(encoding="utf-8-sig", errors="replace") as f: + header = f.readline().rstrip("\r\n") + return [column.strip() for column in next(csv.reader([header]))] + + +def _read_text_table(path: Path) -> pd.DataFrame: + """Read a comma-separated table with unique column names.""" + df = pd.read_csv(path, engine="pyarrow") + # pandas renames repeated names (DAPI.1; with the pyarrow engine since 3.1), so the names come from the header and + # are made unique like the marker names + df.columns = pd.Index(dedupe_channel_names(_read_header(path))) + return df + + +def resolve_per_region( + value: str | Path | Mapping[str, str | Path] | None, + regions: Sequence[str], + available: Sequence[str] = (), +) -> dict[str, Path]: + """Resolve a ``cells_table`` or ``labels`` argument to one file per region. + + Parameters + ---------- + value + A path, a mapping from region id to path, or `None`. The keys are sanitised with + :func:`spatialdata.sanitize_name`, like the region ids. + regions + The regions read. + available + All regions of the dataset; mapping keys that are one of them but not read are ignored. + + Returns + ------- + ``{}`` for `None`; ``{region: path}`` for a path and a single region, else ``{"*": path}`` (a file for all + regions); the mapping, in the order of ``regions``. + """ + if value is None: + return {} + if isinstance(value, Mapping): + by_region = {sanitize_name(str(key)): file for key, file in value.items()} + if len(by_region) < len(value): + raise ValueError( + f"Several of the keys {list(value)} name the same region after `spatialdata.sanitize_name`; pass " + "each region once." + ) + unknown = sorted(set(by_region) - set(regions) - set(available)) + if unknown: + raise ValueError(f"Unknown region(s) {unknown}. Available regions: {list(available or regions)}.") + files = {region: Path(by_region[region]) for region in regions if region in by_region} + else: + files = {regions[0] if len(regions) == 1 else "*": Path(value)} + for file in files.values(): + if not file.is_file(): + raise FileNotFoundError(f"{file} does not exist or is not a file.") + return files + + +def read_region_tables( + files: Mapping[str, Path], *, available: Collection[str], selected: Sequence[str] +) -> list[CellTable]: + """Read cell tables, as resolved by :func:`resolve_per_region`, and assign their cells to regions. + + A table given for one region holds the cells of that region; if its region column holds a single value that + names no region, it holds only cells of that region. Cells whose region names no available region are dropped + with a warning; cells of available regions that are not read are dropped silently. + + Parameters + ---------- + files + The table of each region, or ``{"*": path}`` for one table for all regions. + available + All region ids of the dataset. + selected + The region ids read. + + Returns + ------- + The cell tables, one per file. + """ + tables = [] + for key, path in files.items(): + table = read_cell_table(path, default_region=None if key == "*" else key, regions=available) + tables.append(_assign_regions(table, path, key, available, selected)) + return tables + + +def _assign_regions( + table: CellTable, path: Path, key: str, available: Collection[str], selected: Sequence[str] +) -> CellTable: + values = list(table.region.unique()) + unmatched = sorted(value for value in values if value not in available) + if unmatched: + warnings.warn( + f"Ignoring the cells of {unmatched} in {path}: they match none of the regions {sorted(available)}.", + UserWarning, + stacklevel=4, + ) + keep = table.region.isin(selected if key == "*" else [key]).to_numpy() + if not keep.any(): + target = f"region '{key}'" if key != "*" else f"the regions read {list(selected)}" + raise ValueError(f"{path} has no cells of {target}; its region column holds {values}.") + return table.take(np.flatnonzero(keep)) + + +def read_labels( + path: Path, labels_models_kwargs: Mapping[str, Any], transformations: Mapping[str, BaseTransformation] +) -> DataArray | DataTree: + """Lazily read a segmentation mask as a ``Labels2DModel``. + + Parameters + ---------- + path + A TIFF or OME-TIFF with integer instance ids (0 = background). Of a file with several planes, the first is + read. + labels_models_kwargs + Keyword arguments passed to :class:`spatialdata.models.Labels2DModel`; their ``chunks`` also set the + spatial chunks of uncompressed files, which are memory-mapped. + transformations + The transformations of the region's image, by coordinate system. + + Returns + ------- + The labels element. + """ + _, labels_models_kwargs = _initialize_raster_models_kwargs({}, labels_models_kwargs) + _, defaults = _initialize_raster_models_kwargs({}, {}) + image, _, _ = read_tiff(path, spatial_chunks(labels_models_kwargs["chunks"], defaults["chunks"])) + if not np.issubdtype(image.dtype, np.integer): + raise ValueError(f"{path} has dtype {image.dtype}; a segmentation mask must have an integer dtype.") + return Labels2DModel.parse(image[0], dims=("y", "x"), transformations=dict(transformations), **labels_models_kwargs) + + +def read_region_labels( + files: Mapping[str, Path], + transformations: Mapping[str, Mapping[str, BaseTransformation]], + labels_models_kwargs: Mapping[str, Any], + image_shapes: Mapping[str, tuple[int, int]], +) -> dict[str, DataArray | DataTree]: + """Read the mask of each region, as resolved by :func:`resolve_per_region`, with the region's transformations. + + A mask whose ``(y, x)`` shape differs from ``image_shapes`` of its region is scaled onto the image if it is + downsampled by an integer factor, see :func:`mask_downsample_factor`: its transformations become + ``Sequence([Scale([f, f]), t])`` for each region transformation ``t``. Any other mismatch raises a warning: its + pixels would not align with the image. + + Returns + ------- + The labels elements by region id. + """ + if "*" in files: + raise ValueError( + f"A single labels mask {files['*']} was given for several regions; a mask belongs to one region image, " + "so pass a mapping {region: path}." + ) + labels = {} + for region, path in files.items(): + region_transformations = _transformations_of(transformations, region) + labels[region] = read_labels(path, labels_models_kwargs, region_transformations) + shape, image_shape = raster_shape(labels[region]), image_shapes.get(region) + if image_shape is None or shape == tuple(image_shape): + continue + factor = mask_downsample_factor(shape, image_shape) + if factor is None: + warnings.warn( + f"The mask {path} has shape {shape}, but the image of region '{region}' has shape " + f"{tuple(image_shape)} (y, x); the labels will not align with the image.", + UserWarning, + stacklevel=3, + ) + continue + logger.info( + "The mask %s has shape %s, the image of region '%s' %s (y, x); the labels are scaled by %d to align.", + path, + shape, + region, + tuple(image_shape), + factor, + ) + scale = Scale([factor, factor], axes=("y", "x")) + set_transformation( + labels[region], + { + cs: TransformationSequence([scale, transformation]) + for cs, transformation in region_transformations.items() + }, + set_all=True, + ) + return labels + + +def mask_downsample_factor(mask_shape: tuple[int, int], image_shape: tuple[int, int]) -> int | None: + """The integer factor ``f > 1`` by which a mask of ``(y, x)`` shape ``mask_shape`` downsamples an image. + + The mask must be ``image_shape / f`` along both axes, exactly or rounded up (each image size within ``f - 1`` + pixels of ``f`` times the mask size); `None` if there is no such factor. + """ + candidates = {-(-image // mask) for image, mask in zip(image_shape, mask_shape, strict=True) if mask > 0} + for factor in sorted(candidates): + if factor > 1 and all(-(-image // factor) == mask for image, mask in zip(image_shape, mask_shape, strict=True)): + return factor + return None + + +def _transformations_of( + transformations: Mapping[str, Mapping[str, BaseTransformation]], region: str +) -> dict[str, BaseTransformation]: + """A copy of the transformations of ``region``, ``Identity`` to ``"global"`` if it has none.""" + return dict(transformations.get(region, {"global": Identity()})) + + +def raster_shape(element: DataArray | DataTree) -> tuple[int, int]: + """The level-0 ``(y, x)`` shape of a raster element.""" + data = element["scale0"]["image"] if isinstance(element, DataTree) else element + return int(data.sizes["y"]), int(data.sizes["x"]) + + +def build_cells( + tables: Sequence[CellTable], + labels: Mapping[str, DataArray | DataTree], + transformations: Mapping[str, Mapping[str, BaseTransformation]], +) -> tuple[dict[str, ad.AnnData], dict[str, GeoDataFrame]]: + """Build the table elements and, for regions without labels, circle shapes from cell tables. + + Parameters + ---------- + tables + The cell tables, each covering one or more regions. + labels + The labels elements by region id; their regions' cells annotate ``_labels`` instead of the circles + ``_cells`` built here. + transformations + The transformations of each region, see :func:`region_transformations`; missing regions get ``Identity`` + to ``"global"``. + + Returns + ------- + The table elements, ``{"table": ...}`` if all tables have the same ``var`` apart from its ``source_column`` + (dropped if it differs), else ``{"_table": ...}`` per region, and the circle shapes by element name. The + circles have the cells' area (radius 1 px without a positive one). + """ + _check_disjoint_regions(tables) + row_ids = sorted( + region for table in tables if table.row_instance_ids for region in table.region.unique() if region in labels + ) + if row_ids: + warnings.warn( + f"The cell table of region(s) {row_ids} has no instance id column, so its row numbers are used as " + "instance ids; they may not match the values of the mask.", + UserWarning, + stacklevel=3, + ) + # tables with other markers are not concatenated; the source columns only record where a marker was read from + # (e.g. `cyc001_ch001:DAPI1` of a CSV, `cyc001_ch001` of an FCS) + panels = [_panel(table.var) for table in tables] + pieces = _split_by_region(tables) if any(not panel.equals(panels[0]) for panel in panels) else None + + def element(region: str) -> str: + suffix = PhenocyclerKeys.LABELS_SUFFIX if region in labels else PhenocyclerKeys.CELLS_SUFFIX + return f"{region}{suffix}" + + shapes = {} + for table in tables: + for region, rows in _region_rows(table).items(): + if region in labels: + continue + radius = np.ones(len(rows)) if table.radius_px is None else table.radius_px[rows] + if (no_area := ~(radius > 0)).any(): + logger.info( + "%d cells of region '%s' have no positive area; their circles get radius 1.", no_area.sum(), region + ) + radius = np.where(no_area, 1.0, radius) + shapes[element(region)] = ShapesModel.parse( + table.centroids_px[rows], + geometry=0, + radius=radius, + index=pd.Index(table.instance_id.to_numpy()[rows]), + transformations=_transformations_of(transformations, region), + ) + if shapes: + circle_regions = sorted(name.removesuffix(PhenocyclerKeys.CELLS_SUFFIX) for name in shapes) + warnings.warn( + f"No segmentation mask found for region(s) {circle_regions}; cells are represented as circles.", + UserWarning, + stacklevel=3, + ) + + if pieces is None: + return {PhenocyclerKeys.TABLE_KEY.value: _to_anndata(tables, element)}, shapes + return {f"{region}_table": _to_anndata([piece], element) for region, piece in sorted(pieces.items())}, shapes + + +def build_spatialdata( + scans: Sequence[AkoyaScan], + images: Mapping[str, da.Array], + channel_names: Mapping[str, list[str]], + cell_tables: Sequence[CellTable], + mask_files: Mapping[str, Path], + *, + table_paths: Mapping[str, str] | None, + image_models_kwargs: Mapping[str, Any], + labels_models_kwargs: Mapping[str, Any], + duplicate_scans: Mapping[str, AkoyaScan] | None = None, +) -> SpatialData: + """Assemble the regions of an Akoya read into a `SpatialData` with their metadata in ``attrs["akoya"]``. + + Each region gets the image ``_image``, the labels ``_labels`` from its mask (see + :func:`read_region_labels`) and the cells (see :func:`build_cells`), all in the coordinate systems of + :func:`region_transformations`; regions with cells or a mask but no image are unrelated to the scans. The cell + centroids are in image pixels, also when the mask is downsampled. + + Parameters + ---------- + scans + The scans read. + images + The lazy image of each scan, by region id, see :func:`build_image_element`. + channel_names + The channel names of each scan, by region id. + cell_tables + The cell tables read. + mask_files + The mask of each region, as resolved by :func:`resolve_per_region`. + table_paths + The paths of the cell tables, stored in ``attrs["akoya"]["cell_tables"]`` unless `None`. + image_models_kwargs + Keyword arguments passed to :class:`spatialdata.models.Image2DModel`. + labels_models_kwargs + Keyword arguments passed to :class:`spatialdata.models.Labels2DModel`. + duplicate_scans + Further scans of the same physical scan as a region (e.g. its 8-bit file), by that region's id. Each is + stored as its own image ``_image`` with exactly the transformations of the region; their + images and channel names are looked up in ``images`` and ``channel_names`` by ``duplicate.region``. + + Returns + ------- + The `SpatialData` object. + """ + cell_regions = {region for table in cell_tables for region in table.region.unique()} + transformations = region_transformations(scans, cell_regions | (mask_files.keys() - {"*"})) + image_elements = { + f"{scan.region}{PhenocyclerKeys.IMAGE_SUFFIX}": build_image_element( + images[scan.region], scan, channel_names[scan.region], transformations[scan.region], image_models_kwargs + ) + for scan in scans + } + for region, duplicate in (duplicate_scans or {}).items(): + image_elements[f"{duplicate.region}{PhenocyclerKeys.IMAGE_SUFFIX}"] = build_image_element( + images[duplicate.region], + duplicate, + channel_names[duplicate.region], + transformations[region], + image_models_kwargs, + ) + image_shapes = { + scan.region: raster_shape(image_elements[f"{scan.region}{PhenocyclerKeys.IMAGE_SUFFIX}"]) for scan in scans + } + labels = read_region_labels(mask_files, transformations, labels_models_kwargs, image_shapes) + tables, shapes = build_cells(cell_tables, labels, transformations) if cell_tables else ({}, {}) + sdata = SpatialData( + images=image_elements, + labels={f"{region}{PhenocyclerKeys.LABELS_SUFFIX}": element for region, element in labels.items()}, + shapes=shapes, + tables=tables, + ) + attrs = scans_to_attrs([*scans, *(duplicate_scans or {}).values()], channel_names) + if table_paths is not None: + attrs["cell_tables"] = dict(table_paths) + if mask_files: + attrs["labels"] = {region: str(file) for region, file in mask_files.items()} + sdata.attrs[PhenocyclerKeys.ATTRS_KEY.value] = attrs + return sdata + + +def _region_rows(table: CellTable) -> dict[str, np.ndarray]: + """The row indices of each region of ``table``, in order of first appearance.""" + codes, regions = pd.factorize(table.region) + return {str(region): np.flatnonzero(codes == i) for i, region in enumerate(regions)} + + +def _check_disjoint_regions(tables: Sequence[CellTable]) -> None: + seen: set[str] = set() + for table in tables: + regions = set(table.region.unique()) + if overlap := seen & regions: + raise ValueError(f"Several cell tables contain the cells of region(s) {sorted(overlap)}.") + seen |= regions + + +def _panel(var: pd.DataFrame) -> pd.DataFrame: + """The marker metadata of ``var`` without the provenance column ``source_column``.""" + return var.drop(columns="source_column", errors="ignore") + + +def _split_by_region(tables: Sequence[CellTable]) -> dict[str, CellTable]: + return {region: table.take(rows) for table in tables for region, rows in _region_rows(table).items()} + + +def _to_anndata(tables: Sequence[CellTable], element: Callable[[str], str]) -> ad.AnnData: + """Concatenate cell tables with the same markers into one table annotating the elements of their regions.""" + region_key, instance_key = CodexKeys.REGION_KEY.value, CodexKeys.INSTANCE_KEY.value + obs = pd.concat([table.obs for table in tables]) if len(tables) > 1 else tables[0].obs.copy() + regions = pd.concat([table.region for table in tables]) if len(tables) > 1 else tables[0].region + if not obs.index.is_unique: + obs.index = pd.Index(regions.to_numpy(dtype=object) + "_" + obs.index.to_numpy(dtype=object)) + X = np.concatenate([table.X for table in tables]) if len(tables) > 1 else tables[0].X + var = tables[0].var.copy() + if any(not table.var.equals(tables[0].var) for table in tables): + # the same markers read from columns of other names + var = _panel(var) + adata = ad.AnnData(X=X, obs=obs, var=var) + adata.obsm[CodexKeys.SPATIAL_KEY.value] = np.concatenate([table.centroids_px for table in tables]) + # spatialdata only writes alphanumeric column names and `_.-` + if any(sanitize_name(str(column), is_dataframe_column=True) != column for column in adata.obs.columns): + sanitize_table(adata, inplace=True) + elements = regions.map({region: element(region) for region in regions.unique()}).to_numpy() + adata.obs[region_key] = pd.Categorical(elements, categories=sorted(set(elements))) + adata.obs[instance_key] = np.concatenate([table.instance_id.to_numpy() for table in tables]) + return TableModel.parse( + adata, + region=list(adata.obs[region_key].cat.categories), + region_key=region_key, + instance_key=instance_key, + ) diff --git a/src/spatialdata_io/readers/_utils/_qptiff.py b/src/spatialdata_io/readers/_utils/_qptiff.py new file mode 100644 index 00000000..e2c2a2bd --- /dev/null +++ b/src/spatialdata_io/readers/_utils/_qptiff.py @@ -0,0 +1,512 @@ +"""PerkinElmer/Akoya QPTIFF metadata and lazy level-0 image access, shared by the Akoya readers.""" + +from __future__ import annotations + +import json +import re +import warnings +import xml.etree.ElementTree as ET +from collections import Counter +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any + +import dask.array as da +import tifffile +from spatialdata import sanitize_name + +from spatialdata_io._constants._constants import PhenocyclerKeys +from spatialdata_io.readers._utils._akoya import lazy_tiff_array, region_id + +if TYPE_CHECKING: + from collections.abc import Sequence + from pathlib import Path + +__all__ = [ + "AkoyaChannel", + "AkoyaScan", + "qpi_pixel_size_um", + "qptiff_base_name", + "qptiff_regions", + "qptiff_regions_with_duplicates", + "read_qptiff_image", + "read_qptiff_metadata", +] + +_FULL_RESOLUTION = "FullResolution" +_SCAN_SUFFIX = re.compile(r"_Scan\d+$") +_XML_DECLARATION = re.compile(r"^\s*<\?xml[^>]*\?>") +# micrometres per unit of the TIFF ResolutionUnit; XPosition/YPosition are expressed in the same unit +_UM_PER_RESOLUTION_UNIT = {tifffile.RESUNIT.CENTIMETER: 1e4, tifffile.RESUNIT.INCH: 25400.0} + + +@dataclass +class AkoyaChannel: + """One image channel of an Akoya acquisition. + + ``name`` is the marker identity before de-duplication (``Biomarker`` > ``Name`` > ``channel_{i}``). + ``cycle`` and ``channel`` (both 1-based) come from the scan profile and stay ``None`` when it does not + describe this channel. ``extra`` holds further JSON-serialisable metadata. + """ + + name: str + fluor: str | None = None + cycle: int | None = None + channel: int | None = None + exposure_ms: float | None = None + color: str | None = None + extra: dict[str, Any] = field(default_factory=dict) + + +@dataclass +class AkoyaScan: + """Metadata of one Akoya image (a QPTIFF scan or a legacy CODEX Processor region). + + ``generation`` is one of ``"fusion1"``, ``"fusion2"``, ``"cim_qptiff"``, ``"qpi_other"`` or + ``"legacy_processor"``. ``pixel_size_um`` is the calibrated pixel size, ``nominal_pixel_size_um`` the one + the acquisition software reports. ``stage_position_um`` is the ``(x, y)`` stage position of the upper-left + corner. ``raw`` keeps the unparsed source metadata (JSON-serialisable). ``pages`` are the indices of the + full-resolution pages of a QPTIFF, the image's channels. + """ + + path: Path + region: str + generation: str + software: str | None + slide_id: str | None + pixel_size_um: float | None + nominal_pixel_size_um: float | None + stage_position_um: tuple[float, float] | None + channels: list[AkoyaChannel] + rgb: bool + raw: dict[str, Any] + pages: list[int] = field(default_factory=list) + + +def read_qptiff_metadata(path: Path, region: str | None = None) -> AkoyaScan: + """Read the metadata of a PerkinElmer/Akoya QPTIFF. + + Channels are the ``FullResolution`` pages in file order; each page's own XML is authoritative for the + channel identity, while the page-0 ``ScanProfile`` (Fusion JSON or CIM ``ExperimentV4``) only adds cycle, + channel and reagent information. + + Parameters + ---------- + path + Path to the ``.qptiff`` file. + region + Region name; defaults to the file name without the ``.qptiff``/``.er.qptiff`` suffix. + + Returns + ------- + The scan metadata. + """ + with _open_qpi(path) as tf: + page0 = _page(tf, 0) + page_xml = [_parse_description(_page(tf, i).description) for i in range(len(tf.pages))] + full_res = [i for i, xml in enumerate(page_xml) if _text(xml, "ImageType") == _FULL_RESOLUTION] + if not full_res: + raise ValueError(f"{path} has no page with {_FULL_RESOLUTION}.") + rgb = len(full_res) == 1 and _page(tf, full_res[0]).samplesperpixel == 3 + pixel_size_um = qpi_pixel_size_um(page0) + stage_position_um = _stage_position_um(page0) + page0_description = page0.description + + root = page_xml[0] + software = _text(root, "AcquisitionSoftware") + experiment_v4 = root.find("ScanProfile/ExperimentV4") + scan_profile = _scan_profile(root) + channels = [_channel(page_xml[page_index], i, page_index) for i, page_index in enumerate(full_res)] + if isinstance(scan_profile, dict): + _enrich(channels, _fusion_profile_records(scan_profile)) + elif experiment_v4 is not None: + _enrich(channels, _experiment_v4_records(experiment_v4)) + + if pixel_size_um is None: + warnings.warn( + f"Could not determine the pixel size of {path}; pixel_size_um is None.", UserWarning, stacklevel=2 + ) + + return AkoyaScan( + path=path, + region=region if region is not None else qptiff_base_name(path), + generation=_generation(software, experiment_v4 is not None), + software=software, + slide_id=_text(root, "SlideID"), + pixel_size_um=pixel_size_um, + nominal_pixel_size_um=_nominal_pixel_size_um(root), + stage_position_um=stage_position_um, + channels=channels, + rgb=rgb, + raw={ + "page0_xml": page0_description, + "scan_profile": scan_profile, + "image_types": dict(Counter(_text(xml, "ImageType") or "unknown" for xml in page_xml)), + }, + pages=full_res, + ) + + +def read_qptiff_image(scan: AkoyaScan) -> da.Array: + """Lazily read the full-resolution image of a PerkinElmer/Akoya QPTIFF. + + The array keeps the native tile chunks of the file, opens the file only when its chunks are read and can be + pickled, e.g. for ``dask.distributed``. + + Parameters + ---------- + scan + The metadata of the ``.qptiff`` file, see :func:`read_qptiff_metadata`. + + Returns + ------- + ``(c, y, x)`` array, or ``(y, x, 3)`` for a brightfield RGB scan. + """ + with _open_qpi(scan.path) as tf: + page_shape = _page(tf, scan.pages[0]).shape + expected_shape = page_shape if scan.rgb else (len(scan.pages), *page_shape) + + image = lazy_tiff_array(scan.path, series=0) + if image.shape == expected_shape: + return image + # tifffile groups pages into series by shape heuristics; fall back to the pages the XML says are full-res + pages = [lazy_tiff_array(scan.path, page=i) for i in scan.pages] + return pages[0] if scan.rgb else da.stack(pages) + + +def qpi_pixel_size_um(page: tifffile.TiffPage) -> float | None: + """Return the pixel size of a QPTIFF page in micrometres, or ``None`` when it is unknown. + + The calibrated ``XResolution`` is used when ``ResolutionUnit`` is centimetre. Otherwise (e.g. the + meaningless 96 dpi written when the resolution is unknown) the nominal ``PixelSizeMicrons`` or the CIM + ``Resolution_nm`` from the page XML is returned, if present. + """ + unit = page.tags.get("ResolutionUnit") + resolution = page.tags.get("XResolution") + if unit is not None and resolution is not None and unit.value == tifffile.RESUNIT.CENTIMETER: + numerator, denominator = resolution.value + if numerator > 0 and denominator > 0: + return 1e4 * denominator / numerator + return _nominal_pixel_size_um(_parse_description(page.description)) + + +def _open_qpi(path: Path) -> tifffile.TiffFile: + if not path.is_file(): + raise FileNotFoundError(f"{path} does not exist or is not a file.") + tf = tifffile.TiffFile(path) + if not tf.is_qpi: + tf.close() + raise ValueError( + f"{path} is not a PerkinElmer/Akoya QPTIFF (Software tag does not start with 'PerkinElmer-QPI')." + ) + return tf + + +def _page(tf: tifffile.TiffFile, index: int) -> tifffile.TiffPage: + page = tf.pages[index] + if not isinstance(page, tifffile.TiffPage): # only frames when ``useframes`` is set, which we never do + raise TypeError(f"Expected a TiffPage at index {index}, got {type(page).__name__}.") + return page + + +def _parse_description(description: str) -> ET.Element: + # the declaration claims utf-16 but the tag is ASCII, which ElementTree rejects + return ET.fromstring(_XML_DECLARATION.sub("", description, count=1)) + + +def _text(element: ET.Element, tag: str) -> str | None: + """Stripped text of ``tag`` (a path relative to ``element``), ``None`` when missing or empty.""" + found = element.find(tag) + if found is None or found.text is None or not found.text.strip(): + return None + return found.text.strip() + + +def _float(text: str | None) -> float | None: + try: + return None if text is None else float(text) + except ValueError: + return None + + +def _bool(text: str | None) -> bool | None: + return None if text is None else text.lower() == "true" + + +def qptiff_base_name(path: Path) -> str: + """The file name without the ``.er.qptiff``/``.qptiff`` suffix (in any case).""" + for suffix in (PhenocyclerKeys.ER_QPTIFF_SUFFIX.value, PhenocyclerKeys.QPTIFF_SUFFIX.value): + if path.name.lower().endswith(suffix): + return path.name[: -len(suffix)] + return path.stem + + +def _is_extended_range(path: Path) -> bool: + return path.name.lower().endswith(PhenocyclerKeys.ER_QPTIFF_SUFFIX.value) + + +def qptiff_regions(files: Sequence[Path]) -> dict[str, Path]: + """The region id of every qptiff scan among ``files``, the same for both Akoya readers. + + Files of the same scan (folder and base name) are read once, preferring the extended-range ``.er.qptiff``. + The region id is the CODEX region id ``regRRR`` of a base name with ``reg`` and digits not preceded by a letter + (``exp_reg001`` to ``reg001``), else the base name without a trailing ``_Scan``, unless that makes two scans + share a region id. Raises if two scans have the same region id. + + Parameters + ---------- + files + The qptiff files, in the order of the result. + + Returns + ------- + The scan file of every region. + """ + return qptiff_regions_with_duplicates(files)[0] + + +def qptiff_regions_with_duplicates(files: Sequence[Path]) -> tuple[dict[str, Path], dict[str, Path]]: + """Like :func:`qptiff_regions`, and the 8-bit ``.qptiff`` of every region that also has an ``.er.qptiff``. + + Returns + ------- + The scan file of every region (the ``.er.qptiff`` if there is one), and the 8-bit file of the same scan by + region id for the regions that have both. + """ + scans: dict[tuple[Path, str], Path] = {} + eight_bit: dict[tuple[Path, str], Path] = {} + for file in files: + key = (file.parent, qptiff_base_name(file)) + if _is_extended_range(file): + if key in scans and not _is_extended_range(scans[key]): + eight_bit.setdefault(key, scans[key]) + scans[key] = file + elif key not in scans: + scans[key] = file + else: + eight_bit.setdefault(key, file) + # an 8-bit file is only a duplicate if the scan is read from an extended-range file + eight_bit = {key: file for key, file in eight_bit.items() if _is_extended_range(scans[key])} + bases = [base for _, base in scans] + # region ids must be valid element names; inForm names its files e.g. `_[13860,52919]_component_data.tif` + stripped = [sanitize_name(_SCAN_SUFFIX.sub("", base) or base) for base in bases] + # stripping `_Scan` must not merge two scans into one region + regions = [ + region_id(base) or (sanitize_name(base) if stripped.count(region) > 1 else region) + for base, region in zip(bases, stripped, strict=True) + ] + if len(set(regions)) < len(regions): + duplicated = sorted({region for region in regions if regions.count(region) > 1}) + found = [str(file) for file, region in zip(scans.values(), regions, strict=True) if region in duplicated] + raise ValueError( + f"Several qptiff files map to the same region(s) {duplicated}: {found}. Pass the path of a single qptiff " + "file or of a directory that contains each region once." + ) + region_of = dict(zip(scans, regions, strict=True)) + return dict(zip(regions, scans.values(), strict=True)), {region_of[key]: file for key, file in eight_bit.items()} + + +def _generation(software: str | None, has_experiment_v4: bool) -> str: + software = software or "" + if software.startswith("Fusion 1"): + return "fusion1" + if software.startswith("Fusion 2"): + return "fusion2" + if "instrument manager" in software.lower() or has_experiment_v4: + return "cim_qptiff" + return "qpi_other" + + +def _nominal_pixel_size_um(root: ET.Element) -> float | None: + # Vectra/Opal: ScanProfile/.../ScanResolution/PixelSizeMicrons; CIM: ScanProfile/ExperimentV4/Resolution_nm + pixel_size = _float(_text(root, ".//PixelSizeMicrons")) + if pixel_size is not None: + return pixel_size + resolution_nm = _float(_text(root, ".//Resolution_nm")) + return None if resolution_nm is None else resolution_nm / 1000 + + +def _stage_position_um(page: tifffile.TiffPage) -> tuple[float, float] | None: + unit = page.tags.get("ResolutionUnit") + um_per_unit = _UM_PER_RESOLUTION_UNIT.get(unit.value) if unit is not None else None + positions = [page.tags.get(name) for name in ("XPosition", "YPosition")] + if um_per_unit is None or positions[0] is None or positions[1] is None: + return None + (x_num, x_den), (y_num, y_den) = positions[0].value, positions[1].value + if x_den == 0 or y_den == 0: + return None + return x_num / x_den * um_per_unit, y_num / y_den * um_per_unit + + +def _scan_profile(root: ET.Element) -> dict[str, Any] | str | None: + """Page-0 ``ScanProfile``: parsed JSON (Fusion PhenoCycler mode) or the XML string (Opal, Vectra, CIM).""" + element = root.find("ScanProfile") + if element is None: + return None + if len(element): + return "".join(ET.tostring(child, encoding="unicode") for child in element) + text = (element.text or "").strip() + if not text: + return None + try: + profile = json.loads(text) + except json.JSONDecodeError: + return text + return profile if isinstance(profile, dict) else text + + +def _channel(xml: ET.Element, index: int, page_index: int) -> AkoyaChannel: + biomarker = _text(xml, "Biomarker") + fluor = _text(xml, "Name") + exposure_us = _float(_text(xml, "ExposureTime")) + extra: dict[str, Any] = {"page_index": page_index} + for key, tag in ( + ("is_unmixed_component", "IsUnmixedComponent"), + ("autofluorescence_subtracted", "AutofluorescenceSubtracted"), + ): + value = _bool(_text(xml, tag)) + if value is not None: + extra[key] = value + return AkoyaChannel( + name=biomarker or fluor or f"channel_{index}", + fluor=fluor if biomarker else None, + exposure_ms=None if exposure_us is None else exposure_us / 1000, + color=_text(xml, "Color"), + extra=extra, + ) + + +def _fusion_profile_records(profile: dict[str, Any]) -> list[tuple[str, dict[str, Any]]]: + """``(marker, fields)`` per item of the Fusion ``experimentDescription.wells``; one well per cycle. + + Wells and items that are not JSON objects, and ``wells``/``items`` that are not lists, are skipped. + """ + records = [] + description = profile.get("experimentDescription") + wells = description.get("wells") if isinstance(description, dict) else None + for cycle, well in enumerate(_objects(wells), start=1): + for channel, item in enumerate(_objects(well.get("items")), start=1): + marker = item.get("markerName") + if not isinstance(marker, str) or _is_blank(marker): + continue + # Fusion writes the string "None" for reagent fields that do not apply (e.g. DAPI) + extra = { + key: item[key] + for key in ("reporter", "barcode", "clone", "dye") + if item.get(key) not in (None, "", "None") + } + if item.get("channel"): + extra["channel_name"] = item["channel"] + extra["well"] = well.get("wellName") + exposure = item.get("exposure") + records.append( + ( + marker, + { + "cycle": cycle, + "channel": channel, + "exposure_ms": float(exposure) if isinstance(exposure, int | float) else None, + "extra": extra, + }, + ) + ) + return records + + +def _objects(value: Any) -> list[dict[str, Any]]: + """The entries of ``value`` if it is a list, else none; entries that are not objects become empty objects. + + They are kept as empty objects so that the cycle and channel indices still count every entry. + """ + return [entry if isinstance(entry, dict) else {} for entry in value] if isinstance(value, list) else [] + + +def _experiment_v4_records(experiment: ET.Element) -> list[tuple[str, dict[str, Any]]]: + """``(marker, fields)`` per ``Cycles/Cycle/Channels/Channel`` of a CIM ``ExperimentV4``.""" + records = [] + for cycle in experiment.iterfind("Cycles/Cycle"): + cycle_index = _float(_text(cycle, "Index")) + for channel in cycle.iterfind("Channels/Channel"): + marker = _text(channel, "MarkerName") + if marker is None or _is_blank(marker) or _text(channel, "MarkerClass") == "Blank": + continue + channel_index = _float(_text(channel, "Index")) + extra: dict[str, Any] = {} + for key, tag in (("dye", "DyeName"), ("marker_class", "MarkerClass")): + value = _text(channel, tag) + if value is not None: + extra[key] = value + wavelength = _float(_text(channel, "Wavelength_nm")) + if wavelength is not None: + extra["wavelength_nm"] = wavelength + records.append( + ( + marker, + { + "cycle": None if cycle_index is None else int(cycle_index), + "channel": None if channel_index is None else int(channel_index), + "exposure_ms": _float(_text(channel, "ExposureTime_ms")), + "extra": extra, + }, + ) + ) + return records + + +def _is_blank(marker: str) -> bool: + """Whether a scan-profile marker is a blank/empty slot (Fusion ``--``, CIM ``BLANK``/``EMPTY``).""" + return marker.strip() in {"", "--"} or marker.strip().casefold() in {"blank", "empty"} + + +def _name_key(name: str) -> str: + # CIM sanitises page Biomarker names (``F4/80`` -> ``F4-80``), so compare alphanumerics only (in any script) + return "".join(char for char in name.casefold() if char.isalnum()) + + +def _enrich(channels: list[AkoyaChannel], records: list[tuple[str, dict[str, Any]]]) -> None: + """Fill cycle/channel/extra of ``channels`` from the non-blank scan-profile ``records``. + + Record ``k`` enriches page ``k`` when there are as many records as pages and no pair contradicts the names. + Otherwise channels are joined on the normalised name with occurrence index, but only for names that occur + as often in the pages as in the profile; any other channel is left unchanged rather than guessed. The page + XML stays authoritative: names and order never change and the page exposure wins. + """ + for channel, (_, record) in _positional_pairs(channels, records) or _name_pairs(channels, records): + channel.cycle = record["cycle"] + channel.channel = record["channel"] + if channel.exposure_ms is None: + channel.exposure_ms = record["exposure_ms"] + channel.extra.update({k: v for k, v in record["extra"].items() if k not in channel.extra}) + + +def _positional_pairs( + channels: list[AkoyaChannel], records: list[tuple[str, dict[str, Any]]] +) -> list[tuple[AkoyaChannel, tuple[str, dict[str, Any]]]]: + if len(channels) != len(records): + return [] + # Fusion orders the pages of a cycle differently from the well items; a page whose marker sits at another + # position in the profile, or a record naming the marker of another page, means the orders differ and the + # positional pairing would be wrong + page_keys = {_name_key(channel.name) for channel in channels} + record_keys = {_name_key(name) for name, _ in records} + for channel, (name, _) in zip(channels, records, strict=True): + key, record_key = _name_key(channel.name), _name_key(name) + if key != record_key and (key in record_keys or record_key in page_keys): + return [] + return list(zip(channels, records, strict=True)) + + +def _name_pairs( + channels: list[AkoyaChannel], records: list[tuple[str, dict[str, Any]]] +) -> list[tuple[AkoyaChannel, tuple[str, dict[str, Any]]]]: + channels_by_key: dict[str, list[AkoyaChannel]] = {} + records_by_key: dict[str, list[tuple[str, dict[str, Any]]]] = {} + for channel in channels: + channels_by_key.setdefault(_name_key(channel.name), []).append(channel) + for record in records: + records_by_key.setdefault(_name_key(record[0]), []).append(record) + pairs: list[tuple[AkoyaChannel, tuple[str, dict[str, Any]]]] = [] + for key, keyed_channels in channels_by_key.items(): + keyed_records = records_by_key.get(key, []) + if len(keyed_records) == len(keyed_channels): + pairs.extend(zip(keyed_channels, keyed_records, strict=True)) + return pairs diff --git a/src/spatialdata_io/readers/_utils/_utils.py b/src/spatialdata_io/readers/_utils/_utils.py index 428ac900..49b29057 100644 --- a/src/spatialdata_io/readers/_utils/_utils.py +++ b/src/spatialdata_io/readers/_utils/_utils.py @@ -10,7 +10,7 @@ from spatialdata._logging import logger if TYPE_CHECKING: - from collections.abc import Mapping + from collections.abc import Iterator, Mapping from anndata import AnnData from ome_types.model import Pixels @@ -118,6 +118,24 @@ def _set_reader_metadata(sdata: SpatialData, reader: str) -> SpatialData: return sdata +def walk_visible(root: Path) -> Iterator[tuple[Path, list[str], list[str]]]: + """Like :func:`os.walk` of ``root`` in sorted order, without hidden files and folders, following symlinks. + + Hidden are names starting with ``.`` (also macOS AppleDouble ``._*`` files and ``.Trash-*`` folders) and + ``__MACOSX`` folders. A folder reached a second time through a symlink (e.g. a link to an ancestor) is skipped. + The directory names can be pruned in place, as with :func:`os.walk`. + """ + visited = set() + for dirpath, dirnames, filenames in os.walk(root, followlinks=True): + real = os.path.realpath(dirpath) + if real in visited: + dirnames.clear() + continue + visited.add(real) + dirnames[:] = sorted(name for name in dirnames if not name.startswith(".") and name != "__MACOSX") + yield Path(dirpath), dirnames, sorted(name for name in filenames if not name.startswith(".")) + + def parse_physical_size(path: Path | None = None, ome_pixels: Pixels | None = None) -> float: """Parse physical size from OME-TIFF to micrometer.""" from ome_types import from_tiff diff --git a/src/spatialdata_io/readers/codex.py b/src/spatialdata_io/readers/codex.py index 925ed077..3280e033 100644 --- a/src/spatialdata_io/readers/codex.py +++ b/src/spatialdata_io/readers/codex.py @@ -1,41 +1,156 @@ from __future__ import annotations -import os +import codecs import re +import warnings +from collections.abc import Collection, Iterable, Mapping +from dataclasses import dataclass from pathlib import Path from types import MappingProxyType from typing import TYPE_CHECKING, Any -import anndata as ad -import pandas as pd -import readfcs -from dask_image.imread import imread -from spatialdata import SpatialData -from spatialdata._logging import logger -from spatialdata.models import Image2DModel, ShapesModel, TableModel +import dask.array as da +import numpy as np +from spatialdata import SpatialData, sanitize_name from spatialdata_io._constants._constants import CodexKeys from spatialdata_io._docs import inject_docs -from spatialdata_io.readers._utils._utils import _set_reader_metadata +from spatialdata_io.readers._utils._akoya import ( + IMAGE_EXTENSION, + apply_channel_names, + normalize_regions, + read_json_object, + read_tiff, + region_id, + select_regions, + spatial_chunks, +) +from spatialdata_io.readers._utils._cell_tables import ( + CellTable, + build_spatialdata, + is_codex_table, + read_cell_table, + read_region_tables, + resolve_per_region, +) +from spatialdata_io.readers._utils._qptiff import ( + AkoyaChannel, + AkoyaScan, + qptiff_regions, + read_qptiff_image, + read_qptiff_metadata, +) +from spatialdata_io.readers._utils._utils import ( + _initialize_raster_models_kwargs, + _set_reader_metadata, + walk_visible, +) if TYPE_CHECKING: - from collections.abc import Mapping + from collections.abc import Sequence __all__ = ["codex"] +_PER_MARKER = re.compile(r"reg(\d+)_cyc(\d+)_ch(\d+)_(.+)\.tiff?$", re.IGNORECASE) +_MONTAGE = re.compile(r"reg\d+_montage\.tiff?$", re.IGNORECASE) +_TILE_STACK = re.compile(r"reg\d+_X\d+_Y\d+_Z\d+\.tiff?$", re.IGNORECASE) +# raw Keyence acquisition folders, e.g. `cyc001_reg001/1_00001_Z001_CH1.tif` +_RAW_CYCLE_DIR = re.compile(r"cyc\d+_reg\d+", re.IGNORECASE) +_TIFF_SUFFIXES = (".tif", ".tiff") +# the segmentation mask of a region, e.g. HuBMAP's `reg001_mask.ome.tiff` next to `reg001_expr.ome.tiff` +_MASK = re.compile(r"[_-](?:mask|labels)(?:\.ome)?\.tiff?$", re.IGNORECASE) +_TABLE_SUFFIXES = (CodexKeys.FCS_FILE.value, CodexKeys.FCS_FILE_CSV.value) +# deep enough for `/processed_/segm/segm-1/fcs/compensated/` +_MAX_DEPTH = 8 +# `channelNames.txt`/`experiment.json` sit in the experiment folder, up to `processed_/stitched/reg001/` above +_SIDECAR_ANCESTORS = 3 +_PATTERNS = ( + "'regNNN_cycCCC_chHHH_.tif(f)' (per-marker stitched images)", + "'regNNN_montage.tif', 'regNNN_X01_Y01_Z01.tif' or '*regNNN*.tif' (hyperstacks)", + f"'*{CodexKeys.QPTIFF}' (CODEX Instrument Manager)", + f"'*{CodexKeys.FCS_FILE}' or '*{CodexKeys.FCS_FILE_CSV}' (cell tables)", +) + + +@dataclass +class _LegacyLayout: + """Files of a CODEX Processor/Uploader output tree or of CODEX Instrument Manager scans. + + ``images`` maps each region to its per-marker images (sorted by cycle and channel), to its single hyperstack + or to its single ``.qptiff``. ``tables`` maps each region, or ``"*"`` for a table without a region in its + file name, to its ``.fcs``/``.csv`` cell table. ``masks`` maps regions to their segmentation mask. + """ + + images: dict[str, list[Path]] + tables: dict[str, Path] + masks: dict[str, Path] + channel_names_txt: Path | None + experiment_json: Path | None + @inject_docs(vx=CodexKeys) def codex( path: str | Path, fcs: bool = True, imread_kwargs: Mapping[str, Any] = MappingProxyType({}), + *, + regions: Sequence[str] | None = None, + channel_names: Sequence[str] | None = None, + cells_table: str | Path | Mapping[str, str | Path] | None = None, + labels: str | Path | Mapping[str, str | Path] | None = None, + image_models_kwargs: Mapping[str, Any] = MappingProxyType({}), + labels_models_kwargs: Mapping[str, Any] = MappingProxyType({}), ) -> SpatialData: - """Read *CODEX* formatted dataset. + """Read *CODEX* (CODEX Processor, CODEX Uploader and CODEX Instrument Manager) data. - This function reads the following files: + This function reads the following files, searching ``path`` recursively: - - ``_`{vx.FCS_FILE!r}```: Counts and metadata file. - - ``_`{vx.IMAGE_TIF!r}```: High resolution tif image. + - ``regNNN_cycCCC_chHHH_{vx.IMAGE_TIF}(f)``: per-marker stitched images of the CODEX Processor, + stacked lazily into one image per region in (cycle, channel) order; the channel names come from the file + names. + - ``regNNN_montage{vx.IMAGE_TIF}``, ``regNNN_X01_Y01_Z01{vx.IMAGE_TIF}`` or ``*regNNN*{vx.IMAGE_TIF}``: + hyperstacks (CODEX Uploader, older Processor versions). All axes except ``y`` and ``x`` are flattened + into channels, in file order. The channel names come from the ImageJ ``Labels``, else from + ``{vx.CHANNEL_NAMES_FILE!r}`` (one name per line, cycle-major), else they are ``cycCCC_chHHH``. + - ``*{vx.QPTIFF!r}``: scans of the CODEX Instrument Manager, read lazily with their channel metadata. + - ``{vx.EXPERIMENT_FILE!r}``: optional; the pixel size is taken from its ``xyResolution`` (nm). Like + ``{vx.CHANNEL_NAMES_FILE!r}``, it is found in ``path`` or up to three folders above it, ignoring the case + of the file name (e.g. ``channelnames.txt``). + - ``*{vx.FCS_FILE!r}`` / ``*{vx.FCS_FILE_CSV!r}``: cell tables of the CODEX Processor, one per region (or one + for all regions), e.g. ``segm/segm-1/fcs/compensated/regNNN_compensated.fcs``. Compensated tables are + preferred over uncompensated ones; other tables of the same region (e.g. of a second segmentation run) raise + an error. The ``cycCCC_chHHH`` marker columns are stored in ``X`` (dense float32) under their marker names, + taken from the FCS ``$PnS`` keywords or the CSV header ``cycCCC_chHHH:``; all other columns go to + ``obs``, the ``x``/``y`` pixel centroids to ``obsm["spatial"]``. A table without ``cycCCC_chHHH`` columns + reads its numeric columns that are not CODEX Processor columns (``x``, ``size``, ...) as the markers, with a + warning listing them. Other CSV files (without ``x``/``y`` and ``cell_id`` or ``cycCCC_chHHH`` columns, + e.g. antibody lists) and FCS files of other instruments are ignored. + + Raw acquisition and intermediate folders (``cycCCC_regNNN/``, ``tiles/`` and, next to ``stitched/``, ``bestFocus/``) + are ignored, as are hidden files and folders (names starting with ``.``, e.g. macOS ``._*`` files, and + ``__MACOSX/``); symlinked folders are followed. The region id ``regNNN`` is taken from the file names (``reg`` and + digits not preceded by a letter, e.g. ``exp_reg001`` but not ``Treg2``); a ``{vx.QPTIFF!r}`` without one is named as + :func:`spatialdata_io.phenocycler` names it: after its file, without the ``.er.qptiff`` or ``{vx.QPTIFF!r}`` suffix + and a trailing ``_Scan`` (unless that makes two regions share a name), and the extended-range ``.er.qptiff`` of a + scan is preferred, also for one with a region id (``exp_reg001.er.qptiff`` and ``exp_reg001.qptiff`` are region + ``reg001``). Several scans with the same region id raise an error. If a directory contains no such files, its top-level ``{vx.IMAGE_TIF}(f)`` files are read as + hyperstacks named after their files without the extension (``.ome.tif`` included), made valid element names by + :func:`spatialdata.sanitize_name` (e.g. ``Kidney sample.ome.tif`` becomes the region ``Kidney_sample``). Each region is stored as the image ``_image`` + and its metadata in ``sdata.attrs["akoya"]``. + + A single region is stored in the ``"global"`` coordinate system. Several CIM scans of the same sample space + (same slide, with stage positions and equal pixel sizes) are placed at their stage positions in ``"global"``; + otherwise (e.g. CODEX Processor regions, which have no stage positions) each region is stored in its own + coordinate system, named after the region. A region with cells or a mask but no image is unrelated to the + scans, so it also gets its own coordinate system. All elements of a region share its coordinate system. + + The cells of each region are stored as circles ``_cells`` with the cells' area (``size`` column), + annotated by one table ``"table"`` (``obs["region"]`` names the annotated element, ``obs["cell_id"]`` is the + instance key), or by one table ``_table`` per region if the regions have different markers. Columns of + the cell table that spatialdata cannot write are sanitised. A region with a segmentation mask (see ``labels``) + gets the labels ``_labels`` instead of circles, and its cells annotate them. The paths of the tables and + masks read are stored in ``sdata.attrs["akoya"]["cell_tables"]`` and ``["labels"]``, by region id (``"*"`` for a + file of several regions). .. seealso:: @@ -44,60 +159,464 @@ def codex( Parameters ---------- path - Path to the directory containing the data. + Path to a CODEX output directory, or to a single image (``{vx.QPTIFF!r}`` or ``{vx.IMAGE_TIF}(f)``). fcs - Whether a .fcs file is provided. If False, a .csv file is expected. + Whether to prefer ``{vx.FCS_FILE!r}`` cell tables over ``{vx.FCS_FILE_CSV!r}`` ones when both exist. With + `False`, a region that only has a ``{vx.FCS_FILE!r}`` table is read from it, with a warning. imread_kwargs - Keyword arguments passed to :func:`dask_image.imread.imread`. + Ignored, with a warning: every image is read lazily with ``tifffile``. Kept for backward compatibility. + regions + Regions to read, e.g. ``["reg001"]``, including only their cells. If `None`, all regions are read. + channel_names + Channel names to use instead of the discovered ones; must have one entry per channel of every region read. + cells_table + A CODEX cell table to read instead of the discovered ones (explicit beats discovered): a path, or a mapping + from region id to path, e.g. ``{{"reg001": "reg001_cells.csv"}}``. A single path for several regions needs a + ``region`` column. + labels + A segmentation mask (TIFF/OME-TIFF with integer instance ids) per region: a path for a single region, or a + mapping from region id to path; a mask inside ``path`` is not read as an image. If `None`, a TIFF named after + a region image or with its region id, ending in ``_mask`` or ``_labels`` (e.g. ``reg001_mask.ome.tiff``, + ``Kidney_labels.tif`` next to ``Kidney.tif``, or ``Tonsil_mask.tif`` or ``Tonsil_Scan1_mask.tif`` next to + ``Tonsil_Scan1.qptiff``), is the mask of that region; such files are never read as images. + Each mask is stored lazily as + the labels ``_labels`` with the transformation of the region's image; the cells of the region annotate + it instead of circles, and their instance ids must be the mask values. A mask whose shape is the region image's + divided by an integer factor ``f`` (exactly, or rounded up) is scaled onto the image (``Scale([f, f])`` before + the region's transformation); the cell centroids stay in image pixels. Any other shape mismatch raises a + warning. Of a mask file with several planes, the first is read. + image_models_kwargs + Keyword arguments passed to :class:`spatialdata.models.Image2DModel`. + labels_models_kwargs + Keyword arguments passed to :class:`spatialdata.models.Labels2DModel`. Returns ------- :class:`spatialdata.SpatialData` """ + if imread_kwargs: + warnings.warn( + "`imread_kwargs` is ignored: every image is read lazily with tifffile.", UserWarning, stacklevel=2 + ) + regions = normalize_regions(regions) + channel_names = [channel_names] if isinstance(channel_names, str) else channel_names path = Path(path) - patt = re.compile(".*.fcs") if fcs else re.compile(".*.csv") - path_files = [i for i in os.listdir(path) if patt.match(i)] - if path_files and CodexKeys.FCS_FILE or CodexKeys.FCS_FILE_CSV in patt.pattern: - df = ( - readfcs.ReadFCS(path / path_files[0]).data - if CodexKeys.FCS_FILE in path_files[0] - else pd.read_csv(path_files[0], header=0, index_col=None) + # explicit beats discovered: a given table disables table discovery, and a given mask is no image + masks = labels.values() if isinstance(labels, Mapping) else [] if labels is None else [labels] + layout = _discover_legacy( + path, + fcs, + find_tables=cells_table is None, + find_masks=labels is None, + masks={Path(mask).resolve() for mask in masks}, + ) + # the regions of a table without a region in its file name are known only once it is read + regionless = _read_regionless_table(layout.tables["*"], layout) if "*" in layout.tables else None + table_regions = set(regionless.region.unique()) if regionless is not None else set() + available = sorted(layout.images.keys() | layout.tables.keys() - {"*"} | table_regions) + selected = select_regions(available, regions) + _warn_fcs_fallback([file for key, file in layout.tables.items() if key == "*" or key in selected], fcs) + table_files = resolve_per_region(cells_table, selected, available) + if labels is None: + mask_files = {region: file for region, file in layout.masks.items() if region in selected} + else: + mask_files = resolve_per_region(labels, selected, available) + experiment = read_json_object(layout.experiment_json) if layout.experiment_json is not None else None + names_txt = _read_channel_names_txt(layout.channel_names_txt) if layout.channel_names_txt is not None else None + + scans: list[AkoyaScan] = [] + arrays: dict[str, da.Array] = {} + for region in (region for region in selected if region in layout.images): + files = layout.images[region] + if _is_qptiff(files[0]): + scans.append(read_qptiff_metadata(files[0], region=region)) + arrays[region] = read_qptiff_image(scans[-1]) + continue + arrays[region], channels = _read_legacy_image( + files, layout.channel_names_txt, image_models_kwargs=image_models_kwargs + ) + scans.append( + AkoyaScan( + path=files[0] if len(files) == 1 else files[0].parent, + region=region, + generation="legacy_processor", + software=None, + slide_id=None, + pixel_size_um=_pixel_size_um(experiment), + nominal_pixel_size_um=None, + stage_position_um=None, + channels=channels, + rgb=False, + raw={ + "files": [str(file) for file in files], + "experiment_json": experiment, + "channel_names_txt": names_txt, + }, + ) ) + + names = {scan.region: apply_channel_names(scan, channel_names) for scan in scans} + if cells_table is not None: + cell_tables = read_region_tables(table_files, available=list(layout.images), selected=selected) + table_paths = {key: str(file) for key, file in table_files.items()} else: - raise ValueError("Cannot determine data set. Expecting a file with format .fcs or .csv") - - adata = _codex_df_to_anndata(df) - assert isinstance(adata.obs, pd.DataFrame) - - xy = adata.obsm[CodexKeys.SPATIAL_KEY] - shapes = ShapesModel.parse(xy, geometry=0, radius=1, index=adata.obs[CodexKeys.INSTANCE_KEY]) - region = adata.obs[CodexKeys.REGION_KEY].unique()[0] - adata.obs[CodexKeys.REGION_KEY] = adata.obs[CodexKeys.REGION_KEY].astype("category") - table = TableModel.parse(adata, region=region, region_key=CodexKeys.REGION_KEY, instance_key=CodexKeys.INSTANCE_KEY) - - im_patt = re.compile(".*.tif") - path_files = [i for i in os.listdir(path) if im_patt.match(i)] - if path_files and CodexKeys.IMAGE_TIF in path_files[0]: - image = imread(path_files[0], **imread_kwargs) - images = { - "images": Image2DModel.parse( - image, - scale_factors=[2, 2], - rgb=None, + named_files = {key: layout.tables[key] for key in selected if key in layout.tables} + named = read_region_tables(named_files, available=available, selected=selected) + read = dict(zip(named_files, named, strict=True)) + if regionless is not None: + regionless = regionless.take(np.flatnonzero(regionless.region.isin(selected).to_numpy())) + if len(regionless.region): + read = {"*": regionless, **read} + cell_tables = list(read.values()) + table_paths = {} + for key, table in read.items(): + # by region id, like a given table; `"*"` only for a region-less table of several regions + values = table.region.unique() + region = values[0] if key == "*" and len(values) == 1 and values[0] not in read else key + table_paths[region] = str(layout.tables[key]) + sdata = build_spatialdata( + scans, + arrays, + names, + cell_tables, + mask_files, + table_paths=table_paths, + image_models_kwargs=image_models_kwargs, + labels_models_kwargs=labels_models_kwargs, + ) + return _set_reader_metadata(sdata, "codex") + + +def _warn_fcs_fallback(files: Iterable[Path], fcs: bool) -> None: + fallback = [str(file) for file in files if file.suffix.lower() == CodexKeys.FCS_FILE.value] + if not fcs and fallback: + warnings.warn( + f"`fcs=False`, but no CSV cell table was found for {fallback}; reading the FCS table(s) instead.", + UserWarning, + stacklevel=3, + ) + + +def _read_regionless_table(file: Path, layout: _LegacyLayout) -> CellTable: + """Read the discovered cell table ``file`` that has no region in its file name.""" + only_image = next(iter(layout.images)) if len(layout.images) == 1 else None + # next to a single image named after its file (e.g. `Kidney.tif`), a single region value names that image + table = read_cell_table( + file, default_region=only_image, regions=layout.images.keys() | layout.tables.keys() - {"*"} + ) + # values that name no region of the dataset are regions without an image + table.region = table.region.map({value: sanitize_name(value) for value in table.region.unique()}) + return table + + +def _discover_legacy( + path: Path, + fcs: bool, + *, + find_tables: bool = True, + find_masks: bool = True, + masks: Collection[Path] = frozenset(), +) -> _LegacyLayout: + """Find the images, cell tables and sidecar files of a CODEX dataset. + + Parameters + ---------- + path + A CODEX output directory, searched recursively in sorted order, or a single image or table file. + fcs + Whether to prefer ``.fcs`` over ``.csv`` tables of the same region. + find_tables + Whether to discover cell tables. + find_masks + Whether to discover the segmentation masks of the regions. + masks + The resolved paths of the segmentation masks, which are not images. + + Returns + ------- + The discovered files. + """ + if not path.exists(): + raise FileNotFoundError(f"{path} does not exist.") + files = [path] if path.is_file() else _walk(path) + if masks: + files = [file for file in files if file.resolve() not in masks] + # masks named after their region are never images + named_masks = [_is_tiff(file) and _MASK.search(file.name) is not None for file in files] + image_files = [file for file, is_mask in zip(files, named_masks, strict=True) if not is_mask] + mask_files = [file for file, is_mask in zip(files, named_masks, strict=True) if is_mask] + images = _discover_images(image_files) + if not images: + # the original (2019) layout: one stitched `.tif` next to its table, without a region in the file name + top_level = image_files if path.is_file() else [file for file in image_files if file.parent == path] + images = {sanitize_name(IMAGE_EXTENSION.sub("", file.name)): [file] for file in top_level if _is_tiff(file)} + tables = _discover_tables(files, fcs) if find_tables else {} + if not images and not tables: + raise FileNotFoundError(f"No CODEX images or cell tables found in {path}. Expected: {'; '.join(_PATTERNS)}.") + return _LegacyLayout( + images=dict(sorted(images.items())), + tables=dict(sorted(tables.items())), + masks=_discover_masks(mask_files, images) if find_masks else {}, + channel_names_txt=_find_sidecar(path, files, CodexKeys.CHANNEL_NAMES_FILE.value), + experiment_json=_find_sidecar(path, files, CodexKeys.EXPERIMENT_FILE.value), + ) + + +def _read_legacy_image( + files: list[Path], + channel_names_txt: Path | None, + image_models_kwargs: Mapping[str, Any] = MappingProxyType({}), +) -> tuple[da.Array, list[AkoyaChannel]]: + """Lazily read the image of one legacy CODEX region. + + Parameters + ---------- + files + The per-marker images of the region, or its single hyperstack, as returned by :func:`_discover_legacy`. + channel_names_txt + Optional ``channelNames.txt`` naming the hyperstack channels when the file has no ImageJ ``Labels``. + image_models_kwargs + Keyword arguments for :class:`spatialdata.models.Image2DModel`; their ``chunks`` set the spatial chunks of + uncompressed files, which are memory-mapped. + + Returns + ------- + The ``(c, y, x)`` image and its channels. + """ + chunks = _spatial_chunks(image_models_kwargs) + markers = [(match, file) for file in files if (match := _PER_MARKER.search(file.name)) is not None] + if len(markers) == len(files): + return _read_per_marker(markers, chunks) + if len(files) != 1: + raise ValueError(f"Expected per-marker images or a single hyperstack, got {[str(f) for f in files]}.") + + image, leading_shape, labels = read_tiff(files[0], chunks) + n_planes = image.shape[0] + # the last non-spatial axis holds the channels of a cycle; ImageJ axis labels are unreliable (cycles as Z) + n_channels = leading_shape[-1] if leading_shape else 1 + names: list[str] | None = labels if labels is not None and len(labels) == n_planes else None + if names is None and channel_names_txt is not None: + names = _read_channel_names_txt(channel_names_txt) + if len(names) != n_planes: + warnings.warn( + f"{channel_names_txt} lists {len(names)} channel names, but {files[0]} has {n_planes} channels; " + "using 'cycCCC_chHHH' names instead.", + UserWarning, + stacklevel=2, + ) + names = None + channels = [] + for i in range(n_planes): + cycle, channel = i // n_channels + 1, i % n_channels + 1 + name = names[i] if names is not None else f"cyc{cycle:03d}_ch{channel:03d}" + channels.append(AkoyaChannel(name=name, cycle=cycle, channel=channel)) + return image, channels + + +def _read_per_marker( + markers: list[tuple[re.Match[str], Path]], chunks: tuple[int, int] +) -> tuple[da.Array, list[AkoyaChannel]]: + parsed = sorted((int(match.group(2)), int(match.group(3)), match.group(4), file) for match, file in markers) + + planes = [] + for _, _, _, file in parsed: + plane, _, _ = read_tiff(file, chunks) + if plane.shape[0] != 1: + raise ValueError(f"{file} holds {plane.shape[0]} planes; expected a single 2-D marker image.") + planes.append(plane[0]) + if len({(plane.shape, plane.dtype) for plane in planes}) > 1: + found = {file.name: (plane.shape, str(plane.dtype)) for (*_, file), plane in zip(parsed, planes, strict=True)} + raise ValueError(f"The per-marker images of one region must share shape and dtype, got {found}.") + channels = [ + AkoyaChannel(name=marker, cycle=cycle, channel=channel, extra={"file": file.name}) + for cycle, channel, marker, file in parsed + ] + return da.stack(planes), channels + + +def _spatial_chunks(image_models_kwargs: Mapping[str, Any]) -> tuple[int, int]: + """The ``(y, x)`` chunk size requested in ``image_models_kwargs``, with the reader default when unspecified.""" + default, _ = _initialize_raster_models_kwargs({}, {}) + chunks = _initialize_raster_models_kwargs(image_models_kwargs, {})[0]["chunks"] + return spatial_chunks(chunks, default["chunks"][-2:]) + + +def _walk(root: Path) -> list[Path]: + """All visible files below ``root`` in sorted order, skipping raw acquisition and intermediate processing folders.""" + files: list[Path] = [] + for directory, dirnames, filenames in walk_visible(root): + depth = len(directory.relative_to(root).parts) + has_stitched = any(name.casefold() == "stitched" for name in dirnames) + dirnames[:] = [name for name in dirnames if depth < _MAX_DEPTH and not _is_ignored_dir(name, has_stitched)] + files.extend(directory / name for name in filenames) + return files + + +def _is_ignored_dir(name: str, has_stitched: bool) -> bool: + folded = name.casefold() + return folded == "tiles" or (folded == "bestfocus" and has_stitched) or _RAW_CYCLE_DIR.fullmatch(name) is not None + + +def _discover_images(files: Sequence[Path]) -> dict[str, list[Path]]: + per_marker: dict[str, dict[tuple[int, int], Path]] = {} + stacks: dict[str, list[Path]] = {} + all_qptiffs = [] + for file in files: + region = region_id(file.name) + if _is_qptiff(file): + all_qptiffs.append(file) + elif _is_tiff(file) and region is not None: + if (match := _PER_MARKER.search(file.name)) is None: + stacks.setdefault(region, []).append(file) + continue + markers = per_marker.setdefault(region, {}) + key = (int(match.group(2)), int(match.group(3))) + if key in markers: + raise ValueError( + f"Found several images for cycle {key[0]}, channel {key[1]} of region '{region}': " + f"{markers[key]} and {file}. Pass the path of a directory that contains each region once." + ) + markers[key] = file + + # named like `phenocycler()` names them, e.g. `exp_reg001.qptiff` -> `reg001`, `Tonsil_Scan1.er.qptiff` -> `Tonsil` + qptiffs = qptiff_regions(all_qptiffs) + + images: dict[str, list[Path]] = {} + # the Processor's stitched outputs come first: cell tables refer to their pixel space + for region in sorted(per_marker.keys() | stacks.keys() | qptiffs.keys()): + if region in per_marker: + images[region] = [per_marker[region][key] for key in sorted(per_marker[region])] + elif region in stacks: + images[region] = [_single(region, _preferred_stacks(stacks[region]))] + else: + images[region] = [qptiffs[region]] + return images + + +def _discover_masks(files: Sequence[Path], images: Mapping[str, Sequence[Path]]) -> dict[str, Path]: + """The masks among ``files`` (``_mask.tif`` or ``_labels``, also ``.ome.tiff``) of the regions of ``images``. + + A mask belongs to the region in its file name (``regNNN``), else to the region it is named after: the region id + (e.g. ``Tonsil_mask.tif``) or the file name of its image (``Tonsil_Scan1_mask.tif`` of ``Tonsil_Scan1.qptiff``). + """ + by_image = {region: region for region in images} | { + IMAGE_EXTENSION.sub("", file.name): region for region, paths in images.items() for file in paths + } + masks: dict[str, list[Path]] = {} + for file in files: + region = region_id(file.name) or by_image.get(_MASK.sub("", file.name)) + if region is None or region not in images: + warnings.warn( + f"Ignoring the mask {file}: it is named after no region image. Pass it as `labels` to read it.", + UserWarning, + stacklevel=4, + ) + continue + masks.setdefault(region, []).append(file) + for region, paths in masks.items(): + if len(paths) > 1: + raise ValueError( + f"Found several masks for region '{region}': {[str(p) for p in paths]}. Pass the mask to read as " + "`labels`." ) - } - sdata = SpatialData(images=images, shapes={str(region): shapes}, tables={"table": table}) + return {region: paths[0] for region, paths in sorted(masks.items())} + + +def _preferred_stacks(files: list[Path]) -> list[Path]: + for pattern in (_MONTAGE, _TILE_STACK): + matching = [file for file in files if pattern.search(file.name)] + if matching: + return matching + return files + + +def _single(region: str, files: list[Path]) -> Path: + if len(files) > 1: + raise ValueError( + f"Found several images for region '{region}': {[str(f) for f in files]}. Pass the path of a single " + "image or of a directory that contains each region once." + ) + return files[0] + + +def _discover_tables(files: Sequence[Path], fcs: bool) -> dict[str, Path]: + preferred = CodexKeys.FCS_FILE.value if fcs else CodexKeys.FCS_FILE_CSV.value + candidates: dict[str, list[Path]] = {} + for file in files: + if file.suffix.lower() in _TABLE_SUFFIXES and is_codex_table(file): + candidates.setdefault(region_id(file.name) or "*", []).append(file) + + def uncompensated(file: Path) -> bool: + return any("uncompensated" in part.casefold() for part in file.parts) + + tables = {} + for region, paths in candidates.items(): + best = min(map(uncompensated, paths)) + # the same table may be written as both .fcs and .csv: the preferred format is read + same_table: dict[Path, list[Path]] = {} + for path in sorted(paths): + if uncompensated(path) == best: + same_table.setdefault(path.with_suffix(""), []).append(path) + if len(same_table) > 1: + where = "without a region in their file name" if region == "*" else f"of region '{region}'" + raise ValueError( + f"Found several cell tables {where}: {[str(p) for group in same_table.values() for p in group]}; the reader cannot " + "choose between them. Pass the path of a directory that contains only one, or pass the table as " + "`cells_table`." + ) + (formats,) = same_table.values() + tables[region] = min(formats, key=lambda path: path.suffix.lower() != preferred) + return tables + + +def _find_sidecar(path: Path, files: Sequence[Path], name: str) -> Path | None: + """The shallowest file called ``name`` in ``files``, else in the nearest ancestor of ``path`` holding one. + + Names are compared case-insensitively: SPACEc and the HuBMAP codex-pipeline write ``channelnames.txt``. + """ + folded = name.casefold() + matches = sorted( + (file for file in files if file.name.casefold() == folded), key=lambda file: (len(file.parts), str(file)) + ) + if matches: + return matches[0] + # absolute, so that the ancestors of e.g. `codex(".")` are searched too + start = (path.parent if path.is_file() else path).absolute() + for directory in [start, *start.parents][: _SIDECAR_ANCESTORS + 1]: + try: + found = sorted(file for file in directory.iterdir() if file.name.casefold() == folded and file.is_file()) + except OSError: # e.g. an unreadable ancestor directory + continue + if found: + return found[0] + return None + + +def _read_channel_names_txt(path: Path) -> list[str]: + data = path.read_bytes() + # written on Windows: UTF-16 with a byte order mark, UTF-8 (with or without one) or the ANSI code page + if data.startswith((codecs.BOM_UTF16_LE, codecs.BOM_UTF16_BE)): + text = data.decode("utf-16") else: - logger.warning("Cannot find .tif file. Will build spatialdata with shapes and table only.") - sdata = SpatialData(shapes={str(region): shapes}, tables={"table": table}) + try: + text = data.decode("utf-8-sig") + except UnicodeDecodeError: + text = data.decode("cp1252", errors="replace") + return [line.strip() for line in text.splitlines() if line.strip()] - return _set_reader_metadata(sdata, "codex") + +def _pixel_size_um(experiment: Mapping[str, Any] | None) -> float | None: + """Pixel size in µm from the ``experiment.json`` of the CODEX Processor (nm), ``None`` when it is missing.""" + for key in ("xyResolution", "per_pixel_XY_resolution", "Resolution_nm"): + value = (experiment or {}).get(key) + if isinstance(value, int | float) and not isinstance(value, bool) and value > 0: + return value / 1000 + return None + + +def _is_qptiff(path: Path) -> bool: + return path.suffix.lower() == CodexKeys.QPTIFF.value -def _codex_df_to_anndata(df: pd.DataFrame) -> ad.AnnData: - """Convert a codex formatted .fcs dataframe or .csv file to anndata.""" - adata = ad.AnnData(df.filter(regex="cyc.*")) - adata.obs = df[df.columns.drop(list(df.filter(regex="cyc.*")))] - adata.obsm[CodexKeys.SPATIAL_KEY] = df[["x", "y"]].values - adata.var_names_make_unique() - return adata +def _is_tiff(path: Path) -> bool: + return path.suffix.lower() in _TIFF_SUFFIXES diff --git a/src/spatialdata_io/readers/phenocycler.py b/src/spatialdata_io/readers/phenocycler.py new file mode 100644 index 00000000..fdae290a --- /dev/null +++ b/src/spatialdata_io/readers/phenocycler.py @@ -0,0 +1,286 @@ +from __future__ import annotations + +import csv +import re +import warnings +from pathlib import Path +from types import MappingProxyType +from typing import TYPE_CHECKING, Any + +from spatialdata import SpatialData +from spatialdata._logging import logger + +from spatialdata_io._constants._constants import PhenocyclerKeys +from spatialdata_io._docs import inject_docs +from spatialdata_io.readers._utils._akoya import ( + apply_channel_names, + normalize_regions, + read_json_object, + select_regions, +) +from spatialdata_io.readers._utils._cell_tables import build_spatialdata, resolve_per_region +from spatialdata_io.readers._utils._qptiff import ( + _fusion_profile_records, + _is_blank, + _is_extended_range, + qptiff_base_name, + qptiff_regions_with_duplicates, + read_qptiff_image, + read_qptiff_metadata, +) +from spatialdata_io.readers._utils._utils import _set_reader_metadata, walk_visible + +if TYPE_CHECKING: + from collections.abc import Mapping, Sequence + + from spatialdata_io.readers._utils._qptiff import AkoyaScan + +__all__ = ["phenocycler"] + +# the Fusion writes per-cycle acquisition files to ``raw/images/phenocycler/`` (HuBMAP layout) +_ACQUISITION_DIR = ("raw", "images", "phenocycler") +# region id suffix of the 8-bit scan read next to its extended-range scan (``read_all_bit_depths``) +_EIGHT_BIT_SUFFIX = "_8bit" + + +@inject_docs(vx=PhenocyclerKeys) +def phenocycler( + path: str | Path, + *, + regions: Sequence[str] | None = None, + channel_names: Sequence[str] | None = None, + labels: str | Path | Mapping[str, str | Path] | None = None, + read_all_bit_depths: bool = False, + image_models_kwargs: Mapping[str, Any] = MappingProxyType({}), + labels_models_kwargs: Mapping[str, Any] = MappingProxyType({}), +) -> SpatialData: + """Read *Akoya PhenoCycler-Fusion* 1.x/2.x data, and PerkinElmer/Akoya Vectra and PhenoImager QPTIFFs. + + This function reads the following files: + + - ``*{vx.QPTIFF_SUFFIX!r}`` / ``*{vx.ER_QPTIFF_SUFFIX!r}``: the scans, one region each. The images are read + lazily; when both ``X{vx.QPTIFF_SUFFIX}`` and ``X{vx.ER_QPTIFF_SUFFIX}`` exist, the extended-range + (higher bit depth) ``{vx.ER_QPTIFF_SUFFIX!r}`` is used, unless ``read_all_bit_depths`` is set. + - ``*{vx.XPD_SUFFIX!r}``: optional experiment protocol (JSON) next to a scan. + - ``*{vx.CHANNELS_CSV_SUFFIX!r}``: optional channel list next to a scan or, in the HuBMAP layout, in + ``lab_processed/images/``. + + The same QPTIFF container is written by Vectra Polaris and PhenoImager HT (also unmixed, uint16) scans, whose + channels are read like the Fusion's, by brightfield (H&E) scans, read as the three channels ``r``, ``g`` and + ``b``, and by inForm component data (float32, usually a ``.tif`` file, which is read if ``path`` is that single + file). The image keeps the dtype of the file. These files have the generation ``"qpi_other"`` in the metadata. + + Sidecar files are matched to a scan by file-name prefix and stored as metadata only: the channel names always come + from the per-page XML of the qptiff. A sidecar whose number of markers differs from the number of channels raises a + warning. Acquisition files (``*.raw.qptiff``, ``*.qptiff.raw``, ``*.qptiff.intermediate`` and anything in a + ``raw/images/phenocycler/`` folder below ``path``) are ignored, as are hidden files and folders (names starting with + ``.``, e.g. macOS ``._*`` files, and ``__MACOSX/``); symlinked folders are followed. File suffixes are matched + ignoring case. + + The region name is the CODEX region id ``regNNN`` of a file name with ``reg`` and digits not preceded by a letter + (``exp_reg001.qptiff`` is ``reg001``, as in :func:`spatialdata_io.codex`; several scans with the same region id + raise an error), else the file name without the ``{vx.ER_QPTIFF_SUFFIX!r}``/``{vx.QPTIFF_SUFFIX!r}`` suffix and + without a trailing ``_Scan``, unless that makes two regions share a name, made a valid element name by + :func:`spatialdata.sanitize_name` (characters other than letters, digits, ``_``, ``.`` and ``-`` are replaced by + ``_``, e.g. ``A_[1,2]`` becomes ``A__1_2_``). Each region is stored as the image + ``{vx.IMAGE_SUFFIX}``. A single region is stored in the ``"global"`` coordinate system. Several regions + of the same sample space (scanned on the same slide, with stage positions and equal pixel sizes) are placed at + their stage positions in ``"global"``; otherwise the regions are unrelated and each is stored in its own + coordinate system, named after the region. All elements of a region share its coordinate system. The + per-region metadata is stored in ``sdata.attrs[{vx.ATTRS_KEY!r}]``. + + A segmentation mask given in ``labels`` is stored lazily as the labels ``{vx.LABELS_SUFFIX}`` with the + transformation of the region's image; the mask paths are stored in ``sdata.attrs[{vx.ATTRS_KEY!r}]["labels"]``, + by region id. + + .. seealso:: + + - `Akoya PhenoCycler-Fusion `_. + + Parameters + ---------- + path + Path to a ``{vx.QPTIFF_SUFFIX!r}`` file (or an inForm ``.tif`` QPI file), or to a directory that is searched + recursively for ``{vx.QPTIFF_SUFFIX!r}`` files. + regions + Regions to read. If `None`, all regions are read. Region ids, like the keys of the ``labels`` mapping, may + also be given as raw names, which are sanitised like the region ids (e.g. ``"Slide 1"`` selects the region + ``Slide_1``). + channel_names + Channel names to use instead of the ones stored in the qptiff; must have one entry per channel of every + region read. + labels + A segmentation mask (TIFF/OME-TIFF with integer instance ids) per region: a path for a single region, or a + mapping from region id to path. A mask whose shape is the region image's divided by an integer factor ``f`` + (exactly, or rounded up) is scaled onto the image (``Scale([f, f])`` before the region's transformation); + any other shape mismatch raises a warning. Of a mask file with several planes, the first is read. + read_all_bit_depths + If `True`, a scan that exists as both ``X{vx.QPTIFF_SUFFIX}`` (8-bit) and ``X{vx.ER_QPTIFF_SUFFIX}`` + (extended range, e.g. uint16) is read from both files, also if ``path`` is the extended-range file: the + extended-range file as region ``X`` and the 8-bit file as the image ``X_8bit{vx.IMAGE_SUFFIX}`` (region + ``X_8bit`` in ``sdata.attrs[{vx.ATTRS_KEY!r}]``). Both are the same physical scan, so the 8-bit image has + exactly the transformations of region ``X`` (a warning is raised if their shapes differ). Masks belong to + region ``X`` only. If `False`, only the extended-range file is read. + image_models_kwargs + Keyword arguments passed to :class:`spatialdata.models.Image2DModel`. + labels_models_kwargs + Keyword arguments passed to :class:`spatialdata.models.Labels2DModel`. + + Returns + ------- + :class:`spatialdata.SpatialData` + """ + regions = normalize_regions(regions) + channel_names = [channel_names] if isinstance(channel_names, str) else channel_names + path = Path(path) + qptiffs, eight_bit = _discover_qptiffs(path) + selected = select_regions(list(qptiffs), regions) + mask_files = resolve_per_region(labels, selected, list(qptiffs)) + + scans = [read_qptiff_metadata(qptiffs[region], region=region) for region in selected] + duplicates: dict[str, AkoyaScan] = {} + if read_all_bit_depths: + for region in selected: + if region in eight_bit: + duplicates[region] = read_qptiff_metadata(eight_bit[region], region=f"{region}{_EIGHT_BIT_SUFFIX}") + if clash := sorted({scan.region for scan in duplicates.values()} & set(qptiffs)): + raise ValueError( + f"Cannot read the 8-bit scans as region(s) {clash}: a scan with that region id exists. Rename it or " + "pass read_all_bit_depths=False." + ) + names = {} + for scan in [*scans, *duplicates.values()]: + names[scan.region] = apply_channel_names(scan, channel_names) + _read_sidecars(scan) + images = {scan.region: read_qptiff_image(scan) for scan in [*scans, *duplicates.values()]} + for region, duplicate in duplicates.items(): + if images[region].shape[-2:] != images[duplicate.region].shape[-2:]: + warnings.warn( + f"The 8-bit scan {duplicate.path} has shape {images[duplicate.region].shape[-2:]} (y, x), the " + f"extended-range scan {qptiffs[region]} {images[region].shape[-2:]}; they may not be the same scan, " + "but the 8-bit image is stored with its transformations.", + UserWarning, + stacklevel=2, + ) + sdata = build_spatialdata( + scans, + images, + names, + [], + mask_files, + table_paths=None, + image_models_kwargs=image_models_kwargs, + labels_models_kwargs=labels_models_kwargs, + duplicate_scans=duplicates, + ) + return _set_reader_metadata(sdata, "phenocycler") + + +def _discover_qptiffs(path: Path) -> tuple[dict[str, Path], dict[str, Path]]: + """Find the qptiff scans of a PhenoCycler-Fusion dataset: a ``.qptiff`` file, or a directory searched recursively. + + Returns the scan file of every region, sorted by path, and the 8-bit file of the regions that also have an + extended-range one; see :func:`qptiff_regions_with_duplicates`. + """ + if not path.exists(): + raise FileNotFoundError(f"{path} does not exist.") + if path.is_file(): + files = [path] + if _is_extended_range(path): + # its 8-bit file, read with `read_all_bit_depths` + base = qptiff_base_name(path) + files.append(path.with_name(base + path.name[len(base) + len(".er") :])) + return qptiff_regions_with_duplicates([file for file in files if file.is_file()]) + files = sorted( + directory / name + for directory, _, names in walk_visible(path) + for name in names + if name.lower().endswith(PhenocyclerKeys.QPTIFF_SUFFIX) + and (directory / name).is_file() + and not _is_acquisition_file((directory / name).relative_to(path)) + ) + if not files: + raise FileNotFoundError( + f"No .qptiff files found in {path} (searched recursively for '*.qptiff', excluding acquisition files)." + ) + return qptiff_regions_with_duplicates(files) + + +def _is_acquisition_file(path: Path) -> bool: + """Whether ``path``, relative to the searched directory, is a per-cycle acquisition file of the Fusion.""" + if path.name.lower().endswith(".raw" + PhenocyclerKeys.QPTIFF_SUFFIX): + return True + parts = tuple(part.casefold() for part in path.parent.parts) + n = len(_ACQUISITION_DIR) + return any(parts[i : i + n] == _ACQUISITION_DIR for i in range(len(parts) - n + 1)) + + +def _read_sidecars(scan: AkoyaScan) -> None: + """Store the ``.xpd`` and ``*channels.csv`` sidecars of ``scan`` in ``scan.raw`` and check their channel count. + + The sidecars are metadata only, so one that cannot be read is ignored with a warning. + """ + base = qptiff_base_name(scan.path) + xpd = _find_sidecar(scan.path, base, PhenocyclerKeys.XPD_SUFFIX) + protocol = read_json_object(xpd) if xpd is not None else None + if xpd is not None and protocol is not None: + scan.raw["xpd"] = protocol + # an .xpd has the layout of the Fusion scan profile's `experimentDescription` + markers = [marker for marker, _ in _fusion_profile_records({"experimentDescription": protocol})] + if markers: + _check_channel_count(scan, xpd, markers) + elif protocol.get("wells") is not None: + warnings.warn( + f"{xpd.name} has no readable markers in its 'wells', so its channel count is not checked.", + UserWarning, + stacklevel=3, + ) + + channels_csv = _find_sidecar(scan.path, base, PhenocyclerKeys.CHANNELS_CSV_SUFFIX) + if channels_csv is not None: + try: + with channels_csv.open(encoding="utf-8-sig", newline="") as f: + rows = list(csv.DictReader(f)) + except (UnicodeDecodeError, csv.Error) as e: + warnings.warn(f"Ignoring {channels_csv}: it is not a UTF-8 CSV file ({e}).", UserWarning, stacklevel=3) + return + scan.raw["channels_csv"] = rows + header = rows[0] if rows else {} + column = next((key for key in header if isinstance(key, str) and key.strip().casefold() == "marker"), None) + if column is not None: + markers = [row[column] for row in rows if row[column] and not _is_blank(row[column])] + _check_channel_count(scan, channels_csv, markers) + + +def _find_sidecar(qptiff: Path, base: str, suffix: str) -> Path | None: + """First file ending in ``suffix`` whose name starts with ``base``, next to the scan or in HuBMAP processed dir.""" + pattern = re.compile(rf"{re.escape(base)}(?![0-9A-Za-z]).*{re.escape(suffix)}$") + qptiff = qptiff.absolute() # also look above the working directory, e.g. for `phenocycler(".")` + directories = [qptiff.parent] + # HuBMAP keeps `raw/images/.qptiff` and `lab_processed/images/.qptiff.channels.csv` + for ancestor in qptiff.parents: + if (ancestor / "raw").is_dir() and qptiff.is_relative_to(ancestor / "raw"): + directories.append(ancestor / "lab_processed" / "images") + break + for directory in directories: + if not directory.is_dir(): + continue + matches = sorted(p for p in directory.iterdir() if p.is_file() and pattern.match(p.name)) + if matches: + if len(matches) > 1: + logger.info("Found several sidecar files %s for %s; using %s.", matches, qptiff, matches[0]) + return matches[0] + return None + + +def _check_channel_count(scan: AkoyaScan, sidecar: Path, markers: list[str]) -> None: + # the Fusion acquires DAPI in every cycle but stores it as one page, so each marker is counted once + n_markers = len(set(markers)) + n_channels = len({channel.name for channel in scan.channels}) + if n_markers != n_channels: + warnings.warn( + f"{sidecar.name} lists {n_markers} distinct markers, but {scan.path.name} has {n_channels} distinct " + f"channel names ({len(scan.channels)} channels); using the channel names from the qptiff.", + UserWarning, + stacklevel=2, + ) diff --git a/tests/_akoya_data.py b/tests/_akoya_data.py new file mode 100644 index 00000000..b826ce0b --- /dev/null +++ b/tests/_akoya_data.py @@ -0,0 +1,233 @@ +"""Synthetic Akoya (PhenoCycler-Fusion qptiff, CODEX processor) inputs for the reader tests.""" + +from __future__ import annotations + +import json +from typing import TYPE_CHECKING, Any, Literal + +import numpy as np +import pandas as pd +import tifffile + +from spatialdata_io.readers._utils import _akoya + +if TYPE_CHECKING: + from pathlib import Path + + import pytest + + +def fusion_scan_profile_json(wells: list[list[tuple[str, str]]]) -> str: + """A Fusion ``ScanProfile`` JSON; each well is a list of ``(markerName, channel)``.""" + items = [[{"markerName": marker, "channel": channel, "dye": channel} for marker, channel in well] for well in wells] + wells_json = [{"wellName": f"A{i + 1}", "items": well} for i, well in enumerate(items)] + return json.dumps({"experimentDescription": {"formatVersion": 2, "wells": wells_json}}) + + +def cim_experiment_v4_xml(cycles: list[list[str]]) -> str: + """A CODEX Instrument Manager ``ExperimentV4`` document; ``cycles`` lists the marker names per channel.""" + cycle_xml = "".join( + f"{i}" + + "".join( + f"{j}{marker}" + f"{'Blank' if marker == 'BLANK' else 'Other'}" + f"100{358 + 100 * j}" + for j, marker in enumerate(markers, start=1) + ) + + "" + for i, markers in enumerate(cycles, start=1) + ) + return f"377.39{cycle_xml}" + + +def write_qptiff( + path: Path, + *, + biomarkers: list[str | None], + names: list[str] | None = None, + shape: tuple[int, int] = (256, 320), + dtype: str = "uint8", + levels: int = 2, + software: str = "Fusion 2.3.1", + scan_profile: str | None = None, + resolution_px_per_cm: float | None = 1e4 / 0.5069, + resolution_unit: str = "CENTIMETER", + stage_position_cm: tuple[float, float] | None = None, + label_page: bool = True, + rgb: bool = False, + bigtiff: bool = True, + tile: tuple[int, int] | None = (64, 64), + slide_id: str = "SLIDE1", +) -> np.ndarray: + """Write a synthetic QPTIFF and return the level-0 data, ``(c, y, x)`` or ``(y, x, 3)`` when ``rgb``. + + Page order mirrors real files: full-resolution pages, an RGB thumbnail, the reduced-resolution groups (the last + one stripped), an RGB overview and an RGB label. + """ + rng = np.random.default_rng(0) + + def random(size: tuple[int, ...]) -> np.ndarray: + if np.dtype(dtype).kind == "f": # unmixed component data + return (rng.random(size) * 1000).astype(dtype) + return rng.integers(0, np.iinfo(dtype).max, size, dtype=dtype) + + n_pages = 1 if rgb else len(biomarkers) + names = names or [b if b is not None else f"Channel {i + 1}" for i, b in enumerate(biomarkers)] + resolution = None if resolution_px_per_cm is None else (resolution_px_per_cm, resolution_px_per_cm) + + def description(image_type: str, name: str, biomarker: str | None = None, profile: str | None = None) -> str: + fields = { + "AcquisitionSoftware": software, + "ImageType": image_type, + "SlideID": slide_id, + "ExposureTime": "450000", + "Name": name, + "Color": "255,255,0", + "Biomarker": biomarker, + "ScanProfile": profile, + } + tags = "".join(f"<{tag}>{text}\r\n" for tag, text in fields.items() if text is not None) + root = "PerkinElmer-QPI-ImageDescription" + return f'\r\n<{root}>\r\n{tags}' + + def write_page(tw: tifffile.TiffWriter, data: np.ndarray, desc: str, subfiletype: int = 0, **kwargs: Any) -> None: + tw.write( + data, + software="PerkinElmer-QPI", + description=desc, + resolution=resolution, + resolutionunit=resolution_unit if resolution is not None else None, + subfiletype=subfiletype, + photometric="rgb" if data.ndim == 3 else "minisblack", + metadata=None, + **kwargs, + ) + + def rgb_page() -> np.ndarray: + return rng.integers(0, 256, (64, 64, 3), dtype=np.uint8) + + extratags: list[Any] = [] + for code, value in zip((286, 287), stage_position_cm or (), strict=False): + extratags.append((code, 5, 1, (round(value * 1e6), 1_000_000), False)) + + data = random((*shape, 3) if rgb else (n_pages, *shape)) + with tifffile.TiffWriter(path, bigtiff=bigtiff) as tw: + for i in range(n_pages): + desc = description("FullResolution", names[i], biomarkers[i], scan_profile if i == 0 else None) + write_page(tw, data if rgb else data[i], desc, tile=tile, extratags=extratags if i == 0 else None) + write_page(tw, rgb_page(), description("Thumbnail", "Thumbnail")) + for level in range(1, levels): + side = (shape[0] // 2**level, shape[1] // 2**level) + for i in range(n_pages): + page = random((*side, 3) if rgb else side) + desc = description("ReducedResolution", names[i], biomarkers[i]) + write_page(tw, page, desc, subfiletype=1, tile=None if level == levels - 1 else tile) + write_page(tw, rgb_page(), description("Overview", "Overview")) + if label_page: + write_page(tw, rgb_page(), description("Label", "Label")) + return data + + +def write_fcs(path: Path, df: pd.DataFrame, marker_names: dict[str, str] | None = None) -> None: + """Write a minimal FCS 3.0 file: little-endian float32, ``$PnN`` = column name, ``$PnS`` = marker name.""" + marker_names = marker_names or {} + data = np.ascontiguousarray(df.to_numpy(dtype=" bytes: + pairs = [("$BEGINDATA", f"{data_begin:>12}"), ("$ENDDATA", f"{data_end:>12}"), ("$PAR", str(df.shape[1]))] + pairs += [("$TOT", str(len(df))), ("$DATATYPE", "F"), ("$BYTEORD", "1,2,3,4"), ("$MODE", "L")] + pairs.append(("$NEXTDATA", "0")) + for i, col in enumerate(df.columns, start=1): + pairs += [(f"$P{i}N", str(col)), (f"$P{i}B", "32"), (f"$P{i}E", "0,0"), (f"$P{i}R", "262144")] + if col in marker_names: + pairs.append((f"$P{i}S", marker_names[col])) + return ("/" + "/".join(f"{k}/{v}" for k, v in pairs) + "/").encode("ascii") + + text_start = 58 + text_end = text_start + len(text_segment(0, 0)) - 1 # the offsets are fixed-width + data_begin, data_end = text_end + 1, text_end + len(data) + header = f"FCS3.0 {text_start:>8}{text_end:>8}{data_begin:>8}{data_end:>8}{0:>8}{0:>8}".encode("ascii") + path.write_bytes(header.ljust(text_start) + text_segment(data_begin, data_end) + data) + + +def write_legacy_processor_tree( + root: Path, + *, + regions: int = 1, + cycles: int = 2, + channels: int = 4, + shape: tuple[int, int] = (128, 160), + layout: Literal["per_marker", "hyperstack", "old_flat"] = "per_marker", + fcs: bool = True, + csv: bool = True, + n_cells: int = 50, +) -> dict[str, Any]: + """Write a CODEX Processor output tree and return ``{"markers", "images", "cells"}``. + + ``images`` maps region name to ``(n_markers, y, x)`` for ``per_marker`` and ``old_flat`` and to + ``(cycles, channels, y, x)`` for ``hyperstack``. ``cells`` is the table written for the last region. + """ + n_markers = cycles * channels + markers = ["DAPI" if i % channels == 0 else f"MARKER{i}" for i in range(n_markers)] + rng = np.random.default_rng(0) + root.mkdir(parents=True, exist_ok=True) + processed = root / "processed_2020-01-01" + channel_cols = [f"cyc{c + 1:03d}_ch{ch + 1:03d}" for c in range(cycles) for ch in range(channels)] + images: dict[str, np.ndarray] = {} + cells = pd.DataFrame() + + if layout != "old_flat": + (root / "channelNames.txt").write_text("\n".join(markers) + "\n") + (root / "experiment.json").write_text(json.dumps({"numCycles": cycles, "numChannels": channels})) + + for r in range(1, regions + 1): + name = f"reg{r:03d}" + stack = rng.integers(0, 2**16, (n_markers, *shape), dtype=np.uint16) + cells = pd.DataFrame( + { + "cell_id": np.arange(1, n_cells + 1, dtype=float), + "region": float(r), + "x": rng.uniform(0, shape[1], n_cells), + "y": rng.uniform(0, shape[0], n_cells), + "size": rng.uniform(10, 100, n_cells), + **{col: rng.uniform(0, 1000, n_cells) for col in channel_cols}, + } + ) + table_dir = processed / "segm" / "segm-1" / "fcs" / "compensated" + table_stem = f"{name}_compensated" + if layout == "per_marker": + images[name] = stack + image_dir = processed / "stitched" / name + image_dir.mkdir(parents=True, exist_ok=True) + for i, marker in enumerate(markers): + cyc, ch = divmod(i, channels) + tifffile.imwrite(image_dir / f"{name}_cyc{cyc + 1:03d}_ch{ch + 1:03d}_{marker}.tif", stack[i]) + elif layout == "hyperstack": + images[name] = stack.reshape(cycles, channels, *shape) + (processed / "stitched").mkdir(parents=True, exist_ok=True) + tifffile.imwrite( + processed / "stitched" / f"{name}_montage.tif", images[name], imagej=True, metadata={"axes": "TCYX"} + ) + else: + images[name] = stack + tifffile.imwrite(root / f"exp_{name}.tif", stack[0]) + table_dir, table_stem = root, f"exp_{name}" + table_dir.mkdir(parents=True, exist_ok=True) + if fcs: + write_fcs(table_dir / f"{table_stem}.fcs", cells, dict(zip(channel_cols, markers, strict=True))) + if csv and not (layout == "old_flat" and fcs): + cells.to_csv(table_dir / f"{table_stem}.csv", index=False) + return {"markers": markers, "images": images, "cells": cells} + + +def forbid_pixel_reads(monkeypatch: pytest.MonkeyPatch) -> None: + """Make every read of pixel data fail: the chunk reads of the lazy tiff arrays and tifffile's own decoding.""" + + def fail(*args: object, **kwargs: object) -> None: + raise AssertionError("pixel data read eagerly") + + monkeypatch.setattr(_akoya._TiffZarrSource, "__getitem__", fail) + monkeypatch.setattr(_akoya._ContiguousTiff, "__getitem__", fail) + monkeypatch.setattr(tifffile.TiffPage, "decode", property(lambda self: fail)) + monkeypatch.setattr(tifffile.TiffPage, "asarray", fail) + monkeypatch.setattr(tifffile.TiffFile, "asarray", fail) diff --git a/tests/readers/test_akoya_utils.py b/tests/readers/test_akoya_utils.py new file mode 100644 index 00000000..64da6dd3 --- /dev/null +++ b/tests/readers/test_akoya_utils.py @@ -0,0 +1,127 @@ +from __future__ import annotations + +import os +import pickle +from pathlib import Path + +import dask +import dask.array as da +import numpy as np +import pytest +import tifffile +from spatialdata.transformations import Identity, Translation + +from spatialdata_io.readers._utils._akoya import ( + apply_channel_names, + dedupe_channel_names, + read_tiff, + region_transformations, +) +from spatialdata_io.readers._utils._qptiff import AkoyaChannel, AkoyaScan + + +def _scan( + region: str, + stage: tuple[float, float] | None = None, + pixel_size_um: float | None = 0.5, + slide_id: str | None = "slide", +) -> AkoyaScan: + channels = [AkoyaChannel(name=name) for name in ("DAPI", "DAPI", "CD3")] + return AkoyaScan( + Path(f"{region}.qptiff"), region, "fusion2", None, slide_id, pixel_size_um, None, stage, channels, False, {} + ) + + +def test_channel_names() -> None: + # a generated suffix never repeats a name already in the list + assert dedupe_channel_names(["DAPI", "DAPI", "DAPI (2)"]) == ["DAPI", "DAPI (3)", "DAPI (2)"] + scan = _scan("a") + assert apply_channel_names(scan, None) == ["DAPI", "DAPI (2)", "CD3"] + assert apply_channel_names(scan, ["a", "b", "a"]) == ["a", "b", "a (2)"] + with pytest.raises(ValueError, match="expected 3"): + apply_channel_names(scan, ["a", "b"]) + # the keys of a mapping or the order of a set would silently become the channel names + with pytest.raises(TypeError, match="sequence"): + apply_channel_names(scan, {"a", "b", "c"}) # type: ignore[arg-type] + + +_GLOBAL = { + "a": {"global": Translation([0.0, 0.0], axes=("x", "y"))}, + "b": {"global": Translation([200.0, 100.0], axes=("x", "y"))}, +} +_OWN = {"a": {"a": Identity()}, "b": {"b": Identity()}} + + +@pytest.mark.parametrize( + ("scans", "expected"), + [ + ([_scan("a", (0.0, 0.0)), _scan("b", (100.0, 50.0))], _GLOBAL), + # pixel sizes equal up to rounding share the sample space + ([_scan("a", (0.0, 0.0)), _scan("b", (100.0, 50.0), pixel_size_um=0.5004)], None), + ([_scan("a", (0.0, 0.0), slide_id="s1"), _scan("b", (100.0, 50.0), slide_id="s2")], _OWN), + ([_scan("a", (0.0, 0.0)), _scan("b", (1.0, 1.0), pixel_size_um=0.25)], _OWN), + ([_scan("a"), _scan("b")], _OWN), + ([_scan("a", (0.0, 0.0), pixel_size_um=None), _scan("b", (1.0, 1.0))], _OWN), + ([_scan("a", (0.0, 0.0), slide_id=None), _scan("b", (1.0, 1.0), slide_id=None)], _OWN), + ], + ids=[ + "same_slide", + "rounded_pixel_size", + "other_slide", + "other_pixel_size", + "no_stage", + "no_pixel_size", + "no_slide", + ], +) +def test_region_transformations(scans: list[AkoyaScan], expected: dict | None) -> None: + result = region_transformations(scans) + if expected is None: + assert set(result["b"]) == {"global"} + else: + assert result == expected + # a region without an image is unrelated to the scans + assert region_transformations(scans, ["c"])["c"] == {"c": Identity()} + assert region_transformations([_scan("a", (5.0, 5.0))]) == {"a": {"global": Identity()}} + + +@pytest.mark.parametrize("compression", [None, "zlib"]) +def test_read_tiff_is_picklable_and_keys_depend_on_chunks(tmp_path: Path, compression: str | None) -> None: + # compressed tiles are read through tifffile's zarr store, uncompressed data is memory-mapped + path = tmp_path / "x.tif" + data = np.arange(2 * 64 * 96, dtype=np.uint16).reshape(2, 64, 96) + tifffile.imwrite(path, data, compression=compression, tile=(16, 16) if compression else None) + + small, _, _ = read_tiff(path, (16, 16)) + large, _, _ = read_tiff(path, (32, 32)) + + np.testing.assert_array_equal(pickle.loads(pickle.dumps(small)).compute(), data) + # arrays of one file with other chunks must not share dask keys, or computing both mixes their chunks + assert small.name != large.name + np.testing.assert_array_equal(da.stack([small.rechunk(large.chunks), large]).compute(), [data, data]) + + +@pytest.mark.parametrize("imagej", [False, True]) +def test_read_tiff_big_endian(tmp_path: Path, imagej: bool) -> None: + # ImageJ/Fiji write big-endian tiffs by default; the memory-mapped path must use the file's byte order + data = np.array([[[1, 2], [300, 0]]], np.uint16) + tifffile.imwrite(tmp_path / "x.tif", data, byteorder=">", imagej=imagej) + + np.testing.assert_array_equal(read_tiff(tmp_path / "x.tif", (16, 16))[0].compute(), data) + + +@pytest.mark.skipif(not os.path.isdir("/dev/fd"), reason="needs /dev/fd to count open files") +@pytest.mark.parametrize("compression", [None, "zlib"]) +def test_read_tiff_does_not_keep_files_open(tmp_path: Path, compression: str | None) -> None: + # computed chunks and arrays kept alive must not hold a file descriptor each (EMFILE on large reads) + data = np.arange(64 * 96, dtype=np.uint16).reshape(64, 96) + images = [] + for i in range(20): + tifffile.imwrite(tmp_path / f"{i}.tif", data, compression=compression, tile=(16, 16) if compression else None) + images.append(read_tiff(tmp_path / f"{i}.tif", (16, 16))[0]) + n_open = len(os.listdir("/dev/fd")) + + chunks = dask.compute(*[chunk for image in images for chunk in image.to_delayed().ravel()]) + + assert len(os.listdir("/dev/fd")) - n_open <= 2 + np.testing.assert_array_equal(chunks[-1][0], data[48:, 80:]) diff --git a/tests/readers/test_cell_tables.py b/tests/readers/test_cell_tables.py new file mode 100644 index 00000000..185996e6 --- /dev/null +++ b/tests/readers/test_cell_tables.py @@ -0,0 +1,159 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING, Any + +import numpy as np +import pandas as pd +import pytest +import tifffile +from spatialdata import get_extent, read_zarr +from spatialdata.transformations import Identity, Scale, Sequence, Translation, get_transformation +from tests._akoya_data import forbid_pixel_reads, write_fcs, write_qptiff + +from spatialdata_io import phenocycler +from spatialdata_io.readers._utils._cell_tables import read_cell_table, read_labels + +if TYPE_CHECKING: + from pathlib import Path + + +def _cells(columns: tuple[str, ...]) -> pd.DataFrame: + rng = np.random.default_rng(0) + meta = {"cell_id": [1.0, 2.0], "region": 1.0, "x": [1.0, 2.0], "y": [3.0, 4.0], "size": [3.0, 12.0]} + return pd.DataFrame(meta | {column: rng.uniform(0, 1000, 2) for column in columns}) + + +@pytest.mark.parametrize( + ("pns", "names"), + [ + ({"cyc001_ch001": "DAPI", "cyc001_ch002": "CD3", "cyc002_ch001": "DAPI"}, ["DAPI", "CD3", "DAPI (2)"]), + ({"cyc001_ch001": " ", "cyc001_ch002": "CD3"}, ["cyc001_ch001", "CD3", "cyc002_ch001"]), + (None, ["cyc001_ch001", "cyc001_ch002", "cyc002_ch001"]), + ], + ids=["pns", "blank_pns", "no_pns"], +) +def test_fcs_marker_names(tmp_path: Path, pns: dict[str, str] | None, names: list[str]) -> None: + columns = ("cyc001_ch001", "cyc001_ch002", "cyc002_ch001") + cells = _cells(columns) + write_fcs(tmp_path / "reg001.fcs", cells, pns) + + table = read_cell_table(tmp_path / "reg001.fcs", default_region=None) + + assert list(table.var.index) == names + assert list(table.var["cycle"]) == [1, 1, 2] + assert list(table.var["source_column"]) == list(columns) + np.testing.assert_allclose(table.X, cells[list(columns)], rtol=1e-6) + assert list(table.region) == ["reg001", "reg001"] + assert list(table.instance_id) == [1, 2] + np.testing.assert_allclose(table.centroids_px, cells[["x", "y"]]) + assert table.radius_px is not None + np.testing.assert_allclose(table.radius_px, np.sqrt(cells["size"] / np.pi)) + + +def test_csv_header_names_the_markers(tmp_path: Path) -> None: + path = tmp_path / "reg002.csv" + path.write_text( + "cell_id:cell_id,region:region,x:x,y:y,size:size,cyc001_ch001:DAPI1,cyc001_ch002:CD3\n1,2,10,20,30,1,7\n" + ) + + table = read_cell_table(path, default_region=None) + + assert list(table.var.index) == ["DAPI1", "CD3"] + assert list(table.var["source_column"]) == ["cyc001_ch001:DAPI1", "cyc001_ch002:CD3"] + # the instance ids and regions are stored under the keys of the table element + assert list(table.obs.columns) == ["x", "y", "size"] + assert list(table.region) == ["reg002"] + + +@pytest.mark.parametrize("mangled", [False, True]) +def test_csv_without_cycle_columns_reads_numeric_columns( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch, mangled: bool +) -> None: + if mangled: # pandas >= 3.1 renames repeated names (CD3.1) with the pyarrow engine too + read_csv = pd.read_csv + + def mangling_read_csv(*args: Any, **kwargs: Any) -> pd.DataFrame: + df = read_csv(*args, **kwargs) + df.columns = pd.Index([*df.columns[:-2], "CD3.1", df.columns[-1]]) + return df + + monkeypatch.setattr(pd, "read_csv", mangling_read_csv) + path = tmp_path / "cells.csv" + path.write_text("cell_id,x,y,tile_num,CD3,CD8,CD3,label\n1,1.0,2.0,1,5.0,6.0,7.0,T\n") + + with pytest.warns(UserWarning, match=r"no 'cycCCC_chHHH' columns.*\['CD3', 'CD8', 'CD3 \(2\)'\]"): + table = read_cell_table(path, default_region="Kidney") + + assert list(table.var.index) == ["CD3", "CD8", "CD3 (2)"] + assert list(table.obs.columns) == ["x", "y", "tile_num", "label"] + assert list(table.region) == ["Kidney"] + + +@pytest.mark.parametrize( + ("content", "match"), + [ + ("cell_id,cyc001_ch001\n1,2.0\n", "'x' and 'y'"), + ("cell_id,x,y,cyc001_ch001\n1,1.0,2.0,2.0\n", "cannot be assigned to regions"), + ("cell_id,region,x,y,cyc001_ch001\n1,1,1,1,1\n1.5,1,2,2,1\n", r"1 cell\(s\) whose instance id in 'cell_id'"), + ( + "cell_id,region,x,y,cyc001_ch001\n1,1,1,1,1\n1,1,2,2,1\n", + r"1 duplicate instance id\(s\) \[1\] in region 'reg001'", + ), + ], + ids=["no_centroids", "no_region", "fractional_id", "duplicate_id"], +) +def test_bad_tables_raise(tmp_path: Path, content: str, match: str) -> None: + (tmp_path / "cells.csv").write_text(content) + + with pytest.raises(ValueError, match=match): + read_cell_table(tmp_path / "cells.csv", default_region=None) + + +@pytest.mark.parametrize("compression", [None, "zlib"]) +def test_read_labels_is_lazy(tmp_path: Path, compression: str | None) -> None: + mask = np.repeat(np.arange(64, dtype=np.uint16), 48).reshape(64, 48) + path = tmp_path / "mask.ome.tif" + tile = (16, 16) if compression else None + with tifffile.TiffWriter(path) as tw: + tw.write(np.stack([mask, mask + 1]), subifds=1, tile=tile, compression=compression, metadata={"axes": "CYX"}) + tw.write(np.stack([mask, mask + 1])[:, ::2, ::2], subfiletype=1, tile=tile, compression=compression) + translation = Translation([3, 4], axes=("x", "y")) + + with pytest.MonkeyPatch.context() as monkeypatch: + forbid_pixel_reads(monkeypatch) + labels = read_labels(path, {"scale_factors": None}, {"global": translation}) + + # of several planes, the first + np.testing.assert_array_equal(labels.values, mask) + assert get_transformation(labels, "global") == translation + tifffile.imwrite(path, np.zeros((8, 8), dtype=np.float32)) + with pytest.raises(ValueError, match="integer"): + read_labels(path, {}, {"global": Identity()}) + + +@pytest.mark.parametrize(("mask_shape", "factor"), [((32, 40), 2), ((22, 27), 3)]) +def test_downsampled_mask_is_scaled_onto_its_image(tmp_path: Path, mask_shape: tuple[int, int], factor: int) -> None: + # exactly, or with its shape rounded up; before the stage translation of its region + write_qptiff(tmp_path / "S_Scan1.qptiff", biomarkers=["DAPI"], shape=(64, 80), stage_position_cm=(1.0, 2.0)) + tifffile.imwrite(tmp_path / "mask.tif", np.zeros(mask_shape, dtype=np.uint32)) + + sdata = phenocycler(tmp_path, labels=tmp_path / "mask.tif", image_models_kwargs={"scale_factors": [2]}) + + image_transformation = get_transformation(sdata["S_image"], "global") + labels_transformation = get_transformation(sdata["S_labels"], "global") + assert isinstance(labels_transformation, Sequence) + assert labels_transformation.transformations == [Scale([factor, factor], axes=("y", "x")), image_transformation] + image_extent = get_extent(sdata["S_image"], coordinate_system="global") + labels_extent = get_extent(sdata["S_labels"], coordinate_system="global") + for axis, mask_size in zip(("y", "x"), mask_shape, strict=True): + assert labels_extent[axis][0] == pytest.approx(image_extent[axis][0]) + # a mask rounded up overhangs the image by less than one mask pixel + assert labels_extent[axis][1] == pytest.approx(image_extent[axis][0] + mask_size * factor) + + # writing may reorder the axes of the translation, so compare the composed affine maps + sdata.write(tmp_path / "x.zarr") + read = get_transformation(read_zarr(tmp_path / "x.zarr")["S_labels"], "global") + np.testing.assert_allclose( + read.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y")), + labels_transformation.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y")), + ) diff --git a/tests/readers/test_qptiff.py b/tests/readers/test_qptiff.py new file mode 100644 index 00000000..06ebf3c2 --- /dev/null +++ b/tests/readers/test_qptiff.py @@ -0,0 +1,133 @@ +import gc +import pickle +import warnings +from pathlib import Path +from typing import Any + +import dask.array as da +import numpy as np +import pytest +from tests._akoya_data import cim_experiment_v4_xml, fusion_scan_profile_json, write_qptiff + +from spatialdata_io.readers._utils import _qptiff +from spatialdata_io.readers._utils._qptiff import read_qptiff_image, read_qptiff_metadata + +_CIM = "Codex Instrument Manager 1.30.0.12" + + +@pytest.mark.parametrize( + ("kwargs", "names", "generation"), + [ + ({"biomarkers": ["DAPI", "CD3"], "names": ["DAPI", "ATTO 550"]}, ["DAPI", "CD3"], "fusion2"), + # PhenoImager (Opal) mode has no `` + ({"biomarkers": [None, None], "names": ["DAPI", "Opal 480"]}, ["DAPI", "Opal 480"], "fusion2"), + ({"biomarkers": ["DAPI"], "software": "Fusion 1.0.5"}, ["DAPI"], "fusion1"), + ({"biomarkers": [None], "names": ["FITC"], "software": "VectraPolaris 1.0"}, ["FITC"], "qpi_other"), + ], +) +def test_channel_names_and_generation(tmp_path: Path, kwargs: dict[str, Any], names: list[str], generation: str): + write_qptiff(tmp_path / "S_Scan1.qptiff", **kwargs) + + scan = read_qptiff_metadata(tmp_path / "S_Scan1.qptiff") + + assert [c.name for c in scan.channels] == names + assert scan.generation == generation + + +def test_pixel_size_and_stage_position(tmp_path: Path) -> None: + p = tmp_path / "x.qptiff" + write_qptiff(p, biomarkers=["DAPI"], stage_position_cm=(1.5, 0.25)) + scan = read_qptiff_metadata(p) + assert scan.pixel_size_um == pytest.approx(0.5069, rel=1e-6) + assert scan.stage_position_um == pytest.approx((15000.0, 2500.0)) + + # 96 dpi is a screen resolution, not a calibration + write_qptiff(p, biomarkers=["DAPI"], resolution_px_per_cm=96, resolution_unit="INCH") + with pytest.warns(UserWarning, match="Could not determine the pixel size"): + assert read_qptiff_metadata(p).pixel_size_um is None + + +def test_cim_without_label_page_counts_pages_by_image_type(tmp_path: Path) -> None: + p = tmp_path / "cim.qptiff" + biomarkers: list[str | None] = ["DAPI", "CD3", "CD8", "DAPI", "CD20"] + profile = cim_experiment_v4_xml([["DAPI", "CD3", "CD8"], ["DAPI", "CD20"]]) + written = write_qptiff(p, biomarkers=biomarkers, software=_CIM, scan_profile=profile, levels=3, label_page=False) + + scan = read_qptiff_metadata(p) + + assert scan.generation == "cim_qptiff" + assert [(c.name, c.cycle, c.channel) for c in scan.channels][:2] == [("DAPI", 1, 1), ("CD3", 1, 2)] + assert scan.nominal_pixel_size_um == pytest.approx(0.37739) + np.testing.assert_array_equal(read_qptiff_image(scan).compute(), written) + + +@pytest.mark.parametrize( + ("biomarkers", "profile", "expected"), + [ + # by name; a page without a record is left unset + (["DAPI", "CD3", "CD4"], [[("DAPI", "DAPI"), ("CD3", "ATTO 550")]], [(1, 1), (1, 2), (None, None)]), + # DAPI is in several cycles, so it is not paired + ( + ["DAPI", "CD3", "CD8"], + [[("DAPI", "D"), ("CD3", "A")], [("DAPI", "D"), ("CD8", "C")]], + [(None, None), (1, 2), (2, 2)], + ), + # names that differ only in non-ASCII letters, written as character references in the page XML + (["TCRγδ", "TCRαβ"], [[("TCRαβ", "A"), ("TCRγδ", "C")]], [(1, 2), (1, 1)]), + ], +) +def test_scan_profile_sets_cycle_and_channel( + tmp_path: Path, biomarkers: list[str | None], profile: list[list[tuple[str, str]]], expected: list[tuple] +) -> None: + p = tmp_path / "x.qptiff" + write_qptiff(p, biomarkers=biomarkers, scan_profile=fusion_scan_profile_json(profile)) + + assert [(c.cycle, c.channel) for c in read_qptiff_metadata(p).channels] == expected + + +@pytest.mark.parametrize( + "kwargs", + [ + {"biomarkers": ["DAPI", "CD3", "CD8"]}, # three channels are not RGB + {"biomarkers": ["DAPI"], "dtype": "uint16"}, + {"biomarkers": [None], "names": ["BF"], "rgb": True, "software": "VectraPolaris 1.0"}, + {"biomarkers": ["DAPI", "CD3"], "bigtiff": False, "tile": None}, + ], + ids=["multiplex", "uint16", "brightfield", "classic_stripped"], +) +def test_image_is_lazy_and_correct(tmp_path: Path, kwargs: dict[str, Any]) -> None: + p = tmp_path / "x.qptiff" + written = write_qptiff(p, **kwargs) + + scan = read_qptiff_metadata(p) + image = read_qptiff_image(scan) + + assert scan.rgb == kwargs.get("rgb", False) + assert isinstance(image, da.Array) + assert image.dtype == written.dtype + np.testing.assert_array_equal(image.compute(), written) + + +@pytest.mark.parametrize("series_match", [True, False]) +def test_image_is_picklable_and_closes_its_file( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch, series_match: bool +) -> None: + # dask.distributed pickles the graph, so it must not hold open file handles; also for the per-page fallback, + # used when tifffile's series do not match the full-resolution pages of the XML + p = tmp_path / "x.qptiff" + written = write_qptiff(p, biomarkers=["DAPI", "CD3"]) + if not series_match: + real = _qptiff.lazy_tiff_array + + def wrong_series(path: Path, *, series: int = 0, page: int | None = None) -> da.Array: + return real(path, series=series, page=page) if page is not None else da.zeros((1, 1, 1)) + + monkeypatch.setattr(_qptiff, "lazy_tiff_array", wrong_series) + with warnings.catch_warnings(record=True) as record: + warnings.simplefilter("always", ResourceWarning) + image = read_qptiff_image(read_qptiff_metadata(p)) + restored = pickle.loads(pickle.dumps(image)) + np.testing.assert_array_equal(restored.compute(), written) + del image, restored + gc.collect() + assert not [warning for warning in record if issubclass(warning.category, ResourceWarning)] diff --git a/tests/test_akoya_real_data.py b/tests/test_akoya_real_data.py new file mode 100644 index 00000000..31516e2c --- /dev/null +++ b/tests/test_akoya_real_data.py @@ -0,0 +1,180 @@ +"""Tests of the Akoya readers on real public datasets. + +They run only when the environment variable ``AKOYA_TEST_DATA`` points to a directory holding the files below (paths +relative to it). The files are public (Zenodo, Figshare, BioImage Archive, GitHub) but too large to +download in CI. +""" + +from __future__ import annotations + +import os +from dataclasses import dataclass +from pathlib import Path + +import dask.array as da +import numpy as np +import pytest +import tifffile +from spatialdata.models import get_channel_names +from spatialdata.transformations import get_transformation + +from spatialdata_io import codex, phenocycler + +_ROOT = os.environ.get("AKOYA_TEST_DATA") + +pytestmark = [ + pytest.mark.skipif(_ROOT is None, reason="set AKOYA_TEST_DATA to the directory of the real Akoya datasets"), + # upstream noise unrelated to the readers: ome-zarr's deprecated `Scaler`; unclosed files must not be left behind + pytest.mark.filterwarnings( + "ignore:Call to deprecated class Scaler", + "error::ResourceWarning", + "error::pytest.PytestUnraisableExceptionWarning", + ), +] + + +@dataclass(frozen=True) +class Scan: + """What a real dataset is expected to read as: ``shape`` starts ``(c, y, x)``, ``names`` the first channels.""" + + path: str + region: str + generation: str + names: tuple[str, ...] + pixel_size_um: float | None + shape: tuple[int, ...] + dtype: str | None = None + + +# one per generation and file format +_SCANS = [ + Scan( + "figshare_25123538/AML2_183_Scan1.qptiff", "AML2_183", "fusion1", ("DAPI", "NPM1C"), 0.5069, (54, 10800, 7680) + ), + # Fusion 2.x in PhenoImager (Opal) mode: the names come from ``, there is no `` + Scan("bia_S-BIAD3467/Tonsil-ctrl_Scan1.qptiff", "Tonsil-ctrl", "fusion2", ("DAPI", "Opal 690"), 0.4994, (6,)), + # tifffile's QPI series builder finds only 3 of the 5 levels; the full-resolution shape must still be right + Scan( + "zenodo_22107401/Arlene_Tonsil_9x9.qptiff", + "Arlene_Tonsil_9x9", + "cim_qptiff", + ("DAPI-01",), + 0.3774, + (43, 9508, 12672), + ), + # CODEX Uploader montage, TCYX (21, 4, y, x); the names come from `channelNames.txt` + Scan("bia_S-BIAD859/stitched/6371v1", "reg001", "legacy_processor", ("HOECHST", "CollIV"), None, (84, 3528, 3360)), + # CODEX Processor per-marker stitched images `reg001_cycCCC_chHHH_.tiff` + Scan( + "zenodo_7742474/human_pancreas_codex_zip_subset", + "reg001", + "legacy_processor", + ("Hoechst", "CGC"), + None, + (2, 18440, 25816), + ), + # Vectra / PhenoImager (generation `qpi_other`): unmixed uint16, brightfield RGB and inForm float32 component data + Scan( + "zenodo_22107401/Panel33C_Scan2.unmixed.qptiff", + "Panel33C_Scan2.unmixed", + "qpi_other", + ("CD21",), + 0.4965, + (8, 36000, 24960), + "uint16", + ), + Scan( + "ome_vectra_qptiff/HandEcompressed_Scan1.qptiff", + "HandEcompressed", + "qpi_other", + ("r", "g", "b"), + 0.4989, + (3, 26640, 30720), + "uint8", + ), + Scan( + "ome_vectra_qptiff/LuCa-7color_[13860,52919]_1x1component_data.tif", + "LuCa-7color__13860_52919__1x1component_data", + "qpi_other", + ("PDL1 (Opal 520)",), + 0.498, + (8, 1400, 1868), + "float32", + ), +] + + +@pytest.mark.parametrize("scan", _SCANS, ids=[scan.path.split("/")[-1] for scan in _SCANS]) +def test_real_scan(scan: Scan) -> None: + path = Path(_ROOT) / scan.path # type: ignore[arg-type] + if not path.exists(): + pytest.skip(f"{path} is missing") + reader = codex if scan.generation in ("cim_qptiff", "legacy_processor") else phenocycler + + sdata = reader(path) + + image = sdata.images[f"{scan.region}_image"] + data = image["scale0"]["image"].data + assert isinstance(data, da.Array) + assert data.shape[: len(scan.shape)] == scan.shape + assert scan.dtype is None or data.dtype == scan.dtype + names = get_channel_names(image) + assert tuple(names[: len(scan.names)]) == scan.names + assert len(set(names)) == len(names) + metadata = sdata.attrs["akoya"]["regions"][scan.region] + assert metadata["generation"] == scan.generation + # a brightfield scan has one RGB page, which is read as three channels + assert len(metadata["channels"]) == (1 if names == ["r", "g", "b"] else scan.shape[0]) + if scan.pixel_size_um is None: + assert metadata["pixel_size_um"] is None + else: + assert metadata["pixel_size_um"] == pytest.approx(scan.pixel_size_um, abs=1e-3) + assert sdata.attrs["spatialdata_io_reader"] == reader.__name__ + if path.is_file(): + _assert_tile_matches_tifffile(path, data, channel=scan.shape[0] - 1) + + +def _assert_tile_matches_tifffile(path: Path, data: da.Array, channel: int) -> None: + """A full-resolution tile of ``channel`` equals tifffile's decoding of the whole page (not through zarr).""" + with tifffile.TiffFile(path) as tf: + full = [page for page in tf.pages if page.shape[:2] == data.shape[1:]] + # brightfield: one RGB page holds the three channels + page = full[0] if full[0].ndim == 3 else full[channel] + y, x = data.shape[1] // 2, data.shape[2] // 2 + expected = page.asarray()[y : y + 256, x : x + 256] + if expected.ndim == 3: + expected = expected[..., channel] + np.testing.assert_array_equal(data[channel, y : y + 256, x : x + 256].compute(), expected) + + +def test_real_same_qptiff_through_both_readers() -> None: + path = Path(_ROOT) / "figshare_25123538/AML2_183_Scan1.qptiff" # type: ignore[arg-type] + if not path.exists(): + pytest.skip(f"{path} is missing") + + a, b = phenocycler(path), codex(path) + + image_a, image_b = a.images["AML2_183_image"], b.images["AML2_183_image"] + assert get_channel_names(image_a) == get_channel_names(image_b) + assert list(image_a) == list(image_b) + coarsest = list(image_a)[-1] + np.testing.assert_array_equal(image_a[coarsest]["image"].values, image_b[coarsest]["image"].values) + assert get_transformation(image_a, get_all=True) == get_transformation(image_b, get_all=True) + assert a.attrs["akoya"]["regions"] == b.attrs["akoya"]["regions"] + + +def test_real_codex_processor_csv() -> None: + # a 500-cell subset of a CODEX Processor CSV with the `cyc001_ch001:DAPI1` header (github.com/igordot/phenomenalist) + path = Path(_ROOT) / "github_codex_metadata/phenomenalist/tonsil-akoya-2018-500.csv" # type: ignore[arg-type] + if not path.exists(): + pytest.skip(f"{path} is missing") + + with pytest.warns(UserWarning, match="cells are represented as circles"): + sdata = codex(path) + + table = sdata.tables["table"] + assert table.shape == (500, 32) + assert list(table.var_names[:5]) == ["DAPI1", "Blank", "Blank (2)", "Blank (3)", "DAPI2"] + assert list(table.obs["region"].cat.categories) == ["reg001_cells"] + assert table.obsm["spatial"].shape == (500, 2) + assert len(sdata.shapes["reg001_cells"]) == 500 diff --git a/tests/test_cli_alignment.py b/tests/test_cli_alignment.py index f2d6e6b1..5631943a 100644 --- a/tests/test_cli_alignment.py +++ b/tests/test_cli_alignment.py @@ -43,6 +43,7 @@ ("spatialdata_io.readers.macsima", "macsima", "macsima_wrapper"), ("spatialdata_io.readers.mcmicro", "mcmicro", "mcmicro_wrapper"), ("spatialdata_io.readers.merscope", "merscope", "merscope_wrapper"), + ("spatialdata_io.readers.phenocycler", "phenocycler", "phenocycler_wrapper"), ("spatialdata_io.readers.seqfish", "seqfish", "seqfish_wrapper"), ("spatialdata_io.readers.steinbock", "steinbock", "steinbock_wrapper"), ("spatialdata_io.readers.stereoseq", "stereoseq", "stereoseq_wrapper"), diff --git a/tests/test_codex.py b/tests/test_codex.py new file mode 100644 index 00000000..b9ed8c07 --- /dev/null +++ b/tests/test_codex.py @@ -0,0 +1,400 @@ +"""Tests for the CODEX reader on synthetic CODEX Processor/Uploader trees and CIM qptiff files.""" + +from __future__ import annotations + +import json +import shutil +import warnings +from typing import TYPE_CHECKING, Literal + +import numpy as np +import pandas as pd +import pytest +import tifffile +from spatialdata import read_zarr +from spatialdata.models import get_channel_names +from spatialdata.transformations import Identity, Scale, Sequence, get_transformation + +from spatialdata_io import codex, phenocycler +from spatialdata_io.readers._utils._akoya import dedupe_channel_names +from spatialdata_io.readers.codex import _discover_legacy +from tests._akoya_data import ( + forbid_pixel_reads, + write_fcs, + write_legacy_processor_tree, + write_qptiff, +) + +if TYPE_CHECKING: + from pathlib import Path + +_FAST = {"scale_factors": [2]} + + +def _level0(sdata, name): # type: ignore[no-untyped-def] + return sdata.images[name]["scale0"]["image"] + + +def _codex(path: Path | str, **kwargs): # type: ignore[no-untyped-def] + """`codex()` with fast pyramids, ignoring the warning that cells without a mask are circles.""" + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", "No segmentation mask found", UserWarning) + return codex(path, image_models_kwargs=_FAST, labels_models_kwargs=_FAST, **kwargs) + + +@pytest.mark.parametrize("layout", ["per_marker", "hyperstack"]) +def test_codex_images_are_lazy_and_named(tmp_path: Path, layout: Literal["per_marker", "hyperstack"]) -> None: + data = write_legacy_processor_tree(tmp_path, regions=2, layout=layout) + + with pytest.MonkeyPatch.context() as monkeypatch: + forbid_pixel_reads(monkeypatch) + sdata = _codex(tmp_path) + + assert sorted(sdata.images) == ["reg001_image", "reg002_image"] + for region in ("reg001", "reg002"): + assert get_channel_names(sdata.images[f"{region}_image"]) == dedupe_channel_names(data["markers"]) + expected = data["images"][region] + np.testing.assert_array_equal(_level0(sdata, f"{region}_image").values, expected.reshape(-1, 128, 160)) + attrs = sdata.attrs["akoya"]["regions"][region] + assert attrs["generation"] == "legacy_processor" + assert [(c["cycle"], c["channel"]) for c in attrs["channels"][3:5]] == [(1, 4), (2, 1)] + assert sdata.attrs["spatialdata_io_reader"] == "codex" + + +@pytest.mark.parametrize( + ("labels", "names_txt", "expected"), + [ + # ImageJ labels win over channelNames.txt + (["L0", "L1", "L2", "L3"], "o0\no1\no2\no3\n", ["L0", "L1", "L2", "L3"]), + # ImageJ writes an empty label for every unlabelled slice; such labels do not replace channelNames.txt + (["cyc1.tif", "", "", ""], "DAPI\nCD3\nDAPI2\nCD8\n", ["DAPI", "CD3", "DAPI2", "CD8"]), + (None, "DAPI\nIFNµ\nCD3\nCD8\n".encode("utf-16"), ["DAPI", "IFNµ", "CD3", "CD8"]), + (None, "DAPI\nIFNµ\nCD3\nCD8\n".encode("cp1252"), ["DAPI", "IFNµ", "CD3", "CD8"]), + # a channelNames.txt of another length is ignored with a warning + (None, "DAPI\nCD3\n", ["cyc001_ch001", "cyc001_ch002", "cyc002_ch001", "cyc002_ch002"]), + ], + ids=["imagej_labels", "blank_imagej_labels", "utf16", "cp1252", "wrong_length"], +) +def test_codex_hyperstack_channel_names( + tmp_path: Path, labels: list[str] | None, names_txt: str | bytes, expected: list[str] +) -> None: + metadata = {"axes": "TCYX"} if labels is None else {"axes": "TCYX", "Labels": labels} + tifffile.imwrite( + tmp_path / "reg001_montage.tif", np.zeros((2, 2, 16, 16), np.uint16), imagej=True, metadata=metadata + ) + names = tmp_path / "channelNames.txt" + names.write_bytes(names_txt) if isinstance(names_txt, bytes) else names.write_text(names_txt) + + with warnings.catch_warnings(record=True) as record: + warnings.simplefilter("always") + sdata = _codex(tmp_path) + + assert get_channel_names(sdata.images["reg001_image"]) == expected + assert any("channelNames.txt" in str(w.message) for w in record) == (expected[0] == "cyc001_ch001") + + +@pytest.mark.parametrize("from_subdirectory", [False, True]) +def test_codex_sidecars_case_insensitive_and_above_path(tmp_path: Path, from_subdirectory: bool) -> None: + # SPACEc and the HuBMAP codex-pipeline write `channelnames.txt`; older exports `Experiment.json` + data = write_legacy_processor_tree(tmp_path, cycles=2, channels=4, layout="hyperstack", fcs=False, csv=False) + (tmp_path / "channelNames.txt").rename(tmp_path / "channelnames.txt") + (tmp_path / "experiment.json").unlink() + (tmp_path / "Experiment.json").write_text(json.dumps({"xyResolution": 377.5})) + + sdata = _codex(tmp_path / "processed_2020-01-01" / "stitched" if from_subdirectory else tmp_path) + + assert get_channel_names(sdata.images["reg001_image"]) == dedupe_channel_names(data["markers"]) + assert sdata.attrs["akoya"]["regions"]["reg001"]["pixel_size_um"] == pytest.approx(0.3775) + + +def test_codex_regions_channel_names_and_positional_call(tmp_path: Path) -> None: + write_legacy_processor_tree(tmp_path, regions=2, cycles=1, channels=2, fcs=False, csv=False) + + sdata = _codex(tmp_path, regions=["reg002"], channel_names=["a", "b"]) + + assert list(sdata.images) == ["reg002_image"] + assert get_channel_names(sdata.images["reg002_image"]) == ["a", "b"] + assert sdata.coordinate_systems == ["global"] + # `codex(path, fcs, imread_kwargs)` is the signature of earlier releases; `imread_kwargs` is ignored + with pytest.warns(UserWarning, match="imread_kwargs"): + assert sorted(codex(tmp_path, True, {"nframes": 1}).images) == ["reg001_image", "reg002_image"] + for kwargs, error in (({"regions": ["reg003"]}, "reg003"), ({"channel_names": ["a"]}, "channel_names")): + with pytest.raises(ValueError, match=error): + codex(tmp_path, **kwargs) + with pytest.raises(FileNotFoundError, match="reg"): + codex(tmp_path / "processed_2020-01-01" / "segm") + + +def test_codex_and_phenocycler_name_qptiffs_alike(tmp_path: Path) -> None: + write_qptiff(tmp_path / "exp_reg001.qptiff", biomarkers=["DAPI"], dtype="uint8") + write_qptiff(tmp_path / "exp_reg001.er.qptiff", biomarkers=["DAPI"], dtype="uint16") + write_qptiff(tmp_path / "Tonsil_Scan1.er.qptiff", biomarkers=["DAPI"], dtype="uint16") + write_qptiff(tmp_path / "Preg12_Scan1.qptiff", biomarkers=["DAPI"]) + + for reader in (codex, phenocycler): + sdata = reader(tmp_path, image_models_kwargs=_FAST) + assert sorted(sdata.images) == ["Preg12_image", "Tonsil_image", "reg001_image"] + assert sdata.attrs["akoya"]["regions"]["reg001"]["path"].endswith("exp_reg001.er.qptiff") + assert list(reader(tmp_path / "exp_reg001.qptiff").images) == ["reg001_image"] + + write_qptiff(tmp_path / "other_reg1.qptiff", biomarkers=["DAPI"]) + for reader in (codex, phenocycler): + with pytest.raises(ValueError, match=r"Several qptiff files map to the same region\(s\) \['reg001'\]"): + reader(tmp_path) + + +def test_codex_region_ids_and_ignored_files(tmp_path: Path) -> None: + for name in ("reg001_montage.tif", "exp_reg002.tif", "1_reg003.tif", "Treg2.tif", "coreg1.tif"): + tifffile.imwrite(tmp_path / name, np.zeros((1, 4, 4), dtype=np.uint8)) + # raw acquisition and intermediate folders, macOS AppleDouble files and a link back up the tree + for raw in ("cyc001_reg004/1_00001_Z001_CH1.tif", "tiles/reg005_X01_Y01.tif", "__MACOSX/reg006.tif"): + (tmp_path / raw).parent.mkdir() + tifffile.imwrite(tmp_path / raw, np.zeros((4, 4), dtype=np.uint8)) + (tmp_path / "._reg007.tif").write_bytes(b"AppleDouble") + (tmp_path / "loop").symlink_to(tmp_path, target_is_directory=True) + + assert sorted(_discover_legacy(tmp_path, fcs=True).images) == ["reg001", "reg002", "reg003"] + + +def test_codex_discovers_tables(tmp_path: Path) -> None: + cells = pd.DataFrame({"cell_id": [1.0, 2.0], "x": [1.0, 2.0], "y": [3.0, 4.0], "cyc001_ch001": [5.0, 6.0]}) + tifffile.imwrite(tmp_path / "Kidney.tif", np.zeros((2, 16, 16), dtype=np.uint16)) + cells.to_csv(tmp_path / "Kidney.csv", index=False) + # not cell tables: sidecar CSVs (also named after a region), flow cytometry tubes and broken FCS files + (tmp_path / "flow").mkdir() + write_fcs(tmp_path / "flow" / "tube1.fcs", pd.DataFrame({"FSC-A": [1.0], "SSC-A": [2.0]})) + (tmp_path / "flow" / "broken.fcs").write_bytes(b"not an fcs file") + # flowio raises an EOFError and a KeyError for these + write_fcs(tmp_path / "flow" / "reg001.fcs", cells) + content = (tmp_path / "flow" / "reg001.fcs").read_bytes() + (tmp_path / "flow" / "reg001.fcs").write_bytes(content[:100]) + (tmp_path / "flow" / "reg002.fcs").write_bytes(content.replace(b"$DATATYPE", b"$DATATYPX")) + (tmp_path / "reg001_antibodies.csv").write_text("Antibody,Clone\nCD3,UCHT1\n") + (tmp_path / "xy.csv").write_text("x,y\n1.0,2.0\n") + assert _discover_legacy(tmp_path, fcs=True).tables == {"*": tmp_path / "Kidney.csv"} + + # a CODEX gate export next to the table: the reader cannot tell which one holds the cells + write_fcs(tmp_path / "flow" / "CD3pos.fcs", cells.iloc[:1]) + with pytest.raises(ValueError, match="cannot choose"): + _discover_legacy(tmp_path, fcs=True) + + # compensated beats uncompensated before the file format; `fcs` picks the format of one table + fcs_dir = tmp_path / "exp" / "segm" / "segm-1" / "fcs" + for kind in ("compensated", "uncompensated"): + (fcs_dir / kind).mkdir(parents=True) + cells.to_csv(fcs_dir / "compensated" / "reg001_compensated.csv", index=False) + write_fcs(fcs_dir / "compensated" / "reg001_compensated.fcs", cells) + write_fcs(fcs_dir / "uncompensated" / "reg001_uncompensated.fcs", cells) + for fcs, suffix in ((True, ".fcs"), (False, ".csv")): + assert _discover_legacy(tmp_path / "exp", fcs=fcs).tables == { + "reg001": fcs_dir / f"compensated/reg001_compensated{suffix}" + } + + # two segmentation runs of one region are not resolved by their path + segm_10 = tmp_path / "exp" / "segm" / "segm-10" / "fcs" / "compensated" + segm_10.mkdir(parents=True) + cells.to_csv(segm_10 / "reg001_compensated.csv", index=False) + with pytest.raises(ValueError, match=r"segm-1/.*segm-10/|segm-10/.*segm-1/"): + _discover_legacy(tmp_path / "exp", fcs=True) + + +def test_codex_tables_of_several_regions(tmp_path: Path) -> None: + data = write_legacy_processor_tree(tmp_path, regions=2, n_cells=50) + shutil.rmtree(tmp_path / "processed_2020-01-01" / "stitched" / "reg002") + + with pytest.warns(UserWarning, match=r"No segmentation mask found for region\(s\) \['reg001', 'reg002'\]"): + sdata = codex(tmp_path, image_models_kwargs=_FAST) + + table = sdata.tables["table"] + assert table.n_obs == 100 + assert isinstance(table.X, np.ndarray) + assert table.X.dtype == np.float32 + assert list(table.obs["region"].cat.categories) == ["reg001_cells", "reg002_cells"] + assert list(table.var_names) == dedupe_channel_names(data["markers"]) + reg002 = table[table.obs["region"] == "reg002_cells"] + np.testing.assert_allclose(reg002.obsm["spatial"], data["cells"][["x", "y"]], rtol=1e-6) + radius = sdata.shapes["reg002_cells"]["radius"].to_numpy() + np.testing.assert_allclose(radius, np.sqrt(data["cells"]["size"] / np.pi), rtol=1e-6) + assert sdata.attrs["akoya"]["cell_tables"]["reg002"].endswith("reg002_compensated.fcs") + # regions without stage positions, and a region without an image, have their own coordinate systems + assert sorted(sdata.coordinate_systems) == ["reg001", "reg002"] + for name in ("reg001_image", "reg001_cells"): + assert get_transformation(sdata[name], get_all=True) == {"reg001": Identity()} + assert get_transformation(sdata["reg002_cells"], get_all=True) == {"reg002": Identity()} + assert list(_codex(tmp_path, regions=["reg002"]).shapes) == ["reg002_cells"] + + sdata.write(tmp_path / "x.zarr") + read = read_zarr(tmp_path / "x.zarr") + np.testing.assert_array_equal(_level0(read, "reg001_image").values, _level0(sdata, "reg001_image").values) + assert get_transformation(read["reg002_cells"], get_all=True) == {"reg002": Identity()} + assert read.tables["table"].obs.index.name is None + assert read.tables["table"].obs["cell_id"].dtype == np.int64 + assert read.attrs["akoya"] == sdata.attrs["akoya"] + assert read.attrs["spatialdata_io_reader"] == "codex" + + +def test_codex_tables_with_different_panels_or_formats(tmp_path: Path) -> None: + columns = {"cell_id": [1.0], "x": [1.0], "y": [2.0]} + write_fcs(tmp_path / "reg001.fcs", pd.DataFrame(columns | {"cyc001_ch001": [3.0]}), {"cyc001_ch001": "DAPI"}) + pd.DataFrame(columns | {"cyc001_ch001:DAPI": [3.0], "cyc001_ch002:CD3": [4.0]}).to_csv( + tmp_path / "reg002.csv", index=False + ) + + sdata = _codex(tmp_path) + + assert sorted(sdata.tables) == ["reg001_table", "reg002_table"] + assert list(sdata.tables["reg002_table"].var_names) == ["DAPI", "CD3"] + assert sdata.tables["reg001_table"].uns["spatialdata_attrs"]["region"] == ["reg001_cells"] + + # the CSV and FCS tables of one panel share a table; their source columns differ. With `fcs=False`, a region with + # only an FCS table is read from it + pd.DataFrame(columns | {"cyc001_ch001:DAPI": [3.0]}).to_csv(tmp_path / "reg002.csv", index=False) + with pytest.warns(UserWarning, match="no CSV cell table was found"): + sdata = _codex(tmp_path, fcs=False) + assert list(sdata.tables["table"].var_names) == ["DAPI"] + assert "source_column" not in sdata.tables["table"].var + + +@pytest.mark.parametrize( + ("images", "rows", "regions", "expected"), + [ + # the region is the file name without its image extension, `.ome` included + (["Kidney.ome.tif"], ["1,1,1,2"], None, ["Kidney_cells"]), + # region values are matched after `sanitize_name`, like the region ids + ( + ["Tonsil 1.tif", "Tonsil 2.tif"], + ["1,Tonsil 1,1,2", "1,Tonsil 2,3,4"], + None, + ["Tonsil_1_cells", "Tonsil_2_cells"], + ), + # a region without an image is still a region + (["reg001.tif"], ["1,1,1,2", "1,2,3,4"], ["reg002"], ["reg002_cells"]), + ], +) +def test_codex_regionless_table_regions( + tmp_path: Path, images: list[str], rows: list[str], regions: list[str] | None, expected: list[str] +) -> None: + for name in images: + tifffile.imwrite(tmp_path / name, np.zeros((2, 16, 16), dtype=np.uint16)) + (tmp_path / "cells.csv").write_text("cell_id,region,x,y,cyc001_ch001\n" + "".join(f"{row},5\n" for row in rows)) + + sdata = _codex(tmp_path, regions=regions) + + assert list(sdata.tables["table"].obs["region"].cat.categories) == expected + assert sorted(sdata.shapes) == expected + + +@pytest.mark.parametrize( + ("setup", "match"), + [ + ("two_images", "cannot be assigned to regions"), + ("duplicate_ids", r"reg001_compensated.csv has 1 duplicate instance id"), + ("overlapping", r"Several cell tables contain the cells of region\(s\) \['reg001'\]"), + ("other_region", "has no cells of region 'reg002'"), + ], +) +def test_codex_table_errors(tmp_path: Path, setup: str, match: str) -> None: + data = write_legacy_processor_tree(tmp_path, regions=2, cycles=1, channels=1, fcs=False, n_cells=3) + table = next(tmp_path.rglob("reg002_compensated.csv" if setup == "other_region" else "reg001_compensated.csv")) + if setup == "two_images": + table.rename(tmp_path / "cells.csv") + data["cells"].drop(columns="region").to_csv(tmp_path / "cells.csv", index=False) + elif setup == "duplicate_ids": + data["cells"].assign(cell_id=[1.0, 1.0, 2.0], region=1.0).to_csv(table, index=False) + elif setup == "overlapping": + (tmp_path / "all_cells.csv").write_text("cell_id,region,x,y,cyc001_ch001\n1,1,1.0,2.0,3.0\n") + else: + # the table of reg002 whose region column names reg001 + data["cells"].assign(region=1.0).to_csv(table, index=False) + + with pytest.raises(ValueError, match=match): + codex(tmp_path, fcs=False) + + +@pytest.mark.parametrize( + ("image", "mask", "region"), + [ + ("expr/reg001_expr.ome.tiff", "mask/reg001_mask.ome.tiff", "reg001"), + ("Kidney sample.tif", "Kidney sample_labels.tif", "Kidney_sample"), + ("Tonsil_Scan1.er.qptiff", "Tonsil_Scan1_mask.tif", "Tonsil"), + ("Tonsil_Scan1.er.qptiff", "Tonsil_mask.tif", "Tonsil"), + ], +) +def test_codex_discovers_the_mask_of_a_region(tmp_path: Path, image: str, mask: str, region: str) -> None: + for name in (image, mask): + (tmp_path / name).parent.mkdir(exist_ok=True) + if image.endswith("qptiff"): + write_qptiff(tmp_path / image, biomarkers=["DAPI"], shape=(16, 16)) + else: + tifffile.imwrite(tmp_path / image, np.zeros((2, 16, 16), dtype=np.uint16)) + tifffile.imwrite(tmp_path / mask, np.arange(256, dtype=np.uint16).reshape(16, 16)) + (tmp_path / "cells.csv").write_text("cell_id,x,y,cyc001_ch001\n1,0.0,0.0,1.0\n2,1.0,0.0,2.0\n") + + sdata = _codex(tmp_path) + + assert list(sdata.images) == [f"{region}_image"] + assert list(sdata.labels) == [f"{region}_labels"] + assert sdata.tables["table"].uns["spatialdata_attrs"]["region"] == [f"{region}_labels"] + assert sdata.attrs["akoya"]["labels"] == {region: str(tmp_path / mask)} + + +def test_codex_mask_discovery_errors_and_explicit_labels(tmp_path: Path) -> None: + for region in ("reg001", "reg002"): + tifffile.imwrite(tmp_path / f"{region}_montage.tif", np.zeros((2, 16, 16), dtype=np.uint16)) + tifffile.imwrite(tmp_path / "Other_mask.tif", np.zeros((16, 16), dtype=np.uint16)) + with pytest.warns(UserWarning, match="Ignoring the mask .*Other_mask.tif"): + assert not _codex(tmp_path).labels + + for name in ("reg002_mask.tif", "reg002_labels.tif"): + tifffile.imwrite(tmp_path / name, np.ones((16, 16), dtype=np.uint16)) + with pytest.raises(ValueError, match="several masks for region 'reg002'"): + _codex(tmp_path) + + # given labels disable mask discovery and are not read as images + tifffile.imwrite(tmp_path / "reg001_x_mask.tif", np.ones((16, 16), dtype=np.uint16)) + with warnings.catch_warnings(): + warnings.simplefilter("error") + sdata = _codex(tmp_path, labels={"reg001": tmp_path / "reg001_x_mask.tif"}) + assert sorted(sdata.images) == ["reg001_image", "reg002_image"] + assert list(sdata.labels) == ["reg001_labels"] + + +def test_codex_explicit_table_and_downsampled_mask(tmp_path: Path) -> None: + write_legacy_processor_tree(tmp_path, regions=2, n_cells=50) + (tmp_path / "mine.csv").write_text("cell_id,x,y,cyc001_ch001:CD3\n1,1.0,4.0,1.0\n2,3.0,6.0,2.0\n") + mask = np.zeros((64, 80), dtype=np.uint16) + mask[2, 0], mask[3, 1] = 1, 2 + tifffile.imwrite(tmp_path / "mask.tif", mask) + + with warnings.catch_warnings(): + warnings.simplefilter("error") + sdata = _codex( + tmp_path, cells_table={"reg002": tmp_path / "mine.csv"}, labels={"reg002": tmp_path / "mask.tif"} + ) + + table = sdata.tables["table"] + assert list(table.var_names) == ["CD3"] + assert table.uns["spatialdata_attrs"]["region"] == ["reg002_labels"] + assert not sdata.shapes + assert sdata.attrs["akoya"]["cell_tables"] == {"reg002": str(tmp_path / "mine.csv")} + # the mask is scaled onto the image; the centroids stay in image pixels + labels_transformation = get_transformation(sdata["reg002_labels"], "reg002") + assert isinstance(labels_transformation, Sequence) + assert labels_transformation.transformations[0] == Scale([2, 2], axes=("y", "x")) + np.testing.assert_allclose(table.obsm["spatial"], [[1, 4], [3, 6]]) + # an explicit table of an unselected region does not bring the discovered tables back + assert not _codex(tmp_path, regions=["reg001"], cells_table={"reg002": tmp_path / "mine.csv"}).tables + + +@pytest.mark.parametrize("fcs", [True, False]) +def test_codex_old_flat_layout_regression(tmp_path: Path, monkeypatch: pytest.MonkeyPatch, fcs: bool) -> None: + # the original reader read the image (and the CSV) relative to the working directory, and ignored `fcs` + data = write_legacy_processor_tree(tmp_path / "data", layout="old_flat", fcs=fcs) + monkeypatch.chdir(tmp_path) + + sdata = _codex(tmp_path / "data", fcs=fcs) + + np.testing.assert_array_equal(_level0(sdata, "reg001_image").values, data["images"]["reg001"][:1]) + table = sdata.tables["table"] + assert list(table.obs["region"].cat.categories) == ["reg001_cells"] + np.testing.assert_allclose(table.X, data["cells"].filter(regex="^cyc").to_numpy(), rtol=1e-6) diff --git a/tests/test_phenocycler.py b/tests/test_phenocycler.py new file mode 100644 index 00000000..de3ccad1 --- /dev/null +++ b/tests/test_phenocycler.py @@ -0,0 +1,287 @@ +"""Tests for the PhenoCycler-Fusion reader on synthetic qptiff files.""" + +from __future__ import annotations + +import json +import warnings +from pathlib import Path +from typing import Any + +import dask.array as da +import numpy as np +import pytest +import tifffile +from click.testing import CliRunner +from spatialdata import read_zarr +from spatialdata.models import get_channel_names +from spatialdata.transformations import Identity, Scale, Sequence, Translation, get_transformation + +from spatialdata_io import phenocycler +from spatialdata_io.__main__ import phenocycler_wrapper +from tests._akoya_data import forbid_pixel_reads, write_qptiff + +_FAST = {"scale_factors": [2]} + + +def _level0(sdata, name): # type: ignore[no-untyped-def] + return sdata.images[name]["scale0"]["image"] + + +def _write_mask(path: Path, shape: tuple[int, int] = (256, 320)) -> np.ndarray: + mask = np.zeros(shape, dtype=np.uint32) + mask[10:20, 10:20], mask[30:40, 10:20] = 1, 2 + tifffile.imwrite(path, mask) + return mask + + +def test_read_is_lazy_and_round_trips(tmp_path: Path) -> None: + path = tmp_path / "Tonsil_Scan1.qptiff" + data = write_qptiff(path, biomarkers=["DAPI", "CD3", "CD3"]) + + with pytest.MonkeyPatch.context() as monkeypatch: + forbid_pixel_reads(monkeypatch) + sdata = phenocycler(path, image_models_kwargs=_FAST) + + assert isinstance(_level0(sdata, "Tonsil_image").data, da.Array) + assert get_channel_names(sdata.images["Tonsil_image"]) == ["DAPI", "CD3", "CD3 (2)"] + np.testing.assert_array_equal(_level0(sdata, "Tonsil_image").values, data) + assert sdata.attrs["akoya"]["regions"]["Tonsil"]["generation"] == "fusion2" + sdata.write(tmp_path / "x.zarr") + read = read_zarr(tmp_path / "x.zarr") + assert get_channel_names(read.images["Tonsil_image"]) == ["DAPI", "CD3", "CD3 (2)"] + assert read.attrs["akoya"] == sdata.attrs["akoya"] + assert read.attrs["spatialdata_io_reader"] == "phenocycler" + + +@pytest.mark.parametrize( + ("name", "kwargs", "region", "channels"), + [ + ( + "Slide_Scan1.unmixed.qptiff", + {"biomarkers": [None, None], "names": ["PDL1 (Opal 520)", "DAPI"], "dtype": "float32"}, + "Slide_Scan1.unmixed", + ["PDL1 (Opal 520)", "DAPI"], + ), + ("HandE_Scan1.qptiff", {"biomarkers": [None], "names": ["BF"], "rgb": True}, "HandE", ["r", "g", "b"]), + # inForm component data; the brackets and comma of its name are not valid in an element name + ( + "LuCa_[13860,52919]_1x1component_data.tif", + {"biomarkers": [None], "names": ["A"], "dtype": "float32", "levels": 1}, + "LuCa__13860_52919__1x1component_data", + ["A"], + ), + ], + ids=["unmixed", "brightfield", "inform_component"], +) +def test_vectra_and_phenoimager_files( + tmp_path: Path, name: str, kwargs: dict[str, Any], region: str, channels: list[str] +) -> None: + data = write_qptiff(tmp_path / name, software="VectraPolaris 1.0", **kwargs) + + sdata = phenocycler(tmp_path / name, image_models_kwargs=_FAST) + + image = _level0(sdata, f"{region}_image") + assert get_channel_names(sdata.images[f"{region}_image"]) == channels + np.testing.assert_array_equal(image.values, np.moveaxis(data, -1, 0) if kwargs.get("rgb") else data) + assert sdata.attrs["akoya"]["regions"][region]["generation"] == "qpi_other" + + +@pytest.mark.parametrize( + ("files", "match"), + [ + ([], r"No \.qptiff files"), + (["A [1].qptiff", "A_[1].qptiff"], r"same region\(s\) \['A__1_'\]"), + (["plain.tif"], "not a PerkinElmer/Akoya QPTIFF"), + ], +) +def test_bad_inputs_raise(tmp_path: Path, files: list[str], match: str) -> None: + for name in files: + if name.endswith(".tif"): + tifffile.imwrite(tmp_path / name, np.zeros((8, 8), dtype=np.uint8)) + else: + write_qptiff(tmp_path / name, biomarkers=["DAPI"]) + + with pytest.raises((FileNotFoundError, ValueError), match=match): + phenocycler(tmp_path / files[0] if files and files[0].endswith(".tif") else tmp_path) + + +def test_regions_and_coordinate_systems(tmp_path: Path) -> None: + write_qptiff(tmp_path / "A_Scan1.qptiff", biomarkers=["DAPI"], stage_position_cm=(1.0, 2.0)) + write_qptiff(tmp_path / "B_Scan1.qptiff", biomarkers=["DAPI"], stage_position_cm=(1.01, 2.0)) + write_qptiff(tmp_path / "C_Scan1.qptiff", biomarkers=["DAPI"], slide_id="SLIDE2") + write_qptiff(tmp_path / "C_Scan2.qptiff", biomarkers=["DAPI"], slide_id="SLIDE2") + + sdata = phenocycler(tmp_path, regions=["A", "B"], image_models_kwargs=_FAST) + + # scans of the same slide are placed at their stage positions + assert sdata.coordinate_systems == ["global"] + assert get_transformation(sdata["A_image"], get_all=True) == {"global": Translation([0, 0], axes=("x", "y"))} + offset = get_transformation(sdata["B_image"], "global") + assert isinstance(offset, Translation) + np.testing.assert_allclose(offset.translation, [100 / 0.5069, 0], rtol=1e-3) + # unrelated scans get their own coordinate systems; `_Scan` is kept if it tells regions apart + sdata = phenocycler(tmp_path, regions=["A", "C_Scan1", "C_Scan2"], image_models_kwargs=_FAST) + assert sorted(sdata.coordinate_systems) == ["A", "C_Scan1", "C_Scan2"] + assert get_transformation(sdata["C_Scan1_image"], get_all=True) == {"C_Scan1": Identity()} + + +def test_read_all_bit_depths(tmp_path: Path) -> None: + for suffix, dtype in ((".qptiff", "uint8"), (".er.qptiff", "uint16")): + write_qptiff( + tmp_path / f"S{suffix}", biomarkers=["DAPI"], shape=(64, 80), dtype=dtype, stage_position_cm=(1, 2) + ) + write_qptiff(tmp_path / "T.qptiff", biomarkers=["DAPI"], shape=(64, 80), stage_position_cm=(1.1, 2)) + mask = np.zeros((32, 40), dtype=np.uint32) + tifffile.imwrite(tmp_path / "mask.tif", mask) + + assert sorted(phenocycler(tmp_path, image_models_kwargs=_FAST).images) == ["S_image", "T_image"] + sdata = phenocycler( + tmp_path / "S.er.qptiff", labels=tmp_path / "mask.tif", read_all_bit_depths=True, image_models_kwargs=_FAST + ) + + assert sorted(sdata.images) == ["S_8bit_image", "S_image"] + assert _level0(sdata, "S_image").dtype == np.uint16 + assert _level0(sdata, "S_8bit_image").dtype == np.uint8 + assert sdata.attrs["akoya"]["regions"]["S_8bit"]["path"].endswith("S.qptiff") + image_transformations = get_transformation(sdata["S_image"], get_all=True) + assert get_transformation(sdata["S_8bit_image"], get_all=True) == image_transformations + # the mask is scaled onto the extended-range region only + labels_transformation = get_transformation(sdata["S_labels"], "global") + assert isinstance(labels_transformation, Sequence) + assert labels_transformation.transformations == [Scale([2, 2], axes=("y", "x")), image_transformations["global"]] + with pytest.raises(ValueError, match=r"Unknown region\(s\) \['S_8bit'\]"): + phenocycler(tmp_path, labels={"S_8bit": tmp_path / "mask.tif"}, read_all_bit_depths=True) + write_qptiff(tmp_path / "S_8bit.qptiff", biomarkers=["DAPI"]) + with pytest.raises(ValueError, match="S_8bit"): + phenocycler(tmp_path, read_all_bit_depths=True) + + +def test_ignores_acquisition_files(tmp_path: Path) -> None: + images = tmp_path / "raw" / "images" + acquisition = images / "phenocycler" + acquisition.mkdir(parents=True) + write_qptiff(images / "S_Scan1.QPTIFF", biomarkers=["DAPI", "CD3"]) + for name in ("phenocycler/S_Scan1_Cycle1.qptiff", "phenocycler/S.qptiff.raw", "S_Scan1.raw.qptiff"): + write_qptiff(images / name, biomarkers=["DAPI"]) + + assert list(phenocycler(tmp_path, image_models_kwargs=_FAST).images) == ["S_image"] + + +def test_sidecars_hubmap_layout(tmp_path: Path) -> None: + raw_images = tmp_path / "raw" / "images" + processed = tmp_path / "lab_processed" / "images" + raw_images.mkdir(parents=True) + processed.mkdir(parents=True) + write_qptiff(raw_images / "S_Scan1.qptiff", biomarkers=["DAPI", "CD3"]) + wells = [ + {"wellName": "A1", "items": [{"markerName": "DAPI"}]}, + {"wellName": "A2", "items": [{"markerName": "DAPI"}, {"markerName": "CD3"}, {"markerName": "--"}]}, + ] + (raw_images / "S_Scan1.xpd").write_text(json.dumps({"wells": wells}), encoding="utf-8-sig") + (processed / "S_Scan1.qptiff.channels.csv").write_text("Channel,Marker\n1,DAPI\n2,CD3\n", encoding="utf-8-sig") + + # DAPI repeats per cycle and blanks do not count as channels, so the sidecars match the two pages + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + sdata = phenocycler(tmp_path, image_models_kwargs=_FAST) + + raw = sdata.attrs["akoya"]["regions"]["S"]["raw"] + assert raw["xpd"] == {"wells": wells} + assert raw["channels_csv"] == [{"Channel": "1", "Marker": "DAPI"}, {"Channel": "2", "Marker": "CD3"}] + + +@pytest.mark.parametrize( + ("name", "content", "match"), + [ + ("S_Scan1.xpd", json.dumps({"wells": [{"items": [{"markerName": "CD4"}]}]}).encode(), "lists 1 distinct"), + ("S_Scan1.qptiff.channels.csv", b"Marker\nDAPI\nCD3\nCD8\n", "channels.csv lists 3 distinct"), + ("S_Scan1.xpd", b"not json", "not valid JSON"), + ("S_Scan1.xpd", json.dumps({"wells": "x"}).encode(), "no readable markers"), + ("S_Scan1.qptiff.channels.csv", "Marker\nDAPI\n".encode("utf-16"), "Ignoring"), + ], +) +def test_bad_sidecars_warn(tmp_path: Path, name: str, content: bytes, match: str) -> None: + data = write_qptiff(tmp_path / "S_Scan1.qptiff", biomarkers=["DAPI", "CD3"]) + (tmp_path / name).write_bytes(content) + + with pytest.warns(UserWarning, match=match): + sdata = phenocycler(tmp_path / "S_Scan1.qptiff", image_models_kwargs=_FAST) + + # the channel names always come from the qptiff + assert get_channel_names(sdata.images["S_image"]) == ["DAPI", "CD3"] + np.testing.assert_array_equal(_level0(sdata, "S_image").values, data) + + +def test_labels(tmp_path: Path) -> None: + for name in ("Slide 1", "Slide 2"): + write_qptiff(tmp_path / f"{name}.qptiff", biomarkers=["DAPI"]) + mask = _write_mask(tmp_path / "mask.tif") + _write_mask(tmp_path / "small.tif", shape=(200, 300)) + + # keys may be raw names; keys of unselected regions are ignored; other shapes warn + with pytest.warns(UserWarning, match=r"\(200, 300\).*\(256, 320\).*not align"): + sdata = phenocycler( + tmp_path, + regions=["Slide 1"], + labels={"Slide 1": tmp_path / "small.tif", "Slide_2": tmp_path / "mask.tif"}, + image_models_kwargs=_FAST, + ) + assert list(sdata.labels) == ["Slide_1_labels"] + assert get_transformation(sdata["Slide_1_labels"], get_all=True) == {"global": Identity()} + assert sdata.attrs["akoya"]["labels"] == {"Slide_1": str(tmp_path / "small.tif")} + + sdata = phenocycler(tmp_path / "Slide 1.qptiff", labels=tmp_path / "mask.tif", labels_models_kwargs=_FAST) + np.testing.assert_array_equal(sdata.labels["Slide_1_labels"]["scale0"]["image"].values, mask) + assert not sdata.tables + + +@pytest.mark.parametrize( + ("labels", "match"), + [ + ("mask.tif", "mapping"), + ({"Z": "mask.tif"}, "'Z'"), + ({"Slide 1": "mask.tif", "Slide_1": "mask.tif"}, "same region"), + ], +) +def test_bad_labels_raise(tmp_path: Path, labels: str | dict[str, str], match: str) -> None: + for name in ("Slide 1", "Slide 2"): + write_qptiff(tmp_path / f"{name}.qptiff", biomarkers=["DAPI"]) + _write_mask(tmp_path / "mask.tif") + resolved = tmp_path / labels if isinstance(labels, str) else {k: tmp_path / v for k, v in labels.items()} + + with pytest.raises(ValueError, match=match): + phenocycler(tmp_path, labels=resolved) + + +def test_cli(runner: CliRunner, tmp_path: Path) -> None: + write_qptiff(tmp_path / "A_[1,2].qptiff", biomarkers=["DAPI", "HLA-A"]) + write_qptiff(tmp_path / "B.qptiff", biomarkers=["DAPI", "CD3"]) + _write_mask(tmp_path / "mask.tif") + args = ["--input", str(tmp_path), "--image-models-kwargs", json.dumps(_FAST)] + + # region ids and channel names may hold commas + result = runner.invoke( + phenocycler_wrapper, + [ + *args, + "-o", + str(tmp_path / "a.zarr"), + "--regions", + "A_[1,2]", + *("--channel-names", "DAPI", "--channel-names", "HLA-A,B,C"), + *("--labels", json.dumps({"A_[1,2]": str(tmp_path / "mask.tif")})), + ], + ) + assert result.exit_code == 0, result.output + sdata = read_zarr(tmp_path / "a.zarr") + assert get_channel_names(sdata.images["A__1_2__image"]) == ["DAPI", "HLA-A,B,C"] + assert list(sdata.labels) == ["A__1_2__labels"] + + for option, value, match in ( + ("--labels", '{"B": 1}', "labels"), + ("--labels", "{not json", "labels"), + ("--labels-models-kwargs", "null", "JSON object"), + ): + result = runner.invoke(phenocycler_wrapper, [*args, "-o", str(tmp_path / "b.zarr"), option, value]) + assert result.exit_code == 2 + assert match in result.output