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 = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)\r
  70                            0: automatic baseline
  71                            1: decide baseline with a single click and length defined in auto_left_baseline
  72                            2: let user click points of baseline
  73         auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
  74         auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
  75         '''
  76
  77         #MACROS.
  78         #FIXME: to move outside function
  79         def fit_interval_nm(start_index,plot,nm,backwards):
  80             '''
  81             Calculates the number of points to fit, given a fit interval in nm
  82             start_index: index of point
  83             plot: plot to use
  84             backwards: if true, finds a point backwards.
  85             '''
  86             whatset=1 #FIXME: should be decidable
  87             x_vect=plot.vectors[1][0]
  88
  89             c=0
  90             i=start_index
  91             start=x_vect[start_index]
  92             maxlen=len(x_vect)
  93             while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
  94                 if i==0 or i==maxlen-1: #we reached boundaries of vector!
  95                     return c
  96
  97                 if backwards:
  98                     i-=1
  99                 else:
 100                     i+=1
 101                 c+=1
 102             return c
 103
 104         def pickup_contact_point():
 105             '''macro to pick up the contact point by clicking'''
 106             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 107             contact_point_index=contact_point.index
 108             self.wlccontact_point=contact_point
 109             self.wlccontact_index=contact_point.index
 110             self.wlccurrent=self.current.path\r
 111             return contact_point, contact_point_index
 112
 113         def find_current_peaks():
 114             #Find peaks.
 115             defplot=self.current.curve.default_plots()[0]
 116             flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 117             defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 118             peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 119             return peak_location, peak_size
 120
 121         #default fit etc. variables
 122         pl_value=None
 123         T=self.config['temperature']
 124
 125         slope_span=int(self.config['auto_slope_span'])
 126         delta_force=10
 127         rebase=False #if true=we select rebase
 128
 129         #initialize output data vectors
 130         c_lengths=[]
 131         p_lengths=[]
 132         sigma_c_lengths=[]
 133         sigma_p_lengths=[]
 134         forces=[]
 135         slopes=[]
 136
 137         #pick up plot
 138         displayed_plot=self._get_displayed_plot(0)
 139
 140         #COMMAND LINE PARSING
 141         #--Using points instead of nm interval
 142         if 'usepoints' in args.split():
 143             fit_points=int(self.config['auto_fit_points'])
 144             usepoints=True
 145         else:
 146             fit_points=None
 147             usepoints=False
 148         #--Recalculate baseline
 149         if 'rebase' in args or (self.basecurrent != self.current.path):
 150             rebase=True
 151
 152         #--Custom persistent length / custom temperature
 153         for arg in args.split():
 154             #look for a persistent length argument.
 155             if 'pl=' in arg:
 156                 pl_expression=arg.split('=')
 157                 pl_value=float(pl_expression[1]) #actual value
 158             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 159             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 160                 t_expression=arg.split('=')
 161                 T=float(t_expression[1])
 162         #--Contact point arguments
 163         if 'reclick' in args.split():
 164             print 'Click contact point'
 165             contact_point, contact_point_index = pickup_contact_point()
 166         elif 'noauto' in args.split():
 167             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 168                 print 'Click contact point'
 169                 contact_point , contact_point_index = pickup_contact_point()
 170             else:
 171                 contact_point=self.wlccontact_point
 172                 contact_point_index=self.wlccontact_index
 173         else:
 174             #Automatically find contact point
 175             cindex=self.find_contact_point()
 176             contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
 177         #--END COMMAND LINE PARSING--
 178
 179
 180         peak_location, peak_size = find_current_peaks()
 181
 182         if len(peak_location) == 0:
 183             print 'No peaks to fit.'
 184             return
 185
 186         fitplot=copy.deepcopy(displayed_plot)
 187
 188         #Pick up force baseline
 189         if rebase:
 190             clicks=self.config['baseline_clicks']\r
 191             if clicks==0:
 192                 self.basepoints=[]
 193                 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 194                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 195                 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 196                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 197             elif clicks>0:
 198                 print 'Select baseline'
 199                 if clicks==1:
 200                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 201                     base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 202                     self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 203                 else:
 204                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 205
 206             self.basecurrent=self.current.path
 207
 208         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 209         boundaries.sort()
 210         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 211         avg=np.mean(to_average)\r
 212         \r
 213         clicks=self.config['baseline_clicks']\r
 214         if clicks==-1:\r
 215             try:\r
 216                 avg=displayed_plot.vectors[1][1][contact_point_index]
 217             except:\r
 218                 avg=displayed_plot.vectors[1][1][cindex]
 219
 220         for peak in peak_location:
 221             #WLC FITTING
 222             #define fit interval
 223             if not usepoints:
 224                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 225             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 226             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 227
 228             #points for the fit
 229             points=[contact_point, peak_point, other_fit_point]
 230
 231             if abs(peak_point.index-other_fit_point.index) < 2:
 232                 continue
 233
 234             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)
 235
 236
 237             #Measure forces
 238             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 239             y=min(delta_to_measure)
 240             #save force values (pN)
 241             #Measure slopes
 242             slope=self.linefit_between(peak-slope_span,peak)[0]
 243
 244
 245             #check fitted data and, if right, add peak to the measurement
 246             #FIXME: code duplication
 247             if len(params)==1: #if we did choose 1-value fit
 248                 p_lengths.append(pl_value)
 249                 c_lengths.append(params[0]*(1.0e+9))
 250                 sigma_p_lengths.append(0)
 251                 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 252                 forces.append(abs(y-avg)*(1.0e+12))
 253                 slopes.append(slope)
 254                 #Add WLC fit lines to plot
 255                 fitplot.add_set(xfit,yfit)
 256                 if len(fitplot.styles)==0:
 257                     fitplot.styles=[]
 258                     fitplot.colors=[]
 259                 else:
 260                     fitplot.styles.append(None)
 261                     fitplot.colors.append(None)
 262             else: #2-value fit
 263                 p_leng=params[1]*(1.0e+9)
 264                 #check if persistent length makes sense. otherwise, discard peak.
 265                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 266                     p_lengths.append(p_leng)
 267                     c_lengths.append(params[0]*(1.0e+9))
 268                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 269                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 270                     forces.append(abs(y-avg)*(1.0e+12))
 271                     slopes.append(slope)
 272
 273                     #Add WLC fit lines to plot
 274                     fitplot.add_set(xfit,yfit)
 275                     if len(fitplot.styles)==0:
 276                         fitplot.styles=[]
 277                         fitplot.colors=[]
 278                     else:
 279                         fitplot.styles.append(None)
 280                         fitplot.colors.append(None)
 281                 else:
 282                     pass
 283
 284
 285         #add basepoints to fitplot
 286         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]])
 287         fitplot.styles.append('scatter')
 288         fitplot.colors.append(None)
 289
 290         #Show wlc fits and peak locations
 291         self._send_plot([fitplot])
 292         #self.do_peaks('')
 293
 294         print 'Measurements for all peaks detected:'
 295         print 'contour (nm)', c_lengths
 296         print 'sigma contour (nm)',sigma_c_lengths
 297         print 'p (nm)',p_lengths
 298         print 'sigma p (nm)',sigma_p_lengths
 299         print 'forces (pN)',forces
 300         print 'slopes (N/m)',slopes
 301
 302         #Ask the user what peaks to ignore from analysis.
 303         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 304         print 'N to discard measurement'
 305         exclude_raw=raw_input('Input:')
 306         if exclude_raw=='N':
 307             print 'Discarded.'
 308             return
 309         if not exclude_raw=='':
 310             exclude=exclude_raw.split(',')
 311             try:
 312                 exclude=[int(item) for item in exclude]
 313                 for i in exclude:
 314                     c_lengths[i]=None
 315                     p_lengths[i]=None
 316                     forces[i]=None
 317                     slopes[i]=None
 318                     sigma_c_lengths[i]=None
 319                     sigma_p_lengths[i]=None
 320             except:
 321                  print 'Bad input, taking all...'
 322         #Clean data vectors from ignored peaks
 323         #FIXME:code duplication
 324         c_lengths=[item for item in c_lengths if item != None]
 325         p_lengths=[item for item in p_lengths if item != None]
 326         forces=[item for item in forces if item != None]
 327         slopes=[item for item in slopes if item != None]
 328         sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
 329         sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
 330
 331         print 'Measurements for chosen peaks:'
 332         print 'contour (nm)',c_lengths
 333         print 'sigma contour (nm)',sigma_c_lengths
 334         print 'p (nm)',p_lengths
 335         print 'sigma p (nm)',sigma_p_lengths
 336         print 'forces (pN)',forces
 337         print 'slopes (N/m)',slopes
 338
 339         #Save file info
 340         if self.autofile=='':
 341             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 342             if self.autofile=='':
 343                 print 'Not saved.'
 344                 return
 345
 346         if not os.path.exists(self.autofile):
 347             f=open(self.autofile,'w+')
 348             f.write('Analysis started '+time.asctime()+'\n')
 349             f.write('----------------------------------------\n')
 350             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
 351             f.close()
 352
 353         print 'Saving...'
 354         f=open(self.autofile,'a+')
 355
 356         f.write(self.current.path+'\n')
 357         for i in range(len(c_lengths)):
 358             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')
 359
 360         f.close()
 361         self.do_note('autopeak')
 362