Practical 1: Cell segmentation with Cellpose¶
ELIXIR Spatial Omics 2026
In this practical we segment cells in a Xenium morphology image using Cellpose, assign transcripts to the segmented cells, and compare our segmentation against the segmentation provided by 10x Genomics.
The data is the cropped Xenium mouse brain dataset produced in Practical 0
(Xenium_V1_FFPE_TgCRND8_17_9_months_cropped.zarr). It is stored as a
SpatialData object, which bundles the
images, the 10x segmentation masks, the transcripts and the cell-by-gene
table into a single object with a shared coordinate system.
Adpated from Practical 1: Segmentation of spatial transcriptomics data using cellpose Date: 2025-01-22
1. Setup¶
In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import anndata as ad
import spatialdata as sd
import spatialdata_plot
from spatialdata import get_extent
from spatialdata.transformations import get_transformation, Scale, Sequence, set_transformation
from spatialdata.models import PointsModel, Labels2DModel
from cellpose import models, plot
from skimage.segmentation import expand_labels
from skimage.segmentation import relabel_sequential
from skimage.measure import regionprops_table
from skimage.transform import resize
from sklearn.metrics import jaccard_score
sns.set_theme(style="whitegrid")
%matplotlib inline
Load the data¶
In [2]:
from pathlib import Path
ZARR_PATH = Path("../data/Practical1_cropped.zarr")
sdata = sd.read_zarr(ZARR_PATH)
sdata
no parent found for <ome_zarr.reader.Label object at 0x7f584d404fb0>: None no parent found for <ome_zarr.reader.Label object at 0x7f584d417170>: None
Out[2]:
SpatialData object, with associated Zarr store: /home/jovyan/work/data/Practical1_cropped.zarr
├── Images
│ ├── 'morphology_focus': DataTree[cyx] (1, 5559, 7009), (1, 2780, 3505), (1, 1390, 1752), (1, 695, 876), (1, 347, 438)
│ └── 'morphology_mip': DataTree[cyx] (1, 5559, 7009), (1, 2780, 3505), (1, 1390, 1752), (1, 695, 876), (1, 347, 438)
├── Labels
│ ├── 'cell_labels': DataTree[yx] (5559, 7009), (2780, 3505), (1390, 1752), (695, 876), (347, 438)
│ └── 'nucleus_labels': DataTree[yx] (5559, 7009), (2780, 3505), (1390, 1752), (695, 876), (347, 438)
├── Points
│ └── 'transcripts': DataFrame with shape: (<Delayed>, 10) (3D points)
├── Shapes
│ ├── 'cell_boundaries': GeoDataFrame shape: (5569, 1) (2D shapes)
│ └── 'nucleus_boundaries': GeoDataFrame shape: (5422, 1) (2D shapes)
└── Tables
└── 'table': AnnData (5568, 347)
with coordinate systems:
▸ 'global', with elements:
morphology_focus (Images), morphology_mip (Images), cell_labels (Labels), nucleus_labels (Labels), transcripts (Points), cell_boundaries (Shapes), nucleus_boundaries (Shapes)
The object contains:
| Element | What it is |
|---|---|
images/morphology_mip |
maximum-intensity projection of the morphology stain |
images/morphology_focus |
best-focus plane of the same stain |
labels/cell_labels |
10x whole-cell segmentation, as a label image |
labels/nucleus_labels |
10x nuclear segmentation, as a label image |
points/transcripts |
every detected transcript, with its coordinates |
shapes/cell_boundaries |
10x cell outlines, as polygons |
tables/table |
the cell-by-gene count matrix (AnnData) |
Note that the images are stored as multiscale pyramids — several
progressively downsampled copies of the same image (scale0 is full
resolution, scale1 is half, and so on). This lets us work at a coarser
resolution when we want speed, without loading a different file.
2. Choosing a resolution¶
In [3]:
img_multiscale = sdata.images["morphology_mip"]
In [4]:
img_multiscale
Out[4]:
<xarray.DataTree>
Group: /
├── Group: /scale0
│ Dimensions: (c: 1, y: 5559, x: 7009)
│ Coordinates:
│ * c (c) int64 8B 0
│ * y (y) float64 44kB 0.5 1.5 2.5 3.5 ... 5.556e+03 5.558e+03 5.558e+03
│ * x (x) float64 56kB 0.5 1.5 2.5 3.5 ... 7.006e+03 7.008e+03 7.008e+03
│ Data variables:
│ image (c, y, x) uint16 78MB dask.array<chunksize=(1, 4096, 4096), meta=np.ndarray>
├── Group: /scale1
│ Dimensions: (c: 1, y: 2780, x: 3505)
│ Coordinates:
│ * c (c) int64 8B 0
│ * y (y) float64 22kB 0.9998 2.999 4.999 ... 5.556e+03 5.558e+03
│ * x (x) float64 28kB 0.9999 3.0 4.999 ... 7.004e+03 7.006e+03 7.008e+03
│ Data variables:
│ image (c, y, x) uint16 19MB dask.array<chunksize=(1, 2378, 2358), meta=np.ndarray>
├── Group: /scale2
│ Dimensions: (c: 1, y: 1390, x: 1752)
│ Coordinates:
│ * c (c) int64 8B 0
│ * y (y) float64 11kB 2.0 5.999 9.998 ... 5.549e+03 5.553e+03 5.557e+03
│ * x (x) float64 14kB 2.0 6.001 10.0 ... 6.999e+03 7.003e+03 7.007e+03
│ Data variables:
│ image (c, y, x) uint16 5MB dask.array<chunksize=(1, 1390, 1752), meta=np.ndarray>
├── Group: /scale3
│ Dimensions: (c: 1, y: 695, x: 876)
│ Coordinates:
│ * c (c) int64 8B 0
│ * y (y) float64 6kB 3.999 12.0 20.0 ... 5.539e+03 5.547e+03 5.555e+03
│ * x (x) float64 7kB 4.001 12.0 20.0 ... 6.989e+03 6.997e+03 7.005e+03
│ Data variables:
│ image (c, y, x) uint16 1MB dask.array<chunksize=(1, 695, 876), meta=np.ndarray>
└── Group: /scale4
Dimensions: (c: 1, y: 347, x: 438)
Coordinates:
* c (c) int64 8B 0
* y (y) float64 3kB 8.01 24.03 40.05 ... 5.519e+03 5.535e+03 5.551e+03
* x (x) float64 4kB 8.001 24.0 40.01 ... 6.969e+03 6.985e+03 7.001e+03
Data variables:
image (c, y, x) uint16 304kB dask.array<chunksize=(1, 347, 438), meta=np.ndarray>In [5]:
SCALE_LEVEL = 1
SCALE_KEY = f"scale{SCALE_LEVEL}"
# Every pyramid level is stored as a small dataset with one variable, "image".
img_multiscale = sdata.images["morphology_mip"]
image_channel = img_multiscale[SCALE_KEY]["image"].data[0].compute()
image_channel = image_channel.astype(np.uint16)
full_shape = img_multiscale["scale0"]["image"].shape[-2:]
print(f"Full resolution : {full_shape[0]} x {full_shape[1]}")
print(f"Working at {SCALE_KEY} : {image_channel.shape[0]} x {image_channel.shape[1]}")
Full resolution : 5559 x 7009 Working at scale1 : 2780 x 3505
Pixel size¶
The transcripts are stored in microns, but masks is indexed in pixels, so we
need the conversion factor.
Xenium images have a pixel size of 0.2125 µm at full resolution. We are working at a downsampled pyramid level, so the effective pixel size is larger by exactly the downsampling factor.
In [6]:
PIXEL_SIZE_UM = 0.2125 # Xenium pixel size information is from "experiment.xenium"
shrink = img_multiscale["scale0"]["image"].shape[-1] / image_channel.shape[-1]
x_scale = y_scale = 1 / (PIXEL_SIZE_UM * shrink)
print(f"{x_scale:.4f} pixels per micron (scale{SCALE_LEVEL}, {shrink:.0f}x downsampled)")
2.3533 pixels per micron (scale1, 2x downsampled)
Look at the image before segmenting it¶
In [7]:
vmin, vmax = np.percentile(image_channel, [1, 99.5])
plt.figure(figsize=(10, 8))
plt.imshow(image_channel, cmap="gray", vmin=vmin, vmax=vmax)
plt.title("Morphology image (MIP)")
plt.axis("off")
plt.show()
3. Segmentation with Cellpose¶
Cellpose is a deep learning-based tool commonly used for cell segmentation, especially on fluorescence images and it comes with pre-trained models.
There are many different pretrained models available, and you can do a quick test to check which model fits your data/cell types.
In [8]:
print(models.MODEL_NAMES)
['cyto3', 'nuclei', 'cyto2_cp3', 'tissuenet_cp3', 'livecell_cp3', 'yeast_PhC_cp3', 'yeast_BF_cp3', 'bact_phase_cp3', 'bact_fluor_cp3', 'deepbacs_cp3', 'cyto2', 'cyto', 'CPx', 'transformer_cp3', 'neurips_cellpose_default', 'neurips_cellpose_transformer', 'neurips_grayscale_cyto2', 'CP', 'CPx', 'TN1', 'TN2', 'TN3', 'LC1', 'LC2', 'LC3', 'LC4']
We initialize a model instance using models.CellposeModel().
We also set model_type='nuclei', indicating that the model should use Cellpose’s pre-trained “nuclei” model, which is optimized for nuclei segmentation.
In [9]:
model = models.CellposeModel(gpu=False, model_type="nuclei")
100%|██████████| 25.3M/25.3M [00:01<00:00, 18.7MB/s]
Then, we estimate the average diameter in pixels and then use Cellpose to perform segmentation.
We start by setting diameter_um to an estimated diameter in micrometers, which is based on biological knowledge, but this value can be adjusted based on the specific dataset.
To convert this estimate into pixel units, we multiply diameter_um by the scaling factor x_scale (pixels per micrometer), calculated previously. This results in diameter_pixels, an approximation of the diameter in the pixel space of the image. By converting the diameter to pixels, we ensure that Cellpose can interpret the size parameter relative to the image’s resolution.
In [10]:
diameter_um = 5 # µm
diameter_pixels = diameter_um * x_scale
print(f"Estimated diameter: {diameter_pixels:.1f} pixels")
Estimated diameter: 11.8 pixels
The eval() function returns several outputs:
masks: a labeled mask array where each detected nuclei has a unique identifier,
flows: which provides information about nuclei boundary flows,
styles: representing style vectors for detected objects, and
Note: The segmentation code below will take a few minutes to run. You can read about the parameters flow_threshold and cellprob_threshold here
In [11]:
masks, flows, styles = model.eval(
image_channel,
diameter=diameter_pixels,
flow_threshold=0.4,
cellprob_threshold=0
)
print(f"Number of nuclei detected by Cellpose: {masks.max()}")
Number of nuclei detected by Cellpose: 2911
Inspect the result¶
show_segmentation gives you the input, the predicted flows and the resulting
masks side by side. The flows are worth a look: they show how Cellpose
decided where the boundaries are, which often explains a bad segmentation
better than the masks do.
In [12]:
import base64
import io
import plotly.graph_objects as go
from matplotlib.colors import hsv_to_rgb
from PIL import Image
from skimage.segmentation import find_boundaries
def to_data_uri(rgb, quality=50):
"""Encode an RGB uint8 array as a JPEG data URI that plotly can draw directly."""
buf = io.BytesIO()
Image.fromarray(rgb).save(buf, format="JPEG", quality=quality)
return "data:image/jpeg;base64," + base64.b64encode(buf.getvalue()).decode()
def random_colour_fill(label_img, fill_alpha=0.45, edge_alpha=0.9, seed=0):
"""RGBA layer with a random colour per label; the label edges are more opaque than the fill."""
n = label_img.max() + 1
rng = np.random.default_rng(seed)
lut = (hsv_to_rgb(np.column_stack([rng.random(n), np.full(n, 0.8), np.ones(n)])) * 255).astype(np.uint8)
lut[0] = 0 # background stays transparent
rgba = np.zeros(label_img.shape + (4,), dtype=np.uint8)
rgba[..., :3] = lut[label_img]
alpha = rgba[..., 3]
alpha[label_img > 0] = int(fill_alpha * 255)
alpha[find_boundaries(label_img, mode="inner")] = int(edge_alpha * 255)
return rgba
def outlines(label_img, rgb, alpha=0.9):
"""RGBA layer with the label boundaries in one colour, transparent elsewhere."""
rgba = np.zeros(label_img.shape + (4,), dtype=np.uint8)
rgba[..., :3] = rgb
rgba[..., 3] = find_boundaries(label_img, mode="thick") * int(alpha * 255)
return rgba
def blend(base_rgb, *layers):
"""Alpha-blend RGBA layers onto an RGB image (all uint8), one colour channel at a time."""
out = base_rgb.copy()
for layer in layers:
alpha = layer[..., 3].astype(np.float32) / 255
for c in range(3):
out[..., c] = out[..., c] * (1 - alpha) + layer[..., c] * alpha
return out
gray = (np.clip((image_channel - vmin) / (vmax - vmin), 0, 1) * 255).astype(np.uint8)
image_rgb = np.repeat(gray[..., None], 3, axis=-1)
flows_rgb = flows[0] # RGB rendering of the predicted flows, the same one show_segmentation draws
cellpose_layer = random_colour_fill(masks)
views = { # delete the entries you do not need to keep the output small
"Image": image_rgb,
"Image + Cellpose": blend(image_rgb, cellpose_layer),
"Flows": flows_rgb,
}
fig = go.Figure()
for i, (name, rgb) in enumerate(views.items()):
fig.add_trace(go.Image(source=to_data_uri(rgb), name=name, visible=(i == 0), hoverinfo="skip"))
fig.update_layout(
height=750,
margin=dict(l=0, r=0, t=60, b=0),
dragmode="pan",
updatemenus=[dict(
type="buttons", direction="left", x=0, xanchor="left", y=1.0, yanchor="bottom", active=0,
buttons=[
dict(label=name, method="restyle", args=[{"visible": [j == i for j in range(len(views))]}])
for i, name in enumerate(views)
],
)],
)
fig.update_xaxes(title="x (pixels)", showgrid=False)
fig.update_yaxes(title="y (pixels)", showgrid=False, scaleanchor="x")
In [13]:
### Expanding the nuclei to cells
EXPANSION_UM = 15 # µm (In Xenium Onboard Analysis v1.0 - 1.9, the default expansion distance was 15 µm. It is reduced to 5 µm in v2.0 and later.)
expansion_px = round(EXPANSION_UM * x_scale)
cells = expand_labels(masks, distance=expansion_px)
print(f"Expanded {masks.max()} nuclei by {EXPANSION_UM} µm ({expansion_px} px)")
Expanded 2911 nuclei by 15 µm (35 px)
4. Assigning transcripts to cells¶
Now we connect the imaging to the transcriptomics. Each transcript has a
position in microns; each pixel of masks carries a cell ID. Converting the
transcript coordinates into pixel indices lets us simply look up which cell —
if any — a transcript falls inside.
In [14]:
#Save the coordinate transformation info
transformations = get_transformation(sdata.points["transcripts"], get_all=True)
#Load into memory and reset the index
pts_df = sdata.points["transcripts"].compute().reset_index(drop=True)
#Repackage for spatialdata
sdata.points["transcripts"] = PointsModel.parse(pts_df, transformations=transformations)
In [15]:
transcripts = sdata.transform_element_to_coordinate_system(
"transcripts", "global"
).compute()
print(transcripts["x"].min(), transcripts["x"].max())
5899.00448069853 12907.996323529413
In [16]:
transcripts
Out[16]:
| feature_name | cell_id | transcript_id | qv | fov_name | nucleus_distance | overlaps_nucleus | x | y | z | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Cntnap4 | lifokoje-1 | 281539401463113 | 33.907513 | B3 | 2.367967 | 0 | 5901.669922 | 8299.000460 | 21.963037 |
| 1 | Apoe | liempdpj-1 | 281539401463118 | 23.799747 | B3 | 0.000000 | 1 | 5901.749770 | 8081.004136 | 20.539789 |
| 2 | Nrep | bajngdok-1 | 281539401463122 | 40.000000 | B3 | 1.743729 | 0 | 5901.528608 | 7965.936351 | 20.591566 |
| 3 | Neto2 | kamlgfpp-1 | 281539401463125 | 40.000000 | B3 | 1.502710 | 0 | 5901.668773 | 7790.989200 | 21.805622 |
| 4 | Apoe | kgaoibkm-1 | 281539401463127 | 40.000000 | B3 | 4.311846 | 0 | 5899.079733 | 7269.482422 | 20.703056 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1261743 | Unc13c | UNASSIGNED | 281595235878634 | 40.000000 | C5 | 18.305643 | 0 | 12833.463925 | 8505.659467 | 31.441938 |
| 1261744 | Unc13c | obkpjmge-1 | 281595235878637 | 20.607004 | C5 | 0.000000 | 1 | 12838.781020 | 8613.837316 | 31.096588 |
| 1261745 | Unc13c | UNASSIGNED | 281595235878641 | 40.000000 | C5 | 19.928606 | 0 | 12840.970818 | 8496.756089 | 31.529413 |
| 1261746 | Ctsb | obmifhpe-1 | 281595235878649 | 7.263649 | C5 | 0.576466 | 0 | 12853.004366 | 8958.276654 | 31.497055 |
| 1261747 | Ctsb | obmifhpe-1 | 281595235878657 | 36.340111 | C5 | 0.000000 | 1 | 12870.541131 | 8944.841452 | 31.550989 |
1261748 rows × 10 columns
Why we can't use these coordinates directly¶
masks (from Cellpose) is a plain pixel grid, indexed from (0, 0) at the
corner of the crop. So before we can look up masks[row, col] for a
transcript, we have to convert its real-world micrometer position into a pixel
position relative to the crop
In [17]:
# Where does the crop sit? (in micrometers, on the original slide)
extent = get_extent(sdata.images["morphology_mip"], coordinate_system="global")
x0, x1 = extent["x"]
y0, y1 = extent["y"]
print(f"crop spans x: {x0:.1f} - {x1:.1f} um")
print(f"crop spans y: {y0:.1f} - {y1:.1f} um")
height, width = cells.shape # pixel size of the image we segmented
# shift to crop-relative microns, then scale to pixels
x_indices = ((transcripts["x"].values - x0) / (x1 - x0) * width).astype(int)
y_indices = ((transcripts["y"].values - y0) / (y1 - y0) * height).astype(int)
crop spans x: 5899.0 - 12908.0 um crop spans y: 3436.0 - 8995.0 um
In [18]:
# Keep only transcripts that land inside the image after conversion
inside = (
(x_indices >= 0) & (x_indices < width) &
(y_indices >= 0) & (y_indices < height)
)
print(f"{inside.sum():,} of {len(transcripts):,} transcripts fall inside the image")
transcripts = transcripts[inside].copy()
x_indices = x_indices[inside]
y_indices = y_indices[inside]
1,261,748 of 1,261,748 transcripts fall inside the image
If that count came back near zero, the shift-and-scale above did not run, or
masks is not the one you just produced with Cellpose. If it came back near
the full transcript count, you're good to continue.
In [19]:
transcripts["cellpose_id"] = cells[y_indices, x_indices]
assigned = transcripts[transcripts["cellpose_id"] > 0]
print(f"Assigned to a object: {len(assigned):,} ({100 * len(assigned) / len(transcripts):.1f}%)")
Assigned to a object: 1,012,435 (80.2%)
In [20]:
# Interactive version of the plot above, built like the one in section 3 (it reuses image_rgb,
# cellpose_layer, blend and to_data_uri from there). The assigned transcripts are a WebGL point
# cloud drawn on top of the image: hover over a dot to see its pixel coordinates.
is_assigned = transcripts["cellpose_id"].values > 0
backgrounds = {
"Image": image_rgb,
"Image + Cellpose": blend(image_rgb, cellpose_layer),
}
fig = go.Figure()
for i, (name, rgb) in enumerate(backgrounds.items()):
fig.add_trace(go.Image(source=to_data_uri(rgb, quality=100), name=name, visible=(i == 0), hoverinfo="skip"))
fig.add_trace(go.Scattergl(
x=x_indices[is_assigned].astype(np.uint16), # 16-bit integers keep a million points small in the output
y=y_indices[is_assigned].astype(np.uint16),
mode="markers",
marker=dict(size=4, color="red", opacity=1),
name="assigned transcripts",
visible=False,
hovertemplate="x: %{x} px<br>y: %{y} px<extra></extra>",
))
# Which of the three traces (image, image + cellpose, transcripts) each button shows
views = {
"Image": (True, False, False),
"Image + Mask + Transcript": (False, True, True),
}
fig.update_layout(
height=750,
margin=dict(l=0, r=0, t=60, b=0),
dragmode="pan",
showlegend=False,
updatemenus=[dict(
type="buttons", direction="left", x=0, xanchor="left", y=1.0, yanchor="bottom", active=0,
buttons=[dict(label=name, method="restyle", args=[{"visible": list(shown)}]) for name, shown in views.items()],
)],
)
height, width = cells.shape
fig.update_xaxes(title="x (pixels)", showgrid=False, range=[-0.5, width - 0.5]) # pin the axes to the image
fig.update_yaxes(title="y (pixels)", showgrid=False, range=[height - 0.5, -0.5], scaleanchor="x")
Build a cell-by-gene matrix¶
Counting transcripts per cell per gene gives us the expression matrix that all downstream analysis is built on. This is, in essence, what the 10x pipeline does — we have just done it with our own segmentation.
In [21]:
expression_matrix = (
assigned
.groupby(["cellpose_id", "feature_name"])
.size()
.unstack(fill_value=0)
)
print(f"Expression matrix: {expression_matrix.shape[0]} cells x {expression_matrix.shape[1]} genes")
expression_matrix.iloc[:5, :5]
/tmp/ipykernel_749/2175698165.py:3: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning. .groupby(["cellpose_id", "feature_name"])
Expression matrix: 2911 cells x 541 genes
Out[21]:
| feature_name | Sox10 | Ctsh | Wfs1 | Osmr | Fign |
|---|---|---|---|---|---|
| cellpose_id | |||||
| 1 | 0 | 0 | 0 | 0 | 0 |
| 2 | 0 | 0 | 0 | 0 | 0 |
| 3 | 0 | 0 | 0 | 0 | 0 |
| 4 | 0 | 0 | 2 | 1 | 0 |
| 5 | 0 | 0 | 0 | 1 | 0 |
Plot a gene in space¶
Each cell is drawn at its centroid, coloured by the expression of one gene.
In [22]:
# Interactive version of the plot above, drawn as a label image instead of dots: every Cellpose cell is
# filled with the Viridis colour of its ztranscript count for the gene, on a black background.
from matplotlib import colormaps
from matplotlib.colors import Normalize
gene_of_interest = expression_matrix.sum().idxmax() # the most-detected gene
print(f"Showing: {gene_of_interest}")
props = pd.DataFrame(regionprops_table(cells, properties=["label", "centroid"]))
props = props.set_index("label")
common = expression_matrix.index.intersection(props.index)
values = expression_matrix.loc[common, gene_of_interest].values
# Look-up table from cell label to a colour level: 0 = background, 1..255 = transcript count
norm = Normalize(vmin=values.min(), vmax=values.max())
levels = np.full(cells.max() + 1, 1, dtype=np.uint8) # every cell starts at the colour for 0 transcripts
levels[0] = 0 # background stays black
levels[common.to_numpy()] = 1 + np.round(np.asarray(norm(values)) * 254).astype(np.uint8)
palette = np.zeros((256, 3), dtype=np.uint8)
palette[1:] = (colormaps["viridis"](np.linspace(0, 1, 255))[:, :3] * 255).astype(np.uint8)
expression_image = Image.fromarray(levels[cells]) # a paletted PNG: sharp cell edges, small file
expression_image.putpalette(palette.ravel().tolist())
buf = io.BytesIO()
expression_image.save(buf, format="PNG")
image_uri = "data:image/png;base64," + base64.b64encode(buf.getvalue()).decode()
fig = go.Figure()
fig.add_trace(go.Image(source=image_uri, hoverinfo="skip"))
fig.add_trace(go.Scatter( # empty trace that only carries the colorbar
x=[None], y=[None], mode="markers", hoverinfo="skip", showlegend=False,
marker=dict(colorscale="Viridis", cmin=values.min(), cmax=values.max(), showscale=True,
colorbar=dict(title="transcript count")),
))
height, width = cells.shape
fig.update_layout(title=f"Expression of {gene_of_interest}", height=750, margin=dict(l=0, r=0, t=60, b=0), dragmode="pan")
fig.update_xaxes(title="x (pixels)", showgrid=False, range=[-0.5, width - 0.5])
fig.update_yaxes(title="y (pixels)", showgrid=False, range=[height - 0.5, -0.5], scaleanchor="x")
Showing: Cst3
5. Comparing against the 10x segmentation¶
The dataset already contains 10x's own segmentation as a label image, on the same pixel grid and with the same pyramid levels as the morphology image. That makes the comparison straightforward.
In [23]:
ref_multiscale = sdata.labels["cell_labels"]
original_masks = ref_multiscale[SCALE_KEY]["image"].data.compute()
# Guard against any off-by-one between pyramid levels
if original_masks.shape != masks.shape:
print(f"Resizing reference {original_masks.shape} -> {masks.shape}")
original_masks = resize(
original_masks, masks.shape, order=0,
preserve_range=True, anti_aliasing=False,
).astype(original_masks.dtype)
# Renumber the 10x labels to 1..N (their IDs come from the whole, uncropped slide)
original_masks = relabel_sequential(original_masks)[0]
def n_objects(label_img):
return len(np.unique(label_img)) - 1 # minus 1 for the background (0)
print(f"Cellpose cells : {n_objects(cells)}")
print(f"10x cells : {n_objects(original_masks)}")
Cellpose cells : 2911 10x cells : 5568
In [24]:
# Extent of the morphology_mip
extent_mip = get_extent(sdata.images["morphology_mip"], coordinate_system="global")
print("morphology_mip extent:")
print(f" x: {extent_mip['x']}")
print(f" y: {extent_mip['y']}")
# Extent of the cell_labels
extent_labels = get_extent(sdata.labels["cell_labels"], coordinate_system="global")
print("\nnucleus_labels extent:")
print(f" x: {extent_labels['x']}")
print(f" y: {extent_labels['y']}")
morphology_mip extent: x: (np.float64(5899.0), np.float64(12908.0)) y: (np.float64(3436.0), np.float64(8995.0)) nucleus_labels extent: x: (np.float64(5899.0), np.float64(12908.0)) y: (np.float64(3436.0), np.float64(8995.0))
Overlaying the segmentations from cellpose and 10X segmentation
In [25]:
# The image with the Cellpose masks filled red and the
# 10x segmentation drawn as blue outlines. Built with the helpers from section 3.
def colour_labels(label_img, alpha=0.4, seed=0):
"""RGBA layer where each label gets its own random colour (background stays transparent)."""
rng = np.random.default_rng(seed)
lut = rng.integers(40, 256, size=(int(label_img.max()) + 1, 3), dtype=np.uint8)
lut[0] = 0 # background
rgba = np.zeros(label_img.shape + (4,), dtype=np.uint8)
rgba[..., :3] = lut[label_img] # vectorised: one colour per label ID
rgba[..., 3] = (label_img > 0) * int(alpha * 255)
return rgba
overlay_rgb = blend(image_rgb, colour_labels(cells), outlines(original_masks, (0, 0, 255)))
#overlay_rgb = blend(image_rgb, colour_fill(cells, (255, 0, 0)), outlines(original_masks, (0, 0, 255)))
fig = go.Figure(go.Image(source=to_data_uri(overlay_rgb, quality=100), hoverinfo="skip"))
fig.update_layout(title="Red: Cellpose masks Blue: 10x outlines", height=750,
margin=dict(l=0, r=0, t=60, b=0), dragmode="pan")
fig.update_xaxes(title="x (pixels)", showgrid=False)
fig.update_yaxes(title="y (pixels)", showgrid=False, scaleanchor="x")
Quantifying the agreement¶
Two standard overlap measures, both computed on binary foreground/background masks:
- Jaccard index (intersection over union) — the fraction of the combined foreground that both methods agree on.
- Dice coefficient — similar, but weights the intersection more heavily. Always the larger of the two.
In [26]:
cellpose_binary = (cells > 0).astype(np.uint8).ravel()
original_binary = (original_masks > 0).astype(np.uint8).ravel()
jaccard = jaccard_score(original_binary, cellpose_binary)
dice = 2 * jaccard / (1 + jaccard) # Dice and Jaccard are related exactly
print(f"Jaccard index : {jaccard:.4f}")
print(f"Dice coefficient : {dice:.4f}")
Jaccard index : 0.8702 Dice coefficient : 0.9306
In [27]:
cellpose_areas = pd.DataFrame(regionprops_table(cells, properties=["area"]))
cellpose_areas["method"] = "Cellpose"
original_areas = pd.DataFrame(regionprops_table(original_masks, properties=["area"]))
original_areas["method"] = "10x"
areas_df = pd.concat([cellpose_areas, original_areas], ignore_index=True)
plt.figure(figsize=(10, 5))
sns.kdeplot(data=areas_df, x="area", hue="method", common_norm=False, fill=True, alpha=0.3)
plt.xlabel("Object area (pixels)")
plt.title("Size distributions")
plt.show()
areas_df.groupby("method")["area"].describe()
Out[27]:
| count | mean | std | min | 25% | 50% | 75% | max | |
|---|---|---|---|---|---|---|---|---|
| method | ||||||||
| 10x | 5568.0 | 1445.999461 | 1146.498299 | 1.0 | 501.75 | 1155.5 | 2045.0 | 6507.0 |
| Cellpose | 2911.0 | 2414.027138 | 1224.877083 | 194.0 | 1481.50 | 2206.0 | 3201.5 | 6762.0 |
6. Saving your results¶
Uncomment to write the masks and metrics out. Keeping the segmentation means you can come back to it without re-running Cellpose.
In [28]:
# import tifffile
# from pathlib import Path
#
# out = Path("results")
# out.mkdir(exist_ok=True)
#
# tifffile.imwrite(out / "cellpose_masks.tif", masks.astype(np.uint16))
# expression_matrix.to_csv(out / "cellpose_expression_matrix.csv")
# pd.DataFrame({
# "metric": ["jaccard", "dice", "n_cellpose", "n_10x"],
# "value": [jaccard, dice, masks.max(), original_masks.max()],
# }).to_csv(out / "segmentation_metrics.csv", index=False)
A better way to save: back into the SpatialData object¶
Writing the masks as a labels element keeps them aligned with everything else in the dataset.
In [29]:
# ---------------------------------------------------------------------------
# 1. Build the transformation: our pixels -> full-resolution pixels -> global
# ---------------------------------------------------------------------------
# How much bigger is scale0 than the level we segmented at? (per axis, ~2 for scale1)
#full_h, full_w = img_multiscale["scale0"]["image"].shape[-2:]
#sy, sx = full_h / masks.shape[0], full_w / masks.shape[1]
#print(f"Scale factor to full resolution: y={sy:.4f}, x={sx:.4f}")
# full-resolution pixels -> global (same as the 10x labels; includes the crop offset)
#to_global = get_transformation(sdata.labels["nucleus_labels"], to_coordinate_system="global")
# our pixels -> full-resolution pixels -> global
#transform = Sequence([Scale([sy, sx], axes=("y", "x")), to_global])
# ---------------------------------------------------------------------------
# 2. Add the Cellpose nuclei and expanded cells as label layers (in memory)
# ---------------------------------------------------------------------------
#for name, lab in [("cellpose_nuclei", masks), ("cellpose_cells", cells)]:
#element = Labels2DModel.parse(lab, dims=("y", "x"))
#set_transformation(element, transform, to_coordinate_system="global")
#sdata.labels[name] = element
# ---------------------------------------------------------------------------
# 3. Check the alignment BEFORE writing: both extents should match (within ~1 px)
# ---------------------------------------------------------------------------
#print("Image extent :", get_extent(sdata.images["morphology_mip"], coordinate_system="global"))
#print("Labels extent :", get_extent(sdata.labels["cellpose_nuclei"], coordinate_system="global"))
# ---------------------------------------------------------------------------
# 4. Write everything to a NEW zarr (the original stays untouched)
# ---------------------------------------------------------------------------
# Newer pandas stores text as a new string type that anndata only writes when allowed.
#ad.settings.allow_write_nullable_strings = True
#OUT_PATH = Path("../data/Practical1_cellpose.zarr")
#assert OUT_PATH.resolve() != ZARR_PATH.resolve(),
#sdata.write(OUT_PATH, overwrite=True) # full copy: original data + Cellpose labels
# ---------------------------------------------------------------------------
# 5. Reload from disk to confirm the new layers are there
# ---------------------------------------------------------------------------
#sdata_check = sd.read_zarr(OUT_PATH)
#print(sdata_check)
To try:
- Try altering the
diameter,flow_thresholdandcellprob_thresholdto see how it impacts the segmentation - Try segmenting the images by choosing a different model in cellpose
In [ ]: