hooke/plugin/autopeak.py

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