generalvclamp.py

   1 #!/usr/bin/env python
   2
   3 '''
   4 generalvclamp.py
   5
   6 Plugin regarding general velocity clamp measurements
   7 '''
   8
   9 from libhooke import WX_GOOD, ClickedPoint
  10 import wxversion
  11 wxversion.select(WX_GOOD)
  12 from wx import PostEvent
  13 import numpy as np
  14 import scipy as sp
  15 import copy
  16 import os.path
  17 import time
  18
  19 import warnings
  20 warnings.simplefilter('ignore',np.RankWarning)
  21
  22
  23 class generalvclampCommands:
  24
  25     def _plug_init(self):
  26         self.basecurrent=None
  27         self.basepoints=None
  28         self.autofile=''
  29
  30     def do_distance(self,args):
  31         '''
  32         DISTANCE
  33         (generalvclamp.py)
  34         Measure the distance (in nm) between two points.
  35         For a standard experiment this is the delta X distance.
  36         For a force clamp experiment this is the delta Y distance (actually becomes
  37         an alias of zpiezo)
  38         -----------------
  39         Syntax: distance
  40         '''
  41         if self.current.curve.experiment == 'clamp':
  42             print 'You wanted to use zpiezo perhaps?'
  43             return
  44         else:
  45             dx,unitx,dy,unity=self._delta(set=1)
  46             print str(dx*(10**9))+' nm'
  47
  48
  49     def do_force(self,args):
  50         '''
  51         FORCE
  52         (generalvclamp.py)
  53         Measure the force difference (in pN) between two points
  54         ---------------
  55         Syntax: force
  56         '''
  57         if self.current.curve.experiment == 'clamp':
  58             print 'This command makes no sense for a force clamp experiment.'
  59             return
  60         dx,unitx,dy,unity=self._delta(set=1)
  61         print str(dy*(10**12))+' pN'
  62
  63
  64     def do_forcebase(self,args):
  65         '''
  66         FORCEBASE
  67         (generalvclamp.py)
  68         Measures the difference in force (in pN) between a point and a baseline
  69         took as the average between two points.
  70
  71         The baseline is fixed once for a given curve and different force measurements,
  72         unless the user wants it to be recalculated
  73         ------------
  74         Syntax: forcebase [rebase]
  75                 rebase: Forces forcebase to ask again the baseline
  76                 max: Instead of asking for a point to measure, asks for two points and use
  77                      the maximum peak in between
  78         '''
  79         rebase=False #if true=we select rebase
  80         maxpoint=False #if true=we measure the maximum peak
  81
  82         plot=self._get_displayed_plot()
  83         whatset=1 #fixme: for all sets
  84         if 'rebase' in args or (self.basecurrent != self.current.path):
  85             rebase=True
  86         if 'max' in args:
  87             maxpoint=True
  88
  89         if rebase:
  90             print 'Select baseline'
  91             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
  92             self.basecurrent=self.current.path
  93
  94         if maxpoint:
  95             print 'Select two points'
  96             points=self._measure_N_points(N=2, whatset=whatset)
  97             boundpoints=[points[0].index, points[1].index]
  98             boundpoints.sort()
  99             try:
 100                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
 101             except ValueError:
 102                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
 103         else:
 104             print 'Select point to measure'
 105             points=self._measure_N_points(N=1, whatset=whatset)
 106             #whatplot=points[0].dest
 107             y=points[0].graph_coords[1]
 108
 109         #fixme: code duplication
 110         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 111         boundaries.sort()
 112         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
 113
 114         avg=np.mean(to_average)
 115         forcebase=abs(y-avg)
 116         print str(forcebase*(10**12))+' pN'
 117
 118
 119     def plotmanip_flatten(self, plot, current, customvalue=False):
 120         '''
 121         Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
 122         the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
 123         given by the configuration file or by the customvalue.
 124
 125         customvalue= int (>0) --> starts the function even if config says no (default=False)
 126         '''
 127
 128         #not a smfs curve...
 129         if current.curve.experiment != 'smfs':
 130             return plot
 131
 132         #only one set is present...
 133         if len(self.plots[0].vectors) != 2:
 134             return plot
 135
 136         #config is not flatten, and customvalue flag is false too
 137         if (not self.config['flatten']) and (not customvalue):
 138             return plot
 139
 140         max_exponent=12
 141         delta_contact=0
 142
 143         if customvalue:
 144             max_cycles=customvalue
 145         else:
 146             max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
 147
 148         contact_index=self.find_contact_point()
 149
 150         valn=[[] for item in range(max_exponent)]
 151         yrn=[0.0 for item in range(max_exponent)]
 152         errn=[0.0 for item in range(max_exponent)]
 153
 154         for i in range(int(max_cycles)):
 155
 156             x_ext=plot.vectors[0][0][contact_index+delta_contact:]
 157             y_ext=plot.vectors[0][1][contact_index+delta_contact:]
 158             x_ret=plot.vectors[1][0][contact_index+delta_contact:]
 159             y_ret=plot.vectors[1][1][contact_index+delta_contact:]
 160             for exponent in range(max_exponent):
 161                 try:
 162                     valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
 163                     yrn[exponent]=sp.polyval(valn[exponent],x_ret)
 164                     errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
 165                 except TypeError:
 166                     print 'Cannot flatten!'
 167                     return plot
 168
 169             best_exponent=errn.index(min(errn))
 170
 171             #extension
 172             ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
 173             yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
 174             #retraction
 175             ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
 176             yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
 177
 178             ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
 179             ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
 180
 181             plot.vectors[0][1]=list(ycorr_ext)
 182             plot.vectors[1][1]=list(ycorr_ret)
 183
 184         return plot
 185
 186     #---SLOPE---
 187     def do_slope(self,args):
 188         '''
 189         SLOPE
 190         (generalvclamp.py)
 191         Measures the slope of a delimited chunk on the return trace.
 192         The chunk can be delimited either by two manual clicks, or have
 193         a fixed width, given as an argument.
 194         ---------------
 195         Syntax: slope [width]
 196                 The facultative [width] parameter specifies how many
 197                 points will be considered for the fit. If [width] is
 198                 specified, only one click will be required.
 199         (c) Marco Brucale, Massimo Sandal 2008
 200         '''
 201
 202         # Reads the facultative width argument
 203         try:
 204             fitspan=int(args)
 205         except:
 206             fitspan=0
 207
 208         # Decides between the two forms of user input, as per (args)
 209         if fitspan == 0:
 210             # Gets the Xs of two clicked points as indexes on the current curve vector
 211             print 'Click twice to delimit chunk'
 212             clickedpoints=[]
 213             points=self._measure_N_points(N=2,whatset=1)
 214             clickedpoints=[points[0].index,points[1].index]
 215             clickedpoints.sort()
 216         else:
 217             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 218             clickedpoints=[]
 219             points=self._measure_N_points(N=1,whatset=1)
 220             clickedpoints=[points[0].index-fitspan,points[0].index]
 221
 222         # Calls the function linefit_between
 223         parameters=[0,0,[],[]]
 224         parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 225
 226         # Outputs the relevant slope parameter
 227         print 'Slope:'
 228         print str(parameters[0])
 229
 230         # Makes a vector with the fitted parameters and sends it to the GUI
 231         xtoplot=parameters[2]
 232         ytoplot=[]
 233         x=0
 234         for x in xtoplot:
 235             ytoplot.append((x*parameters[0])+parameters[1])
 236
 237         clickvector_x, clickvector_y=[], []
 238         for item in points:
 239             clickvector_x.append(item.graph_coords[0])
 240             clickvector_y.append(item.graph_coords[1])
 241
 242         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 243
 244         lineplot.add_set(xtoplot,ytoplot)
 245         lineplot.add_set(clickvector_x, clickvector_y)
 246
 247         if lineplot.styles==[]:
 248             lineplot.styles=[None,None,None,'scatter']
 249         else:
 250             lineplot.styles+=[None,'scatter']
 251
 252         self._send_plot([lineplot])
 253
 254     def linefit_between(self,index1,index2):
 255         '''
 256         Creates two vectors (xtofit,ytofit) slicing out from the
 257         current return trace a portion delimited by the two indexes
 258         given as arguments.
 259         Then does a least squares linear fit on that slice.
 260         Finally returns [0]=the slope, [1]=the intercept of the
 261         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 262         used for the fit.
 263         (c) Marco Brucale, Massimo Sandal 2008
 264         '''
 265         # Translates the indexes into two vectors containing the x,y data to fit
 266         xtofit=self.plots[0].vectors[1][0][index1:index2]
 267         ytofit=self.plots[0].vectors[1][1][index1:index2]
 268
 269         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 270         linefit=[]
 271         linefit=np.polyfit(xtofit,ytofit,1)
 272
 273         return (linefit[0],linefit[1],xtofit,ytofit)
 274
 275 #====================
 276 #AUTOMATIC ANALYSES
 277 #====================
 278     def do_autopeak(self,args):
 279         '''
 280         AUTOPEAK
 281         (generalvclamp.py)
 282         Automatically performs a number of analyses on the peaks of the given curve.
 283         Currently it automatically:
 284         - fits peaks with WLC function
 285         - measures peak maximum forces with a baseline
 286         - measures slope in proximity of peak maximum
 287         Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
 288
 289         Syntax:
 290         autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
 291
 292         rebase : Re-asks baseline interval
 293
 294         pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
 295                      the fit will be a 2-variable
 296                      fit. DO NOT put spaces between 'pl', '=' and the value.
 297                      The value must be in meters.
 298                      Scientific notation like 0.35e-9 is fine.
 299
 300         t=[value] : Use a user-defined temperature. The value must be in
 301                     kelvins; by default it is 293 K.
 302                     DO NOT put spaces between 't', '=' and the value.
 303
 304         noauto : allows for clicking the contact point by
 305                  hand (otherwise it is automatically estimated) the first time.
 306                  If subsequent measurements are made, the same contact point
 307                  clicked the first time is used
 308
 309         reclick : redefines by hand the contact point, if noauto has been used before
 310                   but the user is unsatisfied of the previously choosen contact point.
 311
 312         When you first issue the command, it will ask for the filename. If you are giving the filename
 313         of an existing file, autopeak will resume it and append measurements to it. If you are giving
 314         a new filename, it will create the file and append to it until you close Hooke.
 315
 316         Useful variables (to set with SET command):
 317
 318         temperature= temperature of the system for wlc fit (in K)
 319         auto_fit_points = number of points to fit before the peak maximum, for wlc
 320         auto_slope_span = number of points on which measure the slope, for slope
 321         '''
 322
 323         pl_value=None
 324         T=self.config['temperature']
 325
 326         fit_points=int(self.config['auto_fit_points'])
 327         slope_span=int(self.config['auto_slope_span'])
 328
 329         delta_force=10
 330         rebase=False #if true=we select rebase
 331
 332         #Pick up plot
 333         displayed_plot=self._get_displayed_plot(0)
 334
 335         if self.current.curve.experiment != 'smfs' or len(displayed_plot.vectors) < 2:
 336             print 'Cannot work on this curve.'
 337             return
 338
 339         #Look for custom persistent length / custom temperature
 340         for arg in args.split():
 341             #look for a persistent length argument.
 342             if 'pl=' in arg:
 343                 pl_expression=arg.split('=')
 344                 pl_value=float(pl_expression[1]) #actual value
 345             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 346             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 347                 t_expression=arg.split('=')
 348                 T=float(t_expression[1])
 349
 350         #Handle contact point arguments
 351         #FIXME: code duplication
 352         if 'reclick' in args.split():
 353             print 'Click contact point'
 354             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 355             contact_point_index=contact_point.index
 356             self.wlccontact_point=contact_point
 357             self.wlccontact_index=contact_point.index
 358             self.wlccurrent=self.current.path
 359         elif 'noauto' in args.split():
 360             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 361                 print 'Click contact point'
 362                 contact_point=self._measure_N_points(N=1, whatset=1)[0]
 363                 contact_point_index=contact_point.index
 364                 self.wlccontact_point=contact_point
 365                 self.wlccontact_index=contact_point.index
 366                 self.wlccurrent=self.current.path
 367             else:
 368                 contact_point=self.wlccontact_point
 369                 contact_point_index=self.wlccontact_index
 370         else:
 371             #Automatically find contact point
 372             cindex=self.find_contact_point()
 373             contact_point=ClickedPoint()
 374             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 375             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 376             contact_point.is_marker=True
 377
 378
 379         #Pick up force baseline
 380         whatset=1 #fixme: for all sets
 381         if 'rebase' in args or (self.basecurrent != self.current.path):
 382             rebase=True
 383         if rebase:
 384             print 'Select baseline'
 385             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 386             self.basecurrent=self.current.path
 387         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 388         boundaries.sort()
 389         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 390         avg=np.mean(to_average)
 391
 392         #Find peaks.
 393         defplot=self.current.curve.default_plots()[0]
 394         flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 395         defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 396         peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 397
 398         #Create a new plot to send
 399         fitplot=copy.deepcopy(displayed_plot)
 400
 401         #Initialize data vectors
 402         c_lengths=[]
 403         p_lengths=[]
 404         forces=[]
 405         slopes=[]
 406
 407         #Cycle between peaks and do analysis.
 408         for peak in peak_location:
 409             #Do WLC fits.
 410             #FIXME: clean wlc fitting, to avoid this clickedpoint hell
 411             #-create a clicked point for the peak point
 412             peak_point=ClickedPoint()
 413             peak_point.absolute_coords=displayed_plot.vectors[1][0][peak], displayed_plot.vectors[1][1][peak]
 414             peak_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 415             #-create a clicked point for the other fit point
 416             other_fit_point=ClickedPoint()
 417             other_fit_point.absolute_coords=displayed_plot.vectors[1][0][peak-fit_points], displayed_plot.vectors[1][1][peak-fit_points]
 418             other_fit_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 419             #do the fit
 420             points=[contact_point, peak_point, other_fit_point]
 421             params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
 422             #save wlc values (nm)
 423             c_lengths.append(params[0]*(1.0e+9))
 424             if len(params)==2: #if we did choose 2-value fit
 425                 p_lengths.append(params[1]*(1.0e+9))
 426             else:
 427                 p_lengths.append(pl_value)
 428             #Add WLC fit lines to plot
 429             fitplot.add_set(xfit,yfit)
 430             if len(fitplot.styles)==0:
 431                 fitplot.styles=[]
 432             else:
 433                 fitplot.styles.append(None)
 434
 435             #Measure forces
 436             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 437             y=min(delta_to_measure)
 438             #save force values (pN)
 439             forces.append(abs(y-avg)*(1.0e+12))
 440
 441             #Measure slopes
 442             slope=self.linefit_between(peak-slope_span,peak)[0]
 443             slopes.append(slope)
 444
 445         #Show wlc fits and peak locations
 446         self._send_plot([fitplot])
 447         self.do_peaks('')
 448
 449         #Ask the user what peaks to ignore from analysis.
 450         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 451         print 'N to discard measurement'
 452         exclude_raw=raw_input('Input:')
 453         if exclude_raw=='N':
 454             print 'Discarded.'
 455             return
 456         if not exclude_raw=='':
 457             exclude=exclude_raw.split(',')
 458             try:
 459                 exclude=[int(item) for item in exclude]
 460                 for i in exclude:
 461                     c_lengths[i]=None
 462                     p_lengths[i]=None
 463                     forces[i]=None
 464                     slopes[i]=None
 465             except:
 466                  print 'Bad input, taking all...'
 467         #Clean data vectors from ignored peaks
 468         c_lengths=[item for item in c_lengths if item != None]
 469         p_lengths=[item for item in p_lengths if item != None]
 470         forces=[item for item in forces if item != None]
 471         slopes=[item for item in slopes if item != None]
 472         print 'contour (nm)',c_lengths
 473         print 'p (nm)',p_lengths
 474         print 'forces (pN)',forces
 475         print 'slopes (N/m)',slopes
 476
 477         #Save file info
 478         if self.autofile=='':
 479             self.autofile=raw_input('Filename? (return to ignore) ')
 480             if self.autofile=='':
 481                 print 'Not saved.'
 482                 return
 483
 484         if not os.path.exists(self.autofile):
 485             f=open(self.autofile,'w+')
 486             f.write('Analysis started '+time.asctime()+'\n')
 487             f.write('----------------------------------------\n')
 488             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) \n')
 489             f.close()
 490
 491         print 'Saving...'
 492         f=open(self.autofile,'a+')
 493
 494         f.write(self.current.path+'\n')
 495         for i in range(len(c_lengths)):
 496             f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
 497
 498         f.close()
 499         self.do_note('autopeak')
 500