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