Feature extraction and measurements

Duration ~40 min  ·  Session Day 2, Python notebooks

Everything so far has been about getting a segmentation. None of it was the point. Nobody publishes a binary image; what gets published is what the binary image made it possible to measure.

This notebook turns a label image into a table of numbers, which is where image analysis stops being about images and becomes ordinary data analysis.

This notebook covers

• turning a label image into a pandas table with regionprops_table

• morphology features versus intensity features, and when each is useful

• measuring one channel through another channel’s binary image

• comparing a measured population between conditions

Data: data/bbbc020/, the nuclei and cells segmented so far. These images depict bone-marrow derived macrophages from C57BL/6 mice. Based on this experiment, the role of toll-like receptor 4 (TLR 4), Myd88 and MAPK for macrophage spreading was investigated. DAPI labels the nuclei and CD11b/APC the cell surface.

1. From labels to a table of measurements

rp = io.imread(DATA / 'misc' / 'sphx_glr_plot_regionprops_001.png')
imshow(rp)
# (image source: https://scikit-image.org/docs/stable/auto_examples/segmentation/plot_regionprops.html)

import re
import tifffile
import numpy as np
import pandas as pd
from course import DATA, show
import matplotlib.pyplot as plt
from skimage.color import rgb2gray
from skimage import io
from iaf.plot import imshow, show_labels
from skimage.measure import regionprops_table
import os; os.environ["DISPLAY"] = ":1"

cells = tifffile.imread(DATA / "bbbc020" / "images" / "jw-24h 2_c1.tif")
nuclei = tifffile.imread(DATA / "bbbc020" / "images" / "jw-24h 2_c5.tif")
labels = tifffile.imread(DATA / "bbbc020" / "gt" / "jw-24h 2_c5 [1].tif")

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 5))
imshow(cells, ax=ax1, title="Cell intensities (raw data, ch1=CD11b/APC)")
imshow(nuclei, ax=ax2, title="Nuceli intensities (raw data, ch5=DAPI)")
show_labels(labels, ax=ax3, title="Nuceli GT labels")
fig.subplots_adjust(wspace=0.3) 
print("label image:", labels.shape, labels.dtype, "|", labels.max(), "objects")

label image: (384, 512) uint16 | 33 objects

We are using the annotated labels (ground truth) here rather than our own segmentation, so that the measurements are not confounded by segmentation error. In practice this would be the output of a segmentation pipeline or manual / expert ground truth.

table = regionprops_table(
    labels,
    properties=("label", "area", "perimeter", "eccentricity", "solidity", "extent"),
)
features = pd.DataFrame(table)
print(features.shape)
features.head() # we only show the first few elements of the table

(33, 6)

	label	area	perimeter	eccentricity	solidity	extent
0	1	367.0	74.284271	0.741365	0.933842	0.764583
1	2	687.0	98.325902	0.625378	0.963534	0.748366
2	3	334.0	73.012193	0.817711	0.922652	0.607273
3	4	647.0	100.911688	0.831878	0.932277	0.660204
4	5	373.0	72.183766	0.664171	0.953964	0.706439

One row per object, one column per measurement. That is the whole idea, and from here on nothing is image-specific: this is a table, and tables are familiar ground.

regionprops_table takes the property names as strings. There are dozens available; the full list is in the regionprops documentation.

2. Morphology features

These describe shape only. They can be computed from the label image alone: no intensities involved.

feature	meaning	intuition
area	pixels in the object	size
perimeter	pixels around its edge	:
eccentricity	0 = circle, →1 = line	how elongated
solidity	area ÷ convex-hull area	how “dented”; 1 = convex
extent	area ÷ bounding-box area	how well it fills its box

The last three are dimensionless ratios, which makes them far more portable than area: they do not change with magnification.

Working with the measurement results table lets us summarize the data and calc the stats of all nuceli per feature.

features[["area", "perimeter", "eccentricity", "solidity", "extent"]].describe().round(2)

	area	perimeter	eccentricity	solidity	extent
