Update otoferlin analysis

Author: constantinpape

scripts/otoferlin/.gitignore

@@ -1 +1,4 @@
 data/
+sync_segmentation.sh
+segmentation/
+results/

scripts/otoferlin/common.py

@@ -36,8 +36,7 @@ def get_folders():
         return INPUT_ROOT, OUTPUT_ROOT
     root_in = "./data/tomograms"
     assert os.path.exists(root_in)
-    root_out = "./data/segmentation"
-    return root_in, root_out
+    return root_in, OUTPUT_ROOT
 
 
 def get_all_tomograms():
@@ -89,6 +88,12 @@ def load_segmentations(seg_path):
     return segmentations
 
 
+def to_condition(mrc_path):
+    fname = os.path.basename(mrc_path)
+    # TODO: Is this correct, or is it the otherway round?
+    return "MUT" if fname.startswith("Otof") else "WT"
+
+
 if __name__ == "__main__":
     tomos = get_all_tomograms()
     print("We have", len(tomos), "tomograms")

scripts/otoferlin/correct_vesicle_pools.py

@@ -18,8 +18,10 @@ def _create_pool_layer(seg, assignment_path):
 
     pool_colors = get_colormaps()["pools"]
     colormap = {}
-    for pool_id, pool_name in enumerate(pool_names, 1):
-        pool_vesicle_ids = assignments[assignments.pool == pool_name].vesicle_id
+    for pool_id, pool_name in enumerate(pool_names):
+        if not isinstance(pool_name, str) and np.isnan(pool_name):
+            continue
+        pool_vesicle_ids = assignments[assignments.pool == pool_name].vesicle_id.values
         pool_mask = np.isin(seg, pool_vesicle_ids)
         pools[pool_mask] = pool_id
         colormap[pool_id] = pool_colors[pool_name]
@@ -37,7 +39,7 @@ def _update_assignments(vesicles, pool_correction, assignment_path):
         correction_val = prop.max_intensity
         if correction_val == 0:
             continue
-        new_assignments[new_assignments.vesicle_id == prop.label] = val_to_pool[correction_val]
+        new_assignments[new_assignments.vesicle_id == prop.label].pool = val_to_pool[correction_val]
 
     new_assignments.to_csv(assignment_path, index=False)
 
@@ -48,6 +50,9 @@ def correct_vesicle_pools(mrc_path):
 
     output_folder = os.path.split(seg_path)[0]
     assignment_path = os.path.join(output_folder, "vesicle_pools.csv")
+    if not os.path.exists(assignment_path):
+        print("Skip", seg_path, "due to missing assignments")
+        return
 
     data, _ = read_mrc(mrc_path)
     segmentations = load_segmentations(seg_path)
@@ -62,6 +67,7 @@ def correct_vesicle_pools(mrc_path):
     vesicle_pools, pool_colors = _create_pool_layer(vesicles, assignment_path)
 
     pool_correction_path = os.path.join(output_folder, "correction", "pool_correction.tif")
+    os.makedirs(os.path.join(output_folder, "correction"), exist_ok=True)
     if os.path.exists(pool_correction_path):
         pool_correction = imageio.imread(pool_correction_path)
     else:
@@ -81,9 +87,10 @@ def update_pools(viewer: napari.Viewer):
         vesicles = viewer.layers["vesicles"].data
         pool_correction = viewer.layers["pool_correction"].data
         _update_assignments(vesicles, pool_correction, assignment_path)
-        imageio.imwrite(pool_correction_path, pool_correction, compression="zlib")
+        # imageio.imwrite(pool_correction_path, pool_correction, compression="zlib")
         pool_data, pool_colors = _create_pool_layer(vesicles, assignment_path)
         viewer.layers["vesicle_pools"].data = pool_data
+        viewer.layers["vesicle_pools"].colormap = pool_colors
 
     v.window.add_dock_widget(update_pools)
 

scripts/otoferlin/export_results.py

