# spatialproteomics > spatialproteomics is a Python toolbox for analyzing highly multiplexed fluorescence images (e.g. CODEX, MIBI, IMC): image processing, cell segmentation, marker quantification, cell type prediction, neighborhood analysis and plotting. It stores each sample as an xarray.Dataset whose layers (image, segmentation, expression matrix, labels) stay synchronized when subsetting. --- # Core concepts and cheatsheet This page summarizes how spatialproteomics is organized and how its API is meant to be used. It is a compact reference for both humans and AI coding assistants; the tutorials cover each topic in depth. ## The data model A spatialproteomics object is a plain `xarray.Dataset`. All data belonging to one sample lives in this single object as named data variables ("layers") that share dimensions: | Layer (`sp.Layers.*`) | Dimensions | Content | |----|----|----| | `_image` | `(channels, y, x)` | Multiplexed image; one channel per marker. | | `_segmentation` | `(y, x)` | Integer cell mask; `0` is background, `n` is cell `n`. | | `_obs` | `(cells, features)` | Per-cell observations (centroids, area, cell type, neighborhood). | | `_intensity` | `(cells, channels)` | Expression matrix produced by `pp.add_quantification`. | | `_la_properties` | `(labels, la_props)` | Cell type names and colors. | | `_neighborhoods` | `(cells, labels)` | Cell type composition of each cell's neighborhood. | | `_nh_properties` | `(neighborhoods, nh_props)` | Neighborhood names and colors. | You can add layers of your own with `ds.pp.add_layer(array, key_added="_name")`. A common one is `_raw_image`: an untouched copy of `_image` kept before thresholding, because `pp.threshold` overwrites `_image` (see the workflow and examples below). Because these layers share coordinates, subsetting is synchronized: cropping the image also crops the segmentation and drops cells outside the crop from `_obs` and `_intensity`; selecting a cell type removes the other cells from the mask and all tables. Use the constants in `sp.Layers`, `sp.Dims`, `sp.Features` and `sp.Props` instead of hard-coding strings, e.g. `sp.Features.LABELS == "_labels"` (the cell type column in `_obs`), `sp.Features.X == "centroid-1"`. ## Accessors Importing spatialproteomics registers five accessors on every `xarray.Dataset`: | Name | Class | Responsibility | |----|----|----| | `pp` | PreprocessingAccessor | Subsetting, image processing, segmentation handling, quantification, tables. | | `tl` | ToolAccessor | Wrappers around external tools (cellpose, stardist, instanseg, astir, spotiflow) and export. | | `la` | LabelAccessor | Cell type labels: prediction, naming, colors, subtype hierarchies. | | `nh` | NeighborhoodAccessor | Cellular neighborhoods and graph features. | | `pl` | PlotAccessor | Plotting of images, masks, labels, neighborhoods and observations. | Every accessor method returns a **new** dataset; nothing is modified in place. Always reassign (`ds = ds.pp.threshold(...)`) or chain calls: ``` python ds = ds.pp.threshold(quantile=0.5).pp.add_quantification().pp.transform_expression_matrix() ``` ## Indexing ``` python ds.pp["DAPI"] # one channel ds.pp[["CD4", "CD8"]] # several channels ds.pp[100:200, 300:400] # spatial crop (x slice, y slice); cells outside are removed ds.pp[["CD4", "CD8"], 100:200, 300:400] ds.la["B cell"] # cells of one type (by name or integer label id) ds.la[["B cell", "T cell"]] ds.la.deselect(["B cell"]) # all cells except these types ds.nh["Neighborhood 1"] # cells of one neighborhood ``` ## Canonical workflow (xarray backend) ``` python import spatialproteomics as sp import xarray as xr from skimage.io import imread # 1. create the object from a (channels, y, x) array image = imread("input.tiff") ds = sp.load_image_data(image, channel_coords=["DAPI", "PAX5", "CD3", "CD4", "CD8"]) # 2. segmentation (requires `pip install cellpose`); adds `_segmentation` and `_obs` ds = ds.tl.cellpose(channel="DAPI") ds = ds.pp.add_observations("area") ds = ds.pp.filter_by_obs("area", func=lambda x: (x > 20) & (x < 150)) ds = ds.pp.grow_cells(iterations=2) # 3. image processing: per-channel thresholding (overwrites `_image`) # `quantile=` is relative to each channel's intensity distribution, `intensity=` is an absolute value; # pass one value per channel, or one value plus `channels=[...]` to threshold only some channels. # Thresholds require human inspection: plot the raw image, the thresholded image and the resulting # cell types (step 5) side by side, as in the "Cell type QC plot" example below. ds = ds.pp.add_layer(ds["_image"].values.copy(), key_added="_raw_image") # keep the original for comparison ds = ds.pp.threshold(quantile=[0.2, 0.6, 0.6, 0.6, 0.6]) # 4. quantification; adds `_intensity` ds = ds.pp.add_quantification(func="intensity_mean").pp.transform_expression_matrix(method="arcsinh") # 5. cell types; maps marker -> cell type ds = ds.la.predict_cell_types_argmax({"PAX5": "B cell", "CD3": "T cell"}) ds = ds.la.set_label_colors(["B cell", "T cell"], ["red", "green"]) # 6. binarization and subtypes; maps marker -> subtype # quantify the fraction of positive (non-zero) pixels per cell on the thresholded image, then threshold # that fraction; this adds a binary `{marker}_binarized` column to `_obs` (e.g. `CD4_binarized`). # See the "Binarization QC plot" example below for checking the result visually. ds = ds.pp.add_quantification(func=sp.percentage_positive, key_added="_percentage_positive") ds = ds.la.threshold_labels({"CD4": 0.5, "CD8": 0.5}, layer_key="_percentage_positive") ds = ds.la.predict_cell_subtypes({"T cell": {"CD4": "T_h", "CD8": "T_tox"}}) ds = ds.la.set_label_colors(["B cell", "T cell", "T_h", "T_tox"], ["red", "green", "yellow", "blue"]) # 7. plotting (all plotting functions accept `ax=`) # plotting an image ds.pp[["CD4", "CD8"]].pl.colorize(["red", "green"]).pl.show() # plotting the raw image instead of the thresholded one ds.pp[["CD4", "CD8"]].pl.colorize(["red", "green"], layer_key="_raw_image").pl.show() # plotting a segmentation mask and the labels ds.pl.show(render_image=False, render_segmentation=True, render_labels=True) # 8. tables and export obs = ds.pp.get_layer_as_df("_obs") # cell types decoded to names expr = ds.pp.get_layer_as_df("_intensity") ds.drop_encoding().to_zarr("sample.zarr") # reload with xr.open_zarr("sample.zarr") adata = ds.tl.convert_to_anndata() # requires anndata sdata = ds.tl.convert_to_spatialdata() # requires spatialdata ``` Other frequently used operations: - Segmentation from elsewhere: `ds.pp.add_segmentation(mask)`; alternatives to cellpose: `ds.tl.stardist()` and `ds.tl.instanseg()` (segments all selected channels jointly) - Labels from elsewhere: `ds.la.add_labels_from_dataframe(df, cell_col="cell", label_col="label")`. - Marker binarization and subtypes: `ds.la.threshold_labels({...})` then `ds.la.predict_cell_subtypes({...})`; switch levels with `ds.la.set_label_level(...)`. - Neighborhoods across samples: `sp.ImageContainer({"s1": ds1, "s2": ds2}).compute_neighborhoods(k=5)` computes compositions, clusters them with k-means and returns a dict of datasets with neighborhoods assigned (`ds.nh["Neighborhood 0"]`). - Composable plotting: `ds.pl.colorize([...]).pl.render_segmentation().pl.render_labels().pl.imshow()`; scatter plots with `ds.pl.scatter_labels()` / `ds.pl.scatter(feature)`. - Large images: `ds.pp.downsample(rate=8)` before plotting, `xr.open_zarr` for lazy loading. ## The spatialdata backend If you start from a `spatialdata.SpatialData` object (`pip install "spatialproteomics[spatialdata]"`), use the module-level functions instead of the accessors. They follow scverse conventions: they take the object as first argument and **modify it in place** unless `copy=True`. ``` python import spatialproteomics as sp sp.tl.cellpose(sdata, channel="DAPI") # writes sdata.labels["segmentation"] sp.pp.threshold(sdata, quantile=[0.2, 0.6, 0.6, 0.6, 0.6]) sp.pp.add_quantification(sdata, func="intensity_mean") # writes sdata.tables["table"] sp.pp.transform_expression_matrix(sdata, method="arcsinh") sp.la.predict_cell_types_argmax(sdata, {"CD4": "T_h", "CD8": "T_tox"}) # marker -> cell type; adds obs["celltype"] ``` Convert between the two backends with `sp.read_from_spatialdata(sdata)` and `ds.tl.convert_to_spatialdata()`. The xarray backend offers the full feature set; the spatialdata backend covers the core pipeline. ## Common pitfalls - The accessors only exist after `import spatialproteomics`, including on objects loaded with `xr.open_zarr`. - Accessor methods return new objects. `ds.pp.threshold(quantile=0.5)` without reassignment has no effect. - `pp.threshold` overwrites `_image`. Keep a copy (`_raw_image`) first if you want to compare raw and thresholded images. - `la.predict_cell_types_argmax` and `la.predict_cell_subtypes` take **marker -\> cell type** mappings, e.g. `{"CD4": "T_h", "CD8": "T_tox", "FOXP3": "T_reg", "BCL6": "T_fh"}`, not the reverse. - `ds.pp[a:b, c:d]` slices **x first, then y**, while the image array is stored as `(channels, y, x)`. Slices select by coordinate value (like `xarray.Dataset.sel`), so the end point is included. - Call `drop_encoding()` before `to_zarr()` to avoid xarray encoding errors. - Optional tools must be installed separately: `cellpose`, `stardist`, `instanseg-torch`, `astir` and `spotiflow` (or all at once with `pip install "spatialproteomics[all]"`). - Many methods require earlier steps: e.g. `la.predict_cell_types_argmax` needs `_intensity`, which needs `_segmentation`. Inspect the object (`print(ds)`) to see which layers are present. ## Examples ### Cell type QC plot Thresholding, argmax cell type prediction and a side-by-side plot of raw image, thresholded image and resulting cell types per marker, for a human to check the thresholds. ``` python import matplotlib.pyplot as plt import xarray as xr import spatialproteomics # noqa: F401 (registers the .pp / .la / .pl accessors) # --- Inputs (replace with your own) ----------------------------------------- # marker -> display color marker_colors = {"CD3": "red", "CD20": "green", "CD68": "blue"} # marker -> cell type it identifies marker_to_celltype = {"CD3": "T cell", "CD20": "B cell", "CD68": "Macrophage"} # cell type -> display color celltype_colors = {"T cell": "red", "B cell": "green", "Macrophage": "blue"} # channel -> absolute intensity threshold thresholds = {"CD3": 0.3, "CD20": 0.25, "CD68": 0.4} def load_dataset(zarr_path, markers): """Open a segmented zarr and keep an untouched copy of the image.""" ds = xr.open_zarr(zarr_path).pp[["DAPI"] + list(markers)] return ds.pp.add_layer(ds["_image"].values.copy(), key_added="_raw_image") def predict_cell_types(ds): """Threshold, quantify, transform and assign one cell type per cell.""" ds = ( ds.pp.threshold(intensity=list(thresholds.values()), channels=list(thresholds.keys())) .pp.add_quantification(func="intensity_mean") .pp.transform_expression_matrix(method="arcsinh") .la.predict_cell_types_argmax(marker_to_celltype) ) return ds.la.set_label_colors(celltype_colors.keys(), celltype_colors.values()).pl.autocrop() def plot_markers_and_celltypes(ds, scaling=5): """Overview row, then one row per marker: raw + segmentation, raw, thresholded, cell type.""" markers, colors = list(marker_colors.keys()), list(marker_colors.values()) nrows, ncols = len(markers) + 1, 4 fig, axes = plt.subplots(nrows, ncols, figsize=(ncols * scaling, nrows * scaling)) # Overview row ds.pp["DAPI"].pl.colorize("gold").pl.show(render_segmentation=True, ax=axes[0, 0]) ds.pp[markers].pl.colorize(colors, layer_key="_raw_image").pl.show(ax=axes[0, 1]) ds.pp[markers].pl.colorize(colors).pl.show(ax=axes[0, 2]) ds.pl.show(render_image=False, render_labels=True, ax=axes[0, 3]) for ax, title in zip(axes[0], ["Segmentation", "Raw", "Thresholded", "Cell types"]): ax.set_title(title) # One row per marker for row, (marker, color) in zip(axes[1:], marker_colors.items()): celltype = marker_to_celltype[marker] raw = ds.pp[marker].pl.colorize(color, layer_key="_raw_image") try: # outline only the cells of the matching type raw.la[celltype].pl.show(render_segmentation=True, ax=row[0]) except ValueError: # no cells of this type raw.pl.show(ax=row[0]) raw.pl.show(ax=row[1]) ds.pp[marker].pl.colorize(color).pl.show(ax=row[2]) try: ds.la[celltype].pl.show(render_image=False, render_labels=True, ax=row[3]) except ValueError: pass # leave the panel empty for ax, title in zip(row, [marker, marker, marker, celltype]): ax.set_title(title) for ax in axes.flat: ax.axis("off") plt.tight_layout() plt.show() if __name__ == "__main__": ds = load_dataset("path/to/sample.zarr", marker_colors.keys()) ds = predict_cell_types(ds) plot_markers_and_celltypes(ds) ``` ### Binarization QC plot Binarizing markers per cell and plotting each step (raw image, pixel-thresholded image, cells called positive) for a human to check the two thresholds. ``` python import matplotlib.pyplot as plt import xarray as xr import spatialproteomics as sp # --- Inputs (replace with your own) ----------------------------------------- # marker -> (absolute pixel intensity threshold, minimum fraction of positive pixels per cell) # Plotting order follows the order of this dict. thresholds = { "Ki-67": (0.3, 0.1), "CD8": (0.25, 0.2), "PD-1": (0.4, 0.15), } def load_dataset(zarr_path): """Open a segmented zarr and keep an untouched copy of the image.""" ds = xr.open_zarr(zarr_path) ds = ds.pp.add_layer(ds["_image"].values.copy(), key_added="_raw_image") return ds.pl.autocrop() def binarize_marker(ds, marker, intensity_threshold, perc_pos_threshold): """Threshold pixels, compute % positive pixels per cell, then binarize cells.""" return ( ds.pp[marker] .pp.threshold(intensity=intensity_threshold) .pp.add_quantification(func=sp.percentage_positive, key_added="_percentage_positive") .la.threshold_labels({marker: perc_pos_threshold}, layer_key="_percentage_positive") ) def plot_binarization(ds, scaling=5): """One row per marker: raw image, thresholded image, positive cells outlined.""" nrows, ncols = len(thresholds), 3 fig, axes = plt.subplots(nrows, ncols, figsize=(ncols * scaling, nrows * scaling), squeeze=False) for row, (marker, (intensity_thr, perc_pos_thr)) in zip(axes, thresholds.items()): ds.pp[marker].pl.colorize("white", layer_key="_raw_image").pl.show(ax=row[0]) ds_bin = binarize_marker(ds, marker, intensity_thr, perc_pos_thr) ds_bin.pp[marker].pl.colorize("white").pl.show(ax=row[1]) # keep only cells called positive (`{marker}_binarized == 1`) and outline them ( ds_bin.pp[marker] .pp.filter_by_obs(f"{marker}_binarized", func=lambda x: x == 1) .pl.colorize("white") .pl.show( render_segmentation=True, segmentation_kwargs={"alpha": 0.5, "colors": ["gold"]}, ax=row[2], ) ) for ax, title in zip(row, ["Raw", "Thresholded", "Binarized"]): ax.set_title(f"{marker} – {title}") for ax in axes.flat: ax.axis("off") plt.tight_layout() plt.show() if __name__ == "__main__": ds = load_dataset("path/to/sample.zarr") plot_binarization(ds) ``` --- # Installation Install `spatialproteomics` with ``` pip install spatialproteomics ``` ## Extras Some of the functionalities depend on external tools (e. g. `stardist`, `cellpose`, `spatialdata`). Since having lots of external tools also leads to lots of dependencies, installing all of them into one environment is prone to cause dependency conflicts. Depending on your needs, you can choose which tools you want to install together with the `spatialproteomics` core like this: ``` pip install "spatialproteomics[stardist]" pip install "spatialproteomics[cellpose]" pip install "spatialproteomics[spatialdata]" pip install "spatialproteomics[astir]" pip install "spatialproteomics[spotiflow]" pip install "spatialproteomics[all]" ``` The latter one will install all external packages automatically. You can also combine dependencies like this: ``` pip install "spatialproteomics[cellpose,spatialdata]" ``` ## Deprecated Packages `Mesmer` only supported `python<=3.10`, which has reached its end of life in October 2026. It has hence been removed from `spatialproteomics`. If you want to run segmentation using `mesmer`, please downgrade to `python==3.10` and `spatialproteomics==0.8.6`. --- # Example Workflow with spatialproteomics Welcome to spatialproteomics! In this notebook, we will go through an example workflow by looking at the following steps: 1. Reading in a highly multiplexed image and creating a spatialproteomics object 2. Performing basic image processing steps to boost the signal-to-noise ratio 3. Performing cell segmentation using _cellpose_ 4. Quantifying protein expression per cell 5. Predicting cell types with a simple argmax technique 6. Plotting the results *[Figure: Spatialproteomics orchestrates analysis workflows for highly multiplexed fluorescence images.]* ```python import spatialproteomics as sp from skimage.io import imread import matplotlib.pyplot as plt from scipy.signal import medfilt2d import pandas as pd import time start_time = time.time() ``` ## 1. Getting Started Before we can get started with _spatialproteomics_, we first need to read in the image. This can be done with the _skimage_ or _tifffile_ libraries. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python # reading in the image and displaying its shape image = imread("../../data/input.tiff") image.shape ``` Output: ```text (5, 500, 500) ``` We can see that our image has 5 channels and has a size of 500 by 500 pixels. Next, we want to create a _spatialproteomics_ object and store the image within the object. To do this, we also need to feed in the names of the channels. ```python sp_object = sp.load_image_data(image, channel_coords=["DAPI", "CD4", "CD8", "FOXP3", "BCL6"]) ``` Instead of manually typing out all channels in the form of a list, you can also simply read it in from a file, and then provide the resulting list to spatialproteomics. ```python markers = pd.read_csv("../../data/markerlist.txt", header=None)[0].tolist() markers ``` Output: ```text ['DAPI', 'CD4', 'CD8', 'FOXP3', 'BCL6'] ``` ```python # loading the markers into the object sp_object = sp.load_image_data(image, channel_coords=markers) ``` `sp.load_image_data()` returns an xarray object that we can simply inspect by calling in an jupyter cell. Note that the image is stored as the data variable `_image`. ```python sp_object ``` Output: ```text Size: 3MB Dimensions: (channels: 5, y: 500, x: 500) Coordinates: * channels (channels) Size: 3MB Dimensions: (channels: 5, y: 500, x: 500, cells: 2278, features: 2) Coordinates: * channels (channels) 20) & (x < 150)) # plotting the result _ = sp_object.pp["DAPI"].pp[100:200, 100:200].pl.show(render_segmentation=True) ``` ```python # expanding the masks sp_object = sp_object.pp.grow_cells(iterations=2) ``` ```python # plotting the resulting segmentation masks fig, ax = plt.subplots(1, 2, figsize=(10, 5)) _ = sp_object.pp["DAPI"].pp[100:200, 100:200].pl.show(render_segmentation=True, ax=ax[0]) _ = sp_object.pp[["CD4", "CD8", "FOXP3", "BCL6"]].pp[100:200, 100:200].pl.show(render_segmentation=True, ax=ax[1]) ``` ## 4. Quantifying Protein Expression per Cell Now that we have sensible segmentation masks, we can quantify the protein expression in each cell. There are multiple ways to do this, but taking the median intensity and then applying an arcsinh-transform has been proven to work pretty well. ```python sp_object = sp_object.pp.add_quantification(func="intensity_mean").pp.transform_expression_matrix(method="arcsinh") sp_object ``` Output: ```text Size: 3MB Dimensions: (channels: 5, y: 500, x: 500, cells: 2260, features: 2) Coordinates: * channels (channels) Size: 3MB Dimensions: (labels: 4, la_props: 2, channels: 5, y: 500, x: 500, cells: 2260, features: 3) Coordinates: * labels (labels) int64 32B 1 2 3 4 * la_props (la_props) # Example Workflow with spatialproteomics and spatialdata Welcome to spatialproteomics, the toolbox to analyze highly multiplexed fluorescence image data! During the processing of such data, you will generate different data structures, which all share some common dimensions. For example, the multiplexed imaging data (channels, x, y) share spatial dimensions with their corresponding segmentation masks (x, y). Similarly, the instance labels in a segmentation mask can be viewed as a dimension that aligns with the expression matrix (cells, channels) quantifying protein expression for each cell. Therefore, a desirable data structure for multiplexed imaging data should keep such shared dimensions consistent, ensuring that changes to one component automatically update the others to maintain consistency. For example, when subsetting a region of an image, the segmentation mask should be subset accordingly, and the expression matrix should only retain cells within this region. To enable all of this, the core functionality in `spatialproteomics` is based on `xarray`. The image below illustrates these concepts. *[Figure: The spatialproteomics data structure enables synchronized subsetting across shared dimensions.]* Recently, the `spatialdata` format has emerged as a popular alternative, promising standardized storage solution and hence better interoperability with other tools in the field. `Spatialproteomics` integrates with `spatialdata` in two ways: 1. If you want to use the full functionality of `spatialproteomics`, it is recommended that you use the `xarray` backend first (which is used in all tutorials but this one), and then export to `spatialdata` using `tl.convert_to_spatialdata()`. Refer to the FAQ for more information on this. 2. Alternatively, `spatialproteomics` also offers the possibility to work with `spatialdata` objects directly. The syntax for this closely resembles the one commonly used in `squidpy` and other tools which are part of the `scverse` ecosystem. In this notebook, we will go through an example workflow that is fully based on this `spatialdata`-based backend. We will be looking at the following steps: 1. Reading in a highly multiplexed image and creating a spatialdata object 2. Performing basic image processing steps to boost the signal-to-noise ratio 3. Performing cell segmentation using _cellpose_ 4. Quantifying protein expression per cell 5. Predicting cell types with a simple argmax technique 6. Plotting the results with spatialdata-plot **Important:** In order to run `spatialproteomics` with the `spatialdata` backend, it is important that `spatialdata` is installed. You can either do this manually, or install `spatialproteomics` with the additional backend by using `pip install "spatialproteomics[spatialdata]"`. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics as sp import spatialdata as sd import spatialdata_plot from skimage.io import imread import matplotlib import matplotlib.pyplot as plt from scipy.signal import medfilt2d import pandas as pd import numpy as np print(f"Spatialdata version: {sd.__version__}") print(f"Spatialdata-plot version: {spatialdata_plot.__version__}") ``` Output: ```text Spatialdata version: 0.7.2 Spatialdata-plot version: 0.3.3 ``` ## 1. Getting Started Before we can get started with _spatialproteomics_, we first need to read in the image. You can do this using the `spatialdata-io` package. Alternatively, you can construct your own _spatialdata_ object from scratch, which is what we will do here. ```python # reading in the image and displaying its shape image = imread("../../data/input.tiff") image.shape ``` Output: ```text (5, 500, 500) ``` We can see that our image has 5 channels and has a size of 500 by 500 pixels. Next, we want to create a _spatialdata_ object and store the image within the object. To do this, we also need to feed in the names of the channels. ```python image = sd.models.Image2DModel.parse( image, transformations=None, dims=("c", "x", "y"), c_coords=["DAPI", "CD4", "CD8", "FOXP3", "BCL6"] ) sdata = sd.SpatialData(images={"image": image}) sdata ``` Output: ```text SpatialData object └── Images └── 'image': DataArray[cyx] (5, 500, 500) with coordinate systems: ▸ 'global', with elements: image (Images) ``` ## 2. Image Processing Highly multiplexed fluorescence imaging techniques frequently suffer from poor signal-to-noise ratio. To alleviate this problem, you can threshold out low intensity pixels, thereby boosting the contrast of the image. While there are automated methods to determine the thresholds for such operations, it is difficult to come up with one that works in all cases. Here, we therefore set the thresholds based on manual inspection. For more details, check the notebook on `Image Processing`. Let's start by looking at the data set in the current form. ```python # setting the same colors as in the other notebook for consistency colors = ["#e6194B", "#3cb44b", "#ffe119", "#4363d8", "#f58231"] ``` ```python # here, we use the spatialdata-plot library for visualization sdata.pl.render_images(palette=colors, channels_as_legend=True).pl.show() ``` This is a bit much to look at. First, we want to check only the DAPI channel, to get a feel for the signal-to-noise ratio. ```python sdata.pl.render_images(channel="DAPI").pl.show() ``` This looks promising. What about the other channels? ```python sdata.pl.render_images(channel=["CD4", "CD8", "FOXP3", "BCL6"], palette=colors[:-1], channels_as_legend=True).pl.show() ``` Looking good, but we can make the signal a bit clearer by performing some image processing. In the following codeblock, we first threshold the channels by some percentile (so every value below that percentile gets set to 0). We then apply a 2D median filter with a kernel size of 3 to apply some smoothing. ```python # the percentiles by which the channels will be thresholded percentiles = [0.2, 0.6, 0.6, 0.6, 0.6] # thresholding and smoothing the data # note that by default, this performs in-place operations, you can avoid this by setting copy=True in any of the methods sp.pp.threshold(sdata, quantile=percentiles) sp.pp.apply(sdata, medfilt2d, kernel_size=3) ``` Note that `spatialdata` warned us that the key `image` already exists and will be overwritten in memory. This is because by default, the operations on the `spatialdata` object are performed in-place to conserve memory. If you do not wish to perform operations in-place, you can simply set `copy=True` in any of the methods. This will keep the original object unaltered, and return a new, modified object instead. Now that we got that out of the way, let's plot the new object. ```python fig, ax = plt.subplots(1, 2, figsize=(10, 5)) sdata.pl.render_images(channel="DAPI").pl.show(ax=ax[0]) sdata.pl.render_images(channel=["CD4", "CD8", "FOXP3", "BCL6"], palette=colors[:-1], channels_as_legend=True).pl.show( ax=ax[1] ) ``` ## 3. Cell Segmentation Great, this is our preprocessing done. Next, we can perform cell segmentation. Since we only have a universal nuclear marker at hand (and no universal membrane marker), we will segment the nuclei and then simply extend the segmentation masks by two pixels in every direction. We are going to use _cellpose_ for this purpose, which is implemented in the tool (tl) module. Note that this requires you to install cellpose first, for example by running `pip install cellpose` in your terminal. ```python sp.tl.cellpose(sdata, channel="DAPI") sdata ``` Output: ```text SpatialData object ├── Images │ └── 'image': DataArray[cyx] (5, 500, 500) └── Labels └── 'segmentation': DataArray[yx] (500, 500) with coordinate systems: ▸ 'global', with elements: image (Images), segmentation (Labels) ``` Looking at the object, you will realize that a new layer called `segmentation` has appeared. We can plot this over the original image with `spatialdata-plot`. Let's also zoom in a little bit to get a clearer picture. ```python sdata_zoomed = sdata.query.bounding_box( axes=["x", "y"], min_coordinate=[300, 400], max_coordinate=[400, 500], target_coordinate_system="global", ) sdata_zoomed.pl.render_images(channel=["DAPI"]).pl.render_labels( outline_alpha=1, fill_alpha=0, outline_color="white" ).pl.show() ``` ## 4. Quantifying Protein Expression per Cell There are a couple of issues which can arise from segmentation. One is that sometimes very small cells get segmented, which are likely artifacts. We can hence filter cells that are too small or too big. In addition, we will grow the masks by two pixels in each direction to try to capture cytoplasm and membrane. With the `spatialdata` backend, it is important that we first quantify the protein expression in each cell. There are multiple ways to do this, but taking the median intensity and then applying an arcsinh-transform has been proven to work pretty well. ```python sp.pp.add_quantification(sdata, func="intensity_mean") sp.pp.transform_expression_matrix(sdata, method="arcsinh") sdata ``` Output: ```text SpatialData object ├── Images │ └── 'image': DataArray[cyx] (5, 500, 500) ├── Labels │ └── 'segmentation': DataArray[yx] (500, 500) └── Tables └── 'table': AnnData (2346, 5) with coordinate systems: ▸ 'global', with elements: image (Images), segmentation (Labels) ``` As you can see, this introduced a new layer called `table`. This is an `anndata` object that contains the quantification of each protein for each cell. ```python sdata.tables["table"] ``` Output: ```text AnnData object with n_obs × n_vars = 2346 × 5 obs: 'id', 'region' uns: 'spatialdata_attrs' ``` ```python sp.pp.add_observations(sdata, "area") ``` ```python # checking the obs slot in our anndata object sdata.tables["table"].obs.head() ``` Output: ```text id region area Cell_1 1 segmentation 38.0 Cell_2 2 segmentation 66.0 Cell_3 3 segmentation 73.0 Cell_4 4 segmentation 46.0 Cell_5 5 segmentation 34.0 ``` ```python # filtering out cells with less than 20 or more than 150 pixels sp.pp.filter_by_obs(sdata, "area", func=lambda x: (x > 20) & (x < 150)) ``` After the filtering, you can see that there are less cells both in the segmentation mask as well as in the anndata object. ```python print( len(np.unique(sdata.labels["segmentation"].values)) - 1 ) # -1 because the background (0) is also counted otherwise print(sdata.tables["table"]) ``` Output: ```text 2330 AnnData object with n_obs × n_vars = 2330 × 5 obs: 'id', 'region', 'area' uns: 'spatialdata_attrs' ``` Let's also note down the average cell area before and after growing the masks. ```python sdata.tables["table"].obs["area"].values.mean() ``` Output: ```text np.float64(75.60858369098712) ``` ```python # expanding the masks sp.pp.grow_cells(sdata, iterations=2) # adding the area back into the object sp.pp.add_observations(sdata, "area") ``` ```python sdata.tables["table"].obs["area"].values.mean() ``` Output: ```text np.float64(93.99656652360515) ``` The average area has indeed increased. Let's also check what the new segmentation masks look like on top of DAPI or our other markers. ```python # plotting the resulting segmentation masks sdata_zoomed = sdata.query.bounding_box( axes=["x", "y"], min_coordinate=[300, 400], max_coordinate=[400, 500], target_coordinate_system="global", ) fig, ax = plt.subplots(1, 2, figsize=(10, 5)) sdata_zoomed.pl.render_images(channel=["DAPI"]).pl.render_labels( outline_alpha=1, fill_alpha=0, outline_color="white" ).pl.show(ax=ax[0]) sdata_zoomed.pl.render_images(channel=["CD4", "CD8", "FOXP3", "BCL6"], channels_as_legend=True).pl.render_labels( outline_alpha=1, fill_alpha=0, outline_color="white" ).pl.show(ax=ax[1]) ``` ## 5. Cell Type Prediction Finally, we can move on to cell type prediction. There are several ways to predict cell types. Since we thresholded our data beforehand, we can simply take the argmax of the cell type specific channels to get an idea of the cell types we are looking at. Methods related to cell type prediction are all implemented in the label (la) module. ```python # this dictionary maps from cell types to markers marker_dict = {"CD4": "T_h", "CD8": "T_tox", "FOXP3": "T_reg", "BCL6": "T_fh"} sp.la.predict_cell_types_argmax(sdata, marker_dict) ``` ```python # checking what the "obs" in the anndata object look like sdata.tables["table"].obs.head() ``` Output: ```text id region area celltype Cell_1 1 segmentation 39.0 T_tox Cell_2 2 segmentation 75.0 T_h Cell_3 3 segmentation 116.0 T_h Cell_4 4 segmentation 49.0 T_tox Cell_5 5 segmentation 38.0 T_tox ``` As you can see, this added a column called `celltype` to our `anndata.obs`. ## 6. Plotting Let's do some plotting of the predicted cell types. ```python fig, ax = plt.subplots(1, 2, figsize=(10, 5)) sdata.pl.render_images(channel=["CD4", "CD8", "FOXP3", "BCL6"], channels_as_legend=True).pl.show( legend_loc="upper left", ax=ax[0] ) sdata.pl.render_labels(color="celltype").pl.show(legend_loc="upper left", ax=ax[1]) ``` ## 7. Saving Results to Disk And this is how easy it can be to perform analysis of highly multiplexed immunofluorescence images! To store your processed object to disk, use `sdata.write("your_file_path")`. If you have any additional questions, check out the other notebooks for details, or write an issue on [GitHub](https://github.com/sagar87/spatialproteomics/issues). --- # Plotting Spatialproteomics has a variety of plotting functions, enabling the plotting of intensities, segmentation masks, and predicted cell type labels. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics import matplotlib.pyplot as plt import xarray as xr ``` ```python # loading in a data set ds_unprocessed = xr.open_zarr("../../data/LN_24_1_unprocessed.zarr") ``` ```python # having a look at the dataset ds_unprocessed ``` Output: ```text Size: 576MB Dimensions: (channels: 56, y: 3000, x: 3000, cells: 18134, features: 2) Coordinates: * channels (channels) _obs (cells, features) float64 290kB dask.array _segmentation (y, x) int64 72MB dask.array ``` ## Plotting Channel Intensities The slicing capabilities of `xarray` make data visualization very easy. The general pattern is here to first select the channels, image regions and labels of interest, and then to use the plotting functions that are available via the `.pl` accessor. Images can be shown by using the `pl.show()` method. ```python # selecting the DAPI channel and showing the result with pl.show() _ = ds_unprocessed.pp["DAPI"].pl.show() ``` If we want different colors, we can specify those using the `pl.colorize` method. ```python _ = ds_unprocessed.pp["DAPI"].pl.colorize(["orange"]).pl.show() ``` The `pl.colorize` function allows us to make image overlays quickly by first selecting the channels of interest. We illustrate this by first selecting the CD4 and CD8 channels and then assigning the colors. ```python _ = ds_unprocessed.pp[["CD3", "PAX5"], 1000:1500, 500:1000].pl.colorize(colors=["orange", "deepskyblue"]).pl.show() ``` All plotting functions accept `matplotlib.Axes` via the `ax` argument. This enables to easily modify plots using the normal `matplotlib` API. For example, let's say we want to plot multiple channels in different subplots. ```python # creating two subplots fig, ax = plt.subplots(1, 2, figsize=(16, 8)) # populating the first subplot # note how we set ax=ax[0] to tell the plotting function to use the first subplot _ = ds_unprocessed.pp["DAPI", 1000:1500, 500:1000].pl.colorize(colors=["red"]).pl.show(ax=ax[0]) # populating the second subplot _ = ( ds_unprocessed.pp[["CD3", "PAX5"], 1000:1500, 500:1000] .pl.colorize(colors=["orange", "deepskyblue"]) .pl.show(ax=ax[1]) ) ``` By default, a legend showing the colors corresponding to the channels is shown. You can disable this with the `legend_image` argument. In addition, matplotlib functions can be used to further alter the plot, e. g. to remove ticklabels. The code snipped below also highlights the `pl.autocrop()` function, which can be useful to remove empty space around a tissue microarray. ```python fig = plt.figure(figsize=(8, 8)) _ = ds_unprocessed.pp["DAPI"].pl.autocrop().pl.colorize(colors=["red"]).pl.show(legend_image=False) plt.xticks([]) plt.yticks([]) plt.tight_layout() ``` ## Rendering Segmentation Masks If your spatialproteomics object contains a segmentation mask, you can simple plot it by setting `render_segmentation` to true. ```python fig = plt.figure(figsize=(8, 8)) _ = ( ds_unprocessed.pp[["CD3", "PAX5"], 1000:1500, 500:1000] .pl.colorize(colors=["orange", "deepskyblue"]) .pl.show(render_segmentation=True) ) ``` ## Rendering Obs Sometimes, you might want to visualize specific features, such as the cell size. With `spatialproteomics`, you can do this while retaining the cell shapes obtained by the segmentation. ```python fig = plt.figure(figsize=(8, 8)) _ = ( ds_unprocessed.pp.add_observations("area") .pl.autocrop() .pl.render_obs(feature="area", cmap="plasma") .pl.imshow(legend_obs=True) ) ``` ```python fig = plt.figure(figsize=(8, 8)) _ = ds_unprocessed.pl.autocrop().pl.render_obs(feature="centroid-0", cmap="plasma").pl.imshow(legend_obs=True) ``` ## Rendering Labels Similar to `render_segmentation`, `pl.show()` also has an argument called `render_labels`, which can be used to render cell type labels. We can simultaneously deactivate the rendering of the intensities to get an image of only the labels. ```python # reading in a dataset containing labels ds = xr.open_zarr("../../data/LN_24_1.zarr") # setting colors for nicer visualizations celltype_colors = { "B cell": "#5799d1", "T cell": "#ebc850", "Myeloid cell": "#de6866", "Dendritic cell": "#4cbcbd", "Macrophage": "#bb7cb4", "Stromal cell": "#62b346", "Endothelial cell": "#bf997d", } ds = ds.la.set_label_level("labels_0", ignore_neighborhoods=True).la.set_label_colors( celltype_colors.keys(), celltype_colors.values() ) ``` ```python fig = plt.figure(figsize=(8, 8)) _ = ds.pp[1000:1500, 500:1000].pl.show(render_image=False, render_labels=True) ``` We can combine this with rendering the segmentation masks for a clearer picture of where exactly the cell boundaries lie. ```python fig = plt.figure(figsize=(8, 8)) _ = ds.pp[1000:1500, 500:1000].pl.show(render_image=False, render_segmentation=True, render_labels=True) ``` Using the `.la[]` accessor we can also subset cell labels, rather than displaying all of them. We can also use the `la.set_label_colors()` to set custom colors for our labels. ```python fig, ax = plt.subplots(1, 2, figsize=(16, 8)) _ = ds.pp[["CD3", "PAX5"], 1000:1500, 500:1000].pl.colorize(["orange", "deepskyblue"]).pl.show(ax=ax[0]) _ = ( ds.pp[1000:1500, 500:1000] .la["T cell", "B cell"] .la.set_label_colors(["T cell", "B cell"], ["orange", "deepskyblue"]) .pl.show(render_image=False, render_segmentation=True, render_labels=True, ax=ax[1]) ) ``` ## Rendering Neighborhoods We can also render neighborhoods. How to obtain such neighborhoods is shown in the `Neighborhoods` notebook. ```python # obtaining an object that already contains neighborhoods ds_neighborhood = xr.open_zarr("../../data/LN_24_1.zarr") ``` ```python # setting neighborhood colors ds_neighborhood = ds_neighborhood.nh.set_neighborhood_colors( [f"Neighborhood {x+1}" for x in range(6)], ["lightgray", "deepskyblue", "purple", "gold", "gray", "darkgray"] ) # plotting the neighborhoods fig = plt.figure(figsize=(8, 8)) _ = ds_neighborhood.pl.show(render_image=False, render_neighborhoods=True) ``` ## Annotation of Regions or Cells We can also annotate regions or cells within the plots. ```python # highlighting a region with pl.add_box() fig = plt.figure(figsize=(8, 8)) _ = ( ds_unprocessed.pp[["CD3", "PAX5"], 1000:1500, 500:1000] .pl.colorize(colors=["orange", "deepskyblue"]) .pl.add_box([1230, 1480], [640, 840]) .pl.show() ) ``` ```python # selecting all B cells and annotating them with pl.annotate() fig = plt.figure(figsize=(8, 8)) _ = ( ds.pp[["CD3", "PAX5"], 1500:2000, 1200:1700] .la["B cell"] .pl.colorize(colors=["orange", "deepskyblue"]) .la.set_label_colors(["B cell"], ["deepskyblue"]) .pl.annotate() .pl.show(render_labels=True, label_kwargs={"alpha": 0.0}) ) ``` ## Expert Plotting To gain more control over your plots, you can use the functions `colorize`, `render_segmentation`, and `render_labels` in combination with `imshow` (instead of `show`, which provides a high-level wrapper around these). Those methods allow for more nuanced adjustments of the final plot. Here are a few examples, for more details please consult the documentation of the corresponding methods. ```python _ = ds_unprocessed.pp[["CD3", "PAX5"], 1000:1500, 500:1000].pl.colorize(colors=["orange", "deepskyblue"]).pl.imshow() ``` ```python _ = ds_unprocessed.pl.render_segmentation().pl.imshow() ``` ```python _ = ds.pl.render_labels().pl.imshow() ``` ```python _ = ( ds_unprocessed.pp[["CD3", "PAX5"], 1000:1500, 500:1000] .pl.colorize(colors=["orange", "deepskyblue"]) .pl.render_segmentation() .pl.imshow() ) ``` You can also use scatterplots instad of always plotting the images. ```python # creating a scatter plot of the labels _ = ds.pl.scatter_labels(legend=False) ``` ```python fig, ax = plt.subplots(1, 2, figsize=(10, 20)) _ = ds.pp["CD8"].pl.colorize("blue").pl.show(ax=ax[0]) _ = ds.pl.scatter(feature="CD8_binarized", size=0.1, palette={0: "gray", 1: "blue"}, ax=ax[1]) ``` ```python # creating a scatterplot of continuous variables # adding area of each cell and plotting it ds_tmp = ds.pp.add_observations("area") _ = ds_tmp.pl.scatter(feature="area", scatter_kws={"cmap": "YlOrRd"}) ``` If you have multiple segmentation masks, you can show all of them at the same time like this. ```python _ = ds_unprocessed.pp[["DAPI", "CD11c"]].pp[2000:2500, 2000:2500].pl.colorize(["orange", "deepskyblue"]).pl.show() ``` ```python ds_multisegmentation = ( ds_unprocessed.pp[["DAPI", "CD11c"]] .pp[2000:2500, 2000:2500] .tl.cellpose(key_added="_cellpose_segmentation", channel_settings=[0, 0]) ) ``` ```python plt.figure(figsize=(16, 16)) _ = ds_multisegmentation.pl.colorize(["gold", "deepskyblue"]).pl.show( render_segmentation=True, segmentation_kwargs={"layer_key": "_cellpose_segmentation", "colors": ["white", "red"]} ) ``` We can also plot histograms of any table, for example the observations (e. g. cell area) or the marker distribution. ```python # quick check what our object looks like ds_tmp ``` ```python # plotting a histogram for a feature in the obs table df = ds_tmp.pp.get_layer_as_df("_obs") _ = plt.hist(df["area"], bins=50) ``` ```python # plotting a histogram for a feature in the expression matrix df = ds_tmp.pp.get_layer_as_df("_intensity") _ = plt.hist(df["DAPI"], bins=50) ``` Sidenote: you can also create a Vornoi diagram of your data by simply growing segmentation masks by a specified distance. ```python plt.figure(figsize=(15, 15)) _ = ( ds_tmp.pp["DAPI"] .pp.grow_cells(iterations=20) .pl.autocrop() .pl.colorize("gold", layer_key="_image_raw") .pl.show(render_segmentation=True) ) ``` For visualization, it is also possible to adjust the background color of the image. However, caution is advised when plotting fluorescence intensities on light backgrounds, as the background colors get mixed with the actual fluorescence values. ```python fig, ax = plt.subplots(1, 2, figsize=(20, 10)) # you can pass the background either in the pl.colorize() method or in pl.show() _ = ds_tmp.pp["DAPI"].pl.autocrop().pl.colorize("gold", background="#2e2e2e", layer_key="_image_raw").pl.show(ax=ax[0]) _ = ds_tmp.pp["DAPI"].pl.autocrop().pl.show(render_image=False, render_labels=True, background="lightgray", ax=ax[1]) ``` --- # Segmentation Segmentation describes the process of finding cells in an image. This can be done either on the nucleus or on whole cell level. Spatialproteomics provides wrappers for `stardist` and `cellpose` (`cellpose 3` or `cellpose-SAM`). If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics import pandas as pd import matplotlib.pyplot as plt import xarray as xr ``` ```python # loading the input data ds = xr.open_zarr("../../data/segmentation_example.zarr").pp[1000:1300, 1000:1300].pp.drop_layers(keep="_image") ``` ```python # checking what the object initially looks like ds ``` Output: ```text Size: 367kB Dimensions: (channels: 4, y: 301, x: 301) Coordinates: * channels (channels) ``` ## Whole-Cell Segmentation with Cellpose Cellpose can either perform segmentation on nuclei only or on whole cells, using a nuclear and a membrane marker. In this dataset, we have DAPI as a nuclear marker and Na/K ATPase as a membrane marker. In cellpose 4, the order of these does not matter. You can refer to the [cellpose documentation](https://cellpose.readthedocs.io/en/latest/settings.html) for more details to tweak the model. ```python # running the segmentation on the selected markers ds_cellpose = ds.tl.cellpose(channel=["Na/K ATPase", "DAPI"]) ds_cellpose ``` Output: ```text Size: 555kB Dimensions: (channels: 4, y: 301, x: 301, cells: 281, features: 2) Coordinates: * channels (channels) _segmentation (y, x) uint16 181kB 0 0 0 0 0 0 0 ... 273 273 273 273 273 0 _obs (cells, features) float64 4kB 1.009e+03 ... 1.004e+03 ``` ```python # for plotting, we select the two markers from before, but you could also show the other markers _ = ds_cellpose.pp[["Na/K ATPase", "DAPI"]].pl.colorize(["deepskyblue", "gold"]).pl.show(render_segmentation=True) ``` ## Whole-Cell Segmentation with InstanSeg InstanSeg can use multiple fluorescent markers to segment cells. In this example, we use DAPI and a couple of membrane markers to obtain a whole-cell segmentation. You can refer to the [InstanSeg documentation](https://instanseg.github.io/instanseg/index.html) for more details to tweak the model. ```python # running the segmentation on the selected markers ds_instanseg = ds.tl.instanseg(channel=["DAPI", "CD68", "CD11c"]) ds_instanseg ``` Output: ```text Model fluorescence_nuclei_and_cells version 0.1.1 already downloaded in /g/huber/users/meyerben/notebooks/codex_analysis/2024-03-18_spatialproteomics_package/spatialproteomics/.venv/lib/python3.12/site-packages/instanseg/utils/../bioimageio_models/, loading Requesting default device: cpu Size: 1MB Dimensions: (channels: 4, y: 301, x: 301, cells: 349, features: 2) Coordinates: * channels (channels) _segmentation (y, x) int64 725kB 0 0 0 0 0 13 13 ... 355 355 355 355 355 0 _obs (cells, features) float64 6kB 1.002e+03 ... 1.288e+03 ``` ```python # for plotting, we select the three markers used in the segmentation, but you could also show other markers _ = ( ds_instanseg.pp[["DAPI", "CD68", "CD11c"]] .pl.colorize(["gold", "purple", "deepskyblue"]) .pl.show(render_segmentation=True) ) ``` ## Whole-Cell Segmentation with Mesmer Segmentation with `mesmer` is deprecated, since `mesmer` only supported `python<=3.10`, which has reached its end of life in October 2026. If you want to use it, please use `spatialproteomics<=0.8.6`. ## Nuclear Segmentation with StarDist Sometimes, you may not have a membrane stain available. In that case, you can perform nuclear segmentation. This is possible with cellpose (simply call `tl.cellpose(channel='DAPI')`) and StarDist. The following shows and example of how StarDist performs on such a task. ```python # running stardist ds_stardist = ds.tl.stardist(channel="DAPI") ``` Output: ```text Found model '2D_versatile_fluo' for 'StarDist2D'. Loading network weights from 'weights_best.h5'. Loading thresholds from 'thresholds.json'. Using default values: prob_thresh=0.479071, nms_thresh=0.3. ``` ```python # plotting the segmentation fix, ax = plt.subplots(1, 2, figsize=(16, 8)) _ = ds_stardist.pp["DAPI"].pl.colorize("orange").pl.show(ax=ax[0]) _ = ds_stardist.pp["DAPI"].pl.colorize("orange").pl.show(render_segmentation=True, ax=ax[1]) ``` You might have noticed that StarDist tends to segment a lot of small cells. Let's look at the distribution of cell sizes, and then remove all nuclei under a certain size. ```python # getting a histogram of the areas _ = plt.hist(ds_stardist.pp.add_observations("area").pp.get_layer_as_df()["area"], bins=100) ``` ```python # removing cells which have an area below 50 pixels ds_stardist = ds_stardist.pp.add_observations("area").pp.filter_by_obs("area", lambda x: x > 50) # plotting the segmentation fix, ax = plt.subplots(1, 2, figsize=(16, 8)) _ = ds_stardist.pp["DAPI"].pl.colorize("orange").pl.show(ax=ax[0]) _ = ds_stardist.pp["DAPI"].pl.colorize("orange").pl.show(render_segmentation=True, ax=ax[1]) ``` To try to capture full cells, we can now grow the segmentation masks, e. g. by two pixels into each direction. If two masks collide, the masks will stop growing into that direction. ```python ds_stardist_grown = ds_stardist.pp.grow_cells(2) ``` ```python # plotting the segmentation fix, ax = plt.subplots(1, 2, figsize=(16, 8)) _ = ds_stardist.pp["DAPI"].pl.colorize("orange").pl.show(render_segmentation=True, ax=ax[0]) _ = ds_stardist_grown.pp["DAPI"].pl.colorize("orange").pl.show(render_segmentation=True, ax=ax[1]) ax[0].set_title("StarDist segmentation") ax[1].set_title("Grown masks") plt.tight_layout() ``` ## Segmenting Multiple Channels with Cellpose Next to predicting nuclei or whole cells using a universal membrane marker, you can also use cellpose to try to predict other cell types. In the following, we will see how we can use cellpose to segment three different channels independently and then merge the segmentation masks into one mask. There are some hyperparameters in cellpose which can help the prediction. ```python # segmenting on DAPI, CD11c (dendritic cells) and CD68 (macrophages) # note that when merging the masks later, the order matters (channels take precedence by position in the subsetting) # so you should have your largest cells first, and then go down by expected size # we use the channel_settings argument here to tell cellpose that it should perform an independent segmentation for each of our markers # without this argument, cellpose-sam would try to segment the cells based on all three markers combined ds_cellpose = ds.pp[["CD11c", "CD68", "DAPI"]].tl.cellpose(key_added="_cellpose_segmentation", channel_settings=[0, 0]) ``` ```python ds_cellpose ``` Output: ```text Size: 820kB Dimensions: (channels: 3, y: 301, x: 301) Coordinates: * channels (channels) _cellpose_segmentation (channels, y, x) uint16 544kB 0 0 0 ... 311 311 311 ``` ```python # merging the segmentation masks into a single segmentation # the threshold value tells you how much of a mask need to overlap for the less important one to be removed # for more information on this, see next section ds_cellpose_merged = ds_cellpose.pp.merge_segmentation( "_cellpose_segmentation", labels=["Dendritic Cell", "Macrophage", "Lymphocyte"] ) ``` ```python ds_cellpose_merged ``` Output: ```text Size: 1MB Dimensions: (channels: 3, y: 301, x: 301) Coordinates: * channels (channels) _cellpose_segmentation (channels, y, x) uint16 544kB 0 0 0 ... 311 311 311 _merged_segmentation (y, x) uint16 181kB 0 11 11 11 ... 310 310 310 310 ``` ```python # adding the merged segmentation as the main segmentation in the object # this is done so that we can visualize the result with pl.show(render_segmentation=True), without needing to specify which segmentation we want to show ds_cellpose_merged = ds_cellpose_merged.pp.add_segmentation("_merged_segmentation") ``` ```python ds_cellpose_merged ``` Output: ```text Size: 1MB Dimensions: (labels: 3, la_props: 2, channels: 3, x: 301, y: 301, cells: 315, features: 3) Coordinates: * labels (labels) int64 24B 1 2 3 * la_props (la_props) _cellpose_segmentation (channels, y, x) uint16 544kB 0 0 0 ... 311 311 311 ... ``` ```python fig, ax = plt.subplots(1, 2, figsize=(16, 8)) # adding colors to match the channel colors colors = ["orange", "deepskyblue", "darkred"] ds_cellpose_merged = ds_cellpose_merged.la.set_label_colors(["Lymphocyte", "Dendritic Cell", "Macrophage"], colors) # merged masks _ = ( ds_cellpose_merged.pp[["DAPI", "CD11c", "CD68"]] .pl.colorize(colors=colors) .pl.show(render_segmentation=True, ax=ax[0]) ) _ = ds_cellpose_merged.pl.show( render_image=False, render_segmentation=True, render_labels=True, legend_label=True, ax=ax[1] ) ``` ## Fine-Tuning Cellpose Cellpose often benefits from fine-tuning, which is detailed in their [documentation](https://cellpose.readthedocs.io/en/latest/gui.html#training-your-own-cellpose-model). You can load fine-tuned cellpose models into `spatialproteomics` by using `ds.tl.cellpose(channel="your_marker", pretrained_model=model_path)`. ## Merging Segmentation Masks from Different Layers In the previous section, we looked at how to merge masks from a single layer. But what do we do if we have multiple segmentations in different layers? Let us consider an example. We can run cellpose on our membrane stain, and then independently on the macrophages and dendritic cells. However, now our segmentation masks are stored in different layers. The good thing is, we can combine them just as easily by using `pp.merge_segmentation()`. ```python # running cellpose to obtain a bunch of different segmentation masks # these are stored in the layers dapi_seg, macro_seg, and dendritic_seg ds_cellpose = ( ds.tl.cellpose(channel=["DAPI", "Na/K ATPase"], key_added="dapi_seg") .tl.cellpose(channel="CD68", key_added="macro_seg") .tl.cellpose(channel="CD11c", key_added="dendritic_seg") ) ds_cellpose ``` Output: ```text Size: 911kB Dimensions: (channels: 4, y: 301, x: 301) Coordinates: * channels (channels) dapi_seg (y, x) uint16 181kB 0 0 0 0 0 0 0 ... 281 281 281 281 281 0 macro_seg (y, x) uint16 181kB 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 dendritic_seg (y, x) uint16 181kB 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 ``` ```python # plotting the individual masks fig, ax = plt.subplots(1, 3, figsize=(15, 5)) _ = ( ds_cellpose.pp[["DAPI", "Na/K ATPase"]] .pl.colorize(["gold", "deepskyblue"]) .pl.show(render_segmentation=True, segmentation_kwargs={"layer_key": "dapi_seg"}, ax=ax[0]) ) _ = ( ds_cellpose.pp[["DAPI", "CD68"]] .pl.colorize(["gold", "red"]) .pl.show(render_segmentation=True, segmentation_kwargs={"layer_key": "macro_seg"}, ax=ax[1]) ) _ = ( ds_cellpose.pp[["DAPI", "CD11c"]] .pl.colorize(["gold", "deepskyblue"]) .pl.show(render_segmentation=True, segmentation_kwargs={"layer_key": "dendritic_seg"}, ax=ax[2]) ) ``` In brief, `pp.merge_segmentation()` works like this: 1. Assume two segmentation masks, $S_1$ (higher priority, e.g., macrophages) and $S_2$ (lower priority, e.g., lymphocytes) 2. We start by copying the segmentation regions from $S_2$ into $S_1$, but only for pixels where $S_1$ has background 3. Some $S_2$ cells may now be partially overwritten by $S_1$ cells — we drop any $S_2$ cell where less than `threshold1` fraction of its original area survived (e.g. with `threshold1=0.9`, an $S_2$ cell that lost more than 10% of its area to $S_1$ is discarded) 4. Of the remaining $S_2$ cells, we additionally drop any that overlap with `threshold2` or more of any single $S_1$ cell's area — this prevents $S_2$ cells from being retained when they largely encapsulate an $S_1$ cell (e.g. with `threshold2=0.9`, an $S_2$ cell is discarded if it covers more than 90% of any $S_1$ cell) 5. The surviving $S_1$ and $S_2$ regions are combined into a final mask with sequential labels If you have more than two segmentation mask, this process is run iteratively until all segmentations have been merged into one. Depending on what segmentation masks you have available, you might want to tweak `threshold1` and `threshold2` to be suitable to your data. ```python # merging the segmentation masks and making the merged segmentation the default one ds_cellpose = ds_cellpose.pp.merge_segmentation( layer_key=["dendritic_seg", "macro_seg", "dapi_seg"], labels=["Dendritic Cell", "Macrophage", "Lymphocyte"], threshold1=0.6, threshold2=0.9, ).pp.add_segmentation("_merged_segmentation") ``` ```python fig, ax = plt.subplots(1, 2, figsize=(16, 8)) # adding colors to match the channel colors colors = ["orange", "deepskyblue", "darkred"] ds_cellpose = ds_cellpose.la.set_label_colors(["Lymphocyte", "Dendritic Cell", "Macrophage"], colors) # merged masks _ = ds_cellpose.pp[["DAPI", "CD11c", "CD68"]].pl.colorize(colors=colors).pl.show(render_segmentation=True, ax=ax[0]) _ = ds_cellpose.pl.show(render_image=False, render_segmentation=True, render_labels=True, legend_label=True, ax=ax[1]) ``` --- # Image Processing This notebook shows some exemplary image processing algorithms. Note that plotting methods should always be called after preprocessing modules. `ds.pp.threshold(0.99).pl.colorize('red').pl.show()` will work, whereas `ds.pl.colorize('red').pp.threshold(0.99).pl.show()` will not show the thresholded image correctly. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics import matplotlib.pyplot as plt import xarray as xr ``` ```python # loading in a data set and performing some formatting for convenience ds = xr.open_zarr("../../data/LN_24_1_unprocessed.zarr") # focusing on just one channel to illustrate the concepts better ds_single_channel = ds.pp["IKZF3"] ds_multichannel = ds.pp[["PAX5", "CD3"]] ds ``` Output: ```text Size: 576MB Dimensions: (channels: 56, y: 3000, x: 3000, cells: 18134, features: 2) Coordinates: * cells (cells) int64 145kB 1 2 3 4 5 ... 18131 18132 18133 18134 * channels (channels) _obs (cells, features) float64 290kB dask.array _segmentation (y, x) int64 72MB dask.array ``` ```python plt.figure(figsize=(10, 10)) _ = ds_single_channel.pl.colorize("gold").pl.show() ``` ## Thresholding The marker above seems to have worked fairly well, but you can observe some unspecific binding. To boost the signal-to-noise ratio, we can threshold the image. This basically means that every value below that threshold gets set to zero. In addition, the other intensities get shifted downwards by your specified value (you can deactivate this behavior by setting `pp.threshold(shift=False)`, which will only set values below the threshold to 0 and keep the other values unchanged). You can either threshold using an absolute intensity, or using a quantile value. ```python # thresholding with absolute intensity fig, ax = plt.subplots(1, 5, figsize=(50, 10)) for i in range(5): _ = ds_single_channel.pp.threshold(intensity=10 * i).pl.colorize("gold").pl.show(ax=ax[i]) ax[i].set_title(f"Absolute threshold: {10 * i}") ``` ```python # thresholding with quantiles fig, ax = plt.subplots(1, 5, figsize=(50, 10)) for i, quantile in enumerate([0.8, 0.9, 0.95, 0.99, 0.999]): _ = ds_single_channel.pp.threshold(quantile=quantile).pl.colorize("gold").pl.show(ax=ax[i]) ax[i].set_title(f"Quantile threshold: {quantile}") ``` You can also use different thresholds for different channels. ```python fig = plt.figure(figsize=(10, 10)) _ = ds_multichannel.pp.threshold(quantile=[0.9, 0.8]).pl.colorize(["gold", "deepskyblue"]).pl.show() ``` ```python fig = plt.figure(figsize=(10, 10)) _ = ds_multichannel.pp.threshold(intensity=[15, 20]).pl.colorize(["gold", "deepskyblue"]).pl.show() ``` After thresholding, a medianfilter can be useful to reduce the noise further. We can use the `pp.apply()` method to apply any preprocessing function of our choice. ```python from scipy.signal import medfilt2d fig, ax = plt.subplots(1, 3, figsize=(30, 10)) _ = ds_single_channel.pl.colorize("gold").pl.show(ax=ax[0]) _ = ds_single_channel.pp.threshold(quantile=0.9).pl.colorize("gold").pl.show(ax=ax[1]) _ = ( ds_single_channel.pp.threshold(quantile=0.9) .pp.apply(medfilt2d, kernel_size=3) .pl.colorize("gold") .pl.show(ax=ax[2]) ) ax[0].set_title("Raw image") ax[1].set_title("Thresholded") ax[2].set_title("Thresholded and median filtered") plt.show() ``` ## Combining Channels Sometimes, you might want to combine different channels into a single one. For example, if you have different membrane channels, you could merge them into a single channel and then perform segmentation on this new channel. This can be done with `pp.merge_channels()`. You can specify if you want to normalize channels before merging with `normalize` (defaults to True) and set a method using the `method` argument. ```python # let's consider the channels CD68 (macrophages) and CD11c (dendritic cells) _ = ds.pp[["CD68", "CD11c"]].pl.autocrop().pl.colorize(["gold", "deepskyblue"]).pl.show() ``` ```python # merging the channels into one single channel ds_merged = ds.pp.merge_channels(["CD68", "CD11c"], key_added="combined_channel", method="max") _ = ds_merged.pp["combined_channel"].pl.autocrop().pl.show() ``` ## Custom Image Processing With the `pp.apply()` method, you can provide existing (e. g. from `scipy` or `skimage`) or custom functions to your multichannel images. ```python # application to a single channel from scipy.ndimage import gaussian_filter def multiply_array(arr, factor=10): return (arr * factor).clip(0, 255).astype("uint8") fig, ax = plt.subplots(1, 3, figsize=(30, 10)) _ = ds_single_channel.pl.colorize("gold").pl.show(ax=ax[0]) _ = ds_single_channel.pp.apply(func=multiply_array).pl.colorize("gold").pl.show(ax=ax[1]) _ = ds_single_channel.pp.apply(func=gaussian_filter, sigma=5).pl.colorize("gold").pl.show(ax=ax[2]) ax[0].set_title("Raw image") ax[1].set_title("Custom method") ax[2].set_title("Gaussian filter") ``` Output: ```text Text(0.5, 1.0, 'Gaussian filter') ``` ```python # application to multiple channels fig, ax = plt.subplots(1, 3, figsize=(30, 10)) _ = ds_multichannel.pl.colorize(["gold", "deepskyblue"]).pl.show(ax=ax[0]) _ = ds_multichannel.pp.apply(func=multiply_array).pl.colorize(["gold", "deepskyblue"]).pl.show(ax=ax[1]) _ = ds_multichannel.pp.apply(func=gaussian_filter, sigma=5).pl.colorize(["gold", "deepskyblue"]).pl.show(ax=ax[2]) ax[0].set_title("Raw image") ax[1].set_title("Custom method") ax[2].set_title("Gaussian filter") ``` Output: ```text Text(0.5, 1.0, 'Gaussian filter') ``` ```python # total variation filter. The idea is that areas with large gradients are smoothed out, as this implies rapid fluctuation of pixel intensity values (noise). from skimage.restoration import denoise_tv_chambolle fig, ax = plt.subplots(1, 2, figsize=(20, 10)) # zooming in for a better view _ = ds_single_channel.pp[1500:2000, 1500:2000].pl.colorize("gold").pl.show(ax=ax[0]) _ = ( ds_single_channel.pp.apply(func=denoise_tv_chambolle, weight=0.001) .pp[1500:2000, 1500:2000] .pl.colorize("gold") .pl.show(ax=ax[1]) ) ax[0].set_title("Raw") ax[1].set_title("Total variation filter") plt.show() ``` We can of course also combine various processing methods. Let's have a look at some more. ```python import skimage import numpy as np # example for local thresholding def local_threshold_single_channel(im, threshold=3): selem = skimage.morphology.disk(100) # mean filter im_mean = skimage.filters.rank.mean(im, footprint=selem) # threshold based on mean filter out = np.array(im > threshold * im_mean, dtype=np.uint8) return out fig, ax = plt.subplots(1, 2, figsize=(20, 10)) # zooming in for a better view _ = ds_single_channel.pp[1500:2000, 1500:2000].pl.colorize("gold").pl.show(ax=ax[0]) _ = ( ds_single_channel.pp.apply(func=local_threshold_single_channel) .pp[1500:2000, 1500:2000] .pl.colorize("gold") .pl.show(ax=ax[1]) ) ax[0].set_title("Raw") ax[1].set_title("Local Thresholding") plt.show() ``` We can also try to automatically find a good threshold. This follows the approach demonstrated [here](http://bebi103.caltech.edu.s3-website-us-east-1.amazonaws.com/2019a/content/lessons/lesson_05/basic_filtering_and_thresholding.html#Thresholding), where we look for the highest rate of change. ```python def find_optimal_threshold(im, min_k=0, max_k=5, step=100): selem = skimage.morphology.disk(100) arr = im[0, :, :] # mean filter arr_mean = skimage.filters.rank.mean(arr, footprint=selem) k = np.linspace(min_k, max_k, step) n_pix = np.empty_like(k) for i in range(len(k)): n_pix[i] = (arr > k[i] * arr_mean).sum() plt.plot(k, n_pix) plt.show() # Compute rough second derivative dn_pix_dk2 = np.diff(np.diff(n_pix)) # Find index of maximal second derivative max_ind = np.argmax(dn_pix_dk2) # Use this index to set k k_opt = k[max_ind - 2] return k_opt k_opt = find_optimal_threshold(ds_single_channel.pp[1500:2000, 1500:2000]["_image"]) print(f"Optimal threshold: {k_opt}") ``` Output: ```text Optimal threshold: 0.8585858585858586 ``` ```python fig, ax = plt.subplots(1, 2, figsize=(20, 10)) _ = ds_single_channel.pp[1500:2000, 1500:2000].pl.colorize("gold").pl.show(ax=ax[0]) _ = ( ds_single_channel.pp.apply(func=local_threshold_single_channel, threshold=k_opt) .pp[1500:2000, 1500:2000] .pl.colorize("gold") .pl.show(ax=ax[1]) ) ax[0].set_title("Raw") ax[1].set_title("Local Thresholding") ``` Output: ```text Text(0.5, 1.0, 'Local Thresholding') ``` In this particular instance, the automated finding of a cutoff does not correspond to the best one, and manually looking at the parameters appears more promising. However, the example above illustrates one way to automate thresholding, which may work in certain situations. --- # Cell Type Prediction The goal of cell type prediction is to assign a label to each cell. While methods such as clustering are frequently used in single cell transcriptomics data, these methods do not necessarily translate well to highly multiplexed fluorescence images. Instead, this notebook follows a two-step approach. In the first step, major cell types (such as B cells, T cells, macrophages, etc.) are predicted by using a set of mutually exclusive markers. These initial labels are predicted with methods such as `argmax` or `astir` and build the foundation for step two, in which we will look at functional markers. Here, we use a binarization-like approach. This means that for our functional markers, each cell gets either a value of 1 (if the cell expresses that marker) or 0 (if the cell does not express that marker). While this is a simplification of the underlying biology, this abstraction can be useful to counteract the often poor signal-to-noise ratio inherent in the data. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics as sp import pandas as pd import matplotlib.pyplot as plt import xarray as xr from scipy.signal import medfilt2d ``` ```python ds = xr.open_zarr("../../data/LN_24_1_unprocessed.zarr") ``` ## Image Processing Before predicting cell types, it is typically helpful to perform some preprocessing on the images. Here, we simply perform some custom thresholding on the markers to clean up the images and boost the signal-to-noise ratio. ```python # preprocessing: thresholding by percentiles and applying a median filter channels = [ "PAX5", "CD3", "CD11b", "CD11c", "CD15", "CD68", "Podoplanin", "CD31", "CD34", "CD90", "CD56", "CD4", "CD8", "Ki-67", "CD45RA", "CD45RO", ] quantiles = [0.8, 0.5, 0.8, 0.8, 0.8, 0.8, 0.95, 0.95, 0.95, 0.95, 0.8, 0.95, 0.95, 0.95, 0.95, 0.95] colors = [ "#e6194B", "#3cb44b", "#ffe119", "#4363d8", "#ffd8b1", "#f58231", "#911eb4", "#fffac8", "#469990", "#fabed4", "#9A6324", ] channels_celltypes = channels[:-5] # processing and recomputing quantification ds_processed = ( ds.pp[channels].pp.threshold(quantiles).pp.add_quantification().pp.transform_expression_matrix(method="arcsinh") ) ``` ```python # plotting the ds and ds processed next to one another fig, ax = plt.subplots(1, 2, figsize=(10, 5)) _ = ds.pp[channels_celltypes].pl.colorize(colors).pl.show(ax=ax[0]) _ = ds_processed.pp[channels_celltypes].pl.colorize(colors).pl.show(ax=ax[1]) ax[0].set_title("Raw") ax[1].set_title("Processed") plt.show() ``` ## Cell Type Prediction with Argmax Given the thresholded image, the simplest way to derive cell type labels is to define associations between cell types and markers, and then assign the cell type whose marker is most highly expressed (i. e. taking the argmax). This technique, also termed `argmax`, is demonstrated below. ```python marker_ct_dict = { "PAX5": "B", "CD3": "T", "CD11b": "Myeloid", "CD11c": "Dendritic", "CD15": "Granulo", "CD68": "Macro", "Podoplanin": "Stromal PDPN", "CD31": "Stromal CD31", "CD34": "Stromal CD34", "CD90": "Stromal CD90", "CD56": "NK", } ``` ```python # predicting the cell type using the argmax ds_with_ct_predictions = ds_processed.la.predict_cell_types_argmax( marker_ct_dict, key="_intensity", overwrite_existing_labels=True ) # adding colors to match the markers ds_with_ct_predictions = ds_with_ct_predictions.la.set_label_colors(list(marker_ct_dict.values()), colors) ``` ```python # plotting the ct predictions next to the processed image fig, ax = plt.subplots(1, 2, figsize=(10, 5)) _ = ds_with_ct_predictions.pp[channels_celltypes].pl.colorize(colors).pl.show(ax=ax[0]) _ = ds_with_ct_predictions.pl.show(render_image=False, render_labels=True, ax=ax[1]) ax[0].set_title("Markers") ax[1].set_title("Cell Type Predictions") ``` Output: ```text Text(0.5, 1.0, 'Cell Type Predictions') ``` ```python fig, ax = plt.subplots(1, 3, figsize=(21, 7)) _ = ds_processed.pp[channels_celltypes].pl.colorize(colors).pl.autocrop(padding=100).pl.show(ax=ax[0]) _ = ds.pp["DAPI"].pl.autocrop(padding=100).pl.show(render_segmentation=True, ax=ax[1]) _ = ds_with_ct_predictions.pl.autocrop(padding=100).pl.show(render_image=False, render_labels=True, ax=ax[2]) # removing the x and y ticks ax[0].set_xticks([]) ax[0].set_yticks([]) ax[1].set_xticks([]) ax[1].set_yticks([]) ax[2].set_xticks([]) ax[2].set_yticks([]) # setting titles ax[0].set_title("Markers") ax[1].set_title("Segmentation") ax[2].set_title("Cell Types") plt.tight_layout() ``` Getting observations from the xarray object is very easy. We can simply use the method `pp.get_layer_as_df()`, which can also parse our cell types into human-readable format automatically. ```python # getting the obs with get_layer_as_df() to get human readable cell type annotations in a pandas data frame ds_with_ct_predictions.pp.get_layer_as_df("_obs", celltypes_to_str=True) ``` Output: ```text _labels centroid-0 centroid-1 1 Stromal CD34 343.043825 1492.916335 2 B 341.740506 1517.082278 3 B 342.250000 1450.492857 4 B 345.150602 1559.259036 5 T 346.454128 1575.885321 ... ... ... ... 18130 B 2662.788136 1437.889831 18131 B 2662.327160 1570.314815 18132 B 2663.756757 1628.040541 18133 B 2662.770833 1656.781250 18134 B 2665.926829 1535.097561 [18134 rows x 3 columns] ``` Depending on which markers you have chosen, it could happen that some cells have no or very low expression for those markers. In such cases, it can be desirable to set cells with no expression to be "Unlabeled". You can control this using the `min_intensity` argument in `la.predict_cell_types_argmax()`. By default this is set to 0, and increasing it will lead to less cells being labeled. The example below illustrates this. ```python # let's just consider two markers # not all cells in this excerpt are B or T cells, hence they will express neither of the two markers ds_tmp = ds_processed.pp[["PAX5", "CD3"], 1000:1200, 1000:1200] # just predicting the maximum ds_all = ds_tmp.la.predict_cell_types_argmax({"PAX5": "B", "CD3": "T"}).la.set_label_colors( ["B", "T"], ["red", "green"] ) # cells where neither marker is above min_intensity are not labeled ds_with_unlabeled = ds_tmp.la.predict_cell_types_argmax( {"PAX5": "B", "CD3": "T"}, min_intensity=1.5 ).la.set_label_colors(["B", "T"], ["red", "green"]) fig, ax = plt.subplots(1, 3, figsize=(15, 5)) _ = ds_tmp.pl.show(render_segmentation=True, ax=ax[0]) _ = ds_all.pl.show(render_image=False, render_segmentation=True, render_labels=True, ax=ax[1]) _ = ds_with_unlabeled.pl.show(render_image=False, render_segmentation=True, render_labels=True, ax=ax[2]) ``` ## Cell Type Prediction with Astir There also exist more involved methods of predicting cell types, such as astir. Here is a quick example on how one could apply astir to the previous image. ```python ds_processed ``` Output: ```text Size: 219MB Dimensions: (cells: 18134, features: 2, channels: 16, x: 3000, y: 3000) Coordinates: * cells (cells) int64 145kB 1 2 3 4 5 ... 18131 18132 18133 18134 * channels (channels) _segmentation (y, x) int64 72MB 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 _image (channels, y, x) uint8 144MB 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 _intensity (cells, channels) float64 2MB 0.0007968 0.1034 ... 0.0 ``` ```python # minor reformatting seg = ds_processed["_segmentation"].values ds_for_astir = ds_processed.pp.drop_layers("_obs").pp.add_segmentation(seg) ``` ```python ds_for_astir ``` Output: ```text Size: 219MB Dimensions: (channels: 16, x: 3000, y: 3000, cells: 18134, features: 2) Coordinates: * channels (channels) Size: 219MB Dimensions: (labels: 12, la_props: 2, channels: 16, x: 3000, y: 3000, cells: 18134, features: 3) Coordinates: * labels (labels) int64 96B 1 2 3 4 5 6 7 8 9 10 11 12 * la_props (la_props) # Neighborhood Analysis To find structures within our tissues, we can use the concept of cellular neighborhoods. The idea is that we start by defining neighborhoods for each cell, which consist of the relative abundances of the surrounding cell types. Afterwards, we can use clustering across multiple samples to define the neighborhoods. Finally, we can add these neighborhood labels back into the `spatialproteomics` object and use them to select specific regions of the tissue. Methods relating to neighborhoods are contained in the `nh` module. For convenience, you can use the `sp.ImageContainer` to perform this analysis with a single function call. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics as sp import matplotlib.pyplot as plt import xarray as xr import numpy as np import pandas as pd from sklearn.cluster import KMeans import seaborn as sns ``` ## Neighborhood Analysis with the ImageContainer As a first step, we open three datasets and store them inside of a dictionary. ```python sp_obj_1 = xr.open_zarr("../../data/LN_24_1.zarr") sp_obj_2 = xr.open_zarr("../../data/LN_11_1.zarr") sp_obj_3 = xr.open_zarr("../../data/LN_13_1.zarr") # === OPTIONAL === # setting the label level to the top level, so that we only consider major cell types (B, T, etc.) # this is done for demonstration purposes, you can just as well compute the neighborhood on more fine-grained cell types celltype_colors = { "B cell": "#5799d1", "T cell": "#ebc850", "Myeloid cell": "#de6866", "Dendritic cell": "#4cbcbd", "Macrophage": "#bb7cb4", "Stromal cell": "#62b346", "Endothelial cell": "#bf997d", } sp_obj_1 = ( sp_obj_1.pp.drop_layers(["_neighborhoods", "_nh_properties", "_adjacency_matrix"]) .la.set_label_level("labels_0") .la.set_label_colors(celltype_colors.keys(), celltype_colors.values()) ) sp_obj_2 = ( sp_obj_2.pp.drop_layers(["_neighborhoods", "_nh_properties", "_adjacency_matrix"]) .la.set_label_level("labels_0") .la.set_label_colors(celltype_colors.keys(), celltype_colors.values()) ) sp_obj_3 = ( sp_obj_3.pp.drop_layers(["_neighborhoods", "_nh_properties", "_adjacency_matrix"]) .la.set_label_level("labels_0") .la.set_label_colors(celltype_colors.keys(), celltype_colors.values()) ) ``` ```python sp_dict = {"1": sp_obj_1, "2": sp_obj_2, "3": sp_obj_3} ``` Next, we put this dictionary into an `ImageContainer` and run the method `compute_neighborhoods()`. This will compute the neighborhoods for each object, and also run k-means clustering over all samples. If you want to do this yourself, please refer to the section at the end of this notebook. ```python image_container = sp.ImageContainer(sp_dict) # this method returns a dict in the same format as we provided as input sp_dict = image_container.compute_neighborhoods(k=3) ``` Let's plot the cell type labels and the resulting neighborhoods. ```python fig, ax = plt.subplots(2, 3, figsize=(15, 10)) ax = ax.flatten() # setting custom colors to be more in line with the cell types neighborhoods, colors = [f"Neighborhood {x}" for x in np.arange(3)], ["#5799d1", "#99a3b5", "#ebc850"] # plotting labels _ = sp_dict["1"].pl.autocrop().pl.show(render_image=False, render_labels=True, ax=ax[0]) _ = sp_dict["2"].pl.autocrop().pl.show(render_image=False, render_labels=True, ax=ax[1]) _ = sp_dict["3"].pl.autocrop().pl.show(render_image=False, render_labels=True, ax=ax[2]) # plotting neighborhoods _ = ( sp_dict["1"] .pl.autocrop() .nh.set_neighborhood_colors(neighborhoods, colors) .pl.show(render_image=False, render_neighborhoods=True, ax=ax[3]) ) _ = ( sp_dict["2"] .pl.autocrop() .nh.set_neighborhood_colors(neighborhoods, colors) .pl.show(render_image=False, render_neighborhoods=True, ax=ax[4]) ) _ = ( sp_dict["3"] .pl.autocrop() .nh.set_neighborhood_colors(neighborhoods, colors) .pl.show(render_image=False, render_neighborhoods=True, ax=ax[5]) ) plt.tight_layout() ``` By using the function `get_neighborhood_composition()` from the image container, we can easily plot a heatmap showing the neighborhood composition. By default, this composition dataframe is standardized per cell type. ```python # obtaining a data frame containing the neighborhood compositions nh_composition = image_container.get_neighborhood_composition() # plotting the neighborhood composition as a seaborn clustermap cluster_map = sns.clustermap( nh_composition, cmap="coolwarm", cbar_pos=(1.02, 0.2, 0.03, 0.77), # Position the colorbar to the right of the heatmap dendrogram_ratio=(0.00001, 0), # Remove row dendrogram figsize=(9, 7), # Adjust the size annot=False, # No annotation center=0, # Setting the center of the colorbar to 0 ) # customizations of the plot cluster_map.ax_heatmap.set_ylabel("") plt.show() ``` We can also rename the neighborhoods for convenience. ```python neighborhood_names = { "Neighborhood 0": "B high", "Neighborhood 1": "Other", "Neighborhood 2": "T high", } sp_dict["1"] = sp_dict["1"].nh.set_neighborhood_name(neighborhood_names.keys(), neighborhood_names.values()) ``` ## Comparison of different neighborhood methods There are three main ways how a neighborhood can be computed: via a radius, a k-nearest neighborh graph, or a Delaunay triangulation. In the `radius` method, a circle with a certain radius is drawn around a cell, and the cell type frequencies are aggregated. With `knn`, the k nearest neighbors of a cell are counted as the neighborhood. Finally, `delaunay` creates a Delaunay triangulation from the cell centroids and counts cells as part of each other's neighborhood if they are connected in the resulting graph. Let's look at how these three methods perform on the example data from above. To illustrate the differences better, we also cluster with a higher k, resulting in more neighborhoods. Note that you could do all of this with the `ImageContainer` as well. ```python # putting the strategy from above into a single method def compute_and_plot_clusters(sp_dict, k=5, method="radius"): # computing the neighborhoods image_container = sp.ImageContainer(sp_dict) sp_dict = image_container.compute_neighborhoods(neighborhood_method=method) # ===== PLOTTING ===== fig, ax = plt.subplots(2, 3, figsize=(15, 10)) ax = ax.flatten() _ = sp_dict["1"].pl.autocrop().pl.show(render_image=False, render_neighborhoods=True, ax=ax[0]) _ = sp_dict["2"].pl.autocrop().pl.show(render_image=False, render_neighborhoods=True, ax=ax[1]) _ = sp_dict["3"].pl.autocrop().pl.show(render_image=False, render_neighborhoods=True, ax=ax[2]) _ = sp_dict["1"].pl.autocrop().pl.show(render_image=False, render_labels=True, ax=ax[3]) _ = sp_dict["2"].pl.autocrop().pl.show(render_image=False, render_labels=True, ax=ax[4]) _ = sp_dict["3"].pl.autocrop().pl.show(render_image=False, render_labels=True, ax=ax[5]) for axis in ax: axis.set_xticks([]) axis.set_yticks([]) plt.suptitle(method) plt.tight_layout() plt.show() ``` Note that the colors for the neighborhoods were chosen arbitrarily here. ```python tmp_dict = {"1": sp_obj_1, "2": sp_obj_2, "3": sp_obj_3} compute_and_plot_clusters(tmp_dict, k=5, method="radius") ``` ```python tmp_dict = {"1": sp_obj_1, "2": sp_obj_2, "3": sp_obj_3} compute_and_plot_clusters(tmp_dict, k=5, method="knn") ``` ```python tmp_dict = {"1": sp_obj_1, "2": sp_obj_2, "3": sp_obj_3} compute_and_plot_clusters(tmp_dict, k=5, method="delaunay") ``` As a small sidenote, let's quickly see how we can plot the network. ```python plt.figure(figsize=(7, 7)) _ = tmp_dict["1"].pp[1500:2000, 1500:2000].pl.scatter_labels(size=15, render_edges=True) ``` ## Computing network features Neighborhoods naturally lend themselves to being defined via networks. Each cell represents a node in the network, and two cells are connected if they are in the same neighborhood. For example, if we define neighborhoods with the distance-based `radius` method, we can connect two cells if their centroids are less than `radius` units apart. This means that we can now employ techniques from network analysis to further analyse structures within the tissue. This step can be performed simply by calling `ds.nh.add_neighborhood_obs()`. For more details on which features you can compute, refer to the documentation of the method. Let's compute a couple of metrics. In the following codeblock, we compute: 1. Node degree: how many cells are connected to a cell. This can give you a hint about how dense the tissue is at any region (given you used the radius method for constructing your neighborhoods). 2. Homophily: ratio of cells of the same cell type as the center cell within the neighborhood. 3. Inter-Label Connectivity: ratio of cells of a different cell type as the center cell within the neighborhood. 4. Diversity Index: Shannon's entropy. Measures the diversity of cells within a cells neighborhood. ```python # setting label colors ds = sp_dict["1"].nh.set_neighborhood_colors(["T high", "Other", "B high"], ["#ebc850", "#99a3b5", "#5799d1"]) # computing the neighborhood features ds = ds.nh.add_neighborhood_obs( features=["degree", "homophily", "inter_label_connectivity", "diversity_index"] ).pl.autocrop() ``` ```python fig, ax = plt.subplots(2, 3, figsize=(15, 10)) ax = ax.flatten() _ = ds.pl.show(render_image=False, render_labels=True, ax=ax[0]) _ = ds.pl.show(render_image=False, render_neighborhoods=True, ax=ax[1]) _ = ds.pl.render_obs(feature="degree").pl.imshow(ax=ax[2]) _ = ds.pl.render_obs(feature="homophily").pl.imshow(ax=ax[3]) _ = ds.pl.render_obs(feature="inter_label_connectivity").pl.imshow(ax=ax[4]) _ = ds.pl.render_obs(feature="diversity_index").pl.imshow(ax=ax[5]) titles = ["Labels", "Neighborhoods", "Node degree", "Homophily", "Inter-Label Connectivity", "Diversity Index"] for i, axis in enumerate(ax): axis.set_title(titles[i]) plt.tight_layout() ``` Note that we can also compute these on subsets of the sample. For example, we might want to select only the B and T cells within the B-rich neighborhood, and then compute the measures on this data alone. ```python # subselecting and recomputing the features ds = ( ds.nh["B high"] .la[["B cell", "T cell"]] .nh.add_neighborhood_obs( features=["degree", "closeness_centrality", "homophily", "inter_label_connectivity", "diversity_index"] ) ) ``` ```python fig, ax = plt.subplots(2, 3, figsize=(15, 10)) ax = ax.flatten() _ = ds.pl.show(render_image=False, render_labels=True, ax=ax[0]) _ = ds.pl.show(render_image=False, render_neighborhoods=True, ax=ax[1]) _ = ds.pl.render_obs(feature="degree").pl.imshow(ax=ax[2]) _ = ds.pl.render_obs(feature="homophily").pl.imshow(ax=ax[3]) _ = ds.pl.render_obs(feature="inter_label_connectivity").pl.imshow(ax=ax[4]) _ = ds.pl.render_obs(feature="diversity_index").pl.imshow(ax=ax[5]) titles = ["Labels", "Neighborhoods", "Node degree", "Homophily", "Inter-Label Connectivity", "Diversity Index"] for i, axis in enumerate(ax): axis.set_title(titles[i]) plt.tight_layout() ``` --- # Interoperability This notebook shows some way that you can import and export data from `spatialproteomics`. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics as sp import pandas as pd import xarray as xr import os import shutil import anndata ``` ## Exporting Data Once you are happy with your analysis, you will likely want to export the results. The easiest way to do this is by using the `zarr` format, but `csv`, `anndata`, and `spatialdata` are also supported. ```python # loading a test file which we will export later # notice how easy it is to load the file from a zarr using xarray ds = xr.open_zarr("../../data/LN_11_1.zarr") ``` ## Exporting to Zarr This is the easiest file format to work with. It allows you to store and load the xarray objects with a single line of code. It is highly recommended to call `drop_encoding()` before exporting to zarr. There are several open issues linked to encoding problems, and this is the easiest way to circumvent them. For more references, refer to these issues: [issue 1](https://github.com/pydata/xarray/issues/3476), [issue 2](https://github.com/pydata/xarray/issues/9037). ```python zarr_path = "tmp.zarr" # removing the zarr if it exists if os.path.exists(zarr_path): shutil.rmtree(zarr_path) # exporting as zarr ds.drop_encoding().to_zarr("tmp.zarr") ``` Output: ```text ``` ## Exporting Tables to CSV Let's say you want to export some tables as csvs. This can be done with pandas. ```python df = ds.pp.get_layer_as_df("_obs") df.head() ``` Output: ```text BCL-2_binarized BCL-6_binarized CD163_binarized CD45RA_binarized \ 1 0.0 0.0 0.0 0.0 2 0.0 0.0 0.0 0.0 3 0.0 0.0 0.0 0.0 4 0.0 0.0 0.0 0.0 5 0.0 0.0 0.0 1.0 CD45RO_binarized CD8_binarized FOXP3_binarized Helios_binarized \ 1 0.0 0.0 0.0 0.0 2 0.0 0.0 0.0 0.0 3 0.0 0.0 0.0 0.0 4 0.0 0.0 0.0 0.0 5 0.0 0.0 0.0 0.0 ICOS_binarized Ki-67_binarized ... TCF7-TCF1_binarized Tim3_binarized \ ... ``` ```python # exporting as csv df.to_csv("tmp.csv") ``` ## Exporting to AnnData AnnData is a format used by scanpy, which can be useful to create interesting plots and downstream analyses. For this reason, you can export the xarray object as an AnnData object. Note that this object will only store the tabular data, but not the image or the segmentation layer. ```python # adding an additional quantification to demonstrate how those can be exported as well ds = ds.pp.add_quantification(func=sp.percentage_positive, key_added="_percentage_positive") ``` ```python # putting the expression matrix into an anndata object adata = ds.tl.convert_to_anndata( expression_matrix_key="_intensity", additional_layers={"percentage_positive": "_percentage_positive"}, additional_uns={"label_colors": "_la_properties"}, ) adata ``` Output: ```text AnnData object with n_obs × n_vars = 16871 × 56 obs: 'BCL-2_binarized', 'BCL-6_binarized', 'CD163_binarized', 'CD45RA_binarized', 'CD45RO_binarized', 'CD8_binarized', 'FOXP3_binarized', 'Helios_binarized', 'ICOS_binarized', 'Ki-67_binarized', 'PD-1_binarized', 'TCF7-TCF1_binarized', 'Tim3_binarized', '_labels', '_neighborhoods', 'centroid-0', 'centroid-1', 'degree', 'diversity_index', 'homophily', 'inter_label_connectivity' uns: '_labels_colors', 'label_colors' obsm: 'spatial' layers: 'percentage_positive' ``` ```python # writing to disk as hdf5 adata.write("tmp.h5ad") ``` ## Exporting to SpatialData SpatialData is a data format which is commonly used for spatial omics analysis and combines the power of zarr with anndata. You can export to this data format as well. ```python spatialdata_object = ds.tl.convert_to_spatialdata(expression_matrix_key="_intensity") spatialdata_object ``` Output: ```text SpatialData object ├── Images │ └── 'image': DataArray[cyx] (56, 3000, 3000) ├── Labels │ └── 'segmentation': DataArray[yx] (3000, 3000) └── Tables └── 'table': AnnData (16871, 56) with coordinate systems: ▸ 'global', with elements: image (Images), segmentation (Labels) ``` ```python # storing as zarr file spatialdata_object.write("tmp.zarr") ``` Output: ```text INFO  The Zarr backing store has been changed from None the new file path: tmp.zarr ``` ## Importing from Spatialdata In the example workflow, you have already seen how to read data from a tiff file. If you already have your data in `spatialdata` format, you can also read it in from there. Reading in the data like this will convert the data from `spatialdata` format to `xarray` format, so that you can use the `xarray` backend of `spatialproteomics`. The example data was obtained from [here](https://spatialdata.scverse.org/en/latest/api/datasets.html#spatialdata.datasets.raccoon). ```python ds = sp.read_from_spatialdata("../../data/spatialdata_example.zarr", image_key="raccoon") ds ``` Output: ```text Size: 9MB Dimensions: (channels: 3, y: 768, x: 1024, cells: 71, features: 2) Coordinates: * channels (channels) int64 24B 0 1 2 * y (y) float64 6kB 0.0 1.0 2.0 3.0 ... 764.0 765.0 766.0 767.0 * x (x) float64 8kB 0.0 1.0 2.0 ... 1.021e+03 1.022e+03 1.023e+03 * cells (cells) int64 568B 1 2 3 4 5 6 7 8 ... 65 66 67 68 69 70 71 * features (features) _segmentation (y, x) int64 6MB 1 1 1 1 1 1 1 1 ... 70 70 70 70 70 70 70 70 _obs (cells, features) float64 1kB 39.8 198.7 ... 736.5 890.5 ``` ## Exporting to Seurat `Seurat` is widely used in the R community to analyze spatial omics data. Here, we show how you can create a `Seurat` object from a `spatialproteomics` one. It is most sensible to store the data frame to disk first, and then read it in with R. For demonstration purposes, we run R directly in this notebook. ```python # setup and data loading import os import sys # setup for rpy2 to use the correct R binary env_r_bin = "/g/huber/users/meyerben/.conda/envs/tmp_env_3/bin" os.environ["PATH"] = env_r_bin + os.pathsep + os.environ.get("PATH", "") if "R_HOME" in os.environ: del os.environ["R_HOME"] import xarray as xr import spatialproteomics as sp celltype_colors = { "B cell": "#5799d1", "T cell": "#ebc850", "Myeloid cell": "#de6866", "Dendritic cell": "#4cbcbd", "Macrophage": "#bb7cb4", "Stromal cell": "#62b346", "Endothelial cell": "#bf997d", } # loading in a data set ds = xr.open_zarr("../../data/LN_24_1.zarr") # for clearer visualizations, we set the cell types to the broadest level (B cells, T cells, ...) ds = ds.la.set_label_level("labels_0", ignore_neighborhoods=True).la.set_label_colors( celltype_colors.keys(), celltype_colors.values() ) ``` ```python # this gives us a data frame of cells by markers intensity_df = ds.pp.get_layer_as_df("_intensity") # getting the spatial information and cell types from the spatialproteomics object spatial_df = ds.pp.get_layer_as_df()[["centroid-0", "centroid-1", "_labels"]] # renaming the columns for better compatibility with R spatial_df.columns = ["centroid.0", "centroid.1", "labels"] ``` ```python # send the data to R, so that we can use it directly ``` ```python library(Seurat) library(dplyr) # intensity_df: each row = cell, each col = marker # transpose so markers become "genes" mat <- t(as.matrix(intensity_df)) # Make sure row/col names align colnames(mat) <- rownames(intensity_df) # Create Seurat object seu <- CreateSeuratObject(counts = mat) # Add metadata (spatial info + cell labels) seu <- AddMetaData(seu, metadata = spatial_df) ``` Output: ```text WARNING: The R package "reticulate" only fixed recently an issue that caused a segfault when used with rpy2: https://github.com/rstudio/reticulate/pull/1188 Make sure that you use a version of that package that includes the fix. Loading required package: SeuratObject Loading required package: sp ‘SeuratObject’ was built with package ‘Matrix’ 1.7.3 but the current version is 1.7.4; it is recomended that you reinstall ‘SeuratObject’ as the ABI for ‘Matrix’ may have changed Attaching package: ‘SeuratObject’ The following objects are masked from ‘package:base’: ... ``` ```python # Make sure the rownames exactly match the Seurat cells spatial_coords <- spatial_df[, c("centroid.0", "centroid.1")] # The cell names in the Seurat object cell_names <- Cells(seu) # Assign rownames rownames(spatial_coords) <- cell_names colnames(spatial_coords) <- c(1, 2) # Convert to matrix embeddings_mat <- as.matrix(spatial_coords) # Now create the spatial reduction seu[["spatial"]] <- CreateDimReducObject( embeddings = embeddings_mat, key = "spatial_", assay = DefaultAssay(seu) ) ``` ```python library(ggplot2) p <- DimPlot( seu, reduction = "spatial", group.by = "labels", pt.size = 1 ) + coord_fixed() # ensures x and y axes have equal scale print(p) ``` Output: ```text Learn more about the underlying theory at https://ggplot2-book.org/ ``` --- # Downstream Analysis After performing all of the preprocessing with `spatialproteomics`, you typically want to use other packages to extract further information from the data. Here, we show two examples, one using `squidpy` and the other one using the R package `spatstat`. `Spatialproteomics` provides several ways to export your data to facilitate interoperability with external tools. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import os import sys # setup for rpy2 to use the correct R binary env_r_bin = "/g/huber/users/meyerben/.conda/envs/tmp_env_3/bin" os.environ["PATH"] = env_r_bin + os.pathsep + os.environ.get("PATH", "") if "R_HOME" in os.environ: del os.environ["R_HOME"] import xarray as xr import spatialproteomics as sp import pandas as pd import squidpy as sq celltype_colors = { "B cell": "#5799d1", "T cell": "#ebc850", "Myeloid cell": "#de6866", "Dendritic cell": "#4cbcbd", "Macrophage": "#bb7cb4", "Stromal cell": "#62b346", "Endothelial cell": "#bf997d", } ``` ```python # loading in a data set ds = xr.open_zarr("../../data/LN_24_1.zarr") # for clearer visualizations, we set the cell types to the broadest level (B cells, T cells, ...) ds = ds.la.set_label_level("labels_0", ignore_neighborhoods=True).la.set_label_colors( celltype_colors.keys(), celltype_colors.values() ) ``` ## Neighborhood Enrichment with Squidpy We can use the `pp.convert_to_anndata()` function to turn the `spatialproteomics` object into an `anndata` object that is accepted by `squidpy`. We can then follow the [squidpy documentation](https://squidpy.readthedocs.io) to perform an enrichment analysis. ```python adata = ds.tl.convert_to_anndata() adata ``` Output: ```text AnnData object with n_obs × n_vars = 18134 × 56 obs: 'BCL-2_binarized', 'BCL-6_binarized', 'CD163_binarized', 'CD45RA_binarized', 'CD45RO_binarized', 'CD8_binarized', 'DAPI_binarized', 'FOXP3_binarized', 'Helios_binarized', 'ICOS_binarized', 'PD-1_binarized', 'TCF7-TCF1_binarized', 'Tim3_binarized', '_labels', '_neighborhoods', 'centroid-0', 'centroid-1', 'degree', 'diversity_index', 'homophily', 'inter_label_connectivity' uns: '_labels_colors' obsm: 'spatial' ``` ```python sq.gr.spatial_neighbors(adata, coord_type="generic", radius=140) sq.gr.nhood_enrichment(adata, cluster_key="_labels") sq.pl.nhood_enrichment(adata, cluster_key="_labels") ``` Output: ```text 0%| | 0/1000 [00:00 # Interactivity with napari *Spatialproteomics* has a variety of functions to ensure interoperability with existing tools. For example, we can export the *spatialproteomics* object into *spatialdata*, and the use *napari* for interactive visualizations. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics as sp import xarray as xr from napari_spatialdata import Interactive import numpy as np ``` ```python # loading in a data set and performing some formatting for convenience ds = xr.open_zarr("../../data/LN_11_1.zarr") ``` ```python # turning the spatialproteomics object into a spatialdata object sd_obj = ds.tl.convert_to_spatialdata() sd_obj ``` Output: ```text INFO  Transposing `data` of type: <class 'dask.array.core.Array'> to ('c', 'y', 'x'). INFO  Transposing `data` of type: <class 'dask.array.core.Array'> to ('y', 'x'). SpatialData object ├── Images │ └── 'image': DataArray[cyx] (56, 3000, 3000) ├── Labels │ └── 'segmentation': DataArray[yx] (3000, 3000) └── Tables └── 'table': AnnData (16871, 56) with coordinate systems: ▸ 'global', with elements: image (Images), segmentation (Labels) ``` ```python interactive = Interactive(sd_obj) interactive.run() ``` Running this code in a python script will open `napari`, which you can then use to interactively look at your data. You can look at different markers, visualize cell types, neighborhoods, and many more. Below are some examples of what you can visualize with `napari`. #### Markers ![Visualization of a marker with napari](../_static/img/napari_marker.png) #### Segmentation ![Visualization of a segmentation with napari](../_static/img/napari_segmentation.png) #### Cell Type Labels ![Visualization of a cell type labels with napari](../_static/img/napari_labels.png) #### Binarization ![Visualization of a marker binarization with napari](../_static/img/napari_binarization.png) #### Neighborhoods ![Visualization of neighborhoods with napari](../_static/img/napari_neighborhoods.png) #### Features ![Visualization of homophily with napari](../_static/img/napari_homophily.png) --- # Customizability `Spatialproteomics` provides wrappers around a variety of external tools, such as `cellpose` or `astir`. However, as technology progresses, new tools emerge. To accomodate for this, `spatialproteomics` is fully customizable, meaning that you can for example load your own segmentation masks into the object, load custom cell type annotations, etc. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics import matplotlib.pyplot as plt import numpy as np import xarray as xr import cv2 import pandas as pd # loading in a data set and only keeping the image ds = xr.open_zarr("../../data/LN_24_1_unprocessed.zarr").pp.drop_layers(keep="_image") ds ``` Output: ```text Size: 504MB Dimensions: (channels: 56, y: 3000, x: 3000) Coordinates: * channels (channels) ``` ## Custom Segmentations If you have a custom segmentation, you can easily add it to the object using `pp.add_segmentation()`. ```python # dummy array. Instead of this, you could run your own segmentation method here and obtain a numpy array like that # note that this array should consist of integers, and 0 should mark the background custom_segmentation = np.zeros((3000, 3000), dtype=int) custom_segmentation[500:1500, 500:1500] = 1 ds = ds.pp.add_segmentation(custom_segmentation) ds ``` Output: ```text Size: 576MB Dimensions: (channels: 56, y: 3000, x: 3000, cells: 1, features: 2) Coordinates: * channels (channels) _segmentation (y, x) int64 72MB 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 _obs (cells, features) float64 16B 999.5 999.5 ``` Note that our object now contains the `_segmentation` layer. It also computed the centroids of the cells and stored them in `_obs`. ## Custom Image Processing To apply custom image processing methods, you can make use of the `pp.apply()` method, into which you can pass any function that transforms an array. For example, let's consider one that blurs the image and subtracts it from the original. ```python def simple_background_subtraction(image, kernel_size=51): # ensure input is numpy array img = np.array(image) # convert dtype for OpenCV (expects uint8 or float32) if img.dtype != np.uint8 and img.dtype != np.float32: img = cv2.normalize(img, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8) # compute background with Gaussian blur background = cv2.GaussianBlur(img, (kernel_size, kernel_size), 0) # subtract background subtracted = cv2.subtract(img, background) # normalize result result = cv2.normalize(subtracted, None, 0, 255, cv2.NORM_MINMAX) return result ds_processed = ds.pp.apply(func=simple_background_subtraction) ``` ```python fig, ax = plt.subplots(1, 2, figsize=(10, 5)) _ = ds.pp["DAPI"].pl.show(ax=ax[0]) _ = ds_processed.pp["DAPI"].pl.show(ax=ax[1]) for axis in ax: axis.axis("off") plt.tight_layout() ``` ## Custom Protein Quantification Typically, the mean or median intensity of each marker are aggregated per cell to obtain an expression matrix. However, this is also fully customizable. You could for example look at how spatially distributed a marker is within each cell mask. To do this, you can define a custom function that takes a segmentation mask and an image as input. ```python # loading the dataset including a proper cell segmentation # also subsetting for easier visualization ds = xr.open_zarr("../../data/LN_24_1_unprocessed.zarr").pp[1000:1300, 1000:1300] ``` ```python # this is a very simple method that looks at the variance of a marker within a cell def spatial_variance(regionmask, intensity_image): intensities = intensity_image[regionmask] return np.var(intensities) ds = ds.pp.add_quantification(func=spatial_variance) ds ``` Output: ```text Size: 6MB Dimensions: (channels: 56, y: 301, x: 301, cells: 416, features: 2) Coordinates: * cells (cells) int64 3kB 4304 4311 4312 4318 ... 7269 7270 7271 7306 * channels (channels) _obs (cells, features) float64 7kB dask.array _segmentation (y, x) int64 725kB 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 _intensity (cells, channels) float64 186kB 265.9 0.0 ... 7.35 0.2526 ``` Note that the computed values are now stored in the `_intensity` layer. For a quick visualization, we can put them into `_obs`. We do this by using `pp.add_feature()`. ```python dapi_variance = ds.pp.get_layer_as_df("_intensity")["DAPI"].values ds = ds.pp.add_feature(feature_name="dapi_variance", feature_values=dapi_variance) ``` ```python # visualization fig, ax = plt.subplots(1, 2, figsize=(10, 5)) _ = ds.pp["DAPI"].pl.show(render_segmentation=True, ax=ax[0]) _ = ds.pl.render_obs("dapi_variance").pl.imshow(legend_obs=True, ax=ax[1]) for axis in ax: axis.axis("off") plt.tight_layout() ``` ## Custom Cell Type Prediction If you have custom labels per cell, you can easily load them into the `spatialproteomics` object with the `la.add_labels_from_dataframe()` method. The input dataframe simply needs a column identifying the cells and one identifying the labels. ```python num_cells = ds.sizes["cells"] cell_ids = ds["cells"].values # define some example cell types and randomly assign them # substitute this with your cell type prediction method of choice cell_types = ["B cell", "T cell", "Macrophage"] cell_labels = np.random.choice(cell_types, size=num_cells, replace=True) label_df = pd.DataFrame({"cell": cell_ids, "label": cell_labels}) label_df.head() ``` Output: ```text cell label 0 4304 B cell 1 4311 Macrophage 2 4312 B cell 3 4318 T cell 4 4327 T cell ``` ```python # adding the labels to the spatialproteomics object ds = ds.la.add_labels_from_dataframe(label_df) # visualization _ = ds.pl.show(render_image=False, render_labels=True) ``` ## Custom Cell Features Depending on your individual workflow, you might have some additional data for each cell. For example, you could train a graph neural network (GNN) to predict cell embeddings based on their neighborhoods. You can easily load such things into the `spatialproteomics` object using the `pp.add_obs_from_dataframe()` function. ```python # Simulated GNN embeddings in the form of a dataframe df = pd.DataFrame( { "Embedding1": np.random.randn(ds.sizes["cells"]), "Embedding2": np.random.randn(ds.sizes["cells"]), "Embedding3": np.random.randn(ds.sizes["cells"]), } ) df.head() ``` Output: ```text Embedding1 Embedding2 Embedding3 0 -0.834064 -0.218655 0.412385 1 0.275729 -1.343602 -1.060857 2 0.962022 -1.457948 -0.463753 3 -0.297840 -0.571965 -0.329516 4 0.004459 0.321857 -1.582663 ``` ```python ds = ds.pp.add_obs_from_dataframe(df) _ = ds.pl.render_obs("Embedding1").pl.imshow(legend_obs=True) ``` ## Spot Detection with Spotiflow `Spatialproteomics` provides a wrapper around `spotiflow` to count the number of spots per cell. By default, only the final number of spots is written to the object, but you can set `return_spots=True` to have a more in-depth look into the spot predictions. Alternatively, you can work with a `spatialdata` object directly, which supports storing the points alongside cells, images, and other associated data. ```python # predicting the number of spots using spotiflow ds = ds.tl.spotiflow(channel="CD62L") # quick look at the obs. This now contains a column called CD62L_spots ds.pp.get_layer_as_df("_obs").head() ``` Output: ```text CD62L_spots Embedding1 Embedding2 Embedding3 _labels \ 4304 0.0 -0.834064 -0.218655 0.412385 B cell 4311 0.0 0.275729 -1.343602 -1.060857 Macrophage 4312 1.0 0.962022 -1.457948 -0.463753 B cell 4318 0.0 -0.297840 -0.571965 -0.329516 T cell 4327 1.0 0.004459 0.321857 -1.582663 T cell centroid-0 centroid-1 dapi_variance 4304 1000.016216 1094.767568 265.855670 4311 1002.259036 1126.759036 567.466987 4312 1001.149660 1256.564626 495.403975 4318 1001.702703 1080.144144 183.464188 4327 1002.888889 1049.494444 451.865198 ``` ```python fig, ax = plt.subplots(1, 2, dpi=200) _ = ds.pp["CD62L"].pl.show(render_segmentation=True, ax=ax[0]) _ = ds.pl.render_obs("CD62L_spots").pl.imshow(legend_obs=True, ax=ax[1]) ``` --- # Subselecting Data You can select parts of your data based on channels, spatial coordinates, cell type labels, and neighborhoods. If you want to follow along with this tutorial, you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import spatialproteomics import pandas as pd import xarray as xr ``` One of the key features of `spatialproteomics` is the ability to slice our image data quickly and intuitively. We start by loading our _spatialproteomics_ object. ```python ds = xr.open_zarr("../../data/LN_24_1.zarr") ``` ## Slicing Channels and Spatial Coordinates To slice specific channels of the image we simply use `.pp` accessor together with the familiar bracket `[]` indexing. ```python _ = ds.pp["CD4"].pl.show() ``` We can also select multiple channels by simply passing a list to the `.pp` accessor. As we will see later, this makes visualising image overlays easy. ```python _ = ds.pp[["CD4", "CD8"]].pl.show() ``` The `.pp` accessor also understands `x` and `y` coordinates. When `x` and `y` coordinates are sliced, we get rid of all cells that do not belong to the respective image slice. ```python _ = ds.pp[500:1500, 500:1500].pp["PAX5"].pl.show() ``` Note that we can also pass `channels` and `x, y` coordinates at the same time. ```python _ = ds.pp[["CD4", "CD8"], 500:1500, 500:1500].pl.show() ``` ## Slicing Labels The labels accessor `.la` allows to select specific cell types by their label number or name. ```python _ = ds.la[1].pl.show(render_image=False, render_labels=True) ``` ```python _ = ds.la["B cell"].pl.show(render_image=False, render_labels=True) ``` Again it is possible to pass multiple cell labels. ```python _ = ds.la[4, 5, 6].pl.show(render_image=False, render_labels=True) ``` Finally, we can select all cells except a cell type using `la.deselect`. ```python _ = ds.la.deselect(["B cell"]).pl.show(render_image=False, render_labels=True) ``` ## Slicing Neighborhoods We can also select by neighborhoods with the `nh` accessor. The syntax is identical to the one in the label subsetting. ```python # subsetting only neighborhood 1 _ = ds.nh[1].pl.render_neighborhoods().pl.imshow() ``` --- # Interactive Removal of Artifacts Sometimes, it can happen that only a part of the tissue is usable for analysis, for example due to artifacts or poor tissue quality in a certain region. In such cases, it can make sense to mask out problematic regions. This way, the image can still be used for downstream analysis, but the artifact will not influence results in a negative way. `Spatialproteomics` allows users to upload a binary mask, which indicates which part of the image is suitable for further analysis and which part is not. The example below shows how one could for example use [TissueTag](https://github.com/Teichlab/TissueTag/tree/main) to obtain such a mask, load it into the spatialproteomics and use the mask to remove the artifact. Note that TissueTag runs in jupyter notebooks or jupyter lab, but has no implementation for jupyterhub at the moment. A link to this exact sample can be found in the spatialproteomics manuscript, if you want to follow along with a different sample you can download the data [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). ```python import xarray as xr import spatialproteomics as sp import matplotlib.pyplot as plt import numpy as np from PIL import Image ``` ### Loading the image from a zarr file ```python ds = xr.open_zarr("4_B4.zarr") # quick plotting of some channels plt.figure(figsize=(10, 10)) ds.pp[["PAX5", "CD3", "Podoplanin"]].pl.autocrop().pl.show() ``` Output: ```text Dimensions: (cells: 8395, celltype_levels: 3, channels: 3, y: 1811, x: 1721, features: 9, labels: 9, props: 2, rgba: 4) Coordinates: * cells (cells) int64 1 2 3 4 5 ... 8391 8392 8393 8394 8395 * celltype_levels (celltype_levels) ``` ```python labels = tt.update_annotator(imarray=rgba_image, labels=labels, anno_dict=annodict, render_dict=render_dict) ``` Output: ```text Artifact ``` ```python # inverting the labels so that 0 becomes 1 and 1 becomes 0 # the region where the artifact is should be 0, everything else should be 1 labels = 1 - labels ``` ### Putting the mask into the spatialproteomics object ```python # by default, this method adds a layer with the key '_mask' ds = ds.pp.add_layer(labels) ``` ```python # plotting the mask (by using the render_segmentation method, since the data structure is conceptually similar to the segmentation) _ = ds.pl.render_segmentation(layer_key="_mask", alpha=0.7).pl.imshow() ``` ### Using the mask to set every pixel intensity to 0 in that region ```python ds_masked_intensities = ds.pp.mask_region() _ = ds_masked_intensities.pp[["PAX5", "CD3", "Podoplanin"]].pl.show() ``` You can see that this set all of the marker intensities at 0 in the problematic region. However, if we plot the segmentation masks on top, you can see that we have not removed those. ```python _ = ds_masked_intensities.pp[["PAX5", "CD3", "Podoplanin"]].pl.show(render_segmentation=True) ``` Most of the time, what we want to do instead is to remove all cells that were segmented in the area. This way we still stick to the original data, but remove the problematic region from analysis. ### Using the mask to remove cells in that region ```python ds_masked_cells = ds.pp.mask_cells() _ = ds_masked_cells.pp[["PAX5", "CD3", "Podoplanin"]].pl.show(render_segmentation=True) ``` You can of course also combine the two. ```python _ = ds.pp.mask_region().pp.mask_cells().pp[["PAX5", "CD3", "Podoplanin"]].pl.show(render_segmentation=True) ``` --- # Integration of Multiple Technologies There are different technologies to obtain multi-channel images: CODEX, Imaging Mass Cytometry (IMC), MACSima imaging cyclic staining (MICS), MIBI, etc. This notebook explores how you can use `spatialproteomics` to jointly analyse multiple modalities. ## Data Origin The IMC data was downloaded from [here](https://figshare.com/articles/dataset/Raw_images_for_DLBCL_samples/20010146?file=35633237). It contains IMC images of DLBCL TMAs. MICS data was downloaded from [here](https://zenodo.org/records/10057717). It contains a tonsil sample. The CODEX data can be downloaded from [here](https://oc.embl.de/index.php/s/XzEa9po1tjiDzzJ). For this example, we consider three lymph node biopsies. ## Data Preparation For all three data modalities, we follow the same general workflow of segmentation, protein quantification, and cell type prediction. Note that we do not need to use exactly the same tools, so if a certain segmentation method works better on one dataset than on the other, we can just apply different ones. Our downstream analysis is then carried out on the cell type level. This means that our data does not even need to have the same markers or cell types. The important part is that our cell type annotations should be sensible. We can then compare cell type abundances, interactions etc. across modalities and datasets. ```python import xarray as xr import spatialproteomics as sp from glob import glob import tifffile import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns celltype_colors = { "B cell": "#5799d1", "T cell": "#ebc850", "Myeloid cell": "#de6866", "Dendritic cell": "#4cbcbd", "Macrophage": "#bb7cb4", "Macrophage (M2)": "#bb7cb4", "Stromal cell": "#62b346", "Endothelial cell": "#bf997d", } ``` ### Imaging Mass Cytometry (IMC) For the IMC data, we have a diffuse large B cell lymphoma (DLBCL) sample. The data is comprised of individual tiff files (one per marker), as well as a tiff file containing a segmentation that was performed with `ilastik`. The first step is to read this data into `spatialproteomics`. ```python def read_imc_data(path): # === DATA READING === images = [] markers = [] for file in glob(path + "*.tiff"): if "ilastik" in file: continue marker_name = file.split("/")[-1].replace(".tiff", "") images.append(tifffile.imread(file)) markers.append(marker_name) image = np.stack(images) segmentation = tifffile.imread(path + "DLBCL TMA #5_1_3964_1_a0_ilastik_s2_Probabilities_mask.tiff") # === CONVERTING TO SPATIALPROTEOMICS === ds_imc = sp.load_image_data(image, channel_coords=markers, segmentation=segmentation) return ds_imc ds_imc = read_imc_data("../../data/raw-ablated-files/3964_DLBCL TMA #5_1_3964_1/") ``` Now that we have a dataset, we can look at the markers. ```python ds_imc["channels"].values ``` Output: ```text array(['Histone3.Yb176Di', 'DNA.Ir193Di', 'CD45RA.Gd155Di', 'Vimentin.Nd143Di', 'Ki67.Er168Di', 'Membrane.In115Di', 'Tim3.Sm154Di', 'ICOS.Nd148Di', 'PD1.NAT105..Lu175Di', 'EphrinB2.Er166Di', 'CD20.Dy161Di', 'BCL2.Nd146Di', 'BCL6.Sm147Di', 'FOXP3.Dy163Di', 'Vista.Dy160Di', 'Lag3.Eu153Di', 'C.Myc.Dy164Di', 'CD3.Er170Di', 'p.Stat3.Gd158Di', 'CD134.Eu151Di', 'CCR4.Sm149Di', 'CD206.Tm169Di', 'DNA2.Ir191Di', 'CXCR3.Nd142Di', 'CD68.Tb159Di', 'Granzym.B.Er167Di', 'CD45RO.Yb173Di', 'CD8.Dy162Di', 'PDL1.Nd150Di', 'CD4.Gd156Di', 'CD30.Ho165Di', 'Tbet.Nd145Di', 'CD31.Nd144Di', 'HLADR.Yb174Di', 'PDL2.Yb172Di'], dtype=' # FAQ Here are some common questions and answers. If you can't find what you're looking for here, please file an issue on our [GitHub page](https://github.com/sagar87/spatialproteomics). ## How can I optimize my memory usage? Running out of memory is a common problem when dealing with large images. Here are a couple of things you could consider to make your workflow more memory efficient. - **Lazy loading with zarr:** if you already have your objects stored in `zarr` format, it is usually more memory-efficient to load them using `xr.open_zarr(file)` than `xr.load_dataset(file)`. This is especially effective when you are subsetting your data, e. g. only looking at 10 channels when plotting a sample. - **Subsetting:** if you already know that you will require only a subset of your data, e. g. looking at certain channels, it is advised to perform subsetting as early as possible. This can be done with `ds.pp[channels]`. - **Deleting objects which are not required anymore:** spatialproteomics deliberately does not perform in-place operations, but rather copies the existing object to return a new one. This can be heavy on memory, if you do not remove intermediate variables once you do not need them anymore. You could consider removing them with `del` like this: ``` ds_new = ds.pp.your_workflow() del ds ``` - **Downsampling:** when looking at large images, you can downsample the image before plotting using `ds.pp.downsample(rate=8)`. When zooming into a specific area, you can omit the downsampling again. - **Garbage collection:** this is especially relevant if you perform operations in a for loop. If you try to store the dataset in the same variable, python's garbage collection might have some troubles freeing up memory due to cyclical references (for more information, please refer to [this post](https://dev.to/karishmashukla/how-python-uses-garbage-collection-for-efficient-memory-management-270h)). In this case, calling the garbage collector manually can help alleviate some of those issues: ``` import gc for ds in [...]: ds = ds.pp.your_workflow() gc.collect() # manually calling the garbage collector after each iteration ``` ## I have a lot of samples. Can I use parallelization to speed up the analysis? Yes, this is exactly what `dask` is useful for. Their [documentation](https://examples.dask.org/applications/embarrassingly-parallel.html) provides a good starting point for this. Since `spatialproteomics` is built on top of `xarray`, parallelization with `dask` is natively supported. If you want to apply the same workflow to a bunch of files, you could do it like this: ``` from dask import delayed, compute def process_file(path): import spatialproteomics as sp ds = (xr.open_zarr(path) .pp.add_quantification(func="intensity_mean") .pp.transform_expression_matrix(method="arcsinh")) return ds tasks = [delayed(process_file)(f) for f in files] results = compute(*tasks, scheduler="processes", num_workers=16) ``` This can result in substantial speedup compared to sequential execution with for loops, as illustrated below. *[Figure: Dask enables the parallelization of workflows, substantially speeding up execution times.]* ## When do I apply method directly on the object, and when do use sp.method()? Spatialproteomics has two distinct backends: a `xarray` backend and a `spatialdata` backend. These follow slightly different philosophies. The `xarray` backend is based on a functional programming design. This means that you can use it to call methods directly on your object, allowing you to pipe data from one step to the next. For example, this could look like `my_data.pp.segment().la.predict_cell_types().pl.show()`. Internally, `spatialproteomics` takes care of synchronizing shared dimensions across your data. The `spatialdata` backend is for when you want to use `spatialdata` objects from the start. In this case, the syntax is slightly different, but more similar to syntax you might be used to from scverse packages (such as scanpy or squidpy). Here, your code would look like this: ``` import spatialproteomics as sp sp.pp.segment(my_data) sp.pp.predict_cell_types(my_data) ``` These operations modify your object in-place, unless you set `copy=True` in the method signature. --- # The preprocessing (`pp`) accessor The preprocessing accessor provides several methods to subset and process image data. ## `ds.pp` (PreprocessingAccessor) The image accessor enables fast indexing and preprocessing of the spatialproteomics object. Available on any `xarray.Dataset` after `import spatialproteomics as sp`, as `ds.pp`. ### `ds.pp[...]` Fast subsetting the image container. The following examples show how the user can subset the image container. Slices of x and y coordinates are label-based (like xarray's ``.sel``), i.e. both start and stop are inclusive: ``ds.pp[0:50, 0:50]`` returns the pixels with coordinates 0 to 50 (51 x 51 pixels for integer coordinates). Cells are kept if their centroid lies within these bounds. Subset the image container using x and y coordinates: >>> ds.pp[0:50, 0:50] Subset the image container using x and y coordinates and channels: >>> ds.pp['Hoechst', 0:50, 0:50] Subset the image container using channels: >>> ds.pp['Hoechst'] Multiple channels can be selected by passing a list of channels: >>> ds.pp[['Hoechst', 'CD4']] **Parameters** indices : str, slice, list, tuple The indices to subset the image container. **Returns** xr.Dataset The subsetted image container. ### `ds.pp.get_bbox(x_slice: slice, y_slice: slice) -> xr.Dataset` Returns the bounds of the image container. **Parameters** x_slice : slice The slice representing the x-coordinates for the bounding box. Start and stop are inclusive. y_slice : slice The slice representing the y-coordinates for the bounding box. Start and stop are inclusive. **Returns** xr.Dataset The updated image container. ### `ds.pp.get_channels(channels: Union[List[str], str]) -> xr.Dataset` Retrieve the specified channels from the dataset. **Parameters** channels : Union[List[str], str] The channels to retrieve. Can be a single channel name or a list of channel names. **Returns** xr.Dataset The dataset containing the specified channels. ### `ds.pp.add_channel(channels: Union[str, list], array: np.ndarray, layer_key: str='_image') -> xr.Dataset` Adds channel(s) to an existing spatialproteomics object. **Parameters** channels : Union[str, list] The name of the channel or a list of channel names to be added. array : np.ndarray The numpy array representing the channel(s) to be added. layer_key : str The layer key where the channel(s) should be added. Default is '_image'. **Returns** xr.Dataset The updated spatialproteomics object with added channel(s). ### `ds.pp.drop_channel(channels: Union[str, list]) -> xr.Dataset` Drops channel(s) from an existing spatialproteomics object. **Parameters** channels : Union[str, list] The name of the channel or a list of channel names to be dropped. **Returns** xr.Dataset The updated spatialproteomics object with dropped channel(s). ### `ds.pp.add_segmentation(segmentation: Union[str, np.ndarray]=None, reindex: bool=False, keep_labels: bool=True, add_obs: bool=True) -> xr.Dataset` Adds a segmentation mask field to the xarray dataset. This will be stored in the '_segmentation' layer. **Parameters** segmentation : str or np.ndarray A segmentation mask, i.e., a np.ndarray with image.shape = (x, y), that indicates the location of each cell, or a layer key. mask_growth : int The number of pixels by which the segmentation mask should be grown. reindex : bool If true the segmentation mask is relabeled to have continuous numbers from 1 to n. keep_labels : bool When using cellpose on multiple channels, you may already get some initial celltype annotations from those. If you want to keep those annotations, set this to True. Default is True. add_obs : bool If True, centroids are added to the xarray. Default is True. **Returns** xr.Dataset The amended xarray. ### `ds.pp.add_layer(array: np.ndarray, key_added: str='_mask') -> xr.Dataset` Adds a layer (such as a mask highlighting artifacts) to the xarray dataset. **Parameters** array : np.ndarray The array representing the layer to be added. Can either be 2D or 3D (in this case, the first dimension should be the number of channels). key_added : str, optional The name of the added layer in the xarray dataset. Default is '_mask'. **Returns** xr.Dataset The updated dataset with the added layer. **Raises** AssertionError If the array is not 2-dimensional or its shape does not match the image shape. **Notes** This method adds a layer to the xarray dataset, where the layer has the same shape as the image field. The array should be a 2-dimensional numpy array representing the segmentation mask or layer to be added. The layer is created as a DataArray with the same coordinates and dimensions as the image field. The name of the added layer in the xarray dataset can be specified using the `key_added` parameter. The amended xarray dataset is returned after merging the original dataset with the new layer. ### `ds.pp.add_layer_from_dataframe(df: pd.DataFrame, key_added: str='_la_layers') -> xr.Dataset` Adds a dataframe as a layer to the xarray object. This is similar to add_obs, with the only difference that it can be used to add any kind of data to the xarray object. Useful to add things like string-based labels or other metadata. **Parameters** df : pd.DataFrame A dataframe with the observation values. **Returns** xr.Dataset The amended image container. ### `ds.pp.add_observations(properties: Union[str, list, tuple]=('label', 'centroid'), layer_key: str='_segmentation', return_xarray: bool=False) -> xr.Dataset` Adds properties derived from the segmentation mask to the image container. **Parameters** properties : Union[str, list, tuple] A list of properties to be added to the image container. See skimage.measure.regionprops_table for a list of available properties. layer_key : str The key of the layer that contains the segmentation mask. return_xarray : bool If true, the function returns an xarray.DataArray with the properties instead of adding them to the image container. **Returns** xr.Dataset The amended image container. ### `ds.pp.drop_observations(properties: Union[str, list, tuple], key: str='features') -> xr.Dataset` ### `ds.pp.add_feature(feature_name: str, feature_values: Union[list, np.ndarray])` Adds a feature to the image container. **Parameters** feature_name : str The name of the feature to be added. feature_values : The values of the feature to be added. **Returns** xr.Dataset The updated image container with the added feature. ### `ds.pp.add_obs_from_dataframe(df: pd.DataFrame) -> xr.Dataset` Adds an observation table to the image container. Each column of the dataframe is added as a feature to the obs layer. If the index of the dataframe contains exactly the cell ids of the image container, the rows are matched to the cells by their id (in any order). Otherwise, the rows are assumed to be in the same order as the cell coordinates of the image container. **Parameters** df : pd.DataFrame A dataframe with the observation values (one row per cell). **Returns** xr.Dataset The amended image container. ### `ds.pp.add_quantification(func: Union[str, Callable]='intensity_mean', key_added: str='_intensity', layer_key: str='_image', return_xarray=False) -> xr.Dataset` Quantify channel intensities over the segmentation mask. **Parameters** func : Callable or str, optional The function used for quantification. Can either be a string to specify a function from skimage.measure.regionprops_table or a custom function. Default is 'intensity_mean'. key_added : str, optional The key under which the quantification data will be stored in the image container. Default is '_intensity'. layer_key : str, optional The key of the layer to be quantified. Default is '_image'. return_xarray : bool, optional If True, the function returns an xarray.DataArray with the quantification data instead of adding it to the image container. **Returns** xr.Dataset or xr.DataArray The updated image container with added quantification data or the quantification data as a separate xarray.DataArray. ### `ds.pp.add_quantification_from_dataframe(df: pd.DataFrame, key_added: str='_intensity') -> xr.Dataset` Adds an observation table to the image container. Columns of the dataframe have to match the channel coordinates of the image container, and the index of the dataframe has to match the cell coordinates of the image container. **Parameters** df : pd.DataFrame A dataframe with the quantification values. key_added : str, optional The key under which the quantification data will be added to the image container. **Returns** xr.Dataset The amended image container. ### `ds.pp.drop_layers(layers: Union[str, list, NoneType]=None, keep: Union[str, list, NoneType]=None, drop_obs: bool=True, suppress_warnings: bool=False) -> xr.Dataset` Drops layers from the image container. Can either drop all layers specified in layers or drop all layers but the ones specified in keep. **Parameters** layers : Union[str, list] The name of the layer or a list of layer names to be dropped. keep : Union[str, list] The name of the layer or a list of layer names to be kept. drop_obs : bool If True, the observations are removed when the label or neighborhood properties are dropped. Default is True. suppress_warnings : bool If True, warnings are suppressed. Default is False. **Returns** xr.Dataset The updated image container with dropped layers. ### `ds.pp.threshold(quantile: Union[float, list]=None, intensity: Union[int, list]=None, key_added: Optional[str]=None, channels: Union[str, list, NoneType]=None, shift: bool=True, **kwargs)` Apply thresholding to the image layer of the object. By default, shift is set to true. This means that the threshold value is subtracted from the image, and all negative values are set to 0. If you instead want to set all values below the threshold to 0 while retaining the rest of the image at the original values, set shift to False. **Parameters** quantile : float The quantile value used for thresholding. If provided, the pixels below this quantile will be set to 0. intensity : int The absolute intensity value used for thresholding. If provided, the pixels below this intensity will be set to 0. key_added : Optional[str]) The name of the new image layer after thresholding. If not provided, the original image layer will be replaced. channels : Optional[Union[str, list]]) The channels to apply the thresholding to. If None, the thresholding will be applied to all channels. shift : bool If True, the thresholded image will be shifted so that values do not start at an arbitrary value. Default is True. **Returns** xr.Dataset The object with the thresholding applied to the image layer. **Raises** ValueError If both quantile and intensity are None or if both quantile and intensity are provided. ### `ds.pp.apply(func: Callable, key: str='_image', key_added: str='_image', **kwargs)` Apply a function to each channel independently. **Parameters** func : Callable The function to apply to the layer. key : str The key of the layer to apply the function to. Default is '_image'. key_added : str The key under which the updated layer will be stored. Default is '_image' (i. e. the original image will be overwritten). **kwargs : dict, optional Additional keyword arguments to pass to the function. **Returns** xr.Dataset The updated image container with the applied function. ### `ds.pp.normalize()` Performs a percentile normalization on each channel using the 3- and 99.8-percentile. Resulting values are in the range of 0 to 1. **Returns** xr.Dataset The image container with the normalized image stored in '_plot'. ### `ds.pp.downsample(rate: int)` Downsamples the entire dataset by selecting every `rate`-th element along the x and y dimensions. **Parameters** rate : int The downsampling rate. Only every `rate`-th pixel (or coordinate) is kept. **Returns** xr.Dataset The downsampled dataset with updated x and y coordinates. ### `ds.pp.rescale(scale: float)` Rescales the image and segmentation mask in the object by a given scale. The image is interpolated, while the segmentation mask is rescaled with nearest-neighbor interpolation so that the cell ids are preserved. The x and y coordinates keep their origin and are spaced so that they still cover the same region as before (analogous to `pp.downsample`). Cells that disappear when downscaling are removed from all layers, and observations that can be derived from the segmentation mask (e.g. centroids or area) are recomputed on the rescaled mask. **Parameters** scale : float The scale factor by which to rescale the image and segmentation mask. **Returns** xr.Dataset The rescaled dataset with updated image and segmentation mask. ### `ds.pp.filter_by_obs(col: str, func: Callable, segmentation_key: str='_segmentation', reindex: bool=False)` Filter the object by observations based on a given feature and filtering function. Parameters: col (str): The name of the feature to filter by. func (Callable): A filtering function that takes in the values of the feature and returns a boolean array. segmentation_key (str): The key of the segmentation mask in the object. Default is Layers.SEGMENTATION. reindex (bool): Whether to reindex the cell IDs after filtering. Default is False. Returns: xr.Dataset: The filtered object with the selected cells and updated segmentation mask. Raises: AssertionError: If the feature does not exist in the object's observations. Notes: - This method filters the object by selecting only the cells that satisfy the filtering condition. Example: To filter the object by the feature "area" and keep only the cells with an area greater than 70px: `obj = obj.pp.add_observations('area').pp.filter_by_obs('area', lambda x: x > 70)` ### `ds.pp.remove_outlying_cells(dilation_size: int=25, threshold: int=5, segmentation_key: str='_segmentation', reindex: bool=False)` Removes outlying cells from the image container. It does so by dilating the segmentation mask and removing cells that belong to a connected component with less than 'threshold' cells. **Parameters** dilation_size : int The size of the dilation kernel. Default is 25. threshold : int The minimum number of cells in a connected component required for the cells to be kept. Default is 5. segmentation_key : str The key of the segmentation mask in the object. Default is '_segmentation'. reindex : bool Whether to reindex the cell IDs after filtering. Default is False. **Returns** xr.Dataset The updated image container with the outlying cells removed. **Raises** ValueError If the object does not contain a segmentation mask. ### `ds.pp.grow_cells(iterations: int=2, suppress_warning: bool=False) -> xr.Dataset` Grows the segmentation masks by expanding the labels in the object. **Parameters** iterations : int The number of iterations to grow the segmentation masks. Default is 2. suppress_warning :bool Whether to suppress the warning about recalculating the observations. Used internally, default is False. **Raises** ValueError If the object does not contain a segmentation mask. **Returns** xr.Dataset The object with the grown segmentation masks and updated observations. ### `ds.pp.merge_segmentation(layer_key: str, key_added: str='_merged_segmentation', labels: Optional[List[str]]=None, threshold1: float=0.7, threshold2: float=0.8)` Merge segmentation masks. Merging works as follows: 1. Assume two segmentation masks, S1 (higher priority, e.g., macrophages) and S2 (lower priority, e.g., lymphocytes) 2. We start by copying the segmentation regions from S2 into S1, but only for pixels where S1 has background 3. Some S2 cells may now be partially overwritten by S1 cells — we drop any S2 cell where less than threshold1 fraction of its original area survived (e.g. with threshold1=0.7, an S2 cell that lost more than 30% of its area to S1 is discarded) 4. Of the remaining S2 cells, we additionally drop any that overlap with threshold2 or more of any single S1 cell's area — this prevents S2 cells from being retained when they largely encapsulate an S1 cell (e.g. with threshold2=0.8, an S2 cell is discarded if it covers more than 80% of any S1 cell) 5. The surviving S1 and S2 regions are combined into a final mask with sequential labels **Parameters** layer_key : Union[str, List[str]] The key(s) of the segmentation mask(s) to merge. Can be a single key (must be 3D) or a list of keys (each 2D). key_added : str The name of the new segmentation mask to be added to the xarray object. Default is "_merged_segmentation". labels : Optional[List[str]] Optional. Labels corresponding to each segmentation mask. If provided, must match number of arrays. By default, the channel names are used for a single 3D mask, and the layer keys for a list of 2D masks. threshold1: float Optional. Threshold for merging cells. Default is 0.7. threshold2: float Optional. Threshold for merging cells. Default is 0.8. **Returns** xr.Dataset The xarray object with the merged segmentation mask. **Raises** AssertionError If specified keys are not found or other input inconsistencies exist. **Notes** - If the input array is 2D, it will be expanded to 3D. - If labels are provided, they need to match the number of arrays. - The merging process starts with merging the biggest cells first, then the smaller ones. - Disconnected cells in the input are handled based on the specified method. ### `ds.pp.get_layer_as_df(layer: str='_obs', celltypes_to_str: bool=True, neighborhoods_to_str: bool=True, idx_to_str: bool=False) -> pd.DataFrame` Returns the specified layer as a pandas DataFrame. Parameters: layer (str): The name of the layer to retrieve. Defaults to Layers.OBS. celltypes_to_str (bool): Whether to convert celltype labels to strings. Defaults to True. neighborhoods_to_str (bool): Whether to convert neighborhood labels to strings. Defaults to True. idx_to_str (bool): Whether to convert the index to strings. Defaults to False. Returns: pandas.DataFrame: The layer data as a DataFrame. ### `ds.pp.get_disconnected_cell() -> int` Returns the first disconnected cell from the segmentation layer. Returns: np.ndarray: The first disconnected cell from the segmentation layer. ### `ds.pp.transform_expression_matrix(method: str='arcsinh', key: str='_intensity', key_added: str='_intensity', cofactor: float=5.0, min_percentile: float=1.0, max_percentile: float=99.0, **kwargs)` Transforms the expression matrix based on the specified mode. Parameters: method (str): The transformation method. Available options are "arcsinh", "zscore", "minmax", "double_zscore", and "clip". key (str): The key of the expression matrix in the object. key_added (str): The key to assign to the transformed matrix in the object. cofactor (float): The cofactor to use for the "arcsinh" transformation. min_percentile (float): The minimum percentile value to use for the "clip" transformation. max_percentile (float): The maximum percentile value to use for the "clip" transformation. Returns: xr.Dataset: The object with the transformed matrix added. Raises: ValueError: If an unknown transformation mode is specified. AssertionError: If no expression matrix is found at the specified layer. ### `ds.pp.mask_region(key: str='_mask', image_key='_image', key_added='_image') -> xr.Dataset` Mask a region in the image. Parameters: key (str): The key of the region to mask. image_key (str): The key of the image layer in the object. Default is Layers.IMAGE. key_added (str): The key to assign to the masked image in the object. Default is Layers.IMAGE, which overwrites the original image. Returns: xr.Dataset: The object with the masked region in the image. ### `ds.pp.mask_cells(mask_key: str='_mask', segmentation_key='_segmentation', reindex: bool=False) -> xr.Dataset` Mask cells in the segmentation mask. Parameters: mask_key (str): The key of the mask to use for masking. segmentation_key (str): The key of the segmentation mask in the object. Default is Layers.SEGMENTATION. reindex (bool): Whether to reindex the cell IDs after filtering. Default is False. Returns: xr.Dataset: The object with the masked cells in the segmentation mask. ### `ds.pp.convert_to_8bit(key: str='_image', key_added: str='_image')` Convert the image to 8-bit. Parameters: key (str): The key of the image layer in the object. Default is '_image'. key_added (str): The key to assign to the 8-bit image in the object. Default is '_image', which overwrites the original image. Returns: xr.Dataset: The object with the image converted to 8-bit. ### `ds.pp.merge_channels(channels: List[str], key_added: str='merged_channel', normalize: bool=True, method: Union[str, Callable]='max', layer_key: str='_image') -> xr.Dataset` Merge specified channels into a single channel by summing their values. **Parameters** channels : List[str] The list of channel names to merge. key_added : str The name of the new channel. normalize : bool Whether to normalize the images before merging. Default is True. method : Union[str, Callable] The method to use for merging. Can be "max", "sum", "mean", or a custom callable function. Default is "max". layer_key : str The key of the image layer in the object. Default is Layers.IMAGE. **Returns** xr.Dataset The object with the merged channel added to the layer given by layer_key. ### `ds.pp.find_bright_cells(channels: Union[str, list, NoneType]=None, frac_channels: float=0.5, quantile: float=0.95, key: str='_intensity', key_added: str='is_bright') -> xr.Dataset` Flags cells that are bright across multiple channels (e.g., due to a lot of antibodies binding to that cell). A cell is considered bright in a channel if its intensity exceeds the given quantile threshold for that channel. Cells that are bright in at least frac_channels of the considered channels are flagged. **Parameters** channels : Optional[Union[str, list]] The channels to consider. If None, all channels are used. frac_channels : float The fraction of channels in which a cell must be bright to be flagged. Default is 0.5. quantile : float The quantile threshold above which a cell is considered bright in a given channel. Computed per channel. Default is 0.95. key : str The key of the intensity layer to use. Default is Layers.INTENSITY. key_added : str The name of the feature added to obs. Default is 'is_bright'. **Returns** xr.Dataset The updated image container with the 'is_bright' flag added to obs. **Raises** ValueError If no intensity matrix is found at the specified layer. **Example** Flag cells that are bright in at least 50% of all channels at the 95th percentile: >>> ds = ds.pp.find_bright_cells() Flag cells that are bright in at least 2 out of 3 specific channels: >>> ds = ds.pp.find_bright_cells(channels=['CD3', 'CD8', 'CD20'], frac_channels=0.66) ## Functions for the spatialdata backend (`sp.pp`) These functions take a `spatialdata.SpatialData` object as first argument and modify it in place unless `copy=True`. ### `sp.pp.add_quantification(sdata, func: Union[str, Callable]='intensity_mean', key_added: str='table', image_key: str='image', segmentation_key: str='segmentation', layer_key: Optional[str]=None, data_key: Optional[str]=None, copy: bool=False, **kwargs)` This function computes the quantification of the image data based on the provided segmentation masks. It extracts the image data and segmentation masks from the spatialdata object, applies the quantification function, and adds the quantification results to the spatialdata object. The quantification results are stored in an AnnData object, which is added to the tables attribute of the spatialdata object. The quantification function can be specified as a string or a callable function. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the image data and segmentation masks. func (Union[str, Callable], optional): The quantification function to be applied. Defaults to "intensity_mean". Can be a string or a callable function. key_added (str, optional): The key under which the quantification results will be stored in the tables attribute of the spatialdata object. Defaults to table. image_key (str, optional): The key for the image data in the spatialdata object. Defaults to image. segmentation_key (str, optional): The key for the segmentation masks in the spatialdata object. Defaults to segmentation. layer_key (Optional[str], optional): The key for the quantification results in the AnnData object. If None, a new layer will be created. Defaults to None. data_key (Optional[str], optional): The key for the image data in the spatialdata object. If None, the image_key will be used. Defaults to None. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. ### `sp.pp.add_observations(sdata, properties: Union[str, list, tuple]=('label', 'centroid'), segmentation_key: str='segmentation', table_key: str='table', copy: bool=False, **kwargs)` This function computes the observations for each region in the segmentation masks. It extracts the segmentation masks from the spatialdata object, computes the region properties, and adds the observations to the AnnData object stored in the tables attribute of the spatialdata object. The observations are computed using the regionprops_table function from skimage.measure. The properties to be computed can be specified as a string or a list/tuple of strings. The default properties are "label" and "centroid", but other properties can be added as well. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the segmentation masks. properties (Union[str, list, tuple], optional): The properties to be computed for each region. Defaults to ("label", "centroid"). segmentation_key (str, optional): The key for the segmentation masks in the spatialdata object. Defaults to segmentation. table_key (str, optional): The key under which the AnnData object is stored in the tables attribute of the spatialdata object. Defaults to table. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. ### `sp.pp.apply(sdata, func: Callable, key_added: str='image', image_key: str='image', data_key: Optional[str]=None, copy: bool=False, **kwargs)` This function applies a given function to the image data in the spatialdata object. It extracts the image data from the spatialdata object, applies the function, and adds the processed image back to the spatialdata object. The processed image is stored in the images attribute of the spatialdata object. The function can be any callable function that takes an image as input and returns a processed image. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the image data. func (Callable): The function to be applied to the image data. It should take an image as input and return a processed image. image_key (str, optional): The key for the image data in the spatialdata object. Defaults to image. data_key (Optional[str], optional): The key for the image data in the spatialdata object. If None, the image_key will be used. Defaults to None. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. **kwargs: Additional keyword arguments to be passed to the function. ### `sp.pp.threshold(sdata, image_key: str='image', quantile: Union[float, list]=None, intensity: Union[int, list]=None, key_added: str='image', channels: Union[str, list, NoneType]=None, shift: bool=True, copy: bool=False, **kwargs)` This function applies a threshold to the image data in the spatialdata object. It extracts the image data from the spatialdata object, applies the thresholding function, and adds the processed image back to the spatialdata object. The processed image is stored in the images attribute of the spatialdata object. The thresholding function can be specified using the quantile or intensity parameters. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the image data. image_key (str, optional): The key for the image data in the spatialdata object. Defaults to image. quantile (Union[float, list], optional): The quantile value(s) to be used for thresholding. If None, the intensity parameter will be used. Defaults to None. intensity (Union[int, list], optional): The intensity value(s) to be used for thresholding. If None, the quantile parameter will be used. Defaults to None. key_added (str, optional): The key under which the processed image will be stored in the images attribute of the spatialdata object. Defaults to image. channels (Optional[Union[str, list]], optional): The channel(s) to be used for thresholding. If None, all channels will be used. Defaults to None. shift (bool, optional): Whether to shift the intensities towards 0 after thresholding. Defaults to True. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. **kwargs: Additional keyword arguments to be passed to the thresholding function. ### `sp.pp.transform_expression_matrix(sdata, method: str='arcsinh', table_key: str='table', cofactor: float=5.0, min_percentile: float=1.0, max_percentile: float=99.0, copy: bool=False, **kwargs)` This function applies a transformation to the expression matrix in the spatialdata object. It extracts the expression matrix from the spatialdata object, applies the transformation function, and adds the transformed expression matrix back to the spatialdata object. The transformed expression matrix is stored in the tables attribute of the spatialdata object. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the expression matrix. method (str, optional): The transformation method to be applied. Defaults to "arcsinh". table_key (str, optional): The key under which the expression matrix is stored in the tables attribute of the spatialdata object. Defaults to "table". cofactor (float, optional): The cofactor to be used for the transformation. Defaults to 5.0. min_percentile (float, optional): The minimum percentile to be used for the transformation. Defaults to 1.0. max_percentile (float, optional): The maximum percentile to be used for the transformation. Defaults to 99.0. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. **kwargs: Additional keyword arguments to be passed to the transformation function. ### `sp.pp.filter_by_obs(sdata, col: str, func: Callable, segmentation_key: str='segmentation', table_key: str='table', reindex: bool=False, copy: bool=False)` Filter the object by observations based on a given feature and filtering function. Parameters: sdata (spatialdata.SpatialData): The spatialdata object to filter. col (str): The name of the feature to filter by. func (Callable): A filtering function that takes in the values of the feature and returns a boolean array. segmentation_key (str): The key of the segmentation mask in the object. Default is SDLayers.SEGMENTATION. table_key (str): The key of the table in the object. Default is SDLayers.TABLE. reindex (bool): Whether to reindex the obs after filtering. Default is False. copy (bool): If True, a copy of the object is returned. Default is False. ### `sp.pp.grow_cells(sdata, iterations: int=2, segmentation_key: str='segmentation', table_key: str='table', suppress_warning: bool=False, copy: bool=False) -> xr.Dataset` Grows the segmentation masks by expanding the labels in the object. **Parameters** sdata : spatialdata.SpatialData The spatialdata object containing the segmentation masks. iterations : int The number of iterations to grow the segmentation masks. Default is 2. segmentation_key : str The key of the segmentation mask in the object. Default is segmentation. suppress_warning :bool Whether to suppress the warning about recalculating the observations. Used internally, default is False. copy : bool If True, a copy of the object is returned. Default is False. **Raises** ValueError If the object does not contain a segmentation mask. **Returns** xr.Dataset The object with the grown segmentation masks and updated observations. ### `sp.pp.merge_channels(sdata, channels: List[str], key_added: str, normalize: bool=True, method: str='sum', image_key: str='image', copy: bool=False)` Merges multiple channels into a single channel in the spatialdata object. **Parameters** sdata : spatialdata.SpatialData The spatialdata object containing the image data. channels : List[str] The list of channels to be merged. key_added : str The key under which the merged channel will be stored in the images attribute of the spatialdata object. normalize : bool Whether to normalize the channels before merging. Default is True. method : str The method to use for merging the channels. Options are "sum", "mean", and "max". Default is "sum". image_key : str The key for the image data in the spatialdata object. Default is SDLayers.IMAGE. copy : bool If True, a copy of the object is returned. Default is False. **Raises** AssertionError If key_added already exists in the object. KeyError If any of the specified channels do not exist in the image data. ValueError If less than two channels are provided to merge. If an unknown merging method is provided. **Returns** xr.Dataset The object with the merged channel added. --- # The label (`la`) accessor The label accessor provides several functions associated with cell type annotations and gating. ## `ds.la` (LabelAccessor) Adds functions for cell phenotyping. Available on any `xarray.Dataset` after `import spatialproteomics as sp`, as `ds.la`. ### `ds.la.deselect(indices)` Deselect specific label indices from the data object. This method deselects specific label indices from the data object, effectively removing them from the selection. The deselection can be performed using slices, lists, tuples, or individual integers. **Parameters** indices : slice, list, tuple, or int The label indices to be deselected. Can be a slice, list, tuple, or an individual integer. **Returns** xr.Dataset The updated data object with the deselected label indices. **Notes** - The function uses 'indices' to specify which labels to deselect. - 'indices' can be provided as slices, lists, tuples, or an integer. - The function then updates the data object to remove the deselected label indices. ### `ds.la.add_label_type(name: str, color: str='w') -> xr.Dataset` Add a new label type to the data object. This method adds a new label type with the specified 'name' and 'color' to the data object. The label type is used to identify and categorize cells in the segmentation mask. **Parameters** name : str The name of the new label type to be added. color : str, optional The color code to represent the new label type in the visualization. Default is "white" ("w"). **Returns** xr.Dataset The updated data object with the newly added label type. **Raises** ValueError If the segmentation mask or observation table is not found in the data object. If the provided 'name' already exists as a label type. **Notes** - The function checks for the existence of the segmentation mask and observation table in the data object. - It ensures that the 'name' of the new label type does not already exist in the label types. - The function then adds the new label type with the given 'name' and 'color' to the data object. ### `ds.la.remove_label_type(cell_type: Union[int, str, List[Union[int, str]]]) -> xr.Dataset` Remove specific cell type label(s) from the data object. This method removes specific cell type label(s) identified by 'cell_type' from the data object. The cell type label(s) are effectively removed, and their associated cells are assigned to the 'Unlabeled' category. **Parameters** cell_type : int, str or list of int or str The ID(s) or name(s) of the cell type label(s) to be removed. **Returns** xr.Dataset The updated data object with the specified cell type label(s) removed. **Raises** ValueError If the data object does not contain any cell type labels. If the specified 'cell_type' is not found among the existing cell type labels. **Notes** - The function first checks for the existence of cell type labels in the data object. - It then removes the specified 'cell_type' from the cell type labels, setting their cells to the 'Unlabeled' category. ### `ds.la.add_label_property(array: Union[np.ndarray, list], prop: str)` Add a label property for each unique cell type label. This method adds a property, specified by 'prop', for each unique cell type label in the data object. The property values are taken from the 'array' argument and assigned to each corresponding cell type label. **Parameters** array : numpy.ndarray or list An array or list containing property values to be assigned to each unique cell type label. prop : str The name of the property to be added to the cell type labels. **Returns** xr.Dataset The updated data object with the added label property. ### `ds.la.set_label_name(label, name)` Set the name of a specific cell type label. This method sets the 'name' of a specific cell type label identified by the 'label'. The 'label' can be either a label ID or the name of the cell type label. **Parameters** label : int or str The ID or name of the cell type label whose name will be updated. name : str The new name to be assigned to the specified cell type label. **Returns** None **Notes** - The function converts the 'label' from its name to the corresponding ID for internal processing. - It updates the name of the cell type label in the data object to the new 'name'. ### `ds.la.set_label_colors(labels: Union[str, List[str]], colors: Union[str, List[str]])` Set the color of a specific cell type label. This method sets the 'color' of a specific cell type label identified by the 'label'. The 'label' can be either a label ID or the name of the cell type label. **Parameters** label : int or str The ID or name of the cell type label whose color will be updated. color : any The new color to be assigned to the specified cell type label. **Returns** xr.Dataset **Notes** - The function converts the 'label' from its name to the corresponding ID for internal processing. - It updates the color of the cell type label in the data object to the new 'color'. ### `ds.la.predict_cell_types_argmax(marker_dict: dict, key: str='_intensity', overwrite_existing_labels: bool=False, cell_col: str='cell', label_col: str='label', min_intensity: float=0.0, ignore_neighborhoods: bool=False)` Predicts cell types based on the argmax classification of marker intensities. Parameters: marker_dict (dict): A dictionary mapping markers to cell types. Each marker should be associated to one specific cell type. key (str, optional): The key of the quantification layer to use for classification. Defaults to "_intensity". overwrite_existing_labels (bool, optional): Whether to overwrite existing labels. Defaults to False. cell_col (str, optional): The name of the column to store cell IDs in the output dataframe. Defaults to "cell". label_col (str, optional): The name of the column to store predicted cell types in the output dataframe. Defaults to "label". min_intensity (float, optional): The minimum intensity threshold for assigning cells to a type. Defaults to 0.0. ignore_neighborhoods (bool, optional): Whether to remove existing neighborhoods when predicting cell types. Defaults to False. Returns: xr.Dataset: A new spatialproteomics object with the predicted cell types added as labels. Raises: AssertionError: If the quantification layer with the specified key is not found. AssertionError: If any of the markers specified in the marker dictionary are not present in the quantification layer. ### `ds.la.threshold_labels(threshold_dict: dict, label: Optional[str]=None, layer_key: str='_intensity')` Binarise based on a threshold. If a label is specified, the binarization is only applied to this cell type. **Parameters** threshold_dict : dict A dictionary mapping channels to threshold values. label : str, optional The specified cell type for which the binarization is applied, by default None. layer_key : str, optional The key for the new binary feature layer, by default "_percentage_positive_intensity". **Returns** xr.Dataset A new dataset object with the binary features added. **Notes** - If a label is specified, the binarization is only applied to the cells of that specific cell type. - The binary feature is computed by comparing the intensity values of each channel to the threshold value. - The binary feature is added as a new layer to the dataset object. ### `ds.la.add_labels(labels: Optional[dict]=None) -> xr.Dataset` Add labels from a mapping (cell -> label) to the spatialproteomics object. **Parameters** labels : Union[dict, None] A dictionary containing cell labels as keys and corresponding labels as values. If None, a default labeling will be added. Default is None. **Returns** xr.Dataset The spatialproteomics object with added labels. Notes: ------ This method converts the input dictionary into a pandas DataFrame and then adds the labels to the object using the `la.add_labels_from_dataframe` method. ### `ds.la.add_labels_from_dataframe(df: Optional[pd.DataFrame]=None, cell_col: str='cell', label_col: str='label', colors: Optional[list]=None, names: Optional[list]=None, ignore_neighborhoods: bool=False) -> xr.Dataset` Adds labels to the image container. **Parameters** df : Union[pd.DataFrame, None], optional A dataframe with the cell and label information. If None, a default labeling will be applied. cell_col : str, optional The name of the column in the dataframe representing cell coordinates. Default is "cell". label_col : str, optional The name of the column in the dataframe representing cell labels. Default is "label". colors : Union[list, None], optional A list of colors corresponding to the cell labels. If None, random colors will be assigned. Default is None. names : Union[list, None], optional A list of names corresponding to the cell labels. If None, default names will be assigned. Default is None. ignore_neighborhoods : bool, optional If True, the function will ignore the neighborhoods in the object. Default is False. **Returns** xr.Dataset The updated image container with added labels. ### `ds.la.add_properties(array: Union[np.ndarray, list], prop: str='_labels', return_xarray: bool=False) -> xr.Dataset` Adds properties to the image container. In order to add properties to an already existing property layer, use the la.add_label_property() method. **Parameters** array : Union[np.ndarray, list] An array or list of properties to be added to the image container. prop : str, optional The name of the property. Default is Features.LABELS. return_xarray : bool, optional If True, the function returns an xarray.DataArray with the properties instead of adding them to the image container. **Returns** xr.Dataset or xr.DataArray The updated image container with added properties or the properties as a separate xarray.DataArray. ### `ds.la.predict_cell_subtypes(subtype_dict: Union[dict, str], overwrite_existing_labels: bool=True) -> xr.Dataset` Predict cell subtypes based on the binarized marker intensities. **Parameters** subtype_dict : dict A dictionary mapping cell subtypes to the binarized markers used for prediction. Instead of a dictionary, a path to a yaml file containing the subtype dictionary can be provided. overwrite_existing_labels : bool, optional If True, existing labels will be overwritten by the new, more granular cell type predictions. Default is True. **Returns** xr.Dataset The updated image container with the predicted cell subtypes. ### `ds.la.set_label_level(level: str, ignore_neighborhoods: bool=False) -> xr.Dataset` Set the label level to a specific level. **Parameters** level : str The name of the label level to set. ignore_neighborhoods : bool, optional If True, the function will ignore the neighborhoods in the object. Default is False. **Returns** xr.Dataset The updated image container with the label level set to the specified level. ## Functions for the spatialdata backend (`sp.la`) These functions take a `spatialdata.SpatialData` object as first argument and modify it in place unless `copy=True`. ### `sp.la.threshold_labels(sdata, threshold_dict: dict, key_added: str='binarization', table_key: str='table', layer_key: str='perc_pos', copy: bool=False)` Binarise based on a threshold. If a label is specified, the binarization is only applied to this cell type. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the expression matrix. threshold_dict (dict): A dictionary containing the threshold values for each channel. key_added (str, optional): The key under which the processed expression matrix will be stored in the obsm attribute of the spatialdata object. Defaults to "binarization". table_key (str, optional): The key under which the expression matrix is stored in the tables attribute of the spatialdata object. Defaults to "table". layer_key (str, optional): The key under which the expression matrix is stored in the layers attribute of the spatialdata object. Defaults to "perc_pos". copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. ### `sp.la.predict_cell_types_argmax(sdata, marker_dict: dict, min_intensity: float=0.0, table_key: str='table', label_key: str='celltype', copy: bool=False)` This function predicts cell types based on the expression matrix using the argmax method. It extracts the expression matrix from the spatialdata object, applies the argmax method, and adds the predicted cell types to the spatialdata object. The predicted cell types are stored in the obs attribute of the AnnData object in the tables attribute of the spatialdata object. The argmax method assigns the cell type with the highest expression value to each cell. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the expression matrix. marker_dict (dict): A dictionary containing the marker genes for each cell type. min_intensity (float, optional): The minimum intensity threshold for assigning cells to a type. Cells with all marker intensities below this threshold will be assigned to "Unlabeled". Defaults to 0.0. table_key (str, optional): The key under which the expression matrix is stored in the tables attribute of the spatialdata object. Defaults to "table". label_key (str, optional): The key under which the cell type labels are stored in the obs attribute of the spatialdata object. Defaults to "celltype". copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. ### `sp.la.predict_cell_subtypes(sdata, subtype_dict: Union[dict, str], table_key: str='table', layer_key: str='binarization', label_key: str='celltype', copy: bool=False)` This function predicts cell subtypes based on the expression matrix using a subtype dictionary. It extracts the expression matrix from the spatialdata object, applies the subtype prediction method, and adds the predicted cell subtypes to the spatialdata object. The predicted cell subtypes are stored in the obs attribute of the AnnData object in the tables attribute of the spatialdata object. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the expression matrix. subtype_dict (Union[dict, str]): A dictionary mapping cell subtypes to the binarized markers used for prediction. Instead of a dictionary, a path to a yaml file containing the subtype dictionary can be provided. table_key (str, optional): The key under which the expression matrix is stored in the tables attribute of the spatialdata object. Defaults to "table". layer_key (str, optional): The key under which the binarized expression matrix is stored in the obsm attribute of the spatialdata object. Defaults to "binarization". label_key (str, optional): The key under which the cell type labels are stored in the obs attribute of the spatialdata object. Defaults to "celltype". copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. --- # The neighborhood (`nh`) accessor The neighborhood accessor provides several functions associated with cellular neighborhoods. ## `ds.nh` (NeighborhoodAccessor) Adds functions for cellular neighborhoods. Available on any `xarray.Dataset` after `import spatialproteomics as sp`, as `ds.nh`. ### `ds.nh.deselect(indices)` Deselect specific neighborhood indices from the data object. This method deselects specific neighborhood indices from the data object, effectively removing them from the selection. The deselection can be performed using slices, lists, tuples, or individual integers. **Parameters** indices : slice, list, tuple, or int The neighborhood indices to be deselected. Can be a slice, list, tuple, or an individual integer. **Returns** xr.Dataset The updated data object with the deselected neighborhood indices. **Notes** - The function uses 'indices' to specify which neighborhoods to deselect. - 'indices' can be provided as slices, lists, tuples, or an integer. - The function then updates the data object to remove the deselected neighborhood indices. ### `ds.nh.add_properties(array: Union[np.ndarray, list], prop: str='_neighborhoods') -> xr.Dataset` Adds neighborhood properties to the image container. **Parameters** array : Union[np.ndarray, list] An array or list of properties to be added to the image container. prop : str, optional The name of the property. Default is Features.NEIGHBORHOODS. **Returns** xr.Dataset The updated image container with added properties. ### `ds.nh.add_neighborhoods_from_dataframe(df: pd.DataFrame, neighborhood_col: str='_neighborhoods', colors: Optional[list]=None, names: Optional[list]=None) -> xr.Dataset` Add neighborhoods to the dataset from a DataFrame. **Parameters** df : pd.DataFrame DataFrame containing neighborhood information. neighborhood_col : str, optional Column name in the DataFrame that contains neighborhood labels, by default '_neighborhoods'. colors : list or None, optional List of colors for the neighborhoods, by default None. If None, random colors will be assigned. names : list or None, optional List of names for the neighborhoods, by default None. If None, default names will be assigned. **Returns** xr.Dataset Updated dataset with neighborhood information added. **Raises** AssertionError If neighborhood properties are already present in the object. If the number of neighborhoods does not match the number of cells. If the specified column does not exist in the DataFrame. If neighborhoods contain negative values. If the length of colors does not match the number of unique neighborhoods. If the length of names does not match the number of unique neighborhoods. ### `ds.nh.set_neighborhood_colors(neighborhoods: Union[str, List[str]], colors: Union[str, List[str]], suppress_warnings: bool=False)` Set the color of a specific neighborhood. This method sets the 'color' of a specific neighborhood identified by the 'neighborhood'. The 'neighborhood' can be either a neighborhood ID or the name of the neighborhood. **Parameters** label : int or str The ID or name of the neighborhood whose color will be updated. color : any The new color to be assigned to the specified neighborhood. suppress_warnings : bool, optional Whether to suppress warnings. Default is False. **Returns** xr.Dataset **Notes** - The function converts the 'neighborhood' from its name to the corresponding ID for internal processing. - It updates the color of the neighborhood in the data object to the new 'color'. ### `ds.nh.set_neighborhood_name(neighborhoods: Union[int, str, List], names: Union[str, List])` Set the name of one or more neighborhoods. This method updates the 'name' of specific neighborhoods identified by the 'neighborhood'. The 'neighborhoods' can be a single neighborhood ID or name, or a list of IDs/names. Similarly, the 'names' parameter can be a single string or a list of strings. **Parameters** neighborhoods : int, str, or list The ID(s) or name(s) of the neighborhoods whose names will be updated. names : str or list The new name(s) to be assigned to the specified neighborhoods. **Returns** xr.Dataset **Notes** - When both parameters are lists, their lengths must match. - The function converts each 'neighborhood' from its name to the corresponding ID for internal processing. - It updates the name(s) of the neighborhood(s) in the data object to the new 'name(s)'. ### `ds.nh.compute_neighborhoods_radius(radius: int=100, include_center: bool=True, key_added: str='_neighborhoods', key_adjacency_matrix: str='_adjacency_matrix')` Compute the neighborhoods of each cell based on a specified radius. This method defines a radius around each cell and identifies all cells within this radius. It then examines the predicted cell types of these cells, including the center cell itself, and computes the frequency of each cell type. **Parameters** radius : int, optional The radius around each cell to define the neighborhood (in pixels). Default is 100. include_center : bool, optional Whether to include the center cell in the neighborhood. Default is True. key_added : str, optional The key under which the computed neighborhoods will be stored in the resulting DataArray. Default is '_neighborhoods'. key_adjacency_matrix : str, optional The key under which the computed adjacency matrix will be stored in the resulting DataArray. Default is '_adjacency_matrix'. **Returns** xarray.Dataset A merged xarray Dataset containing the original data and the computed neighborhoods. ### `ds.nh.compute_neighborhoods_knn(k=10, include_center: bool=True, key_added: str='_neighborhoods', key_adjacency_matrix: str='_adjacency_matrix')` Compute the neighborhoods of each cell based on k-nearest neighbors. This method identifies the k-nearest neighbors for each cell and computes the frequency of each cell type within these k neighbors. **Parameters** k : int, optional The number of nearest neighbors to consider, not counting the center cell. Default is 10. include_center : bool, optional Whether to add the center cell to its k nearest neighbors in the neighborhood. Default is True. key_added : str, optional The key under which the computed neighborhoods will be stored in the resulting DataArray. Default is '_neighborhoods'. key_adjacency_matrix : str, optional The key under which the computed adjacency matrix will be stored in the resulting DataArray. Default is '_adjacency_matrix'. **Returns** xarray.Dataset A merged xarray Dataset containing the original data and the computed neighborhoods. ### `ds.nh.compute_neighborhoods_delaunay(include_center: bool=True, key_added: str='_neighborhoods', key_adjacency_matrix: str='_adjacency_matrix')` Compute the neighborhoods of each cell based on a Delaunay triangulation. This method identifies the neighbors for each cell and computes the frequency of each cell type within these k neighbors. **Parameters** include_center : bool, optional Whether to include the center cell in the neighborhood. Default is True. key_added : str, optional The key under which the computed neighborhoods will be stored in the resulting DataArray. Default is '_neighborhoods'. key_adjacency_matrix : str, optional The key under which the computed adjacency matrix will be stored in the resulting DataArray. Default is '_adjacency_matrix'. **Returns** xarray.Dataset A merged xarray Dataset containing the original data and the computed neighborhoods. ### `ds.nh.add_neighborhood_obs(features: Union[str, List[str]]=['degree', 'homophily', 'inter_label_connectivity', 'di ...])` Adds neighborhood observations to the object by computing network features from the adjacency matrix. This method requires the `networkx` package to be installed. It checks if the adjacency matrix is present in the object, constructs a graph from the adjacency matrix, and computes network features. **Parameters** features : str or List[str], optional The network features to compute. Possible features are ['degree', 'closeness_centrality', 'betweenness_centrality', 'homophily', 'inter_label_connectivity', 'diversity_index']. **Raises** ImportError If the `networkx` package is not installed. AssertionError If the adjacency matrix is not found in the object. **Notes** The adjacency matrix should be computed and stored in the object before calling this method. You can compute the adjacency matrix using methods from the `nh` module, such as `nh.compute_neighborhoods_radius()`. ### `ds.nh.compute_graph_features(features: Union[str, List[str]]=['num_nodes', 'num_edges', 'density', 'modularity', 'as ...])` Compute various graph features from the adjacency matrix of the data. **Parameters** features : Union[str, List[str]], optional A single feature or a list of features to compute. Valid features are: "num_nodes", "num_edges", "density", "modularity", "assortativity". Default is ["num_nodes", "num_edges", "density", "modularity", "assortativity"]. **Returns** dict A dictionary where keys are the names of the computed features and values are the corresponding feature values. **Raises** ImportError If the `networkx` package is not installed. AssertionError If required layers (OBS, LA_PROPERTIES, ADJACENCY_MATRIX) are not found in the object. If no valid features are provided. If the LABELS feature is not found in the observation layer. **Notes** This method requires the `networkx` package to be installed. If you want to compute the modularity of the network, you will also need to install the `python-louvain` package. --- # The plot (`pl`) accessor The plotting accessor provides several plotting functions. ## `ds.pl` (PlotAccessor) Adds plotting functions to the image container. Available on any `xarray.Dataset` after `import spatialproteomics as sp`, as `ds.pl`. ### `ds.pl.colorize(colors: List[str]=['#e6194B', '#3cb44b', '#ffe119', '#4363d8', '#f58231', ...], background: str='black', normalize: bool=True, merge: bool=True, layer_key: str='_image', amin: Union[float, List, NoneType]=None, amax: Union[float, List, NoneType]=None) -> xr.Dataset` Colorizes a stack of images. **Parameters** colors : List[str], optional A list of strings that denote the color of each channel. background : str, optional Background color of the colorized image. Default is "black". normalize : bool, optional Normalize the image prior to colorizing it. Default is True. When amin and amax are not specified, this uses the 3- and 99.8-percentile for normalization. If they are specified, it uses them as absolute values. merge : True, optional Merge the channel dimension. Default is True. layer_key : str, optional The key of the layer containing the image. Default is '_image'. amin : float | List, optional The minimum absolute value for normalization (either global or per channel). Default is None. amax : float | List, optional The maximum absolute value for normalization (either global or per channel). Default is None. **Returns** xr.Dataset The image container with the colorized image stored in '_plot'. ### `ds.pl.show(render_image: bool=True, render_segmentation: bool=False, render_labels: bool=False, render_neighborhoods: bool=False, ax=None, legend_image: bool=True, legend_segmentation: bool=True, legend_label: bool=True, legend_neighborhoods: bool=True, background: str='black', downsample: int=1, label_order: Optional[List]=None, neighborhood_order: Optional[List]=None, legend_kwargs: dict={'framealpha': 1}, segmentation_kwargs: dict={}, label_kwargs: dict={}, neighborhood_kwargs: dict={}) -> xr.Dataset` Display an image with optional rendering elements. Can be used to render intensities, segmentation masks and labels, either individually or all at once. **Parameters** render_image : bool Whether to render the image with channel intensities. Default is True. render_segmentation : bool Whether to render segmentation. Default is False. render_labels : bool Whether to render labels. Default is False. render_neighborhoods : bool Whether to render neighborhoods. Default is False. ax The matplotlib axes to plot on. If None, the current axes will be used. legend_image : bool Whether to show the channel legend. Default is True. legend_segmentation : bool Whether to show the segmentation legend (only becomes relevant when dealing with multiple segmentation layers, e. g. when using cellpose). Default is False. legend_label : bool Whether to show the label legend. Default is True. legend_neighborhoods : bool Whether to show the neighborhood legend. Default is True. background : str Background color of the image. Default is "black". downsample : int Downsample factor for the image. Default is 1 (no downsampling). label_order : list A list specifying the order of the label indices. Default is None. neighborhood_order : list A list specifying the order of the neighborhood indices. Default is None. legend_kwargs : dict Keyword arguments for configuring the legend. Default is {"framealpha": 1}. segmentation_kwargs : dict Keyword arguments for rendering the segmentation. Default is {}. label_kwargs : dict Keyword arguments for rendering the labels. Default is {}. neighborhood_kwargs : dict Keyword arguments for rendering the neighborhoods. Default is {}. **Returns** xr.Dataset The dataset object including the plot. **Raises** - AssertionError: If no rendering element is specified. **Notes** - This method displays an image with optional rendering elements such as intensities, labels, and segmentation. - The `render_intensities`, `render_labels`, and `render_segmentation` parameters control which rendering elements are displayed. - The `ax` parameter allows specifying a specific matplotlib axes to plot on. If None, the current axes will be used. - The `show_channel_legend` and `show_label_legend` parameters control whether to show the channel and label legends, respectively. - The `legend_kwargs`, `segmentation_kwargs`, and `label_kwargs` parameters allow configuring the appearance of the legend, segmentation, and labels, respectively. ### `ds.pl.annotate(variable: str='cell', layer_key: str='_obs', highlight: list=[], text_kwargs: dict={'color': 'w', 'fontsize': 12}, highlight_kwargs: dict={'color': 'w', 'fontsize': 16, 'fontweight': 'bold'}, bbox: Optional[List]=None, format_string: str='', ax=None) -> xr.Dataset` Annotates cells with their respective number on the plot. **Parameters** variable : str, optional The feature in the observation table to be used for cell annotation. Default is "cell". layer_key : str, optional The key representing the layer in the data object. Default is Layers.OBS. highlight : list, optional A list containing cell IDs to be highlighted in the plot. text_kwargs : dict, optional Keyword arguments passed to matplotlib's text function for normal cell annotations. highlight_kwargs : dict, optional Similar to 'text_kwargs' but specifically for the cell IDs passed via 'highlight'. bbox : Union[List, None], optional A list containing bounding box coordinates [xmin, xmax, ymin, ymax] to annotate cells only within the box. Default is None, which annotates all cells. format_string : str, optional The format string used to format the cell annotation text. Default is "" (no formatting). ax : matplotlib.axes, optional The matplotlib axis to plot on. If not provided, the current axis will be used. **Returns** xr.Dataset The updated image container. **Notes** - The function annotates cells with their respective values from the selected feature. - You can highlight specific cells in the plot using the 'highlight' parameter. - Bounding box coordinates can be provided via 'bbox' to annotate cells only within the specified box. - 'format_string' can be used to format the cell annotation text (e.g., "{t:.2f}" for float formatting). ### `ds.pl.render_segmentation(layer_key: str='_segmentation', colors: List[str]=['white', '#e6194B', '#3cb44b', '#ffe119', '#4363d8', ' ...], alpha: float=0.0, alpha_boundary: float=1.0, mode: str='inner') -> xr.Dataset` Renders the segmentation mask with optional alpha blending and boundary rendering. Parameters: layer_key (str, optional): The key of the layer containing the segmentation mask. Default is '_segmentation'. colors (List[str], optional): A list of colors to be used for rendering the segmentation mask. Default is ['white']. alpha (float, optional): The alpha value for blending the segmentation mask with the plot. Default is 0.0. alpha_boundary (float, optional): The alpha value for rendering the boundary of the segmentation mask. Default is 1.0. mode (str, optional): The mode for rendering the segmentation mask. Possible values are "inner" and "outer". Default is "inner". Returns: xr.Dataset: The modified xarray Dataset. Raises: AssertionError: If no segmentation layer is found in the object. Note: - The segmentation mask must be added to the object before calling this method. - The segmentation mask is expected to have a single channel with integer labels. - The rendered segmentation mask will be added as a new layer to the object. ### `ds.pl.render_labels(alpha: float=1.0, alpha_boundary: float=1.0, mode: str='inner', override_color: Optional[str]=None) -> xr.Dataset` Renders cell type labels on the plot. Parameters: alpha (float, optional): The transparency of the labels. Defaults to 1.0. alpha_boundary (float, optional): The transparency of the label boundaries. Defaults to 1.0. mode (str, optional): The mode of rendering. Can be "inner" or "outer". Defaults to "inner". override_color (str, optional): The color to override the label colors. Defaults to None. Returns: xr.Dataset: The modified dataset with rendered labels. Raises: AssertionError: If no labels are found in the object. Notes: - This method requires labels to be present in the object. Add labels first using `la.predict_cell_types_argmax()` or `tl.astir()`. - The `mode` parameter determines whether the labels are rendered inside or outside the label boundaries. - The `override_color` parameter can be used to override the label colors with a single color. ### `ds.pl.render_neighborhoods(style: str='neighborhoods', alpha: float=1.0, alpha_boundary: float=1.0, boundary_color: str='dimgray', boundary_thickness: int=3, dilation_strength: int=40, erosion_strength: int=35) -> xr.Dataset` Render neighborhoods on the spatial data. **Parameters** style : str, optional The style of rendering, either 'cells' or 'neighborhoods'. Default is 'neighborhoods'. alpha : float, optional The alpha transparency for the rendered image. Default is 1.0. alpha_boundary : float, optional The alpha transparency for the boundary of the rendered image. Default is 1.0. boundary_color : str, optional The color of the boundary lines. Default is 'dimgray'. boundary_thickness : int, optional The thickness of the boundary lines. Default is 3. dilation_strength : int, optional The strength of the dilation applied to the cells. Default is 40. erosion_strength : int, optional The strength of the erosion applied to the cells. Default is 35. **Returns** xr.Dataset The dataset with the rendered neighborhoods. **Raises** AssertionError If no neighborhoods are found in the object or if the style is not 'cells' or 'neighborhoods'. ### `ds.pl.render_obs(feature: str=None, cmap: str='viridis', background: str='black', alpha: float=1.0, alpha_boundary: float=1.0, mode: str='inner', vmin: Optional[float]=None, vmax: Optional[float]=None) -> xr.Dataset` Render the observation layer with the specified feature and colormap. **Parameters** feature : str, optional The feature to be rendered from the observation layer. Default is None. cmap : str, optional The colormap to be used for rendering. Default is "viridis". background : str, optional The background color of the rendered image. Default is "black". alpha : float, optional The alpha value for the rendered image. Default is 1.0. alpha_boundary : float, optional The alpha value for the boundaries in the rendered image. Default is 1.0. mode : str, optional The mode for rendering. Default is "inner". vmin : float, optional The minimum value for colormap normalization. Default is None. vmax : float, optional The maximum value for colormap normalization. Default is None. **Returns** xr.Dataset The dataset with the rendered observation layer. **Raises** AssertionError If the observation layer or segmentation layer is not found in the object. If the specified feature is not found in the observation layer. ### `ds.pl.imshow(legend_image: bool=False, legend_segmentation: bool=False, legend_label: bool=False, legend_neighborhoods: bool=False, legend_obs: bool=False, downsample: int=1, label_order: Optional[List[str]]=None, neighborhood_order: Optional[List[str]]=None, legend_kwargs: dict={'framealpha': 1}, cbar_kwargs: dict={}, ax=None)` Plots the image after rendering certain layers. Meant to be used in conjunction with pl.colorize(), pl.render_segmentation, pl.render_label() or pl.render_neighborhoods(). For a more high level wrapper, please refer to pl.show() instead. **Parameters** legend_image : bool, optional Show the legend for the channels. Default is False. legend_segmentation : bool, optional Show the legend for the segmentation. Default is False. legend_label : bool, optional Show the label legend. Default is False. legend_neighborhoods : bool, optional Show the neighborhood legend. Default is False. legend_obs: bool, optional Show the observation colorbar. Default is False. downsample : int, optional Downsample factor for the image. Default is 1. legend_kwargs : dict, optional Additional keyword arguments for configuring the legend. Default is {"framealpha": 1}. cbar_kwargs : dict, optional Additional keyword arguments for configuring the colorbar. Default is {}. label_order: Optional[List[str]], optional The order in which the labels should be displayed. Default is None. neighborhood_order: Optional[List[str]], optional The order in which the neighborhoods should be displayed. Default is None. ax : matplotlib.axes, optional The matplotlib axis to plot on. If not provided, the current axis will be used. **Returns** xr.Dataset The updated image container. **Notes** - The function is used to plot images in conjunction with 'im.colorize' and 'la.render_label'. - The appearance of the plot and the inclusion of legends can be controlled using the respective parameters. ### `ds.pl.scatter_labels(legend: bool=True, size: float=1.0, alpha: float=0.9, zorder: int=10, render_edges: bool=False, render_self_edges: bool=False, ax=None, legend_kwargs: dict={'framealpha': 1}, scatter_kwargs: dict={}) -> xr.Dataset` Scatter plot of labeled cells. **Parameters** legend : bool, optional Whether to show the legend. Default is True. size : float, optional Size of the scatter markers. Default is 1.0. alpha : float, optional Transparency of the scatter markers. Default is 0.9. zorder : int, optional The z-order of the scatter markers. Default is 10. render_edges : bool, optional Whether to render the edges between cells within the same neighborhood. Default is False. render_self_edges : bool, optional Whether to render the edges between cells and themselves. Default is False. ax : matplotlib.axes.Axes, optional The axes on which to plot the scatter. If not provided, the current axes will be used. legend_kwargs : dict, optional Additional keyword arguments for configuring the legend. Default is {"framealpha": 1}. scatter_kwargs : dict, optional Additional keyword arguments for configuring the scatter plot. **Returns** xr.Dataset The modified spatialproteomics object. ### `ds.pl.scatter(feature: str, palette: dict=None, legend: bool=True, layer_key: str='_obs', size: float=1.0, alpha: float=0.9, zorder: int=10, ax=None, legend_kwargs: dict={'framealpha': 1}, scatter_kws: dict={}) -> xr.Dataset` Create a scatter plot of some feature. At the moment, only categorical features are supported. **Parameters** feature : str The feature to be plotted. palette : dict, optional A dictionary mapping feature values to colors. If not provided, a default palette will be used. legend : bool, optional Whether to show the legend. Default is True. layer_key : str, optional The key of the layer to be plotted. Default is Layers.OBS. size : float, optional The size of the scatter points. Default is 1.0. alpha : float, optional The transparency of the scatter points. Default is 0.9. zorder : int, optional The z-order of the scatter points. Default is 10. ax : object, optional The matplotlib axes object to plot on. If not provided, the current axes will be used. legend_kwargs : dict, optional Additional keyword arguments for configuring the legend. Default is {"framealpha": 1}. scatter_kws : dict, optional Additional keyword arguments for configuring the scatter plot. Default is {}. **Returns** xr.Dataset The original data object. **Raises** - AssertionError: If the layer_key is not found in the data object. - AssertionError: If the feature is not found in the specified layer. - AssertionError: If the X or Y coordinates are not found in the specified layer. - AssertionError: If the number of unique feature values is greater than 10 and no color_scheme is provided. - AssertionError: If not all unique feature values are present in the provided palette. ### `ds.pl.add_box(xlim: List[int], ylim: List[int], color: str='w', linewidth: float=2, ax=None)` Adds a rectangular box to the current plot. **Parameters** xlim : List[int] The x-bounds of the box [xstart, xstop]. ylim : List[int] The y-bounds of the box [ymin, ymax]. color : str, optional The color of the box. Default is "w" (white). linewidth : float, optional The linewidth of the box. Default is 2. ax : matplotlib.axes, optional The matplotlib axis to plot on. If not provided, the current axis will be used. **Returns** xr.Dataset The updated spatialproteomics object. ### `ds.pl.autocrop(padding: int=50, downsample: int=10, key: str='_image', channel: Optional[str]=None) -> xr.Dataset` Crop the image so that the background surrounding the tissue/TMA gets cropped away. **Parameters** padding : int The padding to be added around the cropped image in pixels. Default is 50. downsample : int The downsampling factor for the image. Default is 10. key : str The key of the image to be cropped. Default is Layers.IMAGE. channel : str, optional The channel used for cropping. Default is None, which defaults to using the first available channel. **Returns** xr.Dataset The cropped image. --- # The tool (`tl`) accessor The tool accessor contains methods that utilize third-party tools to enable operations such as segmentation or cell type prediction. ## `ds.tl` (ToolAccessor) The tool accessor enables the application of external tools such as StarDist or Astir. Available on any `xarray.Dataset` after `import spatialproteomics as sp`, as `ds.tl`. ### `ds.tl.cellpose(channel: Optional[str]=None, key_added: str='_segmentation', diameter: float=None, channel_settings: list=None, num_iterations: int=2000, cellprob_threshold: float=0.0, flow_threshold: float=0.4, batch_size: int=8, gpu: bool=True, model_type: str='cyto3', pretrained_model: str='cpsam', postprocess_func: Callable=, return_diameters: bool=False, **kwargs)` Segment cells using Cellpose. Adds a layer to the spatialproteomics object with dimension (X, Y) or (C, X, Y) dependent on whether channel argument is specified or not. **Parameters** channel : str, optional Channel to use for segmentation. If None, all channels are used for independent segmentation. key_added : str, optional Key to assign to the segmentation results. diameter : float, optional Expected cell diameter in pixels. channel_settings : List[int], optional Channels for Cellpose to use for segmentation. If [0, 0], independent segmentation is performed on all channels. If it is anything else (e. g. [1, 2]), joint segmentation is attempted. num_iterations : int, optional Maximum number of iterations for segmentation. cellprob_threshold : float, optional Threshold for cell probability. flow_threshold : float, optional Threshold for flow. batch_size : int, optional Batch size for segmentation. gpu : bool, optional Whether to use GPU for segmentation. model_type : str, optional Type of Cellpose model to use. pretrained_model : str, optional Pretrained model to use for Cellpose (4+). postprocess_func : Callable, optional Function to apply to the segmentation masks after prediction. return_diameters : bool, optional Whether to return the cell diameters. **Returns** xr.Dataset Dataset containing original data and segmentation mask. **Notes** This method requires the 'cellpose' package to be installed. ### `ds.tl.stardist(channel: Optional[str]=None, key_added: str='_segmentation', scale: float=3, n_tiles: int=12, normalize: bool=True, predict_big: bool=False, postprocess_func: Callable=, **kwargs) -> xr.Dataset` Apply StarDist algorithm to perform instance segmentation on the nuclear image. **Parameters** channel : str, optional Channel to use for segmentation. If None, all channels are used. key_added : str, optional Key to write the segmentation results to. scale : float, optional Scaling factor for the StarDist model (default is 3). n_tiles : int, optional Number of tiles to split the image into for prediction (default is 12). normalize : bool, optional Flag indicating whether to normalize the nuclear image (default is True). nuclear_channel : str, optional Name of the nuclear channel in the image (default is "DAPI"). predict_big : bool, optional Flag indicating whether to use the 'predict_instances_big' method for large images (default is False). postprocess_func : Callable, optional Function to apply to the segmentation masks after prediction (default is lambda x: x). **kwargs : dict, optional Additional keyword arguments to be passed to the StarDist prediction method. **Returns** obj : xr.Dataset Xarray dataset containing the segmentation mask and centroids. **Raises** ValueError If the object already contains a segmentation mask. ### `ds.tl.instanseg(channel: Optional[str]=None, key_added: str='_segmentation', model_type: str='fluorescence_nuclei_and_cells', pixel_size: Optional[float]=None, target: str='cells', normalize: bool=True, method: str='auto', tile_size: int=512, batch_size: int=1, device: Optional[str]=None, verbosity: int=1, postprocess_func: Callable=, **kwargs) -> xr.Dataset` Apply InstanSeg (https://github.com/instanseg/instanseg) to perform instance segmentation of nuclei and/or cells. InstanSeg is channel-invariant, so in contrast to cellpose and stardist, all selected channels are segmented jointly and a single segmentation mask is returned (e. g. a nuclear marker together with membrane markers). Requires `pip install instanseg-torch`. **Parameters** channel : str or list, optional Channel(s) to use for segmentation. If None, all channels are used. key_added : str, optional Key to write the segmentation results to. If target is "all_outputs", the cell mask is written to key_added and the nuclear mask to f"{key_added}_nuclei". model_type : str, optional Name of the pretrained InstanSeg model, e. g. "fluorescence_nuclei_and_cells" (default) or "brightfield_nuclei". pixel_size : float, optional Pixel size of the image in microns. The image is rescaled to the pixel size the model was trained on (0.5 microns for the default model). If None, the image is not rescaled, which may lead to inaccurate results. target : str, optional Which segmentation to return: "cells" (default), "nuclei" or "all_outputs" (both). Models that only predict nuclei always return the nuclear segmentation. normalize : bool, optional Whether InstanSeg should percentile-normalize the image (default is True). method : str, optional "small" processes the whole image at once, "medium" uses tiled inference. "auto" (default) picks "medium" for large images. tile_size : int, optional Tile size for tiled inference (default is 512). Only used if method is "medium". batch_size : int, optional Number of tiles processed simultaneously (default is 1). Only used if method is "medium". device : str, optional Device to run the model on, e. g. "cpu" or "cuda". If None, the device is chosen automatically. verbosity : int, optional Verbosity of InstanSeg: 0 is silent, 1 is normal, 2 is verbose (default is 1). postprocess_func : Callable, optional Function to apply to each segmentation mask after prediction (default is lambda x: x). **kwargs : dict, optional Additional keyword arguments passed to InstanSeg's eval_small_image or eval_medium_image. **Returns** xr.Dataset Xarray dataset containing the segmentation mask(s). **Raises** KeyError If the key the segmentation should be written to already exists. ### `ds.tl.mesmer(key_added: str='_segmentation', channel: Optional[List]=None, postprocess_func: Callable=, **kwargs)` This function is deprecated. ### `ds.tl.astir(marker_dict: dict, key: str='_intensity', threshold: float=0, seed: int=42, learning_rate: float=0.001, batch_size: float=64, n_init: int=5, n_init_epochs: int=5, max_epochs: int=500, cell_id_col: str='cell_id', cell_type_col: str='cell_type', **kwargs)` This method predicts cell types from an expression matrix using the Astir algorithm. **Parameters** marker_dict : dict Dictionary mapping markers to cell types. Can also include cell states. Example: {"cell_type": {'B': ['PAX5'], 'T': ['CD3'], 'Myeloid': ['CD11b']}} key : str, optional Layer to use as expression matrix. threshold : float, optional Certainty threshold for astir to assign a cell type. Defaults to 0. seed : int, optional Random seed. Defaults to 42. learning_rate : float, optional Learning rate. Defaults to 0.001. batch_size : float, optional Batch size. Defaults to 64. n_init : int, optional Number of initializations. Defaults to 5. n_init_epochs : int, optional Number of epochs for each initialization. Defaults to 5. max_epochs : int, optional Maximum number of epochs. Defaults to 500. cell_id_col : str, optional Column name for cell IDs. Defaults to "cell_id". cell_type_col : str, optional Column name for cell types. Defaults to "cell_type". **kwargs Additional keyword arguments passed to astir's fit_type method. **Raises** ValueError If no expression matrix was present or the image is not of type uint8. **Returns** DataArray A DataArray with the assigned cell types. ### `ds.tl.spotiflow(channel: Optional[str]=None, image_key: str='_image', segmentation_key: str='_segmentation', key_added_suffix: str='_spots', pretrained_name: str='general', return_spots: bool=False, model_kwargs: Optional[dict]={'verbose': False}, predict_kwargs: Optional[dict]={'verbose': False})` This method detects spots in the image using the Spotiflow algorithm and counts the number of spots per cell, per channel. **Parameters** channel : str, optional Channel to use for spot detection. If None, all channels are used. image_key : str, optional Key to the image layer in the spatialproteomics object. Default is 'image'. segmentation_key : str, optional Key to the segmentation layer in the spatialproteomics object. Default is 'segmentation'. key_added_suffix : str, optional Suffix to add to the channel names in the resulting counts. Default is '_spots'. pretrained_name : str, optional Name of the pretrained model to use for spot detection. Default is 'general'. return_spots : bool, optional Whether to return the detected spots as a DataFrame. Default is False. model_kwargs : dict, optional Additional keyword arguments to pass to the Spotiflow model. Default is {'verbose': False}. predict_kwargs : dict, optional Additional keyword arguments to pass to the Spotiflow prediction method. Default is {'verbose': False}. **Returns** xr.Dataset Dataset containing the counts of spots per cell, per channel. If return_spots is True, also returns a DataFrame of the detected spots with their coordinates and assigned cell IDs. ### `ds.tl.convert_to_anndata(expression_matrix_key: str='_intensity', obs_key: str='_obs', additional_layers: Optional[dict]=None, additional_obsm: Optional[dict]=None, additional_uns: Optional[dict]=None)` Convert the spatialproteomics object to an anndata.AnnData object. The resulting AnnData object does not store the original image or segmentation mask. **Parameters** expression_matrix_key : str, optional The key of the expression matrix in the spatialproteomics object. Default is '_intensity'. obs_key : str, optional The key of the observation data in the spatialproteomics object. Default is '_obs'. additional_layers : dict, optional Additional layers to include in the anndata.AnnData object. The keys are the names of the layers and the values are the corresponding keys in the spatialproteomics object. additional_obsm : dict, optional Additional obsm data to include in the anndata.AnnData object. The keys are the names of the obsm data and the values are the corresponding keys in the spatialproteomics object. additional_uns : dict, optional Additional uns data to include in the anndata.AnnData object. The keys are the names of the uns data and the values are the corresponding keys in the spatialproteomics object. **Returns** adata : anndata.AnnData The converted anndata.AnnData object. **Raises** AssertionError If additional layers, obsm entries, or uns entries are not found in the spatialproteomics object, or if additional layers or obsm entries have incompatible shapes. **Notes** - The expression matrix is extracted from the spatialproteomics object using the provided expression matrix key. - If additional layers are specified, they are extracted from the spatialproteomics object and added to the anndata.AnnData object. - If obs_key is present in the spatialproteomics object, it is used to create the obs DataFrame of the anndata.AnnData object. - If additional_uns is specified, the corresponding uns data is extracted from the spatialproteomics object and added to the anndata.AnnData object. ### `ds.tl.convert_to_spatialdata(image_key: str='_image', segmentation_key: str='_segmentation', **kwargs)` Convert the spatialproteomics object to a spatialdata object. Parameters: image_key (str): The key of the image data in the object. Defaults to Layers.IMAGE. segmentation_key (str): The key of the segmentation data in the object. Defaults to Layers.SEGMENTATION. **kwargs: Additional keyword arguments to be passed to the convert_to_anndata method. Returns: spatial_data_object (spatialdata.SpatialData): The converted spatialdata object. ## Functions for the spatialdata backend (`sp.tl`) These functions take a `spatialdata.SpatialData` object as first argument and modify it in place unless `copy=True`. ### `sp.tl.cellpose(sdata, channel: Optional[str]=None, image_key: str='image', key_added: str='segmentation', data_key: Optional[str]=None, copy: bool=False, **kwargs)` This function runs the cellpose segmentation algorithm on the provided image data. It extracts the image data from the spatialdata object, applies the cellpose algorithm, and adds the segmentation masks to the spatialdata object. The segmentation masks are stored in the labels attribute of the spatialdata object. The function also handles multiple channels by iterating over the channels and applying the segmentation algorithm to each channel separately. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the image data. channel (Optional[str]): The channel(s) to be used for segmentation. If None, all channels will be used. image_key (str, optional): The key for the image data in the spatialdata object. Defaults to image. key_added (str, optional): The key under which the segmentation masks will be stored in the labels attribute of the spatialdata object. Defaults to segmentation. data_key (Optional[str], optional): The key for the image data in the spatialdata object. If None, the image_key will be used. Defaults to None. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. **kwargs: Additional keyword arguments to be passed to the cellpose algorithm. ### `sp.tl.stardist(sdata, channel: Optional[str]=None, image_key: str='image', key_added: str='segmentation', data_key: Optional[str]=None, copy: bool=False, **kwargs)` This function runs the stardist segmentation algorithm on the provided image data. It extracts the image data from the spatialdata object, applies the stardist algorithm, and adds the segmentation masks to the spatialdata object. The segmentation masks are stored in the labels attribute of the spatialdata object. The function also handles multiple channels by iterating over the channels and applying the segmentation algorithm to each channel separately. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the image data. channel (Optional[str]): The channel(s) to be used for segmentation. If None, all channels will be used. image_key (str, optional): The key for the image data in the spatialdata object. Defaults to image. key_added (str, optional): The key under which the segmentation masks will be stored in the labels attribute of the spatialdata object. Defaults to segmentation. data_key (Optional[str], optional): The key for the image data in the spatialdata object. If None, the image_key will be used. Defaults to None. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. **kwargs: Additional keyword arguments to be passed to the stardist algorithm. ### `sp.tl.instanseg(sdata, channel: Optional[str]=None, image_key: str='image', key_added: str='segmentation', data_key: Optional[str]=None, copy: bool=False, **kwargs)` This function runs the InstanSeg segmentation algorithm on the provided image data. It extracts the image data from the spatialdata object, applies InstanSeg, and adds the segmentation masks to the labels attribute of the spatialdata object. In contrast to cellpose and stardist, InstanSeg is channel-invariant, i. e. all selected channels are segmented jointly and a single segmentation mask is returned. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the image data. channel (Optional[str]): The channel(s) to be used for segmentation. If None, all channels will be used. image_key (str, optional): The key for the image data in the spatialdata object. Defaults to image. key_added (str, optional): The key under which the segmentation mask will be stored in the labels attribute of the spatialdata object. If target="all_outputs", the cell mask is stored under key_added and the nuclear mask under f"{key_added}_nuclei". Defaults to segmentation. data_key (Optional[str], optional): The key for the image data in the spatialdata object. If None, the image_key will be used. Defaults to None. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. **kwargs: Additional keyword arguments passed to the InstanSeg wrapper, e.g. model_type, pixel_size, target, method or device (see ToolAccessor.instanseg). ### `sp.tl.mesmer(sdata, channel: Optional[str]=None, image_key: str='image', key_added: str='segmentation', data_key: Optional[str]=None, copy: bool=False, **kwargs)` This function is deprecated and will be removed in a future version of spatialproteomics. ### `sp.tl.astir(sdata, marker_dict: dict, table_key='table', threshold: float=0, seed: int=42, learning_rate: float=0.001, batch_size: float=64, n_init: int=5, n_init_epochs: int=5, max_epochs: int=500, cell_id_col: str='cell_id', cell_type_col: str='cell_type', copy: bool=False, **kwargs)` This function applies the ASTIR algorithm to predict cell types based on the expression matrix. It extracts the expression matrix from the spatialdata object, applies the ASTIR algorithm, and adds the predicted cell types to the spatialdata object. The predicted cell types are stored in the obs attribute of the AnnData object in the tables attribute of the spatialdata object. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the expression matrix. marker_dict (dict): A dictionary containing the marker genes for each cell type. table_key (str, optional): The key under which the expression matrix is stored in the tables attribute of the spatialdata object. Defaults to "table". threshold (float, optional): The threshold value to be used for the ASTIR algorithm. Defaults to 0. seed (int, optional): The random seed to be used for the ASTIR algorithm. Defaults to 42. learning_rate (float, optional): The learning rate to be used for the ASTIR algorithm. Defaults to 0.001. batch_size (float, optional): The batch size to be used for the ASTIR algorithm. Defaults to 64. n_init (int, optional): The number of initializations to be used for the ASTIR algorithm. Defaults to 5. n_init_epochs (int, optional): The number of initial epochs to be used for the ASTIR algorithm. Defaults to 5. max_epochs (int, optional): The maximum number of epochs to be used for the ASTIR algorithm. Defaults to 500. cell_id_col (str, optional): The name of the column containing the cell IDs in the expression matrix. Defaults to "cell_id". cell_type_col (str, optional): The name of the column containing the cell types in the expression matrix. Defaults to "cell_type". copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. **kwargs: Additional keyword arguments to be passed to the ASTIR algorithm. ### `sp.tl.spotiflow(sdata, channel: Optional[str]=None, image_key: str='image', segmentation_key: str='segmentation', table_key: str='table', key_added: str='spotiflow', key_added_suffix: str='_spots', pretrained_name: str='general', model_kwargs: Optional[dict]={'verbose': False}, predict_kwargs: Optional[dict]={'verbose': False}, copy: bool=False)` This function applies the Spotiflow algorithm to detect spots in the provided image data. It extracts the image data from the spatialdata object, applies the Spotiflow algorithm, and adds the detected spots to the spatialdata object. The detected spots are stored in the points attribute of the spatialdata object. Args: sdata (spatialdata.SpatialData): The spatialdata object containing the image data. channel (Optional[str]): The channel(s) to be used for spot detection. If None, all channels will be used. image_key (str, optional): The key for the image data in the spatialdata object. Defaults to "image". segmentation_key (str, optional): The key for the segmentation mask in the spatialdata object. Defaults to "segmentation". table_key (str, optional): The key for the table in the spatialdata object. Defaults to "table". key_added (str, optional): The key under which the detected spots will be stored in the points attribute of the spatialdata object. Defaults to "spotiflow". key_added_suffix (str, optional): The suffix to be added to the channel names in the detected spots. Defaults to "_spots". pretrained_name (str, optional): The name of the pretrained model to be used for spot detection. Defaults to "general". model_kwargs (Optional[dict], optional): Additional keyword arguments to be passed to the Spotiflow model. Defaults to {'verbose': False}. predict_kwargs (Optional[dict], optional): Additional keyword arguments to be passed to the Spotiflow prediction method. Defaults to {'verbose': False}. copy (bool, optional): Whether to create a copy of the spatialdata object. Defaults to False. Returns: spatialdata.SpatialData: The spatialdata object with the detected spots added to the points attribute. --- # The image container The image container wraps methods that needs to run on multiple images. For example, clustering neighborhoods is typically performed on all samples simultaneously. ## `sp.ImageContainer` This class is used to store multiple SpatialProteomics objects and perform operations on them. ### `sp.ImageContainer.compute_neighborhoods(neighborhood_method: str='radius', radius=100, knn=10, k=5, overwrite: bool=False, seed: int=0)` Compute neighborhoods for spatial proteomics objects using the specified method and perform clustering. **Parameters** neighborhood_method : str, optional The method to use for computing neighborhoods. Must be one of 'radius', 'knn', or 'delaunay'. Default is 'radius'. radius : int, optional The radius to use for the 'radius' neighborhood method. Default is 100. knn : int, optional The number of nearest neighbors to use for the 'knn' neighborhood method. Default is 10. k : int, optional The number of clusters to form using K-Means clustering. Default is 5. overwrite : bool, optional Whether to overwrite existing neighborhoods in the objects. Default is False. seed : int, optional The random seed to use for K-Means clustering. Default is 0. **Returns** dict A dictionary of spatial proteomics objects with computed neighborhoods and clusters. **Raises** ValueError If neighborhoods are already present in the objects and `overwrite` is set to False. If `neighborhood_method` is not one of 'radius', 'knn', or 'delaunay'. If there is an error in the clustering process. ### `sp.ImageContainer.get_neighborhood_composition(standardize: bool=True, neighborhood_col: str='_neighborhoods') -> pd.DataFrame` Get the composition of neighborhoods across all objects in the ImageContainer. **Parameters** standardize : bool, optional Whether to standardize the composition of neighborhoods. Default is True. neighborhood_col : str, optional The column in the object's DataFrame that contains the neighborhood labels. Default is 'neighborhoods'. **Returns** pd.DataFrame A DataFrame containing the composition of neighborhoods across all objects in the ImageContainer.