Update plotting scripts

Author: constantinpape

scripts/figures/plot_fig2.py

@@ -1,19 +1,20 @@
 import argparse
 import os
+from glob import glob
 
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
 
+from util import literature_reference_values
+
 png_dpi = 300
 
 
-def fig_02c(save_path, plot=False):
+def fig_02c(save_path, plot=False, all_versions=False):
     """Scatter plot showing the precision, recall, and F1-score of SGN (distance U-Net, manual),
     IHC (distance U-Net, manual), and synapse detection (U-Net).
     """
-    setting = ["U-Net", "manual", "U-Net", "manual", "U-Net"]
-
     # precision, recall, f1-score
     sgn_unet = [0.887, 0.88, 0.884]
     sgn_annotator = [0.95, 0.849, 0.9]
@@ -25,19 +26,32 @@ def fig_02c(save_path, plot=False):
     ihc_annotator = [0.958, 0.956, 0.957]
     syn_unet = [0.931, 0.905, 0.918]
 
+    # This is the version with IHC v4b segmentation:
+    # 4th version of the network with optimized segmentation params
     version_1 = [sgn_unet, sgn_annotator, ihc_v4b, ihc_annotator, syn_unet]
-    settings_1 = ["SGN_v2", "manual", "IHC_v4b", "manual", "U-Net"]
+    settings_1 = ["automatic", "manual", "automatic", "manual", "automatic"]
 
+    # This is the version with IHC v4c segmentation:
+    # 4th version of the network with optimized segmentation params and split of falsely merged IHCs
     version_2 = [sgn_unet, sgn_annotator, ihc_v4c, ihc_annotator, syn_unet]
-    settings_2 = ["SGN_v2", "manual", "IHC_v4c", "manual", "U-Net"]
+    settings_2 = ["automatic", "manual", "automatic", "manual", "automatic"]
 
+    # This is the version with IHC v4c + filter segmentation:
+    # 4th version of the network with optimized segmentation params and split of falsely merged IHCs
+    # + filtering out IHCs with zero mapped synapses.
     version_3 = [sgn_unet, sgn_annotator, ihc_v4c_filter, ihc_annotator, syn_unet]
-    settings_3 = ["SGN_v2", "manual", "IHC_v4c_filter", "manual", "U-Net"]
+    settings_3 = ["automatic", "manual", "automatic", "manual", "automatic"]
 
-    versions = [version_1, version_2, version_3]
-    settings = [settings_1, settings_2, settings_3]
-    save_suffix = ["v4b", "v4c", "v4c_filter"]
-    save_paths = [save_path + i for i in save_suffix]
+    if all_versions:
+        versions = [version_1, version_2, version_3]
+        settings = [settings_1, settings_2, settings_3]
+        save_suffix = ["_v4b", "_v4c", "_v4c_filter"]
+        save_paths = [save_path + i for i in save_suffix]
+    else:
+        versions = [version_2]
+        settings = [settings_2]
+        save_suffix = ["_v4c"]
+        save_paths = [save_path + i for i in save_suffix]
 
     for version, setting, save_path in zip(versions, settings, save_paths):
         precision = [i[0] for i in version]
@@ -84,22 +98,37 @@ def fig_02c(save_path, plot=False):
             plt.close()
 
 
-def fig_02d_01(save_path, plot=False):
+# Load the synapse counts for all IHCs from the relevant tables.
+def _load_ribbon_synapse_counts():
+    tables = glob("ihc_counts/*M_LR*.tsv")
+    syn_counts = []
+    for tab in tables:
+        x = pd.read_csv(tab, sep="\t")
+        syn_counts.extend(x["synapse_count"].values.tolist())
+    return syn_counts
+
+
+def fig_02d_01(save_path, plot=False, all_versions=False, plot_average_ribbon_synapses=False):
     """Box plot showing the counts for SGN and IHC per (mouse) cochlea in comparison to literature values.
     """
     main_tick_size = 16
     main_label_size = 24
 
     rows = 1
-    columns = 2
+    columns = 3 if plot_average_ribbon_synapses else 2
 
     sgn_values = [11153, 11398, 10333, 11820]
     ihc_v4b_values = [836, 808, 796, 901]
     ihc_v4c_values = [712, 710, 721, 675]
     ihc_v4c_filtered_values = [562, 647, 626, 628]
 
-    ihc_list = [ihc_v4b_values, ihc_v4c_values, ihc_v4c_filtered_values]
-    suffixes = ["_v4b", "_v4c", "_v4c_filtered"]
+    if all_versions:
+        ihc_list = [ihc_v4b_values, ihc_v4c_values, ihc_v4c_filtered_values]
+        suffixes = ["_v4b", "_v4c", "_v4c_filtered"]
+        assert not plot_average_ribbon_synapses
+    else:
+        ihc_list = [ihc_v4c_values]
+        suffixes = ["_v4c"]
 
     for (ihc_values, suffix) in zip(ihc_list, suffixes):
         fig, axes = plt.subplots(rows, columns, figsize=(columns*4, rows*4))
@@ -126,34 +155,46 @@ def fig_02d_01(save_path, plot=False):
         xmin = 0.5
         xmax = 1.5
         ax[0].set_xlim(xmin, xmax)
-        upper_y = 12000
-        lower_y = 10000
+        lower_y, upper_y = literature_reference_values("SGN")
         ax[0].hlines([lower_y, upper_y], xmin, xmax)
-        ax[0].text(1, upper_y + 100, "literature reference (WIP)", color="C0", fontsize=main_tick_size, ha="center")
+        ax[0].text(1, upper_y + 100, "literature", color="C0", fontsize=main_tick_size, ha="center")
         ax[0].fill_between([xmin, xmax], lower_y, upper_y, color="C0", alpha=0.05, interpolate=True)
 
         ylim0 = 500
         ylim1 = 950
         y_ticks = [i for i in range(500, 900 + 1, 100)]
 
         ax[1].set_xticklabels(["IHC"], fontsize=main_label_size)
-
         ax[1].set_yticks(y_ticks)
         ax[1].set_yticklabels(y_ticks, rotation=0, fontsize=main_tick_size)
         ax[1].set_ylim(ylim0, ylim1)
-        ax[1].yaxis.tick_right()
-        ax[1].yaxis.set_ticks_position("right")
+        if not plot_average_ribbon_synapses:
+            ax[1].yaxis.tick_right()
+            ax[1].yaxis.set_ticks_position("right")
 
         # set range of literature values
         xmin = 0.5
         xmax = 1.5
+        lower_y, upper_y = literature_reference_values("IHC")
         ax[1].set_xlim(xmin, xmax)
-        upper_y = 850
-        lower_y = 780
         ax[1].hlines([lower_y, upper_y], xmin, xmax)
-        ax[1].text(1, lower_y - 10, "literature reference (WIP)", color="C0", fontsize=main_tick_size, ha="center")
+        ax[1].text(1, upper_y + 20, "literature", color="C0", fontsize=main_tick_size, ha="center")
         ax[1].fill_between([xmin, xmax], lower_y, upper_y, color="C0", alpha=0.05, interpolate=True)
 
+        if plot_average_ribbon_synapses:
+            ribbon_synapse_counts = _load_ribbon_synapse_counts()
+            # ylim0 = 4.9
+            # ylim1 = 25.1
+            y_ticks = [0, 10, 20, 30, 40, 50]
+
+            ax[2].boxplot(ribbon_synapse_counts)
+            ax[2].set_xticklabels(["Ribbon Syn. per IHC"], fontsize=main_label_size)
+            ax[2].set_yticks(y_ticks)
+            ax[2].set_yticklabels(y_ticks, rotation=0, fontsize=main_tick_size)
+            # ax[2].set_ylim(ylim0, ylim1)
+
+            # TODO range of values from literature.
+
         plt.tight_layout()
         plt.savefig(save_path_new, dpi=png_dpi)
 
@@ -249,12 +290,17 @@ def main():
 
     os.makedirs(args.figure_dir, exist_ok=True)
 
