Prior-knowledge phenotyping¶

Cell typing driven by a table you write, rather than by clustering. You state once which markers define which cell type, and the same table then runs across every sample identically.

It takes more thought up front than clustering, and it is more sensitive and far more reproducible in return — nothing depends on a resolution parameter or on which clusters happen to form today.

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

import anndata as ad adata = ad.read_h5ad(DATA / "adata_scimap.h5ad") adata

Out[2]:

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

Step 1: find the gates¶

A gate is the intensity above which a marker counts as positive. The two plotting functions below are how you pick one by eye.

Single-marker distributions — a bimodal histogram usually has its gate in the valley:

In [3]:

sp.pl.distPlot(adata, layer="log", markers=["CD45", "ECAD", "FOXP3"], ncols=3,
               figsize=(3, 2.5), fontSize=8)

sp.pl.distPlot(adata, layer="log", markers=["CD45", "ECAD", "FOXP3"], ncols=3, figsize=(3, 2.5), fontSize=8)

And the biaxial view, which is often clearer when one marker is dim:

In [4]:

sp.pl.densityPlot2D(adata, markerA="SMA", markerB="CD45", layer="log", figsize=(3, 3))

sp.pl.densityPlot2D(adata, markerA="SMA", markerB="CD45", layer="log", figsize=(3, 3))

!!! note "scimap's interactive gating is not in scimappro"

`sm.pl.gate_finder` and `sm.pl.napariGater` opened a napari window to pick
gates by eye against the raw image. They are not ported. Use the plots above
plus a `manual_gates.csv`, or let
[`sp.pp.rescale`][scimappro.pp.rescale] fit a Gaussian Mixture Model per
marker when you give it no gates at all.

Write the gates you settle on into a CSV. The first column must be called markers; every other column is named after an image id in obs['imageid'], so one file can carry different gates for different slides.

In [5]:

gates = pd.read_csv(DATA / "manual_gates.csv")
gates

gates = pd.read_csv(DATA / "manual_gates.csv") gates

Out[5]:

	markers	quant
0	DNA_6	6
1	ELANE	6
2	CD57	6
3	CD45	6
4	DNA_7	6
5	CD11B	6
6	SMA	6
7	CD16	6
8	DNA_8	6
9	ECAD	6
10	FOXP3	6
11	NCAM	6

Step 2: rescale¶

[sp.pp.rescale][scimappro.pp.rescale] maps each marker's gate to exactly 0.5, so one threshold works for every marker in every image. Markers absent from the CSV get a gate from a Gaussian Mixture Model automatically.

In [6]:

adata = sp.pp.rescale(adata, gate=str(DATA / "manual_gates.csv"), verbose=False)
adata.uns["gates"]

adata = sp.pp.rescale(adata, gate=str(DATA / "manual_gates.csv"), verbose=False) adata.uns["gates"]

Out[6]:

	quant
ELANE	6
CD57	6
CD45	6
CD11B	6
SMA	6
CD16	6
ECAD	6
FOXP3	6
NCAM	6

Check the result: every marker should now straddle 0.5.

In [7]:

sp.pl.distPlot(adata, markers=["CD45", "ECAD", "FOXP3"], ncols=3, vline=0.5,
               figsize=(3, 2.5), fontSize=8)

sp.pl.distPlot(adata, markers=["CD45", "ECAD", "FOXP3"], ncols=3, vline=0.5, figsize=(3, 2.5), fontSize=8)

Step 3: write the workflow¶

The phenotype workflow is a hierarchy: broad classes first, refined into specific cell types below.

In [8]:

workflow = pd.read_csv(DATA / "phenotype_workflow.csv")
workflow.fillna("")

workflow = pd.read_csv(DATA / "phenotype_workflow.csv") workflow.fillna("")

Out[8]:

	Unnamed: 0	Unnamed: 1	ELANE	CD57	CD45	CD11B	SMA	CD16	ECAD	FOXP3	NCAM
0	all	ECAD+							pos
1	all	Immune			pos
2	all	SMA+					pos
3	Immune	NK cells		allpos		neg		allpos
4	Immune	Other myeloid cells				pos
5	Immune	Treg								pos
6	Other myeloid cells	Dendritic cells		allneg				allneg

Reading the table:

1. Column 1 is the parent group. all means every cell; anything else means "only cells already labelled that at a previous level". So Immune is resolved first, and NK cells then refines the cells that got it.
2. Column 2 is the cell type to assign.
3. The rest are markers, holding one of six keywords:

