Skip to content

Commit

Permalink
write total_traces to file
Browse files Browse the repository at this point in the history 
calculate total_traces for dendrite submembrane and spine head
plot region specific total_traces
  • Loading branch information
Avrama Blackwell committed Jul 13, 2021
1 parent c23176a commit f6b40d4
Show file tree
Hide file tree
Showing 3 changed files with 64 additions and 68 deletions.
33 changes: 16 additions & 17 deletions nrd_group.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,22 +9,20 @@

ms_to_sec=1000
class nrdh5_group(object):
def __init__(self,args,tot_species=[]):
self.ftuples,self.parlist,self.params=h5utils.argparse(args)
def __init__(self,fileroot,parameters,tot_species=[]):
self.ftuples,self.parlist,self.params=h5utils.create_filenames(fileroot,parameters)
self.file_set_conc={k[1]:{} for k in self.ftuples}
self.time_set={k[1]:{} for k in self.ftuples}
self.spatial_means={k[1]:{} for k in self.ftuples}
self.regions_means={k[1]:{} for k in self.ftuples}
self.regions_structure_means={k[1]:{} for k in self.ftuples}
self.spine_means={k[1]:{} for k in self.ftuples}
if len(tot_species):
#Expand this to include spine_means and other spatial information
self.file_set_tot={k[1]:{sp:[] for sp in tot_species} for k in self.ftuples}
self.tot_species=tot_species
self.endtime={k[1]:{sp:[] for sp in tot_species} for k in self.ftuples}
self.endtime={k[1]:{sp:[] for sp in self.tot_species} for k in self.ftuples}
else:
self.tot_species=[]
self.file_set_tot={}
self.file_set_tot={'Overall':{}}

def conc_arrays(self,data):
self.molecules=data.molecules
Expand All @@ -45,14 +43,15 @@ def conc_arrays(self,data):
self.spine_means[data.parval][molecule]=data.means['spines'][molecule]
else:
self.spatial_data=None
if len(data.ss_tot):
#Expand this to include spine_means and other spatial information
for imol,sp in enumerate(self.file_set_tot[data.parval].keys()):
self.file_set_tot[data.parval][sp]=data.ss_tot[sp][:,:]
self.sstart[sp]=data.sstart[sp]
self.ssend[sp]=data.ssend[sp]
self.dt[sp]=data.dt[sp]
self.endtime[data.parval][sp]=data.endtime[sp]
if len(self.tot_species):
for region in data.total_trace.keys():
self.file_set_tot[region]={k[1]:{sp:[] for sp in self.tot_species} for k in self.ftuples}
for imol,sp in enumerate(self.file_set_tot[region][data.parval].keys()):
self.file_set_tot[region][data.parval][sp]=data.total_trace[region][sp][:,:]
self.sstart[sp]=data.sstart[sp]
self.ssend[sp]=data.ssend[sp]
self.dt[sp]=data.dt[sp]
self.endtime[data.parval][sp]=data.endtime[sp]

