1 # -*- coding: utf-8 -*-
3 from hooke.libhooke import WX_GOOD, ClickedPoint
6 wxversion.select(WX_GOOD)
7 from wx import PostEvent
15 warnings.simplefilter('ignore',np.RankWarning)
18 class autopeakCommands(object):
20 def do_autopeak(self,args):
24 Automatically performs a number of analyses on the peaks of the given curve.
25 Currently it automatically:
26 - fits peaks with WLC or FJC function (depending on how the fit_function variable is set)
27 - measures peak maximum forces with a baseline
28 - measures slope in proximity of peak maximum
29 Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
32 autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
34 rebase : Re-asks baseline interval
36 pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
37 the fit will be a 2-variable
38 fit. DO NOT put spaces between 'pl', '=' and the value.
39 The value must be in meters.
40 Scientific notation like 0.35e-9 is fine.
42 t=[value] : Use a user-defined temperature. The value must be in
43 kelvins; by default it is 293 K.
44 DO NOT put spaces between 't', '=' and the value.
46 noauto : allows for clicking the contact point by
47 hand (otherwise it is automatically estimated) the first time.
48 If subsequent measurements are made, the same contact point
49 clicked the first time is used
51 reclick : redefines by hand the contact point, if noauto has been used before
52 but the user is unsatisfied of the previously choosen contact point.
54 usepoints : fit interval by number of points instead than by nanometers
56 noflatten : does not use the "flatten" plot manipulator
58 When you first issue the command, it will ask for the filename. If you are giving the filename
59 of an existing file, autopeak will resume it and append measurements to it. If you are giving
60 a new filename, it will create the file and append to it until you close Hooke.
63 Useful variables (to set with SET command):
65 fit_function = type of function to use for elasticity. If "wlc" worm-like chain is used, if "fjc" freely jointed
68 temperature= temperature of the system for wlc/fjc fit (in K)
70 auto_slope_span = number of points on which measure the slope, for slope
72 auto_fit_nm = number of nm to fit before the peak maximum, for WLC/FJC (if usepoints false)
73 auto_fit_points = number of points to fit before the peak maximum, for WLC/FJC (if usepoints true)
75 baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
77 1: decide baseline with a single click and length defined in auto_left_baseline
78 2: let user click points of baseline
79 auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
80 auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
82 auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
83 outside of which the peak is automatically discarded (in nm)
87 #FIXME: to move outside function
88 def fit_interval_nm(start_index,plot,nm,backwards):
90 Calculates the number of points to fit, given a fit interval in nm
91 start_index: index of point
93 backwards: if true, finds a point backwards.
95 whatset=1 #FIXME: should be decidable
96 x_vect=plot.vectors[1][0]
100 start=x_vect[start_index]
102 while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
103 if i==0 or i==maxlen-1: #we reached boundaries of vector!
113 def pickup_contact_point():
114 '''macro to pick up the contact point by clicking'''
115 contact_point=self._measure_N_points(N=1, whatset=1)[0]
116 contact_point_index=contact_point.index
117 self.wlccontact_point=contact_point
118 self.wlccontact_index=contact_point.index
119 self.wlccurrent=self.current.path
120 return contact_point, contact_point_index
122 def find_current_peaks(noflatten):
124 defplot=self.current.curve.default_plots()[0]
126 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
127 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
128 peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
129 return peak_location, peak_size
131 #default fit etc. variables
133 T=self.config['temperature']
135 slope_span=int(self.config['auto_slope_span'])
137 rebase=False #if true=we select rebase
138 noflatten=False #if true=we avoid flattening
140 #initialize output data vectors
149 displayed_plot=self._get_displayed_plot(0)
151 #COMMAND LINE PARSING
152 #--Using points instead of nm interval
153 if 'usepoints' in args.split():
154 fit_points=int(self.config['auto_fit_points'])
159 #--Recalculate baseline
160 if 'rebase' in args or (self.basecurrent != self.current.path):
163 if 'noflatten' in args:
166 #--Custom persistent length / custom temperature
167 for arg in args.split():
168 #look for a persistent length argument.
170 pl_expression=arg.split('=')
171 pl_value=float(pl_expression[1]) #actual value
172 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
173 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
174 t_expression=arg.split('=')
175 T=float(t_expression[1])
176 #--Contact point arguments
177 if 'reclick' in args.split():
178 print 'Click contact point'
179 contact_point, contact_point_index = pickup_contact_point()
180 elif 'noauto' in args.split():
181 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
182 print 'Click contact point'
183 contact_point , contact_point_index = pickup_contact_point()
185 contact_point=self.wlccontact_point
186 contact_point_index=self.wlccontact_index
188 #Automatically find contact point
189 cindex=self.find_contact_point()
190 contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
191 #--END COMMAND LINE PARSING--
194 peak_location, peak_size = find_current_peaks(noflatten)
196 if len(peak_location) == 0:
197 print 'No peaks to fit.'
200 fitplot=copy.deepcopy(displayed_plot)
202 #Pick up force baseline
204 clicks=self.config['baseline_clicks']
207 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
208 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
209 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
210 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
212 print 'Select baseline'
214 self.basepoints=self._measure_N_points(N=1, whatset=whatset)
215 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
216 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
218 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
220 self.basecurrent=self.current.path
222 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
224 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
225 avg=np.mean(to_average)
227 clicks=self.config['baseline_clicks']
230 avg=displayed_plot.vectors[1][1][contact_point_index]
232 avg=displayed_plot.vectors[1][1][cindex]
234 for peak in peak_location:
238 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
239 peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
240 other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
243 points=[contact_point, peak_point, other_fit_point]
245 if abs(peak_point.index-other_fit_point.index) < 2:
248 if self.config['fit_function']=='wlc':
250 params, yfit, xfit, fit_errors = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T, return_errors=True)
251 elif self.config['fit_function']=='fjc':
252 params, yfit, xfit, fit_errors = self.fjc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T, return_errors=True)
254 print 'Unknown fit function'
255 print 'Please set fit_function as wlc or fjc'
260 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
261 y=min(delta_to_measure)
262 #save force values (pN)
264 slope=self.linefit_between(peak-slope_span,peak)[0]
267 #check fitted data and, if right, add peak to the measurement
268 #FIXME: code duplication
269 if len(params)==1: #if we did choose 1-value fit
270 p_lengths.append(pl_value)
271 c_lengths.append(params[0]*(1.0e+9))
272 sigma_p_lengths.append(0)
273 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
274 forces.append(abs(y-avg)*(1.0e+12))
276 #Add WLC fit lines to plot
277 fitplot.add_set(xfit,yfit)
278 if len(fitplot.styles)==0:
282 fitplot.styles.append(None)
283 fitplot.colors.append(None)
285 p_leng=params[1]*(1.0e+9)
286 #check if persistent length makes sense. otherwise, discard peak.
287 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
288 p_lengths.append(p_leng)
289 c_lengths.append(params[0]*(1.0e+9))
290 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
291 sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
292 forces.append(abs(y-avg)*(1.0e+12))
295 #Add WLC fit lines to plot
296 fitplot.add_set(xfit,yfit)
297 if len(fitplot.styles)==0:
301 fitplot.styles.append(None)
302 fitplot.colors.append(None)
307 #add basepoints to fitplot
308 fitplot.add_set([self.basepoints[0].graph_coords[0],self.basepoints[1].graph_coords[0]],[self.basepoints[0].graph_coords[1],self.basepoints[1].graph_coords[1]])
309 fitplot.styles.append('scatter')
310 fitplot.colors.append(None)
312 #Show wlc fits and peak locations
313 self._send_plot([fitplot])
316 print 'Using fit function: ',self.config['fit_function']
317 print 'Measurements for all peaks detected:'
318 print 'contour (nm)', c_lengths
319 print 'sigma contour (nm)',sigma_c_lengths
320 print 'p (nm)',p_lengths
321 print 'sigma p (nm)',sigma_p_lengths
322 print 'forces (pN)',forces
323 print 'slopes (N/m)',slopes
326 while controller==False:
327 #Ask the user what peaks to ignore from analysis.
328 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
329 print 'N to discard measurement'
330 exclude_raw=raw_input('Input:')
336 if not exclude_raw=='':
337 exclude=exclude_raw.split(',')
339 exclude=[int(item) for item in exclude]
345 sigma_c_lengths[i]=None
346 sigma_p_lengths[i]=None
352 #Clean data vectors from ignored peaks
353 #FIXME:code duplication
354 c_lengths=[item for item in c_lengths if item != None]
355 p_lengths=[item for item in p_lengths if item != None]
356 forces=[item for item in forces if item != None]
357 slopes=[item for item in slopes if item != None]
358 sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
359 sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
361 print 'Measurements for chosen peaks:'
362 print 'contour (nm)',c_lengths
363 print 'sigma contour (nm)',sigma_c_lengths
364 print 'p (nm)',p_lengths
365 print 'sigma p (nm)',sigma_p_lengths
366 print 'forces (pN)',forces
367 print 'slopes (N/m)',slopes
370 if self.autofile=='':
371 self.autofile=raw_input('Autopeak filename? (return to ignore) ')
372 if self.autofile=='':
376 if not os.path.exists(self.autofile):
377 f=open(self.autofile,'w+')
378 f.write('Analysis started '+time.asctime()+'\n')
379 f.write('----------------------------------------\n')
380 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) ; Sigma contour (nm) ; Sigma persistence (nm)\n')
384 f=open(self.autofile,'a+')
386 f.write(self.current.path+'\n')
387 for i in range(len(c_lengths)):
388 f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+' ; '+str(sigma_c_lengths[i])+' ; '+str(sigma_p_lengths[i])+'\n')
391 self.do_note('autopeak')