This tutorial will give you an overview of how to analyze neuron morphology. See the :ref:`API reference <api_morph>` for a complete list of available functions/methods.

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. :func:`~navis.tortuosity`) - mostly because they require/allow some parameters to be set.

With that in mind, let's dive right in:

.. code:: ipython3

    import navis
    import matplotlib.pyplot as plt
    
    import seaborn as sns
    sns.set_color_codes('muted')

Accessing attributes for a single neuron:

.. code:: ipython3

    # Load a single skeleton (i.e. a TreeNeuron)
    n = navis.example_neurons(n=1, kind='skeleton')
    n.cable_length




.. parsed-literal::

    266476.88



If your neuron has its ``units`` properties set, you can also combine those:

.. code:: ipython3

    print((n.cable_length * n.units).to('microns'))


.. parsed-literal::

    2131.815 micron


Accessing the same properties via a :class:`~navis.NeuronList` will return an array:

.. code:: ipython3

    # Load a couple skeletons
    nl = navis.example_neurons(kind='skeleton')
    print(nl.cable_length)


.. parsed-literal::

    [266476.88 304332.66 274703.38 286522.47 291265.3 ]


.. code:: ipython3

    print(nl.n_branches)


.. parsed-literal::

    [599 735 633 696 626]


:meth:`navis.NeuronList.summary` is a useful way to collect some of the basic parameters:

.. code:: ipython3

    df = nl.summary()
    df.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>type</th>
          <th>name</th>
          <th>id</th>
          <th>n_nodes</th>
          <th>n_connectors</th>
          <th>n_branches</th>
          <th>n_leafs</th>
          <th>cable_length</th>
          <th>soma</th>
          <th>units</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>navis.TreeNeuron</td>
          <td>1734350788</td>
          <td>1734350788</td>
          <td>4465</td>
          <td>2705</td>
          <td>599</td>
          <td>618</td>
          <td>266476.87500</td>
          <td>[4177]</td>
          <td>8 nanometer</td>
        </tr>
        <tr>
          <th>1</th>
          <td>navis.TreeNeuron</td>
          <td>1734350908</td>
          <td>1734350908</td>
          <td>4847</td>
          <td>3042</td>
          <td>735</td>
          <td>761</td>
          <td>304332.65625</td>
          <td>[6]</td>
          <td>8 nanometer</td>
        </tr>
        <tr>
          <th>2</th>
          <td>navis.TreeNeuron</td>
          <td>722817260</td>
          <td>722817260</td>
          <td>4332</td>
          <td>3136</td>
          <td>633</td>
          <td>656</td>
          <td>274703.37500</td>
          <td>None</td>
          <td>8 nanometer</td>
        </tr>
        <tr>
          <th>3</th>
          <td>navis.TreeNeuron</td>
          <td>754534424</td>
          <td>754534424</td>
          <td>4696</td>
          <td>3010</td>
          <td>696</td>
          <td>726</td>
          <td>286522.46875</td>
          <td>[4]</td>
          <td>8 nanometer</td>
        </tr>
        <tr>
          <th>4</th>
          <td>navis.TreeNeuron</td>
          <td>754538881</td>
          <td>754538881</td>
          <td>4881</td>
          <td>2943</td>
          <td>626</td>
          <td>642</td>
          <td>291265.31250</td>
          <td>[701]</td>
          <td>8 nanometer</td>
        </tr>
      </tbody>
    </table>
    </div>



For ``MeshNeurons`` the available properties look different. For example, you can get its volume:

.. code:: ipython3

    # Load a single MeshNeuron
    m = navis.example_neurons(n=1, kind='mesh')
    print(m.volume)


.. parsed-literal::

    1291610825.1668377


For topological properties, we need to convert to skeleton. The fastest way is to simply access the ``MeshNeuron``'s ``.skeleton`` property:

.. code:: ipython3

    print(m.skeleton.cable_length)


.. parsed-literal::

    168890.5891599271


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 :func:`navis.skeletonize`.

