| OpenMS 2.5.0

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