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         try:
 232             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 233         except:
 234             print 'Cannot fit. Did you click twice the same point?'
 235             return
 236
 237         # Outputs the relevant slope parameter
 238         print 'Slope:'
 239         print str(parameters[0])
 240         to_dump='slope '+self.current.path+' '+str(parameters[0])
 241         self.outlet.push(to_dump)
 242
 243         # Makes a vector with the fitted parameters and sends it to the GUI
 244         xtoplot=parameters[2]
 245         ytoplot=[]
 246         x=0
 247         for x in xtoplot:
 248             ytoplot.append((x*parameters[0])+parameters[1])
 249
 250         clickvector_x, clickvector_y=[], []
 251         for item in points:
 252             clickvector_x.append(item.graph_coords[0])
 253             clickvector_y.append(item.graph_coords[1])
 254
 255         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 256
 257         lineplot.add_set(xtoplot,ytoplot)
 258         lineplot.add_set(clickvector_x, clickvector_y)
 259
 260
 261         if lineplot.styles==[]:
 262             lineplot.styles=[None,None,None,'scatter']
 263         else:
 264             lineplot.styles+=[None,'scatter']
 265         if lineplot.colors==[]:
 266             lineplot.styles=[None,None,None,None]
 267         else:
 268             lineplot.colors+=[None,None]
 269
 270
 271         self._send_plot([lineplot])
 272
 273     def linefit_between(self,index1,index2,whatset=1):
 274         '''
 275         Creates two vectors (xtofit,ytofit) slicing out from the
 276         current return trace a portion delimited by the two indexes
 277         given as arguments.
 278         Then does a least squares linear fit on that slice.
 279         Finally returns [0]=the slope, [1]=the intercept of the
 280         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 281         used for the fit.
 282         (c) Marco Brucale, Massimo Sandal 2008
 283         '''
 284         # Translates the indexes into two vectors containing the x,y data to fit
 285         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
 286         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
 287
 288         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 289         linefit=[]
 290         linefit=np.polyfit(xtofit,ytofit,1)
 291
 292         return (linefit[0],linefit[1],xtofit,ytofit)
 293
 294
 295
 296 #====================
 297 #AUTOMATIC ANALYSES
 298 #====================
 299     '''
 300     def do_autopeak(self,args):
 301         #FIXME: this function is too long. split it and make it rational.
 302         #FIXME: also, *generalize fits* to allow FJC and any other model in the future!
 303
 304         def fit_interval_nm(start_index,plot,nm,backwards):
 305     '''
 306
 307     '''
 308             Calculates the number of points to fit, given a fit interval in nm
 309             start_index: index of point
 310             plot: plot to use
 311             backwards: if true, finds a point backwards.
 312     '''
 313
 314     '''
 315             x_vect=plot.vectors[1][0]
 316
 317             c=0
 318             i=start_index
 319             start=x_vect[start_index]
 320             maxlen=len(x_vect)
 321             while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
 322                 if i==0 or i==maxlen-1: #we reached boundaries of vector!
 323                     return c
 324
 325                 if backwards:
 326                     i-=1
 327                 else:
 328                     i+=1
 329                 c+=1
 330             return c
 331
 332
 333         pl_value=None
 334         T=self.config['temperature']
 335
 336         if 'usepoints' in args.split():
 337             fit_points=int(self.config['auto_fit_points'])
 338             usepoints=True
 339         else:
 340             fit_points=None
 341             usepoints=False
 342
 343         slope_span=int(self.config['auto_slope_span'])
 344         delta_force=10
 345         rebase=False #if true=we select rebase
 346
 347         #Pick up plot
 348         displayed_plot=self._get_displayed_plot(0)
 349
 350         if self.current.curve.experiment != 'smfs' or len(displayed_plot.vectors) < 2:
 351             print 'Cannot work on this curve.'
 352             return
 353
 354         #Look for custom persistent length / custom temperature
 355         for arg in args.split():
 356             #look for a persistent length argument.
 357             if 'pl=' in arg:
 358                 pl_expression=arg.split('=')
 359                 pl_value=float(pl_expression[1]) #actual value
 360             #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
 361             if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
 362                 t_expression=arg.split('=')
 363                 T=float(t_expression[1])
 364
 365         #Handle contact point arguments
 366         def pickup_contact_point():
 367             #macro to pick up the contact point by clicking
 368             contact_point=self._measure_N_points(N=1, whatset=1)[0]
 369             contact_point_index=contact_point.index
 370             self.wlccontact_point=contact_point
 371             self.wlccontact_index=contact_point.index
 372             self.wlccurrent=self.current.path
 373             return contact_point, contact_point_index
 374
 375         if 'reclick' in args.split():
 376             print 'Click contact point'
 377             contact_point, contact_point_index = pickup_contact_point()
 378         elif 'noauto' in args.split():
 379             if self.wlccontact_index==None or self.wlccurrent != self.current.path:
 380                 print 'Click contact point'
 381                 contact_point , contact_point_index = pickup_contact_point()
 382             else:
 383                 contact_point=self.wlccontact_point
 384                 contact_point_index=self.wlccontact_index
 385         else:
 386             #Automatically find contact point
 387             cindex=self.find_contact_point()
 388             contact_point=ClickedPoint()
 389             contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
 390             contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 391             contact_point.is_marker=True
 392
 393
 394         #Find peaks.
 395         defplot=self.current.curve.default_plots()[0]
 396         flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 397         defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 398         peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
 399
 400         #Create a new plot to send
 401         fitplot=copy.deepcopy(displayed_plot)
 402
 403         #Pick up force baseline
 404         whatset=1 #fixme: for all sets
 405         if 'rebase' in args or (self.basecurrent != self.current.path):
 406             rebase=True
 407         if rebase:
 408             clicks=self.config['baseline_clicks']
 409             if clicks==0:
 410                 self.basepoints=[]
 411                 self.basepoints.append(ClickedPoint())
 412                 self.basepoints.append(ClickedPoint())
 413                 self.basepoints[0].index=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 414                 self.basepoints[1].index=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 415                 for point in self.basepoints:
 416                     #for graphing etc. purposes, fill-in with coordinates
 417                     point.absolute_coords=displayed_plot.vectors[1][0][point.index], displayed_plot.vectors[1][1][point.index]
 418                     point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 419             elif clicks>0:
 420                 print 'Select baseline'
 421                 if clicks==1:
 422                     self.basepoints=self._measure_N_points(N=1, whatset=whatset)
 423                     self.basepoints.append(ClickedPoint())
 424                     self.basepoints[1].index=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 425                     #for graphing etc. purposes, fill-in with coordinates
 426                     self.basepoints[1].absolute_coords=displayed_plot.vectors[1][0][self.basepoints[1].index], displayed_plot.vectors[1][1][self.basepoints[1].index]
 427                     self.basepoints[1].find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 428                 else:
 429                     self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 430
 431             self.basecurrent=self.current.path
 432
 433         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 434         boundaries.sort()
 435         to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 436         avg=np.mean(to_average)
 437
 438
 439         #Initialize data vectors
 440         c_lengths=[]
 441         p_lengths=[]
 442         forces=[]
 443         slopes=[]
 444
 445
 446
 447         #Cycle between peaks and do analysis.
 448         for peak in peak_location:
 449             #Do WLC fits.
 450             #FIXME: clean wlc fitting, to avoid this clickedpoint hell
 451             #-create a clicked point for the peak point
 452             peak_point=ClickedPoint()
 453             peak_point.absolute_coords=displayed_plot.vectors[1][0][peak], displayed_plot.vectors[1][1][peak]
 454             peak_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 455
 456             if not usepoints:
 457                 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
 458
 459             #-create a clicked point for the other fit point
 460             other_fit_point=ClickedPoint()
 461             other_fit_point.absolute_coords=displayed_plot.vectors[1][0][peak-fit_points], displayed_plot.vectors[1][1][peak-fit_points]
 462             other_fit_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
 463             #do the fit
 464             points=[contact_point, peak_point, other_fit_point]
 465
 466             #Check if we have enough points for a fit. If not, wlc_fit could crash
 467             if abs(peak_point.index-other_fit_point.index) < 2:
 468                 continue
 469
 470             params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
 471             #save wlc values (nm)
 472             c_lengths.append(params[0]*(1.0e+9))
 473             if len(params)==2: #if we did choose 2-value fit
 474                 p_lengths.append(params[1]*(1.0e+9))
 475             else:
 476                 p_lengths.append(pl_value)
 477             #Add WLC fit lines to plot
 478             fitplot.add_set(xfit,yfit)
 479
 480             if len(fitplot.styles)==0:
 481                 fitplot.styles=[]
 482             else:
 483                 fitplot.styles.append(None)
 484
 485             #Measure forces
 486             delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
 487             y=min(delta_to_measure)
 488             #save force values (pN)
 489             forces.append(abs(y-avg)*(1.0e+12))
 490
 491             #Measure slopes
 492             slope=self.linefit_between(peak-slope_span,peak)[0]
 493             slopes.append(slope)
 494
 495         #--DEBUG STUFF--
 496         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]])
 497         fitplot.styles.append('scatter')
 498
 499         #Show wlc fits and peak locations
 500         self._send_plot([fitplot])
 501         self.do_peaks('')
 502
 503         #Ask the user what peaks to ignore from analysis.
 504         print 'Peaks to ignore (0,1...n from contact point,return to take all)'
 505         print 'N to discard measurement'
 506         exclude_raw=raw_input('Input:')
 507         if exclude_raw=='N':
 508             print 'Discarded.'
 509             return
 510         if not exclude_raw=='':
 511             exclude=exclude_raw.split(',')
 512             try:
 513                 exclude=[int(item) for item in exclude]
 514                 for i in exclude:
 515                     c_lengths[i]=None
 516                     p_lengths[i]=None
 517                     forces[i]=None
 518                     slopes[i]=None
 519             except:
 520                  print 'Bad input, taking all...'
 521         #Clean data vectors from ignored peaks
 522         c_lengths=[item for item in c_lengths if item != None]
 523         p_lengths=[item for item in p_lengths if item != None]
 524         forces=[item for item in forces if item != None]
 525         slopes=[item for item in slopes if item != None]
 526         print 'contour (nm)',c_lengths
 527         print 'p (nm)',p_lengths
 528         print 'forces (pN)',forces
 529         print 'slopes (N/m)',slopes
 530
 531         #Save file info
 532         if self.autofile=='':
 533             self.autofile=raw_input('Autopeak filename? (return to ignore) ')
 534             if self.autofile=='':
 535                 print 'Not saved.'
 536                 return
 537
 538         if not os.path.exists(self.autofile):
 539             f=open(self.autofile,'w+')
 540             f.write('Analysis started '+time.asctime()+'\n')
 541             f.write('----------------------------------------\n')
 542             f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) \n')
 543             f.close()
 544
 545         print 'Saving...'
 546         f=open(self.autofile,'a+')
 547
 548         f.write(self.current.path+'\n')
 549         for i in range(len(c_lengths)):
 550             f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
 551
 552         f.close()
 553         self.do_note('autopeak')
 554         '''