hooke/plugin/autopeak.py

   1 # Copyright (C) 2008-2012 Fabrizio Benedetti <fabrizio.benedetti.82@gmail.com>
   2 #                         Marco Brucale <marco.brucale@unibo.it>
   3 #                         Massimo Sandal <devicerandom@gmail.com>
   4 #                         W. Trevor King <wking@tremily.us>
   5 #
   6 # This file is part of Hooke.
   7 #
   8 # Hooke is free software: you can redistribute it and/or modify it under the
   9 # terms of the GNU Lesser General Public License as published by the Free
  10 # Software Foundation, either version 3 of the License, or (at your option) any
  11 # later version.
  12 #
  13 # Hooke is distributed in the hope that it will be useful, but WITHOUT ANY
  14 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  15 # A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  16 # details.
  17 #
  18 # You should have received a copy of the GNU Lesser General Public License
  19 # along with Hooke.  If not, see <http://www.gnu.org/licenses/>.
  20
  21 """The autopeak module provides :class:`Autopeak`, a
  22 :class:`hooke.plugin.Plugin` for automatically extracting force peaks
  23 (unfolding events) from force curves.
  24 """
  25
  26 from hooke.libhooke import WX_GOOD
  27
  28 import wxversion
  29 wxversion.select(WX_GOOD)
  30 from wx import PostEvent
  31 import numpy as np
  32 import scipy as sp
  33 import copy
  34 import os.path
  35 import time
  36
  37 import warnings
  38 warnings.simplefilter('ignore',np.RankWarning)
  39
  40 #from .. import ui.gui.results as results
  41
  42
  43 class autopeakCommands(object):
  44     '''
  45     Autopeak carries out force curve fitting with a chosen model:
  46         - WLC
  47         - FJC
  48         - FJC-PEG
  49     '''
  50     def do_autopeak(self, args):
  51         '''
  52         AUTOPEAK
  53         (autopeak.py)
  54         Automatically performs a number of analyses on the peaks of the given curve.
  55         Currently it automatically:
  56         - fits peaks with WLC or FJC function (depending on how the fit_function variable is set)
  57         - measures peak maximum forces with a baseline
  58         - measures slope in proximity of peak maximum
  59         Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
  60
  61         Syntax:
  62         autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
  63
  64         rebase : Re-asks baseline interval
  65
  66         pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
  67                      the fit will be a 2-variable
  68                      fit. DO NOT put spaces between 'pl', '=' and the value.
  69                      The value must be in meters.
  70                      Scientific notation like 0.35e-9 is fine.
  71
  72         t=[value] : Use a user-defined temperature. The value must be in
  73                     kelvins; by default it is 293 K.
  74                     DO NOT put spaces between 't', '=' and the value.
  75
  76         noauto : allows for clicking the contact point by
  77                  hand (otherwise it is automatically estimated) the first time.
  78                  If subsequent measurements are made, the same contact point
  79                  clicked the first time is used
  80
  81         reclick : redefines by hand the contact point, if noauto has been used before
  82                   but the user is unsatisfied of the previously choosen contact point.
  83
  84         usepoints : fit interval by number of points instead than by nanometers
  85
  86         noflatten : does not use the "flatten" plot manipulator
  87
  88         When you first issue the command, it will ask for the filename. If you are giving the filename
  89         of an existing file, autopeak will resume it and append measurements to it. If you are giving
  90         a new filename, it will create the file and append to it until you close Hooke.
  91
  92
  93         Useful variables (to set with SET command):
  94         ---
  95         fit_function = type of function to use for elasticity. If "wlc" worm-like chain is used, if "fjc" freely jointed
  96                        chain is used
  97
  98         temperature= temperature of the system for wlc/fjc fit (in K)
  99
 100         auto_slope_span = number of points on which measure the slope, for slope
 101
 102         auto_fit_nm = number of nm to fit before the peak maximum, for WLC/FJC (if usepoints false)
 103         auto_fit_points = number of points to fit before the peak maximum, for WLC/FJC (if usepoints true)
 104
 105         baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
 106                            0: automatic baseline
 107                            1: decide baseline with a single click and length defined in auto_left_baseline
 108                            2: let user click points of baseline
 109         auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
 110         auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
 111
 112         auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
 113                                 outside of which the peak is automatically discarded (in nm)
 114         '''
 115
 116         #default fit etc. variables
 117         pl_value=None
 118         T=self.config['temperature']
 119
 120         slope_span=int(self.config['auto_slope_span'])
 121         delta_force=10
 122         rebase=False #if true=we select rebase
 123         noflatten=False #if true=we avoid flattening
 124
 125         #initialize output data vectors
 126         c_lengths=[]
 127         p_lengths=[]
 128         sigma_c_lengths=[]
 129         sigma_p_lengths=[]
 130         forces=[]
 131         slopes=[]
 132
 133         #pick up plot
 134         displayed_plot=self._get_displayed_plot(0)
 135
 136         #COMMAND LINE PARSING
 137         #--Using points instead of nm interval
 138         if 'usepoints' in args.split():
 139             fit_points=int(self.config['auto_fit_points'])
 140             usepoints=True
 141         else:
 142             fit_points=None
 143             usepoints=False
 144         #--Recalculate baseline
 145         if 'rebase' in args or (self.basecurrent != self.current.path):
 146             rebase=True
 147
 148         if 'noflatten' in args:
 149             noflatten=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 = self.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 = self.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 = self.find_current_peaks(noflatten)
 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             self.basepoints=self.baseline_points(peak_location, displayed_plot)
 190
 191         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 192         boundaries.sort()
 193         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 194         avg=np.mean(to_average)
 195
 196         clicks=self.config['baseline_clicks']
 197         if clicks==-1:
 198             try:
 199                 avg=displayed_plot.vectors[1][1][contact_point_index]
 200             except:
 201                 avg=displayed_plot.vectors[1][1][cindex]
 202
 203         for peak in peak_location:
 204             #WLC FITTING
 205             #define fit interval
 206             if not usepoints:
 207                 fit_points=self.fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 208             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 209             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 210
 211             #points for the fit
 212             points=[contact_point, peak_point, other_fit_point]
 213
 214             if abs(peak_point.index-other_fit_point.index) < 2:
 215                 continue
 216
 217             if self.config['fit_function']=='wlc':
 218
 219                 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)
 220             elif self.config['fit_function']=='fjc':
 221                 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)
 222             else:
 223                 print 'Unknown fit function'
 224                 print 'Please set fit_function as wlc or fjc'
 225                 return
 226
 227
 228             #Measure forces
 229             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 230             y=min(delta_to_measure)
 231             #save force values (pN)
 232             #Measure slopes
 233             slope=self.linefit_between(peak-slope_span,peak)[0]
 234
 235
 236             #check fitted data and, if right, add peak to the measurement
 237             #FIXME: code duplication
 238             if len(params)==1: #if we did choose 1-value fit
 239                 p_lengths.append(pl_value)
 240                 c_lengths.append(params[0]*(1.0e+9))
 241                 sigma_p_lengths.append(0)
 242                 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 243                 forces.append(abs(y-avg)*(1.0e+12))
 244                 slopes.append(slope)
 245                 #Add WLC fit lines to plot
 246                 fitplot.add_set(xfit,yfit)
 247                 if len(fitplot.styles)==0:
 248                     fitplot.styles=[]
 249                     fitplot.colors=[]
 250                 else:
 251                     fitplot.styles.append(None)
 252                     fitplot.colors.append(None)
 253             else: #2-value fit
 254                 p_leng=params[1]*(1.0e+9)
 255                 #check if persistent length makes sense. otherwise, discard peak.
 256                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 257                     p_lengths.append(p_leng)
 258                     c_lengths.append(params[0]*(1.0e+9))
 259                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 260                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 261                     forces.append(abs(y-avg)*(1.0e+12))
 262                     slopes.append(slope)
 263
 264                     #Add WLC fit lines to plot
 265                     fitplot.add_set(xfit,yfit)
 266                     if len(fitplot.styles)==0:
 267                         fitplot.styles=[]
 268                         fitplot.colors=[]
 269                     else:
 270                         fitplot.styles.append(None)
 271                         fitplot.colors.append(None)
 272                 else:
 273                     pass
 274
 275
 276         #add basepoints to fitplot
 277         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]])
 278         fitplot.styles.append('scatter')
 279         fitplot.colors.append(None)
 280
 281         #Show wlc fits and peak locations
 282         self._send_plot([fitplot])
 283
 284         print 'Using fit function: ',self.config['fit_function']
 285         print 'Measurements for all peaks detected:'
 286         print 'contour (nm)', c_lengths
 287         print 'sigma contour (nm)',sigma_c_lengths
 288         print 'p (nm)',p_lengths
 289         print 'sigma p (nm)',sigma_p_lengths
 290         print 'forces (pN)',forces
 291         print 'slopes (N/m)',slopes
 292
 293         controller=False
 294         while controller==False:
 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 exclude_raw=='':
 303               controller=True
 304           if not exclude_raw=='':
 305               exclude=exclude_raw.split(',')
 306               try:
 307                   exclude=[int(item) for item in exclude]
 308                   for i in exclude:
 309                       c_lengths[i]=None
 310                       p_lengths[i]=None
 311                       forces[i]=None
 312                       slopes[i]=None
 313                       sigma_c_lengths[i]=None
 314                       sigma_p_lengths[i]=None
 315                       controller=True
 316               except:
 317                   print 'Bad input.'
 318                   controller=False
 319
 320         #Clean data vectors from ignored peaks
 321         #FIXME:code duplication
 322         c_lengths=[item for item in c_lengths if item != None]
 323         p_lengths=[item for item in p_lengths if item != None]
 324         forces=[item for item in forces if item != None]
 325         slopes=[item for item in slopes if item != None]
 326         sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
 327         sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
 328
 329         print 'Measurements for chosen peaks:'
 330         print 'contour (nm)',c_lengths
 331         print 'sigma contour (nm)',sigma_c_lengths
 332         print 'p (nm)',p_lengths
 333         print 'sigma p (nm)',sigma_p_lengths
 334         print 'forces (pN)',forces
 335         print 'slopes (N/m)',slopes
 336
 337         #Save file info
 338         if self.autofile=='':
 339             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 340             if self.autofile=='':
 341                 print 'Not saved.'
 342                 return
 343
 344         if not os.path.exists(self.autofile):
 345             f=open(self.autofile,'w+')
 346             f.write('Analysis started '+time.asctime()+'\n')
 347             f.write('----------------------------------------\n')
 348             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
 349             f.close()
 350
 351         print 'Saving...'
 352         f=open(self.autofile,'a+')
 353
 354         f.write(self.current.path+'\n')
 355         for i in range(len(c_lengths)):
 356             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')
 357
 358         f.close()
 359         self.do_note('autopeak')
 360