Merge pull request #3 from ttngu207/main

added brain atlas and ontology loader

Author: kabilar

README.md

@@ -33,7 +33,8 @@ fully functional electrode localization pipeline.
 
 + Install `element-interface`
 
-    + `element-interface` is a dependency of `element-electrode-localization`, however it is not contained within `requirements.txt`.
+    + `element-interface` is a dependency of `element-electrode-localization`, however
+      it is not contained within `requirements.txt`.
      
     ```
     pip install "element-interface @ git+https://github.com/datajoint/element-interface"
@@ -53,12 +54,13 @@ To activate the `element-electrode-localization`, ones need to provide:
 
 3. Utility functions. See [example definitions](https://github.com/datajoint/workflow-array-ephys/blob/main/workflow_array_ephys/paths.py).
     + get_ephys_root_data_dir(): Returns your root data directory.
-    + get_session_directory(): Returns the path of the session data relative to the root.
+    + get_session_directory(): Returns the path of the session data relative to the
+      root.
 
 For more details, check the docstring of the `element-electrode-localization`:
-
+```python
     help(electrode.activate)
-
+```
 ### Example usage
 
 See the [workflow-array-ephys project](https://github.com/datajoint/workflow-array-ephys) for an example usage of this Element.
@@ -75,4 +77,4 @@ See the [workflow-array-ephys project](https://github.com/datajoint/workflow-arr
 + DataJoint Elements
     + Yatsenko D, Nguyen T, Shen S, Gunalan K, Turner CA, Guzman R, Sasaki M, Sitonic D, Reimer J, Walker EY, Tolias AS. DataJoint Elements: Data Workflows for Neurophysiology. bioRxiv. 2021 Jan 1. doi: https://doi.org/10.1101/2021.03.30.437358
 
-    + DataJoint Elements ([RRID:SCR_021894](https://scicrunch.org/resolver/SCR_021894)) - Element Electrode Localization (version `<Enter version number>`)
+    + DataJoint Elements ([RRID:SCR_021894](https://scicrunch.org/resolver/SCR_021894)) - Element Electrode Localization (version `<Enter version number>`)

element_electrode_localization/coordinate_framework.py

@@ -1,26 +1,41 @@
+import logging
+import pathlib
+from tqdm import tqdm
 import datajoint as dj
+import pandas as pd
+import numpy as np
+from tqdm import tqdm
+import nrrd
+import re
 
 
+log = logging.getLogger(__name__)
 schema = dj.schema()
 
 
 def activate(schema_name, *, create_schema=True, create_tables=True):
     """
     activate(schema_name, create_schema=True, create_tables=True)
-        :param schema_name: schema name on the database server to activate the `coordinate_framework` element
-        :param create_schema: when True (default), create schema in the database if it does not yet exist.
-        :param create_tables: when True (default), create tables in the database if they do not yet exist.
+        :param schema_name: schema name on the database server to activate the
+                            `coordinate_framework` element
+        :param create_schema: when True (default), create schema in the database if it
+                              does not yet exist.
+        :param create_tables: when True (default), create tables in the database if
+                              they do not yet exist.
     """
-    schema.activate(schema_name, create_schema=create_schema, create_tables=create_tables)
+    schema.activate(schema_name, create_schema=create_schema,
+                    create_tables=create_tables)
+
+    # ----------------------------- Table declarations ----------------------
 
 
 @schema
 class CCF(dj.Lookup):
     definition = """  # Common Coordinate Framework
-    # CCF Dataset Information
-    ccf_id:             int             # CCF ID
+    ccf_id            : int             # CCF ID, a.k.a atlas ID
     ---
-    ccf_version:        int             # Allen CCF Version - e.g. CCFv3
+    ccf_version       : varchar(64)     # Allen CCF Version - e.g. CCFv3
+    ccf_resolution    : float           # voxel resolution in micron
     ccf_description='': varchar(255)    # CCFLabel Description
     """
 
@@ -43,7 +58,7 @@ class BrainRegionAnnotation(dj.Lookup):
     class BrainRegion(dj.Part):
         definition = """
         -> master
-        acronym: varchar(32)
+        acronym: varchar(32)  # CHARACTER SET utf8 COLLATE utf8_bin
         ---
         region_name: varchar(128)
         region_id=null: int
@@ -56,13 +71,124 @@ class Voxel(dj.Part):
         -> CCF.Voxel
         """
 
+    @classmethod
+    def retrieve_acronym(self, acronym):
+        """ Retrieve the DataJoint translation of the CCF acronym"""
+        return re.sub(r'(?<!^)(?=[A-Z])', '_', acronym).lower()
+
+    @classmethod
+    def voxel_query(self, x=None, y=None, z=None):
+        """Given one or more coordinates, return unique brain regions
+        :param x: x coordinate
+        :param y: y coordinate
+        :param z: z coordinate
+        """
+        if not any(x, y, z):
+            raise ValueError('Must specify at least one dimension')
+        # query = self.Voxel  #  TODO: add utility function name lookup
+        raise NotImplementedError('Coming soon')
+
 
 @schema
 class ParentBrainRegion(dj.Lookup):
-    definition = """ # Hierarchical structure between the brain regions
+    definition = """ # Hierarchical structure between the brain regionss
     -> BrainRegionAnnotation.BrainRegion
     ---
     -> BrainRegionAnnotation.BrainRegion.proj(parent='acronym')
     """
 
 
+# ---- HELPERS ----
+
+
+def load_ccf_annotation(ccf_id, version_name, voxel_resolution,
+                        nrrd_filepath, ontology_csv_filepath):
+    """
+    :param ccf_id: unique id to identify a new CCF dataset to be inserted
+    :param version_name: CCF version
+    :param voxel_resolution: voxel resolution in micron
+    :param nrrd_filepath: path to the .nrrd file for the volume data
+    :param ontology_csv_filepath: path to the .csv file for the brain region ontology
+
+    load_ccf_annotation(
+        ccf_id=0, version_name='ccf_2017', voxel_resolution=10,
+        nrrd_filepath='./data/annotation_10.nrrd',
+        ontology_csv_filepath='./data/query.csv')
+
+    For an example Allen brain atlas for mouse, see:
+    http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/annotation/ccf_2017
+
+    For the structure/ontology tree, see:
+    https://community.brain-map.org/t/allen-mouse-ccf-accessing-and-using-related-data-and-tools/359
+    (particularly the ontology file downloadable as CSV)
+    """
+    ccf_key = {'ccf_id': ccf_id}
+    if CCF & ccf_key:
+        print(f'CCF ID {ccf_id} already exists!')
+        return
+
+    nrrd_filepath = pathlib.Path(nrrd_filepath)
+    ontology_csv_filepath = pathlib.Path(ontology_csv_filepath)
+
+    def to_snake_case(s):
+        return re.sub(r'(?<!^)(?=[A-Z])', '_', s).lower()
+
+    ontology = pd.read_csv(ontology_csv_filepath)
+
+    stack, hdr = nrrd.read(nrrd_filepath.as_posix())  # AP (x), DV (y), ML (z)
+
+    log.info('.. loaded atlas brain volume of shape '
+             + f'{stack.shape} from {nrrd_filepath}')
+
+    ccf_key = {'ccf_id': ccf_id}
+    ccf_entry = {**ccf_key,
+                 'ccf_version': version_name,
+                 'ccf_resolution': voxel_resolution,
+                 'ccf_description': (f'Version: {version_name}'
+                                     + f' - Voxel resolution (uM): {voxel_resolution}'
+                                     + f' - Volume file: {nrrd_filepath.name}'
+                                     + ' - Region ontology file: '
+                                     + ontology_csv_filepath.name)
+                 }
+
+    with dj.conn().transaction:
+        CCF.insert1(ccf_entry)
+        BrainRegionAnnotation.insert1(ccf_key)
+        BrainRegionAnnotation.BrainRegion.insert([
+            dict(ccf_id=ccf_id,
+                 acronym=to_snake_case(r.acronym),
+                 region_id=r.id,
+                 region_name=r.safe_name,
+                 color_code=r.color_hex_triplet) for _, r in ontology.iterrows()])
+
+        # Process voxels per brain region
+        for idx, (region_id, r) in enumerate(tqdm(ontology.iterrows())):
+            dj.conn().ping()
+            region_id = int(region_id)
+
+            log.info('.. loading region {} ({}/{}) ({})'
+                     .format(region_id, idx, len(ontology), r.safe_name))
+
+            # extracting filled volumes from stack in scaled [[x,y,z]] shape,
+            vol = (np.array(np.where(stack == region_id)).T * voxel_resolution)
+            vol = pd.DataFrame(vol, columns=['x', 'y', 'z'])
+
+            if not vol.shape[0]:
+                log.info('.. region {} volume: shape {} - skipping'
+                         .format(region_id, vol.shape))
+                continue
+            else:
+                log.info('.. region {} volume: shape {}'.format(
+                    region_id, vol.shape))
+
+            vol['ccf_id'] = [ccf_key['ccf_id']] * len(vol)
+            CCF.Voxel.insert(vol)
+
+            vol['acronym'] = [to_snake_case(r.acronym)] * len(vol)
+            BrainRegionAnnotation.Voxel.insert(vol)
+
+    log.info('.. done.')
+
+
+def load_parent_regions(ccf_id):
+    raise NotImplementedError('Coming soon')