seung-lab/realneuralnetworks.jl

A unified framework for skeletonization, morphological analysis, and connectivity analysis.

3D neuron models extracted from EM image segmentation

Features

• sparse skeletonization. extract the neuron skeletons based on image segmentation (colored labels).
• morphological analysis with a lot of features.
• neural networks including synapses. Most current morphological analysis tools do not have synapses.
• neuron morphology similarity measurement using NBLAST algorithm.

Installation

This package is registered in Julia package management system. Simply following these steps to install it.

• type julia to enter command line interface.
• inside julia package REPL, hit ] to enter package mode, then type add RealNeuralNetworks

Usage

Analysis

All the analysis functions are demonstrated in Jupyter Notebooks.

skeletonization of proofread neurons

The skeletonization could be done using a script. The following command will show the parameters to use it. julia skeletonize.jl -h

Docker

build docker image

sudo docker build . -t realneuralnetworks

The secret files containing cloud authentication keys follows the cloudvolume secret format.

docker run -v /tmp:/tmp -v $HOME/.cloudvolume/secrets:/root/.cloudvolume/secrets --net=host realneuralnetworks julia skeletonize.jl -h

REPL in Julia

using RealNeuralNetworks.NodeNets
using RealNeuralNetworks.Neurons
using RealNeuralNetworks.SWCs

nodeNet = NodeNet(seg::Array{UInt32,3}; obj_id = convert(UInt32,77605))
neuron = Neuron( nodeNet )
swc = SWC(neuron)
SWCs.save(swc, tempname()*".swc")

Morphological Features

features represent a whole neuron

• arbor density
• total path length
• number of segment points
• Median segment length is the median dendritic segment length of all the segments starting and ending at irreducible nodes (in μm). Irreducible nodes are the points of the dendritic arbor corresponding to soma, branching points or terminal points.
• 3D sholl analysis.
• Hull area is the area of the tightest convex hull containing the z-projection of the dendritic arbor (in μm2).
• volume of the convex hull around all neurites
• Average angle is the mean of the positive angle between (parent node, node) and (node, child node) vectors, where node, parent node and child node are irreducible nodes (in radians).
• Average tortuosity is the average value of the ratio of the actual dendritic length to the Euclidean distance between irreducible nodes.
• Asymmetry is the distance of the soma node to the dendritic arbor (skeleton) centre of mass (in nm).
• Typical radius (λ) is the root-mean-square distance of dendritic arbor points to the centre of mass (in nm).
• fractal dimension.
• longest neurite length.
• Strahler number
• distribution of Euclidian distance of segmentes from soma (third principal component)
• distribution of Euclidian distance of segmentes from soma as a function of segment order (third principal component)
• number of segmentes per segment order (second principal component)
• distribution of morphological distance of segmentes from soma along the skeleton as a function of segment order (first principal component)

features represent segmentes in a neuron

• ratio of tail diameter to head. could be useful to identify spines.
• segment order
• segment length
• branching angle. computation using dot product
• tortuosity / curvature. caa4486b501f936743f781782e6561833da7e413
• distance to root path length
• segment asymmetry
• average radius. easy to compute with radius list.

Credit

The skeletonization was originally implemented in Matlab by Alexander Bae using TEASAR algorithm, which was translated to Julia by Nicholas Turner.

This package is developed at Princeton University and Flatiron Institute.

Reference

We have a publication for this repo:

