Configurable group ids region names (#55)

added `group_ids_config.json`, so that ids can be specified at runtime for densities

CHANGELOG.rst

@@ -1,6 +1,13 @@
 Changelog
 =========
 
+Version 0.2.1
+-------------
+* The following *cell-density* sub-commands can now optionally take a ``--group-ids-config-path``:
+   *cell-density*, *glia-cell-densities*, *inhibitory-and-excitatory-neuron-densities*, *fit-average-densities*
+
+   Otherwise they use the default config for the AIBS atlas.
+
 Version 0.2.0
 -------------
 

atlas_densities/app/cell_densities.py

@@ -54,12 +54,18 @@
     common_atlas_options,
     log_args,
     set_verbose,
+    verbose_option,
 )
 from atlas_commons.typing import FloatArray
 from voxcell import RegionMap, VoxelData  # type: ignore
 
 from atlas_densities.app.utils import AD_PATH, DATA_PATH
-from atlas_densities.densities import excitatory_inhibitory_splitting
+from atlas_densities.densities import (
+    excitatory_inhibitory_splitting,
+    inhibitory_neuron_densities_optimization,
+    refined_inhibitory_neuron_densities,
+    utils,
+)
 from atlas_densities.densities.cell_counts import (
     extract_inhibitory_neurons_dataframe,
     glia_cell_counts,
@@ -72,17 +78,11 @@
 )
 from atlas_densities.densities.fitting import linear_fitting
 from atlas_densities.densities.glia_densities import compute_glia_densities
-from atlas_densities.densities.inhibitory_neuron_densities_optimization import (
-    create_inhibitory_neuron_densities as linprog,
-)
 from atlas_densities.densities.inhibitory_neuron_density import compute_inhibitory_neuron_density
 from atlas_densities.densities.measurement_to_density import (
     measurement_to_average_density,
     remove_non_density_measurements,
 )
-from atlas_densities.densities.refined_inhibitory_neuron_densities import (
-    create_inhibitory_neuron_densities as keep_proportions,
-)
 from atlas_densities.exceptions import AtlasDensitiesError
 
 EXCITATORY_SPLIT_CORTEX_ALL_TO_EXC_MTYPES = (
@@ -95,8 +95,8 @@
 LINPROG_PATH = "doc/source/bbpp82_628_linprog.pdf"
 
 ALGORITHMS: Dict[str, Callable] = {
-    "keep-proportions": keep_proportions,
-    "linprog": linprog,
+    "keep-proportions": refined_inhibitory_neuron_densities.create_inhibitory_neuron_densities,
+    "linprog": inhibitory_neuron_densities_optimization.create_inhibitory_neuron_densities,
 }
 
 L = logging.getLogger(__name__)
@@ -155,7 +155,7 @@ def _get_voxel_volume_in_mm3(voxel_data: "VoxelData") -> float:
 
 
 @click.group()
-@click.option("-v", "--verbose", count=True)
+@verbose_option
 def app(verbose):
     """Run the cell densities CLI"""
     set_verbose(L, verbose)
@@ -177,13 +177,14 @@ def app(verbose):
     "A voxel value is a number of cells per mm^3",
 )
 @click.option(
-    "--root-region-name",
-    type=str,
-    default="root",
-    help="Name of the region in the hierarchy",
+    "--group-ids-config-path",
+    type=EXISTING_FILE_PATH,
+    default=utils.GROUP_IDS_PATH,
+    help="Path to density groups ids config",
+    show_default=True,
 )
 @log_args(L)