-    fig_02c(save_path=os.path.join(args.figure_dir, "fig_02c"), plot=args.plot)
-    fig_02d_01(save_path=os.path.join(args.figure_dir, "fig_02d_01"), plot=args.plot)
-    fig_02d_02(save_path=os.path.join(args.figure_dir, "fig_02d_02"), plot=args.plot)
+    # Panel C: Evaluation of the segmentation results:
+    fig_02c(save_path=os.path.join(args.figure_dir, "fig_02c"), plot=args.plot, all_versions=False)
+
+    # Panel D: The number of SGNs, IHCs and average number of ribbon synapses per IHC
+    fig_02d_01(save_path=os.path.join(args.figure_dir, "fig_02d"), plot=args.plot, plot_average_ribbon_synapses=True)
+
+    # Alternative version of synapse distribution for panel D.
+    # fig_02d_02(save_path=os.path.join(args.figure_dir, "fig_02d_02"), plot=args.plot)
     # fig_02d_02(save_path=os.path.join(args.figure_dir, "fig_02d_02_v4c"), filter_zeros=False, plot=plot)
     # fig_02d_02(save_path=os.path.join(args.figure_dir, "fig_02d_02_v4c_filtered"), filter_zeros=True, plot=plot)
-    fig_02d_02(save_path=os.path.join(args.figure_dir, "fig_02d_02_v4b"), filter_zeros=True, plot=args.plot)
+    # fig_02d_02(save_path=os.path.join(args.figure_dir, "fig_02d_02_v4b"), filter_zeros=True, plot=args.plot)
 
 
 if __name__ == "__main__":

scripts/figures/plot_fig3.py

