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