This notebook builds spine-axon touch and synapse graphs for the purposes of testing Peter's Rule


In [1]:
import numpy as np
import ndio.remote.neurodata as neurodata
import ndio.ramon as ramon
from time import time
import ndio

start = time()

token = 'kasthuri2015_ramon_v4'
channel = 'neurons'
res = 3

In [2]:
'''
Identify available objects: 
1600 different neurons within this small region of mammalian brain (several billionths of the whole brain)
'''
nd = neurodata(hostname='brainviz1.cs.jhu.edu')

id_segment = nd.get_ramon_ids(token, channel, ramon_type=ramon.RAMONSegment)
id_neuron = nd.get_ramon_ids(token, channel, ramon_type=ramon.RAMONNeuron)
print 'Number of RAMON segments (unique objects) is: ' + str(len(id_segment))
print 'In the 3 volume cylinder we identified {} (RAMON) neurons'.format(len(id_neuron))


Number of RAMON segments (unique objects) is: 3945
In the 3 volume cylinder we identified 1907 (RAMON) neurons

In [3]:
# Retrieve metadata for all objects

# Get all segments
segAll = nd.get_ramon(token,channel,id_segment)

print '\ntime elapsed so far: ' + str(time()-start)


time elapsed so far: 13.5934951305

In [4]:
# Sanity check

c = 0
for s in segAll:
    if s.kvpairs['is_spine'] == '1':
            c += 1
print c
print len(segAll)


1295
3945

In [5]:
# get all data in 3 cylinders
import numpy as np

# A priori known bounds for cylinders.  Alternatively we could sweep over entire volume - this is more efficient.

xbox = [173,449] #[694,1794]
ybox = [437, 615] #[1750, 2460]
zbox = [1004, 1379];
res = 5
# data are pre-masked
segments = nd.get_cutout(token, channel, xbox[0], xbox[1], ybox[0], ybox[1], zbox[0], zbox[1], resolution = res)

In [6]:
im = nd.get_cutout('kasthuri11cc', 'image', xbox[0], xbox[1], ybox[0], ybox[1], zbox[0], zbox[1], resolution = res)
import ndparse as ndp
%matplotlib inline  

ndp.plot(im, segments, slice =100, cmap1='gray', alpha=0.3)# for each object, loop through and relabel



In [7]:
'''
Identify Spines and Axons
'''

spine_id = []
spine_idx = []
axon_neurites = []
axon_neurons = []
axon_idx = []
c = 0

for x in segAll:
    #print x
    if x.segmentclass == 1:
        axon_neurites.append(int(x.id))
        axon_neurons.append(int(x.neuron))
        axon_idx.append(c)
    if x.kvpairs['is_spine'] == '1':
        spine_id.append(int(x.id))
        spine_idx.append(c)
    c += 1
print 'The total number of axon neurites in the three cylinders is: ' + str(len(axon_neurites))
print 'The total number of axons in the three cylinders is: ' + str(len(np.unique(axon_neurons)))

print 'The total number of spines in the three cylinders is: ' + str(len(spine_id))


The total number of axon neurites in the three cylinders is: 1766
The total number of axons in the three cylinders is: 1423
The total number of spines in the three cylinders is: 1295

In [8]:
# for each spinedil, add edge for each axon touch.neuron, [3, 3, 1]
# for each spine, add edge for each axon edge, which is spine.synapse.otherpartner.neuron
# save both graphs
# for each node in both graphs, add node attribute = spine or axon

import copy
segments2 = copy.copy(segments)
segments2 = np.ravel(segments)

for x in range(0,len(spine_id)):
    segments2[segments2 == id_segment[spine_idx[x]]]
print '\ntime elapsed so far: ' + str(time()-start)


time elapsed so far: 126.299653053

In [9]:
# TODO BETTER

# make spine mask
# make axon mask

spine_vol = np.zeros_like(segments)
for i in spine_id:
    #print i 
    spine_vol[segments == i] = i

