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 TypeError:
 172                     print 'Cannot flatten!'
 173                     return plot
 174
 175             best_exponent=errn.index(min(errn))
 176
 177             #extension
 178             ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
 179             yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
 180             #retraction
 181             ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
 182             yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
 183
 184             ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
 185             ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
 186
 187             plot.vectors[0][1]=list(ycorr_ext)
 188             plot.vectors[1][1]=list(ycorr_ret)
 189
 190         return plot
 191
 192     #---SLOPE---
 193     def do_slope(self,args):
 194         '''
 195         SLOPE
 196         (generalvclamp.py)
 197         Measures the slope of a delimited chunk on the return trace.
 198         The chunk can be delimited either by two manual clicks, or have
 199         a fixed width, given as an argument.
 200         ---------------
 201         Syntax: slope [width]
 202                 The facultative [width] parameter specifies how many
 203                 points will be considered for the fit. If [width] is
 204                 specified, only one click will be required.
 205         (c) Marco Brucale, Massimo Sandal 2008
 206         '''
 207
 208         # Reads the facultative width argument
 209         try:
 210             fitspan=int(args)
 211         except:
 212             fitspan=0
 213
 214         # Decides between the two forms of user input, as per (args)
 215         if fitspan == 0:
 216             # Gets the Xs of two clicked points as indexes on the current curve vector
 217             print 'Click twice to delimit chunk'
 218             clickedpoints=[]
 219             points=self._measure_N_points(N=2,whatset=1)
 220             clickedpoints=[points[0].index,points[1].index]
 221             clickedpoints.sort()
 222         else:
 223             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 224             clickedpoints=[]
 225             points=self._measure_N_points(N=1,whatset=1)
 226             clickedpoints=[points[0].index-fitspan,points[0].index]
 227
 228         # Calls the function linefit_between
 229         parameters=[0,0,[],[]]
 230         parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 231
 232         # Outputs the relevant slope parameter
 233         print 'Slope:'
 234         print str(parameters[0])
 235         to_dump='slope '+self.current.path+' '+str(parameters[0])
 236         self.outlet.push(to_dump)
 237
 238         # Makes a vector with the fitted parameters and sends it to the GUI
 239         xtoplot=parameters[2]
 240         ytoplot=[]
 241         x=0
 242         for x in xtoplot:
 243             ytoplot.append((x*parameters[0])+parameters[1])
 244
 245         clickvector_x, clickvector_y=[], []
 246         for item in points:
 247             clickvector_x.append(item.graph_coords[0])
 248             clickvector_y.append(item.graph_coords[1])
 249
 250         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 251
 252         lineplot.add_set(xtoplot,ytoplot)
 253         lineplot.add_set(clickvector_x, clickvector_y)
 254
 255         if lineplot.styles==[]:
 256             lineplot.styles=[None,None,None,'scatter']
 257         else:
 258             lineplot.styles+=[None,'scatter']
 259
 260         self._send_plot([lineplot])
 261
 262     def linefit_between(self,index1,index2):
 263         '''
 264         Creates two vectors (xtofit,ytofit) slicing out from the
 265         current return trace a portion delimited by the two indexes
 266         given as arguments.
 267         Then does a least squares linear fit on that slice.
 268         Finally returns [0]=the slope, [1]=the intercept of the
 269         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 270         used for the fit.
 271         (c) Marco Brucale, Massimo Sandal 2008
 272         '''
 273         # Translates the indexes into two vectors containing the x,y data to fit
 274         xtofit=self.plots[0].vectors[1][0][index1:index2]
 275         ytofit=self.plots[0].vectors[1][1][index1:index2]
 276
 277         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 278         linefit=[]
 279         linefit=np.polyfit(xtofit,ytofit,1)
 280
 281         return (linefit[0],linefit[1],xtofit,ytofit)
 282
 283 #====================
 284 #AUTOMATIC ANALYSES
 285 #====================
 286     def do_autopeak(self,args):
 287         '''
 288         AUTOPEAK
 289         (generalvclamp.py)
 290         Automatically performs a number of analyses on the peaks of the given curve.
 291         Currently it automatically:
 292         - fits peaks with WLC function
 293         - measures peak maximum forces with a baseline
 294         - measures slope in proximity of peak maximum
 295         Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
 296
 297         Syntax:
 298         autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
 299
 300         rebase : Re-asks baseline interval
 301
 302         pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
 303                      the fit will be a 2-variable
 304                      fit. DO NOT put spaces between 'pl', '=' and the value.
 305                      The value must be in meters.
 306                      Scientific notation like 0.35e-9 is fine.
 307
 308         t=[value] : Use a user-defined temperature. The value must be in
 309                     kelvins; by default it is 293 K.
 310                     DO NOT put spaces between 't', '=' and the value.
 311
 312         noauto : allows for clicking the contact point by
 313                  hand (otherwise it is automatically estimated) the first time.
 314                  If subsequent measurements are made, the same contact point
 315                  clicked the first time is used
 316
 317         reclick : redefines by hand the contact point, if noauto has been used before
 318                   but the user is unsatisfied of the previously choosen contact point.
 319
 320         When you first issue the command, it will ask for the filename. If you are giving the filename
 321         of an existing file, autopeak will resume it and append measurements to it. If you are giving
 322         a new filename, it will create the file and append to it until you close Hooke.
 323
 324         Useful variables (to set with SET command):
 325
 326         temperature= temperature of the system for wlc fit (in K)
 327         auto_fit_points = number of points to fit before the peak maximum, for wlc
 328         auto_slope_span = number of points on which measure the slope, for slope
 329         '''
 330
 331         pl_value=None
 332         T=self.config['temperature']
 333
 334         fit_points=int(self.config['auto_fit_points'])
 335         slope_span=int(self.config['auto_slope_span'])
 336
 337         delta_force=10
 338         rebase=False #if true=we select rebase
 339
 340         #Pick up plot
 341         displayed_plot=self._get_displayed_plot(0)
 342
 343         if self.current.curve.experiment != 'smfs' or len(displayed_plot.vectors) < 2:
 344             print 'Cannot work on this curve.'
 345             return
 346
 347         #Look for custom persistent length / custom temperature
 348         for arg in args.split():
 349             #look for a persistent length argument.
 350             if 'pl=' in arg:
 351                 pl_expression=arg.split('=')
 352                 pl_value=float(pl_expression[1]) #actual value
 353             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 354             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 355                 t_expression=arg.split('=')
 356                 T=float(t_expression[1])
 357
 358         #Handle contact point arguments
 359         #FIXME: code duplication
 360         if 'reclick' in args.split():
 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         elif 'noauto' in args.split():
 368             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 369                 print 'Click contact point'
 370                 contact_point=self._measure_N_points(N=1, whatset=1)[0]
 371                 contact_point_index=contact_point.index
 372                 self.wlccontact_point=contact_point
 373                 self.wlccontact_index=contact_point.index
 374                 self.wlccurrent=self.current.path
 375             else:
 376                 contact_point=self.wlccontact_point
 377                 contact_point_index=self.wlccontact_index
 378         else:
 379             #Automatically find contact point
 380             cindex=self.find_contact_point()
 381             contact_point=ClickedPoint()
 382             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 383             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 384             contact_point.is_marker=True
 385
 386
 387         #Pick up force baseline
 388         whatset=1 #fixme: for all sets
 389         if 'rebase' in args or (self.basecurrent != self.current.path):
 390             rebase=True
 391         if rebase:
 392             print 'Select baseline'
 393             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 394             self.basecurrent=self.current.path
 395         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 396         boundaries.sort()
 397         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 398         avg=np.mean(to_average)
 399
 400         #Find peaks.
 401         defplot=self.current.curve.default_plots()[0]
 402         flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 403         defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 404         peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 405
 406         #Create a new plot to send
 407         fitplot=copy.deepcopy(displayed_plot)
 408
 409         #Initialize data vectors
 410         c_lengths=[]
 411         p_lengths=[]
 412         forces=[]
 413         slopes=[]
 414
 415         #Cycle between peaks and do analysis.
 416         for peak in peak_location:
 417             #Do WLC fits.
 418             #FIXME: clean wlc fitting, to avoid this clickedpoint hell
 419             #-create a clicked point for the peak point
 420             peak_point=ClickedPoint()
 421             peak_point.absolute_coords=displayed_plot.vectors[1][0][peak], displayed_plot.vectors[1][1][peak]
 422             peak_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 423             #-create a clicked point for the other fit point
 424             other_fit_point=ClickedPoint()
 425             other_fit_point.absolute_coords=displayed_plot.vectors[1][0][peak-fit_points], displayed_plot.vectors[1][1][peak-fit_points]
 426             other_fit_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 427             #do the fit
 428             points=[contact_point, peak_point, other_fit_point]
 429             params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
 430             #save wlc values (nm)
 431             c_lengths.append(params[0]*(1.0e+9))
 432             if len(params)==2: #if we did choose 2-value fit
 433                 p_lengths.append(params[1]*(1.0e+9))
 434             else:
 435                 p_lengths.append(pl_value)
 436             #Add WLC fit lines to plot
 437             fitplot.add_set(xfit,yfit)
 438             if len(fitplot.styles)==0:
 439                 fitplot.styles=[]
 440             else:
 441                 fitplot.styles.append(None)
 442
 443             #Measure forces
 444             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 445             y=min(delta_to_measure)
 446             #save force values (pN)
 447             forces.append(abs(y-avg)*(1.0e+12))
 448
 449             #Measure slopes
 450             slope=self.linefit_between(peak-slope_span,peak)[0]
 451             slopes.append(slope)
 452
 453         #Show wlc fits and peak locations
 454         self._send_plot([fitplot])
 455         self.do_peaks('')
 456
 457         #Ask the user what peaks to ignore from analysis.
 458         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 459         print 'N to discard measurement'
 460         exclude_raw=raw_input('Input:')
 461         if exclude_raw=='N':
 462             print 'Discarded.'
 463             return
 464         if not exclude_raw=='':
 465             exclude=exclude_raw.split(',')
 466             try:
 467                 exclude=[int(item) for item in exclude]
 468                 for i in exclude:
 469                     c_lengths[i]=None
 470                     p_lengths[i]=None
 471                     forces[i]=None
 472                     slopes[i]=None
 473             except:
 474                  print 'Bad input, taking all...'
 475         #Clean data vectors from ignored peaks
 476         c_lengths=[item for item in c_lengths if item != None]
 477         p_lengths=[item for item in p_lengths if item != None]
 478         forces=[item for item in forces if item != None]
 479         slopes=[item for item in slopes if item != None]
 480         print 'contour (nm)',c_lengths
 481         print 'p (nm)',p_lengths
 482         print 'forces (pN)',forces
 483         print 'slopes (N/m)',slopes
 484
 485         #Save file info
 486         if self.autofile=='':
 487             self.autofile=raw_input('Filename? (return to ignore) ')
 488             if self.autofile=='':
 489                 print 'Not saved.'
 490                 return
 491
 492         if not os.path.exists(self.autofile):
 493             f=open(self.autofile,'w+')
 494             f.write('Analysis started '+time.asctime()+'\n')
 495             f.write('----------------------------------------\n')
 496             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) \n')
 497             f.close()
 498
 499         print 'Saving...'
 500         f=open(self.autofile,'a+')
 501
 502         f.write(self.current.path+'\n')
 503         for i in range(len(c_lengths)):
 504             f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
 505
 506         f.close()
 507         self.do_note('autopeak')
 508