Excitatory split (#19)

* Write SSCX layer-specific types to all isocortex Changed layer key in csv file, changed layer filter in metadata to make it specific for Isocortex layers, changed code to do isocortex instead of just SSCX * made it isocortex specific This corrects the issue with hippocampus being found by adding an extra filter to make sure returned layers are in isocortex. Still needs to be able to use the native regex and not the accessory one provided, though. There appears to be a function for get_layer_mask that would probably be more elegant to use. Used code snippet from Dimitri to read in metadata and create layers. It seems to filter okay for isocortex. I foresee the metadata.json file being put into the CLL, so for now it's included as a helper file but eventually would be specified. Co-authored-by: dkeller9 <[email protected]>
BlueBrain · Nov 17, 2022 · 6805016 · 6805016
1 parent eca7723
commit 6805016
Show file tree

Hide file tree

Showing 7 changed files with 682 additions and 3 deletions.
diff --git a/.pylintrc b/.pylintrc
@@ -1,5 +1,5 @@
 [MESSAGES CONTROL]
-disable=
+disable=invalid-name
 
 [SIMILARITIES]
 # Minimum lines number of a similarity.

diff --git a/README.rst b/README.rst
@@ -317,8 +317,8 @@ transfer functions computed previously (first density estimates).
         --average-densities-path=data/ccfv2/first_estimates/first_estimates.csv \
         --output-dir=data/ccfv2/density_volumes/
 
-Compute ME-types densities from a probality map
------------------------------------------------
+Compute ME-types densities from a probability map
+-------------------------------------------------
 
 Morphological and Electrical type densities of inhibitory neurons in the isocortex can be estimated
 using Roussel et al.'s pipeline. This pipeline produces a mapping from inhibitory neuron molecular
@@ -337,6 +337,30 @@ mapping csv file (see also `mtypes_probability_map_config.yaml`_).
         --mtypes-config-path=$DATA/mtypes/mtypes_probability_map_config.yaml \
         --output-dir=data/ccfv2/me-types/
 
+Subdivide excitatory files into pyramidal subtypes
+--------------------------------------------------
+
+This should run after the ME mapping from Roussel. To run:
+
+make_new_inhib_exc_outside_SSCX.py
+
+Input support files needed:
+
+annotation atlas (with L23 split): brain_regions.nrrd
+
+annotation hierarchy (with L23 split): hierarchy.json
+
+inhibitory: gad67+_density_v3_MMB.nrrd
+
+total neuron: neuron_density_v3_MMB.nrrd
+
+Output:
+
+excitatory nrrd files in output directory
+
+inhibitory minus cortex in output directory
+
+
 Instructions for developers
 ===========================
 

diff --git a/atlas_densities/app/cell_densities.py b/atlas_densities/app/cell_densities.py
@@ -58,6 +58,7 @@
 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.cell_counts import (
     extract_inhibitory_neurons_dataframe,
     glia_cell_counts,
@@ -84,6 +85,10 @@
 from atlas_densities.densities.utils import zero_negative_values
 from atlas_densities.exceptions import AtlasDensitiesError
 
+EXCITATORY_SPLIT_CORTEX_ALL_TO_EXC_MTYPES = (
+    DATA_PATH / "mtypes" / "mapping_cortex_all_to_exc_mtypes.csv"
+)
+EXCITATORY_SPLIT_METADATA = DATA_PATH / "metadata" / "excitatory-inhibitory-splitting.json"
 HOMOGENOUS_REGIONS_PATH = DATA_PATH / "measurements" / "homogenous_regions.csv"
 HOMOGENOUS_REGIONS_REL_PATH = HOMOGENOUS_REGIONS_PATH.relative_to(AD_PATH)
 MARKERS_README_REL_PATH = (DATA_PATH / "markers" / "README.rst").relative_to(AD_PATH)
@@ -951,3 +956,97 @@ def inhibitory_neuron_densities(
         annotation.with_data(volumetric_density).save_nrrd(
             str(Path(output_dir, f"{cell_type}_density.nrrd"))
         )
+
+
+@app.command()
+@click.option(
+    "--annotation-path",
+    type=EXISTING_FILE_PATH,
+    required=True,
+    help="The path to the whole mouse brain annotation file (nrrd).",
+)
+@click.option(
+    "--hierarchy-path",
+    type=EXISTING_FILE_PATH,
+    required=True,
+    help="The path to the hierarchy file, i.e., AIBS 1.json or BBP hierarchy.json.",
+)
+@click.option(
+    "--neuron-density",
+    type=EXISTING_FILE_PATH,
+    required=True,
+    help="Complete neuron density for full brain",
+)
+@click.option(
+    "--inhibitory-density",
+    type=EXISTING_FILE_PATH,
+    required=True,
+    help="Complete inhibitory density for full brain",
+)
+@click.option(
+    "--cortex-all-to-exc-mtypes",
+    type=EXISTING_FILE_PATH,
+    required=True,
+    default=EXCITATORY_SPLIT_CORTEX_ALL_TO_EXC_MTYPES,
+    help="CSV file with mappings for isocortex mtypes",
+    show_default=True,
+)
+@click.option(
+    "--metadata-path",
+    type=EXISTING_FILE_PATH,
+    required=True,
+    default=EXCITATORY_SPLIT_METADATA,
+    help="CSV file with mappings for isocortex mtypes",
+    show_default=True,
+)
+@click.option("--output-dir", required=True, help="Output path")
+@log_args(L)
+def excitatory_split(
+    annotation_path,
+    hierarchy_path,
+    neuron_density,
+    inhibitory_density,
+    cortex_all_to_exc_mtypes,
+    metadata_path,
+    output_dir,
+):
+    """
+    This program makes exc and inh densities with isocortex cut out
+    It also remaps exc cells to morphological fractions in those regions
+    using the m-type fractions in the csv file. All etype exc types are the same
+    """
+    output_dir = Path(output_dir)
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    region_map = RegionMap.load_json(hierarchy_path)
+    brain_regions = VoxelData.load_nrrd(annotation_path)
+
+    inhibitory_density = VoxelData.load_nrrd(inhibitory_density)
+    excitatory_density = excitatory_inhibitory_splitting.make_excitatory_density(
+        VoxelData.load_nrrd(neuron_density), inhibitory_density
+    )
+
+    layer_ids = excitatory_inhibitory_splitting.gather_isocortex_ids_from_metadata(
+        region_map, metadata_path
+    )
+
+    excitatory_mapping = pd.read_csv(cortex_all_to_exc_mtypes).set_index("layer")
+
+    excitatory_inhibitory_splitting.scale_excitatory_densities(
+        output_dir, brain_regions, excitatory_mapping, layer_ids, excitatory_density
+    )
+
+    remove_ids = sum(layer_ids.values(), [])
+
+    excitatory_inhibitory_splitting.set_ids_to_zero_and_save(
+        str(output_dir / "Generic_Excitatory_Neuron_MType_Generic_Excitatory_Neuron_EType.nrrd"),
+        brain_regions,
+        excitatory_density,
+        remove_ids,
+    )
+    excitatory_inhibitory_splitting.set_ids_to_zero_and_save(
+        str(output_dir / "Generic_Inhibitory_Neuron_MType_Generic_Inhibitory_Neuron_EType.nrrd"),
+        brain_regions,
+        inhibitory_density,
+        remove_ids,
+    )
diff --git a/atlas_densities/app/data/metadata/excitatory-inhibitory-splitting.json b/atlas_densities/app/data/metadata/excitatory-inhibitory-splitting.json
@@ -0,0 +1,15 @@
+{
+    "region": {
+        "name": "Isocortex",
+        "query": "Isocortex",
+        "attribute": "acronym",
+        "with_descendants": true
+
+    },
+    "layers": {
+        "names": ["layer_1", "layer_2", "layer_3", "layer_4", "layer_5", "layer_6"],
+        "queries": ["@.*1[ab]?$", "@.*2[ab]?$", "@.*[^/]3[ab]?$", "@.*4[ab]?$", "@.*5[ab]?$", "@.*6[ab]?$"],
+        "attribute": "name",
+        "with_descendants": true
+    }
+}
diff --git a/atlas_densities/app/data/mtypes/mapping_cortex_all_to_exc_mtypes.csv b/atlas_densities/app/data/mtypes/mapping_cortex_all_to_exc_mtypes.csv
@@ -0,0 +1,6 @@
+layer,L2_IPC,L2_TPC:A,L2_TPC:B,L3_TPC:A,L3_TPC:B,L4_SSC,L4_TPC,L4_UPC,L5_TPC:A,L5_TPC:B,L5_TPC:C,L5_UPC,L6_BPC,L6_HPC,L6_IPC,L6_TPC:A,L6_TPC:C,L6_UPC
+layer_2,0.082612842,0.129525291,0.787861866,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+layer_3,0,0,0,0.776524986,0.223475014,0,0,0,0,0,0,0,0,0,0,0,0,0
+layer_4,0,0,0,0,0,0.090735874,0.641289425,0.267974701,0,0,0,0,0,0,0,0,0,0
+layer_5,0,0,0,0,0,0,0,0,0.478764389,0.373704348,0.069744099,0.077787164,0,0,0,0,0,0
+layer_6,0,0,0,0,0,0,0,0,0,0,0,0,0.181041554,0.165265851,0.226502638,0.194968725,0.129358437,0.102862796
diff --git a/atlas_densities/densities/excitatory_inhibitory_splitting.py b/atlas_densities/densities/excitatory_inhibitory_splitting.py
@@ -0,0 +1,118 @@
+""" This program makes exc and inh densities with SSCX cut out
+
+It also remaps exc cells to morphological fractions in those regions
+using the m-type fractions in the cxv file. All etype exc types are the same.
+"""
+import json
+import logging
+
+import numpy as np
+
+L = logging.getLogger(__name__)
+
+
+def scale_excitatory_densities(output, brain_regions, mapping, layer_ids, density):
+    """Scale density by `mapping`
+
+    Args:
+        output (Path): path to output directory
+        brain_regions (VoxelData): annotations of brain regions
+        mapping (DataFrame): mapping of layers and mtypes to scaling factor ex:
+                       L2_IPC  L2_TPC:A ...
+            layer
+            layer_2  0.082613  0.129525 ...
+            layer_3  0.000000  0.000000 ...
+            layer_4  0.000000  0.000000 ...
+            layer_5  0.000000  0.000000 ...
+            layer_6  0.000000  0.000000 ...
+
+        layer_ids (dict layer_name -> list of ids): ids to scale
+        density (VoxelData): initial density to scale
+    """
+    assert (mapping.to_numpy() >= 0).all(), "Cannot have negative scaling factors"
+
+    for layer, df in mapping.iterrows():
+        L.info("Performing layer: %s", layer)
+        idx = np.nonzero(np.isin(brain_regions.raw, layer_ids[layer]))
+
+        for mtype, scale in df.iteritems():
+            if scale == 0:
+                continue
+
+            L.info("  Performing %s", mtype)
+
+            raw = np.zeros_like(density.raw)
+            raw[idx] = scale * density.raw[idx]
+
+            output_name = output / f"{mtype}_cADpyr.nrrd"
+            density.with_data(raw).save_nrrd(str(output_name))
+
+
+def set_ids_to_zero_and_save(output_name, brain_regions, nrrd, remove_ids):
+    """Take an VoxelData, set some of it's values to 0, and save it
+
+    Args:
+        output_name (str): output filename
+        brain_regions (VoxelData): annotations of brain regions
+        nrrd (VoxelData): data where values should be set to zero
+        remove_ids (list of ids): ids where the value should be set to zero
+    """
+    L.info("Setting region ids to zero, outputting %s", output_name)
+
+    idx = np.nonzero(np.isin(brain_regions.raw, remove_ids))
+    raw = nrrd.raw.copy()
+    raw[idx] = 0.0
+    nrrd.with_data(raw).save_nrrd(str(output_name))
+
+
+def make_excitatory_density(neuron_density, inhibitory_density):
+    """Given neuron and inhibitory cell density, calculate the excitatory density
+
+    Args:
+        neuron_density (VoxelData): total neuron density
+        inhibitory_density (VoxelData): total inhibitory density
+
+    Returns:
+        VoxelData: excitatory density
+    """
+    exc_density = neuron_density.with_data(
+        np.clip(neuron_density.raw - inhibitory_density.raw, a_min=0, a_max=None)
+    )
+    return exc_density
+
+
+def gather_isocortex_ids_from_metadata(region_map, metadata_path):
+    """given metadata, return ids of region interest
+
+    Args:
+        region_map (voxcell.RegionMap): RegionMap object to use
+        metadata_path (str): path to json with metadata
+
+    Returns:
+        dict with layer -> list of ids
+    """
+    with open(metadata_path, encoding="utf-8") as fd:
+        metadata = json.load(fd)
+
+    # set of ids of Isocortex
+    region_ids = region_map.find(
+        metadata["region"]["query"],
+        attr=metadata["region"]["attribute"],
+        with_descendants=metadata["region"].get("with_descendants", False),
+    )
+
+    # set of ids of one layer
+    metadata_layers = metadata["layers"]
+    layer_ids = {
+        f"layer_{i}": list(
+            region_ids
+            & region_map.find(
+                query,
+                attr=metadata_layers["attribute"],
+                with_descendants=metadata_layers.get("with_descendants", False),
+            )
+        )
+        for i, query in enumerate(metadata_layers["queries"], 1)
+    }
+
+    return layer_ids