Update inner ear analysis scripts

Author: constantinpape

scripts/cooper/full_reconstruction/visualize_results.py

@@ -6,11 +6,14 @@
 import numpy as np
 import pandas as pd
 
+from skimage.filters import gaussian
+
 ROOT = "./04_full_reconstruction"
 TABLE = "/home/pape/Desktop/sfb1286/mboc_synapse/draft_figures/full_reconstruction.xlsx"
 
 # Skip datasets for which all figures were already done.
-SKIP_DS = ["20241019_Tomo-eval_MF_Synapse"]
+SKIP_DS = ["20241019_Tomo-eval_MF_Synapse", "20241019_Tomo-eval_PS_Synapse"]
+# SKIP_DS = []
 
 
 def _get_name_and_row(path, table):
@@ -46,20 +49,23 @@ def visualize_result(path, table):
     if ds_name in SKIP_DS:
         return
 
-    # if row["Use for vis"].values[0] == "yes":
-    if row["Use for vis"].values[0] in ("yes", "no"):
+    if row["Use for Vis"].values[0] == "no":
         return
     compartment_ids = _get_compartment_ids(row)
 
     # access = np.s_[:]
-    access = np.s_[::2, ::2, ::2]
+    access = np.s_[::3, ::3, ::3]
 
     with h5py.File(path, "r") as f:
         raw = f["raw"][access]
         vesicles = f["labels/vesicles"][access]
         active_zone = f["labels/active_zone"][access]
         mitos = f["labels/mitochondria"][access]
         compartments = f["labels/compartments"][access]
+    print("Loading done")
+
+    raw = gaussian(raw)
+    print("Gaussian done")
 
     if any(comp_ids is not None for comp_ids in compartment_ids):
         mask = np.zeros(raw.shape, dtype="bool")
@@ -78,12 +84,14 @@ def visualize_result(path, table):
         mitos[~mask] = 0
         compartments = compartments_new
 
+    vesicle_ids = np.unique(vesicles)[1:]
+
     v = napari.Viewer()
     v.add_image(raw)
     v.add_labels(mitos)
-    v.add_labels(vesicles)
-    v.add_labels(compartments)
-    v.add_labels(active_zone)
+    v.add_labels(vesicles, colormap={ves_id: "orange" for ves_id in vesicle_ids})
+    v.add_labels(compartments, colormap={1: "red", 2: "green", 3: "orange"})
+    v.add_labels(active_zone, colormap={1: "blue"})
     v.title = f"{ds_name}/{name}"
     napari.run()
 
@@ -115,6 +123,7 @@ def main():
     paths = sorted(glob(os.path.join(ROOT, "**/*.h5"), recursive=True))
     table = pd.read_excel(TABLE)
     for path in paths:
+        print(path)
         visualize_result(path, table)
         # visualize_only_compartment(path, table)
 

scripts/inner_ear/analysis/analyze_distances.py

@@ -35,7 +35,7 @@ def _plot_all(distances):
 # TODO rename the method names.
 # We only care about the following distances:
 # - MP-V -> PD, AZ (Boundary)
-# - Docked-V -> PD
+# - Docked-V -> PD, AZ
 # - RA-V -> Ribbon
 def _plot_selected(distances, save_path=None):
     fig, axes = plt.subplots(2, 2)
@@ -46,7 +46,7 @@ def _plot_selected(distances, save_path=None):
 
     def _plot(pool_name, distance_col, structure_name, ax):
 
