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

« back to all changes in this revision

Viewing changes to source/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.C

  • 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:
 
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
#include <OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h>
 
36
#include <OpenMS/ANALYSIS/OPENSWATH/DATAACCESS/DataAccessHelper.h>
 
37
 
 
38
bool SortDoubleDoublePairFirst(const std::pair<double, double>& left, const std::pair<double, double>& right)
 
39
{
 
40
  return left.first < right.first;
 
41
}
 
42
 
 
43
namespace OpenMS
 
44
{
 
45
 
 
46
  MRMFeatureFinderScoring::MRMFeatureFinderScoring() :
 
47
    DefaultParamHandler("MRMFeatureFinderScoring"),
 
48
    ProgressLogger()
 
49
  {
 
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);
 
61
 
 
62
    defaults_.insert("TransitionGroupPicker:", MRMTransitionGroupPicker().getDefaults());
 
63
 
 
64
    defaults_.insert("DIAScoring:", DIAScoring().getDefaults());
 
65
 
 
66
    defaults_.insert("EMGScoring:", EmgScoring().getDefaults());
 
67
 
 
68
    // One can turn on / off each score individually
 
69
    Param scores_to_use;
 
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);
 
89
 
 
90
    // write defaults into Param object param_
 
91
    defaultsToParam_();
 
92
 
 
93
    strict_ = true;
 
94
  }
 
95
 
 
96
  MRMFeatureFinderScoring::~MRMFeatureFinderScoring()
 
97
  {
 
98
  }
 
99
 
 
100
  void MRMFeatureFinderScoring::updateMembers_()
 
101
  {
 
102
    /*
 
103
    handle_params();
 
104
  }
 
105
 
 
106
  void MRMFeatureFinderScoring::handle_params()
 
107
  {
 
108
    */
 
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");
 
116
 
 
117
    diascoring_.setParameters(param_.copy("DIAScoring:", true));
 
118
    emgscoring_.setFitterParam(param_.copy("EmgScoring:", true));
 
119
 
 
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();
 
129
  }
 
130
 
 
131
  void MRMFeatureFinderScoring::mapExperimentToTransitionList(OpenSwath::SpectrumAccessPtr input,
 
132
                                                              TargetedExpType& transition_exp, TransitionGroupMapType& transition_group_map,
 
133
                                                              TransformationDescription trafo, double rt_extraction_window)
 
134
  {
 
135
    double rt_min, rt_max, expected_rt;
 
136
    trafo.invert();
 
137
 
 
138
    std::map<String, int> chromatogram_map;
 
139
    Size nr_chromatograms = input->getNrChromatograms();
 
140
    for (Size i = 0; i < input->getNrChromatograms(); i++)
 
141
    {
 
142
      chromatogram_map[input->getChromatogramNativeID(i)] = boost::numeric_cast<int>(i);
 
143
    }
 
144
 
 
145
    // Iterate thorugh all transitions and store the transition with the
 
146
    // corresponding chromatogram in the corresponding transition group
 
147
    Size progress = 0;
 
148
    startProgress(0, nr_chromatograms, "Mapping transitions to chromatograms ");
 
149
    for (Size i = 0; i < transition_exp.getTransitions().size(); i++)
 
150
    {
 
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())
 
154
      {
 
155
        std::cerr << "Error: Transition " + transition->getNativeID() + " from group " +
 
156
        transition->getPeptideRef() + " does not have a corresponding chromatogram" << std::endl;
 
157
        if (strict_)
 
158
        {
 
159
          throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__,
 
160
                                           "Error: Transition " + transition->getNativeID() + " from group " +
 
161
                                           transition->getPeptideRef() + " does not have a corresponding chromatogram");
 
162
        }
 
163
        continue;
 
164
      }
 
165
      MSChromatogram<ChromatogramPeak> chromatogram_old;
 
166
      OpenSwath::ChromatogramPtr cptr = input->getChromatogramById(chromatogram_map[transition->getNativeID()]);
 
167
      OpenSwathDataAccessHelper::convertToOpenMSChromatogram(chromatogram_old, cptr);
 
168
      RichPeakChromatogram chromatogram;
 
169
 
 
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
 
175
      // other way round.
 
176
      rt_max = rt_min = 0;
 
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)
 
182
      {
 
183
        if (rt_extraction_window >= 0 && (it->getRT() < rt_min || it->getRT() > rt_max))
 
184
        {
 
185
          continue;
 
186
        }
 
187
        ChromatogramPeak peak;
 
188
        peak.setMZ(it->getRT());
 
189
        peak.setIntensity(it->getIntensity());
 
190
        chromatogram.push_back(peak);
 
191
      }
 
192
      if (chromatogram.empty())
 
193
      {
 
194
        std::cerr << "Error: Could not find any points for chromatogram " + transition->getNativeID() + \
 
195
        ". Maybe your retention time transformation is off?" << std::endl;
 
196
        if (strict_)
 
197
        {
 
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?");
 
201
        }
 
202
      }
 
203
      chromatogram.setMetaValue("product_mz", transition->getProductMZ());
 
204
      chromatogram.setMetaValue("precursor_mz", transition->getPrecursorMZ());
 
205
      chromatogram.setNativeID(transition->getNativeID());
 
206
 
 
207
      // Create new transition group if there is none for this peptide
 
208
      if (transition_group_map.find(transition->getPeptideRef()) == transition_group_map.end())
 
209
      {
 
210
        MRMTransitionGroupType transition_group;
 
211
        transition_group.setTransitionGroupID(transition->getPeptideRef());
 
212
        transition_group_map[transition->getPeptideRef()] = transition_group;
 
213
      }
 
214
 
 
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());
 
219
 
 
220
      setProgress(++progress);
 
221
    }
 
222
    endProgress();
 
223
 
 
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++)
 
226
    {
 
227
      if (trgroup_it->second.getChromatograms().size() > 0 && (trgroup_it->second.getChromatograms().size() != trgroup_it->second.getTransitions().size()))
 
228
      {
 
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.");
 
232
      }
 
233
    }
 
234
  }
 
235
 
 
236
}