Skip to content

Commit 98ebea9

Browse files
schilling40schilling40
committed
Custom thresholds for SGN subtypes
1 parent 004f0ba commit 98ebea9

File tree

6 files changed

+130
-146
lines changed

6 files changed

+130
-146
lines changed

‎flamingo_tools/segmentation/sgn_subtype_utils.py‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
import os
2+
3+
CUSTOM_THRESHOLDS = {
4+
"M_LR_000098_L": {"Ntng1": {
5+
"0492-0500-0581": {"manual": 2, "annotations": 1.41},
6+
"0561-0179-0502": {"manual": 2.05, "annotations": 1.33},
7+
"0604-0841-0296": {"manual": 2.32, "annotations": 2.32},
8+
"0901-0898-0720": {"manual": 1.82, "annotations": 1.82},
9+
"0929-0749-0290": {"manual": 2.78, "annotations": 2.78},
10+
"0956-0629-0576": {"manual": 3.36, "annotations": 3.36},
11+
}},
12+
"M_LR_000099_L": {"Lypd1": 0.65},
13+
"M_LR_000184_L": {"Prph": 1},
14+
"M_LR_000184_R": {"Prph": {
15+
"0549-0793-0901": {"manual": 1.57, "annotations": 1.57},
16+
"0618-0900-0594": {"manual": 1.71, "annotations": 1.71},
17+
"0786-0525-0928": {"manual": 1.48, "annotations": 1.48},
18+
"0870-0230-0729": {"manual": 1.35, "annotations": 1.08},
19+
"0911-1142-0849": {"manual": 1.2, "annotations": 1.05},
20+
"0970-0883-0666": {"manual": 1.83, "annotations": 1.83},
21+
}},
22+
"M_LR_000260_L": {"Prph": 0.7},
23+
"M_LR_N110_L": {"Ntng1": {
24+
"0592-0600-0527": {"manual": 1.71, "annotations": 1.71},
25+
"0594-0611-0870": {"manual": 1.9, "annotations": 1.64},
26+
"0653-1124-0730": {"manual": 1.9, "annotations": 1.56},
27+
"0840-0887-0878": {"manual": 1.87, "annotations": 1.87},
28+
"0966-0560-0452": {"manual": 2, "annotations": 1.71},
29+
"1213-0301-0584": {"manual": 2.33, "annotations": 2.33},
30+
}},
31+
"M_LR_N152_L": {"Ntng1": {
32+
"0519-0929-0610": {"manual": 2, "annotations": 1.95},
33+
"0752-0997-0849": {"manual": 2, "annotations": 3.51},
34+
"0791-1331-0750": {"manual": 2, "annotations": 1.15},
35+
"0794-0904-0420": {"manual": 2, "annotations": 1.66},
36+
"1013-0330-0643": {"manual": 2, "annotations": 1.45},
37+
"1013-0613-0479": {"manual": 2, "annotations": 3.83},
38+
39+
}},
40+
"M_AMD_N180_L": {"Lypd1": {
41+
"0441-0660-0521": {"manual": 200, "annotations": 186},
42+
"0510-0660-0850": {"manual": 179, "annotations": 179.0},
43+
"0578-1095-0560": {"manual": 200, "annotations": 168.0},
44+
"0728-0809-0463": {"manual": 220, "annotations": 192.0},
45+
}},
46+
}
47+
48+
49+
STAIN_TO_TYPE = {
50+
# Combinations of Calb1 and CR:
51+
"CR+/Calb1+": "Type Ib",
52+
"CR-/Calb1+": "Type IbIc", # Calb1 is expressed at Ic less than Lypd1 but more then CR
53+
"CR+/Calb1-": "Type Ia",
54+
"CR-/Calb1-": "Type II",
55+
56+
# Combinations of Calb1 and Lypd1:
57+
"Calb1+/Lypd1+": "Type IbIc",
58+
"Calb1+/Lypd1-": "Type Ib",
59+
"Calb1-/Lypd1+": "Type Ic",
60+
"Calb1-/Lypd1-": "inconclusive", # Can be Type Ia or Type II
61+
62+
# Combinations of Prph and Tuj1:
63+
"Prph+/Tuj1+": "Type II",
64+
"Prph+/Tuj1-": "Type II",
65+
"Prph-/Tuj1+": "Type I",
66+
"Prph-/Tuj1-": "inconclusive",
67+
68+
# Prph is isolated.
69+
"Prph+": "Type II",
70+
"Prph-": "Type I",
71+
72+
# Combinations of CR and Ntng1
73+
"CR+/Ntng1+": "Type Ib",
74+
"CR+/Ntng1-": "Type Ia",
75+
"CR-/Ntng1+": "Type Ic",
76+
"CR-/Ntng1-": "Type II",
77+
78+
# Combinations of CR and Lypd1
79+
"CR+/Lypd1+": "Type Ib",
80+
"CR+/Lypd1-": "Type Ia",
81+
"CR-/Lypd1+": "Type Ic",
82+
"CR-/Lypd1-": "Type II",
83+
84+
# Combinations of Calb1 and Ntng1
85+
"Calb1+/Ntng1+": "Type Ib",
86+
"Calb1+/Ntng1-": "Type Ia",
87+
"Calb1-/Ntng1+": "Type Ic",
88+
"Calb1-/Ntng1-": "inconclusive",
89+
90+
}

