Soma Detection Analysis of Whole-Brain Fluorescence Images

from brainlit.BrainLine.analyze_results import SomaDistribution
from brainlit.BrainLine.util import (
    json_to_points,
    download_subvolumes,
)
from brainlit.BrainLine.apply_ilastik import (
    ApplyIlastik,
    ApplyIlastik_LargeImage,
    plot_results,
    examine_threshold,
)
from brainlit.BrainLine.parse_ara import *
import xml.etree.ElementTree as ET
from brainlit.BrainLine.imports import *

%gui qt5

1. Before Using this notebook

1a. Install brainlit, and dependencies

1b. Write images to s3 using CloudReg

- e.g. python -m cloudreg.scripts.create_precomputed_volumes --s3_input_paths <path-to-stitched-images> --s3_output_paths  <s3-address-for-output>  --voxel_size <x-resolution> <y-resolution> <z-resolution> --num_procs <num-cpus> --resample_iso <boolean-to-resample-isotropically>

1c. Make point annotations in neuroglancer to identify subvolumes for validation (and possible training)

Instructions. JSON snippet:

,
{
"type":"pointAnnotation",
"name": "soma_val",
"points": []
},
{
"type":"pointAnnotation",
"name": "nonsoma_val",
"points":[]
}

1d. Update soma_data.json file

This file stores information on how to access neuroglancer data.

Data should be stored in the brain2paths dictionary, with entries like:

"<sample ID>": {
    "base_s3": "<Path to directory with layers with CloudVolume prependings (ending with forward slash)>",
    "base_local": "<Path to directory with layers with CloudVolume prependings (ending with forward slash)>",
    "val_info": {
        "url": "<neuroglancer URL>",
        "somas_layer": "<name of layer with coordinates on somas>",
        "nonsomas_layer": "<name of layer with coordinates on non-somas>",
    },
    "somas_atlas_url": "<neuroglancer URL entered after processing with a single annotation layer which contains points of soma detections>",
    "subtype": "<subtype>"
    #Optional:
    "train_info": {
        "url": "<neuroglancer URL>",
        "somas_layer": "<name of layer with coordinates on somas>",
        "nonsomas_layer": "<name of layer with coordinates on non-somas>",
    },
}

brainlit_path = Path(os.path.abspath(""))
brainlit_path = brainlit_path.parents[3]
print(f"Path to brainlit: {brainlit_path}")
data_file = (
    brainlit_path / "docs" / "notebooks" / "pipelines" / "BrainLine" / "soma_data.json"
)

brain = "test"  # brain ID
soma_data_dir = (
    str(
        brainlit_path
        / "docs"
        / "notebooks"
        / "pipelines"
        / "BrainLine"
        / "validation-data"
        / "soma"
    )
    + "/"
)  # path to directory where training/validation data should be stored

# Modify base path of test sample according to your system
with open(data_file, "r") as f:
    data = json.load(f)

base = (
    "precomputed://file://"
    + str(
        brainlit_path
        / "docs"
        / "notebooks"
        / "pipelines"
        / "BrainLine"
        / "example-data"
    )
    + "/"
)
data["brain2paths"]["test"]["base_s3"] = base
data["brain2paths"]["test"]["base_local"] = base

brain2paths = data["brain2paths"]

with open(data_file, "w") as f:
    json.dump(data, f, indent=4)

Validation: Steps 2, 4-5, 7 below can be done via script with: brainlit/BrainLine/scripts/soma_validate.py

2. Download benchmark data

*Inputs*

For the test sample to work, you should be serving it via neuroglancer’s cors_webserver e.g. with the command python cors_webserver.py -d "<path to brainlit>/brainlit/BrainLine/data/example" -p 9010

dataset_to_save = "val"  # train or val

antibody_layer = "antibody"
background_layer = "background"
endogenous_layer = "endogenous"

layer_names = [antibody_layer, background_layer, endogenous_layer]

Download data

download_subvolumes(
    soma_data_dir,
    brain_id=brain,
    data_file=data_file,
    layer_names=layer_names,
    dataset_to_save=dataset_to_save,
)

3. View downloaded data (optional)

*Inputs*

fname = Path(soma_data_dir) / f"brain{brain}" / "val" / "2005_1978_1841_pos.h5"
scale = [1.8, 1.8, 2]  # voxel size in microns

with h5py.File(fname, "r") as f:
    im = f.get("image_3channel")
    image_fg = im[0, :, :, :]
    image_bg = im[1, :, :, :]
    image_endo = im[2, :, :, :]

viewer = napari.Viewer(ndisplay=3)
viewer.add_image(image_fg, scale=scale)
viewer.add_image(image_bg, scale=scale)
viewer.add_image(image_endo, scale=scale)
viewer.scale_bar.visible = True
viewer.scale_bar.unit = "um"

4. Apply ilastik to validation data

Ilastik can be downloaded here.

You can do this programmatically (below), or you can use the ilastik GUI (which is sometimes faster).

* Inputs *

ilastik_project = str(
    brainlit_path
    / Path("experiments/BrainLine/data/models/soma/matt_soma_rabies_pix_3ch.ilp")
)  # path to ilastik model to be used
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)
brains = [brain]

applyilastik = ApplyIlastik(
    ilastik_path=ilastik_path,
    project_path=ilastik_project,
    brains_path=soma_data_dir,
    brains=brains,
)
applyilastik.process_subvols(ncpu=6)
# applyilastik.move_results()

5. Check Results

Image Subvolumes with Two Cell Bodies

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.

doubles = ["2005_1978_1841_pos"]

Validation

fig, _, _ = plot_results(
    data_dir=soma_data_dir,
    brain_ids=[brain],
    object_type="soma",
    positive_channel=0,
    doubles=doubles,
)
plt.show()

If results above are not adequate, improve model and try again

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.

Examine best threshold

examine_threshold(
    data_dir=soma_data_dir,
    brain_id=brain,
    threshold=0.76,
    object_type="soma",
    positive_channel=0,
    doubles=doubles,
)

6. Apply ilastik to whole image

This can be done via a script with: brainlit/BrainLine/soma_detect_image

* Inputs *

threshold = 0.76  # threshold to use for ilastik
data_dir = (
    soma_data_dir + "brainr_temp/"
)  # "/data/tathey1/matt_wright/brainr_temp/"  # directory to store temporary subvolumes for segmentation
results_dir = (
    soma_data_dir + "brainr_results/"
)  # directory to store coordinates of soma detections

min_coords = [2000, 1900, 1800]
max_coords = [2100, 2000, 1900]  # -1 if you want to process the whole dimension

ncpu = 1  # 16  # number of cores to use for detection
chunk_size = [100, 100, 100]  # [256, 256, 300]

layer_names = [antibody_layer, background_layer, endogenous_layer]

ilastik_largeimage = ApplyIlastik_LargeImage(
    ilastik_path=ilastik_path,
    ilastik_project=ilastik_project,
    data_file=data_file,
    results_dir=results_dir,
    ncpu=1,
)

ilastik_largeimage.apply_ilastik_parallel(
    brain_id=brain,
    layer_names=layer_names,
    threshold=threshold,
    data_dir=data_dir,
    chunk_size=chunk_size,
    min_coords=min_coords,
    max_coords=max_coords,
)

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:

python cors_webserver.py -d "<path-to-brainlit>/brainlit/brainlit/BrainLine/data/example" -p 9010

which needs to be run in the neuroglancer folder (git clone from here: https://github.com/google/neuroglancer)

ilastik_largeimage.collect_soma_results(
    brain_id="test"
)  # will not work if you are not serving the data (see neuroglancer cors_webserver info above)

7. Prepare Transformed Layers

atlas_vol = CloudVolume(
    "precomputed://https://open-neurodata.s3.amazonaws.com/ara_2016/sagittal_10um/annotation_10um_2017"
)
for layer in [
    antibody_layer,
    background_layer,
]:  # axon_mask is transformed into an image because nearest interpolation doesnt work well after downsampling
    layer_path = brain2paths[brain]["base_s3"] + layer + "_transformed"
    info = CloudVolume.create_new_info(
        num_channels=1,
        layer_type="image",
        data_type="uint16",  # Channel images might be 'uint8'
        encoding="raw",  # raw, jpeg, compressed_segmentation, fpzip, kempressed
        resolution=atlas_vol.resolution,  # Voxel scaling, units are in nanometers
        voxel_offset=atlas_vol.voxel_offset,
        chunk_size=[32, 32, 32],  # units are voxels
        volume_size=atlas_vol.volume_size,  # e.g. a cubic millimeter dataset
    )
    vol_mask = CloudVolume(layer_path, info=info)
    vol_mask.commit_info()

8. Register volume and transform data to atlas space using CloudReg

8a. You need to find an initial affine alignment using cloudreg.scripts.registration.get_affine_matrix. For example:

A link to the ARA parcellation is:

precomputed://https://open-neurodata.s3.amazonaws.com/ara_2016/sagittal_10um/annotation_10um_2017

And some python commands to help with affine alignment is:

from cloudreg.scripts.registration import get_affine_matrix

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")

8b. Run registration using cloudreg.scripts.registration. For example:

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

8c. Transform data to atlas space using CloudReg

example commands below

Soma coordinates

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

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.

Image

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

Analysis: Steps 9-10 below can be done via a script with: brainlit/BrainLine/scripts/soma_analyze.py

9. View somas in brain space

*Inputs*

colors = {
    "test_type": "red",
}  # colors for different sample types
symbols = ["o", "+", "^", "vbar"]
brain_ids = ["test"]
fold_on = False

ontology_file = (
    brainlit_path / "brainlit" / "Brainline" / "data" / "ara_structure_ontology.json"
)

sd = SomaDistribution(
    brain_ids=brain_ids, data_file=data_file, ontology_file=ontology_file
)

f = sd.napari_coronal_section(
    z=1073, subtype_colors=colors, symbols=symbols, fold_on=fold_on, plot_type="plt"
)
plt.show()

sd.napari_coronal_section(
    z=1073, subtype_colors=colors, symbols=symbols, fold_on=fold_on
)

from vedo import embedWindow

embedWindow(None)  # needed for brainrender in jupyter notebook

sd.brainrender_somas(
    subtype_colors=colors, brain_region="VERM"
)  # might take a while to download the atlas for your first run

10. Display bar charts

regions = [
    512,  # cerebellum
    688,  # cerebral cortex
    698,  # olfactory areas
    1089,  # hippocampal formation
    # 583, # claustrum
    477,  # striatum
    # 803, # pallidum
    351,  # bed nuclei of stria terminalis
    # 703, #cortical subplate
    1097,  # hypothalamus
    549,  # thalamus
    186,  # lateral habenula
    519,  # cerebellar nuclei
    313,  # midbrain
    1065,  # hindbrain
]  # allen atlas region IDs to be shown
# see: https://connectivity.brain-map.org/projection/experiment/480074702?imageId=480075280&initImage=TWO_PHOTON&x=17028&y=11704&z=3

composite_regions = {
    "Amygdalar Nuclei": [131, 295, 319, 780]
}  # Custom composite allen regions where key is region name and value is list of allen regions

sd.region_barchart(regions, composite_regions=composite_regions, normalize_region=512)