-        this_distances = distances[distances["pool"] == pool_name][["approach", distance_col]]
+        this_distances = distances[distances["pool"] == pool_name][["tomogram", "approach", distance_col]]
 
         ax.set_title(f"{pool_name} to {structure_name}")
         sns.histplot(
@@ -56,15 +56,27 @@ def _plot(pool_name, distance_col, structure_name, ax):
 
         if save_path is not None:
             approaches = pd.unique(this_distances["approach"])
-            dist_values = [
-                this_distances[this_distances["approach"] == approach][distance_col].values.tolist()
-                for approach in approaches
-            ]
-            max_len = max([len(vals) for vals in dist_values])
-            save_distances = {
-                approach: dists + [np.nan] * (max_len - len(dists))
-                for approach, dists in zip(approaches, dist_values)
-            }
+            tomo_names = pd.unique(this_distances["tomogram"])
+
+            tomograms = []
+            distance_values = {approach: [] for approach in approaches}
+
+            for tomo in tomo_names:
+                tomo_dists = this_distances[this_distances["tomogram"] == tomo]
+                max_vesicles = 0
+                for approach in approaches:
+                    n_vesicles = len(tomo_dists[tomo_dists["approach"] == approach].values)
+                    if n_vesicles > max_vesicles:
+                        max_vesicles = n_vesicles
+
+                for approach in approaches:
+                    app_dists = tomo_dists[tomo_dists["approach"] == approach][distance_col].values.tolist()
+                    app_dists = app_dists + [np.nan] * (max_vesicles - len(app_dists))
+                    distance_values[approach].extend(app_dists)
+                tomograms.extend([tomo] * max_vesicles)
+
+            save_distances = {"tomograms": tomograms}
+            save_distances.update(distance_values)
             save_distances = pd.DataFrame(save_distances)
 
             sheet_name = f"{pool_name}_{structure_name}"
@@ -74,9 +86,11 @@ def _plot(pool_name, distance_col, structure_name, ax):
             else:
                 save_distances.to_excel(save_path, index=False, sheet_name=sheet_name)
 
+    # NOTE: we over-ride a plot here, should not do this in the actual version
     _plot("MP-V", "pd_distance [nm]", "PD", axes[0, 0])
     _plot("MP-V", "boundary_distance [nm]", "AZ Membrane", axes[0, 1])
     _plot("Docked-V", "pd_distance [nm]", "PD", axes[1, 0])
+    _plot("Docked-V", "boundary_distance [nm]", "AZ Membrane", axes[1, 0])
     _plot("RA-V", "ribbon_distance [nm]", "Ribbon", axes[1, 1])
 
     fig.tight_layout()
@@ -87,17 +101,17 @@ def for_tomos_with_annotation(plot_all=True):
     manual_assignments, semi_automatic_assignments, automatic_assignments = get_measurements_with_annotation()
 
     manual_distances = manual_assignments[
-        ["pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
+        ["tomogram", "pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
     ]
     manual_distances["approach"] = ["manual"] * len(manual_distances)
 
     semi_automatic_distances = semi_automatic_assignments[
-        ["pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
+        ["tomogram", "pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
     ]
     semi_automatic_distances["approach"] = ["semi_automatic"] * len(semi_automatic_distances)
 
     automatic_distances = automatic_assignments[
-        ["pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
+        ["tomogram", "pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
     ]
     automatic_distances["approach"] = ["automatic"] * len(automatic_distances)
 
@@ -113,12 +127,12 @@ def for_all_tomos(plot_all=True):
     semi_automatic_assignments, automatic_assignments = get_all_measurements()
 
     semi_automatic_distances = semi_automatic_assignments[
-        ["pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
+        ["tomogram", "pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
     ]
     semi_automatic_distances["approach"] = ["semi_automatic"] * len(semi_automatic_distances)
 
     automatic_distances = automatic_assignments[
-        ["pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
+        ["tomogram", "pool", "ribbon_distance [nm]", "pd_distance [nm]", "boundary_distance [nm]"]
     ]
     automatic_distances["approach"] = ["automatic"] * len(automatic_distances)
 

scripts/inner_ear/analysis/analyze_vesicle_pools.py

@@ -94,8 +94,7 @@ def for_all_tomos():
         errors.to_excel(writer, sheet_name="StandardDeviation", index=False)
 
 
-# TODO: export the vesicle diameters
-# TODO: export the ribbon and pd stats
+# TODO: export the ribbon and pd stats (first need to discuss this with Fid)
 def main():
     # for_tomos_with_annotation()
     for_all_tomos()

scripts/summarize_data.py

@@ -0,0 +1,27 @@
+import numpy as np
+import pandas as pd
+
+
+az_train = pd.read_excel("data_summary/active_zone_training_data.xlsx")
+compartment_train = pd.read_excel("data_summary/compartment_training_data.xlsx")
+vesicle_train = pd.read_excel("data_summary/vesicle_training_data.xlsx")
+vesicle_da = pd.read_excel("data_summary/vesicle_domain_adaptation_data.xlsx", sheet_name="cryo")
+
+
+def training_resolutions():
+    res_az = np.round(az_train["resolution"].mean(), 2)
+    res_compartment = np.round(compartment_train["resolution"].mean(), 2)
+    res_cryo = np.round(vesicle_da["resolution"].mean(), 2)
+    res_vesicles = np.round(vesicle_train["resolution"].mean(), 2)
+
+    print("Training resolutions for models:")
+    print("active_zone:", res_az)
+    print("compartments:", res_compartment)
+    # TODO
+    print("mitochondria:", 1.0)
+    print("vesicles_2d:", res_vesicles)
+    print("vesicles_3d:", res_vesicles)
+    print("vesicles_cryo:", res_cryo)
+
+
+training_resolutions()