.. _nblast_hemibrain:

NBLAST using light-level data
*****************************

One of the applications of NBLAST is to match neurons across data sets. Here we will illustrate this by taking a light-level, confocal image stack and finding the same neuron in an EM connectome. 

Specifically, we will use an image from Janelia's collection of split-Gal4 driver lines and match it against the neurons in the `hemibrain connectome <https://neuprint.janelia.org>`_. 

Before we get started make sure that:
- ``navis`` is installed and up-to-date 
- `navis-flybrains <https://github.com/navis-org/navis-flybrains>` is installed and you have downloaded the Saalfeld lab's and VFB bridging transforms (see ``flybrains.download_...`` functions)
- download and extract `hemibrain-v1.2-skeletons.tar <https://storage.googleapis.com/hemibrain/v1.2/hemibrain-v1.2-skeletons.tar.gz>`_  (kindly provided by Stuart Berg, Janelia)

Next we need to pick an image stack to use as query. You can browse the expression patterns of the Janelia split-Gal4 lines `here <https://splitgal4.janelia.org/cgi-bin/splitgal4.cgi>`_. Here I picked ``LH1112`` which is a very clean line containing a couple of WED projection neurons. Among other data, you can download these stacks as "gendered" (i.e. female or male) or "unisex" space. Unfortunately, all image stacks are in Janelia's ``.h5j`` format which I haven't figured out how to import straight from Python. Two options:

1. Load them into Fiji and save the GFP signal channel as `.nrrd` file.
2. Go to `VirtualFlyBrain <http://www.virtualflybrain.org/>`_, search for your line of interested LH1112 (not all lines are be available on VFB) and download the "Signal(NRRD)" at bottom of Term Info panel. 

I went for option two here and downloaded a ``VFB_001013cg.nrrd``. This is the neuron we'll be looking for:

.. image:: https://s3.amazonaws.com/janelia-flylight-imagery/Lateral+Horn+2019/LH1112/LH1112-20150313_46_A2-f-20x-brain-Split_GAL4-multichannel_mip.png
   :width: 500
   :align: center
   
Let's get started!

.. code:: ipython3

    import navis

First we need to load the image stack and turn it into dotprops:

.. code:: ipython3

    query = navis.read_nrrd('VFB_001013cg.nrrd', output='dotprops', k=20, threshold=100)
    query.id = 'LH1112'  # manually set the ID to the Janelia identifier
    query




.. 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></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>type</th>
          <td>navis.Dotprops</td>
        </tr>
        <tr>
          <th>name</th>
          <td>VFB_001013cg</td>
        </tr>
        <tr>
          <th>id</th>
          <td>LH1112</td>
        </tr>
        <tr>
          <th>k</th>
          <td>20</td>
        </tr>
        <tr>
          <th>units</th>
          <td>1 micrometer</td>
        </tr>
        <tr>
          <th>n_points</th>
          <td>10375</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig = navis.plot3d(query)

.. raw:: html
   :file: figures/3d_nblast_dotprops.html

|
|
|
|

.. note::
    
   Depending on your image you will have to play around with the threshold to get a decent representation.
   
Next we need to load the hemibrain skeletons and convert them to dotprops:

.. code:: ipython3

    sk = navis.read_swc('/Users/philipps/Downloads/hemibrain-v1.2-skeletons/',
                        include_subdirs=True, fmt='{name,id:int}.swc')






These 97k skeletons include lots of small fragments - there are only ~25k proper neurons or substantial fragments thereof in the hemibrain dataset. So to make our life a little easier, we will keep only the largest 30k skeletons:

.. code:: ipython3

    sk.sort_values('n_nodes')
    sk = sk[:30_000]
    sk




.. raw:: html

    &lt;class 'navis.core.neuronlist.NeuronList'&gt; containing 30000 neurons (est. 3.2GiB)<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>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>425790257</td>
          <td>182514</td>
          <td>None</td>
          <td>20286</td>
          <td>20851</td>
          <td>1.001564e+07</td>
          <td>None</td>
          <td>1 dimensionless</td>
        </tr>
        <tr>
          <th>1</th>
          <td>navis.TreeNeuron</td>
          <td>5813105172</td>
          <td>159749</td>
          <td>None</td>
          <td>22416</td>
          <td>22866</td>
          <td>7.610828e+06</td>
          <td>None</td>
          <td>1 dimensionless</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>29998</th>
          <td>navis.TreeNeuron</td>
          <td>5813035385</td>
          <td>1018</td>
          <td>None</td>
          <td>90</td>
          <td>100</td>
          <td>6.213784e+04</td>
          <td>None</td>
          <td>1 dimensionless</td>
        </tr>
        <tr>
          <th>29999</th>
          <td>navis.TreeNeuron</td>
          <td>2375870447</td>
          <td>1018</td>
          <td>None</td>
          <td>151</td>
          <td>156</td>
          <td>3.505048e+04</td>
          <td>None</td>
          <td>1 dimensionless</td>
        </tr>
      </tbody>
    </table>
    </div>



Next up: turning those skeletons into dotprops:

.. code:: ipython3

    # Note that we are resampling to 1 point for every micron of cable 
    # Because these skeletons are in 8nm voxels we have to use 1000/8
    targets = navis.make_dotprops(sk, k=5, parallel=True, resample=1000/8)






.. note:
    
   Making the dotprops may take a while (mostly because of the resampling). You can dedicate more
   cores via the ``n_cores`` parameter. It may also make sense to save the dotprops for future
   use e.g. by pickling them.
   
Last but not least we need to convert the image's dotprops from their current brain space (``JRC2018U``, ``U`` for "unisex") to hemibrain (``JRCFIB2018Fraw``, ``raw`` for voxels) space.

.. code:: ipython3

    import flybrains
    
    query_xf = navis.xform_brain(query, source='JRC2018U', target='JRCFIB2018Fraw')


.. parsed-literal::

    Transform path: JRC2018U -> JRC2018F -> JRCFIB2018Fum -> JRCFIB2018F -> JRCFIB2018Fraw


Now we can run the NBLAST:

.. code:: ipython3

    # Note that we convert from voxels (8x8x8nm) to microns
    scores = navis.nblast(query_xf / 125, targets / 125, scores='mean')

.. code:: ipython3

    scores.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>425790257</th>
          <th>5813105172</th>
          <th>5813057579</th>
          <th>5813050499</th>
          <th>859839499</th>
          <th>5813039148</th>
          <th>1640909284</th>
          <th>423101189</th>
          <th>1321564092</th>
          <th>5813024698</th>
          <th>...</th>
          <th>2032311172</th>
          <th>2283477509</th>
          <th>1686751696</th>
          <th>2180970798</th>
          <th>1904079437</th>
          <th>1348510839</th>
          <th>5813047115</th>
          <th>516672812</th>
          <th>5813035385</th>
          <th>2375870447</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>LH1112</th>
          <td>-0.793901</td>
          <td>-0.870922</td>
          <td>-0.767066</td>
          <td>-0.627762</td>
          <td>-0.612109</td>
          <td>-0.591567</td>
          <td>-0.827376</td>
          <td>-0.878581</td>
          <td>-0.680821</td>
          <td>-0.832845</td>
          <td>...</td>
          <td>-0.881125</td>
          <td>-0.817378</td>
          <td>-0.88186</td>
          <td>-0.881415</td>
          <td>-0.881596</td>
          <td>-0.549413</td>
          <td>-0.881259</td>
          <td>-0.880987</td>
          <td>-0.88139</td>
          <td>-0.667554</td>
        </tr>
      </tbody>
    </table>
    <p>1 rows × 30000 columns</p>
    </div>



.. note::
    
   You can greatly speed up NBLAST by installing the optional dependency ``pykdtree``:
   ``pip3 install pykdtree``.
   
Now we can sort those scores to find the top matches:

.. code:: ipython3

    scores.loc['LH1112'].sort_values(ascending=False)




.. parsed-literal::

    2030342003    0.218951
    2214504597    0.214873
    1069223047    0.183628
    856131667     0.119994
    2214846055    0.094881
                    ...   
    2056511651   -0.883308
    296501634    -0.883329
    2275102579   -0.883354
    5813091313   -0.883436
    910404284    -0.883766
    Name: LH1112, Length: 30000, dtype: float64



So we did get a couple of hits here. Let's visualize the top 3:

.. code:: ipython3

    fig = navis.plot3d([query_xf, targets.idx[[2030342003, 2214504597, 1069223047]]])

.. raw:: html
   :file: figures/3d_nblast_results2.html

|
|
|
|

On a final note: the scores for those matches are rather low (1 = perfect match). 

The main reason for this is that our image stack contains two neurons (the left and the right version) so half of our ``query`` won't have a match in any of the individual hemibrain neurons. We could have fixed that by subsetting the query to the approximate hemibrain bounding box. This is also a good idea for bilateral neurons that have parts of their arbour outside the hemibrain volume:

.. code:: ipython3

    # Remove the left-hand-side neuron based on the position 
    # along the x-axis (this is just one of the possible approaches)
    # This is the approximate LHS boundary of the volume
    flybrains.JRCFIB2018Fraw.mesh.vertices[:, 0].max()




.. parsed-literal::

    TrackedArray(34499.97265625)



.. code:: ipython3

    query_ss = navis.subset_neuron(query_xf, query_xf.points[:, 0] <= 35_000)
    query_ss




.. 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></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>type</th>
          <td>navis.Dotprops</td>
        </tr>
        <tr>
          <th>name</th>
          <td>VFB_001013cg</td>
        </tr>
        <tr>
          <th>k</th>
          <td>20</td>
        </tr>
        <tr>
          <th>units</th>
          <td>9.999999999999998 nanometer</td>
        </tr>
        <tr>
          <th>n_points</th>
          <td>5093</td>
        </tr>
      </tbody>
    </table>
    </div>



Using ``query_ss`` should yield much improved scores. 

Another potential pitfall is the generation of dotprops from the image itself: if you compare the image- against the skeleton-derived dotprops, you might notice that the latter have fewer and less dense points. That's a natural consequence the image containing multiple individuals of the same cell type but we could have tried to ameliorate this by some pre-processing (e.g. downsampling or thinning the image).