In [1]:
import numpy as np
import ndio.remote.neurodata as neurodata
import ndio.ramon as ramon
import time
import ndio
start = time.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()
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.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 [46]:
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)
t = time.time()
for x in range(0,len(spine_id)):
segments2[segments2 == id_segment[spine_idx[x]]]
print str(time.time()-t)
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)
for i in axon_neurites:
#print i
axon_vol[segments == i] = i
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
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[i]['type'] = 'spine'
for i in axon_neurons:
graph_touch.add_node(i)
graph_touch[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:
# TODO smarter way
r = nd.get_ramon(token,'neurons',v)
graph_touch.add_edge(s,r.neuron)
#print v
#print np.where(r.neuron == axon_neurons)
In [12]:
print graph_touch.number_of_edges()
print graph_touch.number_of_nodes()
In [ ]:
In [13]:
print '\ntime elapsed so far: ' + str(time.time()-start)
In [14]:
'''# Relabel everything
import copy
neurons = copy.copy(segments)
def ismember(a, b):
bind = {}
for i, elt in enumerate(b):
if elt not in bind:
bind[elt] = i
return [bind.get(itm, None) for itm in a] # None can be replaced by any other "not in b" value
for s in range(0,len(id_segment)):
print s
neurons[ismember(id_segment[s],segments)] = segAll[s].neuron
'''
# synapse graph
#ndp.plot(spine_dil,slice=100)
#notaxon = np.setdiff1d(np.asarray(id_segment),np.asarray(axon_neurites))
#notspine = np.setdiff1d(np.asarray(id_segment),np.asarray(spine_id))
#axon_vol = mahotas.labeled.remove_regions(segments, notaxon)
#s2 = mahotas.dilate(segments,np.ones([3,3,1]))
Out[14]:
In [52]:
# Synapse graph
import networkx as nx
graph_synapse = nx.Graph()
for i in spine_id:
graph_synapse.add_node(i)
graph_synapse[i]['type'] = 'spine'
for i in axon_neurons:
graph_synapse.add_node(i)
graph_synapse[i]['type'] = 'axon'
for s in spine_id:
print str(s).zfill(4),
r = nd.get_ramon(token,'neurons',s)
if len(r.synapses) > 0:
syn = r.synapses[0]
ss = nd.get_ramon(token, 'synapses', syn)
axon = ss.kvpairs['presynaptic']
if axon is not None:
axon = int(axon)
graph_synapse.add_edge(s,axon,weight=1)
In [47]:
print graph_synapse.number_of_edges()
print graph_synapse.number_of_nodes()
In [53]:
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')
In [ ]:
print '\ntime elapsed so far: ' + str(time.time()-start)
In [42]:
len(r.synapses)
Out[42]:
In [21]:
syn = r.synapses[0]
In [22]:
syn
Out[22]:
In [ ]: