.. _nblast_fc:

NBLAST against FlyCircuit
-------------------------

If you work with `Drosophila`, chances are you have heard of `FlyCircuit <http://www.flycircuit.tw/>`_. It's a collection of ~16k single neuron clones published by `Chiang et al. (2010) <https://www.cell.com/current-biology/fulltext/S0960-9822(10)01522-8>`_.

For R, there is a package containing dotprops and meta data for these neurons: `nat.flycircuit <https://github.com/natverse/flycircuit>`_. This does not yet exist for Python. However, it's still reasonably straight forward to run an NBLAST against the entire flycircuit dataset in Python:

First you need to download flycircuit dotprops from `Zenodo <https://zenodo.org/record/5205616>`_ (850Mb). Next, unzip the archive containing the dotprops as CSV files. Now we can load them into navis (adjust filepaths as required):

.. code:: ipython3

    import navis 
    import pandas as pd
    
    from pathlib import Path
    from tqdm import tqdm 

.. code:: ipython3

    def load_dotprops_csv(fp):
        """Load dotprops from CSV files in filepath `fp`."""
        # Turn into a path object 
        fp = Path(fp).expanduser()
        
        # Go over each CSV file
        files = list(fp.glob('*.csv'))     
        dotprops = []
        for f in tqdm(files, desc='Loading dotprops', leave=False):
            # Read file 
            csv = pd.read_csv(f)
            
            # Each row is a point with associated vector 
            pts = csv[['pt_x', 'pt_y', 'pt_z']].values
            vect = csv[['vec_x', 'vec_y', 'vec_z']].values
            
            # Turn into a Dotprops        
            dp = navis.Dotprops(points=pts, k=20, vect=vect, units='1 micron')
            
            # Use filename as ID/name
            dp.name = dp.id = f.name[:-4]
            
            # Add this dotprop to the list before moving on to the next
            dotprops.append(dp)
        
        return navis.NeuronList(dotprops)
        
    fc_dps = load_dotprops_csv('~/Downloads/flycircuit_dotprops')
    fc_dps


.. parsed-literal::

                                                                            



.. raw:: html

    &lt;class 'navis.core.neuronlist.NeuronList'&gt; containing 16129 neurons (est. 1.4GiB)<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>k</th>
          <th>units</th>
          <th>n_points</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>navis.Dotprops</td>
          <td>DvGlutMARCM-F1822_seg2</td>
          <td>DvGlutMARCM-F1822_seg2</td>
          <td>20</td>
          <td>1 micrometer</td>
          <td>87</td>
        </tr>
        <tr>
          <th>1</th>
          <td>navis.Dotprops</td>
          <td>DvGlutMARCM-F002216_seg001</td>
          <td>DvGlutMARCM-F002216_seg001</td>
          <td>20</td>
          <td>1 micrometer</td>
          <td>2574</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>16127</th>
          <td>navis.Dotprops</td>
          <td>GadMARCM-F000020_seg001</td>
          <td>GadMARCM-F000020_seg001</td>
          <td>20</td>
          <td>1 micrometer</td>
          <td>477</td>
        </tr>
        <tr>
          <th>16128</th>
          <td>navis.Dotprops</td>
          <td>DvGlutMARCM-F1720_seg1</td>
          <td>DvGlutMARCM-F1720_seg1</td>
          <td>20</td>
          <td>1 micrometer</td>
          <td>941</td>
        </tr>
      </tbody>
    </table>
    </div>



To avoid having to re-generate these dotprops, you could consider pickling them:

.. highlight::python

    import pickle
    with open('flycircuit_dps.pkl', 'wb') as f:
        pickle.dump(fc_dps, f)
        
In the future you can then reload the dotprops like so (much faster):

.. highlight::python

    with open('flycircuit_dps.pkl', 'rb') as f:
        fc_dps = pickle.load(f)

Note: The names/ids correspond to the flycircuit gene + clone names!

.. code:: ipython3

    fc_dps[0]




.. 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>FruMARCM-M002262_seg001</td>
        </tr>
        <tr>
          <th>id</th>
          <td>FruMARCM-M002262_seg001</td>
        </tr>
        <tr>
          <th>k</th>
          <td>5</td>
        </tr>
        <tr>
          <th>units</th>
          <td>1 micrometer</td>
        </tr>
        <tr>
          <th>n_points</th>
          <td>4121</td>
        </tr>
      </tbody>
    </table>
    </div>



In case your query neurons are in a different brain space, you can use `navis-flybrains <https://github.com/navis-org/navis-flybrains>`_ to convert them to flycircuit's ``FCWB`` space.

For demonstration purposes we will use the example neurons - olfactory DA1 projection neurons from the hemibrain connectome - that ship with ``navis``.

.. code:: ipython3

    # Load some of the example neurons
    n = navis.example_neurons(3)

.. code:: ipython3

    # Convert from hemibrain (JRCFIB2018Fraw) to FCWB space
    import flybrains
    
    n_fcwb = navis.xform_brain(n, source='JRCFIB2018Fraw', target='FCWB')


.. parsed-literal::

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


.. code:: ipython3

    # A sanity check to make sure the transform worked
    fig, ax = navis.plot2d([n_fcwb, flybrains.FCWB])
    
    # Rotate to front view
    ax.elev = -90
    
    # Zoom in a little
    ax.dist = 4







