#!/usr/bin/env python # coding: utf-8 # In order to build the masks for KlustaKwik, SpikeDetekt2 needs to know the graph of adjacencies between all of the recording sites on a recording array. However, a bad channel for a recording session can drastically change the requirements for a probe geometry. # # Luckily, there's some math out there that can help us to build these adjacency matrices based solely on the geometry of the recording sites. I've [already demonstrated how to use Delauney triangulation to do this for electrode arrays that are very large](http://www.justinkiggins.com/blog/defining-large-array-geometries-for-klustakwik/), where it would be a burden to define the geometry by hand. That was fun, but the technique has an application that I've already implemented into my own pipeline to extract population spiking activity. # # ## We start with site maps in terms of Neuronexus recording sites # # We are using we NeuroNexus probes. The recording sites are 1-indexed, whereas KlustaKwik uses 0-indexed channels (where each channel is a column in a 2D array). Further, due to our hardware configuration, NeuroNexus "Site 1" doesn't necessarily correspond to KlustaKwik's "Channel 0". So we are going to define a site map that gives us the channel/index for each neuronexus site. The dictionary `s` is composed of site:channel pairs, where 'site' is the Neuronexus recording site and 'channel' is the index of that site's recording in the Klustakwik `*.raw.kwd` file. # In[1]: # site:channel s = {11: 0, 12: 1, 22: 2, 10: 3, 21: 4, 23: 5, 9: 6, 13: 7, 24: 8, 8: 9, 20: 10, 25: 11, 7: 12, 14: 13, 26: 14, 6: 15, 19: 16, 27: 17, 5: 18, 15: 19, 28: 20, 4: 21, 18: 22, 29: 23, 1: 24, 16: 25, 32: 26, 2: 27, 17: 28, 31: 29, 3: 30, 30: 31, } # ## Sometimes you lose channels # # Sometimes you don't have all of your channels. Maybe a channel is bad on your headstage or you're using one channel as a reference. Things happen. That shouldn't prevent us from using the most complete adjacency matrix that we can. # # We're going pretend like we've have some dead (or unrecorded) sites by setting their channel to `None` # In[2]: from numpy import random from pprint import pprint for dead_site in random.choice(s.values(),5): s[dead_site]=None pprint(s) # ## Define our probe geometry in terms of the site map # # If I were to manually define the adjacencies for each site on the NeuroNexus A1x32-Poly3-6mm-50 probe, this is what I would sketch: # # # (3) (30) # | \ / | # (2)-(17)-(31) # | / \ | # (1)-(16)-(32) # | \ / | # (4)-(18)-(29) # | / \ | # (5)-(15)-(28) # | \ / | # (6)-(19)-(27) # | / \ | # (7)-(14)-(26) # | \ / | # (8)-(20)-(25) # | / \ | # (9)-(13)-(24) # | \ / | # (10)-(21)-(23) # | / \ | # (11)-(12)-(22) # # Each number is the NeuroNexus site number and I've drawn my edges between them. This particular probe is a grid, so I have each square drawn and (for the sake of completeness) each square is split into two triangles. For KlustaKwik, I would need to write these out as a list of tuples. # # Instead, we're just going to define the probe *geometry*, giving the x,y coordinates of each site, where `s[1]` indicates site 1. If a different electrophysiology rig has a different hardware configuration, this will change, so it's more robust to define this geometry in terms of the NeuroNexus sites. # # Even though KlustaKwik doesn't care about the geometry dict (except to plot the channels in Klustaviewa), we do. We are going to generate the graph on the fly from this geometry information, so I've defined it in microns, where the lower left recording site (Site 11) is (0,0). # In[3]: channel_groups = { # Shank index. 0: { # List of channels to keep for spike detection. 'channels': s.values(), # 2D positions of the channels # channel: (x,y) 'geometry': { s[3]: (0,500), # column 0 s[2]: (0,450), s[1]: (0,400), s[4]: (0,350), s[5]: (0,300), s[6]: (0,250), s[7]: (0,200), s[8]: (0,150), s[9]: (0,100), s[10]: (0,50), s[11]: (0,0), s[17]: (50,450), # column 1 s[16]: (50,400), s[18]: (50,350), s[15]: (50,300), s[19]: (50,250), s[14]: (50,200), s[20]: (50,150), s[13]: (50,100), s[21]: (50,50), s[12]: (50,0), s[30]: (100,500), # column 2 s[31]: (100,450), s[32]: (100,400), s[29]: (100,350), s[28]: (100,300), s[27]: (100,250), s[26]: (100,200), s[25]: (100,150), s[23]: (100,100), s[24]: (100,50), s[22]: (100,0), } } } pprint(channel_groups) # ## Clean up the dead channels # # We have a lot of `None`s in there, indicating dead channels. Let's strip those out. # In[4]: new_group = {} for gr, group in channel_groups.iteritems(): new_group[gr] = { 'channels': [], 'geometry': {} } new_group[gr]['channels'] = [ch for ch in group['channels'] if ch is not None] new_group[gr]['geometry'] = {ch:xy for (ch,xy) in group['geometry'].iteritems() if ch is not None} channel_groups = new_group pprint(channel_groups) # ## Build the adjacency graph from the probe geometry # # Next, we'll use SciPy's built-in Delaunay Triangulation to construct a graph from the probe geometry # In[7]: from scipy import spatial from scipy.spatial.qhull import QhullError def get_graph_from_geometry(geometry): # let's transform the geometry into lists of channel names and coordinates chans,coords = zip(*[(ch,xy) for ch,xy in geometry.iteritems()]) # we'll perform the triangulation and extract the try: tri = spatial.Delaunay(coords) except QhullError: # oh no! we probably have a linear geometry. chans,coords = list(chans),list(coords) x,y = zip(*coords) # let's add a dummy channel and try again coords.append((max(x)+1,max(y)+1)) tri = spatial.Delaunay(coords) # then build the list of edges from the triangulation indices, indptr = tri.vertex_neighbor_vertices edges = [] for k in range(indices.shape[0]-1): for j in indptr[indices[k]:indices[k+1]]: try: edges.append((chans[k],chans[j])) except IndexError: # let's ignore anything connected to the dummy channel pass return edges def build_geometries(channel_groups): for gr, group in channel_groups.iteritems(): group['graph'] = get_graph_from_geometry(group['geometry']) return channel_groups channel_groups = build_geometries(channel_groups) print(channel_groups[0]['graph']) # ## That was way easier than doing it by hand # Let's see what it looks like... # In[6]: get_ipython().run_line_magic('pylab', 'inline') def plot_channel_groups(channel_groups): n_shanks = len(channel_groups) f,ax = plt.subplots(1,n_shanks,squeeze=False) for sh in range(n_shanks): coords = [xy for ch,xy in channel_groups[sh]['geometry'].iteritems()] x,y = zip(*coords) ax[sh,0].scatter(x,y,color='0.2') for pr in channel_groups[sh]['graph']: points = [channel_groups[sh]['geometry'][p] for p in pr] ax[sh,0].plot(*zip(*points),color='k',alpha=0.2) ax[sh,0].set_xlim(min(x)-10,max(x)+10) ax[sh,0].set_ylim(min(y)-10,max(y)+10) ax[sh,0].set_xticks([]) ax[sh,0].set_yticks([]) ax[sh,0].set_title('group %i'%sh) axis('equal') plot_channel_groups(channel_groups) # # Fantastic # # If we had just dropped our missing sites, we'd have huge gaping holes in our graph. Instead, we now have a bunch of edges linking past the dead sites, keeping our graph in tact.