Grouped subtype plots; Manual thresholds

Author: schilling40

scripts/assign_subtypes.py

@@ -12,6 +12,8 @@
     "M_LR_000184_L": {"seg_data": "SGN_v2b", "subtype": ["Prph"]},
     "M_LR_000184_R": {"seg_data": "SGN_v2b", "subtype": ["Prph"]},
     "M_LR_000260_L": {"seg_data": "SGN_v2", "subtype": ["Prph", "Tuj1"]},
+    "M_LR_N110_L": {"seg_data": "SGN_v2", "subtype": ["Calb1", "Ntng1"]},
+    "M_LR_N110_R": {"seg_data": "SGN_v2", "subtype": ["Calb1", "Ntng1"]},
     "M_LR_N152_L": {"seg_data": "SGN_v2", "subtype": ["CR", "Ntng1"]},
     "M_AMD_N180_L": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"]},
     "M_AMD_N180_R": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"]},
@@ -46,6 +48,12 @@
     "CR+/Ntng1-": "Type Ia",
     "CR-/Ntng1+": "Type Ic",
     "CR-/Ntng1-": "inconclusive",
+
+    # Combinations of Calb1 and Ntng1
+    "Calb1+/Ntng1+": "Type Ib",
+    "Calb1+/Ntng1-": "inconclusive",
+    "Calb1-/Ntng1+": "Type Ic",
+    "Calb1-/Ntng1-": "inconclusive",
 }
 
 
@@ -110,14 +118,14 @@ def filter_subtypes(cochlea, seg_name, subtype, stains=None):
     if len(stain_dict) == 0:
         raise ValueError("The dictionary containing stain information must have at least one entry. Check parameters.")
 
-    subset = table_seg.copy()
-
+    label_ids_subtype = []
     for dic in stain_dict:
+        subset = table_seg.copy()
         for stain in dic.keys():
             expression_value = 1 if dic[stain] == "+" else 2
-            subset = subset.loc[subset[f"marker_{stain}"] == expression_value]
+            subset = subset[subset[f"marker_{stain}"] == expression_value]
 
-    label_ids_subtype = list(subset["label_id"])
+        label_ids_subtype.extend(list(subset["label_id"]))
     return label_ids_subtype
 
 
@@ -145,6 +153,7 @@ def export_lower_resolution(args):
     for subtype in subtypes:
 
         label_ids_subtype = filter_subtypes(cochlea, seg_name=seg_name, subtype=subtype, stains=subtype_stains)
+        print(f"Subtype '{subtype}' with {len(label_ids_subtype)} instances.")
         table.loc[table["label_id"].isin(label_ids_subtype), "subtype_label"] = subtype
 
     table.to_csv(out_path, sep="\t", index=False)

scripts/figures/SGNsub_thresholds.py scripts/figures/plot_SGNsub_thresholds.pyscripts/figures/SGNsub_thresholds.py renamed to scripts/figures/plot_SGNsub_thresholds.py

@@ -12,17 +12,24 @@
 png_dpi = 300
 
 COCHLEAE = {
-    "M_LR_000099_L": {"seg_data": "PV_SGN_v2", "subtype": ["Calb1", "Lypd1"], "intensity": "ratio"},
-    "M_AMD_N180_L": {"seg_data": "SGN_merged", "subtype": ["CR", "Lypd1", "Ntng1"], "intensity": "absolute"},
-    "M_AMD_N180_R": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"], "intensity": "absolute"},
     "M_LR_000098_L": {"seg_data": "SGN_v2", "subtype": ["CR", "Ntng1"], "intensity": "ratio"},
+    "M_LR_000099_L": {"seg_data": "PV_SGN_v2", "subtype": ["Calb1", "Lypd1"], "intensity": "ratio"},
     "M_LR_000184_L": {"seg_data": "SGN_v2", "subtype": ["Prph"], "output_seg": "SGN_v2b", "intensity": "ratio"},
     "M_LR_000184_R": {"seg_data": "SGN_v2", "subtype": ["Prph"], "output_seg": "SGN_v2b", "intensity": "ratio"},
     "M_LR_000214_L": {"seg_data": "PV_SGN_v2", "subtype": ["Calb1"], "intensity": "ratio"},
     "M_LR_000260_L": {"seg_data": "SGN_v2", "subtype": ["Prph", "Tuj1"], "intensity": "ratio"},
+    "M_LR_N110_L": {"seg_data": "SGN_v2", "subtype": ["Calb1", "Ntng1"], "intensity": "ratio"},
+    "M_LR_N110_R": {"seg_data": "SGN_v2", "subtype": ["Calb1", "Ntng1"], "intensity": "ratio"},
     "M_LR_N152_L": {"seg_data": "SGN_v2", "subtype": ["CR", "Ntng1"], "intensity": "ratio"},
+    "M_AMD_N180_L": {"seg_data": "SGN_merged", "subtype": ["CR", "Lypd1", "Ntng1"], "intensity": "absolute"},
+    "M_AMD_N180_R": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"], "intensity": "absolute"},
 }
 
+CUSTOM_THRESHOLDS = {
+    "M_LR_000099_L": {"Lypd1": 0.65},
+    "M_LR_000184_L": {"Prph": 1},
+    "M_LR_000260_L": {"Prph": 0.7},
+}
 
 def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False):
     """Plot histograms for positive and negative populations of subtype markers based on thresholding.
@@ -55,9 +62,15 @@ def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False):
             table_measurement_path = f"{cochlea}/tables/{data_name}/{stain}_{seg_str}_object-measures.tsv"
             column = "median"
 
-        rows = 1
-        columns = number_plots // rows
-        fig, ax = plt.subplots(rows, columns, figsize=(columns*4, rows*4), sharex=True)
+        # check for custom threshold
+        if CUSTOM_THRESHOLDS.get(cochlea, {}).get(stain) is not None:
+            rows = 2
+        else:
+            rows = 1
+
+        columns = number_plots
+        fig, axes = plt.subplots(rows, columns, figsize=(columns*4, rows*4), sharex=True)
+        ax = axes.flatten()
         table_path_s3, fs = get_s3_path(table_measurement_path)
         with fs.open(table_path_s3, "r") as f:
             table_measurement = pd.read_csv(f, sep="\t")
@@ -72,8 +85,8 @@ def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False):
             neg_values = list(subset_neg[column])
             pos_values = list(subset_pos[column])
 
-            bins = np.linspace(min(min(neg_values), min(pos_values)),
-                               max(max(neg_values), max(pos_values)), 30)
+            bins = np.linspace(min(neg_values + pos_values),
+                               max(neg_values + pos_values), 30)
 
             ax[num].hist(pos_values, bins=bins, alpha=0.6, label='Positive', color='tab:blue')
             ax[num].hist(neg_values, bins=bins, alpha=0.6, label='Negative', color='tab:orange')
@@ -82,6 +95,27 @@ def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False):
             ax[num].legend()
             ax[num].set_title(center_str)
 
+        if rows  == 2:
+            for num, center_str in enumerate(center_strs):
+                seg_ids = param_dicts[center_str]["seg_ids"]
+                threshold = CUSTOM_THRESHOLDS[cochlea][stain]
+
+                subset = table_measurement[table_measurement["label_id"].isin(seg_ids)]
+                subset_neg = subset[(subset[column] < threshold)]
+                subset_pos = subset[(subset[column] > threshold)]
+                neg_values = list(subset_neg[column])
+                pos_values = list(subset_pos[column])
+
+                bins = np.linspace(min(neg_values + pos_values),
+                                   max(neg_values + pos_values), 30)
+
+                ax[columns + num].hist(pos_values, bins=bins, alpha=0.6, label='Positive (custom)', color='tab:blue')
+                ax[columns + num].hist(neg_values, bins=bins, alpha=0.6, label='Negative (custom)', color='tab:orange')
+                ax[columns + num].set_ylabel('Count')
+                ax[columns + num].set_xlabel('Intensity')
+                ax[columns + num].legend()
+                ax[columns + num].set_title(center_str)
+
         fig.suptitle(f"{cochlea} - {stain}", fontsize=30)
 
         plt.tight_layout()
@@ -101,14 +135,14 @@ def main():
         description="Assign each segmentation instance a marker based on annotation thresholds."
     )
     parser.add_argument("-c", "--cochlea", type=str, nargs="+", default=COCHLEAE, help="Cochlea(e) to process.")
+    parser.add_argument("--figure_dir", "-f", type=str, required=True, help="Output directory for plots.")
     parser.add_argument("-i", "--input_dir", type=str,
-                        default="/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/mobie_project/cochlea-lightsheet/tables/Subtype_marker", help="Input directory.") #noqa
-    parser.add_argument("-o", "--output_dir", type=str,
-                        default="/user/schilling40/u15000/plots/SGNsub_thresholds", help="Output directory.")
+                        default="/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/mobie_project/cochlea-lightsheet/tables/Subtype_marker",
+                        help="Directory containing object measures of the cochleae.") #noqa
 
     args = parser.parse_args()
     for cochlea in args.cochlea:
-        plot_intensity_thresholds(args.input_dir, args.output_dir, cochlea)
+        plot_intensity_thresholds(args.input_dir, args.figure_dir, cochlea)
 
 
 if __name__ == "__main__":

scripts/measurements/evaluate_marker_annotations_subtype.py

@@ -13,15 +13,23 @@
 # The cochlea for the CHReef analysis.
 
 COCHLEAE = {
-    "M_LR_000099_L": {"seg_data": "PV_SGN_v2", "subtype": ["Calb1", "Lypd1"], "intensity": "ratio"},
-    "M_AMD_N180_L": {"seg_data": "SGN_merged", "subtype": ["CR", "Lypd1", "Ntng1"], "intensity": "absolute"},
-    "M_AMD_N180_R": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"], "intensity": "absolute"},
     "M_LR_000098_L": {"seg_data": "SGN_v2", "subtype": ["CR", "Ntng1"], "intensity": "ratio"},
+    "M_LR_000099_L": {"seg_data": "PV_SGN_v2", "subtype": ["Calb1", "Lypd1"], "intensity": "ratio"},
     "M_LR_000184_L": {"seg_data": "SGN_v2", "subtype": ["Prph"], "output_seg": "SGN_v2b", "intensity": "ratio"},
     "M_LR_000184_R": {"seg_data": "SGN_v2", "subtype": ["Prph"], "output_seg": "SGN_v2b", "intensity": "ratio"},
     "M_LR_000214_L": {"seg_data": "PV_SGN_v2", "subtype": ["Calb1"], "intensity": "ratio"},
     "M_LR_000260_L": {"seg_data": "SGN_v2", "subtype": ["Prph", "Tuj1"], "intensity": "ratio"},
+    "M_LR_N110_L": {"seg_data": "SGN_v2", "subtype": ["Calb1", "Ntng1"], "intensity": "ratio"},
+    "M_LR_N110_R": {"seg_data": "SGN_v2", "subtype": ["Calb1", "Ntng1"], "intensity": "ratio"},
     "M_LR_N152_L": {"seg_data": "SGN_v2", "subtype": ["CR", "Ntng1"], "intensity": "ratio"},
+    "M_AMD_N180_L": {"seg_data": "SGN_merged", "subtype": ["CR", "Lypd1", "Ntng1"], "intensity": "absolute"},
+    "M_AMD_N180_R": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"], "intensity": "absolute"},
+}
+
+CUSTOM_THRESHOLDS = {
+    "M_LR_000099_L": {"Lypd1": 0.65},
+    "M_LR_000184_L": {"Prph": 1},
+    "M_LR_000260_L": {"Prph": 0.7},
 }
 
 
@@ -44,7 +52,7 @@ def get_length_fraction_from_center(table, center_str):
     return length_fraction
 
 
-def apply_nearest_threshold(intensity_dic, table_seg, table_measurement, column="median", suffix="labels"):
+def apply_nearest_threshold(intensity_dic, table_seg, table_measurement, column="median", suffix="labels", custom_threshold=None):
     """Apply threshold to nearest segmentation instances.
     Crop centers are transformed into the "length fraction" parameter of the segmentation table.
     This avoids issues with the spiral shape of the cochlea and maps the assignment onto the Rosenthal"s canal.
@@ -54,7 +62,11 @@ def apply_nearest_threshold(intensity_dic, table_seg, table_measurement, column=
     for key in intensity_dic.keys():
         length_fraction = get_length_fraction_from_center(table_seg, key)
         intensity_dic[key]["length_fraction"] = length_fraction
-        lf_intensity[length_fraction] = {"threshold": intensity_dic[key]["median_intensity"]}
+        if custom_threshold is None:
+            lf_intensity[length_fraction] = {"threshold": intensity_dic[key]["median_intensity"]}
+        else:
+            print(f"Using custom threshold {custom_threshold} for crop {key}.")
+            lf_intensity[length_fraction] = {"threshold": custom_threshold}
 
     # get limits for checking marker thresholds
     lf_intensity = dict(sorted(lf_intensity.items()))
@@ -273,6 +285,7 @@ def evaluate_marker_annotation(
             else:
                 annot_table = pd.concat([annot_table, get_annotation_table(annot_dic, subtype)], ignore_index=True)
 
+            # create dictionary containing median intensity and segmentation ids for every crop
             om_dic = get_object_measures(annot_dic, intensity_dic, intensity_mode, subtype)
             om_out_path = os.path.join(output_dir, f"{cochlea_str}_{subtype}_om.json")
             with open(om_out_path, "w") as f:
@@ -293,8 +306,13 @@ def evaluate_marker_annotation(
                     table_measurement = pd.read_csv(f, sep="\t")
 
             # Apply the threshold to all SGNs.
+            if CUSTOM_THRESHOLDS.get(cochlea, {}).get(subtype) is not None:
+                custom_threshold = CUSTOM_THRESHOLDS[cochlea][subtype]
+            else:
+                custom_threshold = None
             table_seg = apply_nearest_threshold(
                 intensity_dic, table_seg, table_measurement, column=column, suffix=subtype,
+                custom_threshold=custom_threshold
             )
 
         # Save the table with positives / negatives for all SGNs.