-
Notifications
You must be signed in to change notification settings - Fork 1
/
AllRampAnal.py
147 lines (129 loc) · 6.3 KB
/
AllRampAnal.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
#This is the python program that reads in test pulse traces
#extracts baseline values, peak PSP, etc
#input is the full filename (including directory) of one of the traces, then python will find all the rest of the traces
#Outputs single file with each column containing a separate measure
#Can use Zbyszek's python program to read in IF or IV curves and extract spike measures
#run the program from within python by typing: execfile('PSPanalSA.py')
import os
from igor import binarywave
from numpy import * # whats the asterisk doing again? importing All things in numpy? Yes.
from pylab import *
from matplotlib import pyplot
import sys # if above is true, whats the dif here? ie why not from sys import *, or vice versa above?
from pprint import pprint as pp
import glob
######################### Experiment specific parameters
#MUST BE SPECIFIED:
DirPattern="ramp*_Waves"
Alldirectories=glob.glob(DirPattern)
#some parameters common to ALL experiments, From RE analylsis
#perhaps the base needs to be changed - go from the max to account for leak?
basestarttime=0.028
baseEndtime=0.046
RampStart=0.05
RampVslope=0.5 #mV/ms = V/sec
RampEnd=0.31
PeakStart=0.150
PeakEnd=0.300
maxEnd=0.15
dt=50e-6 #50 microsecond sample rate
#estimate leak from the ramp around -60 mV:
#V=-65 mV at t=0.08 sec and -55 mv at t=0.1 sec)
rampleak1=int(round(0.08/dt))
rampleak2=int(round(0.10/dt))
DeltaV=-10e-3 #10 mV
rampPts=5
RampLeakList=[]
#convert times into points
basestartpnt=basestarttime/dt
baseEndpnt=baseEndtime/dt
PeakStartpnt=PeakStart/dt
PeakEndpnt=PeakEnd/dt
maxEndpt=maxEnd/dt
for dirnum,directory in enumerate(Alldirectories):
dirparts=directory.split('_')
ExperName=dirparts[0]+'_'+dirparts[1]
filenameending="_1_*.ibw"
pattern=directory+'/'+ExperName+filenameending
print "Looking for files using: ", pattern
#########################Below here, no parameters
def sortorder(fname):
parts = fname.split('_')
ans = int(parts[-3]), int(parts[-2])
#print 'here:', fname, '->', ans
return ans
filenames = glob.glob(pattern)
if (len(filenames)==0):
print "Mistyped the filenames. Python found NO files:"
pp(filenames)
else:
filenames = sorted(filenames, key=sortorder) #sorted is alphabetical unless integers given
#key is output of function "sortorder," using argument "filenames"
#pp(filenames)
#initialize the arrays that will contain the data from all series and traces
numtraces=len(filenames)
Ihold=zeros(numtraces)
IMin=zeros(numtraces)
Amp=zeros(numtraces)
IMax=zeros(numtraces)
Amp2=zeros(numtraces)
tracestring=[]
#-0.5e9 converts to nA and leak due to 5 mV (not 10 mV) Vm change
RampLeak=zeros(numtraces)
fig=pyplot.figure(figsize=(6,6))
fig.canvas.set_window_title('Experiment '+ExperName)
axes=fig.add_subplot(111)
#loop through each trace from each series
for fileindex,eachfile in enumerate(filenames):
parts = eachfile.split('_')
tracestring.append(ExperName+'_'+parts[-3]+"_"+parts[-2]+" ") #extra spaces forces column_stack to use more significant figures when converting floats to string. There should be a better way
data = binarywave.load(eachfile) #read in data
trace=data['wave']['wData']
LeakVsTime=zeros(len(trace))
#pp(data) # optional - look at data array to see what else is there
dtfile=data['wave']['wave_header']['hsA'] #verify that the dt we have is correct
if (dt != dtfile):
print "Error, dt of file is different than expected"
#calculate various measures
Ihold[fileindex]=mean(trace[basestartpnt:baseEndpnt])
IMin[fileindex]=min(trace[PeakStartpnt:PeakEndpnt])
mintimePt=trace[PeakStartpnt:PeakEndpnt].argmin()+PeakStartpnt
mintime=mintimePt*dt
maxtimePt=trace[PeakStartpnt:maxEndpt].argmax()+PeakStartpnt
maxtime=maxtimePt*dt
IMax[fileindex]=mean(trace[maxtimePt-10:maxtimePt+10])
Amp[fileindex]=Ihold[fileindex]-IMin[fileindex]
#Calculate slope of I vs V from linear part:
DeltaI=mean(trace[rampleak1-rampPts:rampleak1+rampPts])-mean(trace[rampleak2-rampPts:rampleak2+rampPts])
LeakSlope=DeltaI/DeltaV
#-0.5e-9 converts leak into value in response to 5 mV depol
RampLeak[fileindex]=-0.5e9*(DeltaI)
#Calculate Ramp current versus time
endtime=len(trace)*dtfile
tracetime=arange(0,endtime-dt,dt)
LeakVsTime=LeakSlope*RampVslope*(tracetime-RampStart)+Ihold[fileindex]
LeakVsTime[0:int(RampStart/dt)]=Ihold[fileindex]
LeakVsTime[int(RampEnd/dt):]=Ihold[fileindex]
Amp2[fileindex]=LeakVsTime[mintimePt]-IMin[fileindex]
LeakSubtrRamp=trace-LeakVsTime
print tracestring[fileindex], "ramp Leak {0:.6f}".format( RampLeak[fileindex] )
#Optionally, plot the data to verify the results
axes.plot(tracetime,trace+fileindex*0.1e-9,label=fileindex)
axes.plot(tracetime,LeakSubtrRamp+fileindex*0.1e-9,label=fileindex)
axes.plot(tracetime,LeakVsTime)
axes.plot(mintime,LeakVsTime[mintimePt]+fileindex*0.1e-9,'ko') #label the min
axes.plot(mintime,IMin[fileindex]+fileindex*0.1e-9,'b*') #label the min
axes.legend(fontsize=8, loc='best')
fig.canvas.draw()
fig.show()
endLeak=min(5,len(RampLeak))
print "endleak", tracestring[endLeak-1],"Mean Leak", np.mean(RampLeak[0:endLeak])
RampLeakList.append(column_stack((ExperName,np.mean(RampLeak[0:endLeak]))))
#header contains wave names for igor
header='trace'+ExperName+' AmpLkSub'+ExperName+' Amp'+ExperName+' Leak'+ExperName+'\n'
outputdata=column_stack((tracestring,Amp2*1e9,Amp*1e9,RampLeak))
outfname=ExperName+"RampAmp.txt"
f=open(outfname,'w')
f.write(header)
savetxt(f, outputdata, fmt='%s', delimiter=' ')
f.close()