one 'X' means claim exists, two 'XX' means verified
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
token = 'kasthuri2015_ramon_v4'
channel = 'neurons'
segAll = nd.get_ramon(token,channel,id_segment)
print '\ntime elapsed so far: ' + str(time.time()-start)
In [4]:
z = []
for x in segAll:
z.append(x.kvpairs['segment_subtype'])
print np.unique(z)
In [5]:
'''
Count Dendrite Things
193 dendrites, 92% spiny and purportedly excitatory (177/193), remainder smooth
'''
# TODO: Maybe count difference based on cylinder of interest
dendrite_neurites = []
dendrite_neurons = []
dendrite_excitatory = []
dendrite_smooth = []
spine = 0
spine_id = []
for x in segAll:
if x.segmentclass == 2:
dendrite_neurites.append(x.id)
dendrite_neurons.append(x.neuron)
if x.kvpairs['segment_subtype'] == 'spiny':
dendrite_excitatory.append(x.neuron)
if x.kvpairs['segment_subtype'] == 'smooth':
dendrite_smooth.append(x.neuron)
if x.kvpairs['is_spine'] == '1':
spine += 1
spine_id.append(x.id)
print 'The total number of dendrite neurites in the three cylinders is: ' + str(len(dendrite_neurites))
print 'The total number of dendrites in the three cylinders is: ' + str(len(np.unique(dendrite_neurons)))
print 'The total number of excitatory dendrites in the three cylinders is: {}'.format(len(np.unique(dendrite_excitatory)))
print 'The percentage of excitatory dendrites is: ' + \
str(round(100.0*len(np.unique(dendrite_excitatory))/len(np.unique(dendrite_neurons)),0))
print 'The percentage of smooth dendrites is: ' + \
str(round(100.0*len(np.unique(dendrite_smooth))/len(np.unique(dendrite_neurons)),0))
print 'The total number of spines in the three cylinders is: ' + str(spine)
In [6]:
'''
Count Axon Things
1407 unmyelinated axons, 93% excitatory (1308/1407); most of remaining are inhibitory, 5 have no/ambiguous
'''
axon_neurites = []
axon_neurons = []
axon_excitatory = []
axon_inhibitory = []
for x in segAll:
if x.segmentclass == 1:
axon_neurites.append(x.id)
axon_neurons.append(x.neuron)
if x.kvpairs['segment_subtype'] == 'excitatory':
axon_excitatory.append(x.neuron)
if x.kvpairs['segment_subtype'] == 'inhibitory':
axon_inhibitory.append(x.neuron)
n_axons = len(np.unique(axon_neurons))
n_excitatory_axon = len(np.unique(axon_excitatory))
n_inhibitory_axon = len(np.unique(axon_inhibitory))
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(n_axons)
print 'The total number of excitatory axons in the three cylinders is: ' + str(n_excitatory_axon)
print 'The percentage of excitatory axons is: {}'.format(round(100.0*n_excitatory_axon/n_axons,0))
print 'The total number of inhibitory axons in the three cylinders is: ' + str(n_inhibitory_axon)
print 'The percentage of inhibitory axons is: {}'.format(round(100.0*n_inhibitory_axon/n_axons,0))
In [7]:
'''
Count Other Things
We also observed astrocytic processes, myelinated axons, oligodendrocyte processes and 20 entities we could not easily classify
'''
axon_myelinated_neurites = []
axon_myelinated_neurons = []
glia_astrocytes = []
glia_oligo = []
other_obj = []
c = 0
for x in segAll:
if x.kvpairs['segment_subtype'] == 'myelinated':
# all myelinated objects are axons in this dataset
axon_myelinated_neurites.append(x.id)
axon_myelinated_neurons.append(x.neuron)
if x.segmentclass == 0:
if x.kvpairs['segment_subtype'] == 'glia_oligodendrocyte':
glia_oligo.append(x.id)
if x.kvpairs['segment_subtype'] == 'glia_astrocyte':
glia_astrocytes.append(x.id)
if x.kvpairs['segment_subtype'] == 'other':
other_obj.append(x.id)
c += 1
n_myelinated_neurons = len(np.unique(axon_myelinated_neurons))
n_glia_astrocytes = len(glia_oligo)
n_glia_oligo = len(glia_astrocytes)
n_glia = n_glia_astrocytes + n_glia_oligo
n_other = len(other_obj) + c-n_myelinated_neurons-n_glia_astrocytes-n_glia_oligo
print 'Not counted or other objects in the three cylinders is: ' + str(n_other)
print 'The total number of myelinated axons in the three cylinders is: ' + str(n_myelinated_neurons)
print 'The total number of astrocytes in the three cylinders is: ' + str(n_glia_astrocytes)
print 'The total number of oligodendrocytes in the three cylinders is: ' + str(n_glia_oligo)
print 'The total number of glia in the three cylinders is: ' + str(n_glia)
In [8]:
# A priori known bounds for cylinders. Alternatively we could sweep over entire volume - this is more efficient.
# TODO: assume that all synapses are inside cylinders, which we know to be true - should do with manual masking or a
# RAMONId predicate query
token = 'kasthuri2015_ramon_v4'
channel = 'neurons'
res = 3
xbox = [694,1794];
ybox = [1750, 2460];
zbox = [1004, 1379];
# These calls take about 60 seconds to execute
rcyl = nd.get_volume('kat11redcylinder','annotation', xbox[0], xbox[1], ybox[0], ybox[1], zbox[0], zbox[1], resolution = res)
seg_masked = nd.get_volume(token, channel, xbox[0], xbox[1], ybox[0], ybox[1], zbox[0], zbox[1], resolution = res)
mask = rcyl.cutout > 0
seg = seg_masked.cutout
seg[mask == 0] = 0
In [9]:
'''
Area calculations based on all 3 cylinders
- Neuronal processes (axons and dendrites) occupy 92% of the cellular volume with glial processes occupying much of the remaining 8%
- The non-cellular (extracellular) space accounts for 6% of the total volume, less than half the space estimates from living brains
- Axons extend into a 7 fold greater volume than dendrites on average.
'''
def ismember(A, B):
return [ np.sum(a == B) for a in A ]
import skimage.measure as measure
seg = np.asarray(seg, dtype='uint32')
rp = measure.regionprops(seg,intensity_image=None)
dendrite_neurites = np.asarray(dendrite_neurites,dtype='int')
axon_neurites = np.asarray(axon_neurites,dtype='int')
glia = np.concatenate([glia_astrocytes, glia_oligo])
glia = np.asarray(glia, dtype='int')
glia_area = 0
dendrite_area = 0
axon_area = 0
for x in rp:
#print x
if np.sum(x.label == axon_neurites):
axon_area += x.area
if np.sum(x.label == dendrite_neurites):
dendrite_area += x.area
if np.sum(x.label == glia):
glia_area += x.area
mask_area = np.count_nonzero(seg)
neurite_area = axon_area + dendrite_area
cellular_volume = neurite_area + glia_area
print 'Total voxel area of axons in volume: ' + str(axon_area)
print 'Total voxel area of dendrites in volume: ' + str(dendrite_area)
print 'Total voxel area of neurites in volume: ' + str(neurite_area)
print 'Total voxel area of glia in volume: ' + str(glia_area)
print 'Total voxel area of masked region: ' + str(mask_area)
print 'Percentage of voxels that are neurites in cellular volume: ' + str(round(1.0*neurite_area/cellular_volume,2))
print 'Percentage of voxels that are glia in cellular volume: ' + str(round(1.0*glia_area/cellular_volume,2))
print 'Percentage of voxels that are extracellular space: ' + str(round(1.0-1.0*(cellular_volume)/mask_area,4))
In [10]:
'''
Axons v. Dendrites
Axons extend into a 7 fold greater volume than dendrites on average.
'''
ratio_voxel_axon_dendrite = 1.0*axon_area / dendrite_area
ratio_count_axon_dendrite = 1.0*len(np.unique(axon_neurons))/ len(np.unique(dendrite_neurons))
print 'The voxel ratio of axons to dendrites is: ' + str(ratio_voxel_axon_dendrite)
print 'The neuron ratio of axons to dendrites is: ' + str(ratio_count_axon_dendrite)
In [11]:
'''
Looking for Orphans
In 500 um^3 cylinder, no axonal or dendritic orphans, 568 spines and 601 terminal axon branches
'''
import copy
import scipy.ndimage.morphology as morpho
mask2 = morpho.binary_erosion(mask, structure=np.ones((3,3,1)))
seg2 = copy.deepcopy(seg)
all_inside_obj = np.unique(seg2[mask2])
all_obj = np.unique(seg)
# Orphans are objects that don't hit a boundary
print 'Number of orphans is: ' + str(len(all_obj)-len(all_inside_obj))
# Count spines
spine_rc = 0
spine_id = np.asarray(spine_id,dtype='int')
for x in all_obj:
if x in spine_id:
spine_rc += 1
print 'Total number of spines in volume: ' + str(spine_rc)
# Count objects that are axon neurites
axon_rc = 0
axon_neurites = np.asarray(axon_neurites,dtype='int')
for x in all_obj:
if x in axon_neurites:
axon_rc += 1
print 'Total number of axon neurites in volume: ' + str(axon_rc)
# Count terminal axons TODO: requires more thought
# Take seg2, split into connected components for each id
# If axon has two or more hits assume passing through - 1 hit = terminal
In [12]:
print 'timing for this notebook: ' + str(time.time()-start)
print 'ndio version used: ' + str(ndio.version)
In [ ]: