Update figures

Author: schilling40

scripts/figures/plot_fig2.py

@@ -1,6 +1,5 @@
 import argparse
 import os
-from glob import glob
 
 import numpy as np
 import pandas as pd
@@ -100,7 +99,9 @@ def fig_02c(save_path, plot=False, all_versions=False):
 
 # Load the synapse counts for all IHCs from the relevant tables.
 def _load_ribbon_synapse_counts():
-    tables = glob("ihc_counts/*M_LR*.tsv")
+    ihc_version = "ihc_counts_v4c"
+    synapse_dir = f"/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/predictions/synapses/{ihc_version}"
+    tables = [entry.path for entry in os.scandir(synapse_dir) if "ihc_count_M_LR" in entry.name]
     syn_counts = []
     for tab in tables:
         x = pd.read_csv(tab, sep="\t")
@@ -141,9 +142,9 @@ def fig_02d_01(save_path, plot=False, all_versions=False, plot_average_ribbon_sy
         # Labels and formatting
         ax[0].set_xticklabels(["SGN"], fontsize=main_label_size)
 
-        ylim0 = 9500
+        ylim0 = 8500
         ylim1 = 12500
-        y_ticks = [i for i in range(10000, 12000 + 1, 1000)]
+        y_ticks = [i for i in range(9000, 12000 + 1, 1000)]
 
         ax[0].set_ylabel("Count per cochlea", fontsize=main_label_size)
         ax[0].set_yticks(y_ticks)
@@ -157,12 +158,12 @@ def fig_02d_01(save_path, plot=False, all_versions=False, plot_average_ribbon_sy
         ax[0].set_xlim(xmin, xmax)
         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", color="C0", fontsize=main_tick_size, ha="center")
+        ax[0].text(1.1, (lower_y + upper_y) // 2, "literature", color="C0", fontsize=main_tick_size, ha="left")
         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)]
+        ylim0 = 600
+        ylim1 = 800
+        y_ticks = [i for i in range(600, 800 + 1, 100)]
 
         ax[1].set_xticklabels(["IHC"], fontsize=main_label_size)
         ax[1].set_yticks(y_ticks)
@@ -178,28 +179,28 @@ def fig_02d_01(save_path, plot=False, all_versions=False, plot_average_ribbon_sy
         lower_y, upper_y = literature_reference_values("IHC")
         ax[1].set_xlim(xmin, xmax)
         ax[1].hlines([lower_y, upper_y], xmin, xmax)
-        ax[1].text(1, upper_y + 20, "literature", color="C0", fontsize=main_tick_size, ha="center")
+        ax[1].text(1.1, (lower_y + upper_y) // 2, "literature", color="C0", fontsize=main_tick_size, ha="left")
         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
+            ylim0 = -1
+            ylim1 = 41
             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)
+            ax[2].set_ylim(ylim0, ylim1)
 
             # set range of literature values
             xmin = 0.5
             xmax = 1.5
             lower_y, upper_y = literature_reference_values("synapse")
             ax[2].set_xlim(xmin, xmax)
             ax[2].hlines([lower_y, upper_y], xmin, xmax)
-            ax[2].text(1, upper_y + 2, "literature", color="C0", fontsize=main_tick_size, ha="center")
+            ax[2].text(1.1, (lower_y + upper_y) // 2, "literature", color="C0", fontsize=main_tick_size, ha="left")
             ax[2].fill_between([xmin, xmax], lower_y, upper_y, color="C0", alpha=0.05, interpolate=True)
 
         plt.tight_layout()

scripts/figures/plot_fig4.py

@@ -71,8 +71,9 @@ def get_chreef_data():
         chreef_data[cochlea] = values
 
     with open(cache_path, "wb") as f:
-        chreef_data = pickle.dump(chreef_data, f)
-    return chreef_data
+        pickle.dump(chreef_data, f)
+    with open(cache_path, "rb") as f:
+        return pickle.load(f)
 
 
 def group_lr(names_lr, values):
@@ -143,7 +144,7 @@ def fig_04c(chreef_data, save_path, plot=False, plot_by_side=False):
     plt.xlim(xmin, xmax)
     lower_y, upper_y = literature_reference_values("SGN")
     plt.hlines([lower_y, upper_y], xmin, xmax)
-    plt.text(1, lower_y - 400, "literature", color="C0", fontsize=main_tick_size, ha="center")
+    plt.text(1.5, lower_y - 400, "literature", color="C0", fontsize=main_tick_size, ha="center")
     plt.fill_between([xmin, xmax], lower_y, upper_y, color="C0", alpha=0.05, interpolate=True)
 
     sgn_values = [11153, 11398, 10333, 11820]
@@ -154,12 +155,12 @@ def fig_04c(chreef_data, save_path, plot=False, plot_by_side=False):
     lower_y = sgn_value - 1.96 * sgn_std
 
     plt.hlines([lower_y, upper_y], xmin, xmax, colors=["C1" for _ in range(2)])
-    plt.text(1, upper_y + 100, "healthy cochleae (95% confidence interval)",
+    plt.text(1.5, upper_y + 100, "healthy cochleae (95% confidence interval)",
              color="C1", fontsize=main_tick_size, ha="center")
     plt.fill_between([xmin, xmax], lower_y, upper_y, color="C1", alpha=0.05, interpolate=True)
 
-    plt.savefig(save_path, bbox_inches="tight", pad_inches=0.1, dpi=png_dpi)
     plt.tight_layout()
+    plt.savefig(save_path, bbox_inches="tight", pad_inches=0.1, dpi=png_dpi)
 
     if plot:
         plt.show()
@@ -179,12 +180,10 @@ def fig_04d(chreef_data, save_path, plot=False, plot_by_side=False, intensity=Fa
             intensities = vals["median"].values
             values.append(intensities.mean())
         else:
-            # The marker labels don't make sense yet, they are in
-            # 0: unlabeled
+            # marker labels
+            # 0: unlabeled - no median intensity in object_measures table
             # 1: positive
             # 2: negative
-            # but they should all be either positive or negative.
-            # Or am I missing something?
             marker_labels = vals["marker_labels"].values
             n_pos = (marker_labels == 1).sum()
             n_neg = (marker_labels == 2).sum()
@@ -222,8 +221,8 @@ def fig_04d(chreef_data, save_path, plot=False, plot_by_side=False, intensity=Fa
     if not intensity:
         plt.ylim(0.5, 1.05)
 
-    plt.savefig(save_path, bbox_inches="tight", pad_inches=0.1, dpi=png_dpi)
     plt.tight_layout()
+    plt.savefig(save_path, bbox_inches="tight", pad_inches=0.1, dpi=png_dpi)
 
     if plot:
         plt.show()
@@ -255,30 +254,76 @@ def fig_04e(chreef_data, save_path, plot, intensity=False):
     band_to_x = {band: i for i, band in enumerate(bin_labels)}
     result["x_pos"] = result["octave_band"].map(band_to_x)
 
-    fig, axes = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True)
+    fig, ax = plt.subplots(figsize=(8, 4))
+
+    sub_tick_label_size = 8
+    tick_label_size = 12
+    label_size = 12
+    legend_size = 8
+    if intensity:
+        band_label_offset_y = 0.09
+    else:
+        band_label_offset_y = 0.07
+        ax.set_ylim(0.45, 1.05)
+
+    # Offsets within each octave band
+    offset_map = {"L": -0.15, "R": 0.15}
+    sublabels = {"L": "L", "R": "R"}
+
+    # Assign a color to each cochlea (ignoring side)
+    cochleas = sorted({name_lr[:-1] for name_lr in result["cochlea"].unique()})
+    colors = plt.cm.tab10.colors  # pick a colormap
+    color_map = {cochlea: colors[i % len(colors)] for i, cochlea in enumerate(cochleas)}
+
+    # Track which cochlea names we have already added to the legend
+    legend_added = set()
+
+    all_x_positions = []
+    all_x_labels = []
+
     for name_lr, grp in result.groupby("cochlea"):
         name, side = name_lr[:-1], name_lr[-1]
-        ax, marker = (axes[0], "o") if side == "L" else (axes[1], "x")
-        ax.scatter(grp["x_pos"], grp["value"], label=name, s=60, alpha=0.8, marker=marker)
-
-    for ax in axes:
-        ax.set_xticks(range(len(bin_labels)))
-        ax.set_xticklabels(bin_labels)
-        ax.set_xlabel("Octave band (kHz)")
+        x_positions = grp["x_pos"] + offset_map[side]
+        ax.scatter(
+            x_positions,
+            grp["value"],
+            label=name if name not in legend_added else None,
+            s=60,
+            alpha=0.8,
+            marker="o" if side == "L" else "x",
+            color=color_map[name]
+        )
+        if name not in legend_added:
+            legend_added.add(name)
+
+        # Store for sublabel ticks
+        all_x_positions.extend(x_positions)
+        all_x_labels.extend([sublabels[side]] * len(x_positions))
+
+    # Create combined tick positions & labels
+    main_ticks = range(len(bin_labels))
+    # add a final tick for label '>64k'
+    ax.set_xticks([pos + offset_map["L"] for pos in main_ticks[:-1]] +
+                  [pos + offset_map["R"] for pos in main_ticks[:-1]] +
+                  [pos for pos in main_ticks[-1:]])
+    ax.set_xticklabels(["L"] * len(main_ticks[:-1]) + ["R"] * len(main_ticks[:-1]) + [""], fontsize=sub_tick_label_size)
+
+    # Add main octave band labels above sublabels
+    for i, label in enumerate(bin_labels):
+        ax.text(i, ax.get_ylim()[0] - band_label_offset_y*(ax.get_ylim()[1]-ax.get_ylim()[0]),
+                label, ha='center', va='top', fontsize=tick_label_size, fontweight='bold')
+
+    ax.set_xlabel("Octave band (kHz)", fontsize=label_size)
+    ax.xaxis.set_label_coords(.5, -.16)
 
     if intensity:
-        axes[0].set_ylabel("Marker Intensity")
-        for ax, side in zip(axes, ("Left", "Right")):
-            ax.set_title(f"Intensity per octave band ({side})")
+        ax.set_ylabel("Marker Intensity", fontsize=label_size)
+        ax.set_title("Intensity per octave band (Left/Right)")
     else:
-        axes[0].set_ylabel("Transduction Efficiency")
-        axes[0].set_ylim(0.5, 1.05)
-        for ax, side in zip(axes, ("Left", "Right")):
-            ax.set_title(f"Transduction efficiency per octave band ({side})")
-
-    # FIXME make this uniform across the plots!
-    axes[0].legend(title="Cochlea")
-    axes[1].legend(title="Cochlea")
+        ax.set_ylabel("Transduction Efficiency", fontsize=label_size)
+        ax.set_title("Transduction efficiency per octave band (Left/Right)")
+
+    ax.legend(title="Cochlea", fontsize=legend_size)
     plt.tight_layout()
 
     plt.savefig(save_path, bbox_inches="tight", pad_inches=0.1, dpi=png_dpi)

scripts/figures/plot_fig5.py

@@ -0,0 +1,93 @@
+import argparse
+import os
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+from util import literature_reference_values
+
+png_dpi = 300
+
+
+def _load_ribbon_synapse_counts():
+    ihc_version = "ihc_counts_v4b"
+    synapse_dir = f"/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/predictions/synapses/{ihc_version}"
+    tables = [entry.path for entry in os.scandir(synapse_dir) if "ihc_count_M_AMD" in entry.name]
+    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_05c(save_path, plot=False):
+    """Bar plot showing the distribution of synapse markers per IHC segmentation average over OTOF cochlea.
+    """
+
+    main_label_size = 20
+    main_tick_size = 12
+    htext_size = 10
+
+    ribbon_synapse_counts = _load_ribbon_synapse_counts()
+
+    fig, ax = plt.subplots(figsize=(8, 4))
+
+    ylim0 = -1
+    ylim1 = 24
+    y_ticks = [i for i in range(0, 25, 5)]
+
+    ax.set_ylabel("Ribbon Syn. per IHC", fontsize=main_label_size)
+    ax.set_yticks(y_ticks)
+    ax.set_yticklabels(y_ticks, rotation=0, fontsize=main_tick_size)
+    ax.set_ylim(ylim0, ylim1)
+
+    ax.boxplot(ribbon_synapse_counts)
+    ax.set_xticklabels(["MOTOF1L"], fontsize=main_label_size)
+
+    # set range of literature values
+    xmin = 0.5
+    xmax = 1.5
+    lower_y, upper_y = literature_reference_values("synapse")
+    ax.set_xlim(xmin, xmax)
+    ax.hlines([lower_y, upper_y], xmin, xmax)
+    ax.text(1.25, upper_y + 0.01*ax.get_ylim()[1]-ax.get_ylim()[0], "literature",
+            color="C0", fontsize=htext_size, ha="center")
+    ax.fill_between([xmin, xmax], lower_y, upper_y, color="C0", alpha=0.05, interpolate=True)
+
+    ihc_values = [14.1, 12.7, 13.8, 13.4]  # MLR226L, MLR226R, MLR227L, MLR227R
+
+    ihc_value = np.mean(ihc_values)
+    ihc_std = np.std(ihc_values)
+
+    upper_y = ihc_value + 1.96 * ihc_std
+    lower_y = ihc_value - 1.96 * ihc_std
+
+    plt.hlines([lower_y, upper_y], xmin, xmax, colors=["C1" for _ in range(2)])
+    plt.text(1.25, upper_y + 0.01*ax.get_ylim()[1]-ax.get_ylim()[0], "healthy cochleae (95% confidence interval)",
+             color="C1", fontsize=htext_size, ha="center")
+    plt.fill_between([xmin, xmax], lower_y, upper_y, color="C1", alpha=0.05, interpolate=True)
+
+    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 main():
+    parser = argparse.ArgumentParser(description="Generate plots for Fig 5 of the cochlea paper.")
+    parser.add_argument("--figure_dir", "-f", type=str, help="Output directory for plots.", default="./panels/fig5")
+    parser.add_argument("--plot", action="store_true")
+    args = parser.parse_args()
+
+    os.makedirs(args.figure_dir, exist_ok=True)
+
+    # Panel C: Monitoring of the Syn / IHC loss
+    fig_05c(save_path=os.path.join(args.figure_dir, "fig_05c"), plot=args.plot)
+
+
+if __name__ == "__main__":
+    main()

scripts/figures/util.py

@@ -63,11 +63,11 @@ def sliding_runlength_sum(run_length, values, width):
 # TODO determine these from Aleyna's table!
 def literature_reference_values(structure):
     if structure == "SGN":
-        lower_bound, upper_bound = 10000, 12000
+        lower_bound, upper_bound = 8728, 10325
     elif structure == "IHC":
-        lower_bound, upper_bound = 780, 850
+        lower_bound, upper_bound = 656, 681
     elif structure == "synapse":
-        lower_bound, upper_bound = 10, 25
+        lower_bound, upper_bound = 9.1, 20.7
     else:
         raise ValueError
     return lower_bound, upper_bound