-def cell_density(annotation_path, hierarchy_path, nissl_path, output_path, root_region_name):
+def cell_density(annotation_path, hierarchy_path, nissl_path, output_path, group_ids_config_path):
     """Compute and save the overall mouse brain cell density.
 
     The input Nissl stain volume of AIBS is turned into an actual density field complying with
@@ -210,13 +211,14 @@ def cell_density(annotation_path, hierarchy_path, nissl_path, output_path, root_
     assert_properties([annotation, nissl])
 
     region_map = RegionMap.load_json(hierarchy_path)
+    group_ids_config = utils.load_json(group_ids_config_path)
 
     overall_cell_density = compute_cell_density(
         region_map,
         annotation.raw,
         _get_voxel_volume_in_mm3(annotation),
         nissl.raw,
-        root_region_name=root_region_name,
+        group_ids_config=group_ids_config,
     )
     nissl.with_data(overall_cell_density).save_nrrd(output_path)
 
@@ -268,6 +270,13 @@ def cell_density(annotation_path, hierarchy_path, nissl_path, output_path, root_
     help="Path to the directory where to write the output cell density nrrd files."
     " It will be created if it doesn't exist already.",
 )
+@click.option(
+    "--group-ids-config-path",
+    type=EXISTING_FILE_PATH,
+    default=utils.GROUP_IDS_PATH,
+    help="Path to density groups ids config",
+    show_default=True,
+)
 @log_args(L)
 def glia_cell_densities(
     annotation_path,
@@ -279,6 +288,7 @@ def glia_cell_densities(
     microglia_density_path,
     glia_proportions_path,
     output_dir,
+    group_ids_config_path,
 ):  # pylint: disable=too-many-arguments, too-many-locals
     """Compute and save the glia cell densities.
 
@@ -338,20 +348,20 @@ def glia_cell_densities(
         "microglia": VoxelData.load_nrrd(microglia_density_path),
     }
 
-    atlases = list(glia_densities.values())
-    atlases += [annotation, overall_cell_density]
+    atlases = list(glia_densities.values()) + [annotation, overall_cell_density]
     L.info("Checking input files consistency ...")
     assert_properties(atlases)
 
     L.info("Loading hierarchy ...")
     region_map = RegionMap.load_json(hierarchy_path)
-    with open(glia_proportions_path, "r", encoding="utf-8") as file_:
-        glia_proportions = json.load(file_)
+    glia_proportions = utils.load_json(glia_proportions_path)
 
     glia_densities = {
-        glia_cell_type: voxel_data.raw for (glia_cell_type, voxel_data) in glia_densities.items()
+        glia_cell_type: voxel_data.raw for glia_cell_type, voxel_data in glia_densities.items()
     }
 
