pcluster.py

   1 #!/usr/bin/env python
   2
   3 from libhooke import WX_GOOD, ClickedPoint
   4 import wxversion
   5 wxversion.select(WX_GOOD)
   6 from wx import PostEvent
   7 import numpy as np
   8 import scipy as sp
   9 import copy
  10 import os.path
  11 import time
  12
  13 import warnings
  14 warnings.simplefilter('ignore',np.RankWarning)
  15
  16
  17 class pclusterCommands:
  18
  19     def do_pcluster(self,args):
  20         '''
  21         pCLUSTER
  22         (pcluster.py)
  23         Automatically measures peaks and extracts informations for further clustering
  24         (c)Paolo Pancaldi, Massimo Sandal 2009
  25         '''
  26         #--Custom persistent length
  27         pl_value=None
  28         for arg in args.split():
  29             #look for a persistent length argument.
  30             if 'pl=' in arg:
  31                 pl_expression=arg.split('=')
  32                 pl_value=float(pl_expression[1]) #actual value
  33             else:
  34                 pl_value=None
  35
  36         #configuration variables
  37         min_npks = self.convfilt_config['minpeaks']
  38         min_deviation = self.convfilt_config['mindeviation']
  39
  40         pclust_filename=raw_input('Automeasure filename? ')
  41         realclust_filename=raw_input('Coordinates filename? ')
  42
  43         f=open(pclust_filename,'w+')
  44         f.write('Analysis started '+time.asctime()+'\n')
  45         f.write('----------------------------------------\n')
  46         f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
  47         f.close()
  48
  49         f=open(realclust_filename,'w+')
  50         f.write('Analysis started '+time.asctime()+'\n')
  51         f.write('----------------------------------------\n')
  52         f.write('; Peak number ; Mean delta (nm)  ;  Median delta (nm) ;  Mean force (pN)  ;  Median force (pN) ; First peak length (nm) ; Last peak length (nm) ; Max force (pN) ; Min force (pN) ; Max delta (nm) ; Min delta (nm)\n')
  53         f.close()
  54         # ------ FUNCTION ------
  55         def fit_interval_nm(start_index,plot,nm,backwards):
  56             '''
  57             Calculates the number of points to fit, given a fit interval in nm
  58             start_index: index of point
  59             plot: plot to use
  60             backwards: if true, finds a point backwards.
  61             '''
  62             whatset=1 #FIXME: should be decidable
  63             x_vect=plot.vectors[1][0]
  64
  65             c=0
  66             i=start_index
  67             start=x_vect[start_index]
  68             maxlen=len(x_vect)
  69             while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
  70                 if i==0 or i==maxlen-1: #we reached boundaries of vector!
  71                     return c
  72                 if backwards:
  73                     i-=1
  74                 else:
  75                     i+=1
  76                 c+=1
  77             return c
  78
  79         def plot_informations(itplot,pl_value):
  80             '''
  81             OUR VARIABLES
  82             contact_point.absolute_coords               (2.4584142802103689e-007, -6.9647135616234017e-009)
  83             peak_point.absolute_coords                  (3.6047748250571423e-008, -7.7142802788854212e-009)
  84             other_fit_point.absolute_coords     (4.1666139243838867e-008, -7.3759393477579707e-009)
  85             peak_location                                                                               [510, 610, 703, 810, 915, 1103]
  86             peak_size                                                                                           [-1.2729111505202212e-009, -9.1632775347399312e-010, -8.1707438353929907e-010, -8.0335812578148904e-010, -8.7483955226387558e-010, -3.6269619757067322e-009]
  87             params                                                                                                      [2.2433999931959462e-007, 3.3230248825175678e-010]
  88             fit_errors                                                                                  [6.5817195369767644e-010, 2.4415923138871498e-011]
  89             '''
  90             fit_points=int(self.config['auto_fit_points']) # number of points to fit before the peak maximum <50>
  91
  92             T=self.config['temperature'] #temperature of the system in kelvins. By default it is 293 K. <301.0>
  93             cindex=self.find_contact_point() #Automatically find contact point <158, libhooke.ClickedPoint>
  94             contact_point=self._clickize(itplot[0].vectors[1][0], itplot[0].vectors[1][1], cindex)
  95             self.basepoints=[]
  96             base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], itplot[0], self.config['auto_right_baseline'],False)
  97             self.basepoints.append(self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],base_index_0))
  98             base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, itplot[0], self.config['auto_left_baseline'],False)
  99             self.basepoints.append(self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],base_index_1))
 100             self.basecurrent=self.current.path
 101             boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 102             boundaries.sort()
 103             to_average=itplot[0].vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 104             avg=np.mean(to_average)
 105             return fit_points, contact_point, pl_value, T, cindex, avg
 106
 107         def features_peaks(itplot, peak, fit_points, contact_point, pl_value, T, cindex, avg):
 108             '''
 109             calculate informations for each peak and add they in
 110             c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes
 111             '''
 112             c_leng=None
 113             p_leng=None
 114             sigma_c_leng=None
 115             sigma_p_leng=None
 116             force=None
 117             slope=None
 118
 119             delta_force=10
 120             slope_span=int(self.config['auto_slope_span'])
 121
 122             peak_point=self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],peak)
 123             other_fit_point=self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],peak-fit_points)
 124
 125             points=[contact_point, peak_point, other_fit_point]
 126
 127             params, yfit, xfit, fit_errors = self.wlc_fit(points, itplot[0].vectors[1][0], itplot[0].vectors[1][1], pl_value, T, return_errors=True)
 128
 129             #Measure forces
 130             delta_to_measure=itplot[0].vectors[1][1][peak-delta_force:peak+delta_force]
 131             y=min(delta_to_measure)
 132             #Measure slopes
 133             slope=self.linefit_between(peak-slope_span,peak)[0]
 134             #check fitted data and, if right, add peak to the measurement
 135             if len(params)==1: #if we did choose 1-value fit
 136                 p_leng=pl_value
 137                 c_leng=params[0]*(1.0e+9)
 138                 sigma_p_lengths=0
 139                 sigma_c_lengths=fit_errors[0]*(1.0e+9)
 140                 force = abs(y-avg)*(1.0e+12)
 141             else: #2-value fit
 142                 p_leng=params[1]*(1.0e+9)
 143                 #check if persistent length makes sense. otherwise, discard peak.
 144                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 145                     '''
 146                     p_lengths.append(p_leng)
 147                     c_lengths.append(params[0]*(1.0e+9))
 148                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 149                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 150                     forces.append(abs(y-avg)*(1.0e+12))
 151                     slopes.append(slope)
 152                     '''
 153                     c_leng=params[0]*(1.0e+9)
 154                     sigma_c_leng=fit_errors[0]*(1.0e+9)
 155                     sigma_p_leng=fit_errors[1]*(1.0e+9)
 156                     force=abs(y-avg)*(1.0e+12)
 157                 else:
 158                     p_leng=None
 159                     slope=None
 160             #return c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes
 161             return  c_leng, p_leng, sigma_c_leng, sigma_p_leng, force, slope
 162
 163
 164         # ------ PROGRAM -------
 165         c=0
 166         for item in self.current_list:
 167             c+=1
 168             item.identify(self.drivers)
 169             itplot=item.curve.default_plots()
 170             try:
 171                 peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
 172             except:
 173                 #We have troubles with exec_has_peaks (bad curve, whatever).
 174                 #Print info and go to next cycle.
 175                 print 'Cannot process ',item.path
 176                 continue
 177
 178             if len(peak_location)==0:
 179                 continue
 180
 181             fit_points, contact_point, pl_value, T, cindex, avg = plot_informations(itplot,pl_value)
 182             print '\n\nCurve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'
 183
 184             #initialize output data vectors
 185             c_lengths=[]
 186             p_lengths=[]
 187             sigma_c_lengths=[]
 188             sigma_p_lengths=[]
 189             forces=[]
 190             slopes=[]
 191
 192             #loop each peak of my curve
 193             for peak in peak_location:
 194                 c_leng, p_leng, sigma_c_leng, sigma_p_leng, force, slope = features_peaks(itplot, peak, fit_points, contact_point, pl_value, T, cindex, avg)
 195                 for var, vector in zip([c_leng, p_leng, sigma_c_leng, sigma_p_leng, force, slope],[c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes]):
 196                     if var is not None:
 197                         vector.append(var)
 198
 199             #FIXME: We need a dictionary here...
 200             allvects=[c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes]
 201             for vect in allvects:
 202                 if len(vect)==0:
 203                     for i in range(len(c_lengths)):
 204                         vect.append(0)
 205
 206             print 'Measurements for all peaks detected:'
 207             print 'contour (nm)', c_lengths
 208             print 'sigma contour (nm)',sigma_c_lengths
 209             print 'p (nm)',p_lengths
 210             print 'sigma p (nm)',sigma_p_lengths
 211             print 'forces (pN)',forces
 212             print 'slopes (N/m)',slopes
 213
 214             '''
 215             write automeasure text file
 216             '''
 217             print 'Saving automatic measurement...'
 218             f=open(pclust_filename,'a+')
 219             f.write(item.path+'\n')
 220             for i in range(len(c_lengths)):
 221                 f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+' ; '+str(sigma_c_lengths[i])+' ; '+str(sigma_p_lengths[i])+'\n')
 222             f.close()
 223
 224             '''
 225             calculate clustering coordinates
 226             '''
 227             peak_number=len(c_lengths)
 228             if peak_number > 1:
 229                 deltas=[]
 230                 for i in range(len(c_lengths)-1):
 231                     deltas.append(c_lengths[i+1]-c_lengths[i])
 232
 233                 delta_mean=np.mean(deltas)
 234                 delta_median=np.median(deltas)
 235
 236                 force_mean=np.mean(forces)
 237                 force_median=np.median(forces)
 238
 239                 first_peak_cl=c_lengths[0]
 240                 last_peak_cl=c_lengths[-1]
 241
 242                 max_force=max(forces[:-1])
 243                 min_force=min(forces)
 244
 245                 max_delta=max(deltas)
 246                 min_delta=min(deltas)
 247
 248                 print 'Coordinates'
 249                 print 'Peaks',peak_number
 250                 print 'Mean delta',delta_mean
 251                 print 'Median delta',delta_median
 252                 print 'Mean force',force_mean
 253                 print 'Median force',force_median
 254                 print 'First peak',first_peak_cl
 255                 print 'Last peak',last_peak_cl
 256                 print 'Max force',max_force
 257                 print 'Min force',min_force
 258                 print 'Max delta',max_delta
 259                 print 'Min delta',min_delta
 260
 261                 '''
 262                 write clustering coordinates
 263                 '''
 264                 f=open(realclust_filename,'a+')
 265                 f.write(item.path+'\n')
 266                 f.write(' ; '+str(peak_number)+' ; '+str(delta_mean)+' ; '+str(delta_median)+' ; '+str(force_mean)+' ; '+str(force_median)+' ; '+str(first_peak_cl)+' ; '+str(last_peak_cl)+ ' ; '+str(max_force)+' ; '
 267                 +str(min_force)+' ; '+str(max_delta)+' ; '+str(min_delta)+ '\n')
 268                 f.close()
 269             else:
 270                 pass