count	36.00	36.00	36.00	36.00	36.00
mean	2740.83	203.27	0.75	0.96	0.75
std	1282.39	42.80	0.13	0.02	0.06
min	320.00	104.35	0.41	0.92	0.63
25%	2078.75	179.99	0.67	0.96	0.70
50%	2430.00	197.32	0.77	0.97	0.75
75%	2979.00	221.53	0.82	0.98	0.79
max	5610.00	288.31	0.99	0.99	0.85

fig, axes = plt.subplots(1, 4, figsize=(17, 3.2))
for ax, column in zip(axes, ["area", "eccentricity", "solidity", "extent"]):
    ax.hist(features[column], bins=15, color="0.4")
    ax.set_title(column)
    ax.set_ylabel("nuclei" if column == "area" else "")
plt.tight_layout()
plt.show()

It’s good practice to plot a histogram and look at the distribution of the data. For example, area ranges over a factor of five, while solidity sits between 0.93 and 0.98 for every single nucleus, because nuclei are convex blobs and solidity measures exactly that.

A feature with no variance carries no information. Therefore, solidity is useless for telling these nuclei apart.

3. Intensity features

These need a second input: the image whose intensities are to be summarised inside each object. Pass it as intensity_image.

This is the same idea as Redirect to… in Fiji’s Set Measurements: the binary image comes from one channel, the numbers from another.

both = pd.DataFrame(regionprops_table(
    labels,
    intensity_image = rgb2gray(cells),
    properties=("label", "area", "intensity_mean", "intensity_max", "intensity_min"),
))
both.head()

	label	area	intensity_mean	intensity_max	intensity_min
0	1	2683.0	0.351108	0.542185	0.265634
1	2	5070.0	0.199823	0.392993	0.098248
2	3	2519.0	0.144074	0.261995	0.076415
3	4	4832.0	0.146022	0.243801	0.076415
4	5	2806.0	0.234752	0.593128	0.149192

Every row is one nucleus, and intensity_mean is how much CD11b/APC signal of the cell surface belongs to inside that nucleus’s segmentation mask.

That is the shape of a very large fraction of all quantitative microscopy: use a reliable channel to find the objects, then measure the channel of interest inside them.

# The same objects, measured in each channel.
in_nuclei = pd.DataFrame(regionprops_table(
    labels, intensity_image=rgb2gray(nuclei), properties=("label", "intensity_mean")
)).rename(columns={"intensity_mean": "dapi"})

in_cells = both[["label", "intensity_mean"]].rename(columns={"intensity_mean": "cd11b_apc"})

merged = in_nuclei.merge(in_cells, on="label")
print(merged[["dapi", "cd11b_apc"]].mean().round(2).to_string())

plt.figure(figsize=(5, 4.5))
plt.scatter(merged.dapi, merged.cd11b_apc, s=25)
plt.xlabel("mean DAPI in nucleus mask")
plt.ylabel("mean CD11b/APC in nucleus mask")
plt.title("same objects, two channels")
plt.show()

dapi         0.03
cd11b_apc    0.26

Same nuclei, two different numbers each, and no particular relationship between them, which is what one would hope, since DAPI brightness and CD11b/APC expression have no reason to track each other. Remember DAPI stains the nucleus and CD11b/APC the cell surface.

A strong diagonal in this plot would need explaining. It is usually a sign of bleed-through between channels, or of an autofluorescent sample, and it would invalidate any claim that the two measurements are independent.

4. Measuring across a whole experiment

The value of doing this in Python code is that we can automate the signal quantification and run over every image in the experiment of a time course.

def measure_field(field):
    """Every annotated nucleus in one field, measured in both channels."""
    labels = tifffile.imread(DATA / "bbbc020" / "gt" / f"{field}_nuclei_labels.tif")
    nuclei = tifffile.imread(DATA / "bbbc020" / "images" / f"{field}_nuclei.tif")

    table = pd.DataFrame(regionprops_table(
        labels,
        intensity_image=nuclei,
        properties=("label", "area", "eccentricity", "extent", "intensity_mean"),
    ))
    table["field"] = field
    # "15min_1" -> "15min", "Kontrolle1" -> "Kontrolle": strip the trailing
    # replicate number, with or without an underscore before it.
    table["condition"] = re.sub(r"_?\d+$", "", field)
    return table


