add command line script

script/nuclei_quantification.py

+#!/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)