From ab63cbcdd5c7d919441df2ecb9713544933f58c1 Mon Sep 17 00:00:00 2001 From: devicerandom Date: Mon, 27 Oct 2008 15:53:15 +0000 Subject: [PATCH] autopeak now in autopeak.py and eliminated from generalvclamp.py --- autopeak.py | 315 +++++++++++++++++++++++++++++++++++++++++++++++ generalvclamp.py | 63 +--------- 2 files changed, 320 insertions(+), 58 deletions(-) create mode 100644 autopeak.py diff --git a/autopeak.py b/autopeak.py new file mode 100644 index 0000000..44c8a60 --- /dev/null +++ b/autopeak.py @@ -0,0 +1,315 @@ +#!/usr/bin/env python + +from libhooke import WX_GOOD, ClickedPoint +import wxversion +wxversion.select(WX_GOOD) +from wx import PostEvent +import numpy as np +import scipy as sp +import copy +import os.path +import time + +import warnings +warnings.simplefilter('ignore',np.RankWarning) + + +class autopeakCommands: + + def do_autopeak(self,args): + ''' + AUTOPEAK + (autopeak.py) + Automatically performs a number of analyses on the peaks of the given curve. + Currently it automatically: + - fits peaks with WLC function + - measures peak maximum forces with a baseline + - measures slope in proximity of peak maximum + Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt + + Syntax: + autopeak [rebase] [pl=value] [t=value] [noauto] [reclick] + + rebase : Re-asks baseline interval + + pl=[value] : Use a fixed persistent length for the fit. If pl is not given, + the fit will be a 2-variable + fit. DO NOT put spaces between 'pl', '=' and the value. + The value must be in meters. + Scientific notation like 0.35e-9 is fine. + + t=[value] : Use a user-defined temperature. The value must be in + kelvins; by default it is 293 K. + DO NOT put spaces between 't', '=' and the value. + + noauto : allows for clicking the contact point by + hand (otherwise it is automatically estimated) the first time. + If subsequent measurements are made, the same contact point + clicked the first time is used + + reclick : redefines by hand the contact point, if noauto has been used before + but the user is unsatisfied of the previously choosen contact point. + + usepoints : fit interval by number of points instead than by nanometers + + When you first issue the command, it will ask for the filename. If you are giving the filename + of an existing file, autopeak will resume it and append measurements to it. If you are giving + a new filename, it will create the file and append to it until you close Hooke. + + + Useful variables (to set with SET command): + --- + temperature= temperature of the system for wlc fit (in K) + + auto_slope_span = number of points on which measure the slope, for slope + + auto_fit_nm = number of nm to fit before the peak maximum, for WLC (if usepoints false) + auto_fit_points = number of points to fit before the peak maximum, for WLC (if usepoints true) + + baseline_clicks = 0: automatic baseline + 1: decide baseline with a single click and length defined in auto_left_baseline + 2: let user click points of baseline + auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1) + auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0) + ''' + + #MACROS. + #FIXME: to move outside function + def fit_interval_nm(start_index,plot,nm,backwards): + ''' + Calculates the number of points to fit, given a fit interval in nm + start_index: index of point + plot: plot to use + backwards: if true, finds a point backwards. + ''' + x_vect=plot.vectors[1][0] + + c=0 + i=start_index + start=x_vect[start_index] + maxlen=len(x_vect) + while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm: + if i==0 or i==maxlen-1: #we reached boundaries of vector! + return c + + if backwards: + i-=1 + else: + i+=1 + c+=1 + return c + + def pickup_contact_point(): + '''macro to pick up the contact point by clicking''' + contact_point=self._measure_N_points(N=1, whatset=1)[0] + contact_point_index=contact_point.index + self.wlccontact_point=contact_point + self.wlccontact_index=contact_point.index + self.wlccurrent=self.current.path + return contact_point, contact_point_index + + def find_current_peaks(): + #Find peaks. + defplot=self.current.curve.default_plots()[0] + flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip + defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks + peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation']) + return peak_location, peak_size + + + + + #default fit etc. variables + pl_value=None + T=self.config['temperature'] + + slope_span=int(self.config['auto_slope_span']) + delta_force=10 + rebase=False #if true=we select rebase + + #initialize output data vectors + c_lengths=[] + p_lengths=[] + forces=[] + slopes=[] + + #pick up plot + displayed_plot=self._get_displayed_plot(0) + + #COMMAND LINE PARSING + #--Using points instead of nm interval + if 'usepoints' in args.split(): + fit_points=int(self.config['auto_fit_points']) + usepoints=True + else: + fit_points=None + usepoints=False + #--Recalculate baseline + if 'rebase' in args or (self.basecurrent != self.current.path): + rebase=True + + #--Custom persistent length / custom temperature + for arg in args.split(): + #look for a persistent length argument. + if 'pl=' in arg: + pl_expression=arg.split('=') + pl_value=float(pl_expression[1]) #actual value + #look for a T argument. FIXME: spaces are not allowed between 'pl' and value + if ('t=' in arg[0:2]) or ('T=' in arg[0:2]): + t_expression=arg.split('=') + T=float(t_expression[1]) + #--Contact point arguments + if 'reclick' in args.split(): + print 'Click contact point' + contact_point, contact_point_index = pickup_contact_point() + elif 'noauto' in args.split(): + if self.wlccontact_index==None or self.wlccurrent != self.current.path: + print 'Click contact point' + contact_point , contact_point_index = pickup_contact_point() + else: + contact_point=self.wlccontact_point + contact_point_index=self.wlccontact_index + else: + #Automatically find contact point + cindex=self.find_contact_point() + contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex) + ''' + contact_point=ClickedPoint() + contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex] + contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1]) + contact_point.is_marker=True + ''' + #--END COMMAND LINE PARSING-- + + + peak_location, peak_size = find_current_peaks() + + fitplot=copy.deepcopy(displayed_plot) + + #Pick up force baseline + if rebase: + clicks=self.config['baseline_clicks'] + if clicks==0: + self.basepoints=[] + base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False) + self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0)) + base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False) + self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1)) + elif clicks>0: + print 'Select baseline' + if clicks==1: + self.basepoints=self._measure_N_points(N=1, whatset=whatset) + base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False) + self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1)) + else: + self.basepoints=self._measure_N_points(N=2, whatset=whatset) + + self.basecurrent=self.current.path + + boundaries=[self.basepoints[0].index, self.basepoints[1].index] + boundaries.sort() + to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average + avg=np.mean(to_average) + + + for peak in peak_location: + #WLC FITTING + #define fit interval + if not usepoints: + fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True) + peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak) + other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points) + + #points for the fit + points=[contact_point, peak_point, other_fit_point] + + if abs(peak_point.index-other_fit_point.index) < 2: + continue + + params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T) + + #save wlc values (nm) + c_lengths.append(params[0]*(1.0e+9)) + if len(params)==2: #if we did choose 2-value fit + p_lengths.append(params[1]*(1.0e+9)) + else: + p_lengths.append(pl_value) + #Add WLC fit lines to plot + fitplot.add_set(xfit,yfit) + + if len(fitplot.styles)==0: + fitplot.styles=[] + else: + fitplot.styles.append(None) + + #Measure forces + delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force] + y=min(delta_to_measure) + #save force values (pN) + forces.append(abs(y-avg)*(1.0e+12)) + + #Measure slopes + slope=self.linefit_between(peak-slope_span,peak)[0] + slopes.append(slope) + + #add basepoints to fitplot + 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]]) + fitplot.styles.append('scatter') + + #Show wlc fits and peak locations + self._send_plot([fitplot]) + self.do_peaks('') + + + #Ask the user what peaks to ignore from analysis. + print 'Peaks to ignore (0,1...n from contact point,return to take all)' + print 'N to discard measurement' + exclude_raw=raw_input('Input:') + if exclude_raw=='N': + print 'Discarded.' + return + if not exclude_raw=='': + exclude=exclude_raw.split(',') + try: + exclude=[int(item) for item in exclude] + for i in exclude: + c_lengths[i]=None + p_lengths[i]=None + forces[i]=None + slopes[i]=None + except: + print 'Bad input, taking all...' + #Clean data vectors from ignored peaks + c_lengths=[item for item in c_lengths if item != None] + p_lengths=[item for item in p_lengths if item != None] + forces=[item for item in forces if item != None] + slopes=[item for item in slopes if item != None] + print 'contour (nm)',c_lengths + print 'p (nm)',p_lengths + print 'forces (pN)',forces + print 'slopes (N/m)',slopes + + #Save file info + if self.autofile=='': + self.autofile=raw_input('Autopeak filename? (return to ignore) ') + if self.autofile=='': + print 'Not saved.' + return + + if not os.path.exists(self.autofile): + f=open(self.autofile,'w+') + f.write('Analysis started '+time.asctime()+'\n') + f.write('----------------------------------------\n') + f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) \n') + f.close() + + print 'Saving...' + f=open(self.autofile,'a+') + + f.write(self.current.path+'\n') + for i in range(len(c_lengths)): + f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n') + + f.close() + self.do_note('autopeak') + \ No newline at end of file diff --git a/generalvclamp.py b/generalvclamp.py index 4039d50..aadbd63 100644 --- a/generalvclamp.py +++ b/generalvclamp.py @@ -284,73 +284,20 @@ class generalvclampCommands: #==================== #AUTOMATIC ANALYSES #==================== + ''' def do_autopeak(self,args): - ''' - AUTOPEAK - (generalvclamp.py) - Automatically performs a number of analyses on the peaks of the given curve. - Currently it automatically: - - fits peaks with WLC function - - measures peak maximum forces with a baseline - - measures slope in proximity of peak maximum - Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt - - Syntax: - autopeak [rebase] [pl=value] [t=value] [noauto] [reclick] - - rebase : Re-asks baseline interval - - pl=[value] : Use a fixed persistent length for the fit. If pl is not given, - the fit will be a 2-variable - fit. DO NOT put spaces between 'pl', '=' and the value. - The value must be in meters. - Scientific notation like 0.35e-9 is fine. - - t=[value] : Use a user-defined temperature. The value must be in - kelvins; by default it is 293 K. - DO NOT put spaces between 't', '=' and the value. - - noauto : allows for clicking the contact point by - hand (otherwise it is automatically estimated) the first time. - If subsequent measurements are made, the same contact point - clicked the first time is used - - reclick : redefines by hand the contact point, if noauto has been used before - but the user is unsatisfied of the previously choosen contact point. - - usepoints : fit interval by number of points instead than by nanometers - - When you first issue the command, it will ask for the filename. If you are giving the filename - of an existing file, autopeak will resume it and append measurements to it. If you are giving - a new filename, it will create the file and append to it until you close Hooke. - - - Useful variables (to set with SET command): - --- - temperature= temperature of the system for wlc fit (in K) - - auto_slope_span = number of points on which measure the slope, for slope - - auto_fit_nm = number of nm to fit before the peak maximum, for WLC (if usepoints false) - auto_fit_points = number of points to fit before the peak maximum, for WLC (if usepoints true) - - baseline_clicks = 0: automatic baseline - 1: decide baseline with a single click and length defined in auto_left_baseline - 2: let user click points of baseline - auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1) - auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0) - ''' - #FIXME: this function is too long. split it and make it rational. #FIXME: also, *generalize fits* to allow FJC and any other model in the future! def fit_interval_nm(start_index,plot,nm,backwards): + ''' ''' Calculates the number of points to fit, given a fit interval in nm start_index: index of point plot: plot to use backwards: if true, finds a point backwards. ''' + ''' x_vect=plot.vectors[1][0] c=0 @@ -403,7 +350,7 @@ class generalvclampCommands: #Handle contact point arguments def pickup_contact_point(): - '''macro to pick up the contact point by clicking''' + #macro to pick up the contact point by clicking contact_point=self._measure_N_points(N=1, whatset=1)[0] contact_point_index=contact_point.index self.wlccontact_point=contact_point @@ -590,4 +537,4 @@ class generalvclampCommands: f.close() self.do_note('autopeak') - \ No newline at end of file + ''' \ No newline at end of file -- 2.26.2