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