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