@@ -0,0 +1,114 @@
+import os
+from datetime import datetime
+
+import numpy as np
+import pandas as pd
+from common import get_all_tomograms, get_seg_path, to_condition
+
+from synapse_net.distance_measurements import load_distances
+
+
+def get_output_folder():
+    output_root = "./results"
+    date = datetime.now().strftime("%Y%m%d")
+
+    version = 1
+    output_folder = os.path.join(output_root, f"{date}_{version}")
+    while os.path.exists(output_folder):
+        version += 1
+        output_folder = os.path.join(output_root, f"{date}_{version}")
+
+    os.makedirs(output_folder)
+    return output_folder
+
+
+def _export_results(tomograms, result_path, result_extraction):
+    results = {}
+    for tomo in tomograms:
+        condition = to_condition(tomo)
+        res = result_extraction(tomo)
+        if condition in results:
+            results[condition].append(res)
+        else:
+            results[condition] = [res]
+
+    for condition, res in results.items():
+        res = pd.concat(res)
+        if os.path.exists(result_path):
+            with pd.ExcelWriter(result_path, engine="openpyxl", mode="a") as writer:
+                res.to_excel(writer, sheet_name=condition, index=False)
+        else:
+            res.to_excel(result_path, sheet_name=condition, index=False)
+
+
+def export_vesicle_pools(tomograms, result_path):
+
+    def result_extraction(tomo):
+        folder = os.path.split(get_seg_path(tomo))[0]
+        measure_path = os.path.join(folder, "vesicle_pools.csv")
+        measures = pd.read_csv(measure_path).dropna()
+        pool_names, counts = np.unique(measures.pool.values, return_counts=True)
+        res = {"tomogram": [os.path.basename(tomo)]}
+        res.update({k: v for k, v in zip(pool_names, counts)})
+        res = pd.DataFrame(res)
+        return res
+
+    _export_results(tomograms, result_path, result_extraction)
+
+
+def export_distances(tomograms, result_path):
+    def result_extraction(tomo):
+        folder = os.path.split(get_seg_path(tomo))[0]
+        measure_path = os.path.join(folder, "vesicle_pools.csv")
+        measures = pd.read_csv(measure_path).dropna()
+
+        measures = measures[measures.pool.isin(["MP-V", "Docked-V"])][["vesicle_id", "pool"]]
+
+        # Load the distances to PD.
+        pd_distances, _, _, seg_ids = load_distances(os.path.join(folder, "distances", "PD.npz"))
+        pd_distances = {sid: dist for sid, dist in zip(seg_ids, pd_distances)}
+        measures["distance-to-pd"] = [pd_distances[vid] for vid in measures.vesicle_id.values]
+
+        # Load the distances to membrane.
+        mem_distances, _, _, seg_ids = load_distances(os.path.join(folder, "distances", "membrane.npz"))
+        mem_distances = {sid: dist for sid, dist in zip(seg_ids, mem_distances)}
+        measures["distance-to-membrane"] = [mem_distances[vid] for vid in measures.vesicle_id.values]
+
+        measures = measures.drop(columns=["vesicle_id"])
+        measures.insert(0, "tomogram", len(measures) * [os.path.basename(tomo)])
+
+        return measures
+
+    _export_results(tomograms, result_path, result_extraction)
+
+
+def export_diameter(tomograms, result_path):
+    def result_extraction(tomo):
+        folder = os.path.split(get_seg_path(tomo))[0]
+        measure_path = os.path.join(folder, "vesicle_pools.csv")
+        measures = pd.read_csv(measure_path).dropna()
+
+        measures = measures[measures.pool.isin(["MP-V", "Docked-V"])][["pool", "diameter"]]
+        measures.insert(0, "tomogram", len(measures) * [os.path.basename(tomo)])
+
+        return measures
+
+    _export_results(tomograms, result_path, result_extraction)
+
+
+def main():
+    tomograms = get_all_tomograms()
+    result_folder = get_output_folder()
+
+    result_path = os.path.join(result_folder, "vesicle_pools.xlsx")
+    export_vesicle_pools(tomograms, result_path)
+
+    result_path = os.path.join(result_folder, "distances.xlsx")
+    export_distances(tomograms, result_path)
+
+    result_path = os.path.join(result_folder, "diameter.xlsx")
+    export_diameter(tomograms, result_path)
+
+
+if __name__ == "__main__":
+    main()