NBLAST

NBLAST

NBLAST (Costa et al., 2016) is a method to quantify morphological similarity. It works on “dotprops” which represent neurons as tangent vectors. For each tangent vector in the query neuron, NBLAST finds the closest tangent vector in the target neuron and calculates a score from the distance between and the dotproduct of the two vectors. The final NBLAST score is the sum over all query-target vector pairs. Typically, this score is normalized to a self-self comparison (i.e. a perfect match would be 1).

VFB computes NBLAST scores for all (?) neurons in its database. So if all you want is a list of similar neurons, it’s fastest (and easiest) to get those directly from VFB.

# Import libs and initialise API objects
from vfb_connect.cross_server_tools import VfbConnect
import navis.interfaces.neuprint as neu

import pandas as pd
import navis

navis.set_pbars(jupyter=False)

vc = VfbConnect()
client = neu.Client('https://neuprint.janelia.org', dataset='hemibrain:v1.1')

First let’s write a little function to grab VFB neurons and turn them into navis neurons. This is modified from vc.neo_query_wrapper.get_images.

import requests
from tqdm import tqdm

def get_vfb_neurons(vfb_ids, template='JRC2018Unisex'):
    """Load neurons from VFB as navis.TreeNeurons."""
    vfb_ids = list(navis.utils.make_iterable(vfb_ids))
    inds = vc.neo_query_wrapper.get_anatomical_individual_TermInfo(short_forms=vfb_ids)
    nl = []
    # Note: queries should be parallelized in the future 
    # -> for pymaid I use request-futures which works really well
    for i in tqdm(inds, desc='Loading VFB neurons', leave=False):        
        if not ('has_image' in i['term']['core']['types']):
            continue
        label = i['term']['core']['label']
        image_matches = [x['image'] for x in i['channel_image']]
        if not image_matches:
            continue
        for imv in image_matches:
            if imv['template_anatomy']['label'] == template:
                r = requests.get(imv['image_folder'] + '/volume.swc')
                ### Slightly dodgy warning - could mask network errors
                if not r.ok:
                    warnings.warn("No '%s' file found for '%s'." % (image_type, label))
                    continue

                # `id` should ideally be unique but that's not enforced
                n = navis.TreeNeuron(r.text, name=label, id=i['term']['core']['short_form'])   

                # This registers attributes that you want to show in the summary
                # alternatively we could have a generic attribute like so:
                # n.template = template
                n._register_attr("template", template, temporary=False)

                # I assume all available templates are in microns?
                # @Robbie -do we have this metadata on templates?  If not, we should add ASAP.
                n.units = '1 micron'

                nl.append(n)
    return navis.NeuronList(nl)

# Search for similar neurons using the VFB ID of an ellipsoid body neuron in FAFB
# We got that ID from a search on the VFB website
similar_to_EPG6L1 = vc.get_similar_neurons('VFB_001012ay')
similar_to_EPG6L1 

# Read those neurons from VFB 
query = get_vfb_neurons('VFB_001012ay')
matches = get_vfb_neurons(similar_to_EPG6L1.id.values)

matches

# Define colors such that query is black 
import seaborn as sns 
colors = ['k'] + sns.color_palette('muted', len(matches))

navis.plot3d([query, matches], color=colors)

So that looks decent!

Next: Let’s say you have a set of neurons and you want to group them into morphologically similar “types”. To demonstrate this, we will take ellipsoid body (protocerebral bridge) neurons from VFB and cluster them using NBLAST.

# Find instances of EB-PB neurons
eb_pb_ids = pd.DataFrame.from_records(vc.get_instances("eb-pb-vbo neuron"))
eb_pb_ids.head()

# Get the neurons from VFB
eb_pb = get_vfb_neurons(eb_pb_ids.id.values)
eb_pb

Now that we have the skeletons, we need to turn them into dotprops. Couple things to keep in mind:

1. NBLAST is optimized for data in microns. Since our neurons are in the JRC2018Unisex template, they already are.
2. Neurons should have the same sampling rate (i.e. nodes per micron). To ensure this we will resample the neurons to 1 point per micron.
3. We won’t be doing this this time around but we typically simplify the neurons slightly by e.g. pruning small twigs (navis.prune_twigs()) to reduce noise.

eb_dps = navis.make_dotprops(eb_pb, resample=1, k=5)

