plugins/autopeak.py

   1 #!/usr/bin/env python
   2
   3 '''
   4 autopeak.py
   5
   6 Automatic peak detection and analysis.
   7
   8 Copyright ???? by ?
   9 with modifications by Dr. Rolf Schmidt (Concordia University, Canada)
  10
  11 This program is released under the GNU General Public License version 2.
  12 '''
  13
  14 import lib.libhooke as lh
  15 import wxversion
  16 wxversion.select(lh.WX_GOOD)
  17
  18 import copy
  19 from numpy import mean, RankWarning
  20
  21 import warnings
  22 warnings.simplefilter('ignore', RankWarning)
  23
  24 #import config
  25 import lib.plugin
  26 import lib.results
  27
  28 class autopeakCommands:
  29     '''
  30     Autopeak carries out force curve fitting with a chosen model:
  31         - WLC
  32         - FJC
  33         - FJC-PEG
  34     '''
  35
  36     def do_autopeak(self, plot=None):
  37         '''
  38         AUTOPEAK
  39         (autopeak.py)
  40         Automatically performs a number of analyses on the peaks of the given curve.
  41         Currently it automatically:
  42         - fits peaks with WLC or FJC function (depending on how the fit_function variable is set)
  43         - measures peak maximum forces with a baseline
  44         - measures slope in proximity of peak maximum
  45         Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
  46
  47         Syntax:
  48         autopeak [rebase] [persistence_length=value] [t=value] [noauto] [reclick]
  49
  50         rebase : Re-asks baseline interval
  51
  52         persistence_length=[value] : Use a fixed persistent length for the fit. If persistence_length is not given,
  53                      the fit will be a 2-variable
  54                      fit. DO NOT put spaces between 'persistence_length', '=' and the value.
  55                      The value must be in nanometers.
  56                      Scientific notation like 0.35 is fine.
  57
  58         t=[value] : Use a user-defined temperature. The value must be in
  59                     kelvins; by default it is 293 K.
  60                     DO NOT put spaces between 't', '=' and the value.
  61
  62         noauto : allows for clicking the contact point by
  63                  hand (otherwise it is automatically estimated) the first time.
  64                  If subsequent measurements are made, the same contact point
  65                  clicked the first time is used
  66
  67         reclick : redefines by hand the contact point, if noauto has been used before
  68                   but the user is unsatisfied of the previously choosen contact point.
  69
  70         usepoints : fit interval by number of points instead than by nanometers
  71
  72         noflatten : does not use the "flatten" plot manipulator
  73
  74         When you first issue the command, it will ask for the filename. If you are giving the filename
  75         of an existing file, autopeak will resume it and append measurements to it. If you are giving
  76         a new filename, it will create the file and append to it until you close Hooke.
  77
  78
  79         Useful variables (to set with SET command):
  80         ---
  81         fit_function = type of function to use for elasticity. If "wlc" worm-like chain is used, if "fjc" freely jointed
  82                        chain is used
  83
  84         temperature= temperature of the system for wlc/fjc fit (in K)
  85
  86         auto_slope_span = number of points on which measure the slope, for slope
  87
  88         auto_fit_nm = number of nm to fit before the peak maximum, for WLC/FJC (if usepoints false)
  89         auto_fit_points = number of points to fit before the peak maximum, for WLC/FJC (if usepoints true)
  90
  91         baseline_clicks = contact point: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
  92                           automatic:     automatic baseline
  93                           1 point:       decide baseline with a single click and length defined in auto_left_baseline
  94                           2 points:      let user click points of baseline
  95         auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks = automatic , 1 point)
  96         auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = automatic)
  97
  98         auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
  99                                   outside of which the peak is automatically discarded (in nm)
 100         '''
 101
 102         #default fit etc. variables
 103         auto_fit_nm = self.GetFloatFromConfig('autopeak', 'auto_fit_nm')
 104         auto_left_baseline = self.GetFloatFromConfig('autopeak', 'auto_left_baseline')
 105         auto_max_p = self.GetFloatFromConfig('autopeak', 'auto_max_p')
 106         auto_min_p = self.GetFloatFromConfig('autopeak', 'auto_min_p')
 107         auto_right_baseline = self.GetFloatFromConfig('autopeak', 'auto_right_baseline')
 108         baseline_clicks = self.GetStringFromConfig('autopeak', 'baseline_clicks')
 109         fit_function = self.GetStringFromConfig('autopeak', 'fit_function')
 110         fit_points = self.GetIntFromConfig('autopeak', 'auto_fit_points')
 111         noauto = self.GetBoolFromConfig('autopeak', 'noauto')
 112         #noflatten: if true we do not flatten the curve
 113         noflatten = self.GetBoolFromConfig('autopeak', 'noflatten')
 114         peak_show = self.GetBoolFromConfig('autopeak', 'peak_show')
 115         persistence_length = self.GetFloatFromConfig('autopeak', 'persistence_length')
 116         #rebase: redefine the baseline
 117         rebase = self.GetBoolFromConfig('autopeak', 'rebase')
 118         reclick = self.GetBoolFromConfig('autopeak', 'reclick')
 119         slope_span = self.GetIntFromConfig('autopeak', 'auto_slope_span')
 120         T = self.GetFloatFromConfig('autopeak', 'temperature')
 121         usepl = self.GetBoolFromConfig('autopeak', 'usepl')
 122         if not usepl:
 123             pl_value = None
 124         else:
 125             pl_value = persistence_length / 10**9
 126         usepoints = self.GetBoolFromConfig('autopeak', 'usepoints')
 127         whatset_str = self.GetStringFromConfig('autopeak', 'whatset')
 128         if whatset_str == 'extension':
 129             whatset = lh.EXTENSION
 130         if whatset_str == 'retraction':
 131             whatset = lh.RETRACTION
 132
 133         #TODO: should this be variable?
 134         delta_force = 10
 135
 136         #setup header column labels for results
 137         if fit_function == 'wlc':
 138             fit_results = lib.results.ResultsWLC()
 139             segment_str = 'Persistence length'
 140             sigma_segment_str = 'sigma persistence length'
 141         elif fit_function == 'fjc' or fit_function == 'fjcPEG':
 142             fit_results = lib.results.ResultsFJC()
 143             segment_str = 'Kuhn length'
 144             sigma_segment_str = 'sigma Kuhn length'
 145         else:
 146             self.AppendToOutput('Unknown fit function, Please set fit_function as wlc, fjc or fjcPEG')
 147             return
 148
 149         #initialize output data vectors
 150         c_lengths = []
 151         p_lengths = []
 152         sigma_c_lengths = []
 153         sigma_p_lengths = []
 154         forces = []
 155         slopes = []
 156
 157         #pick up plot
 158         if plot is None:
 159             plot = copy.deepcopy(self.GetActivePlot())
 160         filename = self.GetActiveFile().name
 161
 162         #apply all active plotmanipulators and add the 'manipulated' data
 163         for plotmanipulator in self.plotmanipulators:
 164             if self.GetBoolFromConfig('core', 'plotmanipulators', plotmanipulator.name):
 165                 if plotmanipulator.name == 'flatten':
 166                     if not noflatten:
 167                         plot = plotmanipulator.method(plot, self.GetActiveFile())
 168                 else:
 169                     plot = plotmanipulator.method(plot, self.GetActiveFile())
 170
 171         #--Using points instead of nm interval
 172         if not usepoints:
 173             fit_points = None
 174
 175         #--Contact point arguments
 176         if reclick:
 177             contact_point, contact_point_index = self.pickup_contact_point(filename=filename)
 178         elif noauto:
 179             if self.wlccontact_index is None or self.wlccurrent != filename:
 180                 contact_point, contact_point_index = self.pickup_contact_point(filename=filename)
 181             else:
 182                 contact_point = self.wlccontact_point
 183                 contact_point_index = self.wlccontact_index
 184         else:
 185             #Automatically find contact point
 186             cindex = self.find_contact_point(plot)
 187             contact_point = self._clickize(plot.curves[lh.RETRACTION].x, plot.curves[lh.EXTENSION].y, cindex)
 188
 189         #peak_size comes from convolution curve
 190         peak_location, peak_size = self.find_current_peaks(plot=plot, noflatten=noflatten)
 191
 192         if len(peak_location) == 0:
 193             self.AppendToOutput('No peaks to fit.')
 194             return
 195
 196         #Pick up force baseline
 197         if baseline_clicks == 'contact point':
 198             try:
 199                 avg = plot.curves[lh.RETRACTION].y[contact_point_index]
 200             except:
 201                 avg = plot.curves[lh.RETRACTION].y[cindex]
 202
 203         if rebase or (self.basecurrent != filename) or self.basepoints is None:
 204             if baseline_clicks == 'automatic':
 205                 self.basepoints = []
 206                 base_index_0 = peak_location[-1] + self.fit_interval_nm(peak_location[-1], plot.curves[lh.RETRACTION].x, auto_right_baseline, False)
 207                 self.basepoints.append(self._clickize(plot.curves[lh.RETRACTION].x, plot.curves[lh.RETRACTION].y, base_index_0))
 208                 base_index_1 = self.basepoints[0].index + self.fit_interval_nm(self.basepoints[0].index, plot.curves[lh.RETRACTION].x, auto_left_baseline, False)
 209                 self.basepoints.append(self._clickize(plot.curves[lh.RETRACTION].x, plot.curves[lh.RETRACTION].y, base_index_1))
 210             if baseline_clicks == '1 point':
 211                 self.basepoints=self._measure_N_points(N=1, message='Click on 1 point to select the baseline.', whatset=whatset)
 212                 base_index_1 = self.basepoints[0].index + self.fit_interval_nm(self.basepoints[0].index, plot.curves[lh.RETRACTION].x, auto_left_baseline, False)
 213                 self.basepoints.append(self._clickize(plot.curves[lh.RETRACTION].x, plot.curves[lh.RETRACTION].y, base_index_1))
 214             if baseline_clicks == '2 points':
 215                 self.basepoints=self._measure_N_points(N=2, message='Click on 2 points to select the baseline.', whatset=whatset)
 216         if baseline_clicks != 'contact point':
 217             boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 218             boundaries.sort()
 219             to_average = plot.curves[lh.RETRACTION].y[boundaries[0]:boundaries[1]] #y points to average
 220             avg = mean(to_average)
 221             self.basecurrent = filename
 222
 223         x_values = plot.curves[lh.RETRACTION].x
 224         y_values = plot.curves[lh.RETRACTION].y
 225         for index, peak in enumerate(peak_location):
 226             #WLC FITTING
 227             #define fit interval
 228             if not usepoints:
 229                 fit_points = self.fit_interval_nm(peak, plot.curves[lh.RETRACTION].x, auto_fit_nm, True)
 230             peak_point = self._clickize(x_values, y_values, peak)
 231             other_fit_point=self._clickize(x_values, y_values, peak - fit_points)
 232
 233             #points for the fit
 234             points = [contact_point, peak_point, other_fit_point]
 235
 236             if abs(peak_point.index - other_fit_point.index) < 2:
 237                 continue
 238
 239             if fit_function == 'wlc':
 240                 params, yfit, xfit, fit_errors = self.wlc_fit(points, x_values, y_values, pl_value, T, return_errors=True)
 241             elif fit_function == 'fjc':
 242                 params, yfit, xfit, fit_errors = self.fjc_fit(points, x_values, y_values, pl_value, T, return_errors=True)
 243             elif fit_function == 'fjcPEG':
 244                 params, yfit, xfit, fit_errors = self.fjcPEG_fit(points, x_values, y_values, pl_value, T, return_errors=True)
 245
 246             #Measure forces
 247             delta_to_measure = y_values[peak - delta_force:peak + delta_force]
 248             y = min(delta_to_measure)
 249             #save force values (pN)
 250             #Measure slopes
 251             slope = self.linefit_between(peak - slope_span, peak, whatset=lh.RETRACTION)[0]
 252
 253             #check fitted data and, if right, add peak to the measurement
 254             fit_result = lib.results.Result()
 255
 256             fit_result.result['Contour length'] = params[0]
 257             fit_result.result['sigma contour length'] = fit_errors[0]
 258             fit_result.result['Rupture force'] = abs(y - avg)
 259             fit_result.result['Slope'] = slope
 260             active_file = self.GetActiveFile()
 261             if active_file.driver.retract_velocity:
 262                 fit_result.result['Loading rate'] = slope * active_file.driver.retract_velocity
 263             else:
 264                 fit_result.result['Loading rate'] = -1
 265             if len(params) == 1: #if we did choose 1-value fit
 266                 fit_result.result[segment_str] = pl_value
 267                 fit_result.result[sigma_segment_str] = 0
 268
 269                 p_lengths.append(pl_value)
 270                 c_lengths.append(params[0]*(1.0e+9))
 271                 sigma_p_lengths.append(0)
 272                 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 273                 forces.append(abs(y-avg)*(1.0e+12))
 274                 slopes.append(slope)
 275             else: #2-value fit
 276                 p_leng = params[1] * (1.0e+9)
 277                 #check if persistence length makes sense, otherwise discard peak.
 278                 if p_leng > auto_min_p and p_leng < auto_max_p:
 279                     fit_result.result[segment_str] = params[1]
 280                     fit_result.result[sigma_segment_str] = fit_errors[1]
 281
 282                     p_lengths.append(p_leng)
 283                     c_lengths.append(params[0]*(1.0e+9))
 284                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 285                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 286                     forces.append(abs(y-avg)*(1.0e+12))
 287                     slopes.append(slope)
 288                 else:
 289                     fit_result.result = {}
 290
 291             if len(fit_result.result) > 0:
 292                 fit_result.label = fit_function + '_' + str(index)
 293                 fit_result.title = plot.curves[lh.RETRACTION].title
 294                 fit_result.units.x = plot.curves[lh.RETRACTION].units.x
 295                 fit_result.units.y = plot.curves[lh.RETRACTION].units.y
 296                 fit_result.visible = True
 297                 fit_result.x = xfit
 298                 fit_result.y = yfit
 299                 fit_results.results.append(fit_result)
 300
 301         if fit_results.results:
 302             fit_results.set_multipliers(0)
 303             plot = self.GetActivePlot()
 304             plot.results[fit_function] = fit_results
 305             if peak_show:
 306                 plugin = lib.plugin.Plugin()
 307                 plugin.name = 'autopeak'
 308                 plugin.section = 'autopeak'
 309                 plugin.prefix = 'peak_'
 310                 self.do_peaks(plugin=plugin, peak_location=peak_location, peak_size=peak_size)
 311             else:
 312                 self.UpdatePlot()
 313
 314         else:
 315             self.AppendToOutput('No peaks found.')
 316
 317         #TODO:
 318         #self.do_note('autopeak')
 319
 320     def find_current_peaks(self, plot=None, noflatten=True):
 321         if not noflatten:
 322             plot_temp = self.plotmanip_flatten(plot, self.GetActiveFile(), customvalue=1)
 323         peak_location, peak_size = self.has_peaks(plot_temp)
 324         return peak_location, peak_size
 325
 326     def fit_interval_nm(self, start_index, x_vect, nm, backwards):
 327         '''
 328         Calculates the number of points to fit, given a fit interval in nm
 329         start_index: index of point
 330         plot: plot to use
 331         backwards: if true, finds a point backwards.
 332         '''
 333
 334         c = 0
 335         i = start_index
 336         maxlen=len(x_vect)
 337         while abs(x_vect[i] - x_vect[start_index]) * (10**9) < nm:
 338             if i == 0 or i == maxlen-1: #we reached boundaries of vector!
 339                 return c
 340             if backwards:
 341                 i -= 1
 342             else:
 343                 i += 1
 344             c += 1
 345         return c
 346
 347     def pickup_contact_point(self, filename=''):
 348         '''
 349         macro to pick up the contact point by clicking
 350         '''
 351         contact_point = self._measure_N_points(N=1, message='Please click on the contact point.')[0]
 352         contact_point_index = contact_point.index
 353         self.wlccontact_point = contact_point
 354         self.wlccontact_index = contact_point.index
 355         self.wlccurrent = filename
 356         return contact_point, contact_point_index
 357
 358