.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/applications/plot_fluorescence_nuclear_envelope.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_applications_plot_fluorescence_nuclear_envelope.py: ====================================================== Measure fluorescence intensity at the nuclear envelope ====================================================== This example reproduces a well-established workflow in bioimage data analysis for measuring the fluorescence intensity localized to the nuclear envelope, in a time sequence of cell images (each with two channels and two spatial dimensions) which shows a process of protein re-localization from the cytoplasmic area to the nuclear envelope. This biological application was first presented by Andrea Boni and collaborators in [1]_; it was used in a textbook by Kota Miura [2]_ as well as in other works ([3]_, [4]_). In other words, we port this workflow from ImageJ Macro to Python with scikit-image. .. [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 .. GENERATED FROM PYTHON SOURCE LINES 31-42 .. code-block:: 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 .. GENERATED FROM PYTHON SOURCE LINES 43-44 We start with a single cell/nucleus to construct the workflow. .. GENERATED FROM PYTHON SOURCE LINES 44-49 .. code-block:: Python image_sequence = protein_transport() print(f'shape: {image_sequence.shape}') .. rst-class:: sphx-glr-script-out .. code-block:: none shape: (15, 2, 180, 183) .. GENERATED FROM PYTHON SOURCE LINES 50-51 The dataset is a 2D image stack with 15 frames (time points) and 2 channels. .. GENERATED FROM PYTHON SOURCE LINES 51-65 .. code-block:: Python 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) .. raw:: html :file: images/sphx_glr_plot_fluorescence_nuclear_envelope_001.html .. GENERATED FROM PYTHON SOURCE LINES 66-68 To begin with, let us consider the first channel of the first image (step ``a)`` in the figure below). .. GENERATED FROM PYTHON SOURCE LINES 68-71 .. code-block:: Python image_t_0_channel_0 = image_sequence[0, 0, :, :] .. GENERATED FROM PYTHON SOURCE LINES 72-79 Segment the nucleus rim ======================= Let us apply a Gaussian low-pass filter to this image in order to smooth it (step ``b)``). Next, we segment the nuclei, finding the threshold between the background and foreground with Otsu's method: We get a binary image (step ``c)``). We then fill the holes in the objects (step ``c-1)``). .. GENERATED FROM PYTHON SOURCE LINES 79-87 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 88-91 Following the original workflow, let us remove objects which touch the image border (step ``c-2)``). Here, we can see that part of another nucleus was touching the bottom right-hand corner. .. GENERATED FROM PYTHON SOURCE LINES 91-95 .. code-block:: Python clear = segmentation.clear_border(fill) clear.dtype .. rst-class:: sphx-glr-script-out .. code-block:: none dtype('bool') .. GENERATED FROM PYTHON SOURCE LINES 96-98 We compute both the morphological dilation of this binary image (step ``d)``) and its morphological erosion (step ``e)``). .. GENERATED FROM PYTHON SOURCE LINES 98-103 .. code-block:: Python dilate = morphology.binary_dilation(clear) erode = morphology.binary_erosion(clear) .. GENERATED FROM PYTHON SOURCE LINES 104-107 Finally, we subtract the eroded from the dilated to get the nucleus rim (step ``f)``). This is equivalent to selecting the pixels which are in ``dilate``, but not in ``erode``: .. GENERATED FROM PYTHON SOURCE LINES 107-110 .. code-block:: Python mask = np.logical_and(dilate, ~erode) .. GENERATED FROM PYTHON SOURCE LINES 111-112 Let us visualize these processing steps in a sequence of subplots. .. GENERATED FROM PYTHON SOURCE LINES 112-144 .. code-block:: Python 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-sg:: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_002.png :alt: a) Raw, b) Blur, c) Threshold, c-1) Fill in, c-2) Keep one nucleus, d) Dilate, e) Erode, f) Nucleus Rim :srcset: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 145-149 Apply the segmented rim as a mask ================================= Now that we have segmented the nuclear membrane in the first channel, we use it as a mask to measure the intensity in the second channel. .. GENERATED FROM PYTHON SOURCE LINES 149-165 .. code-block:: Python 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() .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_003.png :alt: Second channel (raw), Selection :srcset: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 166-170 Measure the total intensity =========================== The mean intensity is readily available as a region property in a labeled image. .. GENERATED FROM PYTHON SOURCE LINES 170-177 .. code-block:: Python props = measure.regionprops_table( mask.astype(np.uint8), intensity_image=image_t_0_channel_1, properties=('label', 'area', 'intensity_mean'), ) .. GENERATED FROM PYTHON SOURCE LINES 178-179 We may want to check that the value for the total intensity .. GENERATED FROM PYTHON SOURCE LINES 179-182 .. code-block:: Python selection.sum() .. rst-class:: sphx-glr-script-out .. code-block:: none np.uint64(28350) .. GENERATED FROM PYTHON SOURCE LINES 183-184 can be recovered from: .. GENERATED FROM PYTHON SOURCE LINES 184-187 .. code-block:: Python props['area'] * props['intensity_mean'] .. rst-class:: sphx-glr-script-out .. code-block:: none array([28350.]) .. GENERATED FROM PYTHON SOURCE LINES 188-193 Process the entire image sequence ================================= Instead of iterating the workflow for each time point, we process the multidimensional dataset directly (except for the thresholding step). Indeed, most scikit-image functions support nD images. .. GENERATED FROM PYTHON SOURCE LINES 193-200 .. code-block:: Python 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)] .. GENERATED FROM PYTHON SOURCE LINES 201-215 Alternatively, we could compute ``thresh_values`` without using a list comprehension, by reshaping ``smooth_seq`` to make it 2D (where the first dimension still corresponds to time points, but the second and last dimension now contains all pixel values), and applying the thresholding function on the image sequence along its second axis: .. code-block:: python thresh_values = np.apply_along_axis(filters.threshold_otsu, axis=1, arr=smooth_seq.reshape(n_z, -1)) We use the following flat structuring element for morphological computations (``np.newaxis`` is used to prepend an axis of size 1 for time): .. GENERATED FROM PYTHON SOURCE LINES 215-219 .. code-block:: Python footprint = ndi.generate_binary_structure(2, 1)[np.newaxis, ...] footprint .. rst-class:: sphx-glr-script-out .. code-block:: none array([[[False, True, False], [ True, True, True], [False, True, False]]]) .. GENERATED FROM PYTHON SOURCE LINES 220-222 This way, each frame is processed independently (pixels from consecutive frames are never spatial neighbors). .. GENERATED FROM PYTHON SOURCE LINES 222-225 .. code-block:: Python fill_seq = ndi.binary_fill_holes(thresh_seq, structure=footprint) .. GENERATED FROM PYTHON SOURCE LINES 226-250 When clearing objects which touch the image border, we want to make sure that the bottom (first) and top (last) frames are not considered as borders. In this case, the only relevant border is the edge at the greatest (x, y) values. This can be seen in 3D by running the following code: .. code-block:: python 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() .. GENERATED FROM PYTHON SOURCE LINES 250-259 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 260-263 Let us give each mask (corresponding to each time point) a different label, running from 1 to 15. We use ``np.min_scalar_type`` to determine the minimum-size integer dtype needed to represent the number of time points: .. GENERATED FROM PYTHON SOURCE LINES 263-269 .. code-block:: Python 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] .. GENERATED FROM PYTHON SOURCE LINES 270-272 Let us compute the region properties of interest for all these labeled regions. .. GENERATED FROM PYTHON SOURCE LINES 272-296 .. code-block:: Python 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() .. image-sg:: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_004.png :alt: Change in fluorescence intensity at the nuclear envelope :srcset: /auto_examples/applications/images/sphx_glr_plot_fluorescence_nuclear_envelope_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 297-300 Reassuringly, we find the expected result: The total fluorescence intensity at the nuclear envelope increases 1.3-fold in the initial five time points, and then becomes roughly constant. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.661 seconds) .. _sphx_glr_download_auto_examples_applications_plot_fluorescence_nuclear_envelope.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/scikit-image/scikit-image/v0.24.0?filepath=notebooks/auto_examples/applications/plot_fluorescence_nuclear_envelope.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_fluorescence_nuclear_envelope.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_fluorescence_nuclear_envelope.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_fluorescence_nuclear_envelope.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_