@article{wuRealNeuralNetworksJlIntegrated2022,
	title = {{RealNeuralNetworks}.jl: {An} {Integrated} {Julia} {Package} for {Skeletonization}, {Morphological} {Analysis}, and {Synaptic} {Connectivity} {Analysis} of {Terabyte}-{Scale} {3D} {Neural} {Segmentations}},
	volume = {16},
	issn = {1662-5196},
	shorttitle = {{RealNeuralNetworks}.jl},
	url = {https://www.frontiersin.org/article/10.3389/fninf.2022.828169},
	urldate = {2022-03-02},
	journal = {Frontiers in Neuroinformatics},
	author = {Wu, Jingpeng and Turner, Nicholas and Bae, J. Alexander and Vishwanathan, Ashwin and Seung, H. Sebastian},
	year = {2022},
}

realneuralnetworks.jl's Issues

only compute NBLAST score from small neuron to large neuron

The NBLAST score seems not catching the similar cells if one neuron is much bigger than the other. take neuron 76631 and 82169 for example.

One simple solution is only computing the NBLAST score from small neuron to big neuron to reduce the number of mismatching.

Run femtocleaner

use BitArray rather than Array{Bool, 3} for binary image

BitArray used all the bits in a byte to save memory usage.

read swc/swc.bin file will ignore the node class

the node class, such as axon/dendrite, information will lost while loading.
The information was lost in transforming NodeNet to Neuron type.

the nodeClassList was not read out and used.
https://github.com/seung-lab/RealNeuralNetworks.jl/blob/master/src/Neurons.jl#L38

remove blob terminal branches

implemented, but not tested yet.
87f488a

the closest node could be wrong if the gap is large

it looks like a bug, the connected distance looks farther then the correct node.

considering the direction of branch could be helpful to enhance the accuracy.

skeletonization is broken

although the program runs, the output swc is not correct. It is not the same with previous version and there is an error while transforming to Neuron type:

connecting to a segment end...                                                                       
the closest segment is a terminal segment, merge the root segment of the other net                   
need to break the closest segment and rebuild connectivity matrix                                    
parentSegmentIdList = [6, 10]                                                                        
ERROR: LoadError: AssertionError: length(parentSegmentIdList) <= 1                                   
Stacktrace:                                                                                          
 [1] macro expansion at ./show.jl:556 [inlined]                                                      
 [2] merge(::Neuron, ::Neuron, ::Int64, ::Int64) at /usr/people/jingpeng/.julia/dev/RealNeuralNetwork
s/src/Neurons.jl:1121                                                                                
 [3] Neuron(::NodeNet) at /usr/people/jingpeng/.julia/dev/RealNeuralNetworks/src/Neurons.jl:75       
 [4] #trace#3(::String, ::UInt32, ::String, ::Tuple{Int64,Int64,Int64}, ::String, ::Function, ::Int64
) at /usr/people/jingpeng/.julia/dev/RealNeuralNetworks/scripts/skeletonize.jl:42                    
 [5] (::getfield(Main, Symbol("#kw##trace")))(::NamedTuple{(:swcDir, :mip, :voxelSize, :meshName, :se
gmentationLayer),Tuple{String,UInt32,Tuple{Int64,Int64,Int64},String,String}}, ::typeof(trace), ::Int
64) at ./none:0                                                                                      
 [6] macro expansion at ./logging.jl:322 [inlined]                                                   
 [7] main() at /usr/people/jingpeng/.julia/dev/RealNeuralNetworks/scripts/skeletonize.jl:103         
 [8] top-level scope at none:0                                                                       
 [9] include at ./boot.jl:317 [inlined]                                                              
 [10] include_relative(::Module, ::String) at ./loading.jl:1038                                      
 [11] include(::Module, ::String) at ./sysimg.jl:29                                                  
 [12] exec_options(::Base.JLOptions) at ./client.jl:239                                              
 [13] _start() at ./client.jl:432                                                                    
in expression starting at /usr/people/jingpeng/.julia/dev/RealNeuralNetworks/scripts/skeletonize.jl:1
37                                                                                                   

morphological features

features represent a whole neuron

• arbor density 5d39fb9
• total path length
• total frustum-based surface 4d9df08
• total frustum-based volume 4d9df08
• number of branch points c78104a
• Median branch length is the median dendritic segment length of all the segments starting and ending at irreducible nodes (in μm). Irreducible nodes are the points of the dendritic arbor corresponding to soma, branching points or terminal points.
• 3D sholl analysis. c493411
• Hull area is the area of the tightest convex hull containing the z-projection of the dendritic arbor (in μm2).
• volume of the convex hull around all neurites
• Average angle is the mean of the positive angle between (parent node, node) and (node, child node) vectors, where node, parent node and child node are irreducible nodes (in radians). 3f11d24
• Average tortuosity is the average value of the ratio of the actual dendritic length to the Euclidean distance between irreducible nodes. caa4486
• Asymmetry is the distance of the soma node to the dendritic arbor (skeleton) centre of mass (in nm). 6991059
• Typical radius (λ) is the root-mean-square distance of dendritic arbor points to the centre of mass (in nm). 87f488a
• fractal dimension. computed using cumulative-mass method. wiki. paper. paper2. 8923553 , 09f846b
• longest neurite length. 903fedf
• distribution of Euclidian distance of branches from soma (third principal component)
• distribution of Euclidian distance of branches from soma as a function of branch order (third principal component)
• number of branches per branch order (second principal component)
• distribution of morphological distance of branches from soma along the skeleton as a function of branch order (first principal component)

