Reading segmented data with Visium_HD (spaceranger>4)

[ddiez]

The newest version (>4) of spaceranger supports nuclei/cell segmentation of VisiumHD data. I wonder if there are official on going plans to give support to this new feature?

I am trying to get this data into spatialdata object but having some difficulties. This is how far I have come.

I can read the files with the nuclei and cell level segmentation:

sdata_hm["HM_nucleus_segmentation"] = io.geojson("data/HM/outs/segmented_outputs/nucleus_segmentations.geojson", coordinate_system="HM")
sdata_hm["HM_cell_segmentation"] = io.geojson("data/HM/outs/segmented_outputs/cell_segmentations.geojson", coordinate_system="HM")
sdata_hm

SpatialData object
├── Images
│     ├── 'HM_hires_image': DataArray[cyx] (3, 5057, 6000)
│     └── 'HM_lowres_image': DataArray[cyx] (3, 506, 600)
├── Shapes
│     ├── 'HM_cell_segmentation': GeoDataFrame shape: (69604, 1) (2D shapes)
│     ├── 'HM_nucleus_segmentation': GeoDataFrame shape: (69604, 1) (2D shapes)
│     └── 'HM_square_016um': GeoDataFrame shape: (48206, 1) (2D shapes)
└── Tables
      └── 'square_016um': AnnData (48206, 18085)
with coordinate systems:
    ▸ 'HM', with elements:
        HM_hires_image (Images), HM_lowres_image (Images), HM_cell_segmentation (Shapes), HM_nucleus_segmentation (Shapes), HM_square_016um (Shapes)
    ▸ 'HM_downscaled_hires', with elements:
        HM_hires_image (Images), HM_square_016um (Shapes)
    ▸ 'HM_downscaled_lowres', with elements:
        HM_lowres_image (Images), HM_square_016um (Shapes)

I can read and assign the AnnData table.

adata = sc.read_10x_h5("data/HM/outs/segmented_outputs/filtered_feature_cell_matrix.h5", gex_only=False)
adata.var_names_make_unique()
sdata_hm["cell_bins"] = adata
sdata_hm

