Brain Region Connectivity function

[stuarteberg]

The website allows you to ask for aggregate connectivity between ROIs, which (presumably) it extracts from the cached connectivity matrix the server constructed for the front page.

It would be nice to include a python function that can compute that info.

Implementation here:
https://flyemteam.slack.com/archives/CA42GNHUL/p1699586516438849

And here

from neuprint import Client, fetch_neurons, NeuronCriteria as NC
c = Client('neuprint.janelia.org', 'hemibrain:v1.2.1')

# Select the ROIs of interest.
# Right now, we default to all "primary" ROIs.
# Make sure this is sorted!
heatmap_rois = sorted(c.primary_rois)

# Fetch each Neuron's ROI-wise synapse distribution.
_, syndist = fetch_neurons(NC())

# Filter out neurons which don't touch our selected ROIs
syndist = syndist.query('roi in @heatmap_rois')

# Drop neurons which don't have any inputs at all.
has_inputs = syndist.groupby('bodyId')['upstream'].transform(np.any)
syndist = syndist.loc[has_inputs]

# If we eliminated any ROIs entirely, drop them from our list.
heatmap_rois = sorted(pd.unique(syndist['roi']))

# For each neuron, multiply (pairwise) all upstream fractions by all downstream counts.
upstream = syndist.set_index(['bodyId', 'roi'])['upstream'].unstack(-1, 0.0)
downstream = syndist.set_index(['bodyId', 'roi'])['downstream'].unstack(-1, 0.0)

upstream_totals = upstream.sum(axis=1).values[:, None]
heatmap = (upstream.values / upstream_totals).T @ downstream.values
heatmap = pd.DataFrame(heatmap, index=heatmap_rois, columns=heatmap_rois)

# Use hierarchical clustering to find a pleasing column order.
# I seem to get better results if I take the log of the data before clustering.
from scipy.cluster.hierarchy import linkage, dendrogram
Z = linkage(np.log(heatmap + 1), method='single', optimal_ordering=True)
d = dendrogram(Z, no_plot=True)

# hvplot provides a convenient heatmap plot.
# To use the 'logz' colormap mode, I need to avoid
# zeros, hence "heatmap + 1".
import hvplot.pandas
((heatmap + 1).iloc[d['leaves'], d['leaves']]
     .hvplot.heatmap(
         x='columns',
         y='index',
         height=1000,
         width=1000,
         cmap='blues',
         xaxis='top',
         flip_yaxis=True,
         logz=True
).opts(xrotation=45))