{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Soma Detection Analysis of Whole-Brain Fluorescence Images" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from brainlit.BrainLine.analyze_results import SomaDistribution\n", "from brainlit.BrainLine.util import (\n", " json_to_points,\n", " download_subvolumes,\n", ")\n", "from brainlit.BrainLine.apply_ilastik import (\n", " ApplyIlastik,\n", " ApplyIlastik_LargeImage,\n", " plot_results,\n", " examine_threshold,\n", ")\n", "from brainlit.BrainLine.parse_ara import *\n", "import xml.etree.ElementTree as ET\n", "from brainlit.BrainLine.imports import *\n", "\n", "%gui qt5" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Before Using this notebook" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 1a. Install brainlit, and dependencies\n", "### 1b. Write images to s3 using CloudReg\n", " - e.g. python -m cloudreg.scripts.create_precomputed_volumes --s3_input_paths --s3_output_paths --voxel_size --num_procs --resample_iso \n", "### 1c. Make point annotations in neuroglancer to identify subvolumes for validation (and possible training)\n", "[Instructions](https://neurodata.io/help/neuroglancer-pt-annotations/). JSON snippet:\n", "\n", " ,\n", " {\n", " \"type\":\"pointAnnotation\",\n", " \"name\": \"soma_val\",\n", " \"points\": []\n", " },\n", " {\n", " \"type\":\"pointAnnotation\",\n", " \"name\": \"nonsoma_val\",\n", " \"points\":[]\n", " }\n", "### 1d. Update soma_data.json file\n", "\n", "This file stores information on how to access neuroglancer data.\n", "\n", "Data should be stored in the `brain2paths` dictionary, with entries like:\n", "\n", " \"\": {\n", " \"base_s3\": \"\",\n", " \"base_local\": \"\",\n", " \"val_info\": {\n", " \"url\": \"\",\n", " \"somas_layer\": \"\",\n", " \"nonsomas_layer\": \"\",\n", " },\n", " \"somas_atlas_url\": \"\",\n", " \"subtype\": \"\"\n", " #Optional:\n", " \"train_info\": {\n", " \"url\": \"\",\n", " \"somas_layer\": \"\",\n", " \"nonsomas_layer\": \"\",\n", " },\n", " }" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "brainlit_path = Path(os.path.abspath(\"\"))\n", "brainlit_path = brainlit_path.parents[3]\n", "print(f\"Path to brainlit: {brainlit_path}\")\n", "data_file = (\n", " brainlit_path / \"docs\" / \"notebooks\" / \"pipelines\" / \"BrainLine\" / \"soma_data.json\"\n", ")\n", "\n", "brain = \"test\" # brain ID\n", "soma_data_dir = (\n", " str(\n", " brainlit_path\n", " / \"docs\"\n", " / \"notebooks\"\n", " / \"pipelines\"\n", " / \"BrainLine\"\n", " / \"validation-data\"\n", " / \"soma\"\n", " )\n", " + \"/\"\n", ") # path to directory where training/validation data should be stored\n", "\n", "# Modify base path of test sample according to your system\n", "with open(data_file, \"r\") as f:\n", " data = json.load(f)\n", "\n", "base = (\n", " \"precomputed://file://\"\n", " + str(\n", " brainlit_path\n", " / \"docs\"\n", " / \"notebooks\"\n", " / \"pipelines\"\n", " / \"BrainLine\"\n", " / \"example-data\"\n", " )\n", " + \"/\"\n", ")\n", "data[\"brain2paths\"][\"test\"][\"base_s3\"] = base\n", "data[\"brain2paths\"][\"test\"][\"base_local\"] = base\n", "\n", "brain2paths = data[\"brain2paths\"]\n", "\n", "with open(data_file, \"w\") as f:\n", " json.dump(data, f, indent=4)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Validation: Steps 2, 4-5, 7 below can be done via script with: `brainlit/BrainLine/scripts/soma_validate.py`" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Download benchmark data" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### \\*Inputs\\*\n", "\n", "For the test sample to work, you should be serving it via neuroglancer's [cors_webserver](https://github.com/google/neuroglancer) e.g. with the command ``python cors_webserver.py -d \"/brainlit/BrainLine/data/example\" -p 9010``" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dataset_to_save = \"val\" # train or val\n", "\n", "antibody_layer = \"antibody\"\n", "background_layer = \"background\"\n", "endogenous_layer = \"endogenous\"\n", "\n", "layer_names = [antibody_layer, background_layer, endogenous_layer]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Download data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "download_subvolumes(\n", " soma_data_dir,\n", " brain_id=brain,\n", " data_file=data_file,\n", " layer_names=layer_names,\n", " dataset_to_save=dataset_to_save,\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 3. View downloaded data (optional)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### \\*Inputs\\*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fname = Path(soma_data_dir) / f\"brain{brain}\" / \"val\" / \"2005_1978_1841_pos.h5\"\n", "scale = [1.8, 1.8, 2] # voxel size in microns" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "with h5py.File(fname, \"r\") as f:\n", " im = f.get(\"image_3channel\")\n", " image_fg = im[0, :, :, :]\n", " image_bg = im[1, :, :, :]\n", " image_endo = im[2, :, :, :]\n", "\n", "viewer = napari.Viewer(ndisplay=3)\n", "viewer.add_image(image_fg, scale=scale)\n", "viewer.add_image(image_bg, scale=scale)\n", "viewer.add_image(image_endo, scale=scale)\n", "viewer.scale_bar.visible = True\n", "viewer.scale_bar.unit = \"um\"" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Apply ilastik to validation data" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Ilastik can be downloaded [here](https://www.ilastik.org/index.html).\n", "\n", "You can do this programmatically (below), or you can use the ilastik GUI (which is sometimes faster)." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### \\* Inputs \\*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ilastik_project = str(\n", " brainlit_path\n", " / Path(\"experiments/BrainLine/data/models/soma/matt_soma_rabies_pix_3ch.ilp\")\n", ") # path to ilastik model to be used\n", "ilastik_path = \"/Applications/ilastik-1.4.0b21-OSX.app/Contents/ilastik-release/run_ilastik.sh\" # path to ilastik executable (see here: https://www.ilastik.org/documentation/basics/headless.html)\n", "brains = [brain]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "applyilastik = ApplyIlastik(\n", " ilastik_path=ilastik_path,\n", " project_path=ilastik_project,\n", " brains_path=soma_data_dir,\n", " brains=brains,\n", ")\n", "applyilastik.process_subvols(ncpu=6)\n", "# applyilastik.move_results()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Check Results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Image Subvolumes with Two Cell Bodies\n", "\n", "If image subvolumes contain two cell bodies, the performance statistics can be adjusted using the `doubles` variable, which is a list of file name strings. Enter the file name (without extension) of the subvolumes with two cell bodies." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "doubles = [\"2005_1978_1841_pos\"]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Validation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, _, _ = plot_results(\n", " data_dir=soma_data_dir,\n", " brain_ids=[brain],\n", " object_type=\"soma\",\n", " positive_channel=0,\n", " doubles=doubles,\n", ")\n", "plt.show()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### If results above are not adequate, improve model and try again\n", "\n", "In my case, I identify more subvolumes from the sample at hand using the same process as for validation data, and add it as training data to the model and retrain." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Examine best threshold" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "examine_threshold(\n", " data_dir=soma_data_dir,\n", " brain_id=brain,\n", " threshold=0.76,\n", " object_type=\"soma\",\n", " positive_channel=0,\n", " doubles=doubles,\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Apply ilastik to whole image\n", "## This can be done via a script with: `brainlit/BrainLine/soma_detect_image`" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### \\* Inputs \\*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "threshold = 0.76 # threshold to use for ilastik\n", "data_dir = (\n", " soma_data_dir + \"brainr_temp/\"\n", ") # \"/data/tathey1/matt_wright/brainr_temp/\" # directory to store temporary subvolumes for segmentation\n", "results_dir = (\n", " soma_data_dir + \"brainr_results/\"\n", ") # directory to store coordinates of soma detections\n", "\n", "min_coords = [2000, 1900, 1800]\n", "max_coords = [2100, 2000, 1900] # -1 if you want to process the whole dimension\n", "\n", "ncpu = 1 # 16 # number of cores to use for detection\n", "chunk_size = [100, 100, 100] # [256, 256, 300]\n", "\n", "layer_names = [antibody_layer, background_layer, endogenous_layer]\n", "\n", "ilastik_largeimage = ApplyIlastik_LargeImage(\n", " ilastik_path=ilastik_path,\n", " ilastik_project=ilastik_project,\n", " data_file=data_file,\n", " results_dir=results_dir,\n", " ncpu=1,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ilastik_largeimage.apply_ilastik_parallel(\n", " brain_id=brain,\n", " layer_names=layer_names,\n", " threshold=threshold,\n", " data_dir=data_dir,\n", " chunk_size=chunk_size,\n", " min_coords=min_coords,\n", " max_coords=max_coords,\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Before this step you will need to make sure that the original data is being served to neuroglancer. For example, in this case, our data is local so we can serve it with the command:\n", "\n", "`python cors_webserver.py -d \"/brainlit/brainlit/BrainLine/data/example\" -p 9010`\n", "\n", "which needs to be run in the neuroglancer folder (git clone from here: https://github.com/google/neuroglancer)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ilastik_largeimage.collect_soma_results(\n", " brain_id=\"test\"\n", ") # will not work if you are not serving the data (see neuroglancer cors_webserver info above)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 7. Prepare Transformed Layers" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "atlas_vol = CloudVolume(\n", " \"precomputed://https://open-neurodata.s3.amazonaws.com/ara_2016/sagittal_10um/annotation_10um_2017\"\n", ")\n", "for layer in [\n", " antibody_layer,\n", " background_layer,\n", "]: # axon_mask is transformed into an image because nearest interpolation doesnt work well after downsampling\n", " layer_path = brain2paths[brain][\"base_s3\"] + layer + \"_transformed\"\n", " info = CloudVolume.create_new_info(\n", " num_channels=1,\n", " layer_type=\"image\",\n", " data_type=\"uint16\", # Channel images might be 'uint8'\n", " encoding=\"raw\", # raw, jpeg, compressed_segmentation, fpzip, kempressed\n", " resolution=atlas_vol.resolution, # Voxel scaling, units are in nanometers\n", " voxel_offset=atlas_vol.voxel_offset,\n", " chunk_size=[32, 32, 32], # units are voxels\n", " volume_size=atlas_vol.volume_size, # e.g. a cubic millimeter dataset\n", " )\n", " vol_mask = CloudVolume(layer_path, info=info)\n", " vol_mask.commit_info()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 8. Register volume and transform data to atlas space using CloudReg" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 8a. You need to find an initial affine alignment using cloudreg.scripts.registration.get_affine_matrix. For example: " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "A link to the ARA parcellation is:\n", "\n", "`precomputed://https://open-neurodata.s3.amazonaws.com/ara_2016/sagittal_10um/annotation_10um_2017`\n", "\n", "And some python commands to help with affine alignment is:\n", "\n", "```\n", "from cloudreg.scripts.registration import get_affine_matrix\n", "\n", "get_affine_matrix([1,1,1], [15,0,0], \"PIR\", \"RAI\", 1.15, \"precomputed://https://open-neurodata.s3.amazonaws.com/ara_2016/sagittal_10um/annotation_10um_2017\")\n", "```" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 8b. Run registration using cloudreg.scripts.registration. For example:" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "python -m cloudreg.scripts.registration -input_s3_path precomputed://s3://smartspim-precomputed-volumes/2022_02_02/8604/Ch_561 --output_s3_path precomputed://s3://smartspim-precomputed-volumes/2022_02_02/8604/atlas_to_target --atlas_s3_path https://open-neurodata.s3.amazonaws.com/ara_2016/sagittal_50um/average_50um --parcellation_s3_path https://open-neurodata.s3.amazonaws.com/ara_2016/sagittal_10um/annotation_10um_2017 --atlas_orientation PIR -orientation ARS --rotation 0 0 0 --translation 0 0 0 --fixed_scale .7 -log_s3_path precomputed://s3://smartspim-precomputed-volumes/2022_02_02/8604/atlas_to_target --missing_data_correction True --grid_correction False --bias_correction True --regularization 5000.0 --iterations 3000 --registration_resolution 100\n", "```" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 8c. Transform data to atlas space using CloudReg\n", "example commands below" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Soma coordinates\n", "\n", "```\n", "python -m cloudreg.scripts.transform_points --target_viz_link \"https://viz.neurodata.io/?json_url=https://json.neurodata.io/v1?NGStateID=pLuy0SXtVkKT0A\" --atlas_viz_link \"https://ara.viz.neurodata.io/?json_url=https://json.neurodata.io/v1?NGStateID=HvyNDGaPsd1wyg\" --affine_path \"/Users/thomasathey/Documents/mimlab/mouselight/brainlit_parent/downloop_1_A.mat\" --velocity_path \"/Users/thomasathey/Documents/mimlab/mouselight/brainlit_parent/downloop_1_v.mat\" --transformation_direction atlas\n", "```\n", "\n", "This will produce a neuroglancer link with the transformed soma coordinates, which should be added to `soma_data.py` under the `somas_atlas_url` key. Then the code below, or `soma_brainrender.py`, can be used to visualize the data." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Image\n", "\n", "```\n", "python -m cloudreg.scripts.transform_data --target_layer_source precomputed://s3://smartspim-precomputed-volumes/2022_09_20/887/Ch_647 --transformed_layer_source precomputed://s3://smartspim-precomputed-volumes/2022_09_20/887/Ch_647_transformed --affine_path /cis/home/tathey/887_Ch_561_registration/downloop_1_A.mat --velocity_path /cis/home/tathey/887_Ch_561_registration/downloop_1_v.mat\n", "```" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis: Steps 9-10 below can be done via a script with: `brainlit/BrainLine/scripts/soma_analyze.py`" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 9. View somas in brain space" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### \\*Inputs\\*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "colors = {\n", " \"test_type\": \"red\",\n", "} # colors for different sample types\n", "symbols = [\"o\", \"+\", \"^\", \"vbar\"]\n", "brain_ids = [\"test\"]\n", "fold_on = False\n", "\n", "ontology_file = (\n", " brainlit_path / \"brainlit\" / \"Brainline\" / \"data\" / \"ara_structure_ontology.json\"\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sd = SomaDistribution(\n", " brain_ids=brain_ids, data_file=data_file, ontology_file=ontology_file\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = sd.napari_coronal_section(\n", " z=1073, subtype_colors=colors, symbols=symbols, fold_on=fold_on, plot_type=\"plt\"\n", ")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sd.napari_coronal_section(\n", " z=1073, subtype_colors=colors, symbols=symbols, fold_on=fold_on\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from vedo import embedWindow\n", "\n", "embedWindow(None) # needed for brainrender in jupyter notebook\n", "\n", "sd.brainrender_somas(\n", " subtype_colors=colors, brain_region=\"VERM\"\n", ") # might take a while to download the atlas for your first run" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 10. Display bar charts" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regions = [\n", " 512, # cerebellum\n", " 688, # cerebral cortex\n", " 698, # olfactory areas\n", " 1089, # hippocampal formation\n", " # 583, # claustrum\n", " 477, # striatum\n", " # 803, # pallidum\n", " 351, # bed nuclei of stria terminalis\n", " # 703, #cortical subplate\n", " 1097, # hypothalamus\n", " 549, # thalamus\n", " 186, # lateral habenula\n", " 519, # cerebellar nuclei\n", " 313, # midbrain\n", " 1065, # hindbrain\n", "] # allen atlas region IDs to be shown\n", "# see: https://connectivity.brain-map.org/projection/experiment/480074702?imageId=480075280&initImage=TWO_PHOTON&x=17028&y=11704&z=3\n", "\n", "composite_regions = {\n", " \"Amygdalar Nuclei\": [131, 295, 319, 780]\n", "} # Custom composite allen regions where key is region name and value is list of allen regions\n", "\n", "sd.region_barchart(regions, composite_regions=composite_regions, normalize_region=512)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 2 }