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