Postprocesamiento para segmentación celular basada en membranas

En este notebook utilizamos un algoritmo de watershed con semillas de napari-segment-blobs-and-things-with-membranes para segmentar células y posteriormente filtrar las células para eliminar objetos segmentados erróneamente.

import matplotlib.pyplot as plt
import numpy as np
from skimage.io import imread

import pyclesperanto_prototype as cle
import napari_segment_blobs_and_things_with_membranes as nsbatwm

Los siguientes datos de imagen de ejemplo son un timelapse de células marcadas con un marcador de membrana. El conjunto de datos es cortesía de Sascha M. Kuhn, Nadler Lab, MPI-CBG Dresden.

image_timelapse = imread("../../data/membrane_2d_timelapse.tif")
image_timelapse.shape

(5, 256, 256)

Comenzamos extrayendo un ejemplo de un solo canal y un solo punto de tiempo del medio de la pila del timelapse.

membranes_single_slice = image_timelapse[2]
membranes_single_slice.shape

(256, 256)

cle.imshow(membranes_single_slice)

El plugin napari-segment-blobs-and-things-with-membranes es programable. Todos los comandos que se pueden llamar desde el menú, también se pueden llamar desde Python. Por ejemplo, el watershed con semillas de mínimos locales se puede llamar así:

cell_labels_pre = nsbatwm.local_minima_seeded_watershed(membranes_single_slice, spot_sigma=7)

cle.imshow(cell_labels_pre, labels=True)

Obviamente, hay demasiados objetos segmentados y etiquetados después de este paso. Por lo tanto, deberíamos explorar qué propiedades de esas células nos permiten diferenciar células reales y fondo segmentado. Por ejemplo, como las células tienen membrana también en la parte frontal y trasera, percibimos esta señal como una intensidad más alta en las células reales. Para visualizar esto cuantitativamente, dibujamos un mapa de intensidad media de las células. En este mapa, cada píxel que pertenece a una célula dada obtiene la intensidad de señal promedio media de toda la célula.

mean_intensity_map = cle.mean_intensity_map(membranes_single_slice, cell_labels_pre)

cle.imshow(mean_intensity_map, colorbar=True, colormap='jet')

A partir de esta vista y la barra de color al lado, podemos adivinar un umbral, por ejemplo 650, que nos permite separar las células reales del fondo. Si tal método de umbralización simple no funciona y/o un solo parámetro como la intensidad media no permite esta diferenciación, puede tener sentido usar clasificadores de objetos, también conocidos como aprendizaje automático. Por ejemplo, el plugin napari-accelerate-pixel-and-object-classification y la biblioteca subyacente apoc pueden ayudar.

cell_labels = cle.exclude_labels_with_map_values_within_range(mean_intensity_map, cell_labels_pre, maximum_value_range=700)

cle.imshow(cell_labels, labels=True)

Aplicando el flujo de trabajo a todo el timelapse

Ahora tenemos una idea aproximada de cómo se ve un flujo de trabajo para segmentar las células. Reescribimos todo el flujo de trabajo en una función. La función también tiene una documentación adecuada para que el yo futuro aún sepa lo que hace la función.

def cell_segmentation_worflow(membranes_single_slice, cell_segmentation_spot_sigma=7, maximum_mean_intensity_per_cell=700):
    """Segmentación celular basada en señal de membrana
    
    Parámetros
    ----------
    membranes_single_slice
        Imagen a segmentar
    cell_segmentation_spot_sigma: float, opcional
        Permite configurar la segmentación. Cuanto mayor sea sigma, menos células habrá.
    maximum_mean_intensity_per_cell: float, opcional
        Umbral para la intensidad media por célula. Las células con señal por debajo serán excluidas
        
    Retorna
    -------
    imagen de etiquetas celulares
    """
    
    # etiquetar candidatos celulares
    cell_labels_pre = nsbatwm.local_minima_seeded_watershed(
                                membranes_single_slice, 
                                spot_sigma=cell_segmentation_spot_sigma)
    
    # medir intensidad en todos los candidatos
    mean_intensity_map = cle.mean_intensity_map(
                                membranes_single_slice, 
                                cell_labels_pre)
    
    # excluir candidatos con señal baja
    cell_labels = cle.exclude_labels_with_map_values_within_range(
                        mean_intensity_map, cell_labels_pre, 
                        maximum_value_range=maximum_mean_intensity_per_cell)
    
    return cell_labels

Esta función, podemos aplicarla a todo el timelapse y ver si la segmentación también funciona bien en todos los puntos de tiempo.

for t in range(len(image_timelapse)):
    # extraer un solo corte / punto de tiempo
    membranes_single_slice = image_timelapse[t]
    
    # segmentar células 
    cell_labels = cell_segmentation_worflow(membranes_single_slice)

    # imprimir encabezado
    print("t = ", t, ", número de células: ", cell_labels.max())
    
    # mostrar tres imágenes: imagen, etiquetas, superposición
    fix, axs = plt.subplots(1, 3, figsize=(10,10))
    
    cle.imshow(membranes_single_slice, plot=axs[0])
    cle.imshow(cell_labels, plot=axs[1], labels=True)
    cle.imshow(membranes_single_slice, plot=axs[2], alpha=0.5, continue_drawing=True)
    cle.imshow(cell_labels, plot=axs[2], labels=True, alpha=0.5)
    plt.show()
    

t =  0 , number of cells:  17.0

t =  1 , number of cells:  18.0

t =  2 , number of cells:  17.0

t =  3 , number of cells:  18.0

t =  4 , number of cells:  17.0

Ejercicio

En lugar de filtrar células con baja intensidad de señal, filtra células que son grandes y/o tocan el borde de la imagen. Aplica el nuevo flujo de trabajo a todo el timelapse.