注意

转到末尾 下载完整示例代码。 或通过 Binder 在浏览器中运行此示例。

测量核膜处的荧光强度#

此示例再现了生物图像数据分析中一个成熟的工作流程，用于测量定位在核膜处的荧光强度，在细胞图像时间序列（每个图像有两个通道和两个空间维度）中，该时间序列显示了蛋白质从细胞质区域重新定位到核膜的过程。 这一生物学应用最初由 Andrea Boni 及其合作者在 [1] 中提出； 它被 Kota Miura 的教科书 [2] 以及其他著作（[3]，[4]）使用。 换句话说，我们将此工作流程从 ImageJ 宏移植到使用 scikit-image 的 Python 中。

import matplotlib.pyplot as plt
import numpy as np
import plotly.io
import plotly.express as px
from scipy import ndimage as ndi

from skimage import filters, measure, morphology, segmentation
from skimage.data import protein_transport

我们从单个细胞/细胞核开始构建工作流程。

image_sequence = protein_transport()

print(f'shape: {image_sequence.shape}')

shape: (15, 2, 180, 183)

数据集是一个 2D 图像堆栈，包含 15 个帧（时间点）和 2 个通道。

vmin, vmax = 0, image_sequence.max()

fig = px.imshow(
    image_sequence,
    facet_col=1,
    animation_frame=0,
    zmin=vmin,
    zmax=vmax,
    binary_string=True,
    labels={'animation_frame': 'time point', 'facet_col': 'channel'},
)
plotly.io.show(fig)

首先，让我们考虑第一张图像的第一个通道（下图中的步骤 a)）。

image_t_0_channel_0 = image_sequence[0, 0, :, :]

分割细胞核边缘#

让我们对该图像应用高斯低通滤波器以使其平滑（步骤 b)）。 接下来，我们分割细胞核，使用 Otsu 方法找到背景和前景之间的阈值： 我们得到一个二值图像（步骤 c)）。 然后我们填充物体中的孔洞（步骤 c-1)）。

smooth = filters.gaussian(image_t_0_channel_0, sigma=1.5)

thresh_value = filters.threshold_otsu(smooth)
thresh = smooth > thresh_value

fill = ndi.binary_fill_holes(thresh)

遵循原始工作流程，让我们删除接触图像边界的物体（步骤 c-2)）。 在这里，我们可以看到另一个细胞核的一部分接触了右下角。

clear = segmentation.clear_border(fill)
clear.dtype

dtype('bool')

我们计算此二值图像的形态学膨胀（步骤 d)）和其形态学腐蚀（步骤 e)）。

dilate = morphology.binary_dilation(clear)

erode = morphology.binary_erosion(clear)

最后，我们从膨胀中减去腐蚀以获得细胞核边缘（步骤 f)）。 这相当于选择存在于 dilate 中但不存在于 erode 中的像素。

mask = np.logical_and(dilate, ~erode)

让我们在一系列子图中可视化这些处理步骤。

fig, ax = plt.subplots(2, 4, figsize=(12, 6), sharey=True)

ax[0, 0].imshow(image_t_0_channel_0, cmap=plt.cm.gray)
ax[0, 0].set_title('a) Raw')

ax[0, 1].imshow(smooth, cmap=plt.cm.gray)
ax[0, 1].set_title('b) Blur')

ax[0, 2].imshow(thresh, cmap=plt.cm.gray)
ax[0, 2].set_title('c) Threshold')

ax[0, 3].imshow(fill, cmap=plt.cm.gray)
ax[0, 3].set_title('c-1) Fill in')

ax[1, 0].imshow(clear, cmap=plt.cm.gray)
ax[1, 0].set_title('c-2) Keep one nucleus')

ax[1, 1].imshow(dilate, cmap=plt.cm.gray)
ax[1, 1].set_title('d) Dilate')

ax[1, 2].imshow(erode, cmap=plt.cm.gray)
ax[1, 2].set_title('e) Erode')

ax[1, 3].imshow(mask, cmap=plt.cm.gray)
ax[1, 3].set_title('f) Nucleus Rim')

for a in ax.ravel():
    a.set_axis_off()

fig.tight_layout()
a) Raw, b) Blur, c) Threshold, c-1) Fill in, c-2) Keep one nucleus, d) Dilate, e) Erode, f) Nucleus Rim

将分割的边缘用作掩码#

现在我们已经分割了第一个通道中的细胞核膜，我们将其用作掩码来测量第二个通道中的强度。

image_t_0_channel_1 = image_sequence[0, 1, :, :]
selection = np.where(mask, image_t_0_channel_1, 0)

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 6), sharey=True)

ax0.imshow(image_t_0_channel_1)
ax0.set_title('Second channel (raw)')
ax0.set_axis_off()