SpatialData object, with associated Zarr store: [/media/diez/Biodev/projects/yamanaka-kidney/output/sdata_hm.zarr](https://vscode-remote+ssh-002dremote-002btower.vscode-resource.vscode-cdn.net/media/diez/Biodev/projects/yamanaka-kidney/output/sdata_hm.zarr)
├── Images
│     ├── 'HM_hires_image': DataArray[cyx] (3, 5057, 6000)
│     └── 'HM_lowres_image': DataArray[cyx] (3, 506, 600)
├── Shapes
│     ├── 'HM_cell_segmentation': GeoDataFrame shape: (69604, 1) (2D shapes)
│     ├── 'HM_nucleus_segmentation': GeoDataFrame shape: (69604, 1) (2D shapes)
│     └── 'HM_square_016um': GeoDataFrame shape: (48206, 1) (2D shapes)
└── Tables
      ├── 'cell_bins': AnnData (69601, 18085)
      └── 'square_016um': AnnData (48206, 18085)
with coordinate systems:
    ▸ 'HM', with elements:
        HM_hires_image (Images), HM_lowres_image (Images), HM_cell_segmentation (Shapes), HM_nucleus_segmentation (Shapes), HM_square_016um (Shapes)
    ▸ 'HM_downscaled_hires', with elements:
        HM_hires_image (Images), HM_square_016um (Shapes)
    ▸ 'HM_downscaled_lowres', with elements:
        HM_lowres_image (Images), HM_square_016um (Shapes)

With this set up I can visualize the shapes for the nuclei and cells. I can also select a ROI and visualize it and everything looks OK. But I cannot plot the gene expression. This is likely because I did not pass the AnnData through TableModel.parse. However, in order to do so I need to add additional data, one comming from VisiumHDKeys.TISSUE_POSITIONS_FILE that is tissue_positions.parquet and it is located in the spatial folder for the binned_outputs. However this file is missing in the segmented_outputs. I need to think about it because I still have gaps in my knowledge of the spatialdata framework but perhaps this makes sense since the cell shapes are irregular.

If no one is working on this I would like to help getting this supported in spatialdata although I would need some hints on what it is missing at the moment. Thanks.

[ddiez]

For reference, this is the code I am using now that collects most of the info and highlights the parts that are missing from the current approach with bin data. It returns a table that can be assigned to a SpatialData object.

def read_visium_hd_segmented(path):
  from spatialdata_io._constants._constants import VisiumHDKeys
  from spatialdata.models import Image2DModel, ShapesModel, TableModel
  from spatialdata.transformations import Affine, Identity, Scale, set_transformation
  import json
  from pathlib import Path
  import scanpy as sc
  import numpy as np
  

  path = Path(path)

  bin_size = "segmented_outputs"

  #counts_file = VisiumHDKeys.FILTERED_COUNTS_FILE # We may need a new registered name for this.
  counts_file = "filtered_feature_cell_matrix.h5"
  adata = sc.read_10x_h5(path / bin_size / counts_file, gex_only=False)
  adata.var_names_make_unique()

  path_bin = path / bin_size
  path_bin_spatial = path_bin / VisiumHDKeys.SPATIAL

  with open(path_bin_spatial / VisiumHDKeys.SCALEFACTORS_FILE) as file:
    scalefactors = json.load(file)

  shapes_name = "cells"

  #tissue_positions_file = path_bin_spatial / VisiumHDKeys.TISSUE_POSITIONS_FILE
  #coords = pd.read_parquet(tissue_positions_file) # This file does not exist in the segmented outputs.
  #coords.set_index(VisiumHDKeys.BARCODE, inplace=True, drop=True)
  #coords_filtered = coords.loc[adata.obs.index]
  #adata.obsm["spatial"] = data.obs[[VisiumHDKeys.LOCATIONS_X, VisiumHDKeys.LOCATIONS_Y]].values # compatibility with squidpy.
  # adata.obs.drop(
  #     columns=[
  #         VisiumHDKeys.LOCATIONS_X,
  #         VisiumHDKeys.LOCATIONS_Y,
  #     ],
  #     inplace=True,
  # )
  adata.obs[VisiumHDKeys.INSTANCE_KEY] = np.arange(len(adata))

  transform_original = Identity()
  transform_lowres = Scale(
      np.array(
          [
              scalefactors[VisiumHDKeys.SCALEFACTORS_LOWRES],
              scalefactors[VisiumHDKeys.SCALEFACTORS_LOWRES],
          ]
      ),
      axes=("x", "y"),
  )
  transform_hires = Scale(
      np.array(
          [
              scalefactors[VisiumHDKeys.SCALEFACTORS_HIRES],
              scalefactors[VisiumHDKeys.SCALEFACTORS_HIRES],
          ]
      ),
      axes=("x", "y"),
  )
  
  adata.obs[VisiumHDKeys.REGION_KEY] = shapes_name
  adata.obs[VisiumHDKeys.REGION_KEY] = adata.obs[VisiumHDKeys.REGION_KEY].astype("category")
  
  return TableModel.parse(
    adata=adata,
    region=shapes_name,
    region_key=str(VisiumHDKeys.REGION_KEY),
    instance_key=str(VisiumHDKeys.INSTANCE_KEY)
  )

[ddiez]

Thanks to some user that posted a comment here and then, it seems, later deleted it, I could fix the code above. I was able to read it because I got the email notification. The main points are:

1. Use raw_feature_cell_matrix.h5 so that the matrix has the same shape as the boundaries.
2. Use shapes_name = "cell_segmentation" as the name of the shapes and the regions in the AnnData have to match.

With this, it is possible to have the cell masks colored with gene expression values. I was testing myself another code not yet posted, based on the original spatialdata_io.visium_hd code, with the goal to add directly support. In that code I had already started using the raw feature matrix because of the shape issue but could not find the issue with the name of the shape segmentation. I need to polish the code further before considering a PR.

I would like to thank the annonymous user that helped me find this issue.

[chrissikath]

It was me, I took it back because I wanted to be able to use the scanpy plotting functions too.
It should now work for both SpatialData and Scanpy plotting. I'm not sure if this is the right way to do it, though.

def read_visium_hd_segmented(outs_path: str, sample_id: str) -> TableModel:
    outs_path = Path(outs_path)
    bin_size = "segmented_outputs"
    path_bin = outs_path / bin_size
    path_bin_spatial = path_bin / VisiumHDKeys.SPATIAL

    # Load gene expression
    counts_file = "raw_feature_cell_matrix.h5"
    adata = sc.read_10x_h5(path_bin / counts_file, gex_only=False)
    adata.var_names_make_unique()

    # Load scalefactors and images
    with open(path_bin_spatial / VisiumHDKeys.SCALEFACTORS_FILE) as f:
        scalefactors = json.load(f)

    hires_img = np.array(Image.open(path_bin_spatial / "tissue_hires_image.png"))
    lowres_img = np.array(Image.open(path_bin_spatial / "tissue_lowres_image.png"))

    library_id = "visium_hd_segmentation"
    adata.uns["spatial"] = {
        library_id: {
            "images": {
                "hires": hires_img,
                "lowres": lowres_img,
            },
            "scalefactors": scalefactors,
            "metadata": {
                "source_image_path": "tissue_hires_image.png"
            }
        }
    }

    # Assign region info
    shapes_name = "cell_segmentation"
    adata.obs[VisiumHDKeys.INSTANCE_KEY] = np.arange(len(adata))
    adata.obs[VisiumHDKeys.REGION_KEY] = shapes_name
    adata.obs[VisiumHDKeys.REGION_KEY] = adata.obs[VisiumHDKeys.REGION_KEY].astype("category")

    # Load cell shapes
    shapes = spatialdata_io.geojson(outs_path / "segmented_outputs" / "cell_segmentations.geojson", coordinate_system=sample_id)

    # Match shape order to adata
    centroids = shapes.geometry.centroid
    adata.obsm["spatial"] = np.vstack([centroids.x.values, centroids.y.values]).T

    # Estimate spot diameter
    areas = shapes.geometry.area
    mean_area = np.mean(areas)
    diameter_pixels = 2 * np.sqrt(mean_area / np.pi)
    adata.uns["spatial"][library_id]["scalefactors"]["spot_diameter_fullres"] = diameter_pixels

    return TableModel.parse(
        adata=adata,
        region=shapes_name,
        region_key=str(VisiumHDKeys.REGION_KEY),
        instance_key=str(VisiumHDKeys.INSTANCE_KEY),
    )

[ddiez]

@chrissikath finally had time to get into this again. Thanks for sharing the updated code. It works nicely to enable plotting with scanpy out of the box.