Counting nuclei according to expression in multiple channels#

A common bio-image analysis task is counting cells according to their signal expression in multiple channels. In this example we take a two-channel image of nuclei which express Cy3 and eGFP. Visually, we can easily see that some nuclei expressing Cy3 also express eGFP, others don’t. This notebook demonstrates how to count these groups of nuclei.

import pyclesperanto_prototype as cle
import numpy as np
from skimage.io import imread, imshow
import matplotlib.pyplot as plt

cle.get_device()
<Apple M1 Max on Platform: Apple (2 refs)>

We’re using a dataset published by Heriche et al. licensed CC BY 4.0 available in the Image Data Resource.

# load file
raw_image = imread('plate1_1_013 [Well 5, Field 1 (Spot 5)].png')

# visualize
imshow(raw_image)
<matplotlib.image.AxesImage at 0x1796beb20>
../_images/counting_nuclei_multichannel_3_1.png

First, we need to split channels (read more). After that, we can actually see that not all cells marked with Cy3 (channel 0) are also marked with eGFP (channel 1):

# extract channels
channel_0 = raw_image[...,0]
channel_1 = raw_image[...,1]

# visualize
fig, axs = plt.subplots(1, 2, figsize=(15,15))
axs[0].imshow(channel_0, cmap='gray')
axs[1].imshow(channel_1, cmap='gray')
<matplotlib.image.AxesImage at 0x17b6d1160>
../_images/counting_nuclei_multichannel_5_1.png

Segmenting nuclei#

As the staining marks the whole nucleus in both cases, it is reasonable to segmentn nuclei in both images and then process the segmented images further. We use Voronoi-Otsu-Labeling for the segmentation because it is a quick and straightforward approach.

# segmentation
nuclei_cy3 = cle.voronoi_otsu_labeling(channel_0, spot_sigma=20)
nuclei_egfp = cle.voronoi_otsu_labeling(channel_1, spot_sigma=20)

# visualize
fig, axs = plt.subplots(1, 2, figsize=(15,15))
cle.imshow(nuclei_cy3, plot=axs[0], labels=True)
cle.imshow(nuclei_egfp, plot=axs[1], labels=True)
../_images/counting_nuclei_multichannel_7_0.png

The above shown label images have inside nuclei pixel intensity values that correspond to the number of the nucleus. In nucleus 1, all pixels have intensity 1. In nucleus 2, all pixels have intensity 2 and so on. Hence, from these label images, we can already determine the number of nuclei in both channels, by measuring the maximum intensity in the label images:

# determine maximum in both label images
number_of_nuclei_cy3 = nuclei_cy3.max()
number_of_nuclei_egfp = nuclei_egfp.max()

# print out result
print("Nuclei Cy3 positive:", number_of_nuclei_cy3)
print("Nuclei eGFP positive:", number_of_nuclei_egfp)
Nuclei Cy3 positive: 31.0
Nuclei eGFP positive: 23.0

Technically, we haven’t checked yet if all eGFP positive nuclei are also Cy3 positive. We can do this by determining how many eGFP positive nuclei are close by each individual Cy3 positive nucleus. Therefore, we need to set a maximum distance threshold allowing us to specify how far away centroids of nuclei are allowed to be.

maximum_distance = 15 # pixels

# draw a parametric map of cell counts
count_map = cle.proximal_other_labels_count_map(nuclei_cy3, nuclei_egfp)
cle.imshow(count_map, colorbar=True)
../_images/counting_nuclei_multichannel_11_0.png

The count_map is a parametric image. We can identify all the nuclei where the count value >= 1. These are all the Cy3-positive nuclei which have at least one eGFP-positive nucleus with a centroid distance <= 15 pixels.

double_positive_nuclei = cle.exclude_labels_with_map_values_out_of_range(
    count_map, 
    nuclei_cy3, 
    minimum_value_range=1)

cle.imshow(double_positive_nuclei, labels=True)
../_images/counting_nuclei_multichannel_13_0.png

And we can also count those similar to shown above:

number_of_double_positives = double_positive_nuclei.max()
print("Number of Cy3 positives that also express eGFP", number_of_double_positives)
Number of Cy3 positives that also express eGFP 23.0

Visualization#

We can also use the outlines around cells which are double positive and visualize those on the original images of both channels.

# determine outlines
outlines = cle.detect_label_edges(double_positive_nuclei)

# add outlines to original images. As outlines have value 1, 
# we need to multiply them to make them properly visible:
channel_0_with_outlines = cle.maximum_images(channel_0, outlines * channel_0.max())

# visualize result
cle.imshow(channel_0_with_outlines)

# let's zoom in
cle.imshow(channel_0_with_outlines.get()[400:800, 1000:1700])
../_images/counting_nuclei_multichannel_17_0.png ../_images/counting_nuclei_multichannel_17_1.png

For interactive visualization, we can also use napari:

# startup a viewer
import napari
viewer = napari.Viewer()

# add raw images in color to the viewer
viewer.add_image(channel_0, colormap='magenta')
viewer.add_image(channel_1, colormap='green', blending='additive')

# add labels and configure it so that we see the contours as thick lines
labels_layer = viewer.add_labels(double_positive_nuclei)
labels_layer.contour=5

# make a screenshot of the viewer
napari.utils.nbscreenshot(viewer)
Warning: Could not find scikit-tensor which is needed for separable approximations...
If you want to compute separable approximations, please install it with
pip install scikit-tensor-py3
/Users/haase/opt/anaconda3/envs/bio_39/lib/python3.9/site-packages/napari/plugins/_plugin_manager.py:494: UserWarning: Plugin 'napari-accelerated-pixel-and-object-classification' provided a non-callable type to `napari_experimental_provide_function`: <class 'magicgui._magicgui.MagicFactory'>. Function widget ignored.
  warn(message=warn_message)

References#

Some of the functions we used might be uncommon. Hence, we can add their documentation for reference.

print(cle.voronoi_otsu_labeling.__doc__)
Labels objects directly from grey-value images.

    The two sigma parameters allow tuning the segmentation result. Under the hood,
    this filter applies two Gaussian blurs, spot detection, Otsu-thresholding [2] and Voronoi-labeling [3]. The
    thresholded binary image is flooded using the Voronoi tesselation approach starting from the found local maxima.
    
    Parameters
    ----------
    source : Image
        Input grey-value image
    label_image_destination : Image, optional
        Output image
    spot_sigma : float
        controls how close detected cells can be
    outline_sigma : float
        controls how precise segmented objects are outlined.
    
    Returns
    -------
    label_image_destination
    
    Examples
    --------
    >>> import pyclesperanto_prototype as cle
    >>> cle.voronoi_otsu_labeling(source, label_image_destination, 10, 2)
    
    References
    ----------
    .. [1] https://clij.github.io/clij2-docs/reference_voronoiOtsuLabeling
    .. [2] https://ieeexplore.ieee.org/document/4310076
    .. [3] https://en.wikipedia.org/wiki/Voronoi_diagram
    
print(cle.proximal_other_labels_count_map.__doc__)
    Count number of labels within a given radius in an other label image and returns the result as parametric map.

    Parameters
    ----------
    label_image: Image
    other_label_image: Image
    count_map: Image, optional
        parametric image where the values will be written in.
    maximum_distance: Number, optional
        maximum distance in pixels

    Returns
    -------
    count_map