Sgn density and volume analysis (#78)

* Implement SGN density analysis

* Add volume analysis

* Update density analysis

* Update rankovic analysis

* Add IHC density analysis

* Update SGN density plots

* Update alias in density and volume scripts

Author: constantinpape

flamingo_tools/measurements.py

@@ -43,7 +43,7 @@ def _get_bounding_box_and_center(table, seg_id, resolution, shape, dilation):
     if dilation is not None and dilation > 0:
         bb_extension = dilation + 1
     else:
-        bb_extension = 1
+        bb_extension = 2
 
     bb_min = np.array([
         row.bb_min_z.item(), row.bb_min_y.item(), row.bb_min_x.item()
@@ -169,6 +169,21 @@ def _default_object_features(
     return measures
 
 
+def _morphology_features(seg_id, table, image, segmentation, resolution, **kwargs):
+    measures = {"label_id": seg_id}
+
+    bb, center = _get_bounding_box_and_center(table, seg_id, resolution, image.shape, dilation=0)
+    mask = segmentation[bb] == seg_id
+
+    # Hard-coded value for LaVision cochleae. This is a hack for the wrong voxel size in MoBIE.
+    # resolution = (3.0, 0.76, 0.76)
+
+    volume, surface = _measure_volume_and_surface(mask, resolution)
+    measures["volume"] = volume
+    measures["surface"] = surface
+    return measures
+
+
 def _regionprops_features(seg_id, table, image, segmentation, resolution, background_mask=None, dilation=None):
     bb, _ = _get_bounding_box_and_center(table, seg_id, resolution, image.shape, dilation)
 
@@ -222,6 +237,7 @@ def get_object_measures_from_table(arr_seg, table):
     "skimage": _regionprops_features,
     "default_background_norm": partial(_default_object_features, background_radius=75, norm=np.divide),
     "default_background_subtract": partial(_default_object_features, background_radius=75, norm=np.subtract),
+    "morphology": _morphology_features,
 }
 """The different feature functions that are supported in `compute_object_measures` and
 that can be selected via the feature_set argument. Currently this supports:

scripts/figures/plot_fig6.py

@@ -168,6 +168,7 @@ def fig_06e_octave(otof_data, save_path, plot=False, use_alias=True, trendline_m
         marker_neg = len([1 for i in marker_labels if i == 2])
         expression_eff = marker_pos / (marker_pos + marker_neg)
         print(f"Cochlea {name}, average expression efficiency {expression_eff}")
+        print(f"Cochlea {name}, number of IHCs: {len(freq)}")
         octave_binned = frequency_mapping(
             freq, marker_labels, animal="mouse", transduction_efficiency=True,
             bin_edges=bin_edges, bin_labels=bin_labels

scripts/figures/util.py

@@ -168,7 +168,8 @@ def _get_mapping(animal):
 
 
 def frequency_mapping(
-    frequencies, values, animal="mouse", transduction_efficiency=False, bin_edges=None, bin_labels=None
+    frequencies, values, animal="mouse", transduction_efficiency=False,
+    bin_edges=None, bin_labels=None, aggregation="mean",
 ):
     # Get the mapping of frequencies to octave bands for the given species.
     if bin_edges is None:
@@ -186,11 +187,8 @@ def frequency_mapping(
         num_tot = df[df["value"].isin([1, 2])].groupby("octave_band", observed=False).size()
         value_by_band = (num_pos / num_tot).reindex(bin_labels)
     else:  # Otherwise, aggregate the values over the octave band using the mean.
-        value_by_band = (
-            df.groupby("octave_band", observed=True)["value"]
-              .mean()
-              .reindex(bin_labels)   # keep octave order even if a bin is empty
-        )
+        aggregator = getattr(df.groupby("octave_band", observed=True)["value"], aggregation)
+        value_by_band = aggregator().reindex(bin_labels)  # keep octave order even if a bin is empty
     return value_by_band
 
 

scripts/la-vision/otof/README.md

@@ -0,0 +1,19 @@
+# Analysis Reuslts OTOF La-Vision Data
+
+## From Rankovic et al.
+
+- OTOF-23R:
+    - N-IHCS: 340
+    - Expression Efficiency: 0.3058823529411765
+- OTOF-25R:
+    - N-IHCS: 346
+    - Expression Efficiency: 0.21676300578034682
+
+## Our Analysis
+
+- OTOF-23R:
+    - N-IHCS: 643
+    - Expression Efficiency: 0.38202247191011235
+- OTOF-25R:
+    - N-IHCS: 520
+    - Expression Efficiency: 0.37884615384615383

scripts/la-vision/otof/comparison_rankovic.py

@@ -0,0 +1,108 @@
+import sys
+
+import matplotlib.pyplot as plt
+
+sys.path.append("../../figures")
+
+
+def plot_tonotopy(bin_labels, values, names, expression_eff_list):
+    from util import prism_style, prism_cleanup_axes
+
+    color_dict = {"O1": "#9C5027", "O2": "#67279C"}
+
+    prism_style()
+    label_size = 20
+    tick_label_size = 14
+
+    # result = pd.DataFrame(result)
+    # bin_labels = pd.unique(result["octave_band"])
+    # band_to_x = {band: i for i, band in enumerate(bin_labels)}
+    # result["x_pos"] = result["octave_band"].map(band_to_x)
+
+    fig, ax = plt.subplots(figsize=(8, 4))
+
+    offset = 0.08
+    for num in range(len(names)):
+
+        name = names[num]
+        expr_vals = values[num]
+
+        x_positions = range(len(bin_labels))
+        x_positions = [x - len(bin_labels) // 2 * offset + offset * num for x in x_positions]
+
+        x_positions = [pos for i, pos in enumerate(x_positions) if expr_vals[i] is not None]
+        expr_vals = [val for val in expr_vals if val is not None]
+        assert len(x_positions) == len(expr_vals)
+
+        ax.scatter(x_positions, expr_vals, marker="o", label=name, s=80, alpha=1, color=color_dict[name])
+
+    xlim_left, xlim_right = ax.get_xlim()
+    y_offset = [0.01, -0.04]
+    x_offset = 0.5
+    plt.xlim(xlim_left, xlim_right)
+    for num, key in enumerate(color_dict.keys()):
+        color = color_dict[key]
+        expression_eff = expression_eff_list[num]
+
+        ax.text(xlim_left + x_offset, expression_eff + y_offset[num], "mean",
+                color=color, fontsize=tick_label_size, ha="center")
+        trend_r, = ax.plot(
+            [xlim_left, xlim_right],
+            [expression_eff, expression_eff],
+            linestyle="dashed",
+            color=color,
+            alpha=0.7,
+            zorder=0
+        )
+
+    ax.set_xticks(range(len(bin_labels)))
+    ax.set_xticklabels(bin_labels)
+    ax.set_xlabel("Octave band [kHz]", fontsize=label_size)
+
+    ax.set_ylabel("Expression efficiency")
+    plt.tight_layout()
+    prism_cleanup_axes(ax)
+
+    plt.show()
+
+
+def main():
+    frequencies = ["4-8", "8-12", "12-16", "16-24", "24-32"]
+
+    n_ihcs_23R = [51, 107, 77, 105, None]
+    n_pos_23R = [15, 36, 24, 29, None]
+    assert len(n_ihcs_23R) == len(n_pos_23R)
+    assert len(n_ihcs_23R) == len(frequencies)
+
+    n_ihcs_25R = [None, 103, 71, 102, 70]
+    n_pos_25R = [None, 24, 21, 28, 2]
+    assert len(n_ihcs_25R) == len(n_pos_25R)
+    assert len(n_ihcs_25R) == len(frequencies)
+
+    total_23R = sum(n for n in n_ihcs_23R if n is not None)
+    pos_23R = sum(n for n in n_pos_23R if n is not None)
+    eff_23R = float(pos_23R) / total_23R
+    print("N-IHCs 23R:", total_23R)
+    print("Expr. Eff.:", eff_23R)
+    print()
+
+    total_25R = sum(n for n in n_ihcs_25R if n is not None)
+    pos_25R = sum(n for n in n_pos_25R if n is not None)
+    eff_25R = float(pos_25R) / total_25R
+    print("N-IHCs 25R:", total_25R)
+    print("Expr. Eff.:", eff_25R)
+
+    expr_23R = [None if n is None else float(pos) / n
+                for pos, n in zip(n_pos_23R, n_ihcs_23R)]
+    expr_25R = [None if n is None else float(pos) / n
+                for pos, n in zip(n_pos_25R, n_ihcs_25R)]
+    plot_tonotopy(
+        frequencies,
+        values=[expr_23R, expr_25R],
+        names=["O1", "O2"],
+        expression_eff_list=[eff_23R, eff_25R]
+    )
+
+
+if __name__ == "__main__":
+    main()