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 ..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
  48
  49
  50 class FlatFiltPlugin (Plugin):
  51     """Standard-devitiation-based peak recognition and filtering.
  52     """
  53     def __init__(self):
  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).
  58 """.strip()),
  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'),
  62                           ('stable', 0.005),
  63                           ('max cut', 0.2),
  64                           ('min deviation', 9.0),
  65                           ('min points', 4),
  66                           ('see double', 10e-9),
  67                           ]:
  68             argument = [a for a in self._arguments if a.name == key][0]
  69             argument.default = value
  70         self._settings = [
  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)]
  77
  78     def dependencies(self):
  79         return ['vclamp']
  80
  81     def default_settings(self):
  82         return self._settings
  83
  84
  85 class FlatPeaksCommand (Command):
  86     """Detect peaks in velocity clamp data using noise statistics.
  87
  88     Notes
  89     -----
  90
  91     Noise analysis on the retraction curve:
  92
  93     1) A median window filter (using
  94       :func:`scipy.signal.signaltools.medfilt`) smooths the
  95       deflection.
  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`.
 100
 101     The algorithm was originally Francesco Musiani's idea.
 102     """
 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)
 112
 113     def _run(self, hooke, inqueue, outqueue, params):
 114         z_data,d_data,params = self._setup(params)
 115         start_index = 0
 116         while z_data[start_index] < params['bind window']:
 117             start_index += 1
 118         median = medfilt(d_data[start_index:], params['median window']),
 119         deriv = diff(median)
 120         kwargs = dict([(a.name, params[a.name]) for a in find_peaks_arguments])
 121         peaks = find_peaks(deriv, **kwargs)
 122         for peak in peaks:
 123             peak.name = 'flat filter of %s with %s' \
 124                 % (params['deflection column name'], params['convolution'])
 125             peak.index += start_index
 126         outqueue.put(peaks)
 127
 128     def _setup(self, params):
 129         """Setup `params` from config and return the z piezo and
 130         deflection arrays.
 131         """
 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']:
 138                 scale(curve)
 139         data = None
 140         for block in curve.data:
 141             if block.info['name'].startswith('retract'):
 142                 data = block
 143                 break
 144         if data == None:
 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)'
 149         else:
 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
 157
 158 class FlatFilterCommand (FilterCommand):
 159     u"""Create a subset playlist of curves with enough flat peaks.
 160
 161     Notes
 162     -----
 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%.
 168
 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 :math:`\alpha`-Synuclein at the
 172       single molecule level."
 173       PLOS Biology, 2009.
 174       doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
 175
 176     See Also
 177     --------
 178     FlatCommand : Underlying flat-based peak detection.
 179     """
 180     def __init__(self, plugin):
 181         super(FlatFilterCommand, self).__init__(
 182             plugin, name='flat filter playlist')
 183         self.arguments.extend(plugin._arguments)
 184
 185     def filter(self, curve, hooke, inqueue, outqueue, params):
 186         inq = Queue()
 187         outq = Queue()
 188         conv_command = [c for c in self.hooke.commands
 189                         if c.name=='flat filter peaks'][0]
 190         conv_command.run(hooke, inq, outq, **params)
 191         peaks = outq.get()
 192         if not isinstance(peaks[0], Peak):
 193             raise Failure('Expected a list of Peaks, not %s' % peaks)
 194         ret = outq.get()
 195         if not isinstance(ret, Success):
 196             raise ret
 197         return len(peaks) >= params['min peaks']
 198