‎scripts/assign_subtypes.py‎

Lines changed: 12 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
import pandas as pd
55

66
from flamingo_tools.s3_utils import get_s3_path, BUCKET_NAME, SERVICE_ENDPOINT
7+
from flamingo_tools.segmentation.sgn_subtype_utils import STAIN_TO_TYPE
78
# from skimage.segmentation import relabel_sequential
89

910
COCHLEA_DICT = {
@@ -15,48 +16,11 @@
1516
"M_LR_N110_L": {"seg_data": "SGN_v2", "subtype": ["Calb1", "Ntng1"]},
1617
"M_LR_N110_R": {"seg_data": "SGN_v2", "subtype": ["Calb1", "Ntng1"]},
1718
"M_LR_N152_L": {"seg_data": "SGN_v2", "subtype": ["CR", "Ntng1"]},
18-
"M_AMD_N180_L": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"]},
19+
"M_AMD_N180_L": {"seg_data": "SGN_merged", "subtype": ["CR", "Lypd1"]},
1920
"M_AMD_N180_R": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"]},
2021
}
2122

2223

23-
STAIN_TO_TYPE = {
24-
# Combinations of Calb1 and CR:
25-
"CR+/Calb1+": "Type Ib",
26-
"CR-/Calb1+": "Type IbIc", # Calb1 is expressed at Ic less than Lypd1 but more then CR
27-
"CR+/Calb1-": "Type Ia",
28-
"CR-/Calb1-": "Type II",
29-
30-
# Combinations of Calb1 and Lypd1:
31-
"Calb1+/Lypd1+": "Type IbIc",
32-
"Calb1+/Lypd1-": "Type Ib",
33-
"Calb1-/Lypd1+": "Type Ic",
34-
"Calb1-/Lypd1-": "inconclusive", # Can be Type Ia or Type II
35-
36-
# Combinations of Prph and Tuj1:
37-
"Prph+/Tuj1+": "Type II",
38-
"Prph+/Tuj1-": "Type II",
39-
"Prph-/Tuj1+": "Type I",
40-
"Prph-/Tuj1-": "inconclusive",
41-
42-
# Prph is isolated.
43-
"Prph+": "Type II",
44-
"Prph-": "Type I",
45-
46-
# Combinations of CR and Ntng1
47-
"CR+/Ntng1+": "Type Ib",
48-
"CR+/Ntng1-": "Type Ia",
49-
"CR-/Ntng1+": "Type Ic",
50-
"CR-/Ntng1-": "inconclusive",
51-
52-
# Combinations of Calb1 and Ntng1
53-
"Calb1+/Ntng1+": "Type Ib",
54-
"Calb1+/Ntng1-": "inconclusive",
55-
"Calb1-/Ntng1+": "Type Ic",
56-
"Calb1-/Ntng1-": "inconclusive",
57-
}
58-
59-
6024
def types_for_stain(stains):
6125
stains.sort()
6226
assert len(stains) in (1, 2)
@@ -129,14 +93,13 @@ def filter_subtypes(cochlea, seg_name, subtype, stains=None):
12993
return label_ids_subtype
13094

