Note
Click here to download the full example code
Analyzing Neuron Morphology#
This tutorial will give you an overview of how to analyze neuron morphology.
Disclaimer: As you might imagine some properties can be gathered for all/most neuron types while others will only work on specific types. For example, topological properties such as cable length, branch points, etc. are easy to get for a skeleton (i.e. a TreeNeuron) but not for an image (i.e. a VoxelNeuron). Make sure to check the respective function's docstring for details!
Basic Properties#
NAVis provides some basic morphometrics as neuron properties. Others need to be computed using a specific function (e.g. navis.tortuosity) - mostly because they require/allow some parameters to be set.
With that in mind, let's dive right in:
import navis
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_color_codes("muted")
Accessing attributes for a single neuron:
# Load a single skeleton (i.e. a TreeNeuron)
n = navis.example_neurons(n=1, kind="skeleton")
print(f"This neuron has {n.cable_length} of cable")
Out:
This neuron has 266476.875 of cable
If your neuron has its units properties set, you can also combine those:
print(f"This neuron has {(n.cable_length * n.units).to('microns')} of cable")
Out:
This neuron has 2131.815 micron of cable
Accessing the same properties via a navis.NeuronList will return an array:
nl = navis.example_neurons(kind="skeleton")
print(f"Cable lengths for neurons in this list:\n{nl.cable_length}")
Out:
Cable lengths for neurons in this list:
[266476.88 304332.66 274703.38 286522.47 291265.3 ]
navis.NeuronList.summary is a useful way to collect some of the basic parameters:
df = nl.summary()
df.head()
For MeshNeurons the available properties look different. For example, you can get its volume:
# Load a single MeshNeuron
m = navis.example_neurons(n=1, kind="mesh")
print(f"Neuron volume: {m.volume}")
# Again, we can use the `units` property to convert:
print(f"Neuron volume: {(m.volume * n.units **3).to('microns ** 3')}")
Out:
Neuron volume: 1291610825.1668377
Neuron volume: 661.3047424854211 micron ** 3
For topological properties, we need to convert to skeleton. The fastest way is to simply access the MeshNeuron's .skeleton property:
print(f"Mesh cable length: {m.skeleton.cable_length * m.units}")
Out:
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/skeletor/skeletonize/wave.py:198: DeprecationWarning:
Graph.clusters() is deprecated; use Graph.connected_components() instead
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/skeletor/skeletonize/wave.py:228: DeprecationWarning:
Graph.shortest_paths() is deprecated; use Graph.distances() instead
Mesh cable length: 1351124.75 nanometer
Note
The above .skeleton property is simply an automatically generated TreeNeuron representation of the mesh. It uses sensible defaults but as said initially: it's good practice to create and check the skeleton yourself via navis.skeletonize.
Importantly, some NAVis functions (e.g. navis.segment_analysis, see below) that accept MeshNeurons as input, really use this skeleton representation under-the-hood.
The skeleton representation of the mesh lets us access many toplogical properties:
m.skeleton.n_leafs
Out:
155
m.skeleton.n_branches
Out:
122
You may have already noticed here and in other examples the use of n_{property} (e.g. n.n_leafs). These are in fact generic: you can use any n_... and - assuming that property exists - NAVis will return a count:
m.n_vertices
Out:
6309
m.skeleton.n_nodes
Out:
1058
# Illustrate with a random property
m.my_counts = [1, 2, 3, 4, 5]
m.n_my_counts
Out:
5
Segment Analysis#
navis.segment_analysis is a great entry point for collecting a bunch of morphometrics for your neuron(s) of interest. It returns Strahler index, cable length, distance to root, radius and tortuosity for each linear segment:
sa = navis.segment_analysis(m)
sa.head()
# See if segment length correlates with radius
ax = sns.scatterplot(
data=sa, x="length", y="radius_mean", size="strahler_index", alpha=0.7
)
ax.set_xscale("log")
ax.set_yscale("log")
Sholl Analysis#
For an example of a Sholl analyses, check out the MICrONS tutorial.
Geodesic Distances#
Working with Euclidean distances is straight forward and we won't cover this extensively but here is an example where we are measuring the average distances between a node and its parent (= the sampling rate):
import numpy as np
# Get nodes but remove the root (has no parent)
nodes = nl[0].nodes[nl[0].nodes.parent_id > 0]
# Get the x/y/z coordinates of all nodes (except root)
node_locs = nodes[["x", "y", "z"]].values
# For each node, get its parent's location
parent_locs = (
nl[0].nodes.set_index("node_id").loc[nodes.parent_id.values, ["x", "y", "z"]].values
)
# Calculate Euclidean distances
distances = np.sqrt(np.sum((node_locs - parent_locs) ** 2, axis=1))
# Use the neuron's units to convert into nm
distances = distances * nl[0].units
print(
f"Mean distance between nodes: {np.mean(distances):.2f} (+/- {np.std(distances):.2f})"
)
Out:
Mean distance between nodes: 477.56 nanometer (+/- 361.10 nanometer)
What if you wanted to know the distance between the soma and all terminal nodes? In that case Euclidean distance would be insufficient as the neuron is not a straight line. Instead, you need the geodesic, the "along-the-arbor" distance.
NAVis comes with a couple functions that help you get geodesic distances. For single node-to-node queries, navis.dist_between should be sufficient:
n = nl[0]
end = n.nodes[n.nodes.type == "end"].node_id.values[0]
d_geo = navis.dist_between(n, n.soma, end) * n.units
print(f"Euclidean distance between soma and terminal node {end}: {d_geo:.2f}")
Out:
Euclidean distance between soma and terminal node 465: 444096.17 nanometer
Let's visualize this:
import networkx as nx
# First we need to find the path between the soma and the terminal node
path = nx.shortest_path(n.graph.to_undirected(), n.soma, end)
# Get coordinates for the path
path_co = n.nodes.set_index("node_id").loc[path, ["x", "y", "z"]].copy()
# Add a small offset
path_co.x += 500
path_co.y -= 500
# Plot neuron
fig, ax = navis.plot2d(n, c="blue", method="2d", view=("x", "-z"))
# Add geodesic path
ax.plot(path_co.x, path_co.z, c="r", ls="--")
# Add Euclidean path
end_loc = n.nodes.set_index("node_id").loc[end, ["x", "y", "z"]]
soma_loc = n.nodes.set_index("node_id").loc[n.soma, ["x", "y", "z"]]
ax.plot([soma_loc.x, end_loc.x], [soma_loc.z, end_loc.z], c="g", ls="--")
d_eucl = np.sqrt(np.sum((end_loc - soma_loc) ** 2)) * n.units
# Annotate distances
_ = ax.text(
x=0.1,
y=0.3,
s=f"Euclidean distance:\n{d_eucl.to_compact():.0f}",
transform=ax.transAxes,
c="g",
)
_ = ax.text(
x=0.85,
y=0.5,
s=f"Geodesic distance:\n{d_geo.to_compact():.0f}",
transform=ax.transAxes,
c="r",
ha="right",
)
If you want to look at geodesic distances across many nodes, you should consider using navis.geodesic_matrix. To demonstrate, let's generate a geodesic distance matrix between all terminal nodes:
# Calculate distances from all end nodes to all other nodes
ends = n.nodes[n.nodes.type == "end"].node_id.values
m = navis.geodesic_matrix(n, from_=ends)
# Subset to only end-nodes-to-end_nodes
m = m.loc[ends, ends]
m.head()
Let's see if we can visualize any clusters using a seaborn clustermap:
import seaborn as sns
from scipy.cluster.hierarchy import linkage
from scipy.spatial.distance import squareform
# Generate a linkage from the distances
Z = linkage(squareform(m, checks=False), method="ward")
# Plot
cm = sns.clustermap(m, cmap="Greys", col_linkage=Z, row_linkage=Z)
cm.ax_heatmap.set_xticks([])
cm.ax_heatmap.set_yticks([])
Out:
[]
As you can see in the heatmap, the dendrites and the axon nicely separate.
That's it for now! Please see the NBLAST tutorial for morphological comparisons using NBLAST and the API reference for a full list of morphology-related functions.
Total running time of the script: ( 0 minutes 2.417 seconds)
Download Python source code: tutorial_morpho_01_analyze.py