1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2010 Alberto Gomez-Casado
5 # Massimo Sandal <devicerandom@gmail.com>
6 # W. Trevor King <wking@drexel.edu>
8 # This file is part of Hooke.
10 # Hooke is free software: you can redistribute it and/or
11 # modify it under the terms of the GNU Lesser General Public
12 # License as published by the Free Software Foundation, either
13 # version 3 of the License, or (at your option) any later version.
15 # Hooke is distributed in the hope that it will be useful,
16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 # GNU Lesser General Public License for more details.
20 # You should have received a copy of the GNU Lesser General Public
21 # License along with Hooke. If not, see
22 # <http://www.gnu.org/licenses/>.
24 """The ``flatfilt`` module provides :class:`~FlatFiltPlugin` and
25 several associated :class:`~hooke.command.Command`\s for removing flat
26 (featureless) :mod:`~hooke.curve.Curve`\s from
27 :class:`~hooke.playlist.Playlist`\s.
31 :mod:`~hooke.plugin.convfilt` for a convolution-based filter for
32 :class:`~hooke.experiment.VelocityClamp` experiments.
36 from multiprocessing import Queue
38 from numpy import diff
39 from scipy.signal.signaltools import medfilt
41 from ..command import Command, Argument, Success, Failure, UncaughtException
42 from ..config import Setting
43 from ..experiment import VelocityClamp
44 from ..plugin import Plugin, argument_to_setting
45 from ..plugin.curve import CurveArgument
46 from ..plugin.playlist import FilterCommand
47 from ..plugin.vclamp import scale
48 from ..util.fit import PoorFit
49 from ..util.peak import find_peaks, find_peaks_arguments, Peak, _kwargs
52 class FlatFiltPlugin (Plugin):
53 """Standard-devitiation-based peak recognition and filtering.
56 super(FlatFiltPlugin, self).__init__(name='flatfilt')
57 self._arguments = [ # For Command initialization
58 Argument('median window', type='int', default=7, help="""
59 Median window filter size (in points).
61 Argument('blind window', type='float', default=20e-9, help="""
62 Meters after the contact point where we do not count peaks to avoid
63 non-specific surface interaction.
65 Argument('min peaks', type='int', default=4, help="""
66 Minimum number of peaks for curve acceptance.
68 ] + copy.deepcopy(find_peaks_arguments)
69 # Set flat-filter-specific defaults for the fit_peak_arguments.
70 for key,value in [('cut side', 'both'),
73 ('min deviations', 9.0),
75 ('see double', 10e-9),
77 argument = [a for a in self._arguments if a.name == key][0]
78 argument.default = value
80 Setting(section=self.setting_section, help=self.__doc__)]
81 for argument in self._arguments:
82 self._settings.append(argument_to_setting(
83 self.setting_section, argument))
84 argument.default = None # if argument isn't given, use the config.
85 self._commands = [FlatPeaksCommand(self), FlatFilterCommand(self)]
87 def dependencies(self):
90 def default_settings(self):
94 class FlatPeaksCommand (Command):
95 """Detect peaks in velocity clamp data using noise statistics.
100 Noise analysis on the retraction curve:
102 1) A median window filter (using
103 :func:`scipy.signal.signaltools.medfilt`) smooths the
105 2) The deflection derivative is calculated (using
106 :func:`numpy.diff` which uses forward differencing).
107 3) Peaks in the derivative curve are extracted with
108 :func:`~hooke.plugins.peak.find_peaks`.
110 The algorithm was originally Francesco Musiani's idea.
112 def __init__(self, plugin):
113 config_arguments = [a for a in plugin._arguments
114 if a.name != 'min peaks']
115 # Drop min peaks, since we're not filtering with this
116 # function, just detecting peaks.
117 super(FlatPeaksCommand, self).__init__(
118 name='flat filter peaks',
121 ] + config_arguments,
122 help=self.__doc__, plugin=plugin)
124 def _run(self, hooke, inqueue, outqueue, params):
125 z_data,d_data,params = self._setup(hooke, params)
127 while (start_index < len(z_data)
128 and z_data[start_index] < params['blind window']):
130 median = medfilt(d_data[start_index:], params['median window'])
132 peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
133 argument_input_keys=True))
135 peak.name = 'flat filter of %s' % (params['deflection column name'])
136 peak.index += start_index
139 def _setup(self, hooke, params):
140 """Setup `params` from config and return the z piezo and
143 curve = params['curve']
144 if curve.info['experiment'] != VelocityClamp:
145 raise Failure('%s operates on VelocityClamp experiments, not %s'
146 % (self.name, curve.info['experiment']))
147 for col in ['surface distance (m)', 'deflection (N)']:
148 if col not in curve.data[0].info['columns']:
151 for i,block in enumerate(curve.data):
152 if block.info['name'].startswith('retract'):
156 raise Failure('No retraction blocks in %s.' % curve)
157 z_data = data[:,data.info['columns'].index('surface distance (m)')]
158 if 'flattened deflection (N)' in data.info['columns']:
159 params['deflection column name'] = 'flattened deflection (N)'
161 params['deflection column name'] = 'deflection (N)'
162 d_data = data[:,data.info['columns'].index(
163 params['deflection column name'])]
164 for key,value in params.items():
165 if value == None: # Use configured default value.
166 params[key] = self.plugin.config[key]
167 # TODO: better option parser to do this automatically by Argument.type
168 for key in ['blind window', 'median window', 'max cut', 'min deviations', 'min points', 'see double', 'stable']:
169 params[key] = float(params[key])
170 # TODO: convert 'see double' from nm to points
171 return z_data,d_data,params
173 class FlatFilterCommand (FilterCommand):
174 u"""Create a subset playlist of curves with enough flat peaks.
178 This type of filter is of course very raw, and requires relatively
179 conservative settings to safely avoid false negatives (that is, to
180 avoid discarding interesting curves). Using it on the protein
181 unfolding experiments described by Sandal [#sandal2008] it has
182 been found to reduce the data set to analyze by hand by 60-80%.
184 .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
185 F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
186 "Conformational equilibria in monomeric α-Synuclein at the
187 single molecule level."
189 doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
193 FlatCommand : Underlying flat-based peak detection.
195 def __init__(self, plugin):
196 super(FlatFilterCommand, self).__init__(
197 plugin, name='flat filter playlist')
198 self.arguments.extend(plugin._arguments)
200 def filter(self, curve, hooke, inqueue, outqueue, params):
201 params['curve'] = curve
204 filt_command = [c for c in hooke.commands
205 if c.name=='flat filter peaks'][0]
206 filt_command.run(hooke, inq, outq, **params)
208 if isinstance(peaks, UncaughtException) \
209 and isinstance(peaks.exception, PoorFit):
211 if not (isinstance(peaks, list) and (len(peaks) == 0
212 or isinstance(peaks[0], Peak))):
213 raise Failure('Expected a list of Peaks, not %s' % peaks)
215 if not isinstance(ret, Success):
217 if params['min peaks'] == None: # Use configured default value.
218 params['min peaks'] = self.plugin.config['min peaks']
219 return len(peaks) >= int(params['min peaks'])