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         if len(peak_location) == 0:
 182             print 'No peaks to fit.'
 183             return
 184
 185         fitplot=copy.deepcopy(displayed_plot)
 186
 187         #Pick up force baseline
 188         if rebase:
 189             clicks=self.config['baseline_clicks']
 190             if clicks==0:
 191                 self.basepoints=[]
 192                 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 193                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 194                 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 195                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 196             elif clicks>0:
 197                 print 'Select baseline'
 198                 if clicks==1:
 199                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 200                     base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 201                     self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 202                 else:
 203                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 204
 205             self.basecurrent=self.current.path
 206
 207         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 208         boundaries.sort()
 209         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 210         avg=np.mean(to_average)
 211
 212
 213         for peak in peak_location:
 214             #WLC FITTING
 215             #define fit interval
 216             if not usepoints:
 217                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 218             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 219             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 220
 221             #points for the fit
 222             points=[contact_point, peak_point, other_fit_point]
 223
 224             if abs(peak_point.index-other_fit_point.index) < 2:
 225                 continue
 226
 227             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)
 228
 229
 230             #Measure forces
 231             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 232             y=min(delta_to_measure)
 233             #save force values (pN)
 234             #Measure slopes
 235             slope=self.linefit_between(peak-slope_span,peak)[0]
 236
 237
 238             #check fitted data and, if right, add peak to the measurement
 239             #FIXME: code duplication
 240             if len(params)==1: #if we did choose 1-value fit
 241                 p_lengths.append(pl_value)
 242                 c_lengths.append(params[0]*(1.0e+9))
 243                 sigma_p_lengths.append(0)
 244                 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 245                 forces.append(abs(y-avg)*(1.0e+12))
 246                 slopes.append(slope)
 247                 #Add WLC fit lines to plot
 248                 fitplot.add_set(xfit,yfit)
 249                 if len(fitplot.styles)==0:
 250                     fitplot.styles=[]
 251                     fitplot.colors=[]
 252                 else:
 253                     fitplot.styles.append(None)
 254                     fitplot.colors.append(None)
 255             else: #2-value fit
 256                 p_leng=params[1]*(1.0e+9)
 257                 #check if persistent length makes sense. otherwise, discard peak.
 258                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 259                     p_lengths.append(p_leng)
 260                     c_lengths.append(params[0]*(1.0e+9))
 261                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 262                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 263                     forces.append(abs(y-avg)*(1.0e+12))
 264                     slopes.append(slope)
 265
 266                     #Add WLC fit lines to plot
 267                     fitplot.add_set(xfit,yfit)
 268                     if len(fitplot.styles)==0:
 269                         fitplot.styles=[]
 270                         fitplot.colors=[]
 271                     else:
 272                         fitplot.styles.append(None)
 273                         fitplot.colors.append(None)
 274                 else:
 275                     pass
 276
 277
 278         #add basepoints to fitplot
 279         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]])
 280         fitplot.styles.append('scatter')
 281         fitplot.colors.append(None)
 282
 283         #Show wlc fits and peak locations
 284         self._send_plot([fitplot])
 285         #self.do_peaks('')
 286
 287         print 'Measurements for all peaks detected:'
 288         print 'contour (nm)', c_lengths
 289         print 'sigma contour (nm)',sigma_c_lengths
 290         print 'p (nm)',p_lengths
 291         print 'sigma p (nm)',sigma_p_lengths
 292         print 'forces (pN)',forces
 293         print 'slopes (N/m)',slopes
 294
 295         #Ask the user what peaks to ignore from analysis.
 296         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 297         print 'N to discard measurement'
 298         exclude_raw=raw_input('Input:')
 299         if exclude_raw=='N':
 300             print 'Discarded.'
 301             return
 302         if not exclude_raw=='':
 303             exclude=exclude_raw.split(',')
 304             try:
 305                 exclude=[int(item) for item in exclude]
 306                 for i in exclude:
 307                     c_lengths[i]=None
 308                     p_lengths[i]=None
 309                     forces[i]=None
 310                     slopes[i]=None
 311                     sigma_c_lengths[i]=None
 312                     sigma_p_lengths[i]=None
 313             except:
 314                  print 'Bad input, taking all...'
 315         #Clean data vectors from ignored peaks
 316         #FIXME:code duplication
 317         c_lengths=[item for item in c_lengths if item != None]
 318         p_lengths=[item for item in p_lengths if item != None]
 319         forces=[item for item in forces if item != None]
 320         slopes=[item for item in slopes if item != None]
 321         sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
 322         sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
 323
 324         print 'Measurements for chosen peaks:'
 325         print 'contour (nm)',c_lengths
 326         print 'sigma contour (nm)',sigma_c_lengths
 327         print 'p (nm)',p_lengths
 328         print 'sigma p (nm)',sigma_p_lengths
 329         print 'forces (pN)',forces
 330         print 'slopes (N/m)',slopes
 331
 332         #Save file info
 333         if self.autofile=='':
 334             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 335             if self.autofile=='':
 336                 print 'Not saved.'
 337                 return
 338
 339         if not os.path.exists(self.autofile):
 340             f=open(self.autofile,'w+')
 341             f.write('Analysis started '+time.asctime()+'\n')
 342             f.write('----------------------------------------\n')
 343             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
 344             f.close()
 345
 346         print 'Saving...'
 347         f=open(self.autofile,'a+')
 348
 349         f.write(self.current.path+'\n')
 350         for i in range(len(c_lengths)):
 351             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')
 352
 353         f.close()
 354         self.do_note('autopeak')
 355