More plot update

Author: constantinpape

scripts/export_frequency_mapping.py

@@ -36,8 +36,8 @@ def export_frequency_mapping(cochlea, scale, output_folder, source_name, colorma
     seg_path = os.path.join(cochlea, source["imageData"]["ome.zarr"]["relativePath"])
     s3_store, _ = get_s3_path(seg_path, bucket_name=BUCKET_NAME, service_endpoint=SERVICE_ENDPOINT)
     input_key = f"s{scale}"
-    with zarr.open(s3_store, mode="r") as f:
-        seg = f[input_key][:]
+    f = zarr.open(s3_store, mode="r")
+    seg = f[input_key][:]
 
     mapping = {int(seg_id): freq for seg_id, freq in zip(table.label_id, frequencies)}
     lut = np.zeros(max_id + 1, dtype="float32")

scripts/figures/plot_fig3.py

@@ -11,7 +11,14 @@
 
 from util import sliding_runlength_sum, frequency_mapping, SYNAPSE_DIR_ROOT
 
-INPUT_ROOT = "/home/pape/Work/my_projects/flamingo-tools/scripts/M_LR_000227_R/scale3/frequency_mapping"
+INPUT_ROOT = "/home/pape/Work/my_projects/flamingo-tools/scripts/M_LR_000227_R/scale3"
+
+TYPE_TO_CHANNEL = {
+    "Type-Ia": "CR",
+    "Type-Ib": "Calb1",
+    "Type-Ic": "Lypd1",
+    "Type-II": "Prph",
+}
 
 png_dpi = 300
 
@@ -42,35 +49,27 @@ def _plot_colormap(vol, title, plot, save_path):
     plt.tight_layout()
     if plot:
         plt.show()
+
     plt.savefig(save_path)
     plt.close()
 
 
 def fig_03a(save_path, plot, plot_napari):
-    path = os.path.join(INPUT_ROOT, "frequencies_IHC_v4c.tif")
-    vol = imageio.imread(path)
-    _plot_colormap(vol, title="Tonotopic Mapping: IHCs", plot=plot, save_path=save_path)
+    path_ihc = os.path.join(INPUT_ROOT, "frequencies_IHC_v4c.tif")
+    path_sgn = os.path.join(INPUT_ROOT, "frequencies_SGN_v2.tif")
+    sgn = imageio.imread(path_sgn)
+    ihc = imageio.imread(path_ihc)
+    _plot_colormap(sgn, title="Tonotopic Mapping", plot=plot, save_path=save_path)
 
     # Show the image in napari for rendering.
     if plot_napari:
         import napari
         v = napari.Viewer()
-        v.add_image(vol, colormap="viridis")
+        v.add_image(ihc, colormap="viridis")
+        v.add_image(sgn, colormap="viridis")
         napari.run()
 
 
-def fig_03b(save_path, plot, plot_napari):
-    path = os.path.join(INPUT_ROOT, "frequencies_SGN_v2.tif")
-    vol = imageio.imread(path)
-    _plot_colormap(vol, title="Tonotopic Mapping: SGNs", plot=plot, save_path=save_path)
-
-    # Show the image in napari for rendering.
-    if plot_napari:
-        import napari
-        v = napari.Viewer()
-        v.add_image(vol, colormap="viridis")
-
-
 def fig_03c_rl(save_path, plot=False):
     tables = glob("./ihc_counts/ihc_count_M_LR*.tsv")
     fig, ax = plt.subplots(figsize=(8, 4))
@@ -102,9 +101,9 @@ def fig_03c_rl(save_path, plot=False):
 
 
 def fig_03c_octave(save_path, plot=False):
-    # TODO update this table
     ihc_version = "ihc_counts_v4c"
     tables = glob(os.path.join(SYNAPSE_DIR_ROOT, ihc_version, "ihc_count_M_LR*.tsv"))
+    assert len(tables) == 4
 
     result = {"cochlea": [], "octave_band": [], "value": []}
     for tab_path in tables:
@@ -114,7 +113,6 @@ def fig_03c_octave(save_path, plot=False):
         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))
@@ -164,7 +162,11 @@ def fig_03d(save_path, plot, print_stats=True):
                 vals = table[col].values
                 subtype = col[3:]
                 n_subtype = vals.sum()