13195

132-
def export_lower_resolution(args):
96+
def export_lower_resolution(cochlea, output_folder, subtype_column="subtype_label"):
13397

134-
cochlea = args.cochlea
13598
subtype_stains = COCHLEA_DICT[cochlea]["subtype"]
13699
subtype_stains.sort()
137100
seg_name = COCHLEA_DICT[cochlea]["seg_data"]
138101

139-
out_path = os.path.join(args.output_folder, f"{cochlea}_subtypes.tsv")
102+
out_path = os.path.join(output_folder, f"{cochlea}_subtypes.tsv")
140103

141104
table_seg_path = f"{cochlea}/tables/{seg_name}/default.tsv"
142105
table_path_s3, fs = get_s3_path(table_seg_path)
@@ -149,23 +112,26 @@ def export_lower_resolution(args):
149112

150113
# Subtype labels
151114
subtype_labels = ["None" for _ in range(len(table))]
152-
table["subtype_label"] = subtype_labels
115+
table[subtype_column] = subtype_labels
153116
for subtype in subtypes:
154117

155118
label_ids_subtype = filter_subtypes(cochlea, seg_name=seg_name, subtype=subtype, stains=subtype_stains)
156119
print(f"Subtype '{subtype}' with {len(label_ids_subtype)} instances.")
157-
table.loc[table["label_id"].isin(label_ids_subtype), "subtype_label"] = subtype
120+
table.loc[table["label_id"].isin(label_ids_subtype), subtype_column] = subtype
158121

159122
table.to_csv(out_path, sep="\t", index=False)
160123

161124

162125
def main():
163126
parser = argparse.ArgumentParser()
164-
parser.add_argument("--cochlea", "-c", required=True)
165-
parser.add_argument("--output_folder", "-o", required=True)
127+
parser.add_argument("-c", "--cochlea", type=str, nargs="+", required=True, help="Cochlea(e) to process.")
128+
parser.add_argument("-o", "--output_folder", required=True)
129+
parser.add_argument("-s", "--subtype_column", default="subtype_label",
130+
help="Custom name for column in output which features subtypes.")
166131
args = parser.parse_args()
167132

168-
export_lower_resolution(args)
133+
for cochlea in args.cochlea:
134+
export_lower_resolution(cochlea, args.output_folder, args.subtype_column)
169135

170136

171137
if __name__ == "__main__":

‎scripts/export_lower_resolution_subtypes.py‎

Lines changed: 1 addition & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
import zarr
88

99
from flamingo_tools.s3_utils import get_s3_path, BUCKET_NAME, SERVICE_ENDPOINT
10+
from flamingo_tools.segmentation.sgn_subtype_utils import STAIN_TO_TYPE
1011
# from skimage.segmentation import relabel_sequential
1112

1213
COCHLEA_DICT = {
@@ -18,37 +19,6 @@
1819
}
1920

2021

21-
STAIN_TO_TYPE = {
22-
# Combinations of Calb1 and CR:
23-
"CR+/Calb1+": "Type Ib",
24-
"CR-/Calb1+": "Type IbIc", # Calb1 is expressed at Ic less than Lypd1 but more then CR
25-
"CR+/Calb1-": "Type Ia",
26-
"CR-/Calb1-": "Type II",
27-
28-
# Combinations of Calb1 and Lypd1:
29-
"Calb1+/Lypd1+": "Type IbIc",
30-
"Calb1+/Lypd1-": "Type Ib",
31-
"Calb1-/Lypd1+": "Type Ic",
32-
"Calb1-/Lypd1-": "inconclusive", # Can be Type Ia or Type II
33-
34-
# Combinations of Prph and Tuj1:
35-
"Prph+/Tuj1+": "Type II",
36-
"Prph+/Tuj1-": "Type II",
37-
"Prph-/Tuj1+": "Type I",
38-
"Prph-/Tuj1-": "inconclusive",
39-
40-
# Prph is isolated.
41-
"Prph+": "Type II",
42-
"Prph-": "Type I",
43-
44-
# Combinations of CR and Ntng1
45-
"CR+/Ntng1+": "Type Ib",
46-
"CR+/Ntng1-": "Type Ia",
47-
"CR-/Ntng1+": "Type Ic",
48-
"CR-/Ntng1-": "inconclusive",
49-
}
50-
51-
5222
def types_for_stain(stains):
5323
stains.sort()
5424
assert len(stains) in (1, 2)

