diff --git a/.pylintrc b/.pylintrc index bcef07f..d46cfb3 100644 --- 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 index c466b98..9106552 100644 --- 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 index f64bfc8..4781296 100644 --- 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 new file mode 100644 index 0000000..83da595 --- /dev/null +++ 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 new file mode 100644 index 0000000..c83f4cf --- /dev/null +++ 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 new file mode 100644 index 0000000..4474cfe --- /dev/null +++ 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 diff --git a/tests/densities/test_excitatory_inhibitory_splitting.py b/tests/densities/test_excitatory_inhibitory_splitting.py new file mode 100644 index 0000000..b212da5 --- /dev/null +++ b/tests/densities/test_excitatory_inhibitory_splitting.py @@ -0,0 +1,417 @@ +from pathlib import Path +from unittest.mock import patch + +import numpy as np +import numpy.testing as npt +import pandas as pd +import pytest +from voxcell import RegionMap, VoxelData + +import atlas_densities.densities.excitatory_inhibitory_splitting as tested +from atlas_densities.app.utils import DATA_PATH + +TESTS_PATH = Path(__file__).parent.parent +HIERARCHY_PATH = Path(TESTS_PATH, "1.json") + + +def test_scale_excitatory_densities(tmp_path): + shape = ( + 3, + 3, + ) + brain_regions = VoxelData( + np.ones(shape), + voxel_dimensions=(10.0, 10.0), + ) + mapping = pd.DataFrame() + layer_ids = pd.DataFrame() + density = VoxelData( + np.ones(shape), + voxel_dimensions=(10.0, 10.0), + ) + tested.scale_excitatory_densities(tmp_path, brain_regions, mapping, layer_ids, density) + + # empty mapping + output = tmp_path.glob("*_cADpyr.nrrd") + assert len(list(output)) == 0 + + mapping = pd.DataFrame( + { + "L2_IPC": { + "layer_2": 0.1, + "layer_3": 0, + }, + "L2_TPC:A": { + "layer_2": 0.2, + "layer_3": 0, + }, + "L2_TPC:B": { + "layer_2": 0.8, + "layer_3": 0, + }, + "L3_TPC:A": { + "layer_2": 0, + "layer_3": 0.7, + }, + } + ) + mapping.index.name = "layer" + layer_ids = {"layer_2": [], "layer_3": []} + tested.scale_excitatory_densities(tmp_path, brain_regions, mapping, layer_ids, density) + output = list(tmp_path.glob("*_cADpyr.nrrd")) + + filenames = {n.name: n for n in output} + assert set(filenames) == { + "L3_TPC:A_cADpyr.nrrd", + "L2_IPC_cADpyr.nrrd", + "L2_TPC:A_cADpyr.nrrd", + "L2_TPC:B_cADpyr.nrrd", + } + # no ids specified, all densities should be 0 + assert all(VoxelData.load_nrrd(v).raw.sum() == 0.0 for v in filenames.values()) + + layer_ids = {"layer_2": [], "layer_3": [1]} + tested.scale_excitatory_densities(tmp_path, brain_regions, mapping, layer_ids, density) + output = list(tmp_path.glob("*_cADpyr.nrrd")) + values = {k: VoxelData.load_nrrd(v).raw.sum() for k, v in filenames.items()} + assert values["L2_IPC_cADpyr.nrrd"] == 0.0 + assert values["L2_TPC:B_cADpyr.nrrd"] == 0.0 + assert values["L2_TPC:A_cADpyr.nrrd"] == 0.0 + assert values["L3_TPC:A_cADpyr.nrrd"] == pytest.approx(0.7 * 3 * 3) + + +def test_set_ids_to_zero_and_save(tmp_path): + shape = ( + 3, + 3, + ) + raw = np.ones(shape) + path = tmp_path / "zero.nrrd" + brain_regions = VoxelData( + raw.copy(), + voxel_dimensions=(10.0, 10.0), + ) + nrrd = VoxelData( + raw.copy(), + voxel_dimensions=(10.0, 10.0), + ) + tested.set_ids_to_zero_and_save( + path, + brain_regions, + nrrd, + [ + 1, + ], + ) + + output = VoxelData.load_nrrd(path) + assert output.shape == nrrd.shape + assert np.sum(output.raw) == 0.0 + + raw[0, 0] = 10 + raw[0, 1] = 11 + brain_regions = VoxelData(raw, voxel_dimensions=(10.0, 10.0)) + tested.set_ids_to_zero_and_save(path, brain_regions, nrrd, [1, 10]) + output = VoxelData.load_nrrd(path) + + assert np.sum(output.raw) == 1.0 + + +def test_make_excitatory_density(): + shape = ( + 3, + 3, + ) + neuron_density = VoxelData( + 2 * np.ones(shape), + voxel_dimensions=(10.0, 10.0), + ) + inhibitory_density = VoxelData( + np.ones(shape), + voxel_dimensions=(10.0, 10.0), + ) + res = tested.make_excitatory_density(neuron_density, inhibitory_density) + + assert res.shape == neuron_density.shape + assert np.sum(res.raw) == np.product(neuron_density.shape) + + # this would create negative densities; make sure they are clipped to zero + res = tested.make_excitatory_density(inhibitory_density, neuron_density) + assert res.shape == neuron_density.shape + assert np.sum(res.raw) == 0.0 + + +def test_gather_isocortex_ids_from_metadata(): + region_map = RegionMap.load_json(HIERARCHY_PATH) + metadata_path = DATA_PATH / "metadata" / "excitatory-inhibitory-splitting.json" + res = tested.gather_isocortex_ids_from_metadata(region_map, metadata_path) + acronyms = { + "layer_1": [ + "ACA1", + "ACAd1", + "ACAv1", + "AId1", + "AIp1", + "AIv1", + "AUDd1", + "AUDp1", + "AUDpo1", + "AUDv1", + "ECT1", + "FRP1", + "GU1", + "ILA1", + "MO1", + "MOp1", + "MOs1", + "ORB1", + "ORBl1", + "ORBm1", + "ORBvl1", + "PERI1", + "PL1", + "PTLp1", + "RSPagl1", + "RSPd1", + "RSPv1", + "SS1", + "SSp-bfd1", + "SSp-ll1", + "SSp-m1", + "SSp-n1", + "SSp-tr1", + "SSp-ul1", + "SSp-un1", + "SSp1", + "SSs1", + "TEa1", + "VIS1", + "VISC1", + "VISa1", + "VISal1", + "VISam1", + "VISl1", + "VISli1", + "VISlla1", + "VISm1", + "VISmma1", + "VISmmp1", + "VISp1", + "VISpl1", + "VISpm1", + "VISpor1", + "VISrl1", + "VISrll1", + ], + "layer_2": ["PL2", "ORBm2", "RSPv2", "ILA2"], + "layer_3": [], + "layer_4": [ + "AUDd4", + "AUDp4", + "AUDpo4", + "AUDv4", + "GU4", + "PTLp4", + "RSPd4", + "SS4", + "SSp-bfd4", + "SSp-ll4", + "SSp-m4", + "SSp-n4", + "SSp-tr4", + "SSp-ul4", + "SSp-un4", + "SSp4", + "SSs4", + "TEa4", + "VIS4", + "VISC4", + "VISa4", + "VISal4", + "VISam4", + "VISl4", + "VISli4", + "VISlla4", + "VISm4", + "VISmma4", + "VISmmp4", + "VISp4", + "VISpl4", + "VISpm4", + "VISpor4", + "VISrl4", + "VISrll4", + ], + "layer_5": [ + "ACA5", + "ACAd5", + "ACAv5", + "AId5", + "AIp5", + "AIv5", + "AUDd5", + "AUDp5", + "AUDpo5", + "AUDv5", + "ECT5", + "FRP5", + "GU5", + "ILA5", + "MO5", + "MOp5", + "MOs5", + "ORB5", + "ORBl5", + "ORBm5", + "ORBvl5", + "PERI5", + "PL5", + "PTLp5", + "RSPagl5", + "RSPd5", + "RSPv5", + "SS5", + "SSp-bfd5", + "SSp-ll5", + "SSp-m5", + "SSp-n5", + "SSp-tr5", + "SSp-ul5", + "SSp-un5", + "SSp5", + "SSs5", + "TEa5", + "VIS5", + "VISC5", + "VISa5", + "VISal5", + "VISam5", + "VISl5", + "VISli5", + "VISlla5", + "VISm5", + "VISmma5", + "VISmmp5", + "VISp5", + "VISpl5", + "VISpm5", + "VISpor5", + "VISrl5", + "VISrll5", + ], + "layer_6": [ + "ACA6a", + "ACA6b", + "ACAd6a", + "ACAd6b", + "ACAv6a", + "ACAv6b", + "AId6a", + "AId6b", + "AIp6a", + "AIp6b", + "AIv6a", + "AIv6b", + "AUDd6a", + "AUDd6b", + "AUDp6a", + "AUDp6b", + "AUDpo6a", + "AUDpo6b", + "AUDv6a", + "AUDv6b", + "ECT6a", + "ECT6b", + "FRP6a", + "FRP6b", + "GU6a", + "GU6b", + "ILA6a", + "ILA6b", + "MO6a", + "MO6b", + "MOp6a", + "MOp6b", + "MOs6a", + "MOs6b", + "ORB6a", + "ORB6b", + "ORBl6a", + "ORBl6b", + "ORBm6a", + "ORBm6b", + "ORBvl6a", + "ORBvl6b", + "PERI6a", + "PERI6b", + "PL6a", + "PL6b", + "PTLp6a", + "PTLp6b", + "RSPagl6a", + "RSPagl6b", + "RSPd6a", + "RSPd6b", + "RSPv6a", + "RSPv6b", + "SS6a", + "SS6b", + "SSp-bfd6a", + "SSp-bfd6b", + "SSp-ll6a", + "SSp-ll6b", + "SSp-m6a", + "SSp-m6b", + "SSp-n6a", + "SSp-n6b", + "SSp-tr6a", + "SSp-tr6b", + "SSp-ul6a", + "SSp-ul6b", + "SSp-un6a", + "SSp-un6b", + "SSp6a", + "SSp6b", + "SSs6a", + "SSs6b", + "TEa6a", + "TEa6b", + "VIS6a", + "VIS6b", + "VISC6a", + "VISC6b", + "VISa6a", + "VISa6b", + "VISal6a", + "VISal6b", + "VISam6a", + "VISam6b", + "VISl6a", + "VISl6b", + "VISli6a", + "VISli6b", + "VISlla6a", + "VISlla6b", + "VISm6a", + "VISm6b", + "VISmma6a", + "VISmma6b", + "VISmmp6a", + "VISmmp6b", + "VISp6a", + "VISp6b", + "VISpl6a", + "VISpl6b", + "VISpm6a", + "VISpm6b", + "VISpor6a", + "VISpor6b", + "VISrl6a", + "VISrl6b", + "VISrll6a", + "VISrll6b", + ], + } + + for layer, ids in res.items(): + assert set(region_map.get(i, "acronym") for i in ids) == set(acronyms[layer])