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, Failure
42 from ..config import Setting
43 from ..plugin import Plugin, argument_to_setting
44 from ..plugin.curve import CurveArgument
45 from ..plugin.playlist import FilterCommand
46 from ..plugin.vclamp import scale
47 from ..util.peak import find_peaks, find_peaks_arguments, Peak
50 class FlatFiltPlugin (Plugin):
51 """Standard-devitiation-based peak recognition and filtering.
54 super(FlatFiltPlugin, self).__init__(name='flatfilt')
55 self._arguments = [ # For Command initialization
56 Argument('median filter', type='int', default=7, help="""
57 Median window filter size (in points).
59 ] + copy.deepcopy(find_peaks_arguments)
60 # Set flat-filter-specific defaults for the fit_peak_arguments.
61 for key,value in [('cut side', 'both'),
64 ('min deviation', 9.0),
66 ('see double', 10e-9),
68 argument = [a for a in self._arguments if a.name == key][0]
69 argument.default = value
71 Setting(section=self.setting_section, help=self.__doc__)]
72 for argument in self._arguments:
73 self._settings.append(argument_to_setting(
74 self.setting_section, argument))
75 argument.default = None # if argument isn't given, use the config.
76 self._commands = [FlatPeaksCommand(self), FlatFilterCommand(self)]
78 def dependencies(self):
81 def default_settings(self):
85 class FlatPeaksCommand (Command):
86 """Detect peaks in velocity clamp data using noise statistics.
91 Noise analysis on the retraction curve:
93 1) A median window filter (using
94 :func:`scipy.signal.signaltools.medfilt`) smooths the
96 2) The deflection derivative is calculated (using
97 :func:`numpy.diff` which uses forward differencing).
98 3) Peaks in the derivative curve are extracted with
99 :func:`~hooke.plugins.peak.find_peaks`.
101 The algorithm was originally Francesco Musiani's idea.
103 def __init__(self, plugin):
104 config_arguments = [a for a in plugin._arguments
105 if a.name != 'min peaks']
106 # Drop min peaks, since we're not filtering with this
107 # function, just detecting peaks.
108 super(FlatPeaksCommand, self).__init__(
109 name='flat filter peaks',
110 arguments=[CurveArgument] + config_arguments,
111 help=self.__doc__, plugin=plugin)
113 def _run(self, hooke, inqueue, outqueue, params):
114 z_data,d_data,params = self._setup(params)
116 while z_data[start_index] < params['bind window']:
118 median = medfilt(d_data[start_index:], params['median window']),
120 kwargs = dict([(a.name, params[a.name]) for a in find_peaks_arguments])
121 peaks = find_peaks(deriv, **kwargs)
123 peak.name = 'flat filter of %s with %s' \
124 % (params['deflection column name'], params['convolution'])
125 peak.index += start_index
128 def _setup(self, params):
129 """Setup `params` from config and return the z piezo and
132 curve = params['curve']
133 if curve.info['experiment'] != VelocityClamp:
134 raise Failure('%s operates on VelocityClamp experiments, not %s'
135 % (self.name, curve.info['experiment']))
136 for col in ['surface z piezo (m)', 'deflection (N)']:
137 if col not in curve.data[0].info['columns']:
140 for block in curve.data:
141 if block.info['name'].startswith('retract'):
145 raise Failure('No retraction blocks in %s.' % curve)
146 z_data = data[:,data.info['columns'].index('surface z piezo (m)')]
147 if 'flattened deflection (N)' in data.info['columns']:
148 params['deflection column name'] = 'flattened deflection (N)'
150 params['deflection column name'] = 'deflection (N)'
151 d_data = data[:,data.info['columns'].index(
152 params['deflection column name'])]
153 for key,value in params.items():
154 if value == None: # Use configured default value.
155 params[key] = self.plugin.config[key]
156 return z_data,d_data,params
158 class FlatFilterCommand (FilterCommand):
159 u"""Create a subset playlist of curves with enough flat peaks.
163 This type of filter is of course very raw, and requires relatively
164 conservative settings to safely avoid false negatives (that is, to
165 avoid discarding interesting curves). Using it on the protein
166 unfolding experiments described by Sandal [#sandal2008] it has
167 been found to reduce the data set to analyze by hand by 60-80%.
169 .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
170 F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
171 "Conformational equilibria in monomeric α-Synuclein at the
172 single molecule level."
174 doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
178 FlatCommand : Underlying flat-based peak detection.
180 def __init__(self, plugin):
181 super(FlatFilterCommand, self).__init__(
182 plugin, name='flat filter playlist')
183 self.arguments.extend(plugin._arguments)
185 def filter(self, curve, hooke, inqueue, outqueue, params):
188 conv_command = [c for c in hooke.commands
189 if c.name=='flat filter peaks'][0]
190 conv_command.run(hooke, inq, outq, **params)
192 if not isinstance(peaks[0], Peak):
193 raise Failure('Expected a list of Peaks, not %s' % peaks)
195 if not isinstance(ret, Success):
197 return len(peaks) >= params['min peaks']