ax1.imshow(selection)
ax1.set_title('Selection')
ax1.set_axis_off()

fig.tight_layout()
Second channel (raw), Selection

测量总强度#

平均强度作为标记图像中的区域属性立即可用。

props = measure.regionprops_table(
    mask.astype(np.uint8),
    intensity_image=image_t_0_channel_1,
    properties=('label', 'area', 'intensity_mean'),
)

我们可能需要检查总强度值

selection.sum()

np.uint64(28350)

可以从

props['area'] * props['intensity_mean']

array([28350.])

处理整个图像序列#

无需为每个时间点迭代工作流程，我们可以直接处理多维数据集（阈值步骤除外）。 实际上，大多数 scikit-image 函数支持 nD 图像。

n_z = image_sequence.shape[0]  # number of frames

smooth_seq = filters.gaussian(image_sequence[:, 0, :, :], sigma=(0, 1.5, 1.5))
thresh_values = [filters.threshold_otsu(s) for s in smooth_seq[:]]
thresh_seq = [smooth_seq[k, ...] > val for k, val in enumerate(thresh_values)]

或者，我们可以计算 thresh_values，而不使用列表推导，方法是将 smooth_seq 重塑为二维（其中第一维仍然对应于时间点，但第二维和最后一维现在包含所有像素值），并将阈值函数应用于图像序列的第二轴

thresh_values = np.apply_along_axis(filters.threshold_otsu,
                                    axis=1,
                                    arr=smooth_seq.reshape(n_z, -1))

我们使用以下扁平结构元素进行形态学计算（np.newaxis 用于为时间添加一个大小为 1 的轴）

footprint = ndi.generate_binary_structure(2, 1)[np.newaxis, ...]
footprint

array([[[False,  True, False],
        [ True,  True,  True],
        [False,  True, False]]])

这样，每个帧都独立处理（来自连续帧的像素永远不是空间邻居）。

fill_seq = ndi.binary_fill_holes(thresh_seq, structure=footprint)

在清除接触图像边界的物体时，我们希望确保底部（第一个）和顶部（最后一个）帧不被视为边界。 在这种情况下，唯一相关的边界是最大（x，y）值处的边缘。 这可以通过运行以下代码在 3D 中看到

import plotly.graph_objects as go

sample = fill_seq
(n_Z, n_Y, n_X) = sample.shape
Z, Y, X = np.mgrid[:n_Z, :n_Y, :n_X]

fig = go.Figure(
    data=go.Volume(
        x=X.flatten(),
        y=Y.flatten(),
        z=Z.flatten(),
        value=sample.flatten(),
        opacity=0.5,
        slices_z=dict(show=True, locations=[n_z // 2])
    )
)
fig.show()

border_mask = np.ones_like(fill_seq)
border_mask[n_z // 2, -1, -1] = False
clear_seq = segmentation.clear_border(fill_seq, mask=border_mask)

dilate_seq = morphology.binary_dilation(clear_seq, footprint=footprint)
erode_seq = morphology.binary_erosion(clear_seq, footprint=footprint)
mask_sequence = np.logical_and(dilate_seq, ~erode_seq)

让我们为每个掩码（对应于每个时间点）提供一个不同的标签，从 1 到 15 运行。 我们使用 np.min_scalar_type 来确定表示时间点数目所需的最小大小整数类型

label_dtype = np.min_scalar_type(n_z)
mask_sequence = mask_sequence.astype(label_dtype)
labels = np.arange(1, n_z + 1, dtype=label_dtype)
mask_sequence *= labels[:, np.newaxis, np.newaxis]

让我们计算所有这些标记区域感兴趣的区域属性。

props = measure.regionprops_table(
    mask_sequence,
    intensity_image=image_sequence[:, 1, :, :],
    properties=('label', 'area', 'intensity_mean'),
)
np.testing.assert_array_equal(props['label'], np.arange(n_z) + 1)

fluorescence_change = [
    props['area'][i] * props['intensity_mean'][i] for i in range(n_z)
]

fluorescence_change /= fluorescence_change[0]  # normalization

fig, ax = plt.subplots()
ax.plot(fluorescence_change, 'rs')
ax.grid()
ax.set_xlabel('Time point')
ax.set_ylabel('Normalized total intensity')
ax.set_title('Change in fluorescence intensity at the nuclear envelope')
fig.tight_layout()

plt.show()
Change in fluorescence intensity at the nuclear envelope

令人欣慰的是，我们找到了预期结果： 在最初的五个时间点，细胞核膜处的总荧光强度增加了 1.3 倍，然后基本保持恒定。

脚本的总运行时间：（0 分钟 1.726 秒）

由 Sphinx-Gallery 生成的画廊

本页内容
显示源代码