hooke/plugin/autopeak.py

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