hooke/plugin/autopeak.py

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