Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
MRMFeatureFinderScoring.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-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_ANALYSIS_OPENSWATH_MRMFEATUREFINDERSCORING_H
36 #define OPENMS_ANALYSIS_OPENSWATH_MRMFEATUREFINDERSCORING_H
37 
38 #define run_identifier "unique_run_identifier"
39 #define USE_SP_INTERFACE
40 
41 // move to TOPPTool
42 #include <OpenMS/FORMAT/MzMLFile.h>
46 
51 
52 // peak picking & noise estimation
55 
56 // data access
62 
63 // scoring
68 
69 // auxiliary
71 
72 #ifdef _OPENMP
73 #include <omp.h>
74 #endif
75 
76 bool SortDoubleDoublePairFirst(const std::pair<double, double>& left, const std::pair<double, double>& right);
77 
78 namespace OpenMS
79 {
80 
86  {
88  double library_corr;
89  double library_rmsd;
90  double norm_rt_score;
93  double massdev_score;
96  double yseries_score;
97  double log_sn_score;
98 
100  double xcorr_shape_score, double log_sn_score)
101  {
102  // some scores based on manual evaluation of 80 chromatograms
103  // quick LDA average model on 100 2xCrossvalidated runs (0.85 TPR/0.17 FDR)
104  // true: mean 4.2 with sd 1.055
105  // false: mean -0.07506772 with sd 1.055
106  // below -0.5 removes around 30% of the peaks
107  // below 0 removes around 50% of the peaks
108  // below 0.5 removes around 70% of the peaks
109  // below 1.0 removes around 85% of the peaks
110  // below 1.5 removes around 93% of the peaks
111  // below 2.0 removes around 97% of the peaks
112  double lda_quick_score =
113  library_corr * -0.5319046 +
114  library_rmsd * 2.1643962 +
115  norm_rt_score * 8.0353047 +
116  xcorr_coelution_score * 0.1458914 +
117  xcorr_shape_score * -1.6901925 +
118  log_sn_score * -0.8002824;
119  return lda_quick_score;
120  }
121 
123  {
124 
125  // LDA average model on 100 2xCrossvalidated runs (0.91 TPR/0.20 FDR)
126  /*
127  double xx_old_lda_prescore =
128  intensity_score * -2.296679 +
129  library_corr * -0.1223876 +
130  library_rmsd * 2.013638 +
131  nr_peaks_score * 0.01683357 +
132  rt_score * 0.00143999 +
133  sn_score * -0.1619762 +
134  total_xic_score * 0.00000003697898 +
135  xcorr_coelution_score * 0.05909583 +
136  xcorr_shape_score * -0.4699841;
137  */
138 
139  return scores.library_corr * -0.34664267 +
140  scores.library_rmsd * 2.98700722 +
141  scores.norm_rt_score * 7.05496384 +
142  scores.xcorr_coelution_score * 0.09445371 +
143  scores.xcorr_shape_score * -5.71823862 +
144  scores.log_sn_score * -0.72989582 +
145  scores.elution_model_fit_score * 1.88443209;
146  }
147 
149  {
150 
151  // Swath - LDA average model on 100 2xCrossvalidated runs (0.76 TPR/0.20 FDR) [without elution model]
152  /*
153  double xx_old_swath_prescore =
154  intensity_score * -3.148838e+00 +
155  library_corr * -7.562403e-02 +
156  library_rmsd * 1.786286e+00 +
157  nr_peaks_score * -7.674263e-03 +
158  rt_score * 1.748377e-03 +
159  sn_score * -1.372636e-01 +
160  total_xic_score * 7.278437e-08 +
161  xcorr_coelution_score * 1.181813e-01 +
162  weighted_coelution_score * -7.661783e-02 +
163  xcorr_shape_score * -6.903933e-02 +
164  weighted_xcorr_shape * -4.234820e-01 +
165  bseries_score * -2.022380e-02 +
166  massdev_score * 2.844948e-02 +
167  massdev_score_weighted * 1.133209e-02 +
168  yseries_score * -9.510874e-02 +
169  isotope_corr * -1.619902e+00 +
170  isotope_overlap * 2.890688e-01 ;
171  */
172 
173  return scores.library_corr * -0.19011762 +
174  scores.library_rmsd * 2.47298914 +
175  scores.norm_rt_score * 5.63906731 +
176  scores.isotope_correlation * -0.62640133 +
177  scores.isotope_overlap * 0.36006925 +
178  scores.massdev_score * 0.08814003 +
179  scores.xcorr_coelution_score * 0.13978311 +
180  scores.xcorr_shape_score * -1.16475032 +
181  scores.yseries_score * -0.19267813 +
182  scores.log_sn_score * -0.61712054;
183  }
184 
185  };
186 
200  class OPENMS_DLLAPI MRMFeatureFinderScoring :
201  public DefaultParamHandler,
202  public ProgressLogger
203  {
204 
205 public:
206 
208 
209 
210  // All the filters expect MSSpectrum<PeakT>, thus we give it an "MSSpectrum"
211  // but filled with Chromatogram Peaks.
212  typedef MSSpectrum<ChromatogramPeak> RichPeakChromatogram; // this is the type in which we store the chromatograms for this analysis
218  typedef MRMTransitionGroup<MSSpectrum <ChromatogramPeak>, TransitionType> MRMTransitionGroupType; // a transition group holds the MSSpectra with the Chromatogram peaks from above
219  typedef std::map<String, MRMTransitionGroupType> TransitionGroupMapType;
221 
224 
227 
228  // pick features in one experiment containing chromatograms
229  // (easy function for wrapping in Python, only uses OpenMS datastructures
230  // and does not return the map)
231  void pickExperiment(MSExperiment<Peak1D> & chromatograms, FeatureMap<Feature>& output, TargetedExperiment& transition_exp_,
233  {
234  OpenSwath::LightTargetedExperiment transition_exp;
235  OpenSwathDataAccessHelper::convertTargetedExp(transition_exp_, transition_exp);
236  TransitionGroupMapType transition_group_map;
237 
240 
241  pickExperiment(chromatogram_ptr, output, transition_exp, trafo, empty_swath_ptr, transition_group_map);
242  }
243 
244  // pick features in one experiment containing chromatograms
246  TransformationDescription trafo, OpenSwath::SpectrumAccessPtr swath_map, TransitionGroupMapType& transition_group_map)
247  {
248  updateMembers_();
249 
250  //
251  // Step 1
252  //
253  // Store the peptide retention times in an intermediate map
254  PeptideRTMap_.clear();
255  for (Size i = 0; i < transition_exp.getPeptides().size(); i++)
256  {
257  PeptideType pep = transition_exp.getPeptides()[i];
258  PeptideRTMap_[pep.id] = pep.rt;
259  PeptideRefMap_[pep.id] = &transition_exp.getPeptides()[i];
260  }
261 
262  // Store the proteins from the input in the output feature map
263  std::vector<ProteinHit> protein_hits;
264  for (Size i = 0; i < transition_exp.getProteins().size(); i++)
265  {
266  const ProteinType& prot = transition_exp.getProteins()[i];
267  ProteinRefMap_[transition_exp.getProteins()[i].id] = &transition_exp.getProteins()[i];
268  ProteinHit prot_hit = ProteinHit();
269  prot_hit.setSequence(prot.sequence);
270  prot_hit.setAccession(prot.id);
271  protein_hits.push_back(prot_hit);
272  }
273 
275  prot_id.setHits(protein_hits);
276  prot_id.setIdentifier(run_identifier);
277  output.getProteinIdentifications().push_back(prot_id);
278 
279  //
280  // Step 2
281  //
282  // Create all MRM transition groups from the individual transitions.
283  mapExperimentToTransitionList(input, transition_exp, transition_group_map, trafo, rt_extraction_window_);
284  int counter = 0;
285  for (TransitionGroupMapType::iterator trgroup_it = transition_group_map.begin(); trgroup_it != transition_group_map.end(); trgroup_it++)
286  {
287  if (trgroup_it->second.getChromatograms().size() > 0) {counter++; }
288  }
289  std::cout << "Will analyse " << counter << " peptides with a total of " << transition_exp.getTransitions().size() << " transitions " << std::endl;
290 
291  //
292  // Step 3
293  //
294  // Go through all transition groups: first create consensus features, then score them
295  Size progress = 0;
296  startProgress(0, transition_group_map.size(), "picking peaks");
297  for (TransitionGroupMapType::iterator trgroup_it = transition_group_map.begin(); trgroup_it != transition_group_map.end(); trgroup_it++)
298  {
299 
300  setProgress(++progress);
301  MRMTransitionGroupType& transition_group = trgroup_it->second;
302  if (transition_group.getChromatograms().size() == 0 || transition_group.getTransitions().size() == 0)
303  {
304  continue;
305  }
306 
307  MRMTransitionGroupPicker trgroup_picker;
308  trgroup_picker.setParameters(param_.copy("TransitionGroupPicker:", true));
309  trgroup_picker.pickTransitionGroup(transition_group);
310  scorePeakgroups_(trgroup_it->second, trafo, swath_map, output);
311 
312  }
313  endProgress();
314 
315  //output.sortByPosition(); // if the exact same order is needed
316  return;
317  }
318 
319  // Map an input experiment (mzML) and transition list (TraML) onto each other
320  // when they share identifiers, e.g. if the transition id is the same as the
321  // chromatogram native id.
322  void mapExperimentToTransitionList(OpenSwath::SpectrumAccessPtr input, OpenSwath::LightTargetedExperiment& transition_exp,
323  TransitionGroupMapType& transition_group_map, TransformationDescription trafo, double rt_extraction_window);
324 
325  void setStrictFlag(bool f)
326  {
327  strict_ = f;
328  }
329 
330 private:
331 
333  template <typename SpectrumT, typename TransitionT>
336  {
337  //std::vector<SignalToNoiseEstimatorMedian<RichPeakChromatogram> > signal_noise_estimators;
339  std::vector<OpenSwath::ISignalToNoisePtr> signal_noise_estimators;
340  std::vector<MRMFeature> feature_list;
341 
342 
343  DoubleReal sn_win_len_ = (DoubleReal)param_.getValue("TransitionGroupPicker:sn_win_len");
344  DoubleReal sn_bin_count_ = (DoubleReal)param_.getValue("TransitionGroupPicker:sn_bin_count");
345  for (Size k = 0; k < transition_group.getChromatograms().size(); k++)
346  {
347  OpenSwath::ISignalToNoisePtr snptr(new OpenMS::SignalToNoiseOpenMS< PeakT >(transition_group.getChromatograms()[k], sn_win_len_, sn_bin_count_));
348  signal_noise_estimators.push_back(snptr);
349  }
350 
351  // get the expected rt value for this peptide
352  double expected_rt = PeptideRTMap_[transition_group.getTransitionGroupID()];
353  TransformationDescription newtr = trafo;
354  newtr.invert();
355  expected_rt = newtr.apply(expected_rt);
356 
357  // Go through all peak groups (found MRM features) and score them
358  for (std::vector<MRMFeature>::iterator mrmfeature = transition_group.getFeaturesMuteable().begin();
359  mrmfeature != transition_group.getFeaturesMuteable().end(); mrmfeature++)
360  {
361 
362  OpenSwath::IMRMFeature* imrmfeature;
363  imrmfeature = new MRMFeatureOpenMS(*mrmfeature);
364 
365  OpenSwath::ITransitionGroup* itransition_group;
366  itransition_group = new TransitionGroupOpenMS<SpectrumT, TransitionT>(transition_group);
367 
368 #ifdef DEBUG_MRMPEAKPICKER
369  std::cout << "000000000000000000000000000000000000000000000000000000000000000000000000000 " << std::endl;
370  std::cout << "scoring feature " << (*mrmfeature) << " == " << mrmfeature->getMetaValue("PeptideRef") <<
371  "[ expected RT " << PeptideRTMap_[mrmfeature->getMetaValue("PeptideRef")] << " / " << expected_rt << " ]" <<
372  " with " << transition_group.size() << " nr transitions and nr chromats " << transition_group.getChromatograms().size() << std::endl;
373 #endif
374 
375  int group_size = boost::numeric_cast<int>(transition_group.size());
376  if (group_size == 0)
377  {
378  throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__,
379  "Error: Transition group " + transition_group.getTransitionGroupID() + " has no chromatograms.");
380  }
381  if (group_size < 2)
382  {
383  std::cerr << "Error: transition group " << transition_group.getTransitionGroupID() << " has less than 2 chromatograms. It has " << group_size << std::endl;
384  continue;
385  //throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Error: transition group " + transition_group.getTransitionGroupID() + " has less than 2 chromatograms.");
386  }
387 
388  // calculate the normalized library intensity (expected value of the intensities)
389  std::vector<double> normalized_library_intensity;
390  transition_group.getLibraryIntensity(normalized_library_intensity);
391  OpenSwath::Scoring::normalize_sum(&normalized_library_intensity[0], boost::numeric_cast<int>(normalized_library_intensity.size()));
392 
393  // calcxcorr -> for each lag do the correlation, normally use lag 0
394  // xcorr_matrix => correlate chromatogram i with chromatogram j
395  bool normalize = true;
396  mrmscore_.initializeXCorrMatrix(imrmfeature, itransition_group, normalize);
397 
398  // XCorr score (coelution)
399  double xcorr_coelution_score = 0;
400  if (use_coelution_score_)
401  {
402  xcorr_coelution_score = mrmscore_.calcXcorrCoelutionScore();
403  mrmfeature->addScore("var_xcorr_coelution", xcorr_coelution_score);
404  }
405 
406  double weighted_coelution_score = 0;
407  if (use_coelution_score_)
408  {
409  weighted_coelution_score = mrmscore_.calcXcorrCoelutionScore_weighted(normalized_library_intensity);
410  mrmfeature->addScore("var_xcorr_coelution_weighted ", weighted_coelution_score);
411  }
412 
413  // XCorr score (shape)
414  // mean over the intensities at the max of the crosscorrelation
415  // FEATURE : weigh by the intensity as done by mQuest
416  // FEATURE : normalize with the intensity at the peak group apex?
417  double xcorr_shape_score = 0;
418  if (use_shape_score_)
419  {
420  xcorr_shape_score = mrmscore_.calcXcorrShape_score();
421  mrmfeature->addScore("var_xcorr_shape", xcorr_shape_score);
422  }
423 
424  double weighted_xcorr_shape = 0;
425  if (use_shape_score_)
426  {
427  weighted_xcorr_shape = mrmscore_.calcXcorrShape_score_weighted(normalized_library_intensity);
428  mrmfeature->addScore("var_xcorr_shape_weighted", weighted_xcorr_shape);
429  }
430 
431  // FEATURE : how should we best calculate correlation between library and experiment?
432  // FEATURE : spectral angle
433  double library_corr = 0, library_rmsd = 0;
434  double library_manhattan, library_dotprod;
435  if (use_library_score_)
436  {
437  mrmscore_.calcLibraryScore(imrmfeature, transition_group.getTransitions(), library_corr, library_rmsd, library_manhattan, library_dotprod);
438  mrmfeature->addScore("var_library_corr", library_corr);
439  mrmfeature->addScore("var_library_rmsd", library_rmsd);
440  mrmfeature->addScore("var_library_manhattan", library_manhattan); // new score
441  mrmfeature->addScore("var_library_dotprod", library_dotprod); // new score
442  }
443 
444  // Retention time score
445  double rt_score = 0, norm_rt_score = 0;
446  if (use_rt_score_)
447  {
448  // get the id, then get the expected and the experimental retention time
449  String native_id = transition_group.getChromatograms()[0].getNativeID();
450  TransitionType tr = transition_group.getTransition(native_id);
451  const PeptideType* pep = PeptideRefMap_[tr.getPeptideRef()];
452  double experimental_rt = mrmfeature->getFeature(native_id).getRT();
453  double normalized_experimental_rt = trafo.apply(experimental_rt);
454  // rt score is delta iRT
455  rt_score = mrmscore_.calcRTScore(*pep, normalized_experimental_rt);
456  norm_rt_score = rt_score / rt_normalization_factor_;
457  mrmfeature->addScore("delta_rt", mrmfeature->getRT() - expected_rt);
458  mrmfeature->addScore("assay_rt", expected_rt);
459  mrmfeature->addScore("norm_RT", normalized_experimental_rt);
460  mrmfeature->addScore("rt_score", rt_score);
461  mrmfeature->addScore("var_norm_rt_score", norm_rt_score);
462  }
463 
464  // Intensity score
465  double intensity_score = 0;
466  if (use_intensity_score_)
467  {
468  intensity_score = mrmfeature->getIntensity() / (double)mrmfeature->getMetaValue("total_xic");
469  mrmfeature->addScore("var_intensity_score", intensity_score);
470  }
471 
472  double total_xic_score = 0;
473  if (use_total_xic_score_)
474  {
475  total_xic_score = (double)mrmfeature->getMetaValue("total_xic");
476  mrmfeature->addScore("total_xic", total_xic_score);
477  }
478 
479  double nr_peaks_score = 0;
480  if (use_nr_peaks_score_)
481  {
482  nr_peaks_score = group_size;
483  mrmfeature->addScore("nr_peaks", nr_peaks_score);
484  }
485 
486  double sn_score = 0, log_sn_score = 0;
487  if (use_sn_score_)
488  {
489  sn_score = mrmscore_.calcSNScore(imrmfeature, signal_noise_estimators);
490  if (sn_score < 1) // fix to make sure, that log(sn_score = 0) = -inf does not occur
491  {
492  log_sn_score = 0;
493  }
494  else
495  {
496  log_sn_score = std::log(sn_score);
497  }
498  mrmfeature->addScore("sn_ratio", sn_score);
499  mrmfeature->addScore("var_log_sn_score", log_sn_score);
500  }
501 
502  OpenSwath_Scores scores;
503  double quick_lda_dismiss = 0;
504  double lda_quick_score = -scores.get_quick_lda_score(library_corr, library_rmsd, norm_rt_score, xcorr_coelution_score, xcorr_shape_score, log_sn_score);
505 
506  if (lda_quick_score < quick_lda_dismiss)
507  {
508  // continue;
509  }
510 
511  double elution_model_fit_score = 0;
512  if (use_elution_model_score_)
513  {
514  elution_model_fit_score = emgscoring_.calcElutionFitScore((*mrmfeature), transition_group);
515  mrmfeature->addScore("var_elution_model_fit_score", elution_model_fit_score);
516  }
517 
518  double xx_lda_prescore;
519  scores.library_corr = library_corr;
520  scores.library_rmsd = library_rmsd;
521  scores.norm_rt_score = norm_rt_score;
522  scores.elution_model_fit_score = elution_model_fit_score;
523  scores.log_sn_score = log_sn_score;
524  scores.xcorr_coelution_score = xcorr_coelution_score;
525  scores.xcorr_shape_score = xcorr_shape_score;
526  xx_lda_prescore = -scores.calculate_lda_prescore(scores);
527 
528  bool swath_present = (swath_map->getNrSpectra() > 0);
529  if (!swath_present)
530  {
531  mrmfeature->addScore("main_var_xx_lda_prelim_score", xx_lda_prescore);
532  mrmfeature->setOverallQuality(xx_lda_prescore);
533  }
534  else
535  {
536  mrmfeature->addScore("xx_lda_prelim_score", xx_lda_prescore);
537  }
538 
539  if (swath_present)
540  {
541  calculateSwathScores_(transition_group, *mrmfeature, swath_map, normalized_library_intensity, scores);
542  }
543 
544 #if 0
545  if (do_local_fdr_)
546  {
547  calculate_local_fdr_scores(transition_group, *mrmfeature, trafo);
548  }
549 #endif
550 
552  // add the peptide hit information to the feature
554 
555  const PeptideType* pep = PeptideRefMap_[transition_group.getTransitions()[0].getPeptideRef()];
556  const ProteinType* prot = ProteinRefMap_[pep->protein_ref];
557 
559  PeptideHit pep_hit_ = PeptideHit();
560 
561  if (pep->getChargeState() != -1)
562  {
563  pep_hit_.setCharge(pep->getChargeState());
564  }
565  pep_hit_.setScore(xx_lda_prescore);
566  if (swath_present)
567  {
568  pep_hit_.setScore(mrmfeature->getScore("xx_swath_prelim_score"));
569  }
570  pep_hit_.setSequence((String)pep->sequence);
571  pep_hit_.addProteinAccession(prot->id);
572  pep_id_.insertHit(pep_hit_);
573  pep_id_.setIdentifier(run_identifier);
574 
575  mrmfeature->getPeptideIdentifications().push_back(pep_id_);
576  mrmfeature->ensureUniqueId();
577  mrmfeature->setMetaValue("PrecursorMZ", transition_group.getTransitions()[0].getPrecursorMZ());
578  mrmfeature->setSubordinates(mrmfeature->getFeatures()); // add all the subfeatures as subordinates
579  double total_intensity = 0, total_peak_apices = 0;
580  for (std::vector<Feature>::iterator sub_it = mrmfeature->getSubordinates().begin(); sub_it != mrmfeature->getSubordinates().end(); sub_it++)
581  {
582  if (!write_convex_hull_) {sub_it->getConvexHulls().clear(); }
583  sub_it->ensureUniqueId();
584  if (sub_it->getMZ() > quantification_cutoff_)
585  {
586  total_intensity += sub_it->getIntensity();
587  total_peak_apices += (DoubleReal)sub_it->getMetaValue("peak_apex_int");
588  }
589  }
590  // overwrite the reported intensities with those above the m/z cutoff
591  mrmfeature->setIntensity(total_intensity);
592  mrmfeature->setMetaValue("peak_apices_sum", total_peak_apices);
593  feature_list.push_back((*mrmfeature));
594 
595  delete imrmfeature;
596  delete itransition_group;
597  }
598 
599  // Order by quality
600  std::sort(feature_list.begin(), feature_list.end(), OpenMS::Feature::OverallQualityLess());
601  std::reverse(feature_list.begin(), feature_list.end());
602 
603  for (Size i = 0; i < feature_list.size(); i++)
604  {
605  if (stop_report_after_feature_ >= 0 && i >= (Size)stop_report_after_feature_) {break; }
606  output.push_back(feature_list[i]);
607  }
608  }
609 
611  OpenSwath::SpectrumPtr getAddedSpectra_(OpenSwath::SpectrumAccessPtr swath_map, double RT, int nr_spectra_to_add)
612  {
613  std::vector<std::size_t> indices = swath_map->getSpectraByRT(RT, 0.0);
614  int closest_idx = boost::numeric_cast<int>(indices[0]);
615  if (indices[0] != 0 &&
616  std::fabs(swath_map->getSpectrumMetaById(boost::numeric_cast<int>(indices[0]) - 1).RT - RT) <
617  std::fabs(swath_map->getSpectrumMetaById(boost::numeric_cast<int>(indices[0])).RT - RT))
618  {
619  closest_idx--;
620  }
621 
622  if (nr_spectra_to_add == 1)
623  {
624  OpenSwath::SpectrumPtr spectrum_ = swath_map->getSpectrumById(closest_idx);
625  return spectrum_;
626  }
627  else
628  {
629  std::vector<OpenSwath::SpectrumPtr> all_spectra;
630  // always add the spectrum 0, then add those right and left
631  all_spectra.push_back(swath_map->getSpectrumById(closest_idx));
632  for (int i = 1; i <= nr_spectra_to_add / 2; i++) // cast to int is intended!
633  {
634  all_spectra.push_back(swath_map->getSpectrumById(closest_idx - i));
635  all_spectra.push_back(swath_map->getSpectrumById(closest_idx + i));
636  }
637  OpenSwath::SpectrumPtr spectrum_ = SpectrumAddition::addUpSpectra(all_spectra, spacing_for_spectra_resampling_, true);
638  return spectrum_;
639  }
640  }
641 
642  template <typename SpectrumT, typename TransitionT>
644  OpenSwath::SpectrumAccessPtr swath_map, std::vector<double>& normalized_library_intensity, OpenSwath_Scores scores)
645  {
646  MRMFeature* mrmfeature = &mrmfeature_;
647 
648  // parameters
649  int by_charge_state = 1; // for which charge states should we check b/y series
650 
651  // find spectrum that is closest to the apex of the peak using binary search
652  OpenSwath::SpectrumPtr spectrum_ = getAddedSpectra_(swath_map, mrmfeature->getRT(), add_up_spectra_);
653  OpenSwath::SpectrumPtr* spectrum = &spectrum_;
654 
655  // Isotope correlation / overlap score: Is this peak part of an
656  // isotopic pattern or is it the monoisotopic peak in an isotopic
657  // pattern?
658  OpenSwath::IMRMFeature* imrmfeature = new MRMFeatureOpenMS(*mrmfeature);
659  double isotope_corr = 0, isotope_overlap = 0;
660  diascoring_.dia_isotope_scores(transition_group.getTransitions(),
661  (*spectrum), imrmfeature, isotope_corr, isotope_overlap);
662  // Mass deviation score
663  double ppm_score = 0, ppm_score_weighted = 0;
664  diascoring_.dia_massdiff_score(transition_group.getTransitions(),
665  (*spectrum), normalized_library_intensity, ppm_score, ppm_score_weighted);
666 
667  // Presence of b/y series score
668  double bseries_score = 0, yseries_score = 0;
669  OpenMS::AASequence aas;
670  OpenSwathDataAccessHelper::convertPeptideToAASequence(*PeptideRefMap_[transition_group.getTransitions()[0].getPeptideRef()], aas);
671  diascoring_.dia_by_ion_score((*spectrum), aas, by_charge_state, bseries_score, yseries_score);
672  mrmfeature->addScore("var_isotope_correlation_score", isotope_corr);
673  mrmfeature->addScore("var_isotope_overlap_score", isotope_overlap);
674 #ifdef DEBUG_MRMPEAKPICKER
675  cout << "added corr isotope_score " << isotope_corr << endl;
676  cout << "added overlap isotope_score " << isotope_overlap << endl;
677 #endif
678 
679  // FEATURE we should not punish so much when one transition is missing!
680  double massdev_score = ppm_score / transition_group.size();
681  double massdev_score_weighted = ppm_score_weighted;
682  mrmfeature->addScore("var_massdev_score", massdev_score);
683  mrmfeature->addScore("var_massdev_score_weighted", massdev_score_weighted);
684 #ifdef DEBUG_MRMPEAKPICKER
685  cout << "added score massdev_score " << massdev_score << endl;
686  cout << "added score weighted massdev_score " << massdev_score_weighted << endl;
687 #endif
688 
689  mrmfeature->addScore("var_bseries_score", bseries_score);
690  mrmfeature->addScore("var_yseries_score", yseries_score);
691 #ifdef DEBUG_MRMPEAKPICKER
692  cout << "added score bseries_score " << bseries_score << endl;
693  cout << "added score yseries_score " << yseries_score << endl;
694 #endif
695 
696  double dotprod_score_dia;
697  double manhatt_score_dia;
698 
699  diascoring_.score_with_isotopes((*spectrum), transition_group.getTransitions(), dotprod_score_dia, manhatt_score_dia);
700 
701  mrmfeature->addScore("var_dotprod_score", dotprod_score_dia);
702  mrmfeature->addScore("var_manhatt_score", manhatt_score_dia);
703 
704  scores.yseries_score = yseries_score;
705  scores.isotope_correlation = isotope_corr;
706  scores.isotope_overlap = isotope_overlap;
707  scores.massdev_score = massdev_score;
708  double xx_swath_prescore = -scores.calculate_swath_lda_prescore(scores);
709  mrmfeature->addScore("main_var_xx_swath_prelim_score", xx_swath_prescore);
710  mrmfeature->setOverallQuality(xx_swath_prescore);
711 #ifdef DEBUG_MRMPEAKPICKER
712  cout << "added xx_swath_prescore (everything above 2 is good) " << xx_swath_prescore << endl;
713 #endif
714  delete imrmfeature;
715  }
716 
717  // void handle_params();
718 
720  void updateMembers_();
721 
722  // Variables
725 
726  // Which scores to use
736 
740 
741  // bool do_local_fdr_;
743  bool strict_;
744 
746 
747  std::map<OpenMS::String, double> PeptideRTMap_;
748  std::map<OpenMS::String, const PeptideType*> PeptideRefMap_;
749  std::map<OpenMS::String, const ProteinType*> ProteinRefMap_;
750 
754  };
755 }
756 
757 #undef run_identifier
758 #endif
double isotope_overlap
Definition: MRMFeatureFinderScoring.h:92
bool use_sn_score_
Definition: MRMFeatureFinderScoring.h:735
std::string getPeptideRef() const
Definition: TransitionExperiment.h:66
std::vector< LightProtein > & getProteins()
Definition: TransitionExperiment.h:135
const double k
Representation of a protein identification run.
Definition: ProteinIdentification.h:61
std::string id
Definition: TransitionExperiment.h:116
double isotope_correlation
Definition: MRMFeatureFinderScoring.h:91
int getChargeState() const
Definition: TransitionExperiment.h:106
const std::vector< TransitionType > & getTransitions() const
Definition: MRMTransitionGroup.h:122
A more convenient string class.
Definition: String.h:56
void setIdentifier(const String &id)
sets the indentifier
void setScore(DoubleReal score)
sets the score of the peptide hit
void setHits(const std::vector< ProteinHit > &hits)
Sets the protein hits.
Size size() const
Definition: MRMTransitionGroup.h:107
OpenSwath::LightPeptide PeptideType
Definition: MRMFeatureFinderScoring.h:215
bool use_shape_score_
Definition: MRMFeatureFinderScoring.h:728
const std::vector< ProteinIdentification > & getProteinIdentifications() const
non-mutable access to the protein identifications
Definition: FeatureMap.h:410
This class implements different scores for peaks found in SRM/MRM.
Definition: MRMScoring.h:75
void setSequence(const String &sequence)
sets the protein sequence
Scoring of an elution peak using an exponentially modified gaussian distribution model.
Definition: EmgScoring.h:60
MRMTransitionGroup< MSSpectrum< ChromatogramPeak >, TransitionType > MRMTransitionGroupType
Definition: MRMFeatureFinderScoring.h:218
OpenMS::DIAScoring diascoring_
Definition: MRMFeatureFinderScoring.h:752
OpenSwath::SpectrumPtr getAddedSpectra_(OpenSwath::SpectrumAccessPtr swath_map, double RT, int nr_spectra_to_add)
Returns the addition of &quot;nr_spectra_to_add&quot; spectra around the given RT.
Definition: MRMFeatureFinderScoring.h:611
OpenSwath::MRMScoring mrmscore_
Definition: MRMFeatureFinderScoring.h:751
const std::vector< SpectrumType > & getChromatograms() const
Definition: MRMTransitionGroup.h:148
A container for features.
Definition: FeatureMap.h:111
OpenSwath::LightTransition TransitionType
Definition: MRMFeatureFinderScoring.h:213
static void convertPeptideToAASequence(const OpenSwath::LightPeptide &peptide, AASequence &aa_sequence)
convert from the LightPeptide to an OpenMS AASequence (with correct modifications) ...
std::map< OpenMS::String, double > PeptideRTMap_
Definition: MRMFeatureFinderScoring.h:747
An implementation of the OpenSWATH Transition Group Access interface using OpenMS.
Definition: MRMFeatureAccessOpenMS.h:109
The MRMFeatureFinder finds and scores peaks of transitions that coelute.
Definition: MRMFeatureFinderScoring.h:200
bool write_convex_hull_
Definition: MRMFeatureFinderScoring.h:742
int add_up_spectra_
Definition: MRMFeatureFinderScoring.h:738
std::map< String, MRMTransitionGroupType > TransitionGroupMapType
Definition: MRMFeatureFinderScoring.h:219
#define run_identifier
Definition: MRMFeatureFinderScoring.h:38
void getLibraryIntensity(std::vector< double > &result) const
Definition: MRMTransitionGroup.h:189
Representation of a peptide/protein sequence.
Definition: AASequence.h:84
QualityLess OverallQualityLess
Compare by quality.
Definition: Feature.h:97
std::string protein_ref
Definition: TransitionExperiment.h:104
CoordinateType getRT() const
Returns the RT coordinate (index 0)
Definition: Peak2D.h:203
void setParameters(const Param &param)
Sets the parameters.
OpenMS::EmgScoring emgscoring_
Definition: MRMFeatureFinderScoring.h:753
OPENSWATHALGO_DLLAPI void normalize(const std::vector< double > &intensities, double normalization_factor, std::vector< double > &normalized_intensities)
Normalize intensities in vector by normalization_factor.
Definition: ITransition.h:65
bool use_nr_peaks_score_
Definition: MRMFeatureFinderScoring.h:734
The representation of a transition group that has information about the individual chromatograms as w...
Definition: MRMTransitionGroup.h:56
int stop_report_after_feature_
Definition: MRMFeatureFinderScoring.h:737
OPENSWATHALGO_DLLAPI void normalize_sum(double x[], unsigned int n)
divide each element of x by the sum of the vector
OpenSwath::LightProtein ProteinType
Definition: MRMFeatureFinderScoring.h:216
double library_rmsd
Definition: MRMFeatureFinderScoring.h:89
Definition: TransitionExperiment.h:99
void scorePeakgroups_(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group, TransformationDescription &trafo, OpenSwath::SpectrumAccessPtr swath_map, FeatureMap< Feature > &output)
Score all peak groups.
Definition: MRMFeatureFinderScoring.h:334
std::string sequence
Definition: TransitionExperiment.h:117
The representation of a 1D spectrum.
Definition: MSSpectrum.h:67
A method or algorithm argument contains illegal values.
Definition: Exception.h:634
OPENSWATHALGO_DLLAPI typedef boost::shared_ptr< ISpectrumAccess > SpectrumAccessPtr
Definition: ANALYSIS/OPENSWATH/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h:73
std::vector< LightTransition > & getTransitions()
Definition: TransitionExperiment.h:125
DoubleReal spacing_for_spectra_resampling_
Definition: MRMFeatureFinderScoring.h:739
void setCharge(Int charge)
sets the charge of the peptide
The MRMTransitionGroupPicker finds peaks in chromatograms that belong to the same precursors...
Definition: MRMTransitionGroupPicker.h:76
DoubleReal quantification_cutoff_
Definition: MRMFeatureFinderScoring.h:724
OPENSWATHALGO_DLLAPI typedef boost::shared_ptr< Spectrum > SpectrumPtr
Definition: ANALYSIS/OPENSWATH/OPENSWATHALGO/DATAACCESS/DataStructures.h:227
Representation of a peptide hit.
Definition: PeptideHit.h:54
static OpenSwath::SpectrumAccessPtr getSpectrumAccessOpenMSPtr(OpenMS::MSExperiment< OpenMS::Peak1D > &exp)
Simple Factory method to get a SpectrumAccess Ptr from an MSExperiment.
double log_sn_score
Definition: MRMFeatureFinderScoring.h:97
An implementation of the OpenSWATH SignalToNoise Access interface using OpenMS.
Definition: MRMFeatureAccessOpenMS.h:152
Scoring of an spectrum at the peak apex of an chromatographic elution peak.
Definition: DIAScoring.h:76
double massdev_score
Definition: MRMFeatureFinderScoring.h:93
bool use_coelution_score_
Definition: MRMFeatureFinderScoring.h:727
DoubleReal rt_normalization_factor_
Definition: MRMFeatureFinderScoring.h:745
double get_quick_lda_score(double library_corr, double library_rmsd, double norm_rt_score, double xcorr_coelution_score, double xcorr_shape_score, double log_sn_score)
Definition: MRMFeatureFinderScoring.h:99
double yseries_score
Definition: MRMFeatureFinderScoring.h:96
std::map< OpenMS::String, const PeptideType * > PeptideRefMap_
Definition: MRMFeatureFinderScoring.h:748
std::vector< MRMFeature > & getFeaturesMuteable()
Definition: MRMTransitionGroup.h:179
void setStrictFlag(bool f)
Definition: MRMFeatureFinderScoring.h:325
bool use_rt_score_
Definition: MRMFeatureFinderScoring.h:729
DoubleReal apply(DoubleReal value) const
Applies the transformation to value.
static void convertTargetedExp(const OpenMS::TargetedExperiment &transition_exp_, OpenSwath::LightTargetedExperiment &transition_exp)
convert from the OpenMS Targeted experiment to the light Targeted Experiment
std::string id
Definition: TransitionExperiment.h:105
void calculateSwathScores_(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group, MRMFeature &mrmfeature_, OpenSwath::SpectrumAccessPtr swath_map, std::vector< double > &normalized_library_intensity, OpenSwath_Scores scores)
Definition: MRMFeatureFinderScoring.h:643
double xcorr_coelution_score
Definition: MRMFeatureFinderScoring.h:94
void setSequence(const AASequence &sequence)
sets the peptide sequence
void setAccession(const String &accession)
sets the accession of the protein
Representation of a protein hit.
Definition: ProteinHit.h:54
Definition: TransitionExperiment.h:46
std::vector< LightPeptide > & getPeptides()
Definition: TransitionExperiment.h:130
DoubleReal rt_extraction_window_
Definition: MRMFeatureFinderScoring.h:723
Definition: ITransition.h:56
std::map< OpenMS::String, const ProteinType * > ProteinRefMap_
Definition: MRMFeatureFinderScoring.h:749
void pickExperiment(MSExperiment< Peak1D > &chromatograms, FeatureMap< Feature > &output, TargetedExperiment &transition_exp_, TransformationDescription trafo, MSExperiment< Peak1D > &swath_map)
Definition: MRMFeatureFinderScoring.h:231
void setOverallQuality(QualityType q)
Set the overall quality.
A structure to hold the different scores computed by the FeatureFinder.
Definition: MRMFeatureFinderScoring.h:85
OpenSwath::LightTargetedExperiment TargetedExpType
Definition: MRMFeatureFinderScoring.h:214
bool use_intensity_score_
Definition: MRMFeatureFinderScoring.h:732
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
void addScore(const String &score_name, double score)
set a single peakgroup score
double xcorr_shape_score
Definition: MRMFeatureFinderScoring.h:95
void pickTransitionGroup(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group)
Definition: MRMTransitionGroupPicker.h:202
Base class for all classes that want to report their progess.
Definition: ProgressLogger.h:56
double calculate_swath_lda_prescore(OpenSwath_Scores scores)
Definition: MRMFeatureFinderScoring.h:148
This class stores an prediction of an SRM/MRM transition.
Definition: TargetedExperiment.h:53
bool SortDoubleDoublePairFirst(const std::pair< double, double > &left, const std::pair< double, double > &right)
double norm_rt_score
Definition: MRMFeatureFinderScoring.h:90
OPENSWATHALGO_DLLAPI typedef boost::shared_ptr< ISignalToNoise > ISignalToNoisePtr
Definition: ITransition.h:78
OpenSwath::LightModification ModificationType
Definition: MRMFeatureFinderScoring.h:217
Generic description of a coordinate transformation.
Definition: TransformationDescription.h:59
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:90
void insertHit(const PeptideHit &hit)
Appends a peptide hit.
double rt
Definition: TransitionExperiment.h:101
double elution_model_fit_score
Definition: MRMFeatureFinderScoring.h:87
bool use_total_xic_score_
Definition: MRMFeatureFinderScoring.h:733
void addProteinAccession(const String &accession)
adds an accession of a protein which contains this peptide hit
A multi-chromatogram MRM feature.
Definition: MRMFeature.h:50
bool use_elution_model_score_
Definition: MRMFeatureFinderScoring.h:731
Definition: TransitionExperiment.h:114
Definition: TransitionExperiment.h:120
const String & getTransitionGroupID() const
Definition: MRMTransitionGroup.h:112
std::string sequence
Definition: TransitionExperiment.h:103
double library_corr
Definition: MRMFeatureFinderScoring.h:88
MSSpectrum< ChromatogramPeak > RichPeakChromatogram
Type definitions.
Definition: MRMFeatureFinderScoring.h:212
const TransitionType & getTransition(String key)
Definition: MRMTransitionGroup.h:138
void invert()
Computes an (approximate) inverse of the transformation.
void pickExperiment(OpenSwath::SpectrumAccessPtr input, FeatureMap< Feature > &output, OpenSwath::LightTargetedExperiment &transition_exp, TransformationDescription trafo, OpenSwath::SpectrumAccessPtr swath_map, TransitionGroupMapType &transition_group_map)
Definition: MRMFeatureFinderScoring.h:245
static OpenSwath::SpectrumPtr addUpSpectra(std::vector< OpenSwath::SpectrumPtr > all_spectra, double sampling_rate, double filter_zeros)
adds up a list of Spectra by resampling them and then addition of intensities
void setIdentifier(const String &id)
Sets the identifier.
Represents the peptide hits for a spectrum.
Definition: PeptideIdentification.h:63
double calculate_lda_prescore(OpenSwath_Scores scores)
Definition: MRMFeatureFinderScoring.h:122
bool use_library_score_
Definition: MRMFeatureFinderScoring.h:730
double DoubleReal
Double-precision real type.
Definition: Types.h:118
bool strict_
Definition: MRMFeatureFinderScoring.h:743
Definition: TransitionExperiment.h:93
An implementation of the OpenSWATH MRM Feature Access interface using OpenMS.
Definition: MRMFeatureAccessOpenMS.h:84

OpenMS / TOPP release 1.11.1 Documentation generated on Thu Nov 14 2013 11:19:17 using doxygen 1.8.5