@@ -0,0 +1,169 @@
+import argparse
+import os
+import imageio.v3 as imageio
+from glob import glob
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from matplotlib import cm, colors
+
+from util import sliding_runlength_sum, frequency_mapping
+
+INPUT_ROOT = "/home/pape/Work/my_projects/flamingo-tools/scripts/M_LR_000227_R/scale3/frequency_mapping"
+
+png_dpi = 300
+
+
+def fig_03a(save_path):
+    import napari
+
+    path = os.path.join(INPUT_ROOT, "frequencies_IHC_v4c.tif")
+    vol = imageio.imread(path)
+
+    # Create the colormap
+    fig, ax = plt.subplots(figsize=(6, 1.3))
+    fig.subplots_adjust(bottom=0.5)
+
+    freq_min = np.min(np.nonzero(vol))
+    freq_max = vol.max()
+    norm = colors.Normalize(vmin=freq_min, vmax=freq_max, clip=True)
+    cmap = plt.get_cmap("viridis")
+
+    cb = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), cax=ax, orientation="horizontal")
+    cb.set_label("Frequency [kHz]")
+    plt.title("Tonotopic Mapping: IHCs")
+    plt.tight_layout()
+    out_path = os.path.join(save_path)
+    plt.savefig(out_path)
+
+    # Show the image in napari for rendering.
+    v = napari.Viewer()
+    v.add_image(vol, colormap="viridis")
+    napari.run()
+
+
+def fig_03b(save_path):
+    import napari
+
+    path = os.path.join(INPUT_ROOT, "frequencies_SGN_v2.tif")
+    vol = imageio.imread(path)
+
+    # Create the colormap
+    fig, ax = plt.subplots(figsize=(6, 1.3))
+    fig.subplots_adjust(bottom=0.5)
+
+    freq_min = np.min(np.nonzero(vol))
+    freq_max = vol.max()
+    norm = colors.Normalize(vmin=freq_min, vmax=freq_max, clip=True)
+    cmap = plt.get_cmap("viridis")
+
+    cb = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), cax=ax, orientation="horizontal")
+    cb.set_label("Frequency [kHz]")
+    plt.title("Tonotopic Mapping: SGNs")
+    plt.tight_layout()
+    out_path = os.path.join(save_path)
+    plt.savefig(out_path)
+
+    # Show the image in napari for rendering.
+    v = napari.Viewer()
+    v.add_image(vol, colormap="viridis")
+    napari.run()
+
+
+def fig_03c_rl(save_path, plot=False):
+    tables = glob("./ihc_counts/ihc_count_M_LR*.tsv")
+    fig, ax = plt.subplots(figsize=(8, 4))
+
+    width = 50  # micron
+
+    for tab_path in tables:
+        # TODO map to alias
+        alias = os.path.basename(tab_path)[10:-4].replace("_", "").replace("0", "")
+        tab = pd.read_csv(tab_path, sep="\t")
+        run_length = tab["run_length"].values
+        syn_count = tab["synapse_count"].values
+
+        # Compute the running sum of 10 micron.
+        run_length, syn_count_running = sliding_runlength_sum(run_length, syn_count, width=width)
+        ax.plot(run_length, syn_count_running, label=alias)
+
+    ax.set_xlabel("Length [µm]")
+    ax.set_ylabel("Synapse Count")
+    ax.set_title(f"Ribbon Syn. per IHC: Runnig sum @ {width} µm")
+    ax.legend(title="cochlea")
+    plt.tight_layout()
+
+    plt.savefig(save_path, bbox_inches="tight", pad_inches=0.1, dpi=png_dpi)
+    if plot:
+        plt.show()
+    else:
+        plt.close()
+
+
+def fig_03c_octave(save_path, plot=False):
+    tables = glob("./ihc_counts/ihc_count_M_LR*.tsv")
+
+    result = {"cochlea": [], "octave_band": [], "value": []}
+    for tab_path in tables:
+        # TODO map to alias
+        alias = os.path.basename(tab_path)[10:-4].replace("_", "").replace("0", "")
+        tab = pd.read_csv(tab_path, sep="\t")
+        freq = tab["frequency"].values
+        syn_count = tab["synapse_count"].values
+
+        # Compute the running sum of 10 micron.
+        octave_binned = frequency_mapping(freq, syn_count, animal="mouse")
+
+        result["cochlea"].extend([alias] * len(octave_binned))
+        result["octave_band"].extend(octave_binned.axes[0].values.tolist())
+        result["value"].extend(octave_binned.values.tolist())
+
+    result = pd.DataFrame(result)
+    bin_labels = pd.unique(result["octave_band"])
+    band_to_x = {band: i for i, band in enumerate(bin_labels)}
+    result["x_pos"] = result["octave_band"].map(band_to_x)
+
+    fig, ax = plt.subplots(figsize=(8, 4))
+    for name, grp in result.groupby("cochlea"):
+        ax.scatter(grp["x_pos"], grp["value"], label=name, s=60, alpha=0.8)
+
+    ax.set_xticks(range(len(bin_labels)))
+    ax.set_xticklabels(bin_labels)
+    ax.set_xlabel("Octave band (kHz)")
+
+    ax.set_ylabel("Average Ribbon Synapse Count per IHC")
+    ax.set_title("Ribbon synapse count per octave band")
+    ax.legend(title="Cochlea")
+
+    plt.savefig(save_path, bbox_inches="tight", pad_inches=0.1, dpi=png_dpi)
+    if plot:
+        plt.show()
+    else:
+        plt.close()
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Generate plots for Fig 3 of the cochlea paper.")
+    parser.add_argument("--figure_dir", "-f", type=str, help="Output directory for plots.", default="./panels/fig3")
+    parser.add_argument("--plot", action="store_true")
+    args = parser.parse_args()
+
+    os.makedirs(args.figure_dir, exist_ok=True)
+
+    # Panel A: Tonotopic mapping of IHCs (rendering in napari)
+    # fig_03a(save_path=os.path.join(args.figure_dir, "fig_03a.png"))
+
+    # Panel B: Tonotopic mapping of SGNs (rendering in napari)
+    # fig_03b(save_path=os.path.join(args.figure_dir, "fig_03b.png"))
+
+    # Panel C: Spatial distribution of synapses across the cochlea.
+    # We have two options: running sum over the runlength or per octave band
+    fig_03c_rl(save_path=os.path.join(args.figure_dir, "fig_03c_runlength.png"), plot=args.plot)
+    fig_03c_octave(save_path=os.path.join(args.figure_dir, "fig_03c_octave.png"), plot=args.plot)
+
+    # TODO: Panel D: Spatial distribution of SGN sub-types.
+
+
+if __name__ == "__main__":
+    main()