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 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         d_name,d_unit = split_data_label(params['input deflection column'])
 163         for i,peak in enumerate(peaks):
 164             peak.name = 'flat filter peak %d of %s' % (i, d_name)
 165             peak.index += start_index
 166         data.info[params['peak info name']] = peaks
 167
 168         # Add a peak-masked deflection column.
 169         new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
 170         new.info = copy.deepcopy(data.info)
 171         new[:,:-1] = data
 172         d_name,d_unit = split_data_label(params['input deflection column'])
 173         new.info['columns'].append(
 174             join_data_label(params['output peak column'], d_unit))
 175         new[:,-1] = peaks_to_mask(d_data, peaks) * d_data
 176         outqueue.put(peaks)
 177         params['curve'].data[block_index] = new
 178
 179     def _setup(self, hooke, params):
 180         """Setup `params` from config and return the z piezo and
 181         deflection arrays.
 182         """
 183         curve = params['curve']
 184         if curve.info['experiment'] != VelocityClamp:
 185             raise Failure('%s operates on VelocityClamp experiments, not %s'
 186                           % (self.name, curve.info['experiment']))
 187         data = None
 188         for i,block in enumerate(curve.data):
 189             if block.info['name'].startswith(params['retract block']):
 190                 data = block
 191                 block_index = i
 192                 break
 193         if data == None:
 194             raise Failure('No retraction blocks in %s.' % curve)
 195         z_data = data[:,data.info['columns'].index(
 196                 params['input distance column'])]
 197         d_data = data[:,data.info['columns'].index(
 198                 params['input deflection column'])]
 199         for key,value in params.items():
 200             if value == None: # Use configured default value.
 201                 params[key] = self.plugin.config[key]
 202         # TODO: convert 'see double' from nm to points
 203         return (data, block_index, z_data, d_data, params)
 204
 205
 206 class FlatFilterCommand (FilterCommand):
 207     u"""Create a subset playlist of curves with enough flat peaks.
 208
 209     Notes
 210     -----
 211     This type of filter is of course very raw, and requires relatively
 212     conservative settings to safely avoid false negatives (that is, to
 213     avoid discarding interesting curves).  Using it on the protein
 214     unfolding experiments described by Sandal [#sandal2008] it has
 215     been found to reduce the data set to analyze by hand by 60-80%.
 216
 217     .. [#sandal2008] M. Sandal, F. Valle, I. Tessari, S. Mammi, E. Bergantino,
 218       F. Musiani, M. Brucale, L. Bubacco, B. Samorì.
 219       "Conformational equilibria in monomeric α-Synuclein at the
 220       single molecule level."
 221       PLOS Biology, 2009.
 222       doi: `10.1371/journal.pbio.0060006 <http://dx.doi.org/10.1371/journal.pbio.0060006>`_
 223
 224     See Also
 225     --------
 226     FlatPeaksCommand : Underlying flat-based peak detection.
 227     """
 228     def __init__(self, plugin):
 229         super(FlatFilterCommand, self).__init__(
 230             plugin, name='flat filter playlist')
 231         self.arguments.extend(plugin._arguments)
 232
 233     def filter(self, curve, hooke, inqueue, outqueue, params):
 234         params['curve'] = curve
 235         inq = Queue()
 236         outq = Queue()
 237         filt_command = [c for c in hooke.commands
 238                         if c.name=='flat filter peaks'][0]
 239         filt_command.run(hooke, inq, outq, **params)
 240         peaks = outq.get()
 241         if isinstance(peaks, UncaughtException) \
 242                 and isinstance(peaks.exception, PoorFit):
 243             return False
 244         if not (isinstance(peaks, list) and (len(peaks) == 0
 245                                              or isinstance(peaks[0], Peak))):
 246             raise Failure('Expected a list of Peaks, not %s' % peaks)
 247         ret = outq.get()
 248         if not isinstance(ret, Success):
 249             raise ret
 250         if params['min peaks'] == None: # Use configured default value.
 251             params['min peaks'] = self.plugin.config['min peaks']
 252         return len(peaks) >= int(params['min peaks'])