major API simplification (WIP)

Author: quentinblampey

CHANGELOG.md

@@ -11,6 +11,7 @@
 
 ### Changed
 - Import submodules in init (segmentation, io, utils)
+- API simplification in progress (new API + tutorial comming soon)
 
 ## [1.1.2] - 2024-07-24
 

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "sopa"
-version = "1.1.2"
+version = "1.1.3"
 description = "Spatial-omics pipeline and analysis"
 documentation = "https://gustaveroussy.github.io/sopa"
 homepage = "https://gustaveroussy.github.io/sopa"
@@ -76,7 +76,7 @@ testpaths = ["tests"]
 python_files = "test_*.py"
 
 [tool.black]
-line-length = 100
+line-length = 120
 include = '\.pyi?$'
 exclude = '''
 /(

sopa/__init__.py

@@ -12,3 +12,6 @@
     from . import utils
     from . import io
     from . import segmentation
+
+    from .segmentation import tissue_segmentation
+    from ._sdata import get_spatial_image, get_spatial_element, to_intrinsic

sopa/_constants.py

@@ -33,6 +33,13 @@ class SopaKeys:
     GEOMETRY_COUNT = "n_components"
 
 
+class SopaAttrs:
+    CELL_SEGMENTATION = "for_cell_segmentation"
+    TISSUE_SEGMENTATION = "for_tissue_segmentation"
+    BINS_AGGREGATION = "for_bins_aggregation"
+    GENE_COLUMN = "feature_key"
+
+
 VALID_DIMENSIONS = ("c", "y", "x")
 LOW_AVERAGE_COUNT = 0.01
 EPS = 1e-5

sopa/_sdata.py

@@ -12,7 +12,7 @@
 from spatialdata.transformations import Identity, get_transformation, set_transformation
 from xarray import DataArray
 
-from ._constants import SopaKeys
+from ._constants import SopaAttrs, SopaKeys
 
 log = logging.getLogger(__name__)
 
@@ -40,9 +40,7 @@ def get_boundaries(
         if res is not None:
             return res
 
-    error_message = (
-        "sdata object has no valid segmentation boundary. Consider running Sopa segmentation first."
-    )
+    error_message = "sdata object has no valid segmentation boundary. Consider running Sopa segmentation first."
 
     if not warn:
         raise ValueError(error_message)
@@ -51,17 +49,13 @@ def get_boundaries(
     return (None, None) if return_key else None
 
 
-def _try_get_boundaries(
-    sdata: SpatialData, shapes_key: str, return_key: bool
-) -> gpd.GeoDataFrame | None:
+def _try_get_boundaries(sdata: SpatialData, shapes_key: str, return_key: bool) -> gpd.GeoDataFrame | None:
     """Try to get a cell boundaries for a given `shapes_key`"""
     if shapes_key in sdata.shapes:
         return (shapes_key, sdata[shapes_key]) if return_key else sdata[shapes_key]
 
 
-def get_intrinsic_cs(
-    sdata: SpatialData, element: SpatialElement | str, name: str | None = None
-) -> str:
+def get_intrinsic_cs(sdata: SpatialData, element: SpatialElement | str, name: str | None = None) -> str:
     """Gets the name of the intrinsic coordinate system of an element
 
     Args:
@@ -86,9 +80,7 @@ def get_intrinsic_cs(
     return name
 
 
-def to_intrinsic(
-    sdata: SpatialData, element: SpatialElement | str, element_cs: SpatialElement | str
-) -> SpatialElement:
+def to_intrinsic(sdata: SpatialData, element: SpatialElement | str, element_cs: SpatialElement | str) -> SpatialElement:
     """Transforms a `SpatialElement` into the intrinsic coordinate system of another `SpatialElement`
 
     Args:
@@ -105,32 +97,6 @@ def to_intrinsic(
     return sdata.transform_element_to_coordinate_system(element, cs)
 
 
-def get_key(sdata: SpatialData, attr: str, key: str | None = None):
-    if key is not None:
-        return key
-
-    elements = getattr(sdata, attr)
-
-    if not len(elements):
-        return None
-
-    assert (
-        len(elements) == 1
-    ), f"Trying to get an element key of `sdata.{attr}`, but it contains multiple values and no dict key was provided"
-
-    return next(iter(elements.keys()))
-
-
-def get_element(sdata: SpatialData, attr: str, key: str | None = None):
-    key = get_key(sdata, attr, key)
-    return sdata[key] if key is not None else None
-
-
-def get_item(sdata: SpatialData, attr: str, key: str | None = None):
-    key = get_key(sdata, attr, key)
-    return key, sdata[key] if key is not None else None
-
-
 def get_intensities(sdata: SpatialData) -> pd.DataFrame | None:
     """Gets the intensity dataframe of shape `n_obs x n_channels`"""
     assert SopaKeys.TABLE in sdata.tables, f"No '{SopaKeys.TABLE}' found in sdata.tables"
@@ -155,35 +121,87 @@ def iter_scales(image: DataTree) -> Iterator[xr.DataArray]:
     Yields:
         Each scale (as a `xr.DataArray`)
     """
-    assert isinstance(
-        image, DataTree
-    ), f"Multiscale iteration is reserved for type DataTree. Found {type(image)}"
+    assert isinstance(image, DataTree), f"Multiscale iteration is reserved for type DataTree. Found {type(image)}"
 
     for scale in image:
         yield next(iter(image[scale].values()))
 
 
+def get_spatial_element(
+    element_dict: dict[str, SpatialElement],
+    key: str | None = None,
+    valid_attr: str | None = None,
+    return_key: bool = False,
+    as_spatial_image: bool = False,
+) -> SpatialElement | tuple[str, SpatialElement]:
+    """Gets an element from a SpatialData object.
+
+    Args:
+        sdata: SpatialData object.
+        key: Optional element key. If `None`, returns the only element (if only one), or tries to find an element with `valid_attr`.
+        return_key: Whether to also return the key of the element.
+        valid_attr: Attribute that the element must have to be considered valid.
+        as_spatial_image: Whether to return the element as a `SpatialImage` (if it is a `DataTree`)
+
+    Returns:
+        If `return_key` is False, only the element is returned, else a tuple `(element_key, element)`
+    """
+    assert len(element_dict), "No spatial element was found in the dict."
+
+    if key is not None:
+        return _return_element(element_dict, key, return_key, as_spatial_image)
+
+    if len(element_dict) == 1:
+        key = next(iter(element_dict.keys()))
+
+        assert valid_attr is None or element_dict[key].attrs.get(
+            valid_attr, True
+        ), f"Element {key} is not valid for the attribute {valid_attr}."
+
+        return _return_element(element_dict, key, return_key, as_spatial_image)
+
+    assert valid_attr is not None, "Multiple elements found. Provide an element key."
+
+    keys = [key for key, element in element_dict.items() if element.attrs.get(valid_attr)]
+
+    assert len(keys) > 0, f"No element with the attribute {valid_attr}. Provide an element key."
+    assert len(keys) == 1, f"Multiple valid elements found: {keys}. Provide an element key."
+
+    return _return_element(element_dict, keys[0], return_key, as_spatial_image)
+
+
 def get_spatial_image(
-    sdata: SpatialData, key: str | None = None, return_key: bool = False
+    sdata: SpatialData,
+    key: str | None = None,
+    return_key: bool = False,
+    valid_attr: str = SopaAttrs.CELL_SEGMENTATION,
 ) -> DataArray | tuple[str, DataArray]:
     """Gets a DataArray from a SpatialData object (if the image has multiple scale, the `scale0` is returned)
 
     Args:
         sdata: SpatialData object.
-        key: Optional image key. If `None`, returns the only image (if only one), or raises an error.
+        key: Optional image key. If `None`, returns the only image (if only one), or tries to find an image with `valid_attr`.
         return_key: Whether to also return the key of the image.
+        valid_attr: Attribute that the image must have to be considered valid.
 
     Returns:
         If `return_key` is False, only the image is returned, else a tuple `(image_key, image)`
     """
-    key = get_key(sdata, "images", key)
+    return get_spatial_element(
+        sdata.images,
+        key=key,
+        valid_attr=valid_attr,
+        return_key=return_key,
+        as_spatial_image=True,
+    )
+
 
-    assert key is not None, "One image in `sdata.images` is required"
+def _return_element(
+    element_dict: dict[str, SpatialElement], key: str, return_key: bool, as_spatial_image: bool
+) -> SpatialElement | tuple[str, SpatialElement]:
+    element = element_dict[key]
 
-    image = sdata.images[key]
-    if isinstance(image, DataTree):
-        image = next(iter(image["scale0"].values()))
+    if as_spatial_image and isinstance(element, DataTree):
+        element = next(iter(element["scale0"].values()))
 
-    if return_key:
-        return key, image
-    return image
+    return (key, element) if return_key else element

sopa/annotation/tangram/run.py

@@ -124,9 +124,7 @@ def init_obsm(self, level: int):
         )
 
     def get_hard_labels(self, df: pd.DataFrame) -> pd.Series:
-        df = df.clip(
-            df.quantile(1 - self.clip_percentile), df.quantile(self.clip_percentile), axis=1
-        )
+        df = df.clip(df.quantile(1 - self.clip_percentile), df.quantile(self.clip_percentile), axis=1)
         df = (df - df.min()) / (df.max() - df.min())
         return df.idxmax(1)
 
@@ -138,9 +136,7 @@ def pp_adata(self, ad_sp_: AnnData, ad_sc_: AnnData, split: np.ndarray) -> AnnDa
         sc.pp.filter_genes(ad_sp_split, min_cells=1)
 
         # Calculate uniform density prior as 1/number_of_spots
-        ad_sp_split.obs["uniform_density"] = (
-            np.ones(ad_sp_split.X.shape[0]) / ad_sp_split.X.shape[0]
-        )
+        ad_sp_split.obs["uniform_density"] = np.ones(ad_sp_split.X.shape[0]) / ad_sp_split.X.shape[0]
 
         # Calculate rna_count_based density prior as % of rna molecule count
         rna_count_per_spot = np.array(ad_sp_split.X.sum(axis=1)).squeeze()
@@ -157,8 +153,7 @@ def pp_adata(self, ad_sp_: AnnData, ad_sc_: AnnData, split: np.ndarray) -> AnnDa
         )
 
         selection = list(
-            set(ad_sp_split.var_names[ad_sp_split.var.counts > 0])
-            & set(ad_sc_.var_names[ad_sc_.var.counts > 0])
+            set(ad_sp_split.var_names[ad_sp_split.var.counts > 0]) & set(ad_sc_.var_names[ad_sc_.var.counts > 0])
         )
 
         assert len(

sopa/cli/annotate.py

@@ -13,9 +13,7 @@
 def fluorescence(
     sdata_path: str = typer.Argument(help=SDATA_HELPER),
     marker_cell_dict: str = typer.Option(callback=ast.literal_eval),
-    cell_type_key: str = typer.Option(
-        "cell_type", help="Key added in `adata.obs` corresponding to the cell type"
-    ),
+    cell_type_key: str = typer.Option("cell_type", help="Key added in `adata.obs` corresponding to the cell type"),
 ):
     """Simple annotation based on fluorescence, where each provided channel corresponds to one cell type.
 
@@ -39,9 +37,7 @@ def fluorescence(
 @app_annotate.command()
 def tangram(
     sdata_path: str = typer.Argument(help=SDATA_HELPER),
-    sc_reference_path: str = typer.Option(
-        help="Path to the scRNAseq annotated reference, as a `.h5ad` file"
-    ),
+    sc_reference_path: str = typer.Option(help="Path to the scRNAseq annotated reference, as a `.h5ad` file"),
     cell_type_key: str = typer.Option(help="Key of `adata_ref.obs` containing the cell-types"),
     reference_preprocessing: str = typer.Option(
         None,

sopa/cli/app.py

@@ -28,9 +28,7 @@
     name="segmentation",
     help="Perform cell segmentation on patches. NB: for `baysor`, use directly the `baysor` command line.",
 )
-app.add_typer(
-    app_resolve, name="resolve", help="Resolve the segmentation conflicts over patches overlaps"
-)
+app.add_typer(app_resolve, name="resolve", help="Resolve the segmentation conflicts over patches overlaps")
 app.add_typer(
     app_patchify,
     name="patchify",
@@ -77,9 +75,7 @@ def read(
 
     io.standardize._check_can_write_zarr(sdata_path)
 
-    assert (
-        technology is not None or config_path is not None
-    ), "Provide the argument `--technology` or `--config-path`"
+    assert technology is not None or config_path is not None, "Provide the argument `--technology` or `--config-path`"
 
     if config_path is not None:
         assert not kwargs, "Provide either a path to a config, or some kwargs, but not both"
@@ -114,9 +110,7 @@ def crop(
         None,
         help="List of channel names to be displayed. Optional if there are already only 1 or 3 channels",
     ),
-    scale_factor: float = typer.Option(
-        10, help="Resize the image by this value (high value for a lower memory usage)"
-    ),
+    scale_factor: float = typer.Option(10, help="Resize the image by this value (high value for a lower memory usage)"),
     margin_ratio: float = typer.Option(
         0.1, help="Ratio of the image margin on the display (compared to the image size)"
     ),
@@ -163,16 +157,12 @@ def aggregate(
         None,
         help="Column of the transcript dataframe representing the gene names. If not provided, it will not compute transcript count",
     ),
-    average_intensities: bool = typer.Option(
-        False, help="Whether to average the channel intensities inside each cell"
-    ),
+    average_intensities: bool = typer.Option(False, help="Whether to average the channel intensities inside each cell"),
     expand_radius_ratio: float = typer.Option(
         default=0,
         help="Cells polygons will be expanded by `expand_radius_ratio * mean_radius` for channels averaging **only**. This help better aggregate boundary stainings",
     ),
-    min_transcripts: int = typer.Option(
-        0, help="Cells with less transcript than this integer will be filtered"
-    ),
+    min_transcripts: int = typer.Option(0, help="Cells with less transcript than this integer will be filtered"),
     min_intensity_ratio: float = typer.Option(
         0,
         help="Cells whose mean channel intensity is less than `min_intensity_ratio * quantile_90` will be filtered",

sopa/cli/check.py

@@ -16,9 +16,7 @@ def _open_config(path: str) -> dict:
         with open(path, "r") as f:
             return yaml.safe_load(f)
     except:
-        log.warn(
-            f"Config file '{path}' could't be read. Make sure that the file exist and that it is a YAML file"
-        )
+        log.warn(f"Config file '{path}' could't be read. Make sure that the file exist and that it is a YAML file")
         return
 
 
@@ -99,8 +97,7 @@ def _check_dict(config: dict, d: dict, log, prefix: str = "config"):
                     break
             else:
                 display = "\n  - ".join(
-                    element if isinstance(element, str) else " AND ".join(element)
-                    for element in values
+                    element if isinstance(element, str) else " AND ".join(element) for element in values
                 )
                 log.warn(f"One of these element must be in {prefix}['{key}']:\n  - {display}")
 

sopa/cli/explorer.py

@@ -6,9 +6,7 @@
 
 app_explorer = typer.Typer()
 
-PIXELSIZE_DEPRECATED = (
-    "`pixelsize` is deprecated and will be removed in future versions. Use `pixel_size` instead."
-)
+PIXELSIZE_DEPRECATED = "`pixelsize` is deprecated and will be removed in future versions. Use `pixel_size` instead."
 
 
 @app_explorer.command()
@@ -18,9 +16,7 @@ def write(
         None,
         help="Path to a directory where Xenium Explorer's outputs will be saved. By default, writes to the same path as `sdata_path` but with the `.explorer` suffix",
     ),
-    gene_column: str = typer.Option(
-        None, help="Column name of the points dataframe containing the gene names"
-    ),
+    gene_column: str = typer.Option(None, help="Column name of the points dataframe containing the gene names"),
     shapes_key: str = typer.Option(
         None,
         help="Sdata key for the boundaries. By default, uses the baysor boundaires, else the cellpose boundaries",
@@ -134,6 +130,4 @@ def add_aligned(
     sdata = spatialdata.read_zarr(sdata_path)
     image = io.ome_tif(image_path, as_image=True)
 
-    align(
-        sdata, image, transformation_matrix_path, overwrite=overwrite, image_key=original_image_key
-    )
+    align(sdata, image, transformation_matrix_path, overwrite=overwrite, image_key=original_image_key)