Unsupervised clustering¶

The other route to cell types: let the expression data group itself, then work out afterwards what each group is.

Faster to get started than prior-knowledge phenotyping, and useful when you do not yet know what is in your tissue — but the clusters depend on your parameters, and they will not match between datasets without effort.

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)

k-means¶

The simplest option. You choose k, and every cell lands in exactly one cluster.

In [3]:

adata = sp.tl.cluster(adata, method="kmeans", k=8, layer="raw", log=True,
                      label="kmeans", verbose=False)
adata.obs["kmeans"].value_counts()

adata = sp.tl.cluster(adata, method="kmeans", k=8, layer="raw", log=True, label="kmeans", verbose=False) adata.obs["kmeans"].value_counts()

Out[3]:

kmeans
4    3546
0    3450
6    1393
5    1297
2     922
3     511
1      66
7      16
Name: count, dtype: int64

Leiden¶

Graph-based, and the usual choice when you do not want to fix the number of clusters in advance. leidenResolution controls granularity — higher gives more, smaller clusters.

In [4]:

adata = sp.tl.cluster(adata, method="leiden", leidenResolution=0.5,
                      leidenNearestNeighbors=30, layer="raw", log=True,
                      label="leiden", verbose=False)
adata.obs["leiden"].value_counts()

adata = sp.tl.cluster(adata, method="leiden", leidenResolution=0.5, leidenNearestNeighbors=30, layer="raw", log=True, label="leiden", verbose=False) adata.obs["leiden"].value_counts()

Out[4]:

leiden
0    2273
1    2078
2    1654
3    1290
4    1166
5    1150
6     700
7     477
8     413
Name: count, dtype: int64

dbscan is also available, via dbscanEps and dbscanMinSamples. It finds arbitrarily-shaped clusters and labels outliers rather than forcing every cell into a group.

What is each cluster?¶

This is the part clustering does not answer for you. The heatmap of mean expression per cluster is where you read it off.

In [5]:

sp.pl.heatmap(adata, groupBy="leiden", layer="raw", standardScale="column",
              showPrevalence=True, figsize=(7, 5))

sp.pl.heatmap(adata, groupBy="leiden", layer="raw", standardScale="column", showPrevalence=True, figsize=(7, 5))

standardScale="column" z-scores each marker across clusters, so a marker that is dim everywhere still shows its relative pattern. Drop it to see absolute levels instead.

On a UMAP¶

In [6]:

adata = sp.tl.umap(adata, layer="raw", log=True, verbose=False)
sp.pl.clusterPlots(adata, groupBy="leiden", size=1)

adata = sp.tl.umap(adata, layer="raw", log=True, verbose=False) sp.pl.clusterPlots(adata, groupBy="leiden", size=1)

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

Colour the same embedding by a marker to check a cluster's identity:

In [7]:

sp.pl.umap(adata, color=["CD45", "ECAD", "SMA"], layer="raw", log=True,
           s=1, ncols=3, figsize=(3.5, 3))

sp.pl.umap(adata, color=["CD45", "ECAD", "SMA"], layer="raw", log=True, s=1, ncols=3, figsize=(3.5, 3))

In the tissue¶

In [8]:

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

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

Naming the clusters¶

Once you know what a cluster is, give it a name. [sp.tl.rename][scimappro.tl.rename] writes a new column, leaving the numbers in place.

In [9]:

adata = sp.tl.rename(
    adata,
    rename={"Tumour": ["0"], "Stroma": ["1"], "Immune": ["2"]},
    fromColumn="leiden",
    toColumn="cluster_named",
    verbose=False,
)
adata.obs["cluster_named"].value_counts()

adata = sp.tl.rename( adata, rename={"Tumour": ["0"], "Stroma": ["1"], "Immune": ["2"]}, fromColumn="leiden", toColumn="cluster_named", verbose=False, ) adata.obs["cluster_named"].value_counts()

Out[9]:

cluster_named
Tumour    2273
Stroma    2078
Immune    1654
3         1290
4         1166
5         1150
6          700
7          477
8          413
Name: count, dtype: int64

Sub-clustering¶

When one cluster is obviously heterogeneous, split just that one and leave the rest alone. collapseLabels=True keeps the other cells' original labels, so the output stays a complete column.

In [10]:

adata = sp.tl.cluster(
    adata,
    method="kmeans",
    k=3,
    subCluster="leiden",
    subClusterGroup=["0"],
    collapseLabels=True,
    label="leiden_sub",
    verbose=False,
)
adata.obs["leiden_sub"].value_counts()

adata = sp.tl.cluster( adata, method="kmeans", k=3, subCluster="leiden", subClusterGroup=["0"], collapseLabels=True, label="leiden_sub", verbose=False, ) adata.obs["leiden_sub"].value_counts()

Out[10]:

leiden_sub
0-0    2191
1      2078
2      1654
3      1290
4      1166
5      1150
6       700
7       477
8       413
0-1      66
0-2      16
Name: count, dtype: int64

Clustering vs. phenotyping¶

	Clustering	Phenotyping
Needs prior knowledge	no	yes — you write the gating table
Reproducible across datasets	not without work	yes, the same table runs everywhere
Finds unexpected populations	yes	only what you asked for
Sensitive to rare cell types	poorly	well
Depends on parameters	strongly (k, resolution)	on the gates

In practice they complement each other: cluster to see what is there, then write a workflow that captures it reproducibly.

Next¶

• Prior-knowledge phenotyping — the reproducible route.
• Latent motifs — clustering neighbourhoods rather than cells.