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))
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)
In [4]:
# Sanity check
c = 0
for s in segAll:
if s.kvpairs['is_spine'] == '1':
c += 1
print c
print len(segAll)
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))
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)
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)
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)
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)
In [12]:
print graph_touch.number_of_edges()
print graph_touch.number_of_nodes()
print '\ntime elapsed so far: ' + str(time()-start)
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)
In [14]:
print graph_synapse.number_of_edges()
print graph_synapse.number_of_nodes()
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)