Merge pull request #335 from MannLabs/fix_bugs_xenium

Fix bugs xenium

Author: sophiamaedler

src/scportrait/pipeline/featurization.py

@@ -125,7 +125,7 @@ def _detect_automatic_inference_device(self) -> str:
 
         if torch.cuda.is_available():
             inference_device = "cuda"
-        if torch.backends.mps.is_available():
+        elif torch.backends.mps.is_available():
             inference_device = torch.device("mps")
         else:
             inference_device = "cpu"

src/scportrait/pipeline/project.py

@@ -15,6 +15,7 @@
 import re
 import shutil
 import tempfile
+import warnings
 from pathlib import PosixPath
 from typing import TYPE_CHECKING, Literal
 
@@ -198,7 +199,7 @@ def __init__(
         if not os.path.isdir(self.project_location):
             os.makedirs(self.project_location)
         else:
-            Warning("There is already a directory in the location path")
+            warnings.warn("There is already a directory in the location path", stacklevel=2)
 
         # === setup sdata reader/writer ===
         self.filehandler = sdata_filehandler(
@@ -499,8 +500,9 @@ def _check_image_dtype(self, image: np.ndarray) -> None:
         """
 
         if not image.dtype == self.DEFAULT_IMAGE_DTYPE:
-            Warning(
-                f"Image dtype is not {self.DEFAULT_IMAGE_DTYPE} but insteadt {image.dtype}. The workflow expects images to be of dtype {self.DEFAULT_IMAGE_DTYPE}. Proceeding with the incorrect dtype can lead to unexpected results."
+            warnings.warn(
+                f"Image dtype is not {self.DEFAULT_IMAGE_DTYPE} but insteadt {image.dtype}. The workflow expects images to be of dtype {self.DEFAULT_IMAGE_DTYPE}. Proceeding with the incorrect dtype can lead to unexpected results.",
+                stacklevel=2,
             )
             self.log(
                 f"Image dtype is not {self.DEFAULT_IMAGE_DTYPE} but insteadt {image.dtype}. The workflow expects images to be of dtype {self.DEFAULT_IMAGE_DTYPE}. Proceeding with the incorrect dtype can lead to unexpected results."
@@ -638,7 +640,10 @@ def close_interactive_viewer(self):
 
     def _check_for_interactive_session(self):
         if self.interactive is not None:
-            Warning("Interactive viewer is still open. Will automatically close before proceeding with processing.")
+            warnings.warn(
+                "Interactive viewer is still open. Will automatically close before proceeding with processing.",
+                stacklevel=2,
+            )
             self.close_interactive_viewer()
 
     #### Functions to visualize results ####
@@ -1452,8 +1457,9 @@ def load_input_from_sdata(
             table = sdata_input[table_elem]
             rename_columns = {}
             if self.DEFAULT_CELL_ID_NAME in table.obs:
-                Warning(
-                    f"Column {self.DEFAULT_CELL_ID_NAME} already exists in table. Renaming to `f{self.DEFAULT_CELL_ID_NAME}_orig` to preserve compatibility with scPortrait workflow."
+                warnings.warn(
+                    f"Column {self.DEFAULT_CELL_ID_NAME} already exists in table. Renaming to `f{self.DEFAULT_CELL_ID_NAME}_orig` to preserve compatibility with scPortrait workflow.",
+                    stacklevel=2,
                 )
                 rename_columns[self.DEFAULT_CELL_ID_NAME] = f"{self.DEFAULT_CELL_ID_NAME}_orig"
                 self.log(
@@ -1495,11 +1501,31 @@ def load_input_from_sdata(
 
             self.filehandler._add_centers(segmentation_label=self.cyto_seg_name)
 
+        # read the sdata object from file to ensure we have access to all newly written elements
+        sdata = SpatialData.read(self.sdata_path)
+
         # ensure that if both an nucleus and cytosol segmentation mask are loaded that they match
         if self.nuc_seg_status and self.cyto_seg_status:
-            ids_nuc = set(sdata[f"{self.DEFAULT_CENTERS_NAME}_{self.nuc_seg_name}"].index.values)
-            ids_cyto = set(sdata[f"{self.DEFAULT_CENTERS_NAME}_{self.cyto_seg_name}"].index.values)
-            assert ids_nuc == ids_cyto, "Nucleus and cytosol segmentation masks do not match."
+            ids_nuc = set(sdata["centers_seg_all_nucleus"].index.compute().values)
+            ids_cyto = set(sdata["centers_seg_all_cytosol"].index.compute().values)
+
+            if ids_nuc.issubset(ids_cyto):
+                if not ids_nuc == ids_cyto:
+                    warnings.warn(
+                        "Nucleus segmentation mask is a subset of cytosol segmentation mask, but they do not match exactly. \n This means that cells exist which do not have a nucleus mask but do have a cytosol mask. \n Please be careful when configuring your extraction workflow.",
+                        stacklevel=2,
+                    )
+
+            elif ids_cyto.issubset(ids_nuc):
+                if not ids_cyto == ids_nuc:
+                    warnings.warn(
+                        "Cytosol segmentation mask is a subset of nucleus segmentation mask, but they do not match exactly. \n  This means cells exist which do not have a cytosol mask but do have a nucleus mask. \n Please be careful when configuring your extraction workflow.",
+                        stacklevel=2,
+                    )
+            else:
+                raise ValueError(
+                    "Nucleus and cytosol segmentation masks do not match. This is unexpected and should be investigated."
+                )
 
         self.get_project_status()
 

src/scportrait/tools/sdata/write/_write.py

@@ -70,12 +70,22 @@ def image(
             # fix until #https://github.com/scverse/spatialdata/issues/528 is resolved
             Image2DModel().validate(image)
 
-            if channel_names is not None:
-                warnings.warn(
-                    "Channel names are ignored when passing a single scale image in the DataArray format. Channel names are read directly from the DataArray.",
-                    stacklevel=2,
-                )
+            # read channel names from the DataArray if not provided
+            if channel_names is None:
+                channel_names = image.coords["c"].values.tolist()
+            else:
+                if len(channel_names) != image.shape[0]:
+                    raise ValueError(
+                        f"Number of channel names ({len(channel_names)}) does not match the number of channels in the image ({image.shape[0]})."
+                    )
+                channel_names_old = image.coords["c"].values.tolist()
+                if channel_names_old != channel_names:
+                    warnings.warn(
+                        f"Channel names in the DataArray ({channel_names_old}) do not match the provided channel names ({channel_names}). The DataArray will be updated with the provided channel names.",
+                        stacklevel=2,
+                    )
 
+            # if so first validate the model since this means we are getting the image from a spatialdata object already
             if chunks is not None:
                 warnings.warn(
                     "Chunks are ignored when passing a single scale image in the DataArray format. Chunks are read directly from the DataArray.",
@@ -84,6 +94,7 @@ def image(
 
             image = Image2DModel.parse(
                 image,
+                c_coords=channel_names,
                 scale_factors=scale_factors,
                 rgb=rgb,
             )