-                print(subtype, ":", n_subtype, "/", n_sgns, f"({np.round(float(n_subtype) / n_sgns * 100, 2)} %)")
+                channel = TYPE_TO_CHANNEL[subtype]
+                print(
+                    f"{subtype} ({channel}):", n_subtype, "/", n_sgns,
+                    f"({np.round(float(n_subtype) / n_sgns * 100, 2)} %)"
+                )
 
             coexpr = np.logical_and(subtype_table.iloc[:, 0].values, subtype_table.iloc[:, 1].values)
             print("Co-expression:", coexpr.sum())
@@ -178,16 +180,13 @@ def main():
 
     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"), plot=args.plot, plot_napari=False)
-
-    # Panel B: Tonotopic mapping of SGNs (rendering in napari)
-    # fig_03b(save_path=os.path.join(args.figure_dir, "fig_03b.png"), plot=args.plot, plot_napari=False)
+    # Panel A: Tonotopic mapping of SGNs and IHCs (rendering in napari + heatmap)
+    # fig_03a(save_path=os.path.join(args.figure_dir, "fig_03a_cmap.png"), plot=args.plot, plot_napari=True)
 
     # 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)
+    fig_03c_octave(save_path=os.path.join(args.figure_dir, "fig_03c_octave.png"), plot=args.plot)
 
     # Panel D: Spatial distribution of SGN sub-types.
     fig_03d(save_path=os.path.join(args.figure_dir, "fig_03d.png"), plot=args.plot)

scripts/figures/util.py

@@ -2,8 +2,8 @@
 import numpy as np
 
 # Directory with synapse measurement tables
-SYNAPSE_DIR_ROOT = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/predictions/synapses"
-# SYNAPSE_DIR_ROOT = "./synapses"
+# SYNAPSE_DIR_ROOT = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/predictions/synapses"
+SYNAPSE_DIR_ROOT = "./synapses"
 
 
 # Define the animal specific octave bands.

scripts/measurements/sgn_subtypes.py

@@ -10,7 +10,8 @@
     "M_LR_000099_L": ["PV", "Calb1", "Lypd1"],
     "M_LR_000214_L": ["PV", "CR", "Calb1"],
     "M_AMD_N62_L": ["PV", "CR", "Calb1"],
-    "M_AMD_Runx1_L": ["PV", "Lypd1", "Calb1"],
+    # Mutant / some stuff is weird.
+    # "M_AMD_Runx1_L": ["PV", "Lypd1", "Calb1"],
     # This one still has to be stitched:
     # "M_LR_000184_R": {"PV", "Prph"},
     # We don't have PV here, so we exclude these two for now.
@@ -32,6 +33,15 @@
     "Prph": "Type-II",
 }
 
+# TODO
+THRESHOLDS = {
+    "M_LR_000214_L": {
+    },
+    "M_AMD_N62_L": {
+        "Calb1": 380,
+    },
+}
+
 
 def check_processing_status():
     s3 = create_s3_target()
@@ -83,10 +93,6 @@ def analyze_subtypes_intensity_based():
         # Remove the PV channel, which we don't need for analysis.
         channels = channels[1:]
 
-        # FIXME
-        if cochlea != "M_LR_000099_L":
-            continue
-
         content = s3.open(f"{BUCKET_NAME}/{cochlea}/dataset.json", mode="r", encoding="utf-8")
         info = json.loads(content.read())
         sources = info["sources"]
@@ -110,7 +116,11 @@ def analyze_subtypes_intensity_based():
         for channel in channels:
             # Load the intensity table.
             intensity_path = os.path.join(table_folder, f"{channel}_PV-SGN-v2_object-measures.tsv")
-            table_content = s3.open(intensity_path, mode="rb")
+            try:
+                table_content = s3.open(intensity_path, mode="rb")
+            except FileNotFoundError:
+                print(intensity_path, "is missing")
+                continue
             intensities = pd.read_csv(table_content, sep="\t")
             intensities = intensities[intensities.label_id.isin(valid_sgns)]
             assert len(table) == len(intensities)
@@ -120,7 +130,10 @@ def analyze_subtypes_intensity_based():
             medians = intensities["median"].values
 
             # TODO: we need to determine the threshold in a better way / validate it in MoBIE.
-            intensity_threshold = float(threshold_otsu(medians))
+            intensity_threshold = THRESHOLDS.get(cochlea, {}).get(channel, None)
+            if intensity_threshold is None:
+                print("Could not find a threshold for", cochlea, channel, "falling back to OTSU")
+                intensity_threshold = float(threshold_otsu(medians))
             threshold_dict[cochlea][channel] = intensity_threshold
 
             subtype = CHANNEL_TO_TYPE[channel]