.ipynb .pdf

Contents

Interoperability

This notebook shows some way that you can import and export data from spatialproteomics.

%reload_ext autoreload
%autoreload 2

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.

# 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("../../tests/test_files/ds_neighborhoods.zarr")
ds

<xarray.Dataset> Size: 171kB
Dimensions:               (cells: 56, cells_2: 56, channels: 5, y: 101, x: 101,
                           la_features: 2, labels: 4, la_props: 2,
                           neighborhoods: 5, nh_props: 2, features: 6)
Coordinates:
  * cells                 (cells) int64 448B 1 2 3 4 5 6 7 ... 51 52 53 54 55 56
  * cells_2               (cells_2) int64 448B 1 2 3 4 5 6 ... 51 52 53 54 55 56
  * channels              (channels) <U11 220B 'DAPI' 'PAX5' 'CD3' 'CD4' 'CD8'
  * features              (features) <U14 336B 'CD4_binarized' ... 'centroid-1'
  * la_features           (la_features) object 16B 'labels_0' 'labels_1'
  * la_props              (la_props) <U6 48B '_color' '_name'
  * labels                (labels) int64 32B 1 2 3 4
  * neighborhoods         (neighborhoods) int64 40B 1 2 3 4 5
  * nh_props              (nh_props) <U6 48B '_color' '_name'
  * x                     (x) int64 808B 1600 1601 1602 1603 ... 1698 1699 1700
  * y                     (y) int64 808B 2100 2101 2102 2103 ... 2198 2199 2200
Data variables:
    _adjacency_matrix     (cells, cells_2) int64 25kB dask.array<chunksize=(56, 56), meta=np.ndarray>
    _image                (channels, y, x) uint8 51kB dask.array<chunksize=(5, 101, 101), meta=np.ndarray>
    _intensity            (cells, channels) float64 2kB dask.array<chunksize=(56, 5), meta=np.ndarray>
    _la_layers            (cells, la_features) object 896B dask.array<chunksize=(56, 2), meta=np.ndarray>
    _la_properties        (labels, la_props) object 64B dask.array<chunksize=(4, 2), meta=np.ndarray>
    _neighborhoods        (cells, labels) float64 2kB dask.array<chunksize=(56, 4), meta=np.ndarray>
    _nh_properties        (neighborhoods, nh_props) <U14 560B dask.array<chunksize=(5, 2), meta=np.ndarray>
    _obs                  (cells, features) float64 3kB dask.array<chunksize=(56, 6), meta=np.ndarray>
    _percentage_positive  (cells, channels) float64 2kB dask.array<chunksize=(56, 5), meta=np.ndarray>
    _segmentation         (y, x) int64 82kB dask.array<chunksize=(101, 101), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • cells: 56
  • cells_2: 56
  • channels: 5
  • y: 101
  • x: 101
  • la_features: 2
  • labels: 4
  • la_props: 2
  • neighborhoods: 5
  • nh_props: 2
  • features: 6
