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 modify it
  11 # under the terms of the GNU Lesser General Public License as
  12 # published by the Free Software Foundation, either version 3 of the
  13 # License, or (at your option) any later version.
  14 #
  15 # Hooke is distributed in the hope that it will be useful, but WITHOUT
  16 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  17 # or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  18 # 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 Queue import Queue
  37
  38 from numpy import absolute, diff
  39 from scipy.signal.signaltools import medfilt
  40
  41 from ..command import Argument, Success, Failure, UncaughtException
  42 from ..config import Setting
  43 from ..curve import Data
  44 from ..experiment import VelocityClamp
  45 from ..plugin import Plugin, argument_to_setting
  46 from ..plugin.curve import ColumnAddingCommand
  47 from ..plugin.playlist import FilterCommand
  48 from ..util.fit import PoorFit
  49 from ..util.peak import (find_peaks, peaks_to_mask,
  50                          find_peaks_arguments, Peak, _kwargs)
  51 from ..util.si import join_data_label, split_data_label
  52
  53
  54 class FlatFiltPlugin (Plugin):
  55     """Standard-devitiation-based peak recognition and filtering.
  56     """
  57     def __init__(self):
  58         super(FlatFiltPlugin, self).__init__(name='flatfilt')
  59         self._arguments = [ # For Command initialization
  60             Argument('median window', type='int', default=7, help="""
  61 Median window filter size (in points).
  62 """.strip()),
  63             Argument('blind window', type='float', default=20e-9, help="""
  64 Meters after the contact point where we do not count peaks to avoid
  65 non-specific surface interaction.
  66 """.strip()),
  67             Argument('min peaks', type='int', default=4, help="""
  68 Minimum number of peaks for curve acceptance.
  69 """.strip()),
  70             ] + copy.deepcopy(find_peaks_arguments)
  71         # Set flat-filter-specific defaults for the fit_peak_arguments.
  72         for key,value in [('cut side', 'both'),
  73                           ('stable', 0.005),
  74                           ('max cut', 0.2),
  75                           ('min deviations', 9.0),
  76                           ('min points', 4),
  77                           ('see double', 10), # TODO: points vs. meters. 10e-9),
  78                           ]:
  79             argument = [a for a in self._arguments if a.name == key][0]
  80             argument.default = value
  81         self._settings = [
  82             Setting(section=self.setting_section, help=self.__doc__)]
  83         for argument in self._arguments:
  84             self._settings.append(argument_to_setting(
  85                     self.setting_section, argument))
  86             argument.default = None # if argument isn't given, use the config.
  87         self._commands = [FlatPeaksCommand(self), FlatFilterCommand(self)]
  88
  89     def dependencies(self):
  90         return ['vclamp']
  91
  92     def default_settings(self):
  93         return self._settings
  94
  95
  96 class FlatPeaksCommand (ColumnAddingCommand):
  97     """Detect peaks in velocity clamp data using noise statistics.
  98
  99     Notes
 100     -----
 101     Noise analysis on the retraction curve:
 102
 103     1) A median window filter (using
 104       :func:`scipy.signal.signaltools.medfilt`) smooths the
 105       deflection.
 106     2) The deflection derivative is calculated (using
 107       :func:`numpy.diff` which uses forward differencing).
 108     3) Peaks in the derivative curve are extracted with
 109       :func:`~hooke.plugins.peak.find_peaks`.
 110
 111     The algorithm was originally Francesco Musiani's idea.
 112     """
 113     def __init__(self, plugin):
 114         plugin_arguments = [a for a in plugin._arguments
 115                             if a.name != 'min peaks']
 116         # Drop min peaks, since we're not filtering with this
 117         # function, just detecting peaks.
 118         super(FlatPeaksCommand, self).__init__(
 119             name='flat filter peaks',
 120             columns=[
 121                 ('distance column', 'surface distance (m)', """
 122 Name of the column to use as the surface position input.
 123 """.strip()),
 124                 ('deflection column', 'surface deflection (m)', """
 125 Name of the column to use as the deflection input.
 126 """.strip()),
 127                 ],
 128             new_columns=[
 129                 ('output peak column', 'flat filter peaks', """
 130 Name of the column (without units) to use as the peak output.
 131 """.split()),
 132                 ],
 133             arguments=[
 134                 Argument(name='peak info name', type='string',
 135                          default='flat filter peaks',
 136                          help="""
 137 Name (without units) for storing the list of peaks in the `.info`
 138 dictionary.
 139 """.strip()),
 140                 ] + plugin_arguments,
 141             help=self.__doc__, plugin=plugin)
 142
 143     def _run(self, hooke, inqueue, outqueue, params):
 144         params = self.__setup_params(hooke=hooke, params=params)
 145         block = self._block(hooke=hooke, params=params)
 146         dist_data = self._get_column(hooke=hooke, params=params,
 147                                      column_name='distance column')
 148         def_data = self._get_column(hooke=hooke, params=params,
 149                                      column_name='deflection column')
 150         start_index = absolute(dist_data-params['blind window']).argmin()
 151         median = medfilt(def_data[start_index:], params['median window'])
 152         deriv = diff(median)
 153         peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
 154                                             argument_input_keys=True))
 155         for i,peak in enumerate(peaks):
 156             peak.name = self._peak_name(params, i)
 157             peak.index += start_index
 158         block.info[params['peak info name']] = peaks
 159
 160         self._set_column(hooke=hooke, params=params,
 161                          column_name='output peak column',
 162                          values=peaks_to_mask(def_data, peaks) * def_data)
 163         outqueue.put(peaks)
 164
 165     def __setup_params(self, hooke, params):
 166         """Setup `params` from config and return the z piezo and
 167         deflection arrays.
 168         """
 169         curve = self._curve(hooke=hooke, params=params)
 170         if curve.info['experiment'] != VelocityClamp:
 171             raise Failure('%s operates on VelocityClamp experiments, not %s'
 172                           % (self.name, curve.info['experiment']))
 173         for key,value in params.items():
 174             if value == None: # Use configured default value.
 175                 params[key] = self.plugin.config[key]
 176         # TODO: convert 'see double' from nm to points
 177         name,def_unit = split_data_label(params['deflection column'])
 178         params['output peak column'] = join_data_label(
 179             params['output peak column'], def_unit)
 180         params['peak info name'] = join_data_label(
 181             params['peak info name'], def_unit)
 182         return params
 183
 184     def _peak_name(self, params, index):
 185         d_name,d_unit = split_data_label(params['deflection column'])
 186         return 'flat filter peak %d of %s' % (index, d_name)
 187
 188
 189 class FlatFilterCommand (FilterCommand):
 190     u"""Create a subset playlist of curves with enough flat peaks.
 191
 192     Notes
 193     -----
 194     This type of filter is of course very raw, and requires relatively
 195     conservative settings to safely avoid false negatives (that is, to
 196     avoid discarding interesting curves).  Using it on the protein
 197     unfolding experiments described by Sandal [#sandal2008] it has
 198     been found to reduce the data set to analyze by hand by 60-80%.
 199
 200     .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
 201       F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
 202       "Conformational equilibria in monomeric α-Synuclein at the
 203       single molecule level."
 204       PLOS Biology, 2009.
 205       doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
 206
 207     See Also
 208     --------
 209     FlatPeaksCommand : Underlying flat-based peak detection.
 210     """
 211     def __init__(self, plugin):
 212         super(FlatFilterCommand, self).__init__(
 213             plugin, name='flat filter playlist')
 214         self.arguments.extend(plugin._arguments)  # TODO: curve & block arguments
 215
 216     def filter(self, curve, hooke, inqueue, outqueue, params):
 217         params['curve'] = curve
 218         inq = Queue()
 219         outq = Queue()
 220         filt_command = [c for c in hooke.commands
 221                         if c.name=='flat filter peaks'][0]
 222         filt_command.run(hooke, inq, outq, **params)
 223         peaks = outq.get()
 224         if isinstance(peaks, UncaughtException) \
 225                 and isinstance(peaks.exception, PoorFit):
 226             return False
 227         if not (isinstance(peaks, list) and (len(peaks) == 0
 228                                              or isinstance(peaks[0], Peak))):
 229             raise Failure('Expected a list of Peaks, not %s' % peaks)
 230         ret = outq.get()
 231         if not isinstance(ret, Success):
 232             raise ret
 233         if params['min peaks'] == None: # Use configured default value.
 234             params['min peaks'] = self.plugin.config['min peaks']
 235         return len(peaks) >= int(params['min peaks'])