SpatialData workflow¶

scimappro speaks SpatialData natively: the same data argument that takes an AnnData takes a SpatialData object or a path to a .zarr store, and every function behaves identically.

There is no scimap equivalent to any of this — see the SpatialData guide for the full contract.

In [1]:

from pathlib import Path

import matplotlib.pyplot as plt
import pandas as pd

import scimappro as sp

# The demo data is not shipped with the package (it is 247 MB). Point DATA at
# wherever you unpacked it; this looks in the repository root first, then in the
# path the docs are built from.
DATA = next(
    path for path in (Path("example_data"), Path("../../../example_data"))
    if path.exists()
)
DATA

from pathlib import Path import matplotlib.pyplot as plt import pandas as pd import scimappro as sp # The demo data is not shipped with the package (it is 247 MB). Point DATA at # wherever you unpacked it; this looks in the repository root first, then in the # path the docs are built from. DATA = next( path for path in (Path("example_data"), Path("../../../example_data")) if path.exists() ) DATA

Out[1]:

PosixPath('../../../example_data')

In [2]:

import anndata as ad

adata = ad.read_h5ad(DATA / "adata_scimap.h5ad")
adata = sp.pp.rescale(adata, gate=str(DATA / "manual_gates.csv"), verbose=False)
adata = sp.tl.phenotype(adata, phenotype=str(DATA / "phenotype_workflow.csv"), verbose=False)
adata.obs["phenotype"].value_counts()

import anndata as ad adata = ad.read_h5ad(DATA / "adata_scimap.h5ad") adata = sp.pp.rescale(adata, gate=str(DATA / "manual_gates.csv"), verbose=False) adata = sp.tl.phenotype(adata, phenotype=str(DATA / "phenotype_workflow.csv"), verbose=False) adata.obs["phenotype"].value_counts()

Out[2]:

phenotype
ECAD+                  7112
Other myeloid cells    2419
Dendritic cells         863
SMA+                    509
Treg                    216
Unknown                  46
NK cells                 35
Immune                    1
Name: count, dtype: int64

Converting a cell table¶

[sp.pp.toSpatialData][scimappro.pp.toSpatialData] builds one circles element per image, holding that image's cell centroids, with the whole table attached to those elements.

In [3]:

sdata = sp.pp.toSpatialData(adata, verbose=False)
sdata

sdata = sp.pp.toSpatialData(adata, verbose=False) sdata

Out[3]:

SpatialData object
├── Shapes
│     └── 'quant': GeoDataFrame shape: (11201, 2) (2D shapes)
└── Tables
      └── 'table': AnnData (11201, 9)
with coordinate systems:
    ▸ 'global', with elements:
        quant (Shapes)

The radius comes from obs['Area'] by default, as sqrt(area / pi). Pass a number or another column name to override it.

In [4]:

sdata["quant"].head()

sdata["quant"].head()

Out[4]:

	geometry	radius
1	POINT (1767.692 257.291)	6.102643
2	POINT (1107.174 665.87)	5.411516
3	POINT (1116.414 671.069)	4.296740
4	POINT (982.729 677.03)	9.253397
5	POINT (1141.071 680.125)	11.396071

Every obs column is preserved, so nothing is lost even though the geometry is two-dimensional.

In [5]:

sdata.tables["table"]

sdata.tables["table"]

Out[5]:

AnnData object with n_obs × n_vars = 11201 × 9
    obs: 'X_centroid', 'Y_centroid', 'Area', 'MajorAxisLength', 'MinorAxisLength', 'Eccentricity', 'Solidity', 'Extent', 'Orientation', 'CellID', 'imageid', 'phenotype'
    var: 'index'
    uns: 'all_markers', 'gates', 'spatialdata_attrs'
    layers: 'log', None (.X)

Running an analysis on it¶

Exactly the same call as for an AnnData. sdataTable names the table to work on and is always the last named parameter.

In [6]:

sdata = sp.tl.spatialDistance(sdata, phenotype="phenotype", sdataTable="table",
                              verbose=False)
sdata.tables["table"].uns["spatial_distance"].head()

sdata = sp.tl.spatialDistance(sdata, phenotype="phenotype", sdataTable="table", verbose=False) sdata.tables["table"].uns["spatial_distance"].head()

Out[6]:

	Dendritic cells	ECAD+	Immune	NK cells	Other myeloid cells	SMA+	Treg	Unknown
obs
quant_1	554.265017	508.856310	0.000000	2053.512851	532.311179	505.979757	575.266490	2086.992073
quant_2	114.340629	10.602318	776.673240	1287.763826	25.455242	27.967658	225.297772	1311.182079
quant_3	122.911429	10.602318	771.606273	1292.597225	24.200678	34.714545	226.977790	1315.778661
quant_4	11.936186	13.731018	890.139838	1178.924606	39.819156	36.240050	168.802786	1204.626638
quant_5	135.256707	15.790519	755.938553	1308.403525	40.006555	58.347658	215.734114	1331.056574

In [7]:

sdata = sp.tl.spatialCooccurrence(sdata, phenotype="phenotype", permutation=200,
                                  sdataTable="table", verbose=False)
sp.pl.spatialCooccurrence(sdata, sdataTable="table")

sdata = sp.tl.spatialCooccurrence(sdata, phenotype="phenotype", permutation=200, sdataTable="table", verbose=False) sp.pl.spatialCooccurrence(sdata, sdataTable="table")

With one table in the store, sdataTable can be omitted entirely:

In [8]:

list(sdata.tables)

