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
  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.peak import find_peaks, find_peaks_arguments, Peak, _kwargs
  49
  50
  51 class FlatFiltPlugin (Plugin):
  52     """Standard-devitiation-based peak recognition and filtering.
  53     """
  54     def __init__(self):
  55         super(FlatFiltPlugin, self).__init__(name='flatfilt')
  56         self._arguments = [ # For Command initialization
  57             Argument('median window', type='int', default=7, help="""
  58 Median window filter size (in points).
  59 """.strip()),
  60             Argument('blind window', type='float', default=20e-9, help="""
  61 Meters after the contact point where we do not count peaks to avoid
  62 non-specific surface interaction.
  63 """.strip()),
  64             Argument('min peaks', type='int', default=4, help="""
  65 Minimum number of peaks for curve acceptance.
  66 """.strip()),
  67             ] + copy.deepcopy(find_peaks_arguments)
  68         # Set flat-filter-specific defaults for the fit_peak_arguments.
  69         for key,value in [('cut side', 'both'),
  70                           ('stable', 0.005),
  71                           ('max cut', 0.2),
  72                           ('min deviations', 9.0),
  73                           ('min points', 4),
  74                           ('see double', 10e-9),
  75                           ]:
  76             argument = [a for a in self._arguments if a.name == key][0]
  77             argument.default = value
  78         self._settings = [
  79             Setting(section=self.setting_section, help=self.__doc__)]
  80         for argument in self._arguments:
  81             self._settings.append(argument_to_setting(
  82                     self.setting_section, argument))
  83             argument.default = None # if argument isn't given, use the config.
  84         self._commands = [FlatPeaksCommand(self), FlatFilterCommand(self)]
  85
  86     def dependencies(self):
  87         return ['vclamp']
  88
  89     def default_settings(self):
  90         return self._settings
  91
  92
  93 class FlatPeaksCommand (Command):
  94     """Detect peaks in velocity clamp data using noise statistics.
  95
  96     Notes
  97     -----
  98
  99     Noise analysis on the retraction curve:
 100
 101     1) A median window filter (using
 102       :func:`scipy.signal.signaltools.medfilt`) smooths the
 103       deflection.
 104     2) The deflection derivative is calculated (using
 105       :func:`numpy.diff` which uses forward differencing).
 106     3) Peaks in the derivative curve are extracted with
 107       :func:`~hooke.plugins.peak.find_peaks`.
 108
 109     The algorithm was originally Francesco Musiani's idea.
 110     """
 111     def __init__(self, plugin):
 112         config_arguments = [a for a in plugin._arguments
 113                             if a.name != 'min peaks']
 114         # Drop min peaks, since we're not filtering with this
 115         # function, just detecting peaks.
 116         super(FlatPeaksCommand, self).__init__(
 117             name='flat filter peaks',
 118             arguments=[
 119                 CurveArgument,
 120                 ] + config_arguments,
 121             help=self.__doc__, plugin=plugin)
 122
 123     def _run(self, hooke, inqueue, outqueue, params):
 124         z_data,d_data,params = self._setup(params)
 125         start_index = 0
 126         while (start_index < len(z_data)
 127                and z_data[start_index] < params['blind window']):
 128             start_index += 1
 129         median = medfilt(d_data[start_index:], params['median window'])
 130         deriv = diff(median)
 131         peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
 132                                             argument_input_keys=True))
 133         for peak in peaks:
 134             peak.name = 'flat filter of %s' % (params['deflection column name'])
 135             peak.index += start_index
 136         outqueue.put(peaks)
 137
 138     def _setup(self, params):
 139         """Setup `params` from config and return the z piezo and
 140         deflection arrays.
 141         """
 142         curve = params['curve']
 143         if curve.info['experiment'] != VelocityClamp:
 144             raise Failure('%s operates on VelocityClamp experiments, not %s'
 145                           % (self.name, curve.info['experiment']))
 146         for col in ['surface z piezo (m)', 'deflection (N)']:
 147             if col not in curve.data[0].info['columns']:
 148                 scale(curve)
 149         data = None
 150         for block in curve.data:
 151             if block.info['name'].startswith('retract'):
 152                 data = block
 153                 break
 154         if data == None:
 155             raise Failure('No retraction blocks in %s.' % curve)
 156         z_data = data[:,data.info['columns'].index('surface z piezo (m)')]
 157         if 'flattened deflection (N)' in data.info['columns']:
 158             params['deflection column name'] = 'flattened deflection (N)'
 159         else:
 160             params['deflection column name'] = 'deflection (N)'
 161         d_data = data[:,data.info['columns'].index(
 162                 params['deflection column name'])]
 163         for key,value in params.items():
 164             if value == None: # Use configured default value.
 165                 params[key] = self.plugin.config[key]
 166         # TODO: better option parser to do this automatically by Argument.type
 167         for key in ['blind window', 'median window', 'max cut', 'min deviations', 'min points', 'see double', 'stable']:
 168             params[key] = float(params[key])
 169         # TODO: convert 'see double' from nm to points
 170         return z_data,d_data,params
 171
 172 class FlatFilterCommand (FilterCommand):
 173     u"""Create a subset playlist of curves with enough flat peaks.
 174
 175     Notes
 176     -----
 177     This type of filter is of course very raw, and requires relatively
 178     conservative settings to safely avoid false negatives (that is, to
 179     avoid discarding interesting curves).  Using it on the protein
 180     unfolding experiments described by Sandal [#sandal2008] it has
 181     been found to reduce the data set to analyze by hand by 60-80%.
 182
 183     .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
 184       F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
 185       "Conformational equilibria in monomeric α-Synuclein at the
 186       single molecule level."
 187       PLOS Biology, 2009.
 188       doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
 189
 190     See Also
 191     --------
 192     FlatCommand : Underlying flat-based peak detection.
 193     """
 194     def __init__(self, plugin):
 195         super(FlatFilterCommand, self).__init__(
 196             plugin, name='flat filter playlist')
 197         self.arguments.extend(plugin._arguments)
 198
 199     def filter(self, curve, hooke, inqueue, outqueue, params):
 200         inq = Queue()
 201         outq = Queue()
 202         filt_command = [c for c in hooke.commands
 203                         if c.name=='flat filter peaks'][0]
 204         filt_command.run(hooke, inq, outq, **params)
 205         peaks = outq.get()
 206         if not (isinstance(peaks, list) and (len(peaks) == 0
 207                                              or isinstance(peaks[0], Peak))):
 208             raise Failure('Expected a list of Peaks, not %s' % peaks)
 209         ret = outq.get()
 210         if not isinstance(ret, Success):
 211             raise ret
 212         if params['min peaks'] == None: # Use configured default value.
 213             params['min peaks'] = self.plugin.config['min peaks']
 214         return len(peaks) >= int(params['min peaks'])