autopeak.py

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