More measurement updates

Author: constantinpape

scripts/export_lower_resolution.py

@@ -10,19 +10,19 @@
 from skimage.segmentation import relabel_sequential
 
 
-def filter_component(fs, segmentation, cochlea, seg_name):
+def filter_component(fs, segmentation, cochlea, seg_name, components):
     # First, we download the MoBIE table for this segmentation.
     internal_path = os.path.join(BUCKET_NAME, cochlea, "tables",  seg_name, "default.tsv")
     with fs.open(internal_path, "r") as f:
         table = pd.read_csv(f, sep="\t")
 
     # Then we get the ids for the components and us them to filter the segmentation.
-    component_mask = np.isin(table.component_labels.values, [1])
+    component_mask = np.isin(table.component_labels.values, components)
     keep_label_ids = table.label_id.values[component_mask].astype("int64")
     filter_mask = ~np.isin(segmentation, keep_label_ids)
     segmentation[filter_mask] = 0
 
-    segmentation, _, _ = relabel_sequential(segmentation)
+    # segmentation, _, _ = relabel_sequential(segmentation)
     return segmentation
 
 
@@ -42,8 +42,8 @@ def export_lower_resolution(args):
         with zarr.open(s3_store, mode="r") as f:
             data = f[input_key][:]
         print(data.shape)
-        if args.filter_by_component:
-            data = filter_component(fs, data, args.cochlea, channel)
+        if args.filter_by_components is not None:
+            data = filter_component(fs, data, args.cochlea, channel, args.filter_by_components)
         tifffile.imwrite(out_path, data, bigtiff=True, compression="zlib")
 
 
@@ -53,7 +53,7 @@ def main():
     parser.add_argument("--scale", "-s", type=int, required=True)
     parser.add_argument("--output_folder", "-o", required=True)
     parser.add_argument("--channels", nargs="+", default=["PV", "VGlut3", "CTBP2"])
-    parser.add_argument("--filter_by_component", action="store_true")
+    parser.add_argument("--filter_by_components", nargs="+", type=int, default=None)
     args = parser.parse_args()
 
     export_lower_resolution(args)

scripts/measurements/evaluate_otof_therapy.py

@@ -49,8 +49,31 @@ def check_project(save=False):
         print()
 
 
+def plot_distribution():
+    import seaborn as sns
+    import matplotlib.pyplot as plt
+
+    table1 = "./results/otof-measurements/M_AMD_OTOF1_L.csv"
+    table2 = "./results/otof-measurements/M_AMD_OTOF2_L.csv"
+
+    table1 = pd.read_csv(table1)
+    table2 = pd.read_csv(table2)
+
+    fig, axes = plt.subplots(1, 2, figsize=(10, 4))
+    sns.histplot(data=table1, x="mean", bins=32, ax=axes[0])
+    axes[0].set_title("Dual AAV")
+    sns.histplot(data=table2, x="mean", bins=32, ax=axes[1])
+    axes[1].set_title("Overloaded AAV")
+
+    fig.suptitle("OTOF Gene Therapy - Mean AlphaTag Intensity of IHCs")
+    plt.tight_layout()
+
+    plt.show()
+
+
 def main():
-    check_project(save=True)
+    # check_project(save=True)
+    plot_distribution()
 
 
 if __name__ == "__main__":

scripts/measurements/evaluate_sgn_therapy.py

@@ -32,8 +32,8 @@ def check_project(save=False):
             continue
 
         # Get the ihc table folder.
-        ihc = sources[sgn_name]["segmentation"]
-        table_folder = os.path.join(BUCKET_NAME, cochlea, ihc["tableData"]["tsv"]["relativePath"])
+        sgn = sources[sgn_name]["segmentation"]
+        table_folder = os.path.join(BUCKET_NAME, cochlea, sgn["tableData"]["tsv"]["relativePath"])
 
         # For debugging.
         x = s3.ls(table_folder)
@@ -47,6 +47,7 @@ def check_project(save=False):
             os.path.join(table_folder, "GFP-resized_SGN-resized-v2_object-measures.tsv"), mode="rb"
         )
         measurement_table = pd.read_csv(measurement_table, sep="\t")
+        breakpoint()
 
         if save:
             os.makedirs(OUTPUT_FOLDER, exist_ok=True)
@@ -58,8 +59,40 @@ def check_project(save=False):
         print()
 
 
+def plot_distribution():
+    import seaborn as sns
+    import matplotlib.pyplot as plt
+
+    table1 = "./results/sgn-measurements/M_LR_000145_L.csv"
+    table2 = "./results/sgn-measurements/M_LR_000151_R.csv"
+    table3 = "./results/sgn-measurements/M_LR_000155_L.csv"
+
+    table1 = pd.read_csv(table1)
+    table2 = pd.read_csv(table2)
+    table3 = pd.read_csv(table3)
+
+    print(len(table1))
+    print(len(table3))
+
+    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
+
+    sns.histplot(data=table1, x="mean", bins=32, ax=axes[0])
+    axes[0].set_title("M145_L")
+    # Something is wrong here, the values are normalized.
+    # sns.histplot(data=table2, x="mean", bins=32, ax=axes[1])
+    # axes[1].set_title("M151_R")
+    sns.histplot(data=table3, x="mean", bins=32, ax=axes[1])
+    axes[1].set_title("M155_L")
+
+    fig.suptitle("SGN Gene Therapy - Mean GFP Intensity of SGNs")
+    plt.tight_layout()
+
+    plt.show()
+
+
 def main():
-    check_project(save=True)
+    # check_project(save=True)
+    plot_distribution()
 
 
 if __name__ == "__main__":

scripts/measurements/measure_sgns.py

@@ -21,10 +21,22 @@ def open_tsv(fs, path):
 
 
 def main():
+    import argparse
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--cochleae", "-c", nargs="+")
+    args = parser.parse_args()
+
     fs = create_s3_target()
     project_info = open_json(fs, "project.json")
-    for dataset in project_info["datasets"]:
-        if dataset == "fens":
+
+    if args.cochleae is None:
+        cochleae = project_info["datasets"]
+    else:
+        cochleae = args.cochleae
+
+    for dataset in cochleae:
+        if dataset not in project_info["datasets"]:
+            print("Could not find cochleae", dataset)
             continue
         print(dataset)
         dataset_info = open_json(fs, os.path.join(dataset, "dataset.json"))
@@ -36,12 +48,22 @@ def main():
             source_info = source_info["segmentation"]
             table_path = source_info["tableData"]["tsv"]["relativePath"]
             table = open_tsv(fs, os.path.join(dataset, table_path, "default.tsv"))
-            component_labels = table.component_labels.values
-            remaining_sgns = component_labels[component_labels != 0]
-            print(source)
-            print("Number of SGNs (all components)   :", len(remaining_sgns))
-            _, n_per_component = np.unique(remaining_sgns, return_counts=True)
-            print("Number of SGNs (largest component):", max(n_per_component))
+
+            if hasattr(table, "component_labels"):
+                component_labels = table.component_labels.values
+                remaining_sgns = component_labels[component_labels != 0]
+                print(source)
+                print(
+                    "Number of SGNs (all components)   :", len(remaining_sgns), "/", len(table),
+                    "(total number of segmented objects)"
+                )
+                component_ids, n_per_component = np.unique(
+                    remaining_sgns, return_counts=True
+                )
+                print("Number of SGNs (largest component):", max(n_per_component))
+            else:
+                print(source)
+                print("Number of SGNs (no postprocessing):", len(table))
 
 
 if __name__ == "__main__":