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 z_data[start_index] < params['blind window']:
 127             start_index += 1
 128         median = medfilt(d_data[start_index:], params['median window'])
 129         deriv = diff(median)
 130         peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
 131                                             argument_input_keys=True))
 132         for peak in peaks:
 133             peak.name = 'flat filter of %s' % (params['deflection column name'])
 134             peak.index += start_index
 135         outqueue.put(peaks)
 136
 137     def _setup(self, params):
 138         """Setup `params` from config and return the z piezo and
 139         deflection arrays.
 140         """
 141         curve = params['curve']
 142         if curve.info['experiment'] != VelocityClamp:
 143             raise Failure('%s operates on VelocityClamp experiments, not %s'
 144                           % (self.name, curve.info['experiment']))
 145         for col in ['surface z piezo (m)', 'deflection (N)']:
 146             if col not in curve.data[0].info['columns']:
 147                 scale(curve)
 148         data = None
 149         for block in curve.data:
 150             if block.info['name'].startswith('retract'):
 151                 data = block
 152                 break
 153         if data == None:
 154             raise Failure('No retraction blocks in %s.' % curve)
 155         z_data = data[:,data.info['columns'].index('surface z piezo (m)')]
 156         if 'flattened deflection (N)' in data.info['columns']:
 157             params['deflection column name'] = 'flattened deflection (N)'
 158         else:
 159             params['deflection column name'] = 'deflection (N)'
 160         d_data = data[:,data.info['columns'].index(
 161                 params['deflection column name'])]
 162         for key,value in params.items():
 163             if value == None: # Use configured default value.
 164                 params[key] = self.plugin.config[key]
 165         # TODO: better option parser to do this automatically by Argument.type
 166         for key in ['blind window', 'median window', 'max cut', 'min deviations', 'min points', 'see double', 'stable']:
 167             params[key] = float(params[key])
 168         # TODO: convert 'see double' from nm to points
 169         return z_data,d_data,params
 170
 171 class FlatFilterCommand (FilterCommand):
 172     u"""Create a subset playlist of curves with enough flat peaks.
 173
 174     Notes
 175     -----
 176     This type of filter is of course very raw, and requires relatively
 177     conservative settings to safely avoid false negatives (that is, to
 178     avoid discarding interesting curves).  Using it on the protein
 179     unfolding experiments described by Sandal [#sandal2008] it has
 180     been found to reduce the data set to analyze by hand by 60-80%.
 181
 182     .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
 183       F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
 184       "Conformational equilibria in monomeric α-Synuclein at the
 185       single molecule level."
 186       PLOS Biology, 2009.
 187       doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
 188
 189     See Also
 190     --------
 191     FlatCommand : Underlying flat-based peak detection.
 192     """
 193     def __init__(self, plugin):
 194         super(FlatFilterCommand, self).__init__(
 195             plugin, name='flat filter playlist')
 196         self.arguments.extend(plugin._arguments)
 197
 198     def filter(self, curve, hooke, inqueue, outqueue, params):
 199         inq = Queue()
 200         outq = Queue()
 201         filt_command = [c for c in hooke.commands
 202                         if c.name=='flat filter peaks'][0]
 203         filt_command.run(hooke, inq, outq, **params)
 204         peaks = outq.get()
 205         if not (isinstance(peaks, list) and (len(peaks) == 0
 206                                              or isinstance(peaks[0], Peak))):
 207             raise Failure('Expected a list of Peaks, not %s' % peaks)
 208         ret = outq.get()
 209         if not isinstance(ret, Success):
 210             raise ret
 211         return len(peaks) >= int(params['min peaks'])