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.
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.
30
// --------------------------------------------------------------------------
31
// $Maintainer: Hannes Roest $
32
// $Authors: Hannes Roest $
33
// --------------------------------------------------------------------------
35
#include <OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h>
36
#include <OpenMS/ANALYSIS/OPENSWATH/DATAACCESS/DataAccessHelper.h>
38
bool SortDoubleDoublePairFirst(const std::pair<double, double>& left, const std::pair<double, double>& right)
40
return left.first < right.first;
46
MRMFeatureFinderScoring::MRMFeatureFinderScoring() :
47
DefaultParamHandler("MRMFeatureFinderScoring"),
50
defaults_.setValue("stop_report_after_feature", -1, "Stop reporting after feature (ordered by quality; -1 means do not stop).");
51
defaults_.setValue("rt_extraction_window", -1.0, "Only extract RT around this value (-1 means extract over the whole range, a value of 500 means to extract around +/- 500 s of the expected elution). For this to work, the TraML input file needs to contain normalized RT values.");
52
defaults_.setValue("rt_normalization_factor", 1.0, "The normalized RT is expected to be between 0 and 1. If your normalized RT has a different range, pass this here (e.g. it goes from 0 to 100, set this value to 100)");
53
defaults_.setValue("quantification_cutoff", 0.0, "Cutoff below which peaks should not be used for quantification any more", StringList::create("advanced"));
54
defaults_.setMinFloat("quantification_cutoff", 0.0);
55
defaults_.setValue("write_convex_hull", "false", "Whether to write out all points of all features into the featureXML", StringList::create("advanced"));
56
defaults_.setValidStrings("write_convex_hull", StringList::create("true,false"));
57
defaults_.setValue("add_up_spectra", 1, "Add up spectra around the peak apex (needs to be a non-even integer)", StringList::create("advanced"));
58
defaults_.setMinInt("add_up_spectra", 1);
59
defaults_.setValue("spacing_for_spectra_resampling", 0.005, "If spectra are to be added, use this spacing to add them up", StringList::create("advanced"));
60
defaults_.setMinFloat("spacing_for_spectra_resampling", 0.0);
62
defaults_.insert("TransitionGroupPicker:", MRMTransitionGroupPicker().getDefaults());
64
defaults_.insert("DIAScoring:", DIAScoring().getDefaults());
66
defaults_.insert("EMGScoring:", EmgScoring().getDefaults());
68
// One can turn on / off each score individually
70
scores_to_use.setValue("use_shape_score", "true", "Use the shape score", StringList::create("advanced"));
71
scores_to_use.setValidStrings("use_shape_score", StringList::create("true,false"));
72
scores_to_use.setValue("use_coelution_score", "true", "Use the coelution score", StringList::create("advanced"));
73
scores_to_use.setValidStrings("use_coelution_score", StringList::create("true,false"));
74
scores_to_use.setValue("use_rt_score", "true", "Use the retention time score", StringList::create("advanced"));
75
scores_to_use.setValidStrings("use_rt_score", StringList::create("true,false"));
76
scores_to_use.setValue("use_library_score", "true", "Use the library score", StringList::create("advanced"));
77
scores_to_use.setValidStrings("use_library_score", StringList::create("true,false"));
78
scores_to_use.setValue("use_elution_model_score", "true", "Use the elution model (EMG) score", StringList::create("advanced"));
79
scores_to_use.setValidStrings("use_elution_model_score", StringList::create("true,false"));
80
scores_to_use.setValue("use_intensity_score", "true", "Use the intensity score", StringList::create("advanced"));
81
scores_to_use.setValidStrings("use_intensity_score", StringList::create("true,false"));
82
scores_to_use.setValue("use_nr_peaks_score", "true", "Use the number of peaks score", StringList::create("advanced"));
83
scores_to_use.setValidStrings("use_nr_peaks_score", StringList::create("true,false"));
84
scores_to_use.setValue("use_total_xic_score", "true", "Use the total XIC score", StringList::create("advanced"));
85
scores_to_use.setValidStrings("use_total_xic_score", StringList::create("true,false"));
86
scores_to_use.setValue("use_sn_score", "true", "Use the SN (signal to noise) score", StringList::create("advanced"));
87
scores_to_use.setValidStrings("use_sn_score", StringList::create("true,false"));
88
defaults_.insert("Scores:", scores_to_use);
90
// write defaults into Param object param_
96
MRMFeatureFinderScoring::~MRMFeatureFinderScoring()
100
void MRMFeatureFinderScoring::updateMembers_()
106
void MRMFeatureFinderScoring::handle_params()
109
stop_report_after_feature_ = (int)param_.getValue("stop_report_after_feature");
110
rt_extraction_window_ = (DoubleReal)param_.getValue("rt_extraction_window");
111
rt_normalization_factor_ = (DoubleReal)param_.getValue("rt_normalization_factor");
112
quantification_cutoff_ = (DoubleReal)param_.getValue("quantification_cutoff");
113
write_convex_hull_ = param_.getValue("write_convex_hull").toBool();
114
add_up_spectra_ = param_.getValue("add_up_spectra");
115
spacing_for_spectra_resampling_ = param_.getValue("spacing_for_spectra_resampling");
117
diascoring_.setParameters(param_.copy("DIAScoring:", true));
118
emgscoring_.setFitterParam(param_.copy("EmgScoring:", true));
120
use_coelution_score_ = param_.getValue("Scores:use_coelution_score").toBool();
121
use_shape_score_ = param_.getValue("Scores:use_shape_score").toBool();
122
use_rt_score_ = param_.getValue("Scores:use_rt_score").toBool();
123
use_library_score_ = param_.getValue("Scores:use_library_score").toBool();
124
use_elution_model_score_ = param_.getValue("Scores:use_elution_model_score").toBool();
125
use_intensity_score_ = param_.getValue("Scores:use_intensity_score").toBool();
126
use_total_xic_score_ = param_.getValue("Scores:use_total_xic_score").toBool();
127
use_nr_peaks_score_ = param_.getValue("Scores:use_nr_peaks_score").toBool();
128
use_sn_score_ = param_.getValue("Scores:use_sn_score").toBool();
131
void MRMFeatureFinderScoring::mapExperimentToTransitionList(OpenSwath::SpectrumAccessPtr input,
132
TargetedExpType& transition_exp, TransitionGroupMapType& transition_group_map,
133
TransformationDescription trafo, double rt_extraction_window)
135
double rt_min, rt_max, expected_rt;
138
std::map<String, int> chromatogram_map;
139
Size nr_chromatograms = input->getNrChromatograms();
140
for (Size i = 0; i < input->getNrChromatograms(); i++)
142
chromatogram_map[input->getChromatogramNativeID(i)] = boost::numeric_cast<int>(i);
145
// Iterate thorugh all transitions and store the transition with the
146
// corresponding chromatogram in the corresponding transition group
148
startProgress(0, nr_chromatograms, "Mapping transitions to chromatograms ");
149
for (Size i = 0; i < transition_exp.getTransitions().size(); i++)
151
// get the current transition and try to find the corresponding chromatogram
152
const TransitionType* transition = &transition_exp.getTransitions()[i];
153
if (chromatogram_map.find(transition->getNativeID()) == chromatogram_map.end())
155
std::cerr << "Error: Transition " + transition->getNativeID() + " from group " +
156
transition->getPeptideRef() + " does not have a corresponding chromatogram" << std::endl;
159
throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__,
160
"Error: Transition " + transition->getNativeID() + " from group " +
161
transition->getPeptideRef() + " does not have a corresponding chromatogram");
165
MSChromatogram<ChromatogramPeak> chromatogram_old;
166
OpenSwath::ChromatogramPtr cptr = input->getChromatogramById(chromatogram_map[transition->getNativeID()]);
167
OpenSwathDataAccessHelper::convertToOpenMSChromatogram(chromatogram_old, cptr);
168
RichPeakChromatogram chromatogram;
170
// Create the chromatogram information
171
// Get the expected retention time, apply the RT-transformation
172
// (which describes the normalization) and then take the difference.
173
// Note that we inverted the transformation in the beginning because
174
// we want to transform from normalized to real RTs here and not the
177
expected_rt = PeptideRTMap_[transition->getPeptideRef()];
178
double de_normalized_experimental_rt = trafo.apply(expected_rt);
179
rt_max = de_normalized_experimental_rt + rt_extraction_window;
180
rt_min = de_normalized_experimental_rt - rt_extraction_window;
181
for (MSChromatogram<ChromatogramPeak>::const_iterator it = chromatogram_old.begin(); it != chromatogram_old.end(); ++it)
183
if (rt_extraction_window >= 0 && (it->getRT() < rt_min || it->getRT() > rt_max))
187
ChromatogramPeak peak;
188
peak.setMZ(it->getRT());
189
peak.setIntensity(it->getIntensity());
190
chromatogram.push_back(peak);
192
if (chromatogram.empty())
194
std::cerr << "Error: Could not find any points for chromatogram " + transition->getNativeID() + \
195
". Maybe your retention time transformation is off?" << std::endl;
198
throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__,
199
"Error: Could not find any points for chromatogram " + transition->getNativeID() + \
200
". Maybe your retention time transformation is off?");
203
chromatogram.setMetaValue("product_mz", transition->getProductMZ());
204
chromatogram.setMetaValue("precursor_mz", transition->getPrecursorMZ());
205
chromatogram.setNativeID(transition->getNativeID());
207
// Create new transition group if there is none for this peptide
208
if (transition_group_map.find(transition->getPeptideRef()) == transition_group_map.end())
210
MRMTransitionGroupType transition_group;
211
transition_group.setTransitionGroupID(transition->getPeptideRef());
212
transition_group_map[transition->getPeptideRef()] = transition_group;
215
// Now add the transition and the chromatogram to the group
216
MRMTransitionGroupType& transition_group = transition_group_map[transition->getPeptideRef()];
217
transition_group.addTransition(*transition, transition->getNativeID());
218
transition_group.addChromatogram(chromatogram, chromatogram.getNativeID());
220
setProgress(++progress);
224
// The assumption is that for each transition that is in the TargetedExperiment we have exactly one chromatogram
225
for (TransitionGroupMapType::iterator trgroup_it = transition_group_map.begin(); trgroup_it != transition_group_map.end(); trgroup_it++)
227
if (trgroup_it->second.getChromatograms().size() > 0 && (trgroup_it->second.getChromatograms().size() != trgroup_it->second.getTransitions().size()))
229
throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Error: Could not match all transition to all chromatograms:\nFor chromatogram " + \
230
trgroup_it->second.getTransitionGroupID() + " I found " + String(trgroup_it->second.getChromatograms().size()) + \
231
" chromatograms but " + String(trgroup_it->second.getTransitions().size()) + " transitions.");