Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
EmgScoring.h
Go to the documentation of this file.
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-2017.
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
43 
46 
48 
49 
50 namespace OpenMS
51 {
52 
60  class EmgScoring
61  {
62 
63  public :
64 
65  EmgScoring() { }
66 
68 
69  void setFitterParam(Param param)
70  {
72  }
73 
75  {
76  return fitter_emg1D_.getDefaults();
77  }
78 
80  template<typename SpectrumType, class TransitionT>
82  {
83  std::vector<double> fit_scores;
84  double avg_score = 0;
85  bool smooth_data = false;
86  for (Size k = 0; k < transition_group.size(); k++)
87  {
88  // get the id, then find the corresponding transition and features within this peakgroup
89  String native_id = transition_group.getChromatograms()[k].getNativeID();
90  Feature f = mrmfeature.getFeature(native_id);
91  OPENMS_PRECONDITION(f.getConvexHulls().size() == 1, "Convex hulls need to have exactly one hull point structure");
92 
93  // TODO what if score is -1 ?? e.g. if it is undefined
94  double fscore = elutionModelFit(f.getConvexHulls()[0].getHullPoints(), smooth_data);
95  fit_scores.push_back(fscore);
96  avg_score += fscore;
97  }
98 
99  avg_score /= transition_group.size();
100  return avg_score;
101  }
102 
103  // Fxn from FeatureFinderAlgorithmMRM
104  // TODO: check whether we can leave out some of the steps here, e.g. gaussian smoothing
105  double elutionModelFit(ConvexHull2D::PointArrayType current_section, bool smooth_data)
106  {
107  // We need at least 2 datapoints in order to create a fit
108  if (current_section.size() < 2)
109  {
110  return -1;
111  }
112 
113  // local PeakType is a small hack since here we *need* data of type
114  // Peak1D, otherwise our fitter will not accept it.
115  typedef Peak1D LocalPeakType;
116 
117  // -- cut line 301 of FeatureFinderAlgorithmMRM
118  std::vector<LocalPeakType> data_to_fit;
119  prepareFit_(current_section, data_to_fit, smooth_data);
120  InterpolationModel * model_rt = 0;
121  double quality = fitRT_(data_to_fit, model_rt);
122  // cut line 354 of FeatureFinderAlgorithmMRM
123  delete model_rt;
124 
125  return quality;
126 
127  }
128 
129  protected:
130  template<class LocalPeakType>
131  double fitRT_(std::vector<LocalPeakType> & rt_input_data, InterpolationModel * & model)
132  {
133  double quality;
134  //Param param;
135 
136  /*EmgFitter
137  param.setValue( "tolerance_stdev_bounding_box", tolerance_stdev_box_);
138  param.setValue( "statistics:mean", rt_stat_.mean() );
139  param.setValue( "statistics:variance", rt_stat_.variance() );
140  param.setValue( "interpolation_step", interpolation_step_rt_ );
141  param.setValue( "max_iteration", max_iteration_);
142  param.setValue( "deltaAbsError", deltaAbsError_);
143  param.setValue( "deltaRelError", deltaRelError_);
144  */
145 
146  // Set parameter for fitter
147  //fitter_emg1D.setParameters(param);
148  // Construct model for rt
149  quality = fitter_emg1D_.fit1d(rt_input_data, model);
150 
151  // Check quality
152  if (boost::math::isnan(quality)) quality = -1.0;
153  return quality;
154  }
155 
156  // Fxn from FeatureFinderAlgorithmMRM
157  // TODO: check whether we can leave out some of the steps here, e.g. gaussian smoothing
158  template<class LocalPeakType>
159  void prepareFit_(const ConvexHull2D::PointArrayType & current_section, std::vector<LocalPeakType> & data_to_fit, bool smooth_data)
160  {
161  // typedef Peak1D LocalPeakType;
162  PeakSpectrum filter_spec;
163  // first smooth the data to prevent outliers from destroying the fit
164  for (ConvexHull2D::PointArrayType::const_iterator it = current_section.begin(); it != current_section.end(); it++)
165  {
166  LocalPeakType p;
167  p.setMZ(it->getX());
168  p.setIntensity(it->getY());
169  filter_spec.push_back(p);
170  }
171 
172  // add two peaks at the beginning and at the end for better fit
173  // therefore calculate average distance first
174  std::vector<double> distances;
175  for (Size j = 1; j < filter_spec.size(); ++j)
176  {
177  distances.push_back(filter_spec[j].getMZ() - filter_spec[j - 1].getMZ());
178  }
179  double dist_average = std::accumulate(distances.begin(), distances.end(), 0.0) / (double) distances.size();
180 
181  // append peaks
182  Peak1D new_peak;
183  new_peak.setIntensity(0);
184  new_peak.setMZ(filter_spec.back().getMZ() + dist_average);
185  filter_spec.push_back(new_peak);
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 
191  // prepend peaks
192  new_peak.setMZ(filter_spec.front().getMZ() - dist_average);
193  filter_spec.insert(filter_spec.begin(), new_peak);
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 
199  // To get an estimate of the peak quality, we probably should not smooth
200  // and/or transform the data.
201  if (smooth_data)
202  {
203  GaussFilter filter;
204  Param filter_param(filter.getParameters());
205  filter.setParameters(filter_param);
206  filter_param.setValue("gaussian_width", 4 * dist_average);
207  filter.setParameters(filter_param);
208  filter.filter(filter_spec);
209  }
210 
211  // transform the data for fitting and fit RT profile
212  for (Size j = 0; j != filter_spec.size(); ++j)
213  {
214  LocalPeakType p;
215  p.setPosition(filter_spec[j].getMZ());
216  p.setIntensity(filter_spec[j].getIntensity());
217  data_to_fit.push_back(p);
218  }
219  }
220 
222  };
223 }
224 
225 #endif /* EMGSCORING_H_ */
const double k
A more convenient string class.
Definition: String.h:57
Scoring of an elution peak using an exponentially modified gaussian distribution model.
Definition: EmgScoring.h:60
#define OPENMS_PRECONDITION(condition, message)
Precondition macro.
Definition: openms/include/OpenMS/CONCEPT/Macros.h:107
Size size() const
Definition: MRMTransitionGroup.h:126
std::vector< PointType > PointArrayType
Definition: ConvexHull2D.h:77
EmgFitter1D fitter_emg1D_
Definition: EmgScoring.h:221
Exponentially modified gaussian distribution fitter (1-dim.) using Levenberg-Marquardt algorithm (Eig...
Definition: EmgFitter1D.h:47
QualityType fit1d(const RawDataArrayType &range, InterpolationModel *&model)
return interpolation model
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
const std::vector< ConvexHull2D > & getConvexHulls() const
Non-mutable access to the convex hulls.
void setMZ(CoordinateType mz)
Mutable access to m/z.
Definition: Peak1D.h:120
void setIntensity(IntensityType intensity)
Mutable access to the data point intensity (height)
Definition: Peak1D.h:111
void setParameters(const Param &param)
Sets the parameters.
The representation of a group of transitions in a targeted proteomics experiment. ...
Definition: MRMTransitionGroup.h:68
Abstract class for 1D-models that are approximated using linear interpolation.
Definition: InterpolationModel.h:55
const Param & getDefaults() const
Non-mutable access to the default parameters.
The representation of a 1D spectrum.
Definition: MSSpectrum.h:67
const std::vector< ChromatogramType > & getChromatograms() const
Definition: MRMTransitionGroup.h:175
A 1-dimensional raw data point or peak.
Definition: Peak1D.h:55
void setFitterParam(Param param)
Definition: EmgScoring.h:69
An LC-MS feature.
Definition: Feature.h:70
const Param & getParameters() const
Non-mutable access to the parameters.
~EmgScoring()
Definition: EmgScoring.h:67
double elutionModelFit(ConvexHull2D::PointArrayType current_section, bool smooth_data)
Definition: EmgScoring.h:105
Management and storage of parameters / INI files.
Definition: Param.h:75
This class represents a Gaussian lowpass-filter which works on uniform as well as on non-uniform prof...
Definition: GaussFilter.h:75
Feature & getFeature(String key)
get a specified feature
double calcElutionFitScore(MRMFeature &mrmfeature, MRMTransitionGroup< SpectrumType, TransitionT > &transition_group)
calculate the elution profile fit score
Definition: EmgScoring.h:81
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:128
void prepareFit_(const ConvexHull2D::PointArrayType &current_section, std::vector< LocalPeakType > &data_to_fit, bool smooth_data)
Definition: EmgScoring.h:159
EmgScoring()
Definition: EmgScoring.h:65
void filter(MSSpectrum &spectrum)
Smoothes an MSSpectrum containing profile data.
Definition: GaussFilter.h:93
A multi-chromatogram MRM feature.
Definition: MRMFeature.h:50
Param getDefaults()
Definition: EmgScoring.h:74
double fitRT_(std::vector< LocalPeakType > &rt_input_data, InterpolationModel *&model)
Definition: EmgScoring.h:131

OpenMS / TOPP release 2.3.0 Documentation generated on Tue Jan 9 2018 18:22:00 using doxygen 1.8.13