Add scripts for initial nephridia analysis

analysis/nephridia.py

+import os
+import numpy as np
+import pandas as pd
+from scripts.analysis.nephridia import filter_by_size, plot_sizes
+from scripts.analysis.nephridia import match_cilia_to_cells
+
+cell_ids = [24449, 22584, 21904, 21590, 21594, 21595, 21910, 21911, 21915]
+print("Number of cells:", len(cell_ids))
+
+tmp_path = './cilia_ids.npy'
+if os.path.exists(tmp_path):
+    cilia_ids = np.load(tmp_path)
+else:
+    cell_table = '../data/0.5.2/tables/sbem-6dpf-1-whole-segmented-cells-labels/default.csv'
+    cilia_path = '../data/0.5.2/segmentations/sbem-6dpf-1-whole-segmented-cilia-labels.h5'
+    cilia_key = 't00000/s00/2/cells'
+    cilia_res = [0.1, 0.04, 0.04]
+    cilia_ids = match_cilia_to_cells(cell_ids, cell_table,
+                                     cilia_path, cilia_key, cilia_res)
+    np.save(tmp_path, cilia_ids)
+
+table_path = '../data/0.5.2/tables/sbem-6dpf-1-whole-segmented-cilia-labels/default.csv'
+table_path2 = '../data/0.5.2/tables/sbem-6dpf-1-whole-segmented-cilia-labels/cilia.csv'
+
+table1 = pd.read_csv(table_path, sep='\t')
+table2 = pd.read_csv(table_path2, sep='\t')
+table2 = table2[['length', 'diameter_mean']]
+table = pd.concat([table1, table2], axis=1)
+assert len(table1) == len(table2) == len(table)
+
+table.set_index('label_id')
+table = table[3:]
+table = table.loc[table['length'] > 0.1]
+ids = table['label_id'].values
+ids = ids[np.isin(ids, cilia_ids)]
+print(ids)
+table = table.loc[ids]
+
+# plot_sizes(table)
+size_threshold = 5000
+table = filter_by_size(table, size_threshold)
+print("Number of cilia:", len(table))
+
+lens = table["length"].values
+for idd, leng in zip(ids, lens):
+    print(idd, leng)
+print("Average len:", np.mean(lens), "+-", np.std(lens))
+
+diameters = table["diameter_mean"].values
+print("Average diameter:", np.mean(diameters), "+-", np.std(diameters))

scripts/analysis/nephridia.py

+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import h5py
+from concurrent import futures
+
+
+def get_bb(idd, table, res):
+    row = table.loc[idd]
+    bb_min = [row.bb_min_z, row.bb_min_y, row.bb_min_x]
+    bb_max = [row.bb_max_z, row.bb_max_y, row.bb_max_x]
+    return tuple(slice(int(mi / re), int(ma / re))
+                 for mi, ma, re in zip(bb_min, bb_max, res))
+
+
+def match_cilia_to_cells(cell_ids, cell_table,
+                         cilia_seg_path, cilia_seg_key, cilia_res):
+    cell_table = pd.read_csv(cell_table, sep='\t')
+    cell_table.set_index('label_id')
+
+    with h5py.File(cilia_seg_path, 'r') as f:
+        ds_cil = f[cilia_seg_key]
+
+        def match_single_cell(cell_id):
+            # get the bounding box of this cell
+            bb = get_bb(cell_id, cell_table, cilia_res)
+            seg = ds_cil[bb]
+            return np.unique(seg)
+
+        with futures.ThreadPoolExecutor(9) as tp:
+            tasks = [tp.submit(match_single_cell, cell_id) for cell_id in cell_ids]
+            res = [t.result() for t in tasks]
+
+        cilia_ids = np.concatenate(res)
+        cilia_ids = np.unique(cilia_ids)
+    return cilia_ids
+
+
+def plot_sizes(table):
+    sizes = table['n_pixels'].values[1:]
+    print(sizes.max(), sizes.min())
+
+    fig, ax = plt.subplots()
+    _, bins, patches = ax.hist(sizes, 32)
+    ax.set_xlabel("Size in pixel")
+    ax.set_ylabel("Count")
+    plt.show()
+
+    sizes = sizes[sizes <= bins[1]]
+    fig, ax = plt.subplots()
+    _, bins, patches = ax.hist(sizes, 32)
+    print(bins)
+    ax.set_xlabel("Size in pixel")
+    ax.set_ylabel("Count")
+    plt.show()
+
+
+def filter_by_size(table, size_threshold):
+    table = table.loc[table['n_pixels'] > size_threshold]
+    return table
+
+
+def compute_offsets(table):
+    df = table[['anchor_x', 'anchor_y', 'anchor_y']]
+    pos = df.values
+    offsets = np.linalg.norm(pos, axis=1)
+    return offsets
+
+
+def plot_offsets(table):
+    offsets = compute_offsets(table)
+    fig, ax = plt.subplots()
+    ax.hist(offsets, 32)
+    ax.set_xlabel("Offset in microns")
+    ax.set_ylabel("Count")
+    plt.show()
+
+
+def filter_by_offset(table, offset_threshold):
+    offsets = compute_offsets(table)
+    table = table.loc[offsets > offset_threshold]
+    return table