hooke/plugin/autopeak.py

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