• Coordinates: (11)
  • cells
    (cells)
    int64
    1 2 3 4 5 6 7 ... 51 52 53 54 55 56

    array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
           19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
           37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
           55, 56])

  • cells_2
    (cells_2)
    int64
    1 2 3 4 5 6 7 ... 51 52 53 54 55 56

    array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
           19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
           37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
           55, 56])

  • channels
    (channels)
    <U11
    'DAPI' 'PAX5' 'CD3' 'CD4' 'CD8'

    array(['DAPI', 'PAX5', 'CD3', 'CD4', 'CD8'], dtype='<U11')

  • features
    (features)
    <U14
    'CD4_binarized' ... 'centroid-1'

    array(['CD4_binarized', 'CD8_binarized', '_labels', '_neighborhoods',
           'centroid-0', 'centroid-1'], dtype='<U14')

  • la_features
    (la_features)
    object
    'labels_0' 'labels_1'

    array(['labels_0', 'labels_1'], dtype=object)

  • la_props
    (la_props)
    <U6
    '_color' '_name'

    array(['_color', '_name'], dtype='<U6')

  • labels
    (labels)
    int64
    1 2 3 4

    array([1, 2, 3, 4])

  • neighborhoods
    (neighborhoods)
    int64
    1 2 3 4 5

    array([1, 2, 3, 4, 5])

  • nh_props
    (nh_props)
    <U6
    '_color' '_name'

    array(['_color', '_name'], dtype='<U6')

  • x
    (x)
    int64
    1600 1601 1602 ... 1698 1699 1700

    array([1600, 1601, 1602, 1603, 1604, 1605, 1606, 1607, 1608, 1609, 1610, 1611,
           1612, 1613, 1614, 1615, 1616, 1617, 1618, 1619, 1620, 1621, 1622, 1623,
           1624, 1625, 1626, 1627, 1628, 1629, 1630, 1631, 1632, 1633, 1634, 1635,
           1636, 1637, 1638, 1639, 1640, 1641, 1642, 1643, 1644, 1645, 1646, 1647,
           1648, 1649, 1650, 1651, 1652, 1653, 1654, 1655, 1656, 1657, 1658, 1659,
           1660, 1661, 1662, 1663, 1664, 1665, 1666, 1667, 1668, 1669, 1670, 1671,
           1672, 1673, 1674, 1675, 1676, 1677, 1678, 1679, 1680, 1681, 1682, 1683,
           1684, 1685, 1686, 1687, 1688, 1689, 1690, 1691, 1692, 1693, 1694, 1695,
           1696, 1697, 1698, 1699, 1700])

  • y
    (y)
    int64
    2100 2101 2102 ... 2198 2199 2200

    array([2100, 2101, 2102, 2103, 2104, 2105, 2106, 2107, 2108, 2109, 2110, 2111,
           2112, 2113, 2114, 2115, 2116, 2117, 2118, 2119, 2120, 2121, 2122, 2123,
           2124, 2125, 2126, 2127, 2128, 2129, 2130, 2131, 2132, 2133, 2134, 2135,
           2136, 2137, 2138, 2139, 2140, 2141, 2142, 2143, 2144, 2145, 2146, 2147,
           2148, 2149, 2150, 2151, 2152, 2153, 2154, 2155, 2156, 2157, 2158, 2159,
           2160, 2161, 2162, 2163, 2164, 2165, 2166, 2167, 2168, 2169, 2170, 2171,
           2172, 2173, 2174, 2175, 2176, 2177, 2178, 2179, 2180, 2181, 2182, 2183,
           2184, 2185, 2186, 2187, 2188, 2189, 2190, 2191, 2192, 2193, 2194, 2195,
           2196, 2197, 2198, 2199, 2200])