# Note this neuronlist contains dotprops now
eb_dps

To illustrate what dotprops are, lets visualize one side-by-side with its skeleton.

sk = eb_pb[0].copy()
dp = eb_dps[0].copy()

# Slightly offset the dotprops by a micron 
dp.points[:, 0] += 1

navis.plot3d([sk, dp], color=['k', 'r'])

Now we can run the NBLAST. Note that in this case, we run an all-by-all NBLAST which means that we don’t really have to worry about the directionality - i.e. whether we NBLAST A->B or B->A.

scores = navis.nblast_allbyall(eb_dps)
scores.head()

Blasting:   0%|          | 0/62 [00:00<?, ?it/s]
Blasting:   2%|▏         | 1/62 [00:00<00:22,  2.69it/s]
Blasting:   3%|▎         | 2/62 [00:00<00:25,  2.39it/s]
Blasting:   5%|▍         | 3/62 [00:01<00:34,  1.71it/s]
Blasting:   5%|▍         | 3/62 [00:01<00:33,  1.75it/s]
Blasting:   6%|▋         | 4/62 [00:02<00:32,  1.79it/s]
Blasting:   8%|▊         | 5/62 [00:02<00:27,  2.08it/s]
Blasting:  10%|▉         | 6/62 [00:02<00:25,  2.16it/s]
Blasting:  11%|█▏        | 7/62 [00:03<00:24,  2.21it/s]
Blasting:  13%|█▎        | 8/62 [00:03<00:20,  2.67it/s]
Blasting:  13%|█▎        | 8/62 [00:03<00:22,  2.43it/s]
Blasting:  16%|█▌        | 10/62 [00:03<00:13,  3.76it/s]
Blasting:  18%|█▊        | 11/62 [00:04<00:15,  3.31it/s]
Blasting:  19%|█▉        | 12/62 [00:04<00:15,  3.15it/s]
Blasting:  21%|██        | 13/62 [00:04<00:13,  3.51it/s]
Blasting:  24%|██▍       | 15/62 [00:05<00:12,  3.75it/s]
Blasting:  26%|██▌       | 16/62 [00:05<00:14,  3.24it/s]
Blasting:  27%|██▋       | 17/62 [00:06<00:15,  2.83it/s]
Blasting:  29%|██▉       | 18/62 [00:06<00:16,  2.72it/s]
Blasting:  31%|███       | 19/62 [00:07<00:17,  2.44it/s]
Blasting:  32%|███▏      | 20/62 [00:07<00:16,  2.55it/s]
Blasting:  34%|███▍      | 21/62 [00:07<00:16,  2.54it/s]
Blasting:  35%|███▌      | 22/62 [00:08<00:16,  2.44it/s]
Blasting:  37%|███▋      | 23/62 [00:08<00:15,  2.47it/s]
Blasting:  39%|███▊      | 24/62 [00:09<00:16,  2.27it/s]
Blasting:  40%|████      | 25/62 [00:10<00:19,  1.86it/s]
Blasting:  42%|████▏     | 26/62 [00:10<00:19,  1.82it/s]
Blasting:  44%|████▎     | 27/62 [00:11<00:17,  2.00it/s]
Blasting:  45%|████▌     | 28/62 [00:11<00:15,  2.25it/s]
Blasting:  47%|████▋     | 29/62 [00:11<00:14,  2.34it/s]
Blasting:  48%|████▊     | 30/62 [00:12<00:15,  2.08it/s]
Blasting:  48%|████▊     | 30/62 [00:12<00:15,  2.09it/s]
Blasting:  52%|█████▏    | 32/62 [00:12<00:10,  2.80it/s]
Blasting:  53%|█████▎    | 33/62 [00:12<00:08,  3.38it/s]
Blasting:  56%|█████▋    | 35/62 [00:13<00:07,  3.74it/s]
Blasting:  58%|█████▊    | 36/62 [00:13<00:07,  3.44it/s]
Blasting:  60%|█████▉    | 37/62 [00:14<00:07,  3.19it/s]
Blasting:  61%|██████▏   | 38/62 [00:14<00:06,  3.54it/s]
Blasting:  63%|██████▎   | 39/62 [00:14<00:05,  3.90it/s]
Blasting:  65%|██████▍   | 40/62 [00:14<00:06,  3.40it/s]
Blasting:  66%|██████▌   | 41/62 [00:15<00:06,  3.07it/s]
Blasting:  68%|██████▊   | 42/62 [00:15<00:06,  3.19it/s]
Blasting:  69%|██████▉   | 43/62 [00:15<00:05,  3.24it/s]
Blasting:  73%|███████▎  | 45/62 [00:16<00:04,  3.84it/s]
Blasting:  74%|███████▍  | 46/62 [00:16<00:04,  3.57it/s]
Blasting:  76%|███████▌  | 47/62 [00:17<00:05,  2.85it/s]
Blasting:  77%|███████▋  | 48/62 [00:17<00:04,  2.93it/s]
Blasting:  79%|███████▉  | 49/62 [00:17<00:04,  2.67it/s]
Blasting:  81%|████████  | 50/62 [00:18<00:04,  2.64it/s]
Blasting:  82%|████████▏ | 51/62 [00:18<00:04,  2.62it/s]
Blasting:  84%|████████▍ | 52/62 [00:18<00:03,  2.65it/s]
Blasting:  85%|████████▌ | 53/62 [00:19<00:03,  2.66it/s]
Blasting:  87%|████████▋ | 54/62 [00:19<00:03,  2.30it/s]
Blasting:  87%|████████▋ | 54/62 [00:20<00:03,  2.18it/s]
Blasting:  90%|█████████ | 56/62 [00:21<00:02,  2.19it/s]
Blasting:  92%|█████████▏| 57/62 [00:21<00:02,  2.14it/s]
Blasting:  94%|█████████▎| 58/62 [00:22<00:02,  1.93it/s]
Blasting:  95%|█████████▌| 59/62 [00:22<00:01,  1.87it/s]
Blasting:  95%|█████████▌| 59/62 [00:22<00:01,  1.69it/s]
Blasting:  98%|█████████▊| 61/62 [00:23<00:00,  2.01it/s]
Blasting: 100%|██████████| 62/62 [00:23<00:00,  2.27it/s]

Next, we need to take those forward scores (A->B) and symmetrize by calculating mean scores (A<->B).

# Symmetrize scores
scores_mean = (scores + scores.T) / 2

scores_mean.head()

From here on out, we use cookie cutter scipy to produce a hierarchical clustering. First we need to take those similarity scores to distances because that’s what scipy operates on:

dist_mean = (scores_mean - 1) * -1
dist_mean.head()

from scipy.spatial.distance import squareform
from scipy.cluster.hierarchy import linkage, dendrogram, fcluster

# This turns the square distance matrix into 1-d vector form
dist_sq = squareform(dist_mean)

# This does the actual hierarchical clustering
Z = linkage(dist_sq, method='ward', optimal_ordering=True)

# Plot a quick dendrogram
import matplotlib.pyplot as plt 
import seaborn as sns 

# Generate the canvas
fig, ax = plt.subplots(figsize=(12, 5))

# Make some more meaningful labels 
vfb2name = dict(zip(eb_pb.id, eb_pb.name))
labels = dist_mean.index.map(lambda x: vfb2name[x].split(' - ')[0])

# Plot the dendrogram
dn = dendrogram(Z, ax=ax, labels=labels, color_threshold=.65, above_threshold_color='lightgrey')

# Make xticks biggger
ax.set_xticklabels(ax.get_xticklabels(), size=11)

# Clean axes
sns.despine(bottom=True)

# Make cluster matching the color threshold in the dendrogram
clusters = fcluster(Z, t=.65, criterion='distance')
clusters

array([16, 14,  4,  5, 11,  4,  3, 12,  2, 19, 13,  4,  8,  3, 19,  2,  5,
       13,  3, 11, 17, 14, 16,  6, 12, 15, 18, 11, 12, 15, 15, 18,  5, 10,
        9,  2, 19,  8,  1,  9, 19,  3,  1, 10,  2,  9,  7, 18,  9,  2, 17,
       13, 12,  4, 17,  5,  7,  6,  4, 12, 16,  4], dtype=int32)

# We have to downsample the neurons for visualization 
# because DeepNote enforces a limit on the size of plots
eb_pb_ds = eb_pb.downsample(10)

# Make a color per cluster 
palette = sns.color_palette('tab20', max(clusters))
cmap = dict(zip(range(1, max(clusters) + 1), palette))

colors = [cmap[i] for i in clusters]

navis.plot3d(eb_pb_ds, color=colors)