-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnrdh5_analv2.py
206 lines (185 loc) · 9.13 KB
/
nrdh5_analv2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 23 10:58:05 2020
@author: kblackw1
#from outside python, type python nrdh5_analv2,py subdir/fileroot -par par1 par2 -mol mol1 mol2 -start 100 200 -tot tot_species_file
from within python, type
ARGS="subdir/fileroot -par par1 par2 -mol mol1 mol2 -start 100 200 -tot tot_species_file"
exec(open(('path/to/file/nrdh5_anal.py').read())
#-par:
# par1 and optionally par2 are specifications of parameter variations, as follows:
# The filenames to read in are constructed as "subdir/fileroot"+"-"+par1+"*"-"+par2+"*"
# DO NOT use hyphens in filenames except for preceding parameter name
#-mol:
# mol1 mol2, etc are the names of molecules to process
# if mol ommitted, then all molecules processed.
#-start:
# start and end time for calculating basal concentration, in sec
# if sstart ssend are ommitted, then calculates basal from 7.5-10% of runtime
#-tot:
# tot_species_files containes a list of molecule forms to total, e.g. pPDE10 and pPDE10cAMP to calculate total pPDE10
#if no parameters specified, then fileroot needs to be full filename (excluding the extension)
#e.g. ARGS='Model_Cof -par HSJCF4trains -mol Cof pCof Cofactin -tot tot_species'
#e.g. ARGS='Model_Cof-HSJCF4trainsspacedcrtl -mol Cof pCof Cofactin -tot tot_species'
additional parameters lines 27-48
"""
#ARGS='/local/vol00/Users/klblackwell/sigpath/cofilin/Model_Cof-4trains -par massed -mol Cof pCof Cofactin PKAc -tot /local/vol00/Users/klblackwell/sigpath/cofilin/tot_species'
import numpy as np
import sys
import plot_h5V2 as pu5
from nrd_output import nrdh5_output
from nrd_group import nrdh5_group
from h5utilsV2 import parse_args,get_tot
#probably should add most of these to args with defaults
submembname='sub'
dendname="dend"
spinehead="head"
stimspine=['sa1[0]'] #list of stimulated spines
spatial_bins=0 #number of spatial bins to subdivide dendrite to look at spatial gradients
window_size=0.1 #number of msec on either side of peak value to average for maximum
#These control what output is printed or written
show_inject=0
write_output=0#one file per molecules per input file
output_auc=0#one file per molecule per set of input files
showplot=3 #0 for none, 1 for overall average, 2 for spine concentration, 3 for spine and nonspine on seperate graph, or for a region plot when there are no spines
show_mol_totals=0
print_head_stats=0
textsize=8
feature_list=[]#['auc','duration']#['duration','auc']
#these molecules MUST be specified as plot_molecules
mol_pairs=[] #[['pCof','actCof'],['CKpCamCa4','PKAphos']]#[['pCof','RacPAK']]#[['CKpCamCa4','ppERK']]#,['ppERK','pSynGap']]
pairs_timeframe=[100,600]#[200,2000] #units are sec
pair_region=['dend','sa1[0]'] #plot mol pairs in which region?
basestart_time=0#2200 #make this value 0 to calculate AUC using baseline calculated from initial time period
aucend=None#600#end time for calculating auc, or None to calculate end time automatically
save_fig=True
#write_trials=True
############## END OF PARAMETERS #################
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:]
import os
mydir=os.getcwd()
sys.path.append(mydir)
print("commandline =", args)
do_exit = True
params=parse_args(args,do_exit)
if params.mol:
plot_molecules=params.mol
else:
plot_molecules=None
figtitle=params.fileroot.split('/')[-1]
if params.par:
figtitle+=' '.join(params.par)
tot_species,weight,sub_species,signature,thresh,min_max=get_tot(params)
num_LTP_stim=params.num_stim
iti=params.iti
og=nrdh5_group(params.fileroot,params.par,tot_species)
for fnum,ftuple in enumerate(og.ftuples):
data=nrdh5_output(ftuple)
data.rows_and_columns(plot_molecules,params.start,params.end)
data.molecule_population()
#print(data.data['model']['grid'][:])
if data.maxvols>1:
data.region_structures(dendname,submembname,spinehead,stimspine) #stimspine is optional
if spatial_bins>0:
data.spatial_structures(spatial_bins,dendname)
data.average_over_voxels()
# need to add another total array for different regions (to use for signature)
#Default outputset is _main_, can specify outset= something
data.total_subspecies(tot_species,sub_species,params.start,weights=weight)
og.conc_arrays(data)
#Now, print or write some optional outputs
if write_output:
data.write_average()
if 'event_statistics' in data.data['trial0']['output'].keys() and show_inject:
print ("seeds", data.seeds," injection stats:")
print('molecule '+' '.join(data.trials))
for imol,inject_sp in enumerate(data.data['model']['event_statistics'][:]):
inject_num=np.column_stack([data.data[t]['output']['event_statistics'][0] for t in data.trials])
print (inject_sp.split()[-1].rjust(20),inject_num[imol])
if print_head_stats:
data.print_head_stats()
#extract some features from the group of data files
#Default numstim = 1, so that parameter not needed for single pulse
#other parameter defaults: lo_thresh_factor=0.2,hi_thresh_factor=0.8, std_factor=2
#another parameter default: end_baseline_start=0 (uses initial baseline to calculate auc).
#Specify specific sim time near end of sim if initialization not sufficient for a good baseline for auc calculation
og.trace_features(window_size,std_factor=2,numstim=num_LTP_stim,end_baseline_start=basestart_time,filt_length=31,aucend=aucend,iti=iti)
if len(feature_list) and output_auc:
og.write_features(feature_list,params.fileroot,params.write_trials)
#################
#print all the features in nice format - Overall values only
regnum=0 #change this to print other regions
features=[k[0:7] for k in og.feature_dict.keys()]
print("file-params molecule " +' '.join(features)+' ratio CV')
for fnum,ftuple in enumerate(og.ftuples):
for imol,mol in enumerate(list(og.molecules)+og.tot_species):
outputvals=[str(np.round(odict[imol,fnum,regnum],3)) for feat,odict in og.mean_feature.items()]
if og.mean_feature['baseline'][imol,fnum,regnum]>1e-5:
outputvals.append(str(np.round(og.mean_feature['amplitude'][imol,fnum,regnum]/og.mean_feature['baseline'][imol,fnum,regnum],2)).rjust(8))
else:
outputvals.append('inf')
print(ftuple[1],mol.rjust(16),' ',' '.join(outputvals),np.std(og.feature_dict['auc'][imol,fnum])/og.mean_feature['auc'][imol,fnum])
######################### Plots
if showplot:
fig,col_inc,scale=pu5.plot_setup(data.molecules,og,len(data.spinelist),showplot)
if showplot==2 and len(stimspine):
figtitle=figtitle+' '+' '.join(stimspine)
if showplot==3:
for spnum,sp in enumerate(data.spinelist):
fig[spnum].suptitle(figtitle+' '+sp)
else:
fig.canvas.manager.set_window_title(figtitle)
pu5.plottrace(data.molecules,og,fig,col_inc,scale,data.spinelist,showplot,textsize=textsize)
if showplot==3 and data.maxvols>1 and len(data.spinelist)==0:
fig2,col_inc,scale=pu5.plot_setup(data.molecules,og,len(data.region_dict),3)
for regnum,reg in enumerate(data.region_dict):
fig2[regnum].suptitle(figtitle+' '+reg)
pu5.plotregions(data.molecules,og,fig2,col_inc,scale,data.region_dict,textsize=textsize)
#also plot the totaled molecule forms
if len(tot_species):
figtot=pu5.plot_total_mol(tot_species,og,figtitle,col_inc,textsize=textsize,regions=['sa1[0]','dendsub'])
for feat in feature_list:
pu5.plot_features(og,feat,figtitle)
if spatial_bins and data.maxvols>1:
pu5.spatial_plot(data,og)
if len(mol_pairs):
pu5.pairs(og,mol_pairs,pairs_timeframe,pair_region)
if len(signature):
og.norm_sig(signature,thresh,min_max)
figsig=pu5.plot_signature(og,thresh,figtitle,col_inc,textsize=textsize) #plot some feature values
for feature in og.sig_features.keys():
print('FEATURE:',feature)
for key in og.sig_features[feature].keys():
print (key,':', og.sig_features[feature][key])
if write_output:
write_sig(og)
if save_fig==True:
for f,sp in zip(fig,data.spinelist):
f.set_size_inches((13,13))
f.savefig(figtitle+'_'+sp+'.png')
if showplot==3 and data.maxvols>1 and len(data.spinelist)==0:
for f,reg in zip(fig2,data.region_dict):
f.savefig(figtitle+'_'+reg+'.png')
if len(tot_species):
figtot.savefig(figtitle+'tot.png')
if len(signature):
figsig.savefig(figtitle+'sig.png')
'''
2. Possibly bring in signature code from sig.py or sig2.py and eliminate one or both of those.
pu5.plot3D is used for signatures in sig.py. How does this differ from spatial_plot?
3. possibly calculate some feature value relative to a control group (e.g. auc_ratio)
4. possibly try to figure out how to extract number of stimuli
inj_keyword='<onset>'
#lookup - how to decode b'<onset>46000.0</onset>'
#injtime=
#duration=b'<duration>1000.0</duration>'
#stop = onset+duration
#sort by injtime
#check whether injtime<=stop - if so, separate pulse? if not continuous
######### Test signature using weight either in total_species or signature
'''