测量核膜处的荧光强度

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

[1]

Boni A, Politi AZ, Strnad P, Xiang W, Hossain MJ, Ellenberg J (2015) “Live imaging and modeling of inner nuclear membrane targeting reveals its molecular requirements in mammalian cells” J Cell Biol 209(5):705–720. ISSN: 0021-9525. DOI:10.1083/jcb.201409133

[2]

Miura K (2020) “Measurements of Intensity Dynamics at the Periphery of the Nucleus” in: Miura K, Sladoje N (eds) Bioimage Data Analysis Workflows. Learning Materials in Biosciences. Springer, Cham. DOI:10.1007/978-3-030-22386-1_2

[3]

Klemm A (2020) “ImageJ/Fiji Macro Language” NEUBIAS Academy Online Course: https://www.youtube.com/watch?v=o8tfkdcd3DA

[4]

Vorkel D and Haase R (2020) “GPU-accelerating ImageJ Macro image processing workflows using CLIJ” https://arxiv.org/abs/2008.11799

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()

将分割的边缘应用为掩码

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

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()

测量总强度

平均强度在标记图像中很容易作为区域属性获得。

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 使其成为 2D（其中第一维仍然对应于时间点，但第二维和最后一维现在包含所有像素值），然后沿着图像序列的第二个轴应用阈值函数

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()

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