#Python version, i.e. alternative of NRDpostAB #in python, type ARGS="par1 par2,mol1 mol2,subdir/fileroot,sstart ssend" then execfile('neurord_analysis.py') #DO NOT PUT ANY SPACES NEXT TO THE COMMAS, DO NOT USE TABS #e.g. ARGS="Ca GaqGTP,Ca GaqGTP Ip3,../Repo/plc/Model_PLCassay,15 20" #from outside python, type python neurord_analysis [par1 par2] [mol1 mol2] #Assumes that molecule outputs are integers, and the hypens used ONLY for parameters #Can process multiple parameter variations, but all files must use same morphology, and meshfile. #It will provide region averages (each spine, dendrite submembrane, cytosol) and if spatialaverage=1, #will calculate an average of n segments along the dendrite, #or whatever structure name is specified in dend variable import os import numpy as np from matplotlib import pyplot from string import * import sys import glob import header_parse as hparse import plot_utils as pu ####################################################### #indicate the name of the injection spines if you want to exclude them fake = 'FAKE' #indicate the name of submembrane region for totaling molecules that are exclusively submembrane #only relevant for tot_species calculation. this should be name of structure concatenated with sub submembname='sub' #Spatial average (=1 to process) only includes the structure "dend", and subdivides into bins: spatialaverage=0 dend="dend" bins=10 #how much info to print prnvox=1 prnheader=0 prninfo=0 showss=0 #outputavg determines whether output files are written outputavg=0 ##change these endings depending on whether using neurord3.x: meshend="*mesh.txt.out" concend='conc.txt.out' ## or neurord2.x (uncomment these) #meshend="*mesh.txt" #concend='conc.txt' #Example of how to total some molecule forms; turn off with tot_species={} #tot_species={ # "PKAtot":["PKA", "PKAcAMP2", "PKAcAMP4", "PKAr"], # "D1Rtot":["D1R","DaD1R", "GsD1R","DaD1RGs", "pDaD1RGs", "PKAcDaD1RGs"], # "pde10tot":["PDE10","pPDE10", "PDE10cAMP","pPDE10cAMP","PKAcPDE10", "PKAcPDE10cAMP"], # "Gitot":["Giabg","AChm4RGi","Gim4R", "GaiGTP", "GaiGDP", "ACGai", "ACGasGai", "ACGasGaiATP"], # "m4Rtot":["AChm4RGi","Gim4R", "m4R", "AChm4R"]} tot_species={} Avogadro=6.023e14 #to convert to nanoMoles mol_per_nM_u3=Avogadro*1e-15 try: args = ARGS.split(",") print "ARGS =", ARGS, "commandline=", args do_exit = False except NameError: #NameError refers to an undefined variable (in this case ARGS) args = sys.argv[1:] print "commandline =", args do_exit = True pattern=args[2]+'*' if len(args[0]): params=args[0].split(" ") for par in params: pattern=pattern+'-'+par+'*' else: params=[] whole_pattern=pattern+concend #A single mesh file means that all files in your list must use the same morphology meshname=pattern.split('-')[0]+meshend lastslash=rfind(pattern,'/') subdir=pattern[0:lastslash] ################################################### def sortorder(ftuple): ans = ftuple[1] #print 'sort', ftuple, '->', ans return ans fnames = glob.glob(whole_pattern) print "NUM FILES:", len(fnames), "CURRENT DIRECTORY:", os.getcwd(), ", Target directory:", subdir if len(fnames)==0: print "MESHFILES:", os.listdir(subdir+'/'+meshend) ss_tot=np.zeros((len(fnames),len(tot_species.keys()))) parlist=[] if len(args[0]): ftuples,parlist=pu.file_tuple(fnames,params) ftuples = sorted(ftuples, key=lambda x:x[1]) else: ftuples=[(fnames[0],1)] #First, read mesh file to determine how many voxels if len(fnames)>0: meshfile=glob.glob(meshname)[0] else: print "********** no meshfile **************" maxvols,vox_volume,xloc,yloc,TotVol,deltaY=hparse.read_mesh(meshfile) #prepare to plot stuff (instead of calculating averages) #plot_molecules determines what is plotted plot_molecules=args[1].split(' ') fig,axes,col_inc,scale,minpar=pu.plot_setup(plot_molecules,parlist,params) fig.suptitle(pattern.split('/')[-1]) ss=np.zeros((len(fnames),len(plot_molecules))) slope=np.zeros((len(fnames),len(plot_molecules))) peaktime=np.zeros((len(fnames),len(plot_molecules))) baseline=np.zeros((len(fnames),len(plot_molecules))) peakval=np.zeros((len(fnames),len(plot_molecules))) lowval=np.zeros((len(fnames),len(plot_molecules))) parval=[] for fnum,ftuple in enumerate(ftuples): fname=ftuple[0] parval.append(ftuple[1]) if fnum == 0: f = open(fname, 'r+') #parse the header to determine identity/structure of voxels and molecules data=f.readline() if (prnheader==1): print "header",data else: print "header not printed" #UPDATE maxvols, or number of voxels in this function regionID,structType,molecules,volnums,maxvols=hparse.header_parse(data,maxvols,prninfo) print "in neurord_analysis: vox#", volnums print " regions",regionID print " structures",structType print " mols",molecules f.close() #all voxels should be read in now with labels #extract number of unique regions (e.g. dendrite, or sa1[0]), #and create list of subvolumes which contribute to that region if maxvols>1: region_list,region_vox,region_col,region_struct_list,region_struct_vox,region_struct_col=hparse.subvol_list(structType,regionID,volnums,fake) RegVol=hparse.region_volume(region_list,region_vox,vox_volume,prnvox) RegStructVol=hparse.region_volume(region_struct_list,region_struct_vox,vox_volume,prnvox) submembVol=0 for region in region_list: smname=region+submembname if smname in region_struct_list: submembVol+=RegStructVol[region_struct_list.index(smname)] # if spatialaverage: hparse.spatial_average(xloc,yloc,bins,regionID,structType,volnums) # #Lastly, read in the data and output separate files of region averages #Can do all molecules in a list without a batch file alldata=np.loadtxt(fname,skiprows=1) time=alldata[:,0]/1000 dt=time[1] data=alldata[:,1:alldata.shape[1]] del alldata #the above eliminates the time column from the data, so that e.g., column 0 = voxel 0 # #reshape the data to create a separate dimension for each molecule rows=data.shape[0] arrays=len(molecules) if maxvols*arrays == data.shape[1]: molecule_array=np.reshape(data, (rows,arrays,maxvols)) del data else: print "UH OH! voxels:", maxvols, "molecules:", len(molecules), "columns:", data.shape[1] plot_array=np.zeros((rows,len(plot_molecules))) sstart=int(float(args[3].split(" ")[0])/dt) ssend=int(float(args[3].split(" ")[1])/dt) ##now, calculate various averages such as soma and dend, subm vs cyt, #use the above lists and volume of each region, and each region-structure # if maxvols>1: data=np.zeros((rows,maxvols)) for imol in range(arrays): if molecules[imol] in plot_molecules: data=molecule_array[:,imol,:] RegionMeans=np.zeros((len(time),len(region_list))) header='#time' #Header for output file for itime in range(len(time)): for j in range(len(region_list)): for k in region_col[j]: RegionMeans[itime,j]+=data[itime,k] #sum the molecules of the voxels in the structure, divide by Avogadro and volume # for j in range(len(region_list)): RegionMeans[:,j]/=(RegVol[j]*mol_per_nM_u3) header=header+' '+molecules[imol]+region_list[j] #Header for output file # #Repeat for regionStructures and overall mean RegionStructMeans=np.zeros((len(time),len(region_struct_list))) OverallMean=np.zeros(len(time)) # for itime in range(len(time)): for j in range(len(region_struct_list)): for k in region_struct_col[j]: RegionStructMeans[itime,j]+=data[itime,k] for k in range(maxvols): OverallMean[itime]+=data[itime,k] # for j in range(len(region_struct_list)): RegionStructMeans[:,j]/=(RegStructVol[j]*mol_per_nM_u3) header=header+' '+molecules[imol]+region_struct_list[j] #Header for output file # if (data[:,1:-1].all==0): OverallMean[:]/=(submembVol*mol_per_nM_u3) else: OverallMean[:]/=(TotVol*mol_per_nM_u3) header=header+' '+molecules[imol]+'AvgTot\n' # if molecules[imol] in plot_molecules: plot_index=plot_molecules.index(molecules[imol]) plot_array[:,plot_index]=OverallMean ss[fnum,plot_index]=plot_array[sstart:ssend,plot_index].mean() # #Repeat for spatial averages if specified if spatialaverage: SpatialMeans=np.zeros((len(time),bins)) for itime in range(len(time)): for j in range(bins): for k in bincolumns[j]: SpatialMeans[itime,j]+=data[itime,k] for j in range(bins): print "j, vol=", j, SpatialVol[j] if (SpatialVol[j] != 0): SpatialMeans[:,j]/=(SpatialVol[j]*mol_per_nM_u3) print SpatialMeans[1:10,j] # #write averages to separate files if outputavg: outfname=fname[0:-8]+molecules[imol]+'_avg.txt' if molecules[imol] in plot_molecules: print 'output file: ', outfname outdata=np.column_stack((time,RegionMeans,RegionStructMeans,OverallMean)) f=open(outfname, 'w') f.write(header) np.savetxt(f, outdata, fmt='%.4f', delimiter=' ') f.close() # #write space if spatialaverage: outnamespace=fname[0:-8]+'-'+molecules[imol]+'_space.txt' outdata=np.column_stack((time,SpatialMeans)) f=open(outnamespace, 'w') f.write(header+'\n') np.savetxt(f, outdata, fmt='%.4f', delimiter=' ') f.close() else: #no processing needed if only a single voxel. Just extract, calculate ss, and plot specified molecules #0 in 3 index of molecule_array indicates that for 1 voxel structures 0th array has total for imol,mol in enumerate(plot_molecules): plot_array[:,imol]=molecule_array[:,molecules.index(mol),0]/TotVol/mol_per_nM_u3 ss[fnum,imol]=plot_array[int(sstart/time[1]):int(ssend/time[1]),imol].mean() # #in both cases (single voxel and multi-voxel): #total some molecule forms - specified by hand above for now for imol,mol in enumerate(tot_species.keys()): for subspecies in tot_species[mol]: mol_sum=molecule_array[0,molecules.index(subspecies),:].sum() #print imol,mol,subspecies,molecule_array[0,molecules.index(subspecies),:],mol_sum ss_tot[fnum,imol]+=mol_sum/TotVol/mol_per_nM_u3 print imol,mol,ss_tot[fnum,imol],"nM, or in picoSD:", ss_tot[fnum,imol]*(TotVol/submembVol)*deltaY[0] #after main processing, extract a few characteristics of molecule trajectory # print params, parval[fnum] print " molecule baseline peakval ptime slope min ratio" for imol,mol in enumerate(plot_molecules): baseline[fnum,imol]=plot_array[sstart:ssend,imol].mean() peakpt=plot_array[ssend:,imol].argmax()+ssend peaktime[fnum,imol]=peakpt*dt peakval[fnum,imol]=plot_array[peakpt-10:peakpt+10,imol].mean() lowpt=plot_array[ssend:,imol].argmin()+ssend lowval[fnum,imol]=plot_array[lowpt-10:lowpt+10,imol].mean() begin_slopeval=0.2*(peakval[fnum,imol]-baseline[fnum,imol])+baseline[fnum,imol] end_slopeval=0.8*(peakval[fnum,imol]-baseline[fnum,imol])+baseline[fnum,imol] exceedsthresh=np.where(plot_array[ssend:,imol]>begin_slopeval) begin_slopept=0 end_slopept=0 found=0 if len(exceedsthresh[0]): begin_slopept=np.min(exceedsthresh)+ssend found=1 exceedsthresh=np.where(plot_array[begin_slopept:,imol]>end_slopeval) if len(exceedsthresh): end_slopept=np.min(exceedsthresh)+begin_slopept else: found=0 if found and len(plot_array[begin_slopept:end_slopept,imol])>1: slope[fnum,imol]=(peakval[fnum,imol]-baseline[fnum,imol])/((end_slopept-begin_slopept)*dt) else: slope[fnum,imol]=-9999 print mol.rjust(14),"%8.2f" % baseline[fnum,imol],"%8.2f" %peakval[fnum,imol], print "%8.2f" % peaktime[fnum,imol], "%8.3f" %slope[fnum,imol], print "%8.2f" %lowval[fnum,imol], "%8.2f" %(peakval[fnum,imol]/baseline[fnum,imol]) # #Now plot some of these molcules, either single voxel or overall average if multi-voxel # pu.plottrace(plot_molecules,time,plot_array,parval[fnum],axes,fig,col_inc,scale,minpar) # #then plot the steady state versus parameter value for each molecule if len(params)>1: print np.column_stack((parval,ss)) xval=np.zeros(len(parval)) for i,pv in enumerate(parval): if len(parlist[0])>len(parlist[1]): xval[i]=pv[0] else: xval[i]=pv[1] if showss: pu.plotss(plot_molecules,xval,ss) else: if showss: #also plot the totaled molecule forms if len(tot_species.keys()): pu.plotss(plot_molecules+tot_species.keys(),parval,np.hstack((ss,ss_tot))) else: pu.plotss(plot_molecules,parval,ss)