add spot technologies support

Author: quentinblampey

CHANGELOG.md

@@ -1,3 +1,9 @@
+## [0.1.x] - tbd
+
+### Added
+- Support spot-based technologies like Visium
+- Faster table conversion (using native CSR matrix arguments)
+
 ## [0.1.4] - 2024-01-11
 
 ### Fix

spatialdata_xenium_explorer/_constants.py

@@ -24,6 +24,11 @@ class Versions:
     CELL_CATEGORIES = [1, 0]
 
 
+class ShapesConstants:
+    RADIUS = "radius"
+    DEFAULT_POINT_RADIUS = 100
+
+
 def cell_categories_attrs() -> dict:
     return {
         "major_version": Versions.CELL_CATEGORIES[0],

spatialdata_xenium_explorer/converter.py

@@ -8,14 +8,14 @@
 from spatialdata import SpatialData
 
 from . import (
+    utils,
     write_cell_categories,
     write_gene_counts,
     write_image,
     write_polygons,
     write_transcripts,
 )
 from ._constants import FileNames, experiment_dict
-from .utils import explorer_file_path, get_element, get_spatial_image, to_intrinsic
 
 log = logging.getLogger(__name__)
 
@@ -28,6 +28,7 @@ def write(
     points_key: str | None = None,
     gene_column: str | None = None,
     pixelsize: float = 0.2125,
+    spot: bool = False,
     layer: str | None = None,
     polygon_max_vertices: int = 13,
     lazy: bool = True,
@@ -64,6 +65,7 @@ def write(
         points_key: Name of the transcripts (key of `sdata.points`). This argument doesn't need to be provided if there is only one points key.
         gene_column: Column name of the points dataframe containing the gene names.
         pixelsize: Number of microns in a pixel. Invalid value can lead to inconsistent scales in the Explorer.
+        spot: Whether the technology is based on spots
         layer: Layer of `sdata.table` where the gene counts are saved. If `None`, uses `sdata.table.X`.
         polygon_max_vertices: Maximum number of vertices for the cell polygons. A higher value will display smoother cells.
         lazy: If `True`, will not load the full images in memory (except if the image memory is below `ram_threshold_gb`).
@@ -73,7 +75,7 @@ def write(
     path: Path = Path(path)
     _check_explorer_directory(path)
 
-    image_key, image = get_spatial_image(sdata, image_key, return_key=True)
+    image_key, image = utils.get_spatial_image(sdata, image_key, return_key=True)
 
     ### Saving cell categories and gene counts
     if sdata.table is not None:
@@ -94,22 +96,27 @@ def write(
             write_cell_categories(path, adata)
 
     ### Saving cell boundaries
-    shapes_key, geo_df = get_element(sdata, "shapes", shapes_key, return_key=True)
+    shapes_key, geo_df = utils.get_element(sdata, "shapes", shapes_key, return_key=True)
 
     if _should_save(mode, "b") and geo_df is not None:
-        geo_df = to_intrinsic(sdata, geo_df, image_key)
+        geo_df = utils.to_intrinsic(sdata, geo_df, image_key)
 
         if sdata.table is not None:
             geo_df = geo_df.loc[adata.obs[adata.uns["spatialdata_attrs"]["instance_key"]]]
 
+        geo_df = utils._standardize_shapes(geo_df)
+
         write_polygons(path, geo_df.geometry, polygon_max_vertices, pixelsize=pixelsize)
 
     ### Saving transcripts
-    df = get_element(sdata, "points", points_key)
+    if spot and sdata.table is not None:
+        df, gene_column = utils._spot_transcripts_origin(adata)
+    else:
+        df = utils.get_element(sdata, "points", points_key)
+        df = utils.to_intrinsic(sdata, df, image_key)
 
     if _should_save(mode, "t") and df is not None:
         if gene_column is not None:
-            df = to_intrinsic(sdata, df, image_key)
             write_transcripts(path, df, gene_column, pixelsize=pixelsize)
         else:
             log.warn("The argument 'gene_column' has to be provided to save the transcripts")
@@ -172,7 +179,7 @@ def write_metadata(
         is_dir: If `False`, then `path` is a path to a single file, not to the Xenium Explorer directory.
         pixelsize: Number of microns in a pixel. Invalid value can lead to inconsistent scales in the Explorer.
     """
-    path = explorer_file_path(path, FileNames.METADATA, is_dir)
+    path = utils.explorer_file_path(path, FileNames.METADATA, is_dir)
 
     with open(path, "w") as f:
         metadata = experiment_dict(image_key, shapes_key, n_obs, pixelsize)

spatialdata_xenium_explorer/core/points.py

@@ -112,7 +112,7 @@ def write_transcripts(
             GRIDS_ATTRS["grid_number_objects"].append([])
             GRIDS_ATTRS["grid_keys"].append([])
 
-            n_tiles_x, n_tiles_y = ceil(xmax / tile_size), ceil(ymax / tile_size)
+            n_tiles_x, n_tiles_y = max(1, ceil(xmax / tile_size)), max(1, ceil(ymax / tile_size))
 
             for tx in range(n_tiles_x):
                 for ty in range(n_tiles_y):
@@ -179,7 +179,7 @@ def write_transcripts(
                         chunks=chunks,
                     )
 
-            if n_tiles_x * n_tiles_y == 1:
+            if n_tiles_x * n_tiles_y == 1 and level > 0:
                 GRIDS_ATTRS["number_levels"] = level + 1
                 break
 

spatialdata_xenium_explorer/utils.py

@@ -4,15 +4,22 @@
 from pathlib import Path
 
 import dask.array as da
+import dask.dataframe as dd
 import dask_image.ndinterp
+import geopandas as gpd
 import numpy as np
+import pandas as pd
 import xarray as xr
+from anndata import AnnData
 from multiscale_spatial_image import MultiscaleSpatialImage
+from shapely.geometry import MultiPolygon, Point, Polygon
 from spatial_image import SpatialImage
 from spatialdata import SpatialData
 from spatialdata.models import SpatialElement
 from spatialdata.transformations import Identity, get_transformation, set_transformation
 
+from ._constants import ShapesConstants
+
 log = logging.getLogger(__name__)
 
 
@@ -213,3 +220,41 @@ def scale_dtype(arr: np.ndarray, dtype: np.dtype) -> np.ndarray:
 
     factor = np.iinfo(dtype).max / np.iinfo(arr.dtype).max
     return (arr * factor).astype(dtype)
+
+
+def _standardize_shapes(geo_df: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
+    if geo_df.geometry.map(lambda geom: isinstance(geom, Polygon)).all():
+        return geo_df
+
+    if geo_df.geometry.map(lambda geom: isinstance(geom, Point)).all():
+        if not ShapesConstants.RADIUS in geo_df:
+            log.warn(
+                f"GeoDataFrame contains only Point objects, but no column '{ShapesConstants.RADIUS}' was found. Using default {ShapesConstants.RADIUS}={ShapesConstants.DEFAULT_POINT_RADIUS}"
+            )
+            geo_df[ShapesConstants.RADIUS] = ShapesConstants.DEFAULT_POINT_RADIUS
+
+        geo_df.geometry = geo_df.apply(lambda row: row.geometry.buffer(row.radius), axis=1)
+        return geo_df
+
+    if geo_df.geometry.map(lambda geom: isinstance(geom, MultiPolygon)).all():
+        log.warn(
+            "GeoDataFrame contains only MultiPolygon objects. For each MultiPolygon, only the Polygon with the largest area will be shown"
+        )
+
+        geo_df.geometry = geo_df.geometry.map(
+            lambda multi_polygon: max(multi_polygon.geoms, key=lambda geom: geom.area)
+        )
+        return geo_df
+
+    raise ValueError(
+        "The provided shapes contain unsupported types, or contain a mix of multiple types. Supported types: Polygon, Point, MultiPolygon."
+    )
+
+
+def _spot_transcripts_origin(adata: AnnData) -> tuple[dd.DataFrame, str]:
+    gene_column = "gene"
+    df = pd.DataFrame(
+        {"x": [0] * adata.n_vars, "y": [0] * adata.n_vars, gene_column: adata.var_names}
+    )
+    df = dd.from_pandas(df, chunksize=10_000)
+    return df, gene_column