Keyword	Meaning
pos	positive for this marker
neg	negative for this marker
anypos	positive for any of the markers marked anypos in this row
anyneg	negative for any of them
allpos	positive for all of them
allneg	negative for all of them

Step 4: run it¶

In [9]:

adata = sp.tl.phenotype(adata, phenotype=str(DATA / "phenotype_workflow.csv"),
                        label="phenotype", verbose=False)
adata.obs["phenotype"].value_counts()

adata = sp.tl.phenotype(adata, phenotype=str(DATA / "phenotype_workflow.csv"), label="phenotype", verbose=False) adata.obs["phenotype"].value_counts()

Out[9]:

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

Unknown means a cell matched nothing at any level. A large Unknown fraction is a signal that the gates are too strict or that the workflow is missing a population.

Note that no cell is reported as Immune-rest: cells that match a parent group but none of its children collapse back to the parent's own name.

Visualising the result¶

The heatmap is the first check. Because the data is rescaled, values above 0.5 are positive — so each phenotype should light up exactly the markers that define it.

In [10]:

sp.pl.heatmap(adata, groupBy="phenotype", standardScale=None, clusterRows=False,
              clusterColumns=False, showPrevalence=True, figsize=(6, 4), vmin=0, vmax=1)

sp.pl.heatmap(adata, groupBy="phenotype", standardScale=None, clusterRows=False, clusterColumns=False, showPrevalence=True, figsize=(6, 4), vmin=0, vmax=1)

Then where the cell types actually sit in the tissue:

In [11]:

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

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

And on a UMAP, which shows whether the phenotypes are separable in expression space as well as recognisable in the tissue:

In [12]:

adata = sp.tl.umap(adata, verbose=False)
sp.pl.umap(adata, color="phenotype", s=1, figsize=(5, 4))

adata = sp.tl.umap(adata, verbose=False) sp.pl.umap(adata, color="phenotype", s=1, figsize=(5, 4))

/Users/aj/Partners HealthCare Dropbox/Ajit Nirmal/nirmal lab/softwares/dev/scimappro/.venv/lib/python3.12/site-packages/umap/umap_.py:1952: UserWarning: n_jobs value 1 overridden to 1 by setting random_state. Use no seed for parallelism.
  warn(

Suppressing rare calls¶

A cell type called from a handful of cells is usually noise. phenoThresholdPercent relabels anything below a percentage of cells as Unknown.

In [13]:

strict = sp.tl.phenotype(
    adata.copy(),
    phenotype=str(DATA / "phenotype_workflow.csv"),
    phenoThresholdPercent=1,
    label="phenotype_strict",
    verbose=False,
)
strict.obs["phenotype_strict"].value_counts()

strict = sp.tl.phenotype( adata.copy(), phenotype=str(DATA / "phenotype_workflow.csv"), phenoThresholdPercent=1, label="phenotype_strict", verbose=False, ) strict.obs["phenotype_strict"].value_counts()

Out[13]:

phenotype_strict
ECAD+                  7112
Other myeloid cells    2419
Dendritic cells         863
SMA+                    509
Treg                    216
Unknown                  82
Name: count, dtype: int64

Refining one population further¶

[sp.tl.classify][scimappro.tl.classify] applies a simple rule to a subset of the cells you already labelled, leaving everything else untouched.

In [14]:

adata = sp.tl.classify(
    adata,
    pos=["CD16"],
    phenotype="phenotype",
    subclassifyPhenotype=["Other myeloid cells"],
    classifyLabel="CD16+",
    collapseFailed=True,
    showPhenotypeLabel=True,
    label="phenotype_refined",
    verbose=False,
)
adata.obs["phenotype_refined"].value_counts()

adata = sp.tl.classify( adata, pos=["CD16"], phenotype="phenotype", subclassifyPhenotype=["Other myeloid cells"], classifyLabel="CD16+", collapseFailed=True, showPhenotypeLabel=True, label="phenotype_refined", verbose=False, ) adata.obs["phenotype_refined"].value_counts()

Out[14]:

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

Overlaying on the raw image¶

[sp.pl.image_viewer][scimappro.pl.image_viewer] opens the OME-TIFF in napari with the cells on top, which is the real test of a phenotype call.

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

!!! warning "Not executed here" image_viewer opens a window, so it cannot run in the executed notebooks on this site or in any headless environment. It also needs a Qt binding — pip install "scimappro[qt]".

Next¶

• Explore cell types — composition, proportions, and differences between samples.
• Distance measurement — the first spatial statistic.
• Unsupervised clustering — the other route to cell types.