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 ..util.fit import PoorFit
  46 from ..util.peak import (find_peaks, peaks_to_mask,
  47                          find_peaks_arguments, Peak, _kwargs)
  48 from ..util.si import join_data_label, split_data_label
  49 from . import Plugin, argument_to_setting
  50 from .curve import ColumnAddingCommand
  51 from .playlist import FilterCommand
  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)]
  88         # append FlatFilterCommand so it can steal arguments from
  89         # FlatPeaksCommand.
  90         self._commands.append(FlatFilterCommand(self))
  91
  92     def dependencies(self):
  93         return ['vclamp']
  94
  95     def default_settings(self):
  96         return self._settings
  97
  98
  99 class FlatPeaksCommand (ColumnAddingCommand):
 100     """Detect peaks in velocity clamp data using noise statistics.
 101
 102     Notes
 103     -----
 104     Noise analysis on the retraction curve:
 105
 106     1) A median window filter (using
 107       :func:`scipy.signal.signaltools.medfilt`) smooths the
 108       deflection.
 109     2) The deflection derivative is calculated (using
 110       :func:`numpy.diff` which uses forward differencing).
 111     3) Peaks in the derivative curve are extracted with
 112       :func:`~hooke.plugins.peak.find_peaks`.
 113
 114     The algorithm was originally Francesco Musiani's idea.
 115     """
 116     def __init__(self, plugin):
 117         plugin_arguments = [a for a in plugin._arguments
 118                             if a.name != 'min peaks']
 119         # Drop min peaks, since we're not filtering with this
 120         # function, just detecting peaks.
 121         super(FlatPeaksCommand, self).__init__(
 122             name='flat filter peaks',
 123             columns=[
 124                 ('distance column', 'surface distance (m)', """
 125 Name of the column to use as the surface position input.
 126 """.strip()),
 127                 ('deflection column', 'surface deflection (m)', """
 128 Name of the column to use as the deflection input.
 129 """.strip()),
 130                 ],
 131             new_columns=[
 132                 ('output peak column', 'flat filter peaks', """
 133 Name of the column (without units) to use as the peak output.
 134 """.split()),
 135                 ],
 136             arguments=[
 137                 Argument(name='peak info name', type='string',
 138                          default='flat filter peaks',
 139                          help="""
 140 Name (without units) for storing the list of peaks in the `.info`
 141 dictionary.
 142 """.strip()),
 143                 ] + plugin_arguments,
 144             help=self.__doc__, plugin=plugin)
 145
 146     def _run(self, hooke, inqueue, outqueue, params):
 147         self._add_to_command_stack(params)
 148         params = self.__setup_params(hooke=hooke, params=params)
 149         block = self._block(hooke=hooke, params=params)
 150         dist_data = self._get_column(hooke=hooke, params=params,
 151                                      column_name='distance column')
 152         def_data = self._get_column(hooke=hooke, params=params,
 153                                      column_name='deflection column')
 154         start_index = absolute(dist_data-params['blind window']).argmin()
 155         median = medfilt(def_data[start_index:], params['median window'])
 156         deriv = diff(median)
 157         peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments,
 158                                             argument_input_keys=True))
 159         for i,peak in enumerate(peaks):
 160             peak.name = self._peak_name(params, i)
 161             peak.index += start_index
 162         block.info[params['peak info name']] = peaks
 163
 164         self._set_column(hooke=hooke, params=params,
 165                          column_name='output peak column',
 166                          values=peaks_to_mask(def_data, peaks) * def_data)
 167         outqueue.put(peaks)
 168
 169     def __setup_params(self, hooke, params):
 170         """Setup `params` from config and return the z piezo and
 171         deflection arrays.
 172         """
 173         curve = self._curve(hooke=hooke, params=params)
 174         if curve.info['experiment'] != VelocityClamp:
 175             raise Failure('%s operates on VelocityClamp experiments, not %s'
 176                           % (self.name, curve.info['experiment']))
 177         for key,value in params.items():
 178             if value == None and key in self.plugin.config:
 179                 # Use configured default value.
 180                 params[key] = self.plugin.config[key]
 181         # TODO: convert 'see double' from nm to points
 182         name,def_unit = split_data_label(params['deflection column'])
 183         params['output peak column'] = join_data_label(
 184             params['output peak column'], def_unit)
 185         params['peak info name'] = join_data_label(
 186             params['peak info name'], def_unit)
 187         return params
 188
 189     def _peak_name(self, params, index):
 190         d_name,d_unit = split_data_label(params['deflection column'])
 191         return 'flat filter peak %d of %s' % (index, d_name)
 192
 193
 194 class FlatFilterCommand (FilterCommand):
 195     u"""Create a subset playlist of curves with enough flat peaks.
 196
 197     Notes
 198     -----
 199     This type of filter is of course very raw, and requires relatively
 200     conservative settings to safely avoid false negatives (that is, to
 201     avoid discarding interesting curves).  Using it on the protein
 202     unfolding experiments described by Sandal [#sandal2008] it has
 203     been found to reduce the data set to analyze by hand by 60-80%.
 204
 205     .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
 206       F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
 207       "Conformational equilibria in monomeric α-Synuclein at the
 208       single molecule level."
 209       PLOS Biology, 2009.
 210       doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
 211
 212     See Also
 213     --------
 214     FlatPeaksCommand : Underlying flat-based peak detection.
 215     """
 216     def __init__(self, plugin):
 217         flat_peaks = [c for c in plugin._commands
 218                       if c.name == 'flat filter peaks'][0]
 219         flat_peaks_arg_names = [a.name for a in flat_peaks.arguments]
 220         plugin_arguments = [a for a in plugin._arguments
 221                             if a.name not in flat_peaks_arg_names]
 222         arguments = flat_peaks.arguments + plugin_arguments
 223         super(FlatFilterCommand, self).__init__(
 224             plugin, name='flat filter playlist')
 225         self.arguments.extend(arguments)
 226
 227     def filter(self, curve, hooke, inqueue, outqueue, params):
 228         params['curve'] = curve
 229         inq = Queue()
 230         outq = Queue()
 231         filt_command = [c for c in hooke.commands
 232                         if c.name=='flat filter peaks'][0]
 233         filt_command.run(hooke, inq, outq, **params)
 234         peaks = outq.get()
 235         if isinstance(peaks, UncaughtException) \
 236                 and isinstance(peaks.exception, PoorFit):
 237             return False
 238         if not (isinstance(peaks, list) and (len(peaks) == 0
 239                                              or isinstance(peaks[0], Peak))):
 240             raise Failure('Expected a list of Peaks, not %s: %s'
 241                           % (type(peaks), peaks))
 242         ret = outq.get()
 243         if not isinstance(ret, Success):
 244             raise ret
 245         if params['min peaks'] == None: # Use configured default value.
 246             params['min peaks'] = self.plugin.config['min peaks']
 247         return len(peaks) >= int(params['min peaks'])