From 06d022963d3a92292bd2f104f2bd99c7b5fea01d Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Fri, 19 Sep 2025 08:21:01 +0200 Subject: [PATCH 1/2] Add script for checking IHCs without associated synapses --- .../check_empty_ihcs.py | 145 ++++++++++++++++++ 1 file changed, 145 insertions(+) create mode 100644 scripts/synapse_marker_detection/check_empty_ihcs.py diff --git a/scripts/synapse_marker_detection/check_empty_ihcs.py b/scripts/synapse_marker_detection/check_empty_ihcs.py new file mode 100644 index 0000000..7f8d6ee --- /dev/null +++ b/scripts/synapse_marker_detection/check_empty_ihcs.py @@ -0,0 +1,145 @@ +import os +from typing import Tuple + +# import imageio.v3 as imageio +import napari +import numpy as np +import pandas as pd +import s3fs +import zarr + +# +# NOTE: This was copied and adapted from 'flamingo_tools.s3_util' +# because of some import errors I have on the jupyter workstation +# that prevent me from importig from the flamingo_tools package. +# + +SERVICE_ENDPOINT = "https://s3.fs.gwdg.de/" +BUCKET_NAME = "cochlea-lightsheet" +CREDENTIAL_FILE = "/mnt/vast-nhr/home/pape41/u12086/.aws/credentials" + + +def read_s3_credentials(credential_file: str) -> Tuple[str, str]: + access_key, secret_key = None, None + with open(credential_file) as f: + for line in f: + if line.startswith("aws_access_key_id"): + access_key = line.rstrip("\n").strip().split(" ")[-1] + if line.startswith("aws_secret_access_key"): + secret_key = line.rstrip("\n").strip().split(" ")[-1] + if access_key is None or secret_key is None: + raise ValueError(f"Invalid credential file {credential_file}") + return access_key, secret_key + + +def _get_fs(): + client_kwargs = {"endpoint_url": SERVICE_ENDPOINT} + key, secret = read_s3_credentials(CREDENTIAL_FILE) + s3_filesystem = s3fs.S3FileSystem(key=key, secret=secret, client_kwargs=client_kwargs) + return s3_filesystem + + +# +# s3 utils go until here. +# + + +def get_ihcs_wo_syn(cochlea, ihc_name, syn_name): + """This function returns the table of IHCs without synapse. + """ + fs = _get_fs() + + table_path = os.path.join(BUCKET_NAME, cochlea, "tables", ihc_name, "default.tsv") + with fs.open(table_path, "r") as f: + ihc_table = pd.read_csv(f, sep="\t") + ihc_table = ihc_table[ihc_table.component_labels == 1] + ihc_ids = ihc_table.label_id.values + + table_path = os.path.join(BUCKET_NAME, cochlea, "tables", syn_name, "default.tsv") + with fs.open(table_path, "r") as f: + syn_table = pd.read_csv(f, sep="\t") + + matched_ihc = syn_table.matched_ihc.values + matched_ihc = matched_ihc[np.isin(matched_ihc, ihc_ids)] + + ihc_ids_with_syn, syn_count = np.unique(matched_ihc, return_counts=True) + ihc_ids_wo_syn_id = np.setdiff1d(ihc_ids, ihc_ids_with_syn) + ihc_ids_wo_syn = ihc_table[ihc_table.label_id.isin(ihc_ids_wo_syn_id)] + return ihc_ids_wo_syn + + +def _get_mobie_ds(cochlea, data_name): + internal_path = os.path.join(BUCKET_NAME, cochlea, "images", "ome-zarr", f"{data_name}.ome.zarr") + s3_filesystem = _get_fs() + print(internal_path) + s3_path = zarr.storage.FSStore(internal_path, fs=s3_filesystem) + return zarr.open(s3_path, mode="r")["s0"] + + +def check_empty_ihcs(cochlea, ihc_name, syn_name, pred_path=None): + ihc_ids_wo_syn = get_ihcs_wo_syn(cochlea, ihc_name, syn_name) + + halo = (96, 192, 192) + resolution = [0.38] * 3 + + ds_vglut = _get_mobie_ds(cochlea, "Vglut3") + ds_ctbp2 = _get_mobie_ds(cochlea, "CTBP2") + ds_seg = _get_mobie_ds(cochlea, ihc_name) + if pred_path is not None: + ds_pred = zarr.open(pred_path, "r")["prediction"] + + for i, row in ihc_ids_wo_syn.iterrows(): + position = [row.anchor_z, row.anchor_y, row.anchor_x] + center = [int(pos / res) for pos, res in zip(position, resolution)] + bb = tuple(slice(ce - ha, ce + ha) for ce, ha in zip(center, halo)) + + ctbp2 = ds_ctbp2[bb] + vglut3 = ds_vglut[bb] + + seg = ds_seg[bb] + if pred_path is not None: + pred = ds_pred[bb] + + pos_string = "-".join(str(int(pos)) for pos in position) + + v = napari.Viewer() + v.add_image(vglut3) + v.add_image(ctbp2) + v.add_labels(seg) + if pred is not None: + v.add_image(pred) + v.title = f"ID: {row.label_id} @ {pos_string}" + napari.run() + + # TODO enable saving of the crops via a magicgui plugin + # if save_raw: + # imageio.imwrite(f"{cochlea}_VGlut3_{'-'.join(str(int(pos)) for pos in position)}.tif", vglut3) + # imageio.imwrite(f"{cochlea}_CTBP2_{'-'.join(str(int(pos)) for pos in position)}.tif", ctbp2) + # continue + + +def check_gerbil(): + cochlea = "G_LR_000233_R" + ihc_name = "IHC_v6" + syn_name = "synapse_v3_ihc_v6" + + pred_root = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/predictions/G_LR_000233_R/synapses_v3" + pred_path = os.path.join(pred_root, "predictions.zarr") + + check_empty_ihcs(cochlea, ihc_name, syn_name, pred_path=pred_path) + + +# For the Gerbil cochlea: +# Synapses are not detected in both of these areas due to the very low contrast +# of the CTBP2 signal +# positions = [ +# [620.0, 1050.0, 1760.0], +# [1116.2177082702497, 1255.3013776373136, 298.84620835303593] +# ] +def main(): + # TODO create an argparser for this instead, so that we can use it to inspect other cochleae. + check_gerbil() + + +if __name__ == "__main__": + main() From e9652e0fa64a9cd31b559306c90320cd73df0c4d Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Fri, 19 Sep 2025 08:29:29 +0200 Subject: [PATCH 2/2] Update script for empty IHC checking --- scripts/synapse_marker_detection/check_empty_ihcs.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/scripts/synapse_marker_detection/check_empty_ihcs.py b/scripts/synapse_marker_detection/check_empty_ihcs.py index 7f8d6ee..8a1cac7 100644 --- a/scripts/synapse_marker_detection/check_empty_ihcs.py +++ b/scripts/synapse_marker_detection/check_empty_ihcs.py @@ -16,7 +16,7 @@ SERVICE_ENDPOINT = "https://s3.fs.gwdg.de/" BUCKET_NAME = "cochlea-lightsheet" -CREDENTIAL_FILE = "/mnt/vast-nhr/home/pape41/u12086/.aws/credentials" +CREDENTIAL_FILE = os.path.expanduser("~/.aws/credentials") def read_s3_credentials(credential_file: str) -> Tuple[str, str]: @@ -71,13 +71,14 @@ def get_ihcs_wo_syn(cochlea, ihc_name, syn_name): def _get_mobie_ds(cochlea, data_name): internal_path = os.path.join(BUCKET_NAME, cochlea, "images", "ome-zarr", f"{data_name}.ome.zarr") s3_filesystem = _get_fs() - print(internal_path) s3_path = zarr.storage.FSStore(internal_path, fs=s3_filesystem) return zarr.open(s3_path, mode="r")["s0"] +# TODO also get the synapse detections and visualize them. def check_empty_ihcs(cochlea, ihc_name, syn_name, pred_path=None): ihc_ids_wo_syn = get_ihcs_wo_syn(cochlea, ihc_name, syn_name) + print(len(ihc_ids_wo_syn), "IHCs don't have synapses.") halo = (96, 192, 192) resolution = [0.38] * 3 @@ -92,6 +93,7 @@ def check_empty_ihcs(cochlea, ihc_name, syn_name, pred_path=None): position = [row.anchor_z, row.anchor_y, row.anchor_x] center = [int(pos / res) for pos, res in zip(position, resolution)] bb = tuple(slice(ce - ha, ce + ha) for ce, ha in zip(center, halo)) + print("Load ID:", row.label_id, "from", center) ctbp2 = ds_ctbp2[bb] vglut3 = ds_vglut[bb] @@ -100,13 +102,13 @@ def check_empty_ihcs(cochlea, ihc_name, syn_name, pred_path=None): if pred_path is not None: pred = ds_pred[bb] - pos_string = "-".join(str(int(pos)) for pos in position) + pos_string = "-".join(str(int(pos)) for pos in position[::-1]) v = napari.Viewer() v.add_image(vglut3) v.add_image(ctbp2) v.add_labels(seg) - if pred is not None: + if pred_path is not None: v.add_image(pred) v.title = f"ID: {row.label_id} @ {pos_string}" napari.run() @@ -125,6 +127,8 @@ def check_gerbil(): pred_root = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/predictions/G_LR_000233_R/synapses_v3" pred_path = os.path.join(pred_root, "predictions.zarr") + if not os.path.exists(pred_path): + pred_path = None check_empty_ihcs(cochlea, ihc_name, syn_name, pred_path=pred_path)