diff --git a/atlas_densities/app/cell_densities.py b/atlas_densities/app/cell_densities.py index d0a618e..f64bfc8 100644 --- a/atlas_densities/app/cell_densities.py +++ b/atlas_densities/app/cell_densities.py @@ -170,7 +170,7 @@ def cell_density(annotation_path, hierarchy_path, nissl_path, output_path): _get_voxel_volume_in_mm3(annotation), nissl.raw, ) - nissl.with_data(np.asarray(overall_cell_density, dtype=float)).save_nrrd(output_path) + nissl.with_data(overall_cell_density).save_nrrd(output_path) @app.command() diff --git a/atlas_densities/app/combination.py b/atlas_densities/app/combination.py index 31d6c53..2709462 100644 --- a/atlas_densities/app/combination.py +++ b/atlas_densities/app/combination.py @@ -52,7 +52,7 @@ def app(verbose): ) @click.option("--output-path", required=True, help="Path of the nrrd file to write") @log_args(L) -def combine_annotations( +def combine_v2_v3_annotations( hierarchy, brain_annotation_ccfv2, fiber_annotation_ccfv2, @@ -60,7 +60,7 @@ def combine_annotations( output_path, ): # pylint: disable=line-too-long - """Generate and save the combined annotation file. + """Generate and save the v2-v3 combined annotation file. The annotation file `brain_annotation_ccfv3` is the annotation file containing the least complete annotation. There are two use cases: a resolution of 10 or 25 um. @@ -102,6 +102,56 @@ def combine_annotations( combined_annotation.save_nrrd(output_path) +@app.command() +@click.option( + "--brain-annotation-ccfv2", + type=EXISTING_FILE_PATH, + required=True, + help=("This brain annotation file contains the most complete annotation."), +) +@click.option( + "--fiber-annotation-ccfv2", + type=EXISTING_FILE_PATH, + required=True, + help="Fiber annotation is not included in the CCF-v2 2011 annotation files.", +) +@click.option("--output-path", required=True, help="Path of the nrrd file to write") +@log_args(L) +def combine_ccfv2_annotations( + brain_annotation_ccfv2, + fiber_annotation_ccfv2, + output_path, +): + # pylint: disable=line-too-long + """Generate and save the ccfv2 combined annotation file. + + The ccfv2 annotations are split into two volumes. `fiber_annotation_ccfv2` describes solely the + fibers and ventricular related regions while `brain_annotation_ccfv2` contains all other brain + regions. There are two use cases: a resolution of 10 or 25 um. + + For a resolution of 10 um, the path arguments should be the following: + + \b + - `brain_annotation_ccfv2` is the path to a copy of + ``AIBS_ANNOTATION_URL``/mouse_2011/annotation_10.nrrd + + \b + - `fiber_annotation` is the path to a copy of + ``AIBS_ANNOTATION_URL``/mouse_2011/annotationFiber_10_2011.nrrd + + where ``AIBS_ANNOTATION_URL`` is + http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/annotation. + """ + + brain_annotation_ccfv2 = voxcell.VoxelData.load_nrrd(brain_annotation_ccfv2) + fiber_annotation_ccfv2 = voxcell.VoxelData.load_nrrd(fiber_annotation_ccfv2) + + combined_annotation = annotations_combinator.combine_ccfv2_annotations( + brain_annotation_ccfv2, fiber_annotation_ccfv2 + ) + combined_annotation.save_nrrd(output_path) + + @app.command() @common_atlas_options @click.option( diff --git a/atlas_densities/app/data/markers/README.rst b/atlas_densities/app/data/markers/README.rst index 095fcff..700b719 100644 --- a/atlas_densities/app/data/markers/README.rst +++ b/atlas_densities/app/data/markers/README.rst @@ -4,38 +4,56 @@ Description This folder contains data pertaining to gene marker volumes provided by AIBS. -The file `realigned_slices.json` contains the slice indices of several gene marker volumetric files which have been -realigned manually by Csaba Eroe (former BBP PHD student). It has the form of a dictionary whose keys are AIBS dataset identifiers -(see the AIBS find queries of `Gene Search`_) and whose values are lists of integers. +The file `realigned_slices_XXX.json` contains position indices of the original image slices for +several AIBS ISH experiments. It has the form of a dictionary whose keys are AIBS dataset +identifiers (see the AIBS find queries of `Gene Search`_) and whose values are the position indices +of the slice images with respect to the reference volumes. These position values can be extracted +from the ISH experiments metadata provided by the AIBS. -The computation of average marker intensities in AIBS brain regions is restricted to these slices which can be trusted. -Such a computation is the basis of a linear fitting -on a 2D point cloud (average marker intensity, average cell density) used to estimate cell densities in regions where no -measurement is available. +The content of `realigned_slices_XXX.json` depends on which annotation and Nissl nrrd files are +considered (`_XXX` in the filename). Two versions of this file are provided in this folder, one for +the AIBS CCFv2 version and one for the CCFbbp version (see `Rodarie et al. (2021)`_ for more +details). The ISH images from the AIBS are registered to the CCFv2 and CCFv3 brain volume so no +transformation is required for these reference volume versions. For the CCFbbp, however, some ISH +slices position are slightly shifted with respect to the Nisslbbp. Hence, a manual selection of the +matching Nissl slice for each ISH slice is required to obtain the best cell density estimates. -The file `realigned_slices.json` is referred to by the gene marker configuration file, -an input of the command line `atlas-densities cell-densities fit-average-densities`, -see `atlas_densities.app.cell_densities`. The gene marker configuration binds -gene marker volumetric files (nrrd files) with the appropriate slices of `realigned_slices.json`. - -The content of `realigned_slices.json` depends on which annotation and Nissl nrrd files are considered. -For atlas files older than AIBS CCFv2, one should use the file stored in this directory when computing average -densities. This file is the companion file of a manual re-alignment process conducted by Csaba Eroe (former BBP PHD student). +The computation of average marker intensities in AIBS brain regions is restricted to the original +ISH slices, which can be trusted (see `Rodarie et al. (2021)`_). Such a computation is the basis of +a linear fitting on a 2D point cloud (average marker intensity, average cell density) used to +estimate inhibitory neuron densities in regions where no literature measurement is available. -For vanilla AIBS CCFV2 or CCFV3 files, `realigned_slices.json` is filled with the indices inferred the AIBS 2D images metadata -(no BBP processing in this case). +The file `realigned_slices_XXX.json` is referred to by the gene marker configuration file, +an input of the command line `atlas-densities cell-densities fit-average-densities`, +see `atlas_densities.app.cell_densities`. The gene marker configuration binds gene marker volumetric +files (nrrd files) with the appropriate slices of `realigned_slices_XXX.json`. -The file `std_cells.json` contains the standard deviation of every AIBS region volumetric cell -density. This file has been produced by Csabe Eroe, assuming that Nissl stains intensity depends linearily on volumetric cell density across the mouse brain. This file is used -to assign a "standard deviation" value to every region where a pre-computed volumetric neuron density -has been used to assign an average neuron density. This holds for every region whose neurons are inhibitory -only. (Those regions are listed in `atlas_densities/app/data/measurements/homogenous_regions.csv`.) +The file `std_cells.json` (in `atlas_densities/app/data/measurements/`) contains the cell density +standard deviation the BBP cell atlas for every AIBS region. This file has been produced by +Csabe Eroe (see `Eroe et al. (2018)`_), based on multiple Nissl volumes, assuming that Nissl stains +intensity depends linearly on volumetric cell density across the mouse brain. This file is used to +assign a "standard deviation" value to every region where a pre-computed volumetric neuron density +has been used to assign an average neuron density. This holds for every region whose neurons are +inhibitory only. (Those regions are listed in +`atlas_densities/app/data/measurements/homogenous_regions.csv`.) Gene markers configuration files -------------------------------- -The option `--gene-config-path` of the CLI `atlas-densities cell-densities fit-average-densities` expects -a path to a yaml file of the following form: +Combine ISH datasets for glia cells +>>>>>>>>>> + +The option `--config` of the CLI `atlas-densities combination combine-markers` expects a path to a +yaml file describing the location of the ISH datasets files as well weight parameters to combine +different genes expressed by glia cell types. `combine_markers_ccfv2_config.yaml` corresponds to +this configuration file, used to obtain the results of `Eroe et al. (2018)`_, using the CCFv2 +reference volumes. + +Fitting of transfer functions from mean region intensity to neuron density +>>>>>>>>>> + +The option `--gene-config-path` of the CLI `atlas-densities cell-densities fit-average-densities` +expects a path to a yaml file of the following form: .. code:: yaml @@ -49,10 +67,13 @@ a path to a yaml file of the following form: sst: 1001 vip: 77371835 gad67: 479 - realignedSlicesPath: "realigned_slices.json" + realignedSlicesPath: "realigned_slices_XXX.json" cellDensityStandardDeviationsPath: "std_cells.json" -The sectionDataSetID values are AIBS dataset identifiers recorded in `realigned_slices.json`. - +The sectionDataSetID values are AIBS dataset identifiers recorded in `realigned_slices_XXX.json`. +An example of this configuration file (`fit_average_densities_ccfv2_config.yaml`) is provided for the +CCFv2 reference volumes. -.. _`Gene Search`: https://mouse.brain-map.org/ \ No newline at end of file +.. _`Gene Search`: https://mouse.brain-map.org/ +.. _`Rodarie et al. (2021)`: https://www.biorxiv.org/content/10.1101/2021.11.20.469384v2 +.. _`Eroe et al. (2018)`: https://www.frontiersin.org/articles/10.3389/fninf.2018.00084/full \ No newline at end of file diff --git a/atlas_densities/app/data/markers/combine_markers_ccfv2_config.yaml b/atlas_densities/app/data/markers/combine_markers_ccfv2_config.yaml new file mode 100644 index 0000000..02a56c8 --- /dev/null +++ b/atlas_densities/app/data/markers/combine_markers_ccfv2_config.yaml @@ -0,0 +1,61 @@ +cellType: +- oligodendrocyte +- astrocyte +- microglia +brainRegion: +- cerebellum +- striatum +gene: +- cnp +- mbp +- gfap +- s100b +- aldh1l1 +- tmem119 +inputGeneVolumePath: + cnp: data/ccfv2/marker_volumes/CNP.nrrd + mbp: data/ccfv2/marker_volumes/MBP.nrrd + gfap: data/ccfv2/marker_volumes/GFAP.nrrd + s100b: data/ccfv2/marker_volumes/S100b.nrrd + aldh1l1: data/ccfv2/marker_volumes/ALDH1L1.nrrd + tmem119: data/ccfv2/marker_volumes/TMEM119.nrrd +outputCellTypeVolumePath: + oligodendrocyte: data/ccfv2/relative_marker_volumes/oligodendrocyte.nrrd + astrocyte: data/ccfv2/relative_marker_volumes/astrocyte.nrrd + microglia: data/ccfv2/relative_marker_volumes/microglia.nrrd +# Overall glia intensity volume +outputOverallGliaVolumePath: data/ccfv2/relative_marker_volumes/glia.nrrd +# The so-called global celltype scaling factors S_celltype of +# "A Cell Atlas for the Mouse Brain", C. Eroe et al, 2018. +# These factors are simply the proportions of the different glia cell types +# in the mouse brain. +outputCellTypeProportionsPath: data/ccfv2/relative_marker_volumes/glia_proportions.json +# Cell densities in number of cells per mm^3 +cellDensity: + cerebellum: + oligodendrocyte: 13750 + astrocyte: 1512 + microglia: 8624 + striatum: + oligodendrocyte: 9950 + astrocyte: 9867 + microglia: 12100 +combination: +- cellType: oligodendrocyte + gene: cnp + averageExpressionIntensity: 35.962800 +- cellType: oligodendrocyte + gene: mbp + averageExpressionIntensity: 3.304965 +- cellType: astrocyte + gene: gfap + averageExpressionIntensity: 3.2097903 +- cellType: astrocyte + gene: s100b + averageExpressionIntensity: 4.174825 +- cellType: astrocyte + gene: aldh1l1 + averageExpressionIntensity: 1.326080 +- cellType: microglia + gene: tmem119 + averageExpressionIntensity: 0.654761 \ No newline at end of file diff --git a/atlas_densities/app/data/markers/fit_average_densities_ccfv2_config.yaml b/atlas_densities/app/data/markers/fit_average_densities_ccfv2_config.yaml new file mode 100644 index 0000000..7a9d9cd --- /dev/null +++ b/atlas_densities/app/data/markers/fit_average_densities_ccfv2_config.yaml @@ -0,0 +1,12 @@ +inputGeneVolumePath: + pv: data/ccfv2/marker_volumes/pvalb.nrrd + sst: data/ccfv2/marker_volumes/SST.nrrd + vip: data/ccfv2/marker_volumes/VIP.nrrd + gad67: data/ccfv2/marker_volumes/gad1.nrrd +sectionDataSetID: + pv: "868" + sst: "1001" + vip: "77371835" + gad67: "479" +realignedSlicesPath: atlas_densities/app/measurements/markers/realigned_slices_ccfv2.json +cellDensityStandardDeviationsPath: atlas_densities/app/measurements/std_cells.json diff --git a/atlas_densities/app/data/markers/realigned_slices.json b/atlas_densities/app/data/markers/realigned_slices_bbp.json similarity index 100% rename from atlas_densities/app/data/markers/realigned_slices.json rename to atlas_densities/app/data/markers/realigned_slices_bbp.json diff --git a/atlas_densities/app/data/markers/realigned_slices_cccfv2.json b/atlas_densities/app/data/markers/realigned_slices_cccfv2.json new file mode 100644 index 0000000..799ebaa --- /dev/null +++ b/atlas_densities/app/data/markers/realigned_slices_cccfv2.json @@ -0,0 +1,587 @@ +{ + "868": [ + 49, + 57, + 74, + 82, + 91, + 100, + 109, + 117, + 126, + 134, + 144, + 152, + 160, + 169, + 177, + 186, + 195, + 203, + 211, + 220, + 228, + 236, + 245, + 254, + 262, + 270, + 279, + 287, + 295, + 304, + 312, + 321, + 329, + 338, + 346, + 354, + 362, + 371, + 379, + 387, + 395, + 403, + 411, + 419, + 428, + 436, + 445, + 453, + 461, + 470, + 478, + 486, + 494, + 502, + 511 + ], + "1001": [ + 64, + 72, + 80, + 88, + 96, + 104, + 112, + 120, + 128, + 136, + 144, + 152, + 160, + 168, + 176, + 184, + 192, + 200, + 208, + 216, + 224, + 232, + 240, + 248, + 256, + 264, + 272, + 280, + 288, + 296, + 303, + 311, + 319, + 328, + 335, + 343, + 351, + 359, + 367, + 375, + 383, + 391 + ], + "77371835": [ + 56, + 64, + 73, + 80, + 97, + 105, + 113, + 121, + 129, + 145, + 153, + 161, + 167, + 178, + 186, + 194, + 202, + 210, + 218, + 226, + 235, + 259, + 267, + 275, + 283, + 292, + 300, + 308, + 316, + 324, + 341, + 349, + 357, + 365, + 382, + 390, + 398, + 407, + 415, + 423, + 431, + 439, + 447, + 456, + 464, + 472, + 480, + 489, + 497, + 505, + 513, + 522, + 526, + 527 + ], + "479": [ + 69, + 77, + 85, + 93, + 101, + 109, + 117, + 125, + 132, + 140, + 148, + 156, + 164, + 172, + 180, + 188, + 196, + 204, + 212, + 220, + 228, + 236, + 244, + 251, + 259, + 267, + 275, + 283, + 291, + 299, + 307, + 315, + 323, + 331, + 339, + 347, + 355, + 362, + 370, + 378, + 386, + 394, + 402, + 410, + 418, + 426, + 434, + 442, + 450, + 458, + 466, + 473, + 481, + 489, + 497, + 505 + ], + "79556706": [ + 66, + 74, + 82, + 90, + 98, + 106, + 115, + 126, + 131, + 139, + 147, + 171, + 179, + 187, + 195, + 203, + 211, + 219, + 227, + 235, + 243, + 251, + 259, + 267, + 275, + 283, + 291, + 299, + 307, + 315, + 323, + 331, + 339, + 347, + 355, + 363, + 371, + 379, + 387, + 394, + 402, + 410, + 418, + 426, + 434, + 442, + 450, + 458, + 466, + 474, + 482, + 490, + 498, + 506, + 513, + 521, + 526, + 527 + ], + "79591671": [ + 68, + 76, + 83, + 91, + 99, + 107, + 115, + 123, + 131, + 139, + 146, + 154, + 162, + 170, + 178, + 194, + 202, + 209, + 217, + 225, + 233, + 241, + 249, + 257, + 264, + 272, + 280, + 288, + 296, + 304, + 319, + 327, + 335, + 351, + 359, + 367, + 374, + 382, + 390, + 398, + 406, + 414, + 422, + 430, + 437, + 445, + 453, + 461, + 469, + 477, + 484, + 492, + 500, + 508, + 516, + 524, + 527 + ], + "112202838": [ + 74, + 81, + 89, + 97, + 104, + 112, + 120, + 127, + 135, + 143, + 151, + 158, + 166, + 174, + 181, + 189, + 197, + 204, + 212, + 220, + 228, + 243, + 251, + 258, + 266, + 274, + 281, + 289, + 297, + 304, + 312, + 320, + 328, + 335, + 343, + 351, + 358, + 366, + 374, + 381, + 389, + 397, + 405, + 412, + 420, + 428, + 435, + 443, + 451, + 459, + 466, + 474, + 482, + 489, + 497, + 505, + 512, + 520, + 527 + ], + "79591593": [ + 73, + 81, + 89, + 96, + 103, + 111, + 119, + 126, + 134, + 141, + 149, + 157, + 164, + 172, + 180, + 188, + 195, + 203, + 211, + 218, + 226, + 234, + 242, + 257, + 265, + 273, + 280, + 288, + 296, + 319, + 327, + 335, + 343, + 351, + 359, + 367, + 383, + 391, + 399, + 407, + 415, + 422, + 430, + 438, + 446, + 454, + 462, + 470, + 478, + 486, + 494, + 502, + 518, + 526, + 527 + ], + "1175": [ + 67, + 76, + 84, + 91, + 107, + 115, + 123, + 131, + 139, + 147, + 155, + 163, + 171, + 179, + 187, + 195, + 203, + 211, + 219, + 227, + 235, + 243, + 251, + 259, + 267, + 275, + 283, + 291, + 299, + 307, + 315, + 323, + 331, + 339, + 347, + 355, + 363, + 372, + 380, + 388, + 396, + 404, + 412, + 420, + 428, + 436, + 444, + 452, + 460, + 468, + 476, + 484, + 492, + 500, + 508 + ], + "68161453": [ + 204, + 198, + 192, + 185, + 167, + 160, + 154, + 148, + 142, + 135, + 117, + 111, + 104, + 98, + 92, + 86, + 432, + 426, + 420, + 413, + 395, + 388, + 382, + 376, + 370, + 363, + 345, + 339, + 332, + 326, + 320, + 314 + ], + "68631984": [ + 255, + 247, + 240, + 232, + 223, + 215, + 207, + 199, + 191, + 183, + 175, + 167, + 159, + 150, + 142, + 134, + 126, + 117, + 109, + 101, + 201, + 209, + 216, + 224, + 451, + 443, + 435, + 427, + 419, + 411, + 403, + 395, + 387, + 378, + 370, + 362, + 354, + 345, + 337, + 329 + ] +} \ No newline at end of file diff --git a/atlas_densities/app/data/measurements/README.rst b/atlas_densities/app/data/measurements/README.rst index e2a5156..1dd24a8 100644 --- a/atlas_densities/app/data/measurements/README.rst +++ b/atlas_densities/app/data/measurements/README.rst @@ -2,22 +2,36 @@ Description =========== -This folder contains various experiment measurements related to cell densities. -Measurements can be average cell densities in number of cells per mm^3, but also +This folder contains various literature measurements related to cell densities. +Literature measurements can be average cell densities in number of cells per mm^3, but also number of cells or volumes. All measurement types are described in `atlas_densities.app.cell_densities`. Content ------- +The excel file `gaba_papers.xlsx` is a compilation of measurements from the scientific literature +collected in `Rodarie et al. (2021)`_. + +The two files `mmc1.xlsx` and `mmc3.xlsx` are unmodified copies of excel files from the +supplementary materials of `Kim et al. (2017)`_. + +The file `std_cells.json` contains the cell density standard deviation the BBP cell atlas for every +AIBS region. This file has been produced by Csabe Eroe (see `Eroe et al. (2018)`_), based on +multiple Nissl volumes, assuming that Nissl stains intensity depends linearly on volumetric cell +density across the mouse brain. This file is used to assign a "standard deviation" value to every +region where a pre-computed volumetric neuron density has been used to assign an average neuron +density. This holds for every region whose neurons are inhibitory only. (Those regions are listed in +`atlas_densities/app/data/measurements/homogenous_regions.csv`.) The file `measurements.csv` has been produced with the command line -`atlas-densities cell-densities compile-measurements` using `mmc3.xlsx`, `gaba_papers.xslx` and `non_density_measurement` as input. -The file `measurements.csv` encloses every measurement serving as input for the computation of cell densities -in the mouse brain. It is consumed by the command line -`atlas-densities cell-densities measurements-to-average-densities`. +`atlas-densities cell-densities compile-measurements` using `mmc3.xlsx`, `gaba_papers.xslx` and +`non_density_measurement` as input. +The file `measurements.csv` encloses every literature measurement serving as input for the +computation of inhibitory neuron type densities in the mouse brain. It is consumed by the command +line `atlas-densities cell-densities measurements-to-average-densities`. -Every new measurement should be added to `measurements.csv`. +Every new literature measurement should be added to `measurements.csv`. The file `homogenous_measurements.csv` is a manually maintained list of brain regions where neurons are all of the "same type", for instance all inhibitory or all excitatory. It is used in particular to set with 0.0 the average inhibitory neuron density of regions @@ -27,16 +41,10 @@ in regions where neurons are inhibitory only, using a precomputed volumetric neu The two previous files taken apart, the data files of this directory are static and should not be updated. -The excel file `gaba_papers.xlsx` is a compilation of measurements from the scientific literature -collected by Dimitri Rodarie (BBP). - -The two files `mmc1.xlsx` and `mmc3.xlsx` are unmodified copies of excel files from the supplementary materials -of `"Brain-wide Maps Reveal Stereotyped Cell-Type-Based Cortical Architecture and Subcortical Sexual Dimorphism"`_ by Kim et al., 2017. - -The file `non_density_measurement.csv` is a series of measurements which have been manually extracted -from `gaba_papers.xlsx` due to their peculiarities. - - -.. _`"Brain-wide Maps Reveal Stereotyped Cell-Type-Based Cortical Architecture and Subcortical Sexual Dimorphism"`: https://www.sciencedirect.com/science/article/pii/S0092867417310693?via%3Dihub +The file `non_density_measurement.csv` is a series of measurements which have been manually +extracted from `gaba_papers.xlsx` due to their peculiarities. +.. _`Kim et al. (2017)`: https://www.sciencedirect.com/science/article/pii/S0092867417310693?via%3Dihub +.. _Rodarie et al. (2021): https://www.biorxiv.org/content/10.1101/2021.11.20.469384v2 +.. _`Eroe et al. (2018)`: https://www.frontiersin.org/articles/10.3389/fninf.2018.00084/full diff --git a/atlas_densities/app/mtype_densities.py b/atlas_densities/app/mtype_densities.py index 3c5bc3e..28a2fad 100644 --- a/atlas_densities/app/mtype_densities.py +++ b/atlas_densities/app/mtype_densities.py @@ -58,7 +58,7 @@ MTYPES_PROFILES_REL_PATH = (DATA_PATH / "mtypes" / "density_profiles").relative_to(AD_PATH) MTYPES_PROBABILITY_MAP_REL_PATH = (DATA_PATH / "mtypes" / "probability_map").relative_to(AD_PATH) MTYPES_COMPOSITION_REL_PATH = (DATA_PATH / "mtypes" / "composition").relative_to(AD_PATH) -METADATA_PATH = DATA_PATH / "data" / "metadata" +METADATA_PATH = DATA_PATH / "metadata" METADATA_REL_PATH = METADATA_PATH.relative_to(AD_PATH) L = logging.getLogger(__name__) diff --git a/atlas_densities/combination/annotations_combinator.py b/atlas_densities/combination/annotations_combinator.py index 47821dd..5212b62 100644 --- a/atlas_densities/combination/annotations_combinator.py +++ b/atlas_densities/combination/annotations_combinator.py @@ -57,6 +57,30 @@ def is_ancestor_(id_1: int, id_2: int) -> bool: return np.vectorize(is_ancestor_, otypes=[bool])(annotation_1, annotation_2) +def combine_ccfv2_annotations( + brain_annotation_ccfv2: "voxcell.VoxelData", + fiber_annotation_ccfv2: "voxcell.VoxelData", +): + """Combine the ccfv2 annotation main file with its fibers + The ccfv2 fiber annotation file is required because fiber tracts + are missing from the ccfv2 brain annotation file, + These assumptions are based on the use case + annotation 2011 (Mouse CCF v2) + + Each annotation file has a resolution, either 10 um or 25 um. + The input files and the output file should all have the same resolution. + Args: + brain_annotation_ccfv2: reference annotation file. + fiber_annotation_ccfv2: fiber annotation. + + Returns: + VoxelData object holding the combined annotation 3D array. + """ + fiber_mask = fiber_annotation_ccfv2.raw > 0 + brain_annotation_ccfv2.raw[fiber_mask] = fiber_annotation_ccfv2.raw[fiber_mask] + return brain_annotation_ccfv2 + + def combine_annotations( region_map: voxcell.RegionMap, brain_annotation_ccfv2: voxcell.VoxelData, @@ -88,8 +112,9 @@ def combine_annotations( Returns: VoxelData object holding the combined annotation 3D array. """ - fiber_mask = fiber_annotation_ccfv2.raw > 0 - brain_annotation_ccfv2.raw[fiber_mask] = fiber_annotation_ccfv2.raw[fiber_mask] + brain_annotation_ccfv2 = combine_ccfv2_annotations( + brain_annotation_ccfv2, fiber_annotation_ccfv2 + ) brain_annotation_ccfv3_mask = brain_annotation_ccfv3.raw > 0 ccfv3_ids = np.unique(brain_annotation_ccfv3.raw[brain_annotation_ccfv3_mask]) missing_ids = np.isin(brain_annotation_ccfv2.raw, ccfv3_ids, invert=True) diff --git a/atlas_densities/densities/cell_density.py b/atlas_densities/densities/cell_density.py index aacdad4..bb60c48 100644 --- a/atlas_densities/densities/cell_density.py +++ b/atlas_densities/densities/cell_density.py @@ -8,7 +8,12 @@ from voxcell import RegionMap # type: ignore from atlas_densities.densities.cell_counts import cell_counts -from atlas_densities.densities.utils import compensate_cell_overlap, get_group_ids, get_region_masks +from atlas_densities.densities.utils import ( + compensate_cell_overlap, + get_group_ids, + get_region_masks, + normalize_intensity, +) def fix_purkinje_layer_intensity( @@ -95,7 +100,9 @@ def compute_cell_density( layer constraint of a constant number of cells per voxel. """ - nissl = compensate_cell_overlap(nissl, annotation, gaussian_filter_stdv=1.0, copy=False) + nissl = np.asarray(nissl, dtype=np.float64) + nissl = normalize_intensity(nissl, annotation, threshold_scale_factor=1.0, copy=False) + nissl = compensate_cell_overlap(nissl, annotation, gaussian_filter_stdv=-1.0, copy=False) group_ids = get_group_ids(region_map) region_masks = get_region_masks(group_ids, annotation) diff --git a/atlas_densities/densities/glia_densities.py b/atlas_densities/densities/glia_densities.py index 4e09b35..f8aeb4b 100644 --- a/atlas_densities/densities/glia_densities.py +++ b/atlas_densities/densities/glia_densities.py @@ -10,6 +10,7 @@ compensate_cell_overlap, constrain_cell_counts_per_voxel, get_group_ids, + normalize_intensity, ) L = logging.getLogger(__name__) @@ -119,7 +120,7 @@ def compute_glia_densities( # pylint: disable=too-many-arguments The overall glia density field is bounded by the `cell_density` field. It sums up to `glia_cell_count` after multiplication by the voxel volume. In addition, it respects some region-specific hints (fiber tracts and Purkinje layer). - For every glia cell type, the correspondging output density field is bounded by the overall + For every glia cell type, the corresponding output density field is bounded by the overall glia density field and sums up to its prescribed cell count when multiplied with by the voxel volume. @@ -127,11 +128,14 @@ def compute_glia_densities( # pylint: disable=too-many-arguments glia_densities = glia_intensities.copy() # The algorithm constraining cell counts per voxel requires double precision for glia_type in glia_densities: - glia_densities[glia_type] = np.asarray(glia_densities[glia_type], dtype=float) - cell_density = np.asarray(cell_density, dtype=float) + glia_densities[glia_type] = np.asarray(glia_densities[glia_type], dtype=np.float64) + cell_density = np.asarray(cell_density, dtype=np.float64) glia_densities["glia"] = compensate_cell_overlap( - glia_densities["glia"], annotation, gaussian_filter_stdv=1.0, copy=False + np.asarray(glia_densities["glia"], dtype=np.float64), + annotation, + gaussian_filter_stdv=-1.0, + copy=copy, ) L.info( "Computing overall glia density field with a target cell count of %d ...", @@ -145,11 +149,17 @@ def compute_glia_densities( # pylint: disable=too-many-arguments glia_densities["glia"], cell_density * voxel_volume, ) + placed_cells = np.zeros_like(glia_densities["glia"]) for glia_type in ["astrocyte", "oligodendrocyte"]: + glia_densities[glia_type] = normalize_intensity( + np.asarray(glia_densities[glia_type], dtype=np.float64), + annotation, + copy=copy, + ) glia_densities[glia_type] = compensate_cell_overlap( glia_densities[glia_type], annotation, - gaussian_filter_stdv=1.0, + gaussian_filter_stdv=2.0, copy=copy, ) cell_count = glia_cell_count * float(glia_proportions[glia_type]) @@ -161,16 +171,15 @@ def compute_glia_densities( # pylint: disable=too-many-arguments glia_densities[glia_type] = constrain_cell_counts_per_voxel( cell_count, glia_densities[glia_type], - glia_densities["glia"], + glia_densities["glia"] - placed_cells, copy=copy, ) - - # pylint: disable=fixme - # FIXME(Luc): The microglia density can be negative. - glia_densities["microglia"] = ( - glia_densities["glia"] - glia_densities["astrocyte"] - glia_densities["oligodendrocyte"] + placed_cells += glia_densities[glia_type] + L.info( + "Computing microglia density field with a target cell count of %d ...", + np.sum(glia_densities["glia"] - placed_cells), ) + glia_densities["microglia"] = glia_densities["glia"] - placed_cells for glia_type, cell_counts_per_voxel in glia_densities.items(): glia_densities[glia_type] = cell_counts_per_voxel / voxel_volume - return glia_densities diff --git a/atlas_densities/densities/inhibitory_neuron_density.py b/atlas_densities/densities/inhibitory_neuron_density.py index 787b222..41cd568 100644 --- a/atlas_densities/densities/inhibitory_neuron_density.py +++ b/atlas_densities/densities/inhibitory_neuron_density.py @@ -166,6 +166,33 @@ def compute_inhibitory_neuron_density( # pylint: disable=too-many-arguments list(group_ids["Cerebellar cortex"] & group_ids["Molecular layer"]), ), ) + inhibitory_neurons_mask = np.logical_or( + inhibitory_neurons_mask, + np.isin( + annotation, + region_map.find("Striatum", attr="name", with_descendants=True), + ), + ) + inhibitory_neurons_mask = np.logical_or( + inhibitory_neurons_mask, + np.isin( + annotation, + region_map.find( + "Reticular nucleus of the thalamus", attr="name", with_descendants=True + ), + ), + ) + # Cortical L1 regions + inhibitory_neurons_mask = np.logical_or( + inhibitory_neurons_mask, + np.isin( + annotation, + ( + region_map.find("Isocortex", attr="name", with_descendants=True) + & region_map.find("@.*1[ab]?$", attr="acronym", with_descendants=True) + ), + ), + ) assert isinstance(inhibitory_data["neuron_count"], int) diff --git a/atlas_densities/densities/utils.py b/atlas_densities/densities/utils.py index dfe1cec..97a30e3 100644 --- a/atlas_densities/densities/utils.py +++ b/atlas_densities/densities/utils.py @@ -75,7 +75,9 @@ def normalize_intensity( output_intensity = copy_array(marker_intensity, copy=copy) output_intensity -= outside_mean * threshold_scale_factor output_intensity[output_intensity < 0.0] = 0.0 - output_intensity /= np.max(output_intensity) + max_intensity = np.max(output_intensity) + if max_intensity > 0: + output_intensity /= max_intensity return output_intensity @@ -343,8 +345,10 @@ def get_group_ids(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict[s | region_map.find("Basic cell groups and regions", attr="name") | region_map.find("Cerebellum", attr="name") ) - molecular_layer_ids = region_map.find("@.*molecular layer", attr="name", with_descendants=True) cerebellar_cortex_ids = region_map.find("Cerebellar cortex", attr="name", with_descendants=True) + molecular_layer_ids = cerebellar_cortex_ids & region_map.find( + "@.*molecular layer", attr="name", with_descendants=True + ) rest_ids = region_map.find("root", attr="name", with_descendants=True) rest_ids -= cerebellum_group_ids | isocortex_group_ids diff --git a/tests/densities/test_glia_densities.py b/tests/densities/test_glia_densities.py index 66b7328..7b8a74e 100644 --- a/tests/densities/test_glia_densities.py +++ b/tests/densities/test_glia_densities.py @@ -75,6 +75,8 @@ def get_glia_input_data(glia_cell_count): } glia_densities["glia"] = np.random.random_sample(shape) + for glia_type, proportion in glia_proportions.items(): + glia_densities[glia_type][0][0][0] = 1e-5 # the outside voxels' intensity should be low return { "region_map": RegionMap.load_json(Path(TESTS_PATH, "1.json")), "annotation": np.arange(8000).reshape(shape), diff --git a/tests/densities/test_inhibitory_neuron_density.py b/tests/densities/test_inhibitory_neuron_density.py index 3a136a1..1be6e03 100644 --- a/tests/densities/test_inhibitory_neuron_density.py +++ b/tests/densities/test_inhibitory_neuron_density.py @@ -2,7 +2,7 @@ Unit tests for inhibitory cell density computation """ from pathlib import Path -from unittest.mock import patch +from unittest.mock import Mock, patch import numpy as np import numpy.testing as npt @@ -85,7 +85,10 @@ def test_compute_inhibitory_neuron_density(): "atlas_densities.densities.inhibitory_neuron_density.compensate_cell_overlap", side_effect=[data["gad1"], data["nrn1"]], ): - region_map = {} # Fake + # have any empty region map + region_map = Mock() + region_map.find = lambda name, attr, with_descendants: set() + inhibitory_neuron_density = tested.compute_inhibitory_neuron_density( region_map, annotation, diff --git a/tests/densities/test_utils.py b/tests/densities/test_utils.py index da3e9ce..13ff32b 100644 --- a/tests/densities/test_utils.py +++ b/tests/densities/test_utils.py @@ -108,7 +108,7 @@ def test_get_group_ids(): assert len(group_ids["Cerebellum group"] & group_ids["Isocortex group"]) == 0 assert len(group_ids["Cerebellum group"] & group_ids["Molecular layer"]) > 0 assert len(group_ids["Cerebellum group"] & group_ids["Purkinje layer"]) > 0 - assert len(group_ids["Isocortex group"] & group_ids["Molecular layer"]) > 0 + assert len(group_ids["Isocortex group"] & group_ids["Molecular layer"]) == 0 assert len(group_ids["Isocortex group"] & group_ids["Purkinje layer"]) == 0 assert len(group_ids["Isocortex group"] & group_ids["Rest"]) == 0 assert len(group_ids["Cerebellum group"] & group_ids["Rest"]) == 0