• Data variables: (10)
  • _adjacency_matrix
    (cells, cells_2)
    int64
    dask.array<chunksize=(56, 56), meta=np.ndarray>

    Array Chunk Bytes 24.50 kiB 24.50 kiB Shape (56, 56) (56, 56) Dask graph 1 chunks in 2 graph layers Data type int64 numpy.ndarray

  • _image
    (channels, y, x)
    uint8
    dask.array<chunksize=(5, 101, 101), meta=np.ndarray>

    Array Chunk Bytes 49.81 kiB 49.81 kiB Shape (5, 101, 101) (5, 101, 101) Dask graph 1 chunks in 2 graph layers Data type uint8 numpy.ndarray

  • _intensity
    (cells, channels)
    float64
    dask.array<chunksize=(56, 5), meta=np.ndarray>

    Array Chunk Bytes 2.19 kiB 2.19 kiB Shape (56, 5) (56, 5) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • _la_layers
    (cells, la_features)
    object
    dask.array<chunksize=(56, 2), meta=np.ndarray>

    Array Chunk Bytes 896 B 896 B Shape (56, 2) (56, 2) Dask graph 1 chunks in 2 graph layers Data type object numpy.ndarray

  • _la_properties
    (labels, la_props)
    object
    dask.array<chunksize=(4, 2), meta=np.ndarray>

    Array Chunk Bytes 64 B 64 B Shape (4, 2) (4, 2) Dask graph 1 chunks in 2 graph layers Data type object numpy.ndarray

  • _neighborhoods
    (cells, labels)
    float64
    dask.array<chunksize=(56, 4), meta=np.ndarray>

    Array Chunk Bytes 1.75 kiB 1.75 kiB Shape (56, 4) (56, 4) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • _nh_properties
    (neighborhoods, nh_props)
    <U14
    dask.array<chunksize=(5, 2), meta=np.ndarray>

    Array Chunk Bytes 560 B 560 B Shape (5, 2) (5, 2) Dask graph 1 chunks in 2 graph layers Data type

  • _obs
    (cells, features)
    float64
    dask.array<chunksize=(56, 6), meta=np.ndarray>

    Array Chunk Bytes 2.62 kiB 2.62 kiB Shape (56, 6) (56, 6) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • _percentage_positive
    (cells, channels)
    float64
    dask.array<chunksize=(56, 5), meta=np.ndarray>

    Array Chunk Bytes 2.19 kiB 2.19 kiB Shape (56, 5) (56, 5) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • _segmentation
    (y, x)
    int64
    dask.array<chunksize=(101, 101), meta=np.ndarray>

    Array Chunk Bytes 79.70 kiB 79.70 kiB Shape (101, 101) (101, 101) Dask graph 1 chunks in 2 graph layers Data type int64 numpy.ndarray

• Indexes: (11)
  • cells
    PandasIndex

    PandasIndex(Index([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
           19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
           37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
           55, 56],
          dtype='int64', name='cells'))

  • cells_2
    PandasIndex

    PandasIndex(Index([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
           19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
           37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
           55, 56],
          dtype='int64', name='cells_2'))

  • channels
    PandasIndex

    PandasIndex(Index(['DAPI', 'PAX5', 'CD3', 'CD4', 'CD8'], dtype='object', name='channels'))

  • features
    PandasIndex

    PandasIndex(Index(['CD4_binarized', 'CD8_binarized', '_labels', '_neighborhoods',
           'centroid-0', 'centroid-1'],
          dtype='object', name='features'))

  • la_features
    PandasIndex

    PandasIndex(Index(['labels_0', 'labels_1'], dtype='object', name='la_features'))

  • la_props
    PandasIndex

    PandasIndex(Index(['_color', '_name'], dtype='object', name='la_props'))

  • labels
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4], dtype='int64', name='labels'))

  • neighborhoods
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4, 5], dtype='int64', name='neighborhoods'))

  • nh_props
    PandasIndex

    PandasIndex(Index(['_color', '_name'], dtype='object', name='nh_props'))

  • x
    PandasIndex

    PandasIndex(Index([1600, 1601, 1602, 1603, 1604, 1605, 1606, 1607, 1608, 1609,
           ...
           1691, 1692, 1693, 1694, 1695, 1696, 1697, 1698, 1699, 1700],
          dtype='int64', name='x', length=101))

  • y
    PandasIndex

    PandasIndex(Index([2100, 2101, 2102, 2103, 2104, 2105, 2106, 2107, 2108, 2109,
           ...
           2191, 2192, 2193, 2194, 2195, 2196, 2197, 2198, 2199, 2200],
          dtype='int64', name='y', length=101))

• Attributes: (0)

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, issue 2.

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")

<xarray.backends.zarr.ZarrStore at 0x7f1c5e212e40>

Exporting Tables to CSV

Let’s say you want to export some tables as csvs. This can be done with pandas.

df = ds.pp.get_layer_as_df("_obs")
df.head()

	CD4_binarized	CD8_binarized	_labels	_neighborhoods	centroid-0	centroid-1
