Explore cell types¶

Once cells are labelled, the first questions are compositional: what is this tissue made of, and does that differ between samples?

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

Composition as bars¶

[sp.pl.barplot][scimappro.pl.barplot] cross-tabulates a cell type column within an image column. method="percent" normalises each bar to 100%, which is what you want whenever images differ in cell count.

In [3]:

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

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

Out[3]:

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

Raw counts instead:

In [4]:

sp.pl.barplot(adata, xAxis="imageid", yAxis="phenotype", method="absolute",
              figSize=(4, 4))

sp.pl.barplot(adata, xAxis="imageid", yAxis="phenotype", method="absolute", figSize=(4, 4))

Out[4]:

(<Figure size 400x400 with 1 Axes>, <Axes: xlabel='imageid', ylabel='Count'>)

Composition as pies¶

In [5]:

sp.pl.pie(adata, phenotype="phenotype", groupBy="imageid")

sp.pl.pie(adata, phenotype="phenotype", groupBy="imageid")

Marker expression per cell type¶

The heatmap is the check that the phenotypes express what they should. On rescaled data, anything above 0.5 is positive.

In [6]:

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)

Comparing groups¶

Everything below needs at least two groups. The demo data is a single image, so we invent two by splitting the slide down the middle.

In [7]:

# The demo data is a single image, and comparisons need at least two groups.
# Splitting it down the middle gives two synthetic "samples" to compare.
# These are NOT biological groups — they exist only to make the code runnable.
median_x = adata.obs["X_centroid"].median()
adata.obs["region"] = ["left" if x < median_x else "right"
                       for x in adata.obs["X_centroid"]]
adata.obs["region"].value_counts()

# The demo data is a single image, and comparisons need at least two groups. # Splitting it down the middle gives two synthetic "samples" to compare. # These are NOT biological groups — they exist only to make the code runnable. median_x = adata.obs["X_centroid"].median() adata.obs["region"] = ["left" if x < median_x else "right" for x in adata.obs["X_centroid"]] adata.obs["region"].value_counts()

Out[7]:

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

Composition side by side¶

In [8]:

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

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

Out[8]:

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

Fold change¶

[sp.tl.foldChange][scimappro.tl.foldChange] divides each group's normalised composition by a reference group's, and runs a Fisher exact test on the counts.

In [9]:

adata = sp.tl.foldChange(adata, fromGroup="left", imageId="region",
                         phenotype="phenotype", normalize=True, verbose=False)
adata.uns["foldchange_fc"]

adata = sp.tl.foldChange(adata, fromGroup="left", imageId="region", phenotype="phenotype", normalize=True, verbose=False) adata.uns["foldchange_fc"]

Out[9]:

phenotype	Dendritic cells	ECAD+	Immune	NK cells	Other myeloid cells	SMA+	Treg	Unknown
region
right	1.560552	1.037053	inf	0.166637	0.853487	0.927865	0.700662	0.0

In [10]:

adata.uns["foldchange_pval"]

adata.uns["foldchange_pval"]

Out[10]:

	Dendritic cells	ECAD+	Immune	NK cells	Other myeloid cells	SMA+	Treg	Unknown
region
right	2.300900e-11	0.011353	1.0	0.000013	0.000012	0.389165	0.009039	1.289840e-14

In [11]:

sp.pl.foldChange(adata, method="heatmap")

sp.pl.foldChange(adata, method="heatmap")

/Users/aj/Partners HealthCare Dropbox/Ajit Nirmal/nirmal lab/softwares/dev/scimappro/scimappro/pl/foldChange.py:153: RuntimeWarning: divide by zero encountered in log2
  matrix = np.log2(plotData.to_numpy(dtype=float)) if log else plotData.to_numpy(dtype=float)

Do samples resemble each other?¶

[sp.pl.groupCorrelation][scimappro.pl.groupCorrelation] correlates groups by their composition — with many samples, this is how you spot which ones cluster together.

In [12]:

sp.pl.groupCorrelation(adata, groupBy="region", condition="phenotype",
                       normalize=True, overlayValues=True, figsize=(3, 3))

sp.pl.groupCorrelation(adata, groupBy="region", condition="phenotype", normalize=True, overlayValues=True, figsize=(3, 3))

Where the cell types are¶

In [13]:

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)

A Voronoi tessellation fills the space between cells, which makes tissue compartments easier to read than a point scatter. It is slower — one polygon per cell — so work one image at a time.

In [14]:

sp.pl.voronoi(adata, colorBy="phenotype", voronoiLineWidth=0.05, voronoiAlpha=0.8)

sp.pl.voronoi(adata, colorBy="phenotype", voronoiLineWidth=0.05, voronoiAlpha=0.8)

Next¶

• Distance measurement — how far apart the cell types sit.
• Co-occurrence analysis — which pairs sit together more than chance.