Add ROIs¶

A region of interest is any subset of the tissue you want to treat as a unit — tumour core, invasive margin, a lymphoid aggregate. Once cells carry an ROI label it behaves like any other obs column: you can compare composition across regions, search for similar ones, or restrict an analysis to one.

There are three routes into scimappro, in decreasing order of how much they let you draw by hand.

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

Route 1: draw in OMERO, import the polygons¶

The usual route for hand-drawn regions. Draw them in OMERO, export the ROI table, and [sp.helpers.addROI_omero][scimappro.helpers.addROI_omero] labels every cell by the polygon it falls inside.

import pandas as pd

roi = pd.read_csv("rois_from_omero.csv")
adata = sp.helpers.addROI_omero(adata, roi=roi, label="ROI")
adata.obs["ROI"].value_counts()

The ROI table needs a WKT geometry column — found under geometry, polygon, or roi — and a name column, given by namingColumn.

Useful options:

• bufferRoi grows every polygon by a distance, so cells just outside a boundary are still included. Negative shrinks it.
• bufferRegions sets that per region, as {name: distance}.
• overwrite=False keeps existing labels and only fills in unlabelled cells, which lets you build up regions across several exports.
• Cells outside every polygon get 'Other'.

!!! warning "Not executed here" The demo data ships no OMERO export, so this cell is illustrative. It also needs shapely: pip install "scimappro[roi]".

Route 2: define a region by a rule¶

Often a region is definable rather than drawable — "wherever ECAD is positive" is a perfectly good tumour mask. [sp.pl.addRoiScatter][scimappro.pl.addRoiScatter] does exactly that.

In [3]:

adata = sp.pl.addRoiScatter(
    adata, marker="ECAD", threshold=0.5, roiName="tumour", label="roi"
)
adata.obs["roi"].value_counts()

adata = sp.pl.addRoiScatter( adata, marker="ECAD", threshold=0.5, roiName="tumour", label="roi" ) adata.obs["roi"].value_counts()

Out[3]:

roi
tumour    11135
Other        66
Name: count, dtype: int64

In [4]:

sp.pl.spatialScatterPlot(adata, colorBy="roi", s=2, figsize=(5, 5), fontSize=7)

sp.pl.spatialScatterPlot(adata, colorBy="roi", s=2, figsize=(5, 5), fontSize=7)

!!! note "This replaces scimap's lasso tool" scimap's addROI_scatter opened a Matplotlib window for freehand lasso selection. That cannot work in a script, a headless notebook, or CI, so scimappro defines the region by a rule instead. For genuinely freehand regions, use route 1.

For more than one marker, or for a negative one, call [sp.tl.classify][scimappro.tl.classify] directly — addRoiScatter is a thin wrapper over it.

In [5]:

adata = sp.tl.classify(
    adata,
    pos=["CD45"],
    neg=["ECAD"],
    classifyLabel="immune_zone",
    failedLabel="Other",
    label="roi_immune",
    verbose=False,
)
adata.obs["roi_immune"].value_counts()

adata = sp.tl.classify( adata, pos=["CD45"], neg=["ECAD"], classifyLabel="immune_zone", failedLabel="Other", label="roi_immune", verbose=False, ) adata.obs["roi_immune"].value_counts()

Out[5]:

roi_immune
Other          11135
immune_zone       66
Name: count, dtype: int64

Route 3: geometry¶

Nothing stops you writing the column yourself. Anything that produces one label per cell works.

In [6]:

median_x = adata.obs["X_centroid"].median()
adata.obs["half"] = ["left" if x < median_x else "right"
                     for x in adata.obs["X_centroid"]]
adata.obs["half"].value_counts()

median_x = adata.obs["X_centroid"].median() adata.obs["half"] = ["left" if x < median_x else "right" for x in adata.obs["X_centroid"]] adata.obs["half"].value_counts()

Out[6]:

half
right    5601
left     5600
Name: count, dtype: int64

What ROIs are for¶

Comparing composition between regions¶

In [7]:

pd.crosstab(adata.obs["roi"], adata.obs["phenotype"], normalize="index").round(3)

pd.crosstab(adata.obs["roi"], adata.obs["phenotype"], normalize="index").round(3)

Out[7]:

phenotype	Dendritic cells	ECAD+	Immune	NK cells	Other myeloid cells	SMA+	Treg	Unknown
roi
Other	0.000	0.000	0.0	0.348	0.000	0.000	0.000	0.652
tumour	0.078	0.639	0.0	0.001	0.217	0.046	0.019	0.000

In [8]:

sp.pl.barplot(adata, xAxis="roi", yAxis="phenotype", method="percent",
              figSize=(3.5, 4))

sp.pl.barplot(adata, xAxis="roi", yAxis="phenotype", method="percent", figSize=(3.5, 4))

Out[8]:

(<Figure size 350x400 with 1 Axes>, <Axes: xlabel='roi', ylabel='Percentage'>)

Restricting an analysis to one region¶

In [9]:

tumourOnly = adata[adata.obs["roi"] == "tumour"].copy()
tumourOnly = sp.tl.spatialDistance(tumourOnly, phenotype="phenotype", verbose=False)
sp.pl.spatialDistanceHeatmap(tumourOnly, phenotype="phenotype")

tumourOnly = adata[adata.obs["roi"] == "tumour"].copy() tumourOnly = sp.tl.spatialDistance(tumourOnly, phenotype="phenotype", verbose=False) sp.pl.spatialDistanceHeatmap(tumourOnly, phenotype="phenotype")

Searching for more regions like it¶

See Search patterns — an ROI column is exactly what [sp.tl.spatialSimilarityLookup][scimappro.tl.spatialSimilarityLookup] takes as its query.

Drawing on the raw image¶

[sp.pl.addRoiImage][scimappro.pl.addRoiImage] attaches polygons with the image recorded for provenance:

adata = sp.pl.addRoiImage(
    str(DATA / "registration" / "image.tif"), adata, roi=roi
)

And [sp.pl.image_viewer][scimappro.pl.image_viewer] opens the image in napari with the cells overlaid and a lasso tool that writes selections back into obs:

adata, viewer = sp.pl.image_viewer(
    str(DATA / "registration" / "image.tif"), adata, overlay="phenotype"
)

!!! warning "Not executed here" Both open windows, so neither runs headless. image_viewer also needs a Qt binding — pip install "scimappro[qt]".

Next¶

• Search patterns — find regions resembling this one.
• Explore cell types — compare composition between regions properly.