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

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.

# 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"]

%load_ext rpy2.ipython

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

# 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"]

# send the data to R, so that we can use it directly
%R -i intensity_df
%R -i spatial_df

%%R
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)

    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’:

    intersect, t


Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

Warning: Data is of class matrix. Coercing to dgCMatrix.

%%R
# 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)
)

%%R
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)

Learn more about the underlying theory at https://ggplot2-book.org/

Neighborhood Analysis Downstream Analysis