various fixes and additions for IO

.pre-commit-config.yaml

@@ -7,7 +7,7 @@ default_stages:
 minimum_pre_commit_version: 2.16.0
 repos:
     - repo: https://github.com/psf/black
-      rev: 22.12.0
+      rev: 23.1.0
       hooks:
           - id: black
     - repo: https://github.com/pre-commit/mirrors-prettier
@@ -19,11 +19,11 @@ repos:
       hooks:
           - id: blacken-docs
     - repo: https://github.com/PyCQA/isort
-      rev: 5.11.4
+      rev: 5.12.0
       hooks:
           - id: isort
     - repo: https://github.com/pre-commit/mirrors-mypy
-      rev: v0.991
+      rev: v1.0.1
       hooks:
           - id: mypy
             additional_dependencies: [numpy==1.24.0]
@@ -50,7 +50,7 @@ repos:
           - id: trailing-whitespace
           - id: check-case-conflict
     - repo: https://github.com/PyCQA/autoflake
-      rev: v2.0.0
+      rev: v2.0.1
       hooks:
           - id: autoflake
             args:

pyproject.toml

@@ -26,6 +26,7 @@ dependencies = [
     "joblib",
     "imagecodecs",
     "dask-image",
+    "pyarrow",
 ]
 
 [project.optional-dependencies]

src/spatialdata_io/__init__.py

@@ -1,13 +1,17 @@
 from importlib.metadata import version
 
 from spatialdata_io.readers.cosmx import cosmx
+from spatialdata_io.readers.mcmicro import mcmicro
+from spatialdata_io.readers.steinbock import steinbock
 from spatialdata_io.readers.visium import visium
 from spatialdata_io.readers.xenium import xenium
 
 __all__ = [
     "visium",
     "xenium",
     "cosmx",
+    "mcmicro",
+    "steinbock",
 ]
 
 __version__ = version("spatialdata-io")

src/spatialdata_io/_constants/_constants.py

@@ -57,7 +57,7 @@ class XeniumKeys(ModeEnum):
     CELL_METADATA_FILE = "cells.parquet"
     CELL_X = "x_centroid"
     CELL_Y = "y_centroid"
-    CELL_AREA = 'cell_area'
+    CELL_AREA = "cell_area"
 
     # morphology iamges
     MORPHOLOGY_MIP_FILE = "morphology_mip.ome.tif"
@@ -85,3 +85,37 @@ class VisiumKeys(ModeEnum):
     SPOTS_FILE = "spatial/tissue_positions.csv"
     SPOTS_X = "pxl_row_in_fullres"
     SPOTS_Y = "pxl_col_in_fullres"
+
+
+@unique
+class SteinbockKeys(ModeEnum):
+    """Keys for *Steinbock* formatted dataset."""
+
+    # files and directories
+    CELLS_FILE = "cells.h5ad"
+    DEEPCELL_MASKS_DIR = "masks_deepcell"
+    ILASTIK_MASKS_DIR = "masks_ilastik"
+    IMAGES_DIR = "ome"
+
+    # suffixes for images and labels
+    IMAGE_SUFFIX = ".ome.tiff"
+    LABEL_SUFFIX = ".tiff"
+
+
+@unique
+class McmicroKeys(ModeEnum):
+    """Keys for *Mcmicro* formatted dataset."""
+
+    # files and directories
+    CELL_FEATURES_SUFFIX = "--unmicst_cell.csv"
+    QUANTIFICATION_DIR = "quantification"
+    MARKERS_FILE = "markers.csv"
+    IMAGES_DIR = "registration"
+    IMAGE_SUFFIX = ".ome.tif"
+    LABELS_DIR = "segmentation"
+    LABELS_PREFIX = "unmicst-"
+
+    # metadata
+    COORDS_X = "X_centroid"
+    COORDS_Y = "Y_centroid"
+    INSTANCE_KEY = "CellID"

src/spatialdata_io/readers/cosmx.py