Importantly, some ``navis`` functions (e.g. :func:`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:

.. code:: ipython3

    m.skeleton.n_leafs




.. parsed-literal::

    155



.. code:: ipython3

    m.skeleton.n_branches




.. parsed-literal::

    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:

.. code:: ipython3

    m.n_vertices




.. parsed-literal::

    6309



.. code:: ipython3

    m.skeleton.n_nodes




.. parsed-literal::

    1058



.. code:: ipython3

    # Illustrate with a random property
    m.my_counts = [1, 2, 3, 4, 5]
    m.n_my_counts




.. parsed-literal::

    5



Segment Analysis
================

:func:`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:

.. code:: ipython3

    sa = navis.segment_analysis(m)
    sa.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>length</th>
          <th>tortuosity</th>
          <th>root_dist</th>
          <th>strahler_index</th>
          <th>radius_mean</th>
          <th>radius_min</th>
          <th>radius_max</th>
          <th>volume</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>84.610286</td>
          <td>1.320717</td>
          <td>7083.656132</td>
          <td>1</td>
          <td>61.481414</td>
          <td>11.867318</td>
          <td>148.842288</td>
          <td>1.420165e+06</td>
        </tr>
        <tr>
          <th>1</th>
          <td>276.524917</td>
          <td>2.589908</td>
          <td>41381.328436</td>
          <td>1</td>
          <td>58.348518</td>
          <td>3.699802</td>
          <td>165.413996</td>
          <td>6.111648e+06</td>
        </tr>
        <tr>
          <th>2</th>
          <td>141.448954</td>
          <td>1.079153</td>
          <td>12560.860744</td>
          <td>1</td>
          <td>86.103133</td>
          <td>24.133240</td>
          <td>156.117439</td>
          <td>4.873987e+06</td>
        </tr>
        <tr>
          <th>3</th>
          <td>59.935134</td>
          <td>1.122367</td>
          <td>11144.594281</td>
          <td>1</td>
          <td>32.241229</td>
          <td>3.699802</td>
          <td>85.624280</td>
          <td>4.333826e+05</td>
        </tr>
        <tr>
          <th>4</th>
          <td>141.513383</td>
          <td>1.362527</td>
          <td>38527.170130</td>
          <td>1</td>
          <td>44.703457</td>
          <td>9.421692</td>
          <td>139.731548</td>
          <td>1.472218e+06</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    # See if segment length correlates with radius
    ax = sns.scatterplot(data=sa, x='length', y='radius_mean', size='strahler_index', alpha=.7)
    ax.set_xscale('log')
    ax.set_yscale('log')



.. image:: morph_analysis_files/morph_analysis_24_0.png


Sholl Analysis
==============

For an example of a Sholl analyses, check out the :ref:`MICrONS tutorial<sholl>`.

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

.. code:: ipython3

    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 Euclidian 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})')


.. parsed-literal::

    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, :func:`~navis.dist_between` should be sufficient:

.. code:: ipython3

    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}')


.. parsed-literal::

    Euclidean distance between soma and terminal node 465: 444096.17 nanometer


Let's visualise this:

.. code:: ipython3

    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[0], 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 Euclidian 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[0], ['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=.3, s=f'Euclidean distance:\n{d_eucl.to_compact():.0f}', transform=ax.transAxes, c='g')
    _ = ax.text(x=.85, y=.5, s=f'Geodesic distance:\n{d_geo.to_compact():.0f}', transform=ax.transAxes, c='r', ha='right')
    
    plt.show()



.. image:: morph_analysis_files/morph_analysis_30_0.png


If you want to look at geodesic distances across many nodes, you should consider using :func:`~navis.geodesic_matrix`. To demonstrate, let's generate a geodesic distance matrix between all terminal nodes:

.. code:: ipython3

    # 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()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>465</th>
          <th>548</th>
          <th>618</th>
          <th>683</th>
          <th>745</th>
          <th>789</th>
          <th>832</th>
          <th>872</th>
          <th>911</th>
          <th>949</th>
          <th>...</th>
          <th>4456</th>
          <th>4457</th>
          <th>4458</th>
          <th>4459</th>
          <th>4460</th>
          <th>4461</th>
          <th>4462</th>
          <th>4463</th>
          <th>4464</th>
          <th>4465</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>465</th>
          <td>0.000000</td>
          <td>54395.706752</td>
          <td>53556.405732</td>
          <td>54489.728022</td>
          <td>53685.768456</td>
          <td>52679.982779</td>
          <td>53139.882271</td>
          <td>20944.502247</td>
          <td>53065.272557</td>
          <td>53873.631365</td>
          <td>...</td>
          <td>50402.960991</td>
          <td>50363.976517</td>
          <td>49936.090104</td>
          <td>49895.070897</td>
          <td>55185.275744</td>
          <td>54742.764466</td>
          <td>55158.505486</td>
          <td>55950.521492</td>
          <td>56111.790505</td>
          <td>55266.165581</td>
        </tr>
        <tr>
          <th>548</th>
          <td>54395.706752</td>
          <td>0.000000</td>
          <td>8980.969668</td>
          <td>10787.026090</td>
          <td>6228.693973</td>
          <td>6866.540223</td>
          <td>8564.446208</td>
          <td>37619.898765</td>
          <td>8489.836493</td>
          <td>6416.556882</td>
          <td>...</td>
          <td>4589.518435</td>
          <td>5788.540454</td>
          <td>5360.654040</td>
          <td>5319.634833</td>
          <td>12127.772206</td>
          <td>11685.260928</td>
          <td>12101.001948</td>
          <td>12893.017955</td>
          <td>13054.286967</td>
          <td>12208.662044</td>
        </tr>
        <tr>
          <th>618</th>
          <td>53556.405732</td>
          <td>8980.969668</td>
          <td>0.000000</td>
          <td>9947.725069</td>
          <td>8271.031372</td>
          <td>7265.245695</td>
          <td>4698.465309</td>
          <td>36780.597744</td>
          <td>3951.850620</td>
          <td>8458.894281</td>
          <td>...</td>
          <td>4988.223907</td>
          <td>4060.795593</td>
          <td>3828.188293</td>
          <td>3866.356197</td>
          <td>11288.471186</td>
          <td>10845.959908</td>
          <td>11261.700928</td>
          <td>12053.716934</td>
          <td>12214.985947</td>
          <td>11369.361023</td>
        </tr>
        <tr>
          <th>683</th>
          <td>54489.728022</td>
          <td>10787.026090</td>
          <td>9947.725069</td>
          <td>0.000000</td>
          <td>10077.087793</td>
          <td>9071.302116</td>
          <td>9531.201609</td>
          <td>37713.920034</td>
          <td>9456.591894</td>
          <td>10264.950703</td>
          <td>...</td>
          <td>6794.280329</td>
          <td>6755.295855</td>
          <td>6327.409441</td>
          <td>6286.390234</td>
          <td>12221.793476</td>
          <td>11779.282198</td>
          <td>12195.023218</td>
          <td>12987.039225</td>
          <td>13148.308237</td>
          <td>12302.683313</td>
        </tr>
        <tr>
          <th>745</th>
          <td>53685.768456</td>
          <td>6228.693973</td>
          <td>8271.031372</td>
          <td>10077.087793</td>
          <td>0.000000</td>
          <td>6156.601927</td>
          <td>7854.507912</td>
          <td>36909.960468</td>
          <td>7779.898197</td>
          <td>4485.482889</td>
          <td>...</td>
          <td>3879.580139</td>
          <td>5078.602158</td>
          <td>4650.715744</td>
          <td>4609.696537</td>
          <td>11417.833910</td>
          <td>10975.322632</td>
          <td>11391.063652</td>
          <td>12183.079659</td>
          <td>12344.348671</td>
          <td>11498.723747</td>
        </tr>
      </tbody>
    </table>
    <p>5 rows × 618 columns</p>
    </div>



Let's see if we can visualize any clusters using a ``seaborn`` clustermap:

.. code:: ipython3

    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([])
    
    plt.show()



.. image:: morph_analysis_files/morph_analysis_34_0.png


As you can see in the heatmap, the dendrites and the axon nicely separate.

That's it for now! Please see the :ref:`NBLAST tutorial <nblast_intro>` for morphological comparisons using NBLAST and the :ref:`API<api_morph>` for a full list of morphology-related functions.