def trace_features(self,trials,window_size,lo_thresh_factor=0.2,hi_thresh_factor=0.8,std_factor=1,numstim=1,end_baseline_start=0,filt_length=5,aucend=None):
import operator
Expand All @@ -78,7 +77,7 @@ def exceeds_thresh_points(traces,startpoints,thresh,relate,endpoints=[-1 for t i
if len(pointset):
earliest_points[i]=np.min(pointset)
return earliest_points
for ii,(molecules,traces) in enumerate(zip([self.molecules,self.tot_species],[self.file_set_conc,self.file_set_tot])):
for ii,(molecules,traces) in enumerate(zip([self.molecules,self.tot_species],[self.file_set_conc,self.file_set_tot['Overall']])):
for parnum,(fname,par) in enumerate(self.ftuples):
for jmol,mol in enumerate(molecules):
imol=jmol+ii*len(self.molecules)
Expand Down Expand Up @@ -171,14 +170,14 @@ def exceeds_thresh_points(traces,startpoints,thresh,relate,endpoints=[-1 for t i

def write_features(self,feature_list,arg0,write_trials=False):
outfname=arg0+'-'+'analysis'+'-'.join([i for i in self.params])+'.txt' #test whether this works for single file
#print('in write_features',outfname, feature_list)
print('in write_features',outfname, feature_list)
if len(self.ftuples)==1:
outputdata=arg0
header='file '
else:
outputdata=['-'.join([str(p) for p in par[1]]) for par in self.ftuples]
header='-'.join([i for i in self.params])+' '
header+=' '.join([m+'_'+f+'_mean ' for m in self.molecules+self.tot_species for f in feature_list]+ [m+'_'+f+'_std ' for m in self.molecules+self.tot_species for f in feature_list])
header+=' '.join([m+'_'+f+'_mean ' for m in list(self.molecules)+self.tot_species for f in feature_list]+ [m+'_'+f+'_std ' for m in list(self.molecules)+self.tot_species for f in feature_list])
for feat in feature_list:
outputdata=np.column_stack((outputdata,np.round(self.mean_feature[feat].T/ms_to_sec,3),np.round(self.std_feature[feat].T/ms_to_sec,3)))
f=open(outfname, 'w')
Expand Down
60 changes: 29 additions & 31 deletions nrd_output.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ def spatial_structures(self,bins,dendname):
if any(v==0 for v in vol):
print ("**********Too many bins for spatial average****************")
#
def rows_and_columns(self,molecules,args):
def rows_and_columns(self,molecules,start_end):
#which columns in the data set contain counts of molecules of interest
self.all_molecules=h5utils.decode(self.data['model']['species'][:])
if molecules is not None:
Expand All @@ -65,12 +65,8 @@ def rows_and_columns(self,molecules,args):
self.out_location=out_location
self.dt=dt
self.rows=rows
if len(args)>3:
arg=args[3]
else:
arg=None
#Which "rows" should be used for baseline value, specifed in args[3]. If different for each file then will have problems later
sstart,ssend=h5utils.sstart_end(self.molecules,self.out_location,self.dt,self.rows,arg)
#Which "rows" should be used for baseline value. If different for each file then will have problems later
sstart,ssend=h5utils.sstart_end(self.molecules,self.out_location,self.dt,self.rows,start_end)
self.sstart=sstart
self.ssend=ssend

Expand Down Expand Up @@ -107,10 +103,16 @@ def average_over_voxels(self):

def write_average(self):
import os ########### NEED TO PUT SPATIAL STUFF IN SEPARATE FILE?
for mol in self.molecules:
outfilename=os.path.splitext(os.path.basename(self.fname))[0]+mol+'_avg.txt'
output_means=np.mean(self.OverallMean[mol],axis=0) #average over trials
output_std=np.std(self.OverallMean[mol],axis=0)
for mol in list(self.molecules)+list(self.tot_species.keys()):
outfilename=os.path.splitext(os.path.basename(self.fname))[0]+mol+'_avg.txt'
if mol in self.molecules:
output_means=np.mean(self.OverallMean[mol],axis=0) #average over trials
output_std=np.std(self.OverallMean[mol],axis=0)
time=self.time[mol]
elif mol in self.tot_species.keys():
output_means=np.mean(self.total_trace['Overall'][mol],axis=0)
output_std=np.std(self.total_trace['Overall'][mol],axis=0)
time=np.linspace(0,self.endtime[mol], len(output_means))
mean_header=mol+'_'.join([str(q) for q in self.parval])+'_All '
if self.maxvols>1:
for mean_dict in self.means.values():
Expand All @@ -123,18 +125,18 @@ def write_average(self):
#print(outfilename,output_header)
f=open(outfilename, 'w')
f.write(output_header)
np.savetxt(f, np.column_stack((self.time[mol],output_means,output_std)), fmt='%.4f', delimiter=' ')
np.savetxt(f, np.column_stack((time,output_means,output_std)), fmt='%.4f', delimiter=' ')
f.close()
#
#FOR SIGNATURE:
#add option to baseline subtract - separate function?
def total_subspecies(self,tot_species,sub_species,args,outset='__main__',weights={}):
def total_subspecies(self,tot_species,sub_species,start_end,outset='__main__',weights={}):
self.tot_species={t:[] for t in tot_species}
self.ss_tot={}
self.total_trace={'Overall':{}}
if self.dsm_vox:
self.dsm_tot={}
self.total_trace['dsm']={}
if self.spinelist:
self.head_tot={}
self.total_trace['spine']={}
for imol,mol in enumerate(tot_species):
mol_set=[];tot_wt=0
#### first set up arrays of all species (sub_species) that are a form of the molecule
Expand All @@ -144,21 +146,17 @@ def total_subspecies(self,tot_species,sub_species,args,outset='__main__',weights
mol_set=[subsp for subsp in self.all_molecules if mol in subsp] #h5utils.decode(self.data['model']['output'][outset]['species'][:])
self.tot_species[mol]=mol_set
out_location,dt,rows=h5utils.get_mol_info(self.data, mol_set)
if len(args)>3:
arg=args[3]
else:
arg=None
sstart,ssend=h5utils.sstart_end(mol_set,out_location,dt,rows,arg)
sstart,ssend=h5utils.sstart_end(mol_set,out_location,dt,rows,start_end)
if len(np.unique(list(dt.values())))>1:
print('PROBLEM, cannot total subspecies for',mol,' with different dt',dt)
else:
dt=np.mean(list(dt.values()))
self.endtime={mol:(rows[0]-1)*dt for mol in tot_species}
self.ss_tot[mol]=np.zeros((len(self.trials),rows[0]))
self.total_trace['Overall'][mol]=np.zeros((len(self.trials),rows[0]))
if self.dsm_vox:
self.dsm_tot[mol]=np.zeros((len(self.trials),rows[0]))
self.total_trace['dsm'][mol]=np.zeros((len(self.trials),rows[0]))
if self.spinelist:
self.head_tot[mol]=np.zeros((len(self.trials),rows[0]))
self.total_trace['spine'][mol]=np.zeros((len(self.trials),rows[0]))
#### second, find molecule index of the sub_species and total them
for subspecie in mol_set:
self.dt[mol]=dt
Expand All @@ -172,18 +170,18 @@ def total_subspecies(self,tot_species,sub_species,args,outset='__main__',weights
mol_pop=self.data[trial]['output'][outname]['population'][:,:,outset['mol_index']] # FIXME won't work if > 1 outset
else:
mol_pop=self.counts[subspecie][trialnum]
#print('tot_species',mol,subspecie,np.shape(mol_pop),np.shape(self.ss_tot[mol][trialnum]))
self.ss_tot[mol][trialnum]+=wt*np.sum(mol_pop,axis=1)/self.TotVol/mol_per_nM_u3
#print('tot_species',mol,subspecie,np.shape(mol_pop),np.shape(self.total_trace['Overall'][mol][trialnum]))
self.total_trace['Overall'][mol][trialnum]+=wt*np.sum(mol_pop,axis=1)/self.TotVol/mol_per_nM_u3
#then total sub_species in submembrane and spine head, if they exist
if self.dsm_vox:
self.dsm_tot[mol][trialnum]+=wt*np.sum(mol_pop[:,self.dsm_vox['vox']],axis=1)/self.dsm_vox['vol']*self.dsm_vox['depth']/mol_per_nM_u3
self.total_trace['dsm'][mol][trialnum]+=wt*np.sum(mol_pop[:,self.dsm_vox['vox']],axis=1)/self.dsm_vox['vol']*self.dsm_vox['depth']/mol_per_nM_u3
if self.spinelist:
self.head_tot[mol][trialnum]+=wt*np.sum(mol_pop[:,self.region_dict[self.head]['vox']],axis=1)/self.region_dict[self.head]['vol']/mol_per_nM_u3
outputline=str(self.parval)+' TOTAL: '+str(np.round(self.ss_tot[mol][0,0],3))+' nM'
self.total_trace['spine'][mol][trialnum]+=wt*np.sum(mol_pop[:,self.region_dict[self.head]['vox']],axis=1)/self.region_dict[self.head]['vol']/mol_per_nM_u3
outputline=str(self.parval)+' TOTAL: '+str(np.round(self.total_trace['Overall'][mol][0,0],3))+' nM'
if self.dsm_vox:
outputline +=', sp: '+str(np.round(self.head_tot[mol][0,0],3))+' nM'
outputline +=', sp: '+str(np.round(self.total_trace['spine'][mol][0,0],3))+' nM'
if self.spinelist:
outputline +=', dsm: '+str(np.round(self.dsm_tot[mol][0,0],3))+' pSD'
outputline +=', dsm: '+str(np.round(self.total_trace['dsm'][mol][0,0],3))+' pSD'
print(outputline)

def print_head_stats(self):
Expand Down
39 changes: 19 additions & 20 deletions plot_h5V2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,29 +183,28 @@ def plotss(plot_mol,xparval,ss):

def plot_signature(tot_species,dataset,figtitle,colinc,textsize,thresholds=[]):
numcols=len(tot_species)
#Will need to specify whether plotting spine (and non-spine) totals, or perhaps region totals
#add additional total arrays in nrd_group
numrows= 1 #one row for each region For now, there is only an overall sum
row=0 #once multiple rows, will loop over rows, and region as y label, mol as column heading?
#Will need to specify whether plotting spine (and non-spine) totals, currently, regions are overall, dsm and spine head. need non-spine (dendrite)
numrows= len(dataset.file_set_tot.keys())
fig,axes=pyplot.subplots(numrows,numcols,sharex=True)
fig.canvas.set_window_title(figtitle+'Totals')
axis=fig.axes
for i,(param,ss_tot) in enumerate(dataset.file_set_tot.items()):
mycolor,plotlabel,par_index,map_index=get_color_label(dataset.parlist,param,colinc)
for col,(mol,trace) in enumerate(ss_tot.items()):
#print('$$$$$$$$$$ pu.ps',param,mol,np.shape(trace),mycolor,plotlabel)
ax=col+row*numrows
newtime = np.linspace(0,dataset.endtime[param][mol], np.shape(trace)[1]) #convert from ms to sec
axis[ax].plot(newtime,np.mean(trace,axis=0),label=plotlabel,color=mycolor)
axis[ax].set_title(mol+' TOTAL',fontsize=textsize)
axis[ax].set_xlabel('Time (sec)',fontsize=textsize)
axis[ax].tick_params(labelsize=textsize)
axis[row*numrows].set_ylabel('Conc (nM)',fontsize=textsize)
axis[0].legend(fontsize=legtextsize, loc='best')
if len(thresholds): #this needs to be fixed. Determine how to match thresholds to mol/sig
r=(1,0)[row==0]
axis[ax].plot([0,newtime[-1]],[thresholds[r],thresholds[r]],color='k',linestyle= 'dashed')
axis[ax+1].plot([0,newtime[-1]],[thresholds[r+numcols],thresholds[r+numcols]],color='k',linestyle= 'dashed')
for row,region in enumerate(dataset.file_set_tot.keys()):
for i,(param,total_trace) in enumerate(dataset.file_set_tot[region].items()):
mycolor,plotlabel,par_index,map_index=get_color_label(dataset.parlist,param,colinc)
for col,(mol,trace) in enumerate(total_trace.items()):
#print('$$$$$$$$$$ pu.ps',param,mol,np.shape(trace),mycolor,plotlabel)
ax=col+row*numrows
newtime = np.linspace(0,dataset.endtime[param][mol], np.shape(trace)[1]) #convert from ms to sec
axis[ax].plot(newtime,np.mean(trace,axis=0),label=plotlabel,color=mycolor)
axis[ax].set_title(mol+' TOTAL',fontsize=textsize)
axis[ax].set_xlabel('Time (sec)',fontsize=textsize)
axis[ax].tick_params(labelsize=textsize)
axis[row*numrows].set_ylabel('Conc (nM)',fontsize=textsize)
axis[0].legend(fontsize=legtextsize, loc='best')#for now put legend into panel 0
if len(thresholds): #this needs to be fixed. Determine how to match thresholds to mol/sig
r=(1,0)[row==0]
axis[ax].plot([0,newtime[-1]],[thresholds[r],thresholds[r]],color='k',linestyle= 'dashed')
axis[ax+1].plot([0,newtime[-1]],[thresholds[r+numcols],thresholds[r+numcols]],color='k',linestyle= 'dashed')
fig.canvas.draw()

def tweak_fig(fig,yrange,legendloc,legendaxis,legtextsize):
Expand Down

0 comments on commit f6b40d4

Please sign in to comment.