← Back to branch summary

~ubuntu-branches/ubuntu/wily/openms/wily

« back to all changes in this revision

Viewing changes to include/OpenMS/TRANSFORMATIONS/FEATUREFINDER/EmgScoring.h

  • Committer: Package Import Robot
  • Author(s): Filippo Rusconi
  • Date: 2013-12-20 11:30:16 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20131220113016-wre5g9bteeheq6he
Tags: 1.11.1-3
* remove version number from libbost development package names;
* ensure that AUTHORS is correctly shipped in all packages.

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/OpenMS/TRANSFORMATIONS/FEATUREFINDER/EmgScoring.h
 
1
// --------------------------------------------------------------------------
 
2
//                   OpenMS -- Open-Source Mass Spectrometry               
 
3
// --------------------------------------------------------------------------
 
4
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
 
5
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
 
6
// 
 
7
// This software is released under a three-clause BSD license:
 
8
//  * Redistributions of source code must retain the above copyright
 
9
//    notice, this list of conditions and the following disclaimer.
 
10
//  * Redistributions in binary form must reproduce the above copyright
 
11
//    notice, this list of conditions and the following disclaimer in the
 
12
//    documentation and/or other materials provided with the distribution.
 
13
//  * Neither the name of any author or any participating institution 
 
14
//    may be used to endorse or promote products derived from this software 
 
15
//    without specific prior written permission.
 
16
// For a full list of authors, refer to the file AUTHORS. 
 
17
// --------------------------------------------------------------------------
 
18
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 
19
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 
20
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 
21
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING 
 
22
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
 
23
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
 
24
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
 
25
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
 
26
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
 
27
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF 
 
28
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
29
// 
 
30
// --------------------------------------------------------------------------
 
31
// $Maintainer: Hannes Roest $
 
32
// $Authors: Hannes Roest $
 
33
// --------------------------------------------------------------------------
 
34
 
 
35
#ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_EMGSCORING_H
 
36
#define OPENMS_TRANSFORMATIONS_FEATUREFINDER_EMGSCORING_H
 
37
 
 
38
#include <vector>
 
39
#include <boost/math/special_functions/fpclassify.hpp> // for isnan
 
40
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/EmgFitter1D.h>
 
41
#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/EmgModel.h>
 
42
#include <OpenMS/FILTERING/SMOOTHING/GaussFilter.h>
 
43
 
 
44
#include <OpenMS/KERNEL/MRMFeature.h>
 
45
#include <OpenMS/KERNEL/MRMTransitionGroup.h>
 
46
 
 
47
#include <OpenMS/KERNEL/StandardTypes.h>
 
48
 
 
49
 
 
50
namespace OpenMS
 
51
{
 
52
 
 
53
  /**
 
54
    @brief Scoring of an elution peak using an exponentially modified gaussian distribution model.
 
55
    
 
56
    This class uses the original ideas from FeatureFinderAlgorithmMRM to
 
57
    construct an interface that allows scoring of chromatographic peaks.
 
58
 
 
59
  */
 
60
  class EmgScoring
 
61
  {
 
62
 
 
63
  public :
 
64
 
 
65
    EmgScoring() { }
 
66
 
 
67
    ~EmgScoring() { }
 
68
 
 
69
    void setFitterParam(Param param)
 
70
    {
 
71
      fitter_emg1D_.setParameters(param);
 
72
    }
 
73
 
 
74
    Param getDefaults()
 
75
    {
 
76
      return fitter_emg1D_.getDefaults();
 
77
    }
 
78
 
 
79
    /// calculate the elution profile fit score
 
80
    template<typename SpectrumType, class TransitionT>
 
81
    double calcElutionFitScore(MRMFeature & mrmfeature, MRMTransitionGroup<SpectrumType, TransitionT> & transition_group)
 
82
    {
 
83
 
 
84
      std::vector<double> fit_scores;
 
85
      double avg_score = 0;
 
86
      bool smooth_data = false;
 
87
      for (Size k = 0; k < transition_group.size(); k++)
 
88
      {
 
89
        // get the id, then find the corresponding transition and features within this peakgroup
 
90
        String native_id = transition_group.getChromatograms()[k].getNativeID();
 
91
        Feature f = mrmfeature.getFeature(native_id);
 
92
        OPENMS_PRECONDITION(f.getConvexHulls().size() == 1, "Convex hulls need to have exactly one hull point structure");
 
93
 
 
94
        // TODO what if score is -1 ?? e.g. if it is undefined
 
95
        double fscore = elutionModelFit(f.getConvexHulls()[0].getHullPoints(), smooth_data);
 
96
        fit_scores.push_back(fscore);
 
97
        avg_score += fscore;
 
98
      }
 
99
 
 
100
      avg_score /= transition_group.size();
 
101
      return avg_score;
 
102
    }
 
103
 
 
104
    // Fxn from FeatureFinderAlgorithmMRM
 
105
    // TODO: check whether we can leave out some of the steps here, e.g. gaussian smoothing
 
106
    double elutionModelFit(ConvexHull2D::PointArrayType current_section, bool smooth_data)
 
107
    {
 
108
      // We need at least 2 datapoints in order to create a fit
 
109
      if (current_section.size() < 2)
 
110
      {
 
111
        return -1;
 
112
      }
 
113
 
 
114
      // local PeakType is a small hack since here we *need* data of type
 
115
      // Peak1D, otherwise our fitter will not accept it.
 
116
      typedef Peak1D LocalPeakType;
 
117
 
 
118
 
 
119
      // -- cut line 301 of FeatureFinderAlgorithmMRM
 
120
      std::vector<LocalPeakType> data_to_fit;
 
121
      prepareFit_(current_section, data_to_fit, smooth_data);
 
122
      InterpolationModel * model_rt = 0;
 
123
      DoubleReal quality = fitRT_(data_to_fit, model_rt);
 
124
      // cut line 354 of FeatureFinderAlgorithmMRM
 
125
      delete model_rt;
 
126
 
 
127
      return quality;
 
128
 
 
129
    }
 
130
 
 
131
  protected:
 
132
    template<class LocalPeakType>
 
133
    double fitRT_(std::vector<LocalPeakType> & rt_input_data, InterpolationModel * & model)
 
134
    {
 
135
      DoubleReal quality;
 
136
      //Param param;
 
137
 
 
138
      /*EmgFitter
 
139
       param.setValue( "tolerance_stdev_bounding_box", tolerance_stdev_box_);
 
140
       param.setValue( "statistics:mean", rt_stat_.mean() );
 
141
       param.setValue( "statistics:variance", rt_stat_.variance() );
 
142
       param.setValue( "interpolation_step", interpolation_step_rt_ );
 
143
       param.setValue( "max_iteration", max_iteration_);
 
144
       param.setValue( "deltaAbsError", deltaAbsError_);
 
145
       param.setValue( "deltaRelError", deltaRelError_);
 
146
       */
 
147
 
 
148
      // Set parameter for fitter
 
149
      //fitter_emg1D.setParameters(param);
 
150
      // Construct model for rt
 
151
      quality = fitter_emg1D_.fit1d(rt_input_data, model);
 
152
 
 
153
      // Check quality
 
154
      if (boost::math::isnan(quality)) quality = -1.0;
 
155
      return quality;
 
156
    }
 
157
 
 
158
    // Fxn from FeatureFinderAlgorithmMRM
 
159
    // TODO: check whether we can leave out some of the steps here, e.g. gaussian smoothing
 
160
    template<class LocalPeakType>
 
161
    void prepareFit_(const ConvexHull2D::PointArrayType & current_section, std::vector<LocalPeakType> & data_to_fit, bool smooth_data)
 
162
    {
 
163
      // typedef Peak1D LocalPeakType;
 
164
      PeakSpectrum filter_spec;
 
165
      // first smooth the data to prevent outliers from destroying the fit
 
166
      for (ConvexHull2D::PointArrayType::const_iterator it = current_section.begin(); it != current_section.end(); it++)
 
167
      {
 
168
        LocalPeakType p;
 
169
        p.setMZ(it->getX());
 
170
        p.setIntensity(it->getY());
 
171
        filter_spec.push_back(p);
 
172
      }
 
173
 
 
174
      // add two peaks at the beginning and at the end for better fit
 
175
      // therefore calculate average distance first
 
176
      std::vector<DoubleReal> distances;
 
177
      for (Size j = 1; j < filter_spec.size(); ++j)
 
178
      {
 
179
        distances.push_back(filter_spec[j].getMZ() - filter_spec[j - 1].getMZ());
 
180
      }
 
181
      DoubleReal dist_average = std::accumulate(distances.begin(), distances.end(), 0.0) / (DoubleReal) distances.size();
 
182
 
 
183
      // append peaks
 
184
      Peak1D new_peak;
 
185
      new_peak.setIntensity(0);
 
186
      new_peak.setMZ(filter_spec.back().getMZ() + dist_average);
 
187
      filter_spec.push_back(new_peak);
 
188
      new_peak.setMZ(filter_spec.back().getMZ() + dist_average);
 
189
      filter_spec.push_back(new_peak);
 
190
      new_peak.setMZ(filter_spec.back().getMZ() + dist_average);
 
191
      filter_spec.push_back(new_peak);
 
192
 
 
193
      // prepend peaks
 
194
      new_peak.setMZ(filter_spec.front().getMZ() - dist_average);
 
195
      filter_spec.insert(filter_spec.begin(), new_peak);
 
196
      new_peak.setMZ(filter_spec.front().getMZ() - dist_average);
 
197
      filter_spec.insert(filter_spec.begin(), new_peak);
 
198
      new_peak.setMZ(filter_spec.front().getMZ() - dist_average);
 
199
      filter_spec.insert(filter_spec.begin(), new_peak);
 
200
 
 
201
      // To get an estimate of the peak quality, we probably should not smooth
 
202
      // and/or transform the data.
 
203
      if (smooth_data)
 
204
      {
 
205
        GaussFilter filter;
 
206
        Param filter_param(filter.getParameters());
 
207
        filter.setParameters(filter_param);
 
208
        filter_param.setValue("gaussian_width", 4 * dist_average);
 
209
        filter.setParameters(filter_param);
 
210
        filter.filter(filter_spec);
 
211
      }
 
212
 
 
213
      // transform the data for fitting and fit RT profile
 
214
      for (Size j = 0; j != filter_spec.size(); ++j)
 
215
      {
 
216
        LocalPeakType p;
 
217
        p.setPosition(filter_spec[j].getMZ());
 
218
        p.setIntensity(filter_spec[j].getIntensity());
 
219
        data_to_fit.push_back(p);
 
220
      }
 
221
    }
 
222
 
 
223
    EmgFitter1D fitter_emg1D_;
 
224
  };
 
225
}
 
226
 
 
227
#endif /* EMGSCORING_H_ */