Implement mapping of synapses to IHC and create draft for marker detection workflow

Author: constantinpape

flamingo_tools/segmentation/marker_detection.py

@@ -0,0 +1,77 @@
+from typing import Optional
+
+import numpy as np
+import pandas as pd
+
+from elf.parallel.distance_transform import map_points_to_objects
+
+
+def map_and_filter_detections(
+    segmentation: np.ndarray,
+    detections: pd.DataFrame,
+    max_distance: float,
+    resolution: float = 0.38,
+    n_threads: Optional[int] = None,
+    verbose: bool = True,
+) -> pd.DataFrame:
+    """Map synapse detections to segmented IHCs and filter out detections above a distance threshold to the IHCs.
+
+    Args:
+        segmentation: The IHC segmentation.
+        detections: The synapse marker detections.
+        max_distance: The maximal distance for a valid match of synapse markers to IHCs.
+        resolution: The resolution / voxel size of the data in micrometer.
+        n_threads: The number of threads for parallelizing the mapping of detections to objects.
+        verbose: Whether to print the progress of the mapping procedure.
+    """
+    # Get the point coordinates.
+    points = detections[["z", "y", "x"]].values.astype("int")
+
+    # Set the block shape (this could also be exposed as a parameter; it should not matter much though).
+    block_shape = (64, 256, 256)
+
+    # Determine the halo. We set it to 2 pixels + the max-distance in pixels, to ensure all distances
+    # that are smaller than the max distance are measured.
+    halo = (2 + int(np.ceil(max_distance / resolution)),) * 3
+
+    # Map the detections to the obejcts in the (IHC) segmentation.
+    object_ids, object_distances = map_points_to_objects(
+        segmentation=segmentation,
+        points=points,
+        block_shape=block_shape,
+        halo=halo,
+        sampling=resolution,
+        n_threads=n_threads,
+        verbose=verbose,
+    )
+    assert len(object_ids) == len(points)
+    assert len(object_distances) == len(points)
+
+    # Add matched ids and distances to the dataframe.
+    detections["matched_ihc"] = object_ids
+    detections["distance_to_ihc"] = object_distances
+
+    # Filter the dataframe by the max distance.
+    detections = detections[detections.distance_to_ihc < max_distance]
+    return detections
+
+
+# TODO implement streamlined workflow for the marker detection, mapping and filtering.
+def marker_detection():
+    """
+    """
+
+    # 1.) Determine mask for inference based on the IHC segmentation.
+    # Best approach: load IHC segmentation at a low scale level, binarize it,
+    # dilate it and use this as mask. It can be mapped back to the full resolution
+    # with `elf.wrapper.ResizedVolume`.
+
+    # 2.) Run inference and detection of maxima.
+    # This can be taken from 'scripts/synapse_marker_detection/run_prediction.py'
+    # (And the run prediction script should then be refactored).
+
+    # 3.) Map the detections to IHC and filter them based on a distance criterion.
+    # Use the function 'map_and_filter_detections' from above.
+
+    # 4.) Add the filtered detections to MoBIE.
+    # IMPORTANT scale the coordinates with the resolution here.

scripts/ihc_analysis/synapse_mapping.py

@@ -0,0 +1,57 @@
+import os
+
+import imageio.v3 as imageio
+import pandas as pd
+
+from flamingo_tools.segmentation.marker_detection import map_and_filter_detections
+
+
+ROOT = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet/croppings/Synapse_crop"
+VGLUT_PATH = os.path.join(ROOT, "M_LR_000226_R_crop_1098-0926-0872_Vglut3.tif")
+CTBP2_PATH = os.path.join(ROOT, "M_LR_000226_R_crop_1098-0926-0872_CTBP2.tif")
+SEG_PATH = os.path.join(ROOT, "M_LR_000226_R_resized_crop_1098-0926-0872_IHC.tif")
+DET_PATH = os.path.join(ROOT, "synapses/synapse_detection.tsv")
+
+
+def check_data():
+    import napari
+
+    detections = pd.read_csv(DET_PATH, sep="\t")
+    detections = detections[["z", "y", "x"]].values
+
+    filtered_path = os.path.join(ROOT, "synapses/synapse_detection_filtered.tsv")
+    filtered_detections = pd.read_csv(filtered_path, sep="\t")
+    filtered_detections = filtered_detections[["z", "y", "x"]].values
+
+    vglut = imageio.imread(VGLUT_PATH)
+    ctbp2 = imageio.imread(CTBP2_PATH)
+    ihcs = imageio.imread(SEG_PATH)
+
+    v = napari.Viewer()
+    v.add_image(vglut)
+    v.add_image(ctbp2)
+    v.add_labels(ihcs)
+    v.add_points(detections)
+    v.add_points(filtered_detections)
+    napari.run()
+
+
+def map_synapses():
+    ihcs = imageio.imread(SEG_PATH)
+    detections = pd.read_csv(DET_PATH, sep="\t")
+    n_detections = len(detections)
+
+    detections = map_and_filter_detections(ihcs, detections, max_distance=2.0, n_threads=8)
+    print("Detections after mapping and fitering:", len(detections), "/", n_detections)
+
+    out_path = os.path.join(ROOT, "synapses/synapse_detection_filtered.tsv")
+    detections.to_csv(out_path, sep="\t", index=False)
+
+
+def main():
+    map_synapses()
+    # check_data()
+
+
+if __name__ == "__main__":
+    main()