hooke/plugin/flatfilt.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2010 Alberto Gomez-Casado
   4 #                         Fabrizio Benedetti
   5 #                         Massimo Sandal <devicerandom@gmail.com>
   6 #                         W. Trevor King <wking@drexel.edu>
   7 #
   8 # This file is part of Hooke.
   9 #
  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.
  14 #
  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.
  19 #
  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/>.
  23
  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.
  28
  29 See Also
  30 --------
  31 :mod:`~hooke.plugin.convfilt` for a convolution-based filter for
  32 :class:`~hooke.experiment.VelocityClamp` experiments.
  33 """
  34
  35 import copy
  36 from multiprocessing import Queue
  37
  38 from numpy import diff
  39 from scipy.signal.signaltools import medfilt
  40
  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
  50
  51
  52 class FlatFiltPlugin (Plugin):
  53     """Standard-devitiation-based peak recognition and filtering.
  54     """
  55     def __init__(self):
  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).
  60 """.strip()),
  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.
  64 """.strip()),
  65             Argument('min peaks', type='int', default=4, help="""
  66 Minimum number of peaks for curve acceptance.
  67 """.strip()),
  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'),
  71                           ('stable', 0.005),
  72                           ('max cut', 0.2),
  73                           ('min deviations', 9.0),
  74                           ('min points', 4),
  75                           ('see double', 10e-9),
  76                           ]:
  77             argument = [a for a in self._arguments if a.name == key][0]
  78             argument.default = value
  79         self._settings = [
  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)]
  86
  87     def dependencies(self):
  88         return ['vclamp']
  89
  90     def default_settings(self):
  91         return self._settings
  92
  93
  94 class FlatPeaksCommand (Command):
  95     """Detect peaks in velocity clamp data using noise statistics.
  96
  97     Notes
  98     -----
  99
 100     Noise analysis on the retraction curve:
 101
 102     1) A median window filter (using
 103       :func:`scipy.signal.signaltools.medfilt`) smooths the
 104       deflection.
 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`.
 109
 110     The algorithm was originally Francesco Musiani's idea.
 111     """
 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',
 119             arguments=[
 120                 CurveArgument,
 121                 ] + config_arguments,
 122             help=self.__doc__, plugin=plugin)
 123
 124     def _run(self, hooke, inqueue, outqueue, params):
 125         z_data,d_data,params = self._setup(hooke, params)
 126         start_index = 0
 127         while (start_index < len(z_data)
 128                and z_data[start_index] < params['blind window']):
 129             start_index += 1
 130         median = medfilt(d_data[start_index:], params['median window'])
 131         deriv = diff(median)
 132         peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
 133                                             argument_input_keys=True))
 134         for peak in peaks:
 135             peak.name = 'flat filter of %s' % (params['deflection column name'])
 136             peak.index += start_index
 137         outqueue.put(peaks)
 138
 139     def _setup(self, hooke, params):
 140         """Setup `params` from config and return the z piezo and
 141         deflection arrays.
 142         """
 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']:
 149                 scale(hooke, curve)
 150         data = None
 151         for i,block in enumerate(curve.data):
 152             if block.info['name'].startswith('retract'):
 153                 data = block
 154                 break
 155         if data == None:
 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)'
 160         else:
 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
 172
 173 class FlatFilterCommand (FilterCommand):
 174     u"""Create a subset playlist of curves with enough flat peaks.
 175
 176     Notes
 177     -----
 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%.
 183
 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."
 188       PLOS Biology, 2009.
 189       doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
 190
 191     See Also
 192     --------
 193     FlatCommand : Underlying flat-based peak detection.
 194     """
 195     def __init__(self, plugin):
 196         super(FlatFilterCommand, self).__init__(
 197             plugin, name='flat filter playlist')
 198         self.arguments.extend(plugin._arguments)
 199
 200     def filter(self, curve, hooke, inqueue, outqueue, params):
 201         params['curve'] = curve
 202         inq = Queue()
 203         outq = Queue()
 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)
 207         peaks = outq.get()
 208         if isinstance(peaks, UncaughtException) \
 209                 and isinstance(peaks.exception, PoorFit):
 210             return False
 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)
 214         ret = outq.get()
 215         if not isinstance(ret, Success):
 216             raise ret
 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'])