hooke/plugin/autopeak.py

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