autopeak.py

   1 #!/usr/bin/env python
   2
   3 from libhooke import WX_GOOD, ClickedPoint
   4 import wxversion
   5 wxversion.select(WX_GOOD)
   6 from wx import PostEvent
   7 import numpy as np
   8 import scipy as sp
   9 import copy
  10 import os.path
  11 import time
  12
  13 import warnings
  14 warnings.simplefilter('ignore',np.RankWarning)
  15
  16
  17 class autopeakCommands:
  18
  19     def do_autopeak(self,args):
  20         '''
  21         AUTOPEAK
  22         (autopeak.py)
  23         Automatically performs a number of analyses on the peaks of the given curve.
  24         Currently it automatically:
  25         - fits peaks with WLC function
  26         - measures peak maximum forces with a baseline
  27         - measures slope in proximity of peak maximum
  28         Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
  29
  30         Syntax:
  31         autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
  32
  33         rebase : Re-asks baseline interval
  34
  35         pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
  36                      the fit will be a 2-variable
  37                      fit. DO NOT put spaces between 'pl', '=' and the value.
  38                      The value must be in meters.
  39                      Scientific notation like 0.35e-9 is fine.
  40
  41         t=[value] : Use a user-defined temperature. The value must be in
  42                     kelvins; by default it is 293 K.
  43                     DO NOT put spaces between 't', '=' and the value.
  44
  45         noauto : allows for clicking the contact point by
  46                  hand (otherwise it is automatically estimated) the first time.
  47                  If subsequent measurements are made, the same contact point
  48                  clicked the first time is used
  49
  50         reclick : redefines by hand the contact point, if noauto has been used before
  51                   but the user is unsatisfied of the previously choosen contact point.
  52
  53         usepoints : fit interval by number of points instead than by nanometers
  54
  55         When you first issue the command, it will ask for the filename. If you are giving the filename
  56         of an existing file, autopeak will resume it and append measurements to it. If you are giving
  57         a new filename, it will create the file and append to it until you close Hooke.
  58
  59
  60         Useful variables (to set with SET command):
  61         ---
  62         temperature= temperature of the system for wlc fit (in K)
  63
  64         auto_slope_span = number of points on which measure the slope, for slope
  65
  66         auto_fit_nm = number of nm to fit before the peak maximum, for WLC (if usepoints false)
  67         auto_fit_points = number of points to fit before the peak maximum, for WLC (if usepoints true)
  68
  69         baseline_clicks = 0: automatic baseline
  70                           1: decide baseline with a single click and length defined in auto_left_baseline
  71                           2: let user click points of baseline
  72         auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
  73         auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
  74         '''
  75
  76         #MACROS.
  77         #FIXME: to move outside function
  78         def fit_interval_nm(start_index,plot,nm,backwards):
  79             '''
  80             Calculates the number of points to fit, given a fit interval in nm
  81             start_index: index of point
  82             plot: plot to use
  83             backwards: if true, finds a point backwards.
  84             '''
  85             whatset=1 #FIXME: should be decidable
  86             x_vect=plot.vectors[1][0]
  87
  88             c=0
  89             i=start_index
  90             start=x_vect[start_index]
  91             maxlen=len(x_vect)
  92             while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
  93                 if i==0 or i==maxlen-1: #we reached boundaries of vector!
  94                     return c
  95
  96                 if backwards:
  97                     i-=1
  98                 else:
  99                     i+=1
 100                 c+=1
 101             return c
 102
 103         def pickup_contact_point():
 104             '''macro to pick up the contact point by clicking'''
 105             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 106             contact_point_index=contact_point.index
 107             self.wlccontact_point=contact_point
 108             self.wlccontact_index=contact_point.index
 109             self.wlccurrent=self.current.path
 110             return contact_point, contact_point_index
 111
 112         def find_current_peaks():
 113             #Find peaks.
 114             defplot=self.current.curve.default_plots()[0]
 115             flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 116             defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 117             peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 118             return peak_location, peak_size
 119
 120         #default fit etc. variables
 121         pl_value=None
 122         T=self.config['temperature']
 123
 124         slope_span=int(self.config['auto_slope_span'])
 125         delta_force=10
 126         rebase=False #if true=we select rebase
 127
 128         #initialize output data vectors
 129         c_lengths=[]
 130         p_lengths=[]
 131         sigma_c_lengths=[]
 132         sigma_p_lengths=[]
 133         forces=[]
 134         slopes=[]
 135
 136         #pick up plot
 137         displayed_plot=self._get_displayed_plot(0)
 138
 139         #COMMAND LINE PARSING
 140         #--Using points instead of nm interval
 141         if 'usepoints' in args.split():
 142             fit_points=int(self.config['auto_fit_points'])
 143             usepoints=True
 144         else:
 145             fit_points=None
 146             usepoints=False
 147         #--Recalculate baseline
 148         if 'rebase' in args or (self.basecurrent != self.current.path):
 149             rebase=True
 150
 151         #--Custom persistent length / custom temperature
 152         for arg in args.split():
 153             #look for a persistent length argument.
 154             if 'pl=' in arg:
 155                 pl_expression=arg.split('=')
 156                 pl_value=float(pl_expression[1]) #actual value
 157             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 158             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 159                 t_expression=arg.split('=')
 160                 T=float(t_expression[1])
 161         #--Contact point arguments
 162         if 'reclick' in args.split():
 163             print 'Click contact point'
 164             contact_point, contact_point_index = pickup_contact_point()
 165         elif 'noauto' in args.split():
 166             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 167                 print 'Click contact point'
 168                 contact_point , contact_point_index = pickup_contact_point()
 169             else:
 170                 contact_point=self.wlccontact_point
 171                 contact_point_index=self.wlccontact_index
 172         else:
 173             #Automatically find contact point
 174             cindex=self.find_contact_point()
 175             contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
 176         #--END COMMAND LINE PARSING--
 177
 178
 179         peak_location, peak_size = find_current_peaks()
 180
 181         fitplot=copy.deepcopy(displayed_plot)
 182
 183         #Pick up force baseline
 184         if rebase:
 185             clicks=self.config['baseline_clicks']
 186             if clicks==0:
 187                 self.basepoints=[]
 188                 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 189                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 190                 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 191                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 192             elif clicks>0:
 193                 print 'Select baseline'
 194                 if clicks==1:
 195                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 196                     base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 197                     self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 198                 else:
 199                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 200
 201             self.basecurrent=self.current.path
 202
 203         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 204         boundaries.sort()
 205         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 206         avg=np.mean(to_average)
 207
 208
 209         for peak in peak_location:
 210             #WLC FITTING
 211             #define fit interval
 212             if not usepoints:
 213                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 214             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 215             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 216
 217             #points for the fit
 218             points=[contact_point, peak_point, other_fit_point]
 219
 220             if abs(peak_point.index-other_fit_point.index) < 2:
 221                 continue
 222
 223             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)
 224
 225
 226             #Measure forces
 227             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 228             y=min(delta_to_measure)
 229             #save force values (pN)
 230             #Measure slopes
 231             slope=self.linefit_between(peak-slope_span,peak)[0]
 232
 233
 234             #check fitted data and, if right, add peak to the measurement
 235             #FIXME: code duplication
 236             if len(params)==1: #if we did choose 1-value fit
 237                 p_lengths.append(pl_value)
 238                 c_lengths.append(params[0]*(1.0e+9))
 239                 sigma_p_lengths.append(0)
 240                 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 241                 forces.append(abs(y-avg)*(1.0e+12))
 242                 slopes.append(slope)
 243                 #Add WLC fit lines to plot
 244                 fitplot.add_set(xfit,yfit)
 245                 if len(fitplot.styles)==0:
 246                     fitplot.styles=[]
 247                 else:
 248                     fitplot.styles.append(None)
 249             else: #2-value fit
 250                 p_leng=params[1]*(1.0e+9)
 251                 #check if persistent length makes sense. otherwise, discard peak.
 252                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 253                     p_lengths.append(p_leng)
 254                     c_lengths.append(params[0]*(1.0e+9))
 255                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 256                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 257                     forces.append(abs(y-avg)*(1.0e+12))
 258                     slopes.append(slope)
 259
 260                     #Add WLC fit lines to plot
 261                     fitplot.add_set(xfit,yfit)
 262                     if len(fitplot.styles)==0:
 263                         fitplot.styles=[]
 264                     else:
 265                         fitplot.styles.append(None)
 266                         fitplot.colors.append(None)
 267                 else:
 268                     pass
 269
 270
 271         #add basepoints to fitplot
 272         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]])
 273         fitplot.styles.append('scatter')
 274
 275
 276         #Show wlc fits and peak locations
 277         self._send_plot([fitplot])
 278         #self.do_peaks('')
 279
 280         print 'Measurements for all peaks detected:'
 281         print 'contour (nm)', c_lengths
 282         print 'sigma contour (nm)',sigma_c_lengths
 283         print 'p (nm)',p_lengths
 284         print 'sigma p (nm)',sigma_p_lengths
 285         print 'forces (pN)',forces
 286         print 'slopes (N/m)',slopes
 287
 288         #Ask the user what peaks to ignore from analysis.
 289         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 290         print 'N to discard measurement'
 291         exclude_raw=raw_input('Input:')
 292         if exclude_raw=='N':
 293             print 'Discarded.'
 294             return
 295         if not exclude_raw=='':
 296             exclude=exclude_raw.split(',')
 297             try:
 298                 exclude=[int(item) for item in exclude]
 299                 for i in exclude:
 300                     c_lengths[i]=None
 301                     p_lengths[i]=None
 302                     forces[i]=None
 303                     slopes[i]=None
 304                     sigma_c_lengths[i]=None
 305                     sigma_p_lengths[i]=None
 306             except:
 307                  print 'Bad input, taking all...'
 308         #Clean data vectors from ignored peaks
 309         #FIXME:code duplication
 310         c_lengths=[item for item in c_lengths if item != None]
 311         p_lengths=[item for item in p_lengths if item != None]
 312         forces=[item for item in forces if item != None]
 313         slopes=[item for item in slopes if item != None]
 314         sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
 315         sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
 316
 317         print 'Measurements for chosen peaks:'
 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
 324
 325         #Save file info
 326         if self.autofile=='':
 327             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 328             if self.autofile=='':
 329                 print 'Not saved.'
 330                 return
 331
 332         if not os.path.exists(self.autofile):
 333             f=open(self.autofile,'w+')
 334             f.write('Analysis started '+time.asctime()+'\n')
 335             f.write('----------------------------------------\n')
 336             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
 337             f.close()
 338
 339         print 'Saving...'
 340         f=open(self.autofile,'a+')
 341
 342         f.write(self.current.path+'\n')
 343         for i in range(len(c_lengths)):
 344             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')
 345
 346         f.close()
 347         self.do_note('autopeak')
 348