features represent branches in a neuron

• ratio of tail diameter to head. could be useful to identify spines. 8af50c0
• branch order
• branch length
• branching angle. 3f11d24 computation using dot product
• tortuosity / curvature. caa4486
• distance to root path length
• branch asymmetry
• average radius. easy to compute.

References:

1. Sümbül U, Song S, McCulloch K, Becker M, Lin B, Sanes JR, Masland RH, Seung HS. A genetic and computational approach to structurally classify neuronal types. Nature communications. 2014 Mar 24;5:3512. link
2. Schierwagen A, Villmann T, Alpár A, Gärtner U. Cluster analysis of cortical pyramidal neurons using som. InIAPR Workshop on Artificial Neural Networks in Pattern Recognition 2010 Apr 11 (pp. 120-130). Springer, Berlin, Heidelberg.
3. Cuntz H, Forstner F, Borst A, H\äusser M. The TREES Toolbox—Probing the Basis of Axonal and Dendritic Branching. Neuroinformatics. 2011;1–6.
4. Schierwagen A. Neuronal morphology: Shape characteristics and models. Neurophysiology. 2008;40(4):310–315.
5. Uylings HB., van Pelt J. Measures for quantifying dendritic arborizations. Network: Computation in Neural Systems. 2002;13(3):397–414. a good review paper
6. Wanner AA, Genoud C, Masudi T, Siksou L, Friedrich RW. Dense EM-based reconstruction of the interglomerular projectome in the zebrafish olfactory bulb. Nature neuroscience. 2016 Jun 1;19(6):816-25. also have clustering methods

use GeometryTypes to replace some type

such as Point3f0 as an alternative of node
and HyperRectangle as an alternative of bounding box.

this will make package more Julia and easier to interact with other Julia packages, such as GLVisualize

the total path length seems wrong.

in the valid cells of 11/30 consensus,
the total path length of all 1288 neurons is 534 micron
the total path length of 413 completed neurons is 207 micron

this result is suspicious and does not look right.

Ashwin had one cell that was ~600um path length. And the total path length of the 22 cells was ~1.2mm.

some branching point position is not accurate

when the neurite is not round (one axis is not as thick as other axies), and the neurite branching angle is small, the branching point is relatively "dragged" to root node.

example

multiple synapses can not attach to same node

this made a lot of synapses overwritten since one big post synapse can receive multiple boutton input. This is very common in the brain.

a merging subnet bug

only one cell (77652) was found with this problem.

distributed skeletonization

use SQS to distribute the skeleton tasks.

reconnecting could be wrong

neuron 76192 was wrong.

link

split out TEASAR part as an individual module

currently, it was merged to NodeNets, and the parameters were not well named following the Clean Code principles.

NBLAST metric with axon and dendrite compared separately

compare axon and dendrite separately with NBLAST

resampling of a branch

currently the nodes inside a branch is pretty dense, should resample to a fixed distance between neighboring nodes.

Reference:

trees toolbox

use Unit package rather than the raw values

https://github.com/ajkeller34/Unitful.jl

skeleton is not aligned with the mesh

since Nico fixed the meshing with an offset by one, now the mesh do not align with skeleton.

seung-lab/igneous#16

remove the branch completely inside the soma

there exist a redundent branch inside some

synaptic features

• distance between synapses. This could be computed from total path length and number of synapses.
• synapse distance to soma center.

use physical coordinate internally

currently, we are using voxel-based coordinate system internally, it might be better to use physical coordinates to compute metrics, such as length.

it is also much easier to get the radius correct.

use TimerOutputs to make the timer print prettier

https://github.com/KristofferC/TimerOutputs.jl

downsampling might not work correctly

I find that downsampling to 1 micron node distance makes some neuron length go to 0!
will need double check.

there exists some point with more than one parent

when I run skeletonization code, the problem was spotted by an assertion:

the closest segment is a terminal segment, merge the root segment of the other net                   
need to break the closest segment and rebuild connectivity matrix                                    
ERROR: LoadError: AssertionError: length(parentSegmentIdList) == 1                                   
Stacktrace:                                                                                          
 [1] get_parent_segment_id(::SparseArrays.SparseMatrixCSC{Bool,Int64}, ::Int64) at /usr/people/jingpe
ng/.julia/dev/RealNeuralNetworks/src/Neurons.jl:589                                                  
 [2] get_parent_segment_id at /usr/people/jingpeng/.julia/dev/RealNeuralNetworks/src/Neurons.jl:579 [
inlined]                                                                                             
 [3] Array{RealNeuralNetworks.SWCs.PointObjs.PointObj,1}(::Neuron) at /usr/people/jingpeng/.julia/dev
/RealNeuralNetworks/src/Neurons.jl:1432                                                              
 [4] #trace#3(::String, ::UInt32, ::String, ::Tuple{Int64,Int64,Int64}, ::String, ::Function, ::Int64
) at /usr/people/jingpeng/.julia/dev/RealNeuralNetworks/scripts/skeletonize.jl:45                    
 [5] (::getfield(Main, Symbol("#kw##trace")))(::NamedTuple{(:swcDir, :mip, :meshName, :voxelSize, :se
gmentationLayer),Tuple{String,UInt32,String,Tuple{Int64,Int64,Int64},String}}, ::typeof(trace), ::Int
64) at ./none:0                                                                                      
 [6] macro expansion at ./logging.jl:314 [inlined]                                                   
 [7] main() at /usr/people/jingpeng/.julia/dev/RealNeuralNetworks/scripts/skeletonize.jl:102         
 [8] top-level scope at none:0                                                                       
 [9] include at ./boot.jl:317 [inlined]                                                              
 [10] include_relative(::Module, ::String) at ./loading.jl:1038                                      
 [11] include(::Module, ::String) at ./sysimg.jl:29                                                  
 [12] exec_options(::Base.JLOptions) at ./client.jl:239                                              
 [13] _start() at ./client.jl:432                                                                    
in expression starting at /usr/people/jingpeng/.julia/dev/RealNeuralNetworks/scripts/skeletonize.jl:1
36                                                                                                   

such as neuron 79951

performance degradation of refactor-nodenet branch

the dbf weight function is about 1 min while it is about 0.2 sec in master branch!

index out of bound error

have already get swc file of more than 10 neurons, so this error happens occasionally.

skeletonization from global point cloud and dbf using TEASAR algorithm...
total number of points: 827520
ERROR: LoadError: ArgumentError: An index is out of bound.
 in sparsevec(::Array{Int64,1}, ::UnitRange{Int64}, ::UInt32, ::Function) at ./sparse/sparsevector.jl:117
 in create_node_lookup(::Array{UInt32,2}) at /home/nico/.julia/v0.5/TEASAR/src/NodeNets.jl:347
 in #NodeNet#3(::Array{Float64,1}, ::TEASAR.NodeNets.#alexs_penalty, ::Tuple{UInt32,UInt32,UInt32}, ::Type{T}, ::Array{UInt
32,2}) at /home/nico/.julia/v0.5/TEASAR/src/NodeNets.jl:75
 in (::Core.#kw#Type)(::Array{Any,1}, ::Type{TEASAR.NodeNets.NodeNet}, ::Array{UInt32,2}) at ./<missing>:0
 in macro expansion at ./util.jl:188 [inlined]
 in trace(::TEASAR.Manifests.Manifest, ::UInt32) at /home/nico/.julia/v0.5/TEASAR/src/Manifests.jl:75
 in #trace#1(::String, ::String, ::UInt32, ::Function, ::UInt32) at /home/nico/.julia/v0.5/TEASAR/scripts/skeletonize.jl:20
 in (::#kw##trace)(::Array{Any,1}, ::#trace, ::UInt32) at ./<missing>:0
 in macro expansion at /home/nico/.julia/v0.5/TEASAR/scripts/skeletonize.jl:90 [inlined]
 in macro expansion at /home/nico/.julia/v0.5/ProgressMeter/src/ProgressMeter.jl:478 [inlined]
 in main() at /home/nico/.julia/v0.5/TEASAR/scripts/skeletonize.jl:89
 in include_from_node1(::String) at ./loading.jl:488
 in process_options(::Base.JLOptions) at ./client.jl:265
 in _start() at ./client.jl:321
while loading /home/nico/.julia/v0.5/TEASAR/scripts/skeletonize.jl, in expression starting on line 99

dense skeletonization

currently, we are doing one neuron a time. the advantages are:

• less boundary effect
• avoid stitching across chunks, which could lead to errors

the disadvantages are:

• the IO is redundant since we only use a very sparse part of each chunk. the IO time is longer than computation.

getting the skeletons around the neuron, we might be able to use skeleton to do agglomeration. the AI need to answer a question, whether this piece should merge to this neuron or not.

this has a very low priority.

visualization in neuroglancer

Problem

currently, we can only save SWC file format and visualize it using some local tools. Since all other data were visualized in neuroglancer, it would be nice to be able to visualize in neuroglancer too.

Solution

chunk-wise storage

neuroglancer chunk manager has binary skeleton support, similar with mesh.

SWC

the Human Brain Project guys added swc support recently.

radius are not accurate for anisotropic images

although we have a voxel_size parameter now, but it only works for (1,1,1) voxel size now.

the mass center was estimated as center of all nodes

currently, there is no guarantee that the nodes are evenly distributed. After, evenly resample the nodes, the mass estimation should be more accurate.

this also will affect the estimation of fractal dimension.

#24

differentiate the score table of axon and dendrite

the axons seems follow the direction, but the dendrite seems not.
The simplest way is only using position for dendrite.
The more complex way is building two new tables for axon and dendrite separately.

update penalty function

Alex have a new penalty function, which can eliminate centerline jitter, and make the center line more smooth.

synapse and neuron connection representation

neurons connect with each other, so we should be able to represent them as connected.

missing branch

there exist some missing branches.

this was probably caused by removing points along the path.
https://github.com/seung-lab/RealNeuralNetworks.jl/blob/master/src/NodeNets.jl#L460

speedup fractal dimension computation

currently, it is pretty expensive to compute fractal dimension.

in the test neuron

get fractal dimension ...                                            
 28.322377 seconds (681.63 M allocations: 20.345 GiB, 9.42% gc time) 
fractalDimension = 1.3603973145963741                                

current implementation is based on nodes, we can change to branch based implementation. we can skip the branches that have starting point closer or farther than r +- path length .

connect the broken pieces

Problem

the segmentation is not one complete connected component, so the neuron is broken. The downsampling also enhanced the broken part.

the broken tree is not usable to post processing, such as feature extraction and clustering.

Solution

• simply find the nearest terminal point and connect them

eliminate the segments with single point

some segments only have one single point!

identify and remove spines

the spines should not be counted in morphological analysis. they should be able to be identified by some metrics.

• it is terminal branch
• inside the terminal branch, the terminating part have larger diameter and the starting part have smaller diameter. the ratio of terminal half radius and starting half radius might be a good metric to identify them.

an example

[Feature] visualize the synaptic distribution of a neuron.

like the one in Adrian's paper.

binary IO

currently, swc string IO is pretty slow (1-8 sec/neuron). It takes a long time to read out all the neurons, making debug slow.
adding the radius information to Neuroglancer binary format, should just work.

reset the root node to soma center

current root nodes were selected randomly, it would be more biologically plausible to make the root node in the soma center, which could be estimated by the largest diameter node.

load secret key from cloudvolume secret file

for task production code,
this is more secure.

upgrade to Julia 1.0

• use property to replace the get functions. This is really a nice feature of Julia 1.0!

transforming from NodeNet to Neuron increases the node number

there are 100 nodes in example.swc, but after transforming to Neuron type, the node number increased to 197.

swc = RealNeuralNetworks.SWCs.load("/usr/people/jingpeng/.julia/dev/RealNeuralNetworks/asset/example.swc")

nodeNet = RealNeuralNetworks.NodeNets.NodeNet(swc)

@show length(nodeNet.nodeList)
neuron = Neurons.Neuron(swc)
@show Neurons.get_num_nodes(neuron)

respect the position of neurons

use arbor density and check their overlap.
this reflect the synaptic connection.

compute DBF directly from segmentation

currently, we first transform pointcloud to binary image, and the DBF was computed from binary image. there is a few more memory allocation there.

should be able to compute DBF directly from the segmentation, could be a little bit more efficient.

the original implementation works fine now, this is very low priority.

binary IO

currently, swc files were stored as ascii text. the IO time is much longer than morphological computation time even though I am using async IO.

compress the ascii text using gzip and save binary files should speedup IO a lot.

out-of-core neuron reconstruction

Problem

a neuron can be pretty big and spans a large space, directly working one the big dense array will probably out of memory.

Solution

use the manifest to load chunks sparsely. record the offset of each chunk and form a list of them.

• load chunks sequentially (or asynchronously), do not load them all to work with large neurons
• binarize the segmentation chunk, and build point cloud and DBF (Distance from Boundary Field).
• concatenate the point clouds and DBFs for the whole neuron
• run remaining algorithms only on point clouds and DBFs. no dense array anymore.