From 00a32d18200d45041ba5d569a360595865609694 Mon Sep 17 00:00:00 2001 From: aliceeeeeeeeee Date: Wed, 19 Nov 2025 19:41:28 +0100 Subject: [PATCH 1/2] adding morphometrics notebook --- src/tapenade/analysis/deformation/__init__.py | 31 +- .../analysis/deformation/morphometrics.py | 296 ++++++++ .../deformation_analysis_notebook.ipynb | 31 +- ...tial_features_morphometrics_notebook.ipynb | 663 ++++++++++++++++++ 4 files changed, 1015 insertions(+), 6 deletions(-) create mode 100644 src/tapenade/analysis/deformation/morphometrics.py create mode 100644 src/tapenade/notebooks/spatial_features_morphometrics_notebook.ipynb diff --git a/src/tapenade/analysis/deformation/__init__.py b/src/tapenade/analysis/deformation/__init__.py index bd0cff0..cd12b6c 100644 --- a/src/tapenade/analysis/deformation/__init__.py +++ b/src/tapenade/analysis/deformation/__init__.py @@ -5,7 +5,23 @@ add_ellipsoidal_coefficients, add_ellipsoidal_nature, add_ellipsoidal_nature_bool, - add_anisotropy_coefficient, + add_anisotropy_coefficient +) + +from .morphometrics import ( + process_genetic_field_cellular, + smooth_cellular_map, + process_radial_distance_cellular, + process_axial_distance_cellular, + process_cell_density_sigma, + process_density_gradient, + process_volume_fraction_sigma, + process_nuclear_volume_cellular, + process_nematic_order, + process_ellipsoidal_coeff_cellular, + process_oblate_prolate_cellular, + process_true_strain_maxeig_cellular, + clear_borders, ) from .deformation_quantification import tensors_to_napari_vectors @@ -19,4 +35,17 @@ "add_ellipsoidal_nature_bool", "add_anisotropy_coefficient", "tensors_to_napari_vectors", + "process_genetic_field_cellular", + "smooth_cellular_map", + "process_radial_distance_cellular", + "process_axial_distance_cellular", + "process_cell_density_sigma", + "process_density_gradient", + "process_volume_fraction_sigma", + "process_nuclear_volume_cellular", + "process_nematic_order", + "process_ellipsoidal_coeff_cellular", + "process_oblate_prolate_cellular", + "process_true_strain_maxeig_cellular", + "clear_borders", ] \ No newline at end of file diff --git a/src/tapenade/analysis/deformation/morphometrics.py b/src/tapenade/analysis/deformation/morphometrics.py new file mode 100644 index 0000000..72f7e4f --- /dev/null +++ b/src/tapenade/analysis/deformation/morphometrics.py @@ -0,0 +1,296 @@ +import numpy as np +from skimage.measure import regionprops, label +from scipy.ndimage import distance_transform_edt +from skimage.morphology import binary_erosion +from tapenade.analysis.deformation.additional_regionprops_properties import add_principal_lengths,add_ellipsoidal_coefficients +from tapenade.preprocessing import ( + masked_gaussian_smoothing) +from tapenade.preprocessing._smoothing import _masked_smooth_gaussian + + + +def clear_borders(labels, buffer_size=0, bgval=0, mask=None, *, out=None): + ## from skimage.segmentation + """Clear objects connected to the label image border. + + Parameters + ---------- + labels : (M[, N[, ..., P]]) array of int or bool + Imaging data labels. + buffer_size : int, optional + The width of the border examined. By default, only objects + that touch the outside of the image are removed. + bgval : float or int, optional + Cleared objects are set to this value. + mask : ndarray of bool, same shape as `image`, optional. + Image data mask. Objects in labels image overlapping with + False pixels of mask will be removed. If defined, the + argument buffer_size will be ignored. + out : ndarray + Array of the same shape as `labels`, into which the + output is placed. By default, a new array is created. + + Returns + ------- + out : (M[, N[, ..., P]]) array + Imaging data labels with cleared borders + + """ + if any(buffer_size >= s for s in labels.shape) and mask is None: + # ignore buffer_size if mask + raise ValueError("buffer size may not be greater than labels size") + + if out is None: + out = labels.copy() + + if mask is not None: + err_msg = ( + f"labels and mask should have the same shape but " + f"are {out.shape} and {mask.shape}" + ) + if out.shape != mask.shape: + raise (ValueError, err_msg) + if mask.dtype != bool: + raise TypeError("mask should be of type bool.") + borders = ~mask + else: + # create borders with buffer_size + borders = np.zeros_like(out, dtype=bool) + ext = buffer_size + 1 + slstart = slice(ext) + slend = slice(-ext, None) + slices = [slice(None) for _ in out.shape] + for d in range(out.ndim): + slices[d] = slstart + borders[tuple(slices)] = True + slices[d] = slend + borders[tuple(slices)] = True + slices[d] = slice(None) + + # Re-label, in case we are dealing with a binary out + # and to get consistent labeling + labels, number = label(out, background=0, return_num=True) + + # determine all objects that are connected to borders + borders_indices = np.unique(labels[borders]) + indices = np.arange(number + 1) + # mask all label indices that are connected to borders + label_mask = np.isin(indices, borders_indices) + # create mask for pixels to clear + mask = label_mask[labels.reshape(-1)].reshape(labels.shape) + + # clear border pixels + out[mask] = bgval + + return out + + + +def smooth_cellular_map( + cellular, + props, + mask, + sigma, + +): + #takes a cellular map (one value per segmented nuclei) and apply a gaussian smoothing + + centroids = np.array([prop.centroid for prop in props]) + centroids_data = np.zeros_like(mask, dtype=np.float32) + centroids_data[tuple(centroids.T.astype(int))] = 1 + + smoothed = masked_gaussian_smoothing( + cellular, + mask_for_volume=centroids_data, + sigmas=sigma, + mask=mask, + n_jobs=-1, + ) + return smoothed + + +def process_genetic_field_cellular(labels,genetic_field): + #return cellular map (one value per segmented nuclei) of mean gene expression in that nucleus + props_genetic_field = regionprops(labels, intensity_image=genetic_field) + cellular_genetic_field = np.zeros_like(labels, dtype=np.float32) + for prop in props_genetic_field: + cellular_genetic_field[prop.slice][prop.image] = prop.mean_intensity + return cellular_genetic_field + +def process_radial_distance_cellular( + mask, labels +): + #return cellular map (one value per segmented nuclei) of radial distance to the tissue surface in 3D + radial_distance_cellular = np.zeros_like(mask, dtype=np.float32) + distance_transform = np.zeros_like(mask, dtype=np.float32) + props = regionprops(labels) + for z in range(len(mask)): + distance_transform[z] = distance_transform_edt( + mask[z].astype(bool)) + + for prop in props: + dist = distance_transform[tuple(np.array(prop.centroid).astype(int))] + radial_distance_cellular[prop.slice][prop.image] = dist + return radial_distance_cellular + +def process_axial_distance_cellular( + mask, labels +): + #return cellular map (one value per segmented nuclei) of axial distance along z to the first plane in 3D + axial_distance_cellular = np.zeros_like(mask, dtype=np.float32) + + props = regionprops(labels) + first_x = np.min([prop.centroid[2] for prop in props]) + for prop in props: + dist = prop.centroid[2]-first_x + axial_distance_cellular[prop.slice][prop.image] = dist + + return axial_distance_cellular + +def process_cell_density_sigma( + props, + mask, + sigma, + +): + centroids = np.array([prop.centroid for prop in props]) + centroids_data = np.zeros_like(mask, dtype=np.float32) + centroids_data[tuple(centroids.T.astype(int))] = 1 + + density = masked_gaussian_smoothing( + centroids_data, + sigmas=sigma, + mask=mask, + n_jobs=-1, + ) + return density + +def process_density_gradient( + mask, + density, + positions_on_grid, +): + gradient_field = np.gradient(density) + gradient_field = np.array(gradient_field).transpose(1, 2, 3, 0) + gradient_field[~binary_erosion(mask)] = 0 + + gradient_magnitude_field = np.linalg.norm(gradient_field, axis=-1) + + + gradient_on_grid = gradient_field[tuple(positions_on_grid.T.astype(int))] + + napari_gradient_on_grid = np.zeros((len(positions_on_grid), 2, 3)) + + napari_gradient_on_grid[:, 0] = positions_on_grid + napari_gradient_on_grid[:, 1] = gradient_on_grid + + angles = np.arctan2(*(napari_gradient_on_grid[:, 1, -2:].reshape(-1, 2).T)) + angles = np.arctan2(np.sin(angles - 1), np.cos(angles - 1)) + #gradient_field can be used to compute dot product between cell density gradient and true strain tensor + return gradient_field, gradient_magnitude_field,napari_gradient_on_grid, angles + +def process_volume_fraction_sigma( + mask, + labels, + sigma, +): + volume_data = labels.astype(bool).astype(np.float32) + volume_fraction = masked_gaussian_smoothing( + volume_data, sigmas=sigma, mask=mask, n_jobs=-1 + ) + return volume_fraction + +def process_nuclear_volume_cellular( + props, + mask, +): + centroids_data_volumes = np.zeros_like(mask, dtype=np.float32) + + for prop in props: + centroids_data_volumes[prop.slice][prop.image] = prop.area + + return centroids_data_volumes + +def process_nematic_order( + mask, labels, sigma +): + props =regionprops(labels) + n_points = len(props) + n_dim_tensor = 9 + + sparse_m_tensor = np.zeros((n_points, n_dim_tensor)) + dense_m_tensor = np.zeros((*mask.shape, n_dim_tensor), dtype=np.float32) + dense_centroids = np.zeros(mask.shape) + + for index_label, prop in enumerate(props): + add_principal_lengths(prop, add_principal_vectors=True, scale=(1, 1, 1)) + orientation_vector = prop.principal_vectors[0] + m = np.outer(orientation_vector, orientation_vector) + sparse_m_tensor[index_label, :] = m.flatten() + dense_m_tensor[tuple(np.array(prop.centroid).astype(int))] = m.flatten() + dense_centroids[tuple(np.array(prop.centroid).astype(int))] = 1 + + smoothed_dense_m_tensor = np.zeros_like(dense_m_tensor) + + for i in range(9): + smoothed_dense_m_tensor[..., i] = _masked_smooth_gaussian( + array=dense_m_tensor[..., i], + sigmas=sigma, + mask=mask, + mask_for_volume=dense_centroids.astype(bool), + ) + + masked_flat = smoothed_dense_m_tensor[mask] + q_matrices = masked_flat.reshape(-1, 3, 3) - (np.eye(3) / 3) + if not np.isnan(q_matrices).any(): + eigenvalues = np.linalg.eigvalsh(q_matrices) + s_values = np.max(eigenvalues, axis=1) * 3 / 2 + + s_dense = np.zeros_like(mask, dtype=np.float32) + s_dense[mask] = s_values + + else: + print("NaNs found in q_matrices") + + return s_dense + +def process_ellipsoidal_coeff_cellular(props,mask) : + + ellipsoidal_coefficients = np.zeros_like(mask, dtype=np.float32) + for prop in props: + if not hasattr(prop, "ellipsoidal_coefficients"): + prop = add_ellipsoidal_coefficients(prop, scale=(1,1,1)) + l1, l2, l3 = prop.principal_lengths + score = np.log(l2**2 / (l1 * l3)) + ellipsoidal_coefficients[prop.slice][prop.image] = score + return ellipsoidal_coefficients + +def process_oblate_prolate_cellular(props,mask) : + + ellipsoidal_coefficients_oblate = np.zeros_like(mask, dtype=np.float32) + ellipsoidal_coefficients_prolate = np.zeros_like(mask, dtype=np.float32) + for prop in props: + if not hasattr(prop, "ellipsoidal_coefficients"): + prop = add_ellipsoidal_coefficients(prop, scale=(1,1,1)) + l1, l2, l3 = prop.principal_lengths + ellipsoidal_coefficients_oblate[prop.slice][prop.image] = (l2-l3) / np.sqrt(l1 * l3) + ellipsoidal_coefficients_prolate[prop.slice][prop.image] = (l1-l2) / np.sqrt(l1 * l3) + return ellipsoidal_coefficients_oblate, ellipsoidal_coefficients_prolate + +def process_true_strain_maxeig_cellular( + props,mask +): + # Maximum eigenvalue of the true strain tensor of each nuclei, cellular map (one value per segmented nuclei) + + true_strain_maxeig_cellular = np.zeros_like(mask, dtype=np.float32) + for prop in props: + if not hasattr(prop, "principal_lengths"): + prop = add_principal_lengths(prop, scale=(1, 1, 1)) + principal_lengths = prop.principal_lengths + denominator = np.power(np.prod(principal_lengths), 1 / 3) + max_eig = np.log(principal_lengths[0] / denominator) + + true_strain_maxeig_cellular[prop.slice][prop.image] = max_eig + + return true_strain_maxeig_cellular + \ No newline at end of file diff --git a/src/tapenade/notebooks/deformation_analysis_notebook.ipynb b/src/tapenade/notebooks/deformation_analysis_notebook.ipynb index 0bbd5e7..5d468b3 100644 --- a/src/tapenade/notebooks/deformation_analysis_notebook.ipynb +++ b/src/tapenade/notebooks/deformation_analysis_notebook.ipynb @@ -19,7 +19,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -73,9 +73,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "🔽 First run – downloading data into C:\\Users\\gros\\Desktop\\CODES\\tapenade\\src\\tapenade\\notebooks\\demo_data …\n" + ] + }, + { + "ename": "FileNotFoundError", + "evalue": "[Errno 2] No such file or directory: 'C:\\\\Users\\\\gros\\\\Desktop\\\\CODES\\\\tapenade\\\\src\\\\tapenade\\\\notebooks\\\\demo_data\\\\payload.zip'", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[2], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m path_to_data \u001b[38;5;241m=\u001b[39m \u001b[43mget_path_to_demo_folder\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 3\u001b[0m data \u001b[38;5;241m=\u001b[39m tifffile\u001b[38;5;241m.\u001b[39mimread(path_to_data \u001b[38;5;241m/\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdapi_isotropized.tif\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m 4\u001b[0m mask \u001b[38;5;241m=\u001b[39m tifffile\u001b[38;5;241m.\u001b[39mimread(path_to_data \u001b[38;5;241m/\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmask_def_corr.tif\u001b[39m\u001b[38;5;124m'\u001b[39m)\n", + "File \u001b[1;32m~\\Desktop\\CODES\\tapenade\\src\\tapenade\\data.py:46\u001b[0m, in \u001b[0;36mget_path_to_demo_folder\u001b[1;34m()\u001b[0m\n\u001b[0;32m 43\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m🔽 First run – downloading data into \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mdata_dir\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m …\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m 44\u001b[0m tmp \u001b[38;5;241m=\u001b[39m data_dir \u001b[38;5;241m/\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpayload.zip\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m---> 46\u001b[0m \u001b[43murllib\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrequest\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43murlretrieve\u001b[49m\u001b[43m(\u001b[49m\u001b[43murl\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtmp\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;66;03m# <- simple, std-lib only\u001b[39;00m\n\u001b[0;32m 48\u001b[0m shutil\u001b[38;5;241m.\u001b[39munpack_archive(tmp, data_dir)\n\u001b[0;32m 49\u001b[0m tmp\u001b[38;5;241m.\u001b[39munlink() \u001b[38;5;66;03m# remove archive after unpack\u001b[39;00m\n", + "File \u001b[1;32mc:\\Users\\gros\\Anaconda3\\envs\\env-tapenade\\lib\\urllib\\request.py:251\u001b[0m, in \u001b[0;36murlretrieve\u001b[1;34m(url, filename, reporthook, data)\u001b[0m\n\u001b[0;32m 249\u001b[0m \u001b[38;5;66;03m# Handle temporary file setup.\u001b[39;00m\n\u001b[0;32m 250\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m filename:\n\u001b[1;32m--> 251\u001b[0m tfp \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mopen\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mwb\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[0;32m 252\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m 253\u001b[0m tfp \u001b[38;5;241m=\u001b[39m tempfile\u001b[38;5;241m.\u001b[39mNamedTemporaryFile(delete\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n", + "\u001b[1;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'C:\\\\Users\\\\gros\\\\Desktop\\\\CODES\\\\tapenade\\\\src\\\\tapenade\\\\notebooks\\\\demo_data\\\\payload.zip'" + ] + } + ], "source": [ "path_to_data = get_path_to_demo_folder()\n", "\n", @@ -354,7 +375,7 @@ ], "metadata": { "kernelspec": { - "display_name": "test-tapenade", + "display_name": "env-tapenade", "language": "python", "name": "python3" }, @@ -368,7 +389,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.18" + "version": "3.10.14" } }, "nbformat": 4, diff --git a/src/tapenade/notebooks/spatial_features_morphometrics_notebook.ipynb b/src/tapenade/notebooks/spatial_features_morphometrics_notebook.ipynb new file mode 100644 index 0000000..aba2d5a --- /dev/null +++ b/src/tapenade/notebooks/spatial_features_morphometrics_notebook.ipynb @@ -0,0 +1,663 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Spatial features analysis: morphometrics measurements\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### After clicking on a cell, press \"Shift+Enter\" to run the code, or click on the \"Run\" button in the toolbar above.
\n", + "\n", + "### Replace \"...\" signs with the appropriate path to your data.\n", + "" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import tifffile\n", + "try:\n", + " import napari\n", + " napari_available = True\n", + "except ImportError:\n", + " print(\"Napari is not installed. Some parts of the notebook will not be available.\")\n", + " napari_available = False\n", + "\n", + "from pathlib import Path\n", + "from tapenade import get_path_to_demo_folder\n", + "from tapenade.analysis.deformation.morphometrics import (\n", + " process_genetic_field_cellular,\n", + " smooth_cellular_map,\n", + " process_radial_distance_cellular,\n", + " process_axial_distance_cellular,\n", + " process_cell_density_sigma,\n", + " process_density_gradient,\n", + " process_volume_fraction_sigma,\n", + " process_nuclear_volume_cellular,\n", + " process_nematic_order,\n", + " process_oblate_prolate_cellular,\n", + " process_ellipsoidal_coeff_cellular,\n", + " process_true_strain_maxeig_cellular,\n", + " clear_borders\n", + "\n", + ")\n", + "from ensure import check\n", + "import numpy as np\n", + "from skimage.measure import regionprops" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook presents typical use cases of morphometric measurements\n", + "\n", + "All functions work on 3D (ZYX convention)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Images need to be processed and have a nuclei segmentation and a tissue binary mask, as computed in the prepreocessing notebook.\n", + "\n", + "From there, you can plot maps of :\n", + "1. **field of gene expression**\n", + "2. **spatial position in the tissue: radial or axial scores**\n", + "3. **cell density, cell density gradient magnitude, cell density gradient orientation**\n", + "4. **volume fraction**: fraction of space occupied by nuclei **and nuclear volume**\n", + "5. **nematic order**: local alignment of nuclei true strain tensors.\n", + "6. **nuclei ellipsoidal coefficients, oblate or prolate scores**\n", + "7. **nuclei elongation**: true strain maximum eigenvalue of each nuclei\n", + "\n", + "True strain tensors can be computed using the notebook deformation_analysis_notebook.ipynb\n", + "\n", + "\n", + "Definitions of all metrics are found in Gros, Vanaret, Dunsing-Eichenhauer et al, eLife 2025.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- This notebook computes and plots different 3D maps on napari. Each feature is shown in an independant napari viewer. If you want to plot all of them in one viewer, just comment or uncomment the lines `viewer=napari.Viewer()` (which creates a new viewer) and `napari.run()`\n", + "- To adjust the min and max values on the 3D maps, add it as a parameter in the `add_image` function : `viewer.add_image(feature, contrast_limits = [vmin,vmax]`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 0. Loading the data and choosing parameters (mandatory before computing any feature)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "path_to_data = get_path_to_demo_folder()\n", + "\n", + "data = tifffile.imread(path_to_data / 'dapi_isotropized.tif')\n", + "mask = tifffile.imread(path_to_data / 'mask_def_corr.tif')\n", + "labels = tifffile.imread(path_to_data / 'labels_def_corr.tif')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Load your data.\n", + "\n", + "For simpler interpretation of 3D fields, your mask and labels should have an isotropic pixel size, for example [1,1,1] µm/pix (ZYX). Please input 3D data with only one channel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "path_to_data = ...\n", + "\n", + "data = tifffile.imread(Path(path_to_data) / \"...tif\")\n", + "mask = tifffile.imread(Path(path_to_data) / \"...tif\")\n", + "labels = tifffile.imread(Path(path_to_data) / \"...tif\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# smoothing parameter\n", + "sigma=20" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Adjust labels" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "labels = clear_borders(labels) #removing labels touching the border because their shape is not reliable\n", + "\n", + "props = regionprops(labels)\n", + "check(labels.shape).equals(mask.shape).or_raise(\n", + " Exception,\n", + " f\"The shape of the selected segmentation {labels.shape} does not match the shape of the mask {mask.shape}. \\n\",\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer = napari.view_image(data, name='raw data')\n", + "viewer.add_image(mask, name='mask')\n", + "viewer.add_labels(labels, name='labels')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def get_napari_angles_cmap():\n", + " cmap_class = napari.utils.colormaps.colormap.Colormap\n", + "\n", + " colors_strs = [\n", + " (195, 163, 75),\n", + " (163, 103, 44),\n", + " (135, 64, 55),\n", + " (115, 57, 87),\n", + " (92, 83, 139),\n", + " (79, 136, 185),\n", + " (116, 187, 205),\n", + " (180, 222, 198),\n", + " (214, 216, 147),\n", + " (195, 163, 75),\n", + " (163, 103, 44),\n", + " (135, 64, 55),\n", + " (115, 57, 87),\n", + " (92, 83, 139),\n", + " (79, 136, 185),\n", + " (116, 187, 205),\n", + " (180, 222, 198),\n", + " (214, 216, 147),\n", + " ]\n", + "\n", + " for i, elem in enumerate(colors_strs):\n", + " colors_strs[i] = tuple([x / 255 for x in elem])\n", + "\n", + " cmap = cmap_class(\n", + " colors_strs,\n", + " controls=np.linspace(0, 1, 1 + len(colors_strs)),\n", + " interpolation=\"zero\",\n", + " )\n", + "\n", + " return cmap" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Field of gene expression" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Images of gene expression markers need to be normalized in 3D to remove intensity gradients linked with optical artifacts, see preprocessing_notebook.ipynb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 1) Sanity checks on image dimensions\n", + "check(data.shape).equals(labels.shape).or_raise(\n", + " Exception,\n", + " f\"The shape of the selected data {data.shape} does not match the shape of the labels {labels.shape}. \\n\",\n", + ")\n", + "check(data.shape).equals(mask.shape).or_raise(\n", + " Exception,\n", + " f\"The shape of the selected data {data.shape} does not match the shape of the mask {mask.shape}. \\n\",\n", + ")\n", + "\n", + "# 2) Check for NaN values in the data\n", + "if np.isnan(data).any():\n", + " data[mask == 0] = 0\n", + "\n", + "check(np.isnan(data).any()).equals(False).or_raise(\n", + " Exception,\n", + " \"The image contains NaN values. \\n\",\n", + ")\n", + "\n", + "# 3) Compute the genetic field\n", + "genetic_field_cellular = process_genetic_field_cellular(\n", + "labels,\n", + "data\n", + ")\n", + "genetic_field_smooth = smooth_cellular_map(\n", + "genetic_field_cellular, props, mask, sigma\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer=napari.view_labels(labels, name='labels')\n", + "viewer.add_image(data)\n", + "viewer.add_image(genetic_field_cellular, name='genetic field cellular',colormap='inferno')\n", + "viewer.add_image(genetic_field_smooth, name='genetic field smooth',colormap='inferno')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Spatial position in the tissue\n", + "\n", + "distance to border (radial) or distance along longest axis = AP distance for gastruloids (axial)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "radial_distance_cellular = process_radial_distance_cellular(\n", + " mask, labels\n", + ")\n", + "axial_distance_cellular = process_axial_distance_cellular(\n", + " mask, labels\n", + ")\n", + "\n", + "axial_distance_smooth = smooth_cellular_map(\n", + " axial_distance_cellular, props, mask, sigma\n", + ")\n", + "radial_distance_smooth = smooth_cellular_map(\n", + " radial_distance_cellular, props, mask, sigma\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer=napari.Viewer()\n", + "viewer.add_image(radial_distance_cellular, name='radial distance cellular',colormap='inferno')\n", + "viewer.add_image(radial_distance_smooth, name='radial distance smooth',colormap='inferno')\n", + "viewer.add_image(axial_distance_cellular, name='axial distance cellular',colormap='inferno')\n", + "viewer.add_image(axial_distance_smooth, name='axial distance smooth',colormap='inferno')\n", + "napari.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Cell density, cell density gradient magnitude and direction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compute grid to plot cell density gradient" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "z_spacing = 10 #to tune\n", + "xy_spacing = 10 #to tune\n", + "\n", + "positions_on_grid = (\n", + " np.mgrid[\n", + " [slice(0, labels.shape[0], z_spacing)]\n", + " + [slice(0, dim, xy_spacing) for dim in labels.shape[1:]]\n", + " ]\n", + " .reshape(labels.ndim, -1)\n", + " .T\n", + " )\n", + "positions_on_grid = positions_on_grid[\n", + " np.where(\n", + " mask[\n", + " tuple(\n", + " [\n", + " positions_on_grid[:, 0],\n", + " positions_on_grid[:, 1],\n", + " positions_on_grid[:, 2],\n", + " ]\n", + " )\n", + " ]\n", + " )\n", + " ]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compute cell density map and gradient" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "centroids = np.array([prop.centroid for prop in props])\n", + "density = process_cell_density_sigma(props,mask,sigma)\n", + "gradient_field, gradient_magnitude_field,napari_gradient_on_grid, density_gradient_angles = process_density_gradient(\n", + " mask,\n", + " density,\n", + " positions_on_grid,\n", + " )\n", + "lengths = np.linalg.norm(napari_gradient_on_grid[:, 1], axis=1)\n", + "napari_gradient_on_grid[:, 1] = (\n", + " napari_gradient_on_grid[:, 1] / np.median(lengths[:, None]) * 5 #*5 is just for visualization purposes\n", + ")\n", + "napari_gradient_on_grid_proj = napari_gradient_on_grid[:, :, 1:]\n", + "cmap_vectors_gradient = get_napari_angles_cmap()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize cell density maps" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer = napari.Viewer()\n", + "viewer.add_image(density,colormap='inferno', name='cell density')\n", + "viewer.add_image(gradient_magnitude_field,colormap='inferno', name='density gradient magnitude')\n", + "napari.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualize cell density gradient" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer = napari.Viewer()\n", + "viewer.add_vectors(\n", + " napari_gradient_on_grid,\n", + " name=\"density_gradient\",\n", + " properties={\"angles\": density_gradient_angles},\n", + " vector_style=\"triangle\",\n", + " blending=\"translucent\",\n", + " opacity=1,\n", + " edge_colormap=cmap_vectors_gradient,\n", + " length=5,\n", + " edge_width=5,\n", + ")\n", + "viewer.add_vectors(\n", + " napari_gradient_on_grid_proj,\n", + " name=\"density_gradient_proj\",\n", + " properties={\"angles\": density_gradient_angles},\n", + " vector_style=\"triangle\",\n", + " blending=\"translucent\",\n", + " opacity=1,\n", + " edge_colormap=cmap_vectors_gradient,\n", + " length=7,\n", + " edge_width=2,\n", + ")\n", + "napari.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Volume fraction and nuclear volume" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "volume_fraction = process_volume_fraction_sigma(\n", + " mask,\n", + " labels,\n", + " sigma\n", + ")\n", + "nuclear_volume_cellular = process_nuclear_volume_cellular(\n", + " props,mask)\n", + "nuclear_volume_smooth = smooth_cellular_map(\n", + " nuclear_volume_cellular, props, mask, sigma\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer = napari.Viewer()\n", + "viewer.add_image(volume_fraction,colormap='inferno', name='volume fraction')\n", + "viewer.add_image(nuclear_volume_cellular,colormap='inferno', name='nuclear volume cellular')\n", + "viewer.add_image(nuclear_volume_smooth,colormap='inferno', name='nuclear volume smooth')\n", + "napari.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 5. Nematic order" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nematic_order = process_nematic_order(\n", + " mask,\n", + " labels,\n", + " sigma\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer=napari.Viewer()\n", + "viewer.add_image(nematic_order, name='nematic order',colormap='inferno')\n", + "napari.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 6. Nuclei ellipsoidal coefficients, oblate or prolate scores" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Two ways of computing it :\n", + "- with two different scores, one oblate, one prolate (>0)\n", + "- with a unique score between -inf and +inf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ellipsoidal_coefficients_oblate, ellipsoidal_coefficients_prolate = process_oblate_prolate_cellular(\n", + " props,mask\n", + ")\n", + "ellipsoidal_coefficients = process_ellipsoidal_coeff_cellular(\n", + " props, mask\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compute min and max values for visualization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "max_oblate_prolate = np.nanmax(np.concatenate((ellipsoidal_coefficients_oblate.flatten(),ellipsoidal_coefficients_prolate.flatten())))\n", + "max_extrema_ellipsoidal = max(abs(np.nanmin(ellipsoidal_coefficients.flatten())),abs(np.nanmax(ellipsoidal_coefficients.flatten()))) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer=napari.Viewer()\n", + "viewer.add_image(ellipsoidal_coefficients_oblate, name='ellipsoidal_coefficients_oblate',colormap='inferno',contrast_limits=(0,max_oblate_prolate))\n", + "viewer.add_image(ellipsoidal_coefficients_prolate, name='ellipsoidal_coefficients_prolate',colormap='inferno',contrast_limits=(0,max_oblate_prolate))\n", + "viewer.add_image(ellipsoidal_coefficients, name='ellipsoidal_coefficients',colormap='PiYG',contrast_limits=(-max_extrema_ellipsoidal,max_extrema_ellipsoidal))\n", + "napari.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 7. Nuclei elongation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "true_strain_maxeig_cellular = process_true_strain_maxeig_cellular(\n", + " props,mask)\n", + "\n", + "true_strain_maxeig_smooth = smooth_cellular_map(\n", + " true_strain_maxeig_cellular, props, mask, sigma)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "viewer = napari.Viewer()\n", + "viewer.add_image(true_strain_maxeig_cellular, name='true_strain_maxeig_cellular',colormap='inferno')\n", + "viewer.add_image(true_strain_maxeig_smooth, name='true_strain_maxeig_smooth',colormap='inferno')\n", + "napari.run()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "env-tapenade", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From d3df14b0d54724fd3781186985659df3c3f387db Mon Sep 17 00:00:00 2001 From: aliceeeeeeeeee Date: Wed, 19 Nov 2025 20:28:25 +0100 Subject: [PATCH 2/2] fix morpho measurements --- ...tial_features_morphometrics_notebook.ipynb | 100 ++++++++++-------- 1 file changed, 58 insertions(+), 42 deletions(-) diff --git a/src/tapenade/notebooks/spatial_features_morphometrics_notebook.ipynb b/src/tapenade/notebooks/spatial_features_morphometrics_notebook.ipynb index aba2d5a..190238d 100644 --- a/src/tapenade/notebooks/spatial_features_morphometrics_notebook.ipynb +++ b/src/tapenade/notebooks/spatial_features_morphometrics_notebook.ipynb @@ -58,7 +58,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This notebook presents typical use cases of morphometric measurements\n", + "This notebook presents typical use cases of morphometric measurements.\n", "\n", "All functions work on 3D (ZYX convention)" ] @@ -106,6 +106,13 @@ "## 0. Loading the data and choosing parameters (mandatory before computing any feature)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Load example data" + ] + }, { "cell_type": "code", "execution_count": null, @@ -125,7 +132,8 @@ "source": [ "Load your data.\n", "\n", - "For simpler interpretation of 3D fields, your mask and labels should have an isotropic pixel size, for example [1,1,1] µm/pix (ZYX). Please input 3D data with only one channel" + "- For simpler interpretation of 3D fields, your mask and labels should have an isotropic pixel size, for example [1,1,1] µm/pix (ZYX)\n", + "- Please input 3D data with only one channel" ] }, { @@ -181,14 +189,10 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "viewer = napari.view_image(data, name='raw data')\n", - "viewer.add_image(mask, name='mask')\n", - "viewer.add_labels(labels, name='labels')" + "Visualized image, mask and segmentation" ] }, { @@ -197,40 +201,9 @@ "metadata": {}, "outputs": [], "source": [ - "def get_napari_angles_cmap():\n", - " cmap_class = napari.utils.colormaps.colormap.Colormap\n", - "\n", - " colors_strs = [\n", - " (195, 163, 75),\n", - " (163, 103, 44),\n", - " (135, 64, 55),\n", - " (115, 57, 87),\n", - " (92, 83, 139),\n", - " (79, 136, 185),\n", - " (116, 187, 205),\n", - " (180, 222, 198),\n", - " (214, 216, 147),\n", - " (195, 163, 75),\n", - " (163, 103, 44),\n", - " (135, 64, 55),\n", - " (115, 57, 87),\n", - " (92, 83, 139),\n", - " (79, 136, 185),\n", - " (116, 187, 205),\n", - " (180, 222, 198),\n", - " (214, 216, 147),\n", - " ]\n", - "\n", - " for i, elem in enumerate(colors_strs):\n", - " colors_strs[i] = tuple([x / 255 for x in elem])\n", - "\n", - " cmap = cmap_class(\n", - " colors_strs,\n", - " controls=np.linspace(0, 1, 1 + len(colors_strs)),\n", - " interpolation=\"zero\",\n", - " )\n", - "\n", - " return cmap" + "viewer = napari.view_image(data, name='raw data')\n", + "viewer.add_image(mask, name='mask')\n", + "viewer.add_labels(labels, name='labels')" ] }, { @@ -398,6 +371,49 @@ "Compute cell density map and gradient" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# colormap for napari vector fields\n", + "def get_napari_angles_cmap():\n", + " cmap_class = napari.utils.colormaps.colormap.Colormap\n", + "\n", + " colors_strs = [\n", + " (195, 163, 75),\n", + " (163, 103, 44),\n", + " (135, 64, 55),\n", + " (115, 57, 87),\n", + " (92, 83, 139),\n", + " (79, 136, 185),\n", + " (116, 187, 205),\n", + " (180, 222, 198),\n", + " (214, 216, 147),\n", + " (195, 163, 75),\n", + " (163, 103, 44),\n", + " (135, 64, 55),\n", + " (115, 57, 87),\n", + " (92, 83, 139),\n", + " (79, 136, 185),\n", + " (116, 187, 205),\n", + " (180, 222, 198),\n", + " (214, 216, 147),\n", + " ]\n", + "\n", + " for i, elem in enumerate(colors_strs):\n", + " colors_strs[i] = tuple([x / 255 for x in elem])\n", + "\n", + " cmap = cmap_class(\n", + " colors_strs,\n", + " controls=np.linspace(0, 1, 1 + len(colors_strs)),\n", + " interpolation=\"zero\",\n", + " )\n", + "\n", + " return cmap" + ] + }, { "cell_type": "code", "execution_count": null,