+    group_ids_config = utils.load_json(group_ids_config_path)
+
     L.info("Compute volumetric glia densities: started")
     glia_densities = compute_glia_densities(
         region_map,
@@ -362,6 +372,7 @@ def glia_cell_densities(
         overall_cell_density.raw,
         glia_proportions,
         copy=False,
+        group_ids_config=group_ids_config,
     )
 
     if not Path(output_dir).exists():
@@ -424,10 +435,11 @@ def glia_cell_densities(
     " It will be created if it doesn't exist already.",
 )
 @click.option(
-    "--root-region-name",
-    type=str,
-    default="root",
-    help="Name of the region in the hierarchy",
+    "--group-ids-config-path",
+    type=EXISTING_FILE_PATH,
+    default=utils.GROUP_IDS_PATH,
+    help="Path to density groups ids config",
+    show_default=True,
 )
 @log_args(L)
 def inhibitory_and_excitatory_neuron_densities(
@@ -438,7 +450,7 @@ def inhibitory_and_excitatory_neuron_densities(
     neuron_density_path,
     inhibitory_neuron_counts_path,
     output_dir,
-    root_region_name,
+    group_ids_config_path,
 ):  # pylint: disable=too-many-arguments
     """Compute and save the inhibitory and excitatory neuron densities.
 
@@ -484,6 +496,7 @@ def inhibitory_and_excitatory_neuron_densities(
 
     region_map = RegionMap.load_json(hierarchy_path)
     inhibitory_df = extract_inhibitory_neurons_dataframe(inhibitory_neuron_counts_path)
+    group_ids_config = utils.load_json(group_ids_config_path)
     inhibitory_neuron_density = compute_inhibitory_neuron_density(
         region_map,
         annotation.raw,
@@ -492,7 +505,7 @@ def inhibitory_and_excitatory_neuron_densities(
         VoxelData.load_nrrd(nrn1_path).raw,
         neuron_density.raw,
         inhibitory_data=inhibitory_data(inhibitory_df),
-        root_region_name=root_region_name,
+        group_ids_config=group_ids_config,
     )
 
     if not Path(output_dir).exists():
@@ -775,6 +788,13 @@ def measurements_to_average_densities(
     help="Path to the json file containing the fitting coefficients and standard deviations"
     "for each region group and each cell type.",
 )
+@click.option(
+    "--group-ids-config-path",
+    type=EXISTING_FILE_PATH,
+    default=utils.GROUP_IDS_PATH,
+    help="Path to density groups ids config",
+    show_default=True,
+)
 @log_args(L)
 def fit_average_densities(
     hierarchy_path,
@@ -786,6 +806,7 @@ def fit_average_densities(
     homogenous_regions_path,
     fitted_densities_output_path,
     fitting_maps_output_path,
+    group_ids_config_path,
 ):  # pylint: disable=too-many-arguments, too-many-locals
     """
     Estimate average cell densities of brain regions in `hierarchy_path` for the cell types
@@ -835,6 +856,7 @@ def fit_average_densities(
     - some regions can have NaN density values for one or more cell types because they are not
     covered by the selected slices of the volumetric gene marker intensities.
     """
+    Path(fitted_densities_output_path).parent.mkdir(parents=True, exist_ok=True)
 
     L.info("Loading annotation ...")
     annotation = VoxelData.load_nrrd(annotation_path)
@@ -857,17 +879,14 @@ def fit_average_densities(
     voxel_data += list(gene_voxeldata.values())
     assert_properties(voxel_data)
 
-    with open(config["realignedSlicesPath"], "r", encoding="utf-8") as file_:
-        slices = json.load(file_)
-
-    with open(config["cellDensityStandardDeviationsPath"], "r", encoding="utf-8") as file_:
-        cell_density_stddev = json.load(file_)
-        cell_density_stddev = {
-            # Use the AIBS name attribute as key (this is a unique identifier in 1.json)
-            # (Ex: former key "|root|Basic cell groups and regions|Cerebrum" -> new key: "Cerebrum")
-            name.split("|")[-1]: stddev
-            for (name, stddev) in cell_density_stddev.items()
-        }
+    slices = utils.load_json(config["realignedSlicesPath"])
+    cell_density_stddev = utils.load_json(config["cellDensityStandardDeviationsPath"])
+    cell_density_stddev = {
+        # Use the AIBS name attribute as key (this is a unique identifier in 1.json)
+        # (Ex: former key "|root|Basic cell groups and regions|Cerebrum" -> new key: "Cerebrum")
+        name.split("|")[-1]: stddev
+        for (name, stddev) in cell_density_stddev.items()
+    }
 
     gene_marker_volumes = {
         gene: {
@@ -877,6 +896,8 @@ def fit_average_densities(
         for (gene, gene_data) in gene_voxeldata.items()
     }
 
+    group_ids_config = utils.load_json(group_ids_config_path)
+
     L.info("Loading average densities dataframe ...")
     average_densities_df = pd.read_csv(average_densities_path)
     homogenous_regions_df = pd.read_csv(homogenous_regions_path)
@@ -891,6 +912,7 @@ def fit_average_densities(
         homogenous_regions_df,
         cell_density_stddev,
         region_name=region_name,
+        group_ids_config=group_ids_config,
     )
 
     # Turn index into column so as to ease off the save and load operations on csv files
@@ -985,14 +1007,15 @@ def inhibitory_neuron_densities(
     neuron_density = VoxelData.load_nrrd(neuron_density_path)
     if np.any(neuron_density.raw < 0.0):
         raise AtlasDensitiesError(f"Negative density value found in {neuron_density_path}.")
+
     L.info("Loading hierarchy ...")
-    with open(hierarchy_path, "r", encoding="utf-8") as file_:
-        hierarchy = json.load(file_)
-        if "msg" in hierarchy:
-            L.warning("Top-most object contains 'msg'; assuming AIBS JSON layout")
-            if len(hierarchy["msg"]) > 1:
-                raise AtlasDensitiesError("Unexpected JSON layout (more than one 'msg' child)")
-            hierarchy = hierarchy["msg"][0]
+
+    hierarchy = utils.load_json(hierarchy_path)
+    if "msg" in hierarchy:
+        L.warning("Top-most object contains 'msg'; assuming AIBS JSON layout")
+        if len(hierarchy["msg"]) > 1:
+            raise AtlasDensitiesError("Unexpected JSON layout (more than one 'msg' child)")
+        hierarchy = hierarchy["msg"][0]
 
     # Consistency check
     L.info("Checking consistency ...")

atlas_densities/app/combination.py

@@ -10,15 +10,21 @@
 import pandas as pd
 import voxcell  # type: ignore
 import yaml  # type: ignore
-from atlas_commons.app_utils import EXISTING_FILE_PATH, common_atlas_options, log_args, set_verbose
+from atlas_commons.app_utils import (
+    EXISTING_FILE_PATH,
+    common_atlas_options,
+    log_args,
+    set_verbose,
+    verbose_option,
+)
 
 from atlas_densities.combination import annotations_combinator, markers_combinator
 
 L = logging.getLogger(__name__)
 
 
 @click.group()
-@click.option("-v", "--verbose", count=True)
+@verbose_option
 def app(verbose):
     """Run the combination CLI"""
     set_verbose(L, verbose)

atlas_densities/app/data/metadata/README.txt

@@ -0,0 +1,47 @@
+The groups_ids.json has the following layout:
+    [
+      {
+        "name": "Cerebellum group",
+        "UNION": [
+            { "name": "Cerebellum", "with_descendants": true },
+            { "name": "arbor vitae", "with_descendants": true }
+        ]
+      },
+      {
+        "name": "Molecular layer",
+        "INTERSECT": [
+          "!Cerebellar group",
+          { "name": "@.*molecular layer", "with_descendants": true }
+        ]
+      },
+      {
+        "name": "Rest",
+        "REMOVE": [
+          { "name": "root", "with_descendants": true },
+          { "UNION": [ "!Cerebellum group", ] }
+        ]
+      }
+      [....]
+    ]
+
+The config is an ordered list of stanzas; each stanza is a dictionary with a "name" key, whose value is the Group Name.
+This is followed by a key with one of the following keywords, and a list of clauses:
+    * `UNION` keyword creates a union of the ids found by the list of clauses.
+    * `INTERSECT` keyword creates an intersection of the ids found by the list of clauses.
+    * `REMOVE` keyword removes the ids in the second clause from the those of the first.
+
+
+A clause is formed of dictionary of the form:
+    {"$attribute_name": "$attribute_value", "with_descendants": true}
+The attribute name is used for the atlas lookup, things like "name" or "acronym" are valid.
+
+Finally, one can refer to a previous stanza by preceding it with a `!`.
+
+The current `group_ids.json was discussed here:
+    https://github.com/BlueBrain/atlas-densities/pull/51#issuecomment-1748445813
+    In very short, in older versions of the annotations, fibers and main regions were in separated files (ccfv2 / ccfbbp).
+    In the main regions, the space allocated to the fibers were annotated to Basic cell groups and regions.
+    For some weird reasons, the fibers do not take the entire space allocated for them.
+    Which means that for most of the voxels annotated to Basic cell groups and regions, they correspond to fibers.
+    I say most of the voxels because there are also voxels at the frontier between two main regions (e.g. Cerebellum) that were not annotated to their closest leaf region.
+    These voxels are gone in ccfv3.

atlas_densities/app/data/metadata/ca1_metadata.json

@@ -13,3 +13,4 @@
         "with_descendants": true
     }
 }
+