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         #MACROS.
  87         #FIXME: to move outside function
  88         def fit_interval_nm(start_index,plot,nm,backwards):
  89             '''
  90             Calculates the number of points to fit, given a fit interval in nm
  91             start_index: index of point
  92             plot: plot to use
  93             backwards: if true, finds a point backwards.
  94             '''
  95             whatset=1 #FIXME: should be decidable
  96             x_vect=plot.vectors[1][0]
  97
  98             c=0
  99             i=start_index
 100             start=x_vect[start_index]
 101             maxlen=len(x_vect)
 102             while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
 103                 if i==0 or i==maxlen-1: #we reached boundaries of vector!
 104                     return c
 105
 106                 if backwards:
 107                     i-=1
 108                 else:
 109                     i+=1
 110                 c+=1
 111             return c
 112
 113         def pickup_contact_point():
 114             '''macro to pick up the contact point by clicking'''
 115             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 116             contact_point_index=contact_point.index
 117             self.wlccontact_point=contact_point
 118             self.wlccontact_index=contact_point.index
 119             self.wlccurrent=self.current.path
 120             return contact_point, contact_point_index
 121
 122         def find_current_peaks(noflatten):
 123             #Find peaks.
 124             defplot=self.current.curve.default_plots()[0]
 125             if not noflatten:
 126                 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 127                 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 128             peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 129             return peak_location, peak_size
 130
 131         #default fit etc. variables
 132         pl_value=None
 133         T=self.config['temperature']
 134
 135         slope_span=int(self.config['auto_slope_span'])
 136         delta_force=10
 137         rebase=False #if true=we select rebase
 138         noflatten=False #if true=we avoid flattening
 139
 140         #initialize output data vectors
 141         c_lengths=[]
 142         p_lengths=[]
 143         sigma_c_lengths=[]
 144         sigma_p_lengths=[]
 145         forces=[]
 146         slopes=[]
 147
 148         #pick up plot
 149         displayed_plot=self._get_displayed_plot(0)
 150
 151         #COMMAND LINE PARSING
 152         #--Using points instead of nm interval
 153         if 'usepoints' in args.split():
 154             fit_points=int(self.config['auto_fit_points'])
 155             usepoints=True
 156         else:
 157             fit_points=None
 158             usepoints=False
 159         #--Recalculate baseline
 160         if 'rebase' in args or (self.basecurrent != self.current.path):
 161             rebase=True
 162
 163         if 'noflatten' in args:
 164             noflatten=True
 165
 166         #--Custom persistent length / custom temperature
 167         for arg in args.split():
 168             #look for a persistent length argument.
 169             if 'pl=' in arg:
 170                 pl_expression=arg.split('=')
 171                 pl_value=float(pl_expression[1]) #actual value
 172             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 173             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 174                 t_expression=arg.split('=')
 175                 T=float(t_expression[1])
 176         #--Contact point arguments
 177         if 'reclick' in args.split():
 178             print 'Click contact point'
 179             contact_point, contact_point_index = pickup_contact_point()
 180         elif 'noauto' in args.split():
 181             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 182                 print 'Click contact point'
 183                 contact_point , contact_point_index = pickup_contact_point()
 184             else:
 185                 contact_point=self.wlccontact_point
 186                 contact_point_index=self.wlccontact_index
 187         else:
 188             #Automatically find contact point
 189             cindex=self.find_contact_point()
 190             contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
 191         #--END COMMAND LINE PARSING--
 192
 193
 194         peak_location, peak_size = find_current_peaks(noflatten)
 195
 196         if len(peak_location) == 0:
 197             print 'No peaks to fit.'
 198             return
 199
 200         fitplot=copy.deepcopy(displayed_plot)
 201
 202         #Pick up force baseline
 203         if rebase:
 204             clicks=self.config['baseline_clicks']
 205             if clicks==0:
 206                 self.basepoints=[]
 207                 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 208                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 209                 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 210                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 211             elif clicks>0:
 212                 print 'Select baseline'
 213                 if clicks==1:
 214                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 215                     base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 216                     self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 217                 else:
 218                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 219
 220             self.basecurrent=self.current.path
 221
 222         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 223         boundaries.sort()
 224         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 225         avg=np.mean(to_average)
 226
 227         clicks=self.config['baseline_clicks']
 228         if clicks==-1:
 229             try:
 230                 avg=displayed_plot.vectors[1][1][contact_point_index]
 231             except:
 232                 avg=displayed_plot.vectors[1][1][cindex]
 233
 234         for peak in peak_location:
 235             #WLC FITTING
 236             #define fit interval
 237             if not usepoints:
 238                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 239             peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
 240             other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
 241
 242             #points for the fit
 243             points=[contact_point, peak_point, other_fit_point]
 244
 245             if abs(peak_point.index-other_fit_point.index) < 2:
 246                 continue
 247
 248             if self.config['fit_function']=='wlc':
 249
 250                 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)
 251             elif self.config['fit_function']=='fjc':
 252                 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)
 253             else:
 254                 print 'Unknown fit function'
 255                 print 'Please set fit_function as wlc or fjc'
 256                 return
 257
 258
 259             #Measure forces
 260             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 261             y=min(delta_to_measure)
 262             #save force values (pN)
 263             #Measure slopes
 264             slope=self.linefit_between(peak-slope_span,peak)[0]
 265
 266
 267             #check fitted data and, if right, add peak to the measurement
 268             #FIXME: code duplication
 269             if len(params)==1: #if we did choose 1-value fit
 270                 p_lengths.append(pl_value)
 271                 c_lengths.append(params[0]*(1.0e+9))
 272                 sigma_p_lengths.append(0)
 273                 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 274                 forces.append(abs(y-avg)*(1.0e+12))
 275                 slopes.append(slope)
 276                 #Add WLC fit lines to plot
 277                 fitplot.add_set(xfit,yfit)
 278                 if len(fitplot.styles)==0:
 279                     fitplot.styles=[]
 280                     fitplot.colors=[]
 281                 else:
 282                     fitplot.styles.append(None)
 283                     fitplot.colors.append(None)
 284             else: #2-value fit
 285                 p_leng=params[1]*(1.0e+9)
 286                 #check if persistent length makes sense. otherwise, discard peak.
 287                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 288                     p_lengths.append(p_leng)
 289                     c_lengths.append(params[0]*(1.0e+9))
 290                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 291                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 292                     forces.append(abs(y-avg)*(1.0e+12))
 293                     slopes.append(slope)
 294
 295                     #Add WLC fit lines to plot
 296                     fitplot.add_set(xfit,yfit)
 297                     if len(fitplot.styles)==0:
 298                         fitplot.styles=[]
 299                         fitplot.colors=[]
 300                     else:
 301                         fitplot.styles.append(None)
 302                         fitplot.colors.append(None)
 303                 else:
 304                     pass
 305
 306
 307         #add basepoints to fitplot
 308         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]])
 309         fitplot.styles.append('scatter')
 310         fitplot.colors.append(None)
 311
 312         #Show wlc fits and peak locations
 313         self._send_plot([fitplot])
 314         #self.do_peaks('')
 315
 316         print 'Using fit function: ',self.config['fit_function']
 317         print 'Measurements for all peaks detected:'
 318         print 'contour (nm)', c_lengths
 319         print 'sigma contour (nm)',sigma_c_lengths
 320         print 'p (nm)',p_lengths
 321         print 'sigma p (nm)',sigma_p_lengths
 322         print 'forces (pN)',forces
 323         print 'slopes (N/m)',slopes
 324
 325         controller=False
 326         while controller==False:
 327           #Ask the user what peaks to ignore from analysis.
 328           print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 329           print 'N to discard measurement'
 330           exclude_raw=raw_input('Input:')
 331           if exclude_raw=='N':
 332               print 'Discarded.'
 333               return
 334           if exclude_raw=='':
 335               controller=True
 336           if not exclude_raw=='':
 337               exclude=exclude_raw.split(',')
 338               try:
 339                   exclude=[int(item) for item in exclude]
 340                   for i in exclude:
 341                       c_lengths[i]=None
 342                       p_lengths[i]=None
 343                       forces[i]=None
 344                       slopes[i]=None
 345                       sigma_c_lengths[i]=None
 346                       sigma_p_lengths[i]=None
 347                       controller=True
 348               except:
 349                   print 'Bad input.'
 350                   controller=False
 351
 352         #Clean data vectors from ignored peaks
 353         #FIXME:code duplication
 354         c_lengths=[item for item in c_lengths if item != None]
 355         p_lengths=[item for item in p_lengths if item != None]
 356         forces=[item for item in forces if item != None]
 357         slopes=[item for item in slopes if item != None]
 358         sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
 359         sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
 360
 361         print 'Measurements for chosen peaks:'
 362         print 'contour (nm)',c_lengths
 363         print 'sigma contour (nm)',sigma_c_lengths
 364         print 'p (nm)',p_lengths
 365         print 'sigma p (nm)',sigma_p_lengths
 366         print 'forces (pN)',forces
 367         print 'slopes (N/m)',slopes
 368
 369         #Save file info
 370         if self.autofile=='':
 371             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 372             if self.autofile=='':
 373                 print 'Not saved.'
 374                 return
 375
 376         if not os.path.exists(self.autofile):
 377             f=open(self.autofile,'w+')
 378             f.write('Analysis started '+time.asctime()+'\n')
 379             f.write('----------------------------------------\n')
 380             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
 381             f.close()
 382
 383         print 'Saving...'
 384         f=open(self.autofile,'a+')
 385
 386         f.write(self.current.path+'\n')
 387         for i in range(len(c_lengths)):
 388             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')
 389
 390         f.close()
 391         self.do_note('autopeak')