hooke/plugin/autopeak.py

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