‎scripts/figures/plot_SGNsub_thresholds.py‎

Lines changed: 13 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
import pandas as pd
99

1010
from flamingo_tools.s3_utils import get_s3_path
11+
from flamingo_tools.segmentation.sgn_subtype_utils import CUSTOM_THRESHOLDS
1112

1213
png_dpi = 300
1314

@@ -25,14 +26,8 @@
2526
"M_AMD_N180_R": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"], "intensity": "absolute"},
2627
}
2728

28-
CUSTOM_THRESHOLDS = {
29-
"M_LR_000099_L": {"Lypd1": 0.65},
30-
"M_LR_000184_L": {"Prph": 1},
31-
"M_LR_000260_L": {"Prph": 0.7},
32-
}
33-
3429

35-
def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False):
30+
def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False, sharex=True):
3631
"""Plot histograms for positive and negative populations of subtype markers based on thresholding.
3732
"""
3833
os.makedirs(output_dir, exist_ok=True)
@@ -55,6 +50,7 @@ def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False):
5550
continue
5651
param_dicts = json.loads(data)
5752
center_strs = param_dicts["center_strings"]
53+
center_strs.sort()
5854
intensity_mode = param_dicts["intensity_mode"]
5955
number_plots = len(center_strs)
6056
if number_plots == 0:
@@ -70,13 +66,13 @@ def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False):
7066
column = "median"
7167

7268
# check for custom threshold
73-
if CUSTOM_THRESHOLDS.get(cochlea, {}).get(stain) is not None:
69+
if cochlea in CUSTOM_THRESHOLDS and stain in CUSTOM_THRESHOLDS[cochlea]:
7470
rows = 2
7571
else:
7672
rows = 1
7773

7874
columns = number_plots
79-
fig, axes = plt.subplots(rows, columns, figsize=(columns*4, rows*4), sharex=True)
75+
fig, axes = plt.subplots(rows, columns, figsize=(columns*4, rows*4), sharex=sharex)
8076
ax = axes.flatten()
8177
table_path_s3, fs = get_s3_path(table_measurement_path)
8278
with fs.open(table_path_s3, "r") as f:
@@ -103,9 +99,13 @@ def plot_intensity_thresholds(input_dir, output_dir, cochlea, plot=False):
10399
ax[num].set_title(center_str)
104100

105101
if rows == 2:
102+
threshold_dic = CUSTOM_THRESHOLDS[cochlea][stain]
106103
for num, center_str in enumerate(center_strs):
107104
seg_ids = param_dicts[center_str]["seg_ids"]
108-
threshold = CUSTOM_THRESHOLDS[cochlea][stain]
105+
if isinstance(threshold_dic, (int, float)):
106+
threshold = threshold_dic
107+
else:
108+
threshold = threshold_dic[center_str]["manual"]
109109

110110
subset = table_measurement[table_measurement["label_id"].isin(seg_ids)]
111111
subset_neg = subset[(subset[column] < threshold)]
@@ -143,13 +143,14 @@ def main():
143143
)
144144
parser.add_argument("-c", "--cochlea", type=str, nargs="+", default=COCHLEAE, help="Cochlea(e) to process.")
145145
parser.add_argument("--figure_dir", "-f", type=str, required=True, help="Output directory for plots.")
146+
parser.add_argument("--sharex", action="store_true", help="Shared x-axis of subplots.")
146147
parser.add_argument("-i", "--input_dir", type=str,
147148
default="/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/mobie_project/cochlea-lightsheet/tables/Subtype_marker", # noqa
148-
help="Directory containing object measures of the cochleae.") # noqa
149+
help="Directory containing object measures of the cochleae.")
149150

