first version working

Author: quentinblampey

README.md

@@ -14,9 +14,17 @@ import spatialdata
 import spatialdata_xenium_explorer
 
 sdata = spatialdata.read_zarr("...")
-image_key = "..." # The name of the MultiscaleSpatialImage to be exported
 
-spatialdata_xenium_explorer.write("/path/to/directory", sdata, image_key)
+spatialdata_xenium_explorer.write("/path/to/directory", sdata, image_key, shapes_key, points_key, gene_column)
 ```
 
-This will create up to 6 files, among which a file called `experiment.xenium`. Double-click on this file to open it on the [Xenium Explorer](https://www.10xgenomics.com/support/software/xenium-explorer/downloads) (make sure you have the latest version of the software).
+For more details about the arguments, see the [function docstrings](https://github.com/quentinblampey/spatialdata_xenium_explorer/blob/master/spatialdata_xenium_explorer/converter.py#L29).
+
+This will create up to 6 files, among which a file called `experiment.xenium`. Double-click on this file to open it on the [Xenium Explorer](https://www.10xgenomics.com/support/software/xenium-explorer/downloads) (make sure you have the latest version of the software).
+
+## Future improvements
+
+- Support all types of images (not just `MultiscaleSpatialImage`)
+- Better user experience (less arguments)
+- Write `.tif` image without loading the spatial image in memory
+- Write transcripts without computing the whole coordinates

spatialdata_xenium_explorer/constants.py spatialdata_xenium_explorer/_constants.pyspatialdata_xenium_explorer/constants.py renamed to spatialdata_xenium_explorer/_constants.py

@@ -10,6 +10,7 @@ class FileNames:
 class ExplorerConstants:
     GRID_SIZE = 250
     QUALITY_SCORE = 40
+    MICRONS_TO_PIXELS = 4.705882
 
 
 class Versions:

spatialdata_xenium_explorer/converter.py

@@ -11,10 +11,10 @@
     write_polygons,
     write_transcripts,
 )
-from .constants import FileNames, experiment_dict
+from ._constants import FileNames, experiment_dict
 
 
-def _order_instances(sdata: SpatialData, shapes_key: str):
+def _reorder_instances(sdata: SpatialData, shapes_key: str):
     adata = sdata.table
 
     instance_key = adata.uns["spatialdata_attrs"]["instance_key"]
@@ -34,6 +34,7 @@ def write(
     points_key: str,
     gene_column: str,
     layer: str | None = None,
+    polygon_max_vertices: int = 13,
 ) -> None:
     """
     Transform a SpatialData object into inputs for the Xenium Explorer.
@@ -47,6 +48,7 @@ def write(
         points_key: Name of the transcripts (key of `sdata.points`).
         gene_column: Column name of the points dataframe containing the gene names.
         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.
     """
     path: Path = Path(path)
     assert (
@@ -55,7 +57,7 @@ def write(
 
     path.mkdir(parents=True, exist_ok=True)
 
-    adata = _order_instances(sdata, shapes_key)
+    adata = _reorder_instances(sdata, shapes_key)
 
     EXPERIMENT = experiment_dict(image_key, shapes_key, adata.n_obs)
     with open(path / FileNames.METADATA, "w") as f:
@@ -68,9 +70,9 @@ def write(
     set_transformation(sdata.images[image_key], Identity(), pixels_cs)
 
     gdf = sdata.transform_element_to_coordinate_system(sdata.shapes[shapes_key], pixels_cs)
-    write_polygons(path / FileNames.SHAPES, gdf.geometry)
-
-    write_multiscale(path / FileNames.IMAGE, sdata.images[image_key])
+    write_polygons(path / FileNames.SHAPES, gdf.geometry, polygon_max_vertices)
 
     df = sdata.transform_element_to_coordinate_system(sdata.points[points_key], pixels_cs)
     write_transcripts(path / FileNames.POINTS, df, gene_column)
+
+    write_multiscale(path / FileNames.IMAGE, sdata.images[image_key])

spatialdata_xenium_explorer/images.py

@@ -2,10 +2,11 @@
 import tifffile as tf
 from multiscale_spatial_image import MultiscaleSpatialImage
 
-from .constants import image_metadata, image_options
+from ._constants import image_metadata, image_options
 
 
 def _astype_uint8(arr: np.ndarray) -> np.ndarray:
+    print(f"   Image of shape {arr.shape}")
     assert np.issubdtype(
         arr.dtype, np.integer
     ), f"The image dtype has to be an integer dtype. Found {arr.dtype}"
@@ -22,6 +23,7 @@ def write_multiscale(
     multiscale: MultiscaleSpatialImage,
     pixelsize: float = 0.2125,
 ):
+    print("Writing multiscale image")
     scale_names = list(multiscale.children)
     channel_names = list(multiscale[scale_names[0]].c.values)
 

spatialdata_xenium_explorer/points.py

@@ -5,7 +5,7 @@
 import numpy as np
 import zarr
 
-from .constants import ExplorerConstants
+from ._constants import ExplorerConstants
 
 
 def subsample_indices(n_samples, factor: int = 4):
@@ -20,33 +20,28 @@ def write_transcripts(
     max_levels: int = 15,
 ):
     # TODO: make everything using dask instead of pandas
+    print(f"Writing {len(df)} transcripts")
     df = df.compute()
 
     num_transcripts = len(df)
     df[gene] = df[gene].astype("category")
 
     location = df[["x", "y"]]
+    location /= ExplorerConstants.MICRONS_TO_PIXELS
     location = np.concatenate([location, np.zeros((num_transcripts, 1))], axis=1)
 
     xmax, ymax = location[:, :2].max(axis=0)
 
-    assert location[:, 0].min() >= 0
-    assert location[:, 1].min() >= 0
-
     gene_names = list(df[gene].cat.categories)
     num_genes = len(gene_names)
 
     codeword_gene_mapping = list(range(num_genes))
 
     valid = np.ones((num_transcripts, 1))
     uuid = np.stack([np.arange(num_transcripts), np.full(num_transcripts, 65535)], axis=1)
-    transcript_id = np.stack(
-        [np.arange(num_transcripts), np.full(num_transcripts, 65535)], axis=1
-    )
+    transcript_id = np.stack([np.arange(num_transcripts), np.full(num_transcripts, 65535)], axis=1)
     gene_identity = df[gene].cat.codes.values[:, None]
-    codeword_identity = np.stack(
-        [gene_identity[:, 0], np.full(num_transcripts, 65535)], axis=1
-    )
+    codeword_identity = np.stack([gene_identity[:, 0], np.full(num_transcripts, 65535)], axis=1)
     status = np.zeros((num_transcripts, 1))
     quality_score = np.full((num_transcripts, 1), ExplorerConstants.QUALITY_SCORE)
 
@@ -55,9 +50,7 @@ def write_transcripts(
         "codeword_gene_mapping": codeword_gene_mapping,
         "codeword_gene_names": gene_names,
         "gene_names": gene_names,
-        "gene_index_map": {
-            name: index for name, index in zip(gene_names, codeword_gene_mapping)
-        },
+        "gene_index_map": {name: index for name, index in zip(gene_names, codeword_gene_mapping)},
         "number_genes": num_genes,
         "spatial_units": "micron",
         "coordinate_space": "refined-final_global_micron",
@@ -85,12 +78,11 @@ def write_transcripts(
         grids = g.create_group("grids")
 
         for level in range(max_levels):
+            print(f"   Level {level}: {len(location)} transcripts")
             level_group = grids.create_group(level)
 
             tile_size = ExplorerConstants.GRID_SIZE * 2**level
 
-            print(f"Level {level}: {len(location)} transcripts")
-
             indices = np.floor(location[:, :2] / tile_size).clip(0).astype(int)
             tiles_str_indices = np.array([f"{tx},{ty}" for (tx, ty) in indices])
 

spatialdata_xenium_explorer/shapes.py

@@ -6,7 +6,7 @@
 import zarr
 from shapely.geometry import Polygon
 
-from .constants import cell_summary_attrs, group_attrs
+from ._constants import ExplorerConstants, cell_summary_attrs, group_attrs
 
 
 def pad_polygon(polygon: Polygon, max_vertices: int, tolerance: float = 1) -> np.ndarray:
@@ -19,18 +19,17 @@ def pad_polygon(polygon: Polygon, max_vertices: int, tolerance: float = 1) -> np
         return coords.flatten()
 
     if n_vertices < max_vertices:
-        return np.pad(
-            coords, ((0, max_vertices - n_vertices), (0, 0)), mode="edge"
-        ).flatten()
+        return np.pad(coords, ((0, max_vertices - n_vertices), (0, 0)), mode="edge").flatten()
 
     # TODO: improve it: how to choose the right tolerance?
     polygon = polygon.simplify(tolerance=tolerance)
     return pad_polygon(polygon, max_vertices, tolerance + 1)
 
 
 def write_polygons(path: Path, polygons: Iterable[Polygon], max_vertices: int) -> None:
+    print(f"Writing {len(polygons)} cell polygons")
     coordinates = np.stack([pad_polygon(p, max_vertices) for p in polygons])
-    coordinates /= 4.705882
+    coordinates /= ExplorerConstants.MICRONS_TO_PIXELS
 
     num_cells = len(coordinates)
     cells_fourth = ceil(num_cells / 4)

spatialdata_xenium_explorer/table.py

@@ -3,10 +3,11 @@
 from anndata import AnnData
 from scipy.sparse import csr_matrix
 
-from .constants import cell_categories_attrs
+from ._constants import cell_categories_attrs
 
 
 def write_gene_counts(path: str, adata: AnnData, layer: str | None) -> None:
+    print(f"Writing table of {adata.n_vars} genes")
     counts = adata.X if layer is None else adata.layers[layer]
     counts = csr_matrix(counts)
 
@@ -74,7 +75,7 @@ def write_cell_categories(path: str, adata: AnnData) -> None:
     # TODO: consider also columns that can be transformed to a categorical column?
     cat_columns = [name for name, cat in adata.obs.dtypes.items() if cat == "category"]
 
-    print(f"Saving {len(cat_columns)} cell categories: {', '.join(cat_columns)}")
+    print(f"Writing {len(cat_columns)} cell categories: {', '.join(cat_columns)}")
 
     ATTRS = cell_categories_attrs()
     ATTRS["number_groupings"] = len(cat_columns)