neurord_analysis.py

#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)