150151
args = parser.parse_args()
151152
for cochlea in args.cochlea:
152-
plot_intensity_thresholds(args.input_dir, args.figure_dir, cochlea)
153+
plot_intensity_thresholds(args.input_dir, args.figure_dir, cochlea, sharex=args.sharex)
153154

154155

155156
if __name__ == "__main__":

‎scripts/measurements/evaluate_marker_annotations_subtype.py‎

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
from flamingo_tools.s3_utils import get_s3_path
99
from flamingo_tools.file_utils import read_image_data
1010
from flamingo_tools.segmentation.chreef_utils import localize_median_intensities, find_annotations
11+
from flamingo_tools.segmentation.sgn_subtype_utils import CUSTOM_THRESHOLDS
1112

1213
MARKER_DIR_SUBTYPE = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/SGN_subtypes"
1314
# The cochlea for the CHReef analysis.
@@ -26,12 +27,6 @@
2627
"M_AMD_N180_R": {"seg_data": "SGN_merged", "subtype": ["CR", "Ntng1"], "intensity": "absolute"},
2728
}
2829

29-
CUSTOM_THRESHOLDS = {
30-
"M_LR_000099_L": {"Lypd1": 0.65},
31-
"M_LR_000184_L": {"Prph": 1},
32-
"M_LR_000260_L": {"Prph": 0.7},
33-
}
34-
3530

3631
def get_length_fraction_from_center(table, center_str):
3732
"""Get 'length_fraction' parameter for center coordinate by averaging nearby segmentation instances.
@@ -52,7 +47,7 @@ def get_length_fraction_from_center(table, center_str):
5247
return length_fraction
5348

5449

55-
def apply_nearest_threshold(intensity_dic, table_seg, table_measurement, column="median", suffix="labels", custom_threshold=None):
50+
def apply_nearest_threshold(intensity_dic, table_seg, table_measurement, column="median", suffix="labels", threshold_dic=None):
5651
"""Apply threshold to nearest segmentation instances.
5752
Crop centers are transformed into the "length fraction" parameter of the segmentation table.
5853
This avoids issues with the spiral shape of the cochlea and maps the assignment onto the Rosenthal"s canal.
@@ -62,9 +57,13 @@ def apply_nearest_threshold(intensity_dic, table_seg, table_measurement, column=
6257
for key in intensity_dic.keys():
6358
length_fraction = get_length_fraction_from_center(table_seg, key)
6459
intensity_dic[key]["length_fraction"] = length_fraction
65-
if custom_threshold is None:
60+
if threshold_dic is None:
6661
lf_intensity[length_fraction] = {"threshold": intensity_dic[key]["median_intensity"]}
6762
else:
63+
if isinstance(threshold_dic, (int, float)):
64+
custom_threshold = threshold_dic
65+
else:
66+
custom_threshold = threshold_dic[key]["manual"]
6867
print(f"Using custom threshold {custom_threshold} for crop {key}.")
6968
lf_intensity[length_fraction] = {"threshold": custom_threshold}
7069

@@ -307,12 +306,13 @@ def evaluate_marker_annotation(
307306

308307
# Apply the threshold to all SGNs.
309308
if CUSTOM_THRESHOLDS.get(cochlea, {}).get(subtype) is not None:
310-
custom_threshold = CUSTOM_THRESHOLDS[cochlea][subtype]
309+
custom_threshold_dic = CUSTOM_THRESHOLDS[cochlea][subtype]
311310
else:
312-
custom_threshold = None
311+
custom_threshold_dic = None
312+
313313
table_seg = apply_nearest_threshold(
314314
intensity_dic, table_seg, table_measurement, column=column, suffix=subtype,
315-
custom_threshold=custom_threshold
315+
threshold_dic=custom_threshold_dic,
316316
)
317317

318318
# Save the table with positives / negatives for all SGNs.

0 commit comments

Comments
 (0)