axon_vol = np.zeros_like(segments)
print '\ntime elapsed so far: ' + str(time()-start)

for i in axon_neurites:
    #print i 
    axon_vol[segments == i] = i

print '\ntime elapsed so far: ' + str(time()-start)


time elapsed so far: 292.498157024

time elapsed so far: 581.035701036

In [10]:
spine_dil = copy.copy(spine_vol)

import mahotas
# dilate all spines
spine_dil = mahotas.dilate(spine_dil, Bc=np.ones([11,11,3]))

# remove all spines in original (to avoid overlap weirdness)
spine_dil[spine_vol > 0] = 0
print '\ntime elapsed so far: ' + str(time()-start)


time elapsed so far: 582.278528929

In [11]:
# touch graph
# scale 5, epsilon, multiedge

import networkx as nx
graph_touch = nx.Graph()
axon_neurons = np.asarray(axon_neurons)

for i in spine_id:
    graph_touch.add_node(i)
    graph_touch.node[i]['type'] = 'spine'

for i in axon_neurons:
    graph_touch.add_node(i)
    graph_touch.node[i]['type'] = 'axon'

for s in spine_id:
    #print str(s).zfill(4),

    val = axon_vol[s == spine_dil]
    val = np.unique(val)
    val = val[val > 0]
    
    for v in val:
        idx = np.where(id_segment == v)[0][0]        
        graph_touch.add_edge(s,segAll[idx].neuron) #avoid retrieving 1 by 1 if possible
print '\ntime elapsed so far: ' + str(time()-start)


time elapsed so far: 618.207873106

In [12]:
print graph_touch.number_of_edges()
print graph_touch.number_of_nodes()
print '\ntime elapsed so far: ' + str(time()-start)


15680
2718

time elapsed so far: 618.235081911

In [13]:
# Synapse graph

# Do with local data cache:
synid = nd.get_ramon_ids(token, 'synapses')
synAll = nd.get_ramon(token, 'synapses',synid)

import networkx as nx
graph_synapse = nx.Graph()

for i in spine_id:
    graph_synapse.add_node(i)
    graph_synapse.node[i]['type'] = 'spine'

for i in axon_neurons:
    graph_synapse.add_node(i)
    graph_synapse.node[i]['type'] = 'axon'

for s in spine_id:
    #print str(s).zfill(4),
    idx = np.where(np.asarray(id_segment) == s)[0][0] # replace with local data cache - WAY faster  #nd.get_ramon(token,'neurons',s)
    r = segAll[idx]
    if len(r.synapses) > 0:
        syn = r.synapses[0]
        #ss = nd.get_ramon(token, 'synapses', syn)
        idx = np.where(np.asarray(synid) == syn)[0][0]
        ss = synAll[idx]
        axon = ss.kvpairs['presynaptic']
    
        if axon is not None:
            axon = int(axon)
            graph_synapse.add_edge(s,axon,weight=1)
print '\ntime elapsed so far: ' + str(time()-start)


time elapsed so far: 625.517513037

In [14]:
print graph_synapse.number_of_edges()
print graph_synapse.number_of_nodes()


1178
3359

In [15]:
nx.write_edgelist(graph_synapse, 'kasthuri2015_ramon_v4_graph_synapse.edgelist')
nx.write_graphml(graph_synapse, 'kasthuri2015_ramon_v4_graph_synapse.graphml')
nx.write_edgelist(graph_touch, 'kasthuri2015_ramon_v4_graph_touch.edgelist')
nx.write_graphml(graph_touch, 'kasthuri2015_ramon_v4_graph_touch.graphml')

In [16]:
import datetime
i = datetime.datetime.now()
print '\ntime elapsed so far: ' + str(time()-start)

print ("Current date & time = %s" % i)


time elapsed so far: 626.392879009
Current date & time = 2016-06-29 10:28:18.646500