autopeak.py

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