list(sdata.tables)

Out[8]:

['table']

In [9]:

result = sp.tl.spatialAggregate(sdata, phenotype="phenotype", radius=30, verbose=False)
result.tables["table"].obs["spatialAggregate"].value_counts()

result = sp.tl.spatialAggregate(sdata, phenotype="phenotype", radius=30, verbose=False) result.tables["table"].obs["spatialAggregate"].value_counts()

Out[9]:

spatialAggregate
ECAD+                  6643
non-significant        2812
Other myeloid cells    1106
Dendritic cells         527
SMA+                     57
Unknown                  34
NK cells                 18
Treg                      4
Name: count, dtype: int64

Writing a store¶

outputDir writes <inputFilename>.zarr there and returns None.

In [10]:

outputDir = Path("tutorial_output")
outputDir.mkdir(exist_ok=True)

sp.pp.toSpatialData(adata, outputDir=str(outputDir), verbose=False)
sorted(p.name for p in outputDir.iterdir())

outputDir = Path("tutorial_output") outputDir.mkdir(exist_ok=True) sp.pp.toSpatialData(adata, outputDir=str(outputDir), verbose=False) sorted(p.name for p in outputDir.iterdir())

Out[10]:

['adata_scimap.zarr']

Reading it back — and working from the path¶

A .zarr path is a valid data argument, so there is nothing to load first.

In [11]:

import spatialdata as sd

storePath = outputDir / "adata_scimap.zarr"
reloaded = sd.read_zarr(storePath)
reloaded

import spatialdata as sd storePath = outputDir / "adata_scimap.zarr" reloaded = sd.read_zarr(storePath) reloaded

Out[11]:

SpatialData object, with associated Zarr store: /Users/aj/Partners HealthCare Dropbox/Ajit Nirmal/nirmal lab/softwares/dev/scimappro/docs/tutorials/nbs/tutorial_output/adata_scimap.zarr
├── Shapes
│     └── 'quant': GeoDataFrame shape: (11201, 2) (2D shapes)
└── Tables
      └── 'table': AnnData (11201, 9)
with coordinate systems:
    ▸ 'global', with elements:
        quant (Shapes)

In [12]:

scored = sp.tl.spatialProximityScore(
    str(storePath), proximity=["Treg", "ECAD+"], radius=20, verbose=False
)
scored.tables["table"].uns["spatial_pscore"]

scored = sp.tl.spatialProximityScore( str(storePath), proximity=["Treg", "ECAD+"], radius=20, verbose=False ) scored.tables["table"].uns["spatial_pscore"]

Out[12]:

spatial_pscore	Treg_ECAD+	All Cells	Proximity Volume	Celltype of interest	Proximity Density
imageid
quant	565	11201	0.050442	7328	0.049809

Updating a store in place¶

Point outputDir at the store the object came from and only that one table is rewritten, rather than the whole store.

In [13]:

sdataOnDisk = sd.read_zarr(storePath)
sp.tl.spatialDistance(sdataOnDisk, phenotype="phenotype", outputDir=str(storePath),
                      verbose=False)

check = sd.read_zarr(storePath)
"spatial_distance" in check.tables["table"].uns

sdataOnDisk = sd.read_zarr(storePath) sp.tl.spatialDistance(sdataOnDisk, phenotype="phenotype", outputDir=str(storePath), verbose=False) check = sd.read_zarr(storePath) "spatial_distance" in check.tables["table"].uns

Out[13]:

True

Coordinates from the elements¶

A SpatialData table often has no X_centroid / Y_centroid / imageid columns — that information lives in the elements the table annotates. scimappro derives them, preserving row order, and stores them on the table. Nothing is required from you.

Watch it happen with verbose=True:

In [14]:

bare = sd.read_zarr(storePath)
bare.tables["table"].obs = bare.tables["table"].obs.drop(
    columns=["X_centroid", "Y_centroid"]
)
derived = sp.tl.spatialDistance(bare, phenotype="phenotype", verbose=True)
derived.tables["table"].obs[["X_centroid", "Y_centroid"]].head()

bare = sd.read_zarr(storePath) bare.tables["table"].obs = bare.tables["table"].obs.drop( columns=["X_centroid", "Y_centroid"] ) derived = sp.tl.spatialDistance(bare, phenotype="phenotype", verbose=True) derived.tables["table"].obs[["X_centroid", "Y_centroid"]].head()

Using SpatialData table 'table'.

Derived obs columns ['X_centroid', 'Y_centroid'] from the SpatialData elements.
Processing Image: quant

  0%|          | 0/8 [00:00<?, ?it/s]

100%|██████████| 8/8 [00:00<00:00, 12368.02it/s]

Out[14]:

	X_centroid	Y_centroid
quant_1	1767.692308	257.290598
quant_2	1107.173913	665.869565
quant_3	1116.413793	671.068966
quant_4	982.728625	677.029740
quant_5	1141.071078	680.125000

Things to know¶

• streamData does not apply. SpatialData tables are already in memory; passing it warns and continues. See Streaming.
• Functions that replace the table — [tl.classify][scimappro.tl.classify], [tl.rename][scimappro.tl.rename], [pp.dropFeatures][scimappro.pp.dropFeatures], [helpers.addROI_omero][scimappro.helpers.addROI_omero] — put the new table back in the container. Keep its region_key and instance_key obs columns intact or the write-back will not validate.
• Several tables without sdataTable raises a ValueError listing them.

Next¶

• The SpatialData guide — the full contract.
• Streaming large files — the other way to handle data that will not fit.