1	0.0	0.0	B	Neighborhood 1	2103.768519	1607.277778
2	0.0	1.0	T_tox	Neighborhood 1	2103.857143	1630.741071
3	1.0	1.0	T_h	Neighborhood 3	2104.837037	1668.733333
4	0.0	1.0	T_tox	Neighborhood 3	2101.750000	1677.000000
5	0.0	1.0	B	Neighborhood 3	2104.416058	1685.627737

# 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.

# 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

AnnData object with n_obs × n_vars = 56 × 5
    obs: 'CD4_binarized', 'CD8_binarized', '_labels', '_neighborhoods', 'centroid-0', 'centroid-1'
    uns: '_labels_colors', 'label_colors'
    obsm: 'spatial'
    layers: 'percentage_positive'

# writing to disk as hdf5
adata.write("tmp.h5ad")

... storing '_labels' as categorical
... storing '_neighborhoods' as categorical

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.

spatialdata_object = ds.tl.convert_to_spatialdata(expression_matrix_key="_intensity")
spatialdata_object

/home/meyerben/meyerben/.conda/envs/tmp_env_3/lib/python3.10/site-packages/dask/dataframe/__init__.py:31: FutureWarning: The legacy Dask DataFrame implementation is deprecated and will be removed in a future version. Set the configuration option `dataframe.query-planning` to `True` or None to enable the new Dask Dataframe implementation and silence this warning.
  warnings.warn(

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] (5, 101, 101)
├── Labels
│     └── 'segmentation': DataArray[yx] (101, 101)
└── Tables
      └── 'table': AnnData (56, 5)
with coordinate systems:
    ▸ 'global', with elements:
        image (Images), segmentation (Labels)

# storing as zarr file
spatialdata_object.write("tmp.zarr")

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.

ds = sp.read_from_spatialdata("../../data/spatialdata_example.zarr", image_key="raccoon")
ds

root_attr: multiscales
root_attr: omero
datasets [{'coordinateTransformations': [{'scale': [1.0, 1.0, 1.0], 'type': 'scale'}], 'path': '0'}]
/home/meyerben/meyerben/.conda/envs/tmp_env_3/lib/python3.10/site-packages/zarr/creation.py:614: UserWarning: ignoring keyword argument 'read_only'
  compressor, fill_value = _kwargs_compat(compressor, fill_value, kwargs)
resolution: 0
 - shape ('c', 'y', 'x') = (3, 768, 1024)
 - chunks =  ['3', '768', '1024']
 - dtype = uint8
root_attr: multiscales
root_attr: omero
Unsupported transform Identity , resetting coordinates for the spatialproteomics object.

<xarray.Dataset> Size: 9MB
Dimensions:        (channels: 3, y: 768, x: 1024, cells: 70, features: 2)
Coordinates:
  * channels       (channels) int64 24B 0 1 2
  * y              (y) int64 6kB 0 1 2 3 4 5 6 7 ... 761 762 763 764 765 766 767
  * x              (x) int64 8kB 0 1 2 3 4 5 6 ... 1018 1019 1020 1021 1022 1023
  * cells          (cells) int64 560B 1 2 3 4 5 6 7 8 ... 64 65 66 67 68 69 70
  * features       (features) <U10 80B 'centroid-0' 'centroid-1'
Data variables:
    _image         (channels, y, x) uint8 2MB dask.array<chunksize=(3, 768, 1024), meta=np.ndarray>
    _segmentation  (y, x) int64 6MB 0 0 0 0 0 0 0 0 ... 69 69 69 69 69 69 69 69
    _obs           (cells, features) float64 1kB 44.79 402.5 ... 736.5 890.5

xarray.Dataset

• Dimensions:
  • channels: 3
  • y: 768
  • x: 1024
  • cells: 70
  • features: 2
• Coordinates: (5)
  • channels
    (channels)
    int64
    0 1 2

    array([0, 1, 2])

  • y
    (y)
    int64
    0 1 2 3 4 5 ... 763 764 765 766 767

    array([  0,   1,   2, ..., 765, 766, 767])

  • x
    (x)
    int64
    0 1 2 3 4 ... 1020 1021 1022 1023

    array([   0,    1,    2, ..., 1021, 1022, 1023])

  • cells
    (cells)
    int64
    1 2 3 4 5 6 7 ... 65 66 67 68 69 70

    array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
           19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
           37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
           55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70])

  • features
    (features)
    <U10
    'centroid-0' 'centroid-1'

    array(['centroid-0', 'centroid-1'], dtype='<U10')

• Data variables: (3)
  • _image
    (channels, y, x)
    uint8
    dask.array<chunksize=(3, 768, 1024), meta=np.ndarray>

    transform :

    {'global': Identity }

    Array Chunk Bytes 2.25 MiB 2.25 MiB Shape (3, 768, 1024) (3, 768, 1024) Dask graph 1 chunks in 2 graph layers Data type uint8 numpy.ndarray

  • _segmentation
    (y, x)
    int64
    0 0 0 0 0 0 0 ... 69 69 69 69 69 69

    array([[ 0,  0,  0, ...,  6,  6,  6],
           [ 0,  0,  0, ...,  6,  6,  6],
           [ 0,  0,  0, ...,  6,  6,  6],
           ...,
           [54, 54, 54, ..., 69, 69, 69],
           [54, 54, 54, ..., 69, 69, 69],
           [54, 54, 54, ..., 69, 69, 69]])

  • _obs
    (cells, features)
    float64
    44.79 402.5 46.1 ... 736.5 890.5

    array([[ 44.78587164, 402.51402555],
           [ 46.09709402, 509.92752137],
           [ 46.69286234, 619.94337501],
           [ 57.15253528, 774.59581809],
           [ 38.4374736 , 864.7918251 ],
           [ 63.55431734, 964.27397468],
           [121.78519217, 418.70400063],
           [106.3700368 ,  54.14516754],
           [113.31814502, 187.80652249],
           [132.80366953, 488.87643827],
           [133.98213592, 716.02457698],
           [153.50623811, 299.81327976],
           [143.60992817, 600.04002052],
           [168.41822767, 855.63922911],
           [123.82931306, 877.34714211],
           [181.28506958, 964.84825122],
           [219.45214906,  58.44334286],
           [230.64467967, 119.52869081],
           [219.6725626 , 215.18721364],
           [245.02390979, 428.43622868],
    ...
           [565.92275762, 158.6677495 ],
           [569.96313057, 227.85944363],
           [560.5407913 , 612.03252639],
           [677.20036704,  69.88212319],
           [680.40545436, 318.10314984],
           [579.23025874, 893.34131384],
           [606.66677498, 584.26271324],
           [698.75664811, 240.51589363],
           [643.95254269, 735.43315407],
           [656.45687871, 896.77599926],
           [703.55210333, 162.2459417 ],
           [676.97454256, 659.92992575],
           [679.89467356, 412.44306841],
           [686.24289934, 358.41354238],
           [733.43646833, 686.88995521],
           [715.61062682, 559.24837415],
           [692.89371859, 863.47663317],
           [730.6795004 , 480.58863298],
           [727.75352895, 974.09770024],
           [736.47787879, 890.45787879]])

• Indexes: (5)
  • channels
    PandasIndex

    PandasIndex(Index([0, 1, 2], dtype='int64', name='channels'))

  • y
    PandasIndex

    PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
           ...
           758, 759, 760, 761, 762, 763, 764, 765, 766, 767],
          dtype='int64', name='y', length=768))

  • x
    PandasIndex

    PandasIndex(Index([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,
           ...
           1014, 1015, 1016, 1017, 1018, 1019, 1020, 1021, 1022, 1023],
          dtype='int64', name='x', length=1024))

  • cells
    PandasIndex

    PandasIndex(Index([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
           19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
           37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
           55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70],
          dtype='int64', name='cells'))

  • features
    PandasIndex

    PandasIndex(Index(['centroid-0', 'centroid-1'], dtype='object', name='features'))

• Attributes: (0)

Neighborhood Analysis Interactivity