From 92b8ebe3960b1e021e71166491c3d7076c8f0caf Mon Sep 17 00:00:00 2001
From: Guillaume Cerutti <guillaume.cerutti@inria.fr>
Date: Wed, 6 Mar 2024 11:46:02 +0100
Subject: [PATCH] add command line script
---
script/nuclei_quantification.py | 141 ++++++++++++++++++++++++++++++++
1 file changed, 141 insertions(+)
create mode 100644 script/nuclei_quantification.py
diff --git a/script/nuclei_quantification.py b/script/nuclei_quantification.py
new file mode 100644
index 0000000..0f6bc41
--- /dev/null
+++ b/script/nuclei_quantification.py
@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import argparse
+import logging
+
+import numpy as np
+import pandas as pd
+
+import matplotlib.pyplot as plt
+
+from timagetk import SpatialImage
+from timagetk.io import imread
+from timagetk.algorithms.peak_detection import detect_nuclei
+from timagetk.algorithms.signal_quantification import quantify_nuclei_signal_intensity
+
+
+if __name__ == "__main__":
+ logging.getLogger().setLevel(logging.INFO)
+
+ parser = argparse.ArgumentParser()
+ parser.add_argument('-d', '--dirname', help='Path to the data directory', default="../data/")
+ parser.add_argument('-f', '--filename', help='Image name identifier', required=True)
+ parser.add_argument('-c', '--channel-names', help='List of channel names in the same order as in the microscopy image', nargs='+', type=str)
+ parser.add_argument('-w', '--channel-weights', help='List of weights associated to each channel to compute the detection image', nargs='+', type=float)
+ parser.add_argument('-q', '--quantification-channels', help='List of channel names on which to quantify nuclei intensities', nargs='+', type=str)
+
+ parser.add_argument('-t', '--detection-threshold', default=100, help="Threshold used for nuclei detection", type=float)
+ parser.add_argument('-rmin', '--detection-radius-min', default=0.8, help="Min radius used for signal quantification", type=float)
+ parser.add_argument('-rmax', '--detection-radius-max', default=3, help="Max radius used for signal quantification", type=float)
+ parser.add_argument('-r', '--quantification-radius', default=2.5, help="Radius used for signal quantification", type=float)
+ args = parser.parse_args()
+
+ dirname = args.dirname # "../data/"
+ filename = args.filename # "230315_GK_86_9_T1_5"
+ channel_names = args.channel_names
+
+ image_filename = f"{dirname}/{filename}.tif"
+ img = imread(image_filename, channel_names=channel_names)
+ channel_names = img.channel_names
+
+ logging.info(f"--> Reading image {image_filename}")
+ if args.quantification_channels is not None:
+ quantification_channels = args.quantification_channels
+ assert all([c in channel_names for c in quantification_channels])
+ else:
+ quantification_channels = channel_names
+
+ detection_threshold = args.detection_threshold # 100
+ detection_radius_range = (args.detection_radius_min, args.detection_radius_max) # (0.8, 3)
+ quantification_radius = args.quantification_radius # 2.5
+
+ logging.info(f"--> Computing merged image for detection")
+ if args.channel_weights is not None:
+ channel_weights = args.channel_weights
+ assert len(channel_weights) == len(channel_names)
+ else:
+ channel_weights = [1. for c in channel_names]
+
+ nuclei_img = SpatialImage(
+ np.sum([
+ w*img[c].get_array()
+ for c, w in zip(channel_names, channel_weights)
+ ], axis=0),
+ voxelsize=img.voxelsize
+ )
+
+ logging.info(f"--> Detecting nuclei")
+ nuclei_points = detect_nuclei(
+ nuclei_img,
+ threshold=detection_threshold,
+ radius_range=detection_radius_range
+ )
+
+ logging.info(f"--> Quantifying nuclei signals for {quantification_channels}")
+ nuclei_signals = {
+ signal_name: quantify_nuclei_signal_intensity(
+ img[signal_name],
+ nuclei_points,
+ nuclei_sigma=quantification_radius
+ )
+ for signal_name in quantification_channels
+ }
+
+ logging.info(f"--> Displaying results")
+ figure = plt.figure(figsize=(10*len(quantification_channels), 10))
+
+ extent = (
+ -img.voxelsize[2]/2,
+ img.extent[2]+img.voxelsize[2]/2,
+ img.extent[1]+img.voxelsize[1]/2,
+ -img.voxelsize[1]/2,
+ )
+
+ for i, signal_name in enumerate(quantification_channels):
+ figure.add_subplot(1, 2, i+1)
+ figure.gca().imshow(
+ img[signal_name].get_array().max(axis=0),
+ extent=extent,
+ )
+
+ figure.gca().scatter(
+ nuclei_points[:, 0],
+ nuclei_points[:, 1],
+ c=nuclei_signals[signal_name],
+ cmap='Reds'
+ )
+
+ figure.gca().set_title(f"{signal_name} nuclei signal intensity", size=18)
+
+ figure.tight_layout()
+ figure_filename = f"{dirname}/{filename}_nuclei_quantification.png"
+ figure.savefig(figure_filename)
+
+ figure = plt.figure(figsize=(10, 10))
+
+ figure.gca().scatter(
+ nuclei_signals[quantification_channels[0]],
+ nuclei_signals[quantification_channels[1]],
+ c=nuclei_points[:, 1],
+ cmap="Reds",
+ )
+ figure.gca().set_xlabel(quantification_channels[0], size=18)
+ figure.gca().set_ylabel(quantification_channels[1], size=18)
+
+ figure.tight_layout()
+ figure_filename = f"{dirname}/{filename}_nuclei_{quantification_channels[0]}_{quantification_channels[1]}_signals.png"
+ figure.savefig(figure_filename)
+
+ nuclei_data = {
+ 'id': range(len(nuclei_points)),
+ }
+ nuclei_data.update({
+ f'center_{dim}': nuclei_points[:, i]
+ for i, dim in enumerate('xyz')
+ })
+ nuclei_data.update(nuclei_signals)
+
+ nuclei_df = pd.DataFrame(nuclei_data)
+ nuclei_filename = f"{dirname}/{filename}_nuclei_data.csv"
+ nuclei_df.to_csv(nuclei_filename, index=False)
--
GitLab