@@ -2,26 +2,29 @@
 
 import os
 import re
-import tempfile
 from collections.abc import Mapping
 from pathlib import Path
 from types import MappingProxyType
 from typing import Any, Optional
-import pyarrow as pa
-import pyarrow.parquet as pq
 
+import dask.array as da
 import numpy as np
 import pandas as pd
+import pyarrow as pa
 from anndata import AnnData
-from dask_image.imread import imread
-import dask.array as da
 from dask.dataframe.core import DataFrame as DaskDataFrame
+from dask_image.imread import imread
 from scipy.sparse import csr_matrix
 
 # from spatialdata._core.core_utils import xy_cs
 from skimage.transform import estimate_transform
 from spatialdata import SpatialData
-from spatialdata._core.models import Image2DModel, Labels2DModel, TableModel, PointsModel
+from spatialdata._core.models import (
+    Image2DModel,
+    Labels2DModel,
+    PointsModel,
+    TableModel,
+)
 
 # from spatialdata._core.ngff.ngff_coordinate_system import NgffAxis  # , CoordinateSystem
 from spatialdata._core.transformations import Affine, Identity
@@ -41,7 +44,6 @@
 def cosmx(
     path: str | Path,
     dataset_id: Optional[str] = None,
-    # shape_size: float | int = 1,
     transcripts: bool = True,
     imread_kwargs: Mapping[str, Any] = MappingProxyType({}),
     image_models_kwargs: Mapping[str, Any] = MappingProxyType({}),
@@ -67,8 +69,8 @@ def cosmx(
         Path to the root directory containing *Nanostring* files.
     dataset_id
         Name of the dataset.
-    shape_size
-        Size of the shape to be used for the centroids of the labels.
+    transcripts
+        Whether to also read in transcripts information.
     imread_kwargs
         Keyword arguments passed to :func:`dask_image.imread.imread`.
     image_models_kwargs
@@ -141,12 +143,6 @@ def cosmx(
 
     fovs_counts = list(map(str, adata.obs.fov.astype(int).unique()))
 
-    # TODO(giovp): uncomment once transform is ready
-    # input_cs = CoordinateSystem("cxy", axes=[c_axis, y_axis, x_axis])
-    # input_cs_labels = CoordinateSystem("cxy", axes=[y_axis, x_axis])
-    # output_cs = CoordinateSystem("global", axes=[c_axis, y_axis, x_axis])
-    # output_cs_labels = CoordinateSystem("global", axes=[y_axis, x_axis])
-
     affine_transforms_to_global = {}
 
     for fov in fovs_counts:
@@ -163,7 +159,10 @@ def cosmx(
 
     table.obsm["global"] = table.obs[[CosmxKeys.X_GLOBAL_CELL, CosmxKeys.Y_GLOBAL_CELL]].to_numpy()
     table.obsm["spatial"] = table.obs[[CosmxKeys.X_LOCAL_CELL, CosmxKeys.Y_LOCAL_CELL]].to_numpy()
-    table.obs.drop(columns=[CosmxKeys.X_LOCAL_CELL, CosmxKeys.Y_LOCAL_CELL, CosmxKeys.X_GLOBAL_CELL, CosmxKeys.Y_GLOBAL_CELL], inplace=True)
+    table.obs.drop(
+        columns=[CosmxKeys.X_LOCAL_CELL, CosmxKeys.Y_LOCAL_CELL, CosmxKeys.X_GLOBAL_CELL, CosmxKeys.Y_GLOBAL_CELL],
+        inplace=True,
+    )
 
     # prepare to read images and labels
     file_extensions = (".jpg", ".png", ".jpeg", ".tif", ".tiff")
@@ -200,7 +199,6 @@ def cosmx(
                 flipped_im = da.flip(im, axis=0)
                 parsed_im = Image2DModel.parse(
                     flipped_im,
-                    name=fov,
                     transformations={
                         fov: Identity(),
                         "global": aff,
@@ -224,7 +222,6 @@ def cosmx(
                 flipped_la = da.flip(la, axis=0)
                 parsed_la = Labels2DModel.parse(
                     flipped_la,
-                    name=fov,
                     transformations={
                         fov: Identity(),
                         "global": aff,
@@ -239,33 +236,27 @@ def cosmx(
 
     points: dict[str, DaskDataFrame] = {}
     if transcripts:
-        # let's convert the .csv to .parquet and let's read it with pyarrow.parquet for faster subsetting
-        with tempfile.TemporaryDirectory() as tmpdir:
-            print("converting .csv to .parquet... ", end="")
-            assert transcripts_file is not None
-            transcripts_data = pd.read_csv(path / transcripts_file, header=0)
-            transcripts_data.to_parquet(Path(tmpdir) / "transcripts.parquet")
-            print("done")
-
-            ptable = pq.read_table(Path(tmpdir) / "transcripts.parquet")
-            for fov in fovs_counts:
-                aff = affine_transforms_to_global[fov]
-                sub_table = ptable.filter(pa.compute.equal(ptable.column(CosmxKeys.FOV), int(fov))).to_pandas()
-                sub_table[CosmxKeys.INSTANCE_KEY] = sub_table[CosmxKeys.INSTANCE_KEY].astype('category')
-                # we rename z because we want to treat the data as 2d
-                sub_table.rename(columns={'z': 'z_raw'}, inplace=True)
-                points[fov] = PointsModel.parse(
-                    sub_table,
-                    coordinates={"x": CosmxKeys.X_LOCAL_TRANSCRIPT, "y": CosmxKeys.Y_LOCAL_TRANSCRIPT},
-                    feature_key=CosmxKeys.TARGET_OF_TRANSCRIPT,
-                    instance_key=CosmxKeys.INSTANCE_KEY,
-                    transformations={
-                        fov: Identity(),
-                        "global": aff,
-                        "global_only_labels": aff,
-                    },
-                )
-
+        assert transcripts_file is not None
+        from pyarrow.csv import read_csv
+
+        ptable = read_csv(path / transcripts_file)  # , header=0)
+        for fov in fovs_counts:
+            aff = affine_transforms_to_global[fov]
+            sub_table = ptable.filter(pa.compute.equal(ptable.column(CosmxKeys.FOV), int(fov))).to_pandas()
+            sub_table[CosmxKeys.INSTANCE_KEY] = sub_table[CosmxKeys.INSTANCE_KEY].astype("category")
+            # we rename z because we want to treat the data as 2d
+            sub_table.rename(columns={"z": "z_raw"}, inplace=True)
+            points[fov] = PointsModel.parse(
+                sub_table,
+                coordinates={"x": CosmxKeys.X_LOCAL_TRANSCRIPT, "y": CosmxKeys.Y_LOCAL_TRANSCRIPT},
+                feature_key=CosmxKeys.TARGET_OF_TRANSCRIPT,
+                instance_key=CosmxKeys.INSTANCE_KEY,
+                transformations={
+                    fov: Identity(),
+                    "global": aff,
+                    "global_only_labels": aff,
+                },
+            )
 
     # TODO: what to do with fov file?
     # if fov_file is not None:

src/spatialdata_io/readers/mcmicro.py

@@ -0,0 +1,130 @@
+from __future__ import annotations
+
+import os
+from collections.abc import Mapping
+from pathlib import Path
+from types import MappingProxyType
+from typing import Any, Union
+
+import numpy as np
+import pandas as pd
+from anndata import AnnData
+from dask_image.imread import imread
+from multiscale_spatial_image.multiscale_spatial_image import MultiscaleSpatialImage
+from spatial_image import SpatialImage
+from spatialdata import Image2DModel, Labels2DModel, SpatialData, TableModel
+
+from spatialdata_io._constants._constants import McmicroKeys
+
+__all__ = ["mcmicro"]
+
+
+def mcmicro(
+    path: str | Path,
+    dataset_id: str,
+    imread_kwargs: Mapping[str, Any] = MappingProxyType({}),
+    image_models_kwargs: Mapping[str, Any] = MappingProxyType({}),
+) -> SpatialData:
+    """
+    Read a *Mcmicro* output into a SpatialData object.
+
+    .. seealso::
+
+        - `Mcmicro pipeline  <https://mcmicro.org/>`_.
+
+    Parameters
+    ----------
+    path
+        Path to the dataset.
+    dataset_id
+        Dataset identifier.
+    imread_kwargs
+        Keyword arguments to pass to the image reader.
+    image_models_kwargs
+        Keyword arguments to pass to the image models.
+
+    Returns
+    -------
+    :class:`spatialdata.SpatialData`
+    """
+    path = Path(path)
+
+    samples = os.listdir(path / McmicroKeys.IMAGES_DIR)
+    if len(samples) > 1:
+        raise ValueError("Only one sample per dataset is supported.")
+    if (dataset_id + McmicroKeys.IMAGE_SUFFIX) not in samples:
+        raise ValueError("Dataset id is not consistent with sample name.")
+
+    images = {}
+    images[dataset_id] = _get_images(
+        path,
+        dataset_id,
+        imread_kwargs,
+        image_models_kwargs,
+    )
+    labels = {}
+    labels[dataset_id] = _get_labels(
+        path,
+        dataset_id,
+        "cell",
+        imread_kwargs,
+        image_models_kwargs,
+    )
+    labels[f"{dataset_id}_nuclei"] = _get_labels(
+        path,
+        dataset_id,
+        "nuclei",
+        imread_kwargs,
+        image_models_kwargs,
+    )
+
+    table = _get_table(path, dataset_id)
+
+    return SpatialData(images=images, labels=labels, table=table)
+
+
+def _get_images(
+    path: Path,
+    sample: str,
+    imread_kwargs: Mapping[str, Any] = MappingProxyType({}),
+    image_models_kwargs: Mapping[str, Any] = MappingProxyType({}),
+) -> Union[SpatialImage, MultiscaleSpatialImage]:
+    image = imread(path / McmicroKeys.IMAGES_DIR / f"{sample}{McmicroKeys.IMAGE_SUFFIX}", **imread_kwargs)
+    return Image2DModel.parse(image, **image_models_kwargs)
+
+
+def _get_labels(
+    path: Path,
+    sample: str,
+    labels_kind: str,
+    imread_kwargs: Mapping[str, Any] = MappingProxyType({}),
+    image_models_kwargs: Mapping[str, Any] = MappingProxyType({}),
+) -> Union[SpatialImage, MultiscaleSpatialImage]:
+    image = imread(
+        path
+        / McmicroKeys.LABELS_DIR
+        / f"{McmicroKeys.LABELS_PREFIX}{sample}"
+        / f"{labels_kind}{McmicroKeys.IMAGE_SUFFIX}",
+        **imread_kwargs,
+    ).squeeze()
+    return Labels2DModel.parse(image, **image_models_kwargs)
+
+
+def _get_table(
+    path: Path,
+    sample: str,
+) -> AnnData:
+    table = pd.read_csv(path / McmicroKeys.QUANTIFICATION_DIR / f"{sample}{McmicroKeys.CELL_FEATURES_SUFFIX}")
+    markers = pd.read_csv(path / McmicroKeys.MARKERS_FILE)
+    markers.index = markers.marker_name
+    var = markers.marker_name.tolist()
+    coords = [McmicroKeys.COORDS_X.value, McmicroKeys.COORDS_Y.value]
+    adata = AnnData(
+        table[var].to_numpy(),
+        obs=table.drop(columns=var + coords),
+        var=markers,
+        obsm={"spatial": table[coords].to_numpy()},
+        dtype=np.float_,
+    )
+
+    return TableModel.parse(adata, region=sample, instance_key=McmicroKeys.INSTANCE_KEY.value)