fields = ["15min_1", "15min_2", "15min_3", "1h_1", "1h_2",
          "2h_1", "2h_2", "24h_1", "24h_2",
          "Kontrolle1", "Kontrolle2", "Kontrolle3"]

population = pd.concat([measure_field(f) for f in fields], ignore_index=True)
print(f"{len(population)} nuclei from {population.field.nunique()} fields")
population.groupby("condition").size().to_frame("nuclei")

309 nuclei from 12 fields

	nuclei
condition
15min	72
1h	61
24h	42
2h	78
Kontrolle	56

Note the two extra columns. field keeps track of which image each object came from, and condition groups the fields: Kontrolle is German for control, and the rest are timepoints after stimulation.

Always keep the identifying columns. Pooling measurements from several images loses the ability to trace a result back to its source, unless the source was carried along.

5. Comparing conditions

summary = population.groupby("condition").agg(
    nuclei=("label", "size"),
    median_area=("area", "median"),
    median_dapi=("intensity_mean", "median"),
).round(1)
summary

	nuclei	median_area	median_dapi
condition
15min	72	629.5	97.3
1h	61	595.0	110.5
24h	42	574.5	143.2
2h	78	540.5	112.8
Kontrolle	56	487.0	122.6

order = ["Kontrolle", "15min", "1h", "2h", "24h"]

fig, axes = plt.subplots(1, 2, figsize=(13, 4))
for ax, column, ylabel in [
    (axes[0], "area", "nucleus area (px)"),
    (axes[1], "intensity_mean", "mean DAPI in nucleus"),
]:
    data = [population.loc[population.condition == c, column] for c in order]
    ax.boxplot(data, tick_labels=order)
    ax.set_ylabel(ylabel)
    ax.set_xlabel("condition")
plt.tight_layout()
plt.show()

A box plot shows the whole distribution rather than just its middle, which matters here: the boxes overlap heavily. Whatever difference exists between conditions is small compared with the spread within each condition. However, there is a trend that the average DAPI signal in the nuclei increases with time

Warning

Conclusions from this need care. Each condition here is one field from one experiment, and the nuclei within a field are not independent observations: they share a dish, a staining, a focus setting. Treating 300 nuclei as 300 independent samples would badly overstate the confidence.

The unit of replication is the experiment, not the cell. This is the single most common statistical error in quantitative microscopy.

6. Which features actually distinguish anything?

# A quick way to see which features differ between two conditions: compare the
# separation between group medians with the spread inside the groups.
start = population[population.condition == "15min"]
late = population[population.condition == "24h"]

print(f"{'feature':16s} {'15min':>9s} {'24h':>9s} {'difference':>11s}")
print("-" * 48)
for column in ["area", "eccentricity", "extent", "intensity_mean"]:
    a, b = start[column], late[column]
    pooled_sd = np.sqrt((a.var() + b.var()) / 2)
    effect = (b.mean() - a.mean()) / pooled_sd if pooled_sd else 0
    print(f"{column:16s} {a.median():9.2f} {b.median():9.2f} {effect:11.2f}")

feature              15min       24h  difference
------------------------------------------------
area                629.50    574.50       -0.15
eccentricity          0.66      0.76        0.65
extent                0.80      0.79       -0.43
intensity_mean       97.26    143.24        1.07

That last column is the difference in means expressed in standard deviations: an effect size. It is far more informative than a p-value on its own, because it says how big the difference is rather than only whether it is detectable. Crucially, standardizing the difference allows you to directly compare completely different measurements (like pixel area versus intensity) on the exact same scale. A rule of thumb: 0.2 is small, 0.5 moderate, 0.8 large.

Notice that the features disagree about which condition differs and by how much. That is normal, and it is why picking a feature after seeing the result is such a dangerous habit: with enough features, something always looks significant. Good experimental design defines the question before looking at the results, preventing researchers from cherry-picking data that matches their assumptions.

Recap

	computed from	examples	portable across magnifications?
morphology	the label image alone	area, eccentricity, solidity, extent	ratios yes, area no
intensity	label image + an intensity image	intensity_mean, intensity_max	yes

Take home:

• keep field and condition columns so pooled measurements stay traceable;

• a feature with no variance carries no information;

• the unit of replication is the experiment, not the cell;

• report an effect size, not just a difference.