.. image:: nblast_flycircuit_files/nblast_flycircuit_9_1.png


.. code:: ipython3

    # Mirror neurons from the right to the left
    # (FlyCircuit neurons are supposed to all be on the left)
    n_fcwb_mirr = navis.mirror_brain(n_fcwb, template='FCWB')

.. code:: ipython3

    # Convert to dotprops
    n_dps = navis.make_dotprops(n_fcwb_mirr, resample=1, k=None)

.. code:: ipython3

    # Run a "smart" NBLAST to get the top hits
    # (ignore the warning about data not being in microns - that's a rounding error from the transform)
    scores = navis.nblast_smart(query=n_dps, target=fc_dps, scores='mean', progress=False)

.. 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>FruMARCM-M002262_seg001</th>
          <th>FruMARCM-M002262_seg002</th>
          <th>FruMARCM-M002278_seg001</th>
          <th>FruMARCM-M002287_seg001</th>
          <th>FruMARCM-M002507_seg001</th>
          <th>FruMARCM-M002522_seg001</th>
          <th>FruMARCM-M002578_seg001</th>
          <th>FruMARCM-M002578_seg002</th>
          <th>FruMARCM-M002579_seg001</th>
          <th>FruMARCM-M002579_seg002</th>
          <th>...</th>
          <th>DvGlutMARCM-F594-X2_seg1</th>
          <th>DvGlutMARCM-F594-X2_seg2</th>
          <th>DvGlutMARCM-F595_seg1</th>
          <th>DvGlutMARCM-F596_seg1</th>
          <th>DvGlutMARCM-F597_seg1</th>
          <th>DvGlutMARCM-F598_seg1</th>
          <th>DvGlutMARCM-F599_seg1</th>
          <th>DvGlutMARCM-F599_seg2</th>
          <th>DvGlutMARCM-F600-X2_seg1</th>
          <th>DvGlutMARCM-F600-X2_seg2</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>1734350788</th>
          <td>-0.260213</td>
          <td>-0.200286</td>
          <td>-0.645926</td>
          <td>-0.639277</td>
          <td>-0.663507</td>
          <td>-0.451059</td>
          <td>-0.229996</td>
          <td>-0.604237</td>
          <td>-0.705247</td>
          <td>-0.307929</td>
          <td>...</td>
          <td>-0.881092</td>
          <td>-0.648820</td>
          <td>-0.660410</td>
          <td>0.118524</td>
          <td>-0.315490</td>
          <td>-0.671891</td>
          <td>-0.499438</td>
          <td>-0.882007</td>
          <td>-0.882510</td>
          <td>-0.880825</td>
        </tr>
        <tr>
          <th>1734350908</th>
          <td>-0.267769</td>
          <td>-0.191842</td>
          <td>-0.655194</td>
          <td>-0.623617</td>
          <td>-0.679471</td>
          <td>-0.459428</td>
          <td>-0.233066</td>
          <td>-0.607806</td>
          <td>-0.715372</td>
          <td>-0.310166</td>
          <td>...</td>
          <td>-0.881365</td>
          <td>-0.654605</td>
          <td>-0.669929</td>
          <td>0.113844</td>
          <td>-0.316001</td>
          <td>-0.682929</td>
          <td>-0.505106</td>
          <td>-0.881421</td>
          <td>-0.881361</td>
          <td>-0.882250</td>
        </tr>
        <tr>
          <th>722817260</th>
          <td>-0.270642</td>
          <td>-0.170228</td>
          <td>-0.680244</td>
          <td>-0.667571</td>
          <td>-0.670850</td>
          <td>-0.431723</td>
          <td>-0.236943</td>
          <td>-0.634743</td>
          <td>-0.704928</td>
          <td>-0.298078</td>
          <td>...</td>
          <td>-0.881398</td>
          <td>-0.649228</td>
          <td>-0.659683</td>
          <td>0.132546</td>
          <td>-0.317336</td>
          <td>-0.627997</td>
          <td>-0.503828</td>
          <td>-0.882124</td>
          <td>-0.881570</td>
          <td>-0.881419</td>
        </tr>
      </tbody>
    </table>
    <p>3 rows × 16129 columns</p>
    </div>



.. code:: ipython3

    # Find the top hits for each of our query neurons
    import numpy as np
    
    for dp in n_dps:
        hit = scores.columns[np.argmax(scores.loc[dp.id].values)]
        sc = scores.loc[dp.id].values.max()
        print(f'Top hit for {dp.id}: {hit} ({sc:.3f})')    


.. parsed-literal::

    Top hit for 1734350788: FruMARCM-F001496_seg001 (0.778)
    Top hit for 1734350908: FruMARCM-F001496_seg001 (0.744)
    Top hit for 722817260: FruMARCM-F001496_seg001 (0.792)


Conveniently all of our query neurons have the same top match (they are all from the same cell type after all): ``FruMARCM-F001496_seg001``

.. code:: ipython3

    # Let's co-visualize:
    # Queries in red, hit in black
    fig, ax = navis.plot2d([n_fcwb_mirr, fc_dps.idx['FruMARCM-F001496_seg001']],
                           color=['r'] * len(n_fcwb_mirr) + ['k'])
    ax.elev = -90







.. image:: nblast_flycircuit_files/nblast_flycircuit_16_1.png


Looking good!