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