Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
OfflinePrecursorIonSelection.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-2015.
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: Alexandra Zerck $
32 // $Authors: Alexandra Zerck $
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_ANALYSIS_TARGETED_OFFLINEPRECURSORIONSELECTION_H
36 #define OPENMS_ANALYSIS_TARGETED_OFFLINEPRECURSORIONSELECTION_H
37 
38 
45 
46 namespace OpenMS
47 {
48  class PeptideIdentification;
49  class ProteinIdentification;
50  class String;
51 
52 
62  class OPENMS_DLLAPI OfflinePrecursorIonSelection :
63  public DefaultParamHandler
64  {
65 public:
67 
70 
80  template <typename InputPeakType>
81  void makePrecursorSelectionForKnownLCMSMap(const FeatureMap& features,
82  const MSExperiment<InputPeakType>& experiment,
84  std::set<Int>& charges_set,
85  bool feature_based);
86 
94  template <typename InputPeakType>
95  void getMassRanges(const FeatureMap& features,
96  const MSExperiment<InputPeakType>& experiment,
97  std::vector<std::vector<std::pair<Size, Size> > >& indices);
98 
99  void createProteinSequenceBasedLPInclusionList(String include, String rt_model_file, String pt_model_file, FeatureMap& precursors);
100 
102  {
103  solver_ = solver;
104  std::cout << " LPSolver set to " << solver_ << std::endl;
105  }
106 
108  {
109  return solver_;
110  }
111 
112 private:
113 
114  typedef std::map<std::pair<double, double>, int, PairComparatorSecondElement<std::pair<double, double> > > ExclusionListType_;
115 
119  template <typename InputPeakType>
120  void calculateXICs_(const FeatureMap& features,
121  const std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
122  const MSExperiment<InputPeakType>& experiment,
123  const std::set<Int>& charges_set,
124  std::vector<std::vector<std::pair<Size, double> > >& xics);
125 
129  template <typename InputPeakType>
130  void checkMassRanges_(std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
131  const MSExperiment<InputPeakType>& experiment);
132 
134  void updateExclusionList_(ExclusionListType_& exclusion_list) const;
135 
137  };
138 
139  template <typename InputPeakType>
141  {
142  bool enclose_hit = false;
143  const std::vector<ConvexHull2D>& hulls = f.getConvexHulls();
144  for (Size i = 0; i < hulls.size(); ++i)
145  {
146  if (hulls[i].getBoundingBox().encloses(rt, mz))
147  {
148  enclose_hit = true;
149  return enclose_hit;
150  }
151  }
152  return enclose_hit;
153  }
154 
155  template <typename InputPeakType>
157  const MSExperiment<InputPeakType>& experiment,
158  std::vector<std::vector<std::pair<Size, Size> > >& indices)
159  {
160  if (experiment.empty())
161  throw Exception::InvalidSize(__FILE__, __LINE__, __PRETTY_FUNCTION__, 0);
162  for (Size f = 0; f < features.size(); ++f)
163  {
164  std::vector<std::pair<Size, Size> > vec;
165 
166  for (Size rt = 0; rt < experiment.size(); ++rt)
167  {
168  // is scan relevant?
169  if (!enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), features[f].getMZ()))
170  continue;
171 
172  std::pair<Size, Size> start;
173  std::pair<Size, Size> end;
174  bool start_found = false;
175  bool end_found = false;
176  typename MSSpectrum<InputPeakType>::ConstIterator mz_iter = experiment[rt].MZBegin(features[f].getMZ());
177  typename MSSpectrum<InputPeakType>::ConstIterator mz_end = mz_iter;
178  if (mz_iter == experiment[rt].end())
179  continue;
180  // check to the left
181  while (enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), mz_iter->getMZ()))
182  {
183  start_found = true;
184  start.first = rt;
185  start.second = distance(experiment[rt].begin(), mz_iter);
186  if (mz_iter == experiment[rt].begin())
187  break;
188  --mz_iter;
189  end_found = true;
190  end.first = rt;
191  end.second = start.second;
192  }
193  // and now to the right
194  while (mz_end != experiment[rt].end() && enclosesBoundingBox<InputPeakType>(features[f], experiment[rt].getRT(), mz_end->getMZ()))
195  {
196  end_found = true;
197  end.first = rt;
198  end.second = distance(experiment[rt].begin(), mz_end);
199  ++mz_end;
200  }
201  if (start_found && end_found)
202  {
203  vec.push_back(start);
204  vec.push_back(end);
205  }
206 #ifdef DEBUG_OPS
207  else if (start_found || end_found)
208  {
209  std::cout << "start " << start_found << " end " << end_found << std::endl;
210  std::cout << "feature: " << f << " rt: " << rt << std::endl;
211  }
212 #endif
213  }
214 #ifdef DEBUG_OPS
215  if (vec.size() > 0)
216  {
217  std::cout << vec.size() << " / 2 scans" << std::endl;
218  for (Size i = 0; i < vec.size(); i += 2)
219  {
220  std::cout << "Feature " << f << " RT : " << vec[i].first
221  << " MZ : " << experiment[vec[i].first][vec[i].second].getMZ() << " "
222  << experiment[vec[i + 1].first][vec[i + 1].second].getMZ() << std::endl;
223  }
224  }
225 #endif
226  if (vec.empty())
227  {
228 #ifdef DEBUG_OPS
229  std::cout << "According to the convex hulls no mass traces found for this feature->estimate!"
230  << features[f].getRT() << " " << features[f].getMZ() << " " << features[f].getCharge() << std::endl;
231 #endif
232  // we estimate the convex hull
233  typename MSExperiment<InputPeakType>::ConstIterator spec_iter = experiment.RTBegin(features[f].getRT());
234  if (spec_iter == experiment.end())
235  --spec_iter;
236 
237  double dist1 = fabs(spec_iter->getRT() - features[f].getRT());
238  double dist2 = std::numeric_limits<double>::max();
239  double dist3 = std::numeric_limits<double>::max();
240  if ((spec_iter + 1) != experiment.end())
241  {
242  dist2 = fabs((spec_iter + 1)->getRT() - features[f].getRT());
243  }
244  if (spec_iter != experiment.begin())
245  {
246  dist3 = fabs((spec_iter - 1)->getRT() - features[f].getRT());
247  }
248  if (dist3 <= dist1 && dist3 <= dist2)
249  {
250  --spec_iter;
251  }
252  else if (dist2 <= dist3 && dist2 <= dist1)
253  {
254  ++spec_iter;
255  }
256  std::pair<Size, Size> start;
257  std::pair<Size, Size> end;
258  start.first = distance(experiment.begin(), spec_iter);
259  end.first = start.first;
260 
261  typename MSSpectrum<InputPeakType>::ConstIterator mz_iter = spec_iter->MZBegin(features[f].getMZ());
262  if (spec_iter->begin() == spec_iter->end())
263  {
264  indices.push_back(vec);
265  continue;
266  }
267  if (mz_iter == spec_iter->end() || (mz_iter->getMZ() > features[f].getMZ() && mz_iter != spec_iter->begin()))
268  --mz_iter;
269  while (mz_iter != spec_iter->begin())
270  {
271  if (fabs((mz_iter - 1)->getMZ() - features[f].getMZ()) < 0.5)
272  --mz_iter;
273  else
274  break;
275  }
276  start.second = distance(spec_iter->begin(), mz_iter);
277  typename MSSpectrum<InputPeakType>::ConstIterator mz_end = mz_iter;
278 #ifdef DEBUG_OPS
279  std::cout << features[f].getMZ() << " Start: " << experiment[start.first].getRT() << " " << experiment[start.first][start.second].getMZ();
280 #endif
281  Int charge = features[f].getCharge();
282  if (charge == 0)
283  charge = 1;
284  while (mz_end + 1 != spec_iter->end())
285  {
286  if (fabs((mz_end + 1)->getMZ() - features[f].getMZ()) < 3.0 / (double)charge)
287  ++mz_end;
288  else
289  break;
290  }
291  end.second = distance(spec_iter->begin(), mz_end);
292 #ifdef DEBUG_OPS
293  std::cout << "\tEnd: " << experiment[end.first].getRT() << " " << experiment[end.first][end.second].getMZ() << std::endl;
294 #endif
295  vec.push_back(start);
296  vec.push_back(end);
297  }
298 
299  indices.push_back(vec);
300  }
301  // eliminate nearby peaks
302  if (param_.getValue("exclude_overlapping_peaks") == "true")
303  checkMassRanges_(indices, experiment);
304  }
305 
306  template <typename InputPeakType>
308  const std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
309  const MSExperiment<InputPeakType>& experiment,
310  const std::set<Int>& charges_set,
311  std::vector<std::vector<std::pair<Size, double> > >& xics)
312  {
313  xics.clear();
314  xics.resize(experiment.size());
315  // for each feature
316  for (Size f = 0; f < mass_ranges.size(); ++f)
317  {
318  // is charge valid
319  if (charges_set.count(features[f].getCharge()) < 1)
320  {
321  continue;
322  }
323  // go through all scans where the feature occurs
324  for (Size s = 0; s < mass_ranges[f].size(); s += 2)
325  {
326  // sum intensity over all raw datapoints belonging to the feature in the current scan
327  double weight = 0.;
328  for (Size j = mass_ranges[f][s].second; j <= mass_ranges[f][s + 1].second; ++j)
329  {
330  weight += experiment[mass_ranges[f][s].first][j].getIntensity();
331  }
332  // enter xic in the vector for scan
333  xics[mass_ranges[f][s].first].push_back(std::make_pair(f, weight));
334  }
335  }
336 
337  for (Size s = 0; s < xics.size(); ++s)
338  {
339  sort(xics[s].begin(), xics[s].end(), PairComparatorSecondElement<std::pair<Size, double> >());
340  }
341  }
342 
343  template <typename InputPeakType>
345  const MSExperiment<InputPeakType>& experiment,
347  std::set<Int>& charges_set, // allowed charges
348  bool feature_based) // (MALDI with ILP) vs. scan_based (ESI)
349  {
350 
351  const double mz_isolation_window = param_.getValue("mz_isolation_window");
352  const double min_mz_peak_distance = param_.getValue("min_mz_peak_distance");
353 
354  double rt_dist = 0.;
355  if (experiment.size() > 1)
356  {
357  rt_dist = experiment[1].getRT() - experiment[0].getRT();
358  }
359 
360  // feature based selection (e.g. with LC-MALDI)
361  if (feature_based)
362  {
363  // create ILP
364  PSLPFormulation ilp_wrapper;
365 
366  // get the mass ranges for each features for each scan it occurs in
367  std::vector<std::vector<std::pair<Size, Size> > > indices;
368  getMassRanges(features, experiment, indices);
369 
370  std::vector<IndexTriple> variable_indices;
371  std::vector<int> solution_indices;
372  ilp_wrapper.createAndSolveILPForKnownLCMSMapFeatureBased(features, experiment, variable_indices,
373  indices, charges_set,
374  param_.getValue("ms2_spectra_per_rt_bin"),
375  solution_indices);
376 
377  sort(variable_indices.begin(), variable_indices.end(), PSLPFormulation::IndexLess());
378 #ifdef DEBUG_OPS
379  std::cout << "best_solution " << std::endl;
380 #endif
381  // print best solution
382  // create inclusion list
383  for (Size i = 0; i < solution_indices.size(); ++i)
384  {
385  Size feature_index = variable_indices[solution_indices[i]].feature;
386  Size feature_scan_idx = variable_indices[solution_indices[i]].scan;
387  typename MSExperiment<InputPeakType>::ConstIterator scan = experiment.begin() + feature_scan_idx;
389  Precursor p;
390  std::vector<Precursor> pcs;
391  p.setIntensity(features[feature_index].getIntensity());
392  p.setMZ(features[feature_index].getMZ());
393  p.setCharge(features[feature_index].getCharge());
394  pcs.push_back(p);
395  ms2_spec.setPrecursors(pcs);
396  ms2_spec.setRT(scan->getRT() + rt_dist / 2.0);
397  ms2_spec.setMSLevel(2);
398  // link ms2 spectrum with features overlapping its precursor
399  // Warning: this depends on the current order of features in the map
400  // Attention: make sure to name ALL features that overlap, not only one!
401  ms2_spec.setMetaValue("parent_feature_ids", ListUtils::create<Int>(String(feature_index)));
402  ms2.addSpectrum(ms2_spec);
403  std::cout << " MS2 spectra generated at: " << scan->getRT() << " x " << p.getMZ() << "\n";
404 
405  }
406 #ifdef DEBUG_OPS
407  std::cout << solution_indices.size() << " out of " << features.size()
408  << " precursors are in best solution.\n";
409 #endif
410  }
411  else // scan based selection (take the x highest signals for each spectrum)
412  {
413 #ifdef DEBUG_OPS
414  std::cout << "scan based precursor selection" << std::endl;
415 #endif
416  // if the highest signals for each scan shall be selected we don't need an ILP formulation
417 
418  //cache the values for each feature
419  std::vector<DoubleList> feature_elution_bounds;
420  std::vector<DoubleList> elution_profile_intensities;
421  std::vector<DoubleList> isotope_intensities;
422 
423  bool meta_values_present = false;
424 
425  if (!features.empty() &&
426  features[0].metaValueExists("elution_profile_bounds") &&
427  features[0].metaValueExists("elution_profile_intensities") &&
428  features[0].metaValueExists("isotope_intensities"))
429  {
430  for (Size feat = 0; feat < features.size(); ++feat)
431  {
432  feature_elution_bounds.push_back(features[feat].getMetaValue("elution_profile_bounds"));
433  elution_profile_intensities.push_back(features[feat].getMetaValue("elution_profile_intensities"));
434  isotope_intensities.push_back(features[feat].getMetaValue("isotope_intensities"));
435  }
436  meta_values_present = true;
437  }
438 
439  // for each feature: cache for which scans it has to be considered
440  std::vector<std::vector<Size> > scan_features(experiment.size());
441 
442  for (Size feat_idx = 0; feat_idx < features.size(); ++feat_idx)
443  {
444  if (charges_set.count(features[feat_idx].getCharge()))
445  {
446  Size lower_rt = features[feat_idx].getConvexHull().getBoundingBox().minX();
447  Size upper_rt = features[feat_idx].getConvexHull().getBoundingBox().maxX();
449  for (it = experiment.RTBegin(lower_rt); it != experiment.RTEnd(upper_rt); ++it)
450  {
451  scan_features[it - experiment.begin()].push_back(feat_idx);
452  }
453  }
454  }
455 
456  bool dynamic_exclusion = param_.getValue("Exclusion:use_dynamic_exclusion").toBool();
457  ExclusionListType_ exclusion_list;
458  int exclusion_scan_count = (int)(floor((double)param_.getValue("Exclusion:exclusion_time") / (double) rt_dist));
459  if (!dynamic_exclusion)
460  {
461  // if the dynamic exclusion if not active we use the exclusion list to guarantee no two peaks within 'min_mz_peak_distance' are selected for single scan
462  exclusion_scan_count = 0;
463  }
464 
465  // cache bounding boxes of features and mass traces (mass trace bb are also widened for effective discovery of enclosing peaks in intervals)
466  std::map<Size, typename OpenMS::DBoundingBox<2> > bounding_boxes_f;
467  std::map<std::pair<Size, Size>, typename OpenMS::DBoundingBox<2> > bounding_boxes;
468  for (Size feature_num = 0; feature_num < features.size(); ++feature_num)
469  {
470  if (charges_set.count(features[feature_num].getCharge()))
471  {
472  bounding_boxes_f.insert(std::make_pair(feature_num, features[feature_num].getConvexHull().getBoundingBox()));
473  const std::vector<ConvexHull2D> mass_traces = features[feature_num].getConvexHulls();
474  for (Size mass_trace_num = 0; mass_trace_num < mass_traces.size(); ++mass_trace_num)
475  {
476  typename OpenMS::DBoundingBox<2> tmp_bbox = mass_traces[mass_trace_num].getBoundingBox();
477  tmp_bbox.setMinY(tmp_bbox.minY() - mz_isolation_window);
478  tmp_bbox.setMaxY(tmp_bbox.maxY() + mz_isolation_window);
479  bounding_boxes.insert(std::make_pair(std::make_pair(feature_num, mass_trace_num), tmp_bbox));
480  }
481  }
482  }
483 
484  Size max_spec = (Int)param_.getValue("ms2_spectra_per_rt_bin");
485  // get best x signals for each scan
486  for (Size i = 0; i < experiment.size(); ++i)
487  {
488 #ifdef DEBUG_OPS
489  std::cout << "scan " << experiment[i].getRT() << ":";
490 #endif
491  // decrease scan counter by 1, and remove if entry on time-out
492  updateExclusionList_(exclusion_list);
493 
494  MSSpectrum<InputPeakType> scan = experiment[i]; // copy & sort
495  scan.sortByIntensity(true);
496  Size selected_peaks = 0;
497 
498  // walk over scan
499  double peak_rt = scan.getRT();
500  for (Size j = 0; j < scan.size(); ++j)
501  {
502  if (selected_peaks >= max_spec) break;
503 
504  double peak_mz = scan[j].getMZ();
505 
506  // is peak on exclusion? (this relies on 'exclusion_list' being sorted by its right-interval!)
507  ExclusionListType_::iterator it_low = exclusion_list.lower_bound(std::make_pair(peak_mz, peak_mz)); // find the first exclusion entry which ENDS(!) behind this peak
508  if (it_low != exclusion_list.end() && it_low->first.first <= peak_mz) // does it START before this peak?
509  { // then this peak is within the window
510  continue;
511  }
512  // --> no, then fragment it
513  ++selected_peaks;
514 
515  // find all features (mass traces that are in the window around peak_mz)
517  std::vector<Precursor> pcs;
518  std::set<std::pair<Size, Size> > selected_mt;
519  IntList parent_feature_ids;
520 
521  double local_mz = peak_mz;
522  // std::cerr<<"MZ pos: "<<local_mz<<std::endl;
523  for (Size scan_feat_id = 0; scan_feat_id < scan_features[i].size(); ++scan_feat_id)
524  {
525  Size feature_num = scan_features[i][scan_feat_id];
526  if (bounding_boxes_f[feature_num].encloses(peak_rt, local_mz))
527  {
528  //find a mass trace enclosing the point
529  double feature_intensity = 0;
530  for (Size mass_trace_num = 0; mass_trace_num < features[feature_num].getConvexHulls().size(); ++mass_trace_num)
531  {
532  if (bounding_boxes[std::make_pair(feature_num, mass_trace_num)].encloses(DPosition<2>(peak_rt, local_mz)))
533  {
534  double elu_factor = 1.0, iso_factor = 1.0;
535  //get the intensity factor for the position in the elution profile
536  if (meta_values_present)
537  {
538  DoubleList xxx = elution_profile_intensities[feature_num];
539  DoubleList yyy = feature_elution_bounds[feature_num];
540 // std::cout << "PEAKRT: " << peak_rt << std::endl;
541 // std::cout << "Max: " << yyy[3] << " vs. " << bounding_boxes_f[feature_num].maxX() << std::endl;
542 // std::cout << "Min: " << yyy[1] << " vs. " << bounding_boxes_f[feature_num].minX() << std::endl;
543  OPENMS_PRECONDITION(i - yyy[0] < xxx.size(), "Tried to access invalid index for elution factor");
544  elu_factor = xxx[i - yyy[0]]; // segfault here: "i-yyy[0]" yields invalid index
545  iso_factor = isotope_intensities[feature_num][mass_trace_num];
546  }
547  feature_intensity += features[feature_num].getIntensity() * iso_factor * elu_factor;
548  }
549  }
550  Precursor p;
551  p.setIntensity(feature_intensity);
552  p.setMZ(features[feature_num].getMZ());
553  p.setCharge(features[feature_num].getCharge());
554  pcs.push_back(p);
555  parent_feature_ids.push_back((Int)feature_num);
556  }
557  }
558 
559  if (!pcs.empty())
560  {
561  //std::cerr<<"scan "<<i<<" added spectrum for features: "<<parent_feature_ids<<std::endl;
562  ms2_spec.setPrecursors(pcs);
563  ms2_spec.setMSLevel(2);
564  ms2_spec.setRT(experiment[i].getRT() + rt_dist / 2.0); //(selected_peaks+1)*rt_dist/(max_spec+1) );
565  ms2_spec.setMetaValue("parent_feature_ids", parent_feature_ids);
566  ms2.addSpectrum(ms2_spec);
567  }
568 
569  // add m/z window to exclusion list
570  exclusion_list.insert(std::make_pair(std::make_pair(peak_mz - min_mz_peak_distance, peak_mz + min_mz_peak_distance), exclusion_scan_count + 1));
571 
572  }
573  }
574  }
575  }
576 
577  template <typename InputPeakType>
578  void OfflinePrecursorIonSelection::checkMassRanges_(std::vector<std::vector<std::pair<Size, Size> > >& mass_ranges,
579  const MSExperiment<InputPeakType>& experiment)
580  {
581  std::vector<std::vector<std::pair<Size, Size> > > checked_mass_ranges;
582  double min_mz_peak_distance = param_.getValue("min_mz_peak_distance");
583  checked_mass_ranges.reserve(mass_ranges.size());
584  for (Size f = 0; f < mass_ranges.size(); ++f)
585  {
586  std::vector<std::pair<Size, Size> > checked_mass_ranges_f;
587  for (Size s_idx = 0; s_idx < mass_ranges[f].size(); s_idx += 2)
588  {
589  Size s = mass_ranges[f][s_idx].first;
590  bool overlapping_features = false;
592  // check if other features overlap with this feature in the current scan
594  const InputPeakType& peak_left_border = experiment[s][mass_ranges[f][s_idx].second];
595  const InputPeakType& peak_right_border = experiment[s][mass_ranges[f][s_idx + 1].second];
596  for (Size fmr = 0; fmr < mass_ranges.size(); ++fmr)
597  {
598  if (fmr == f)
599  continue;
600  for (Size mr = 0; mr < mass_ranges[fmr].size(); mr += 2)
601  {
602  if (mass_ranges[fmr][mr].first == s) // same spectrum
603  {
604  const InputPeakType& tmp_peak_left = experiment[s][mass_ranges[fmr][mr].second];
605  const InputPeakType& tmp_peak_right = experiment[s][mass_ranges[fmr][mr + 1].second];
606 #ifdef DEBUG_OPS
607  std::cout << tmp_peak_left.getMZ() << " < "
608  << peak_left_border.getMZ() - min_mz_peak_distance << " && "
609  << tmp_peak_right.getMZ() << " < "
610  << peak_left_border.getMZ() - min_mz_peak_distance << " ? "
611  << (tmp_peak_left.getMZ() < peak_left_border.getMZ() - min_mz_peak_distance &&
612  tmp_peak_right.getMZ() < peak_left_border.getMZ() - min_mz_peak_distance)
613  << " || "
614  << tmp_peak_left.getMZ() << " > "
615  << peak_right_border.getMZ() + min_mz_peak_distance << " && "
616  << tmp_peak_right.getMZ() << " > "
617  << peak_right_border.getMZ() + min_mz_peak_distance << " ? "
618  << (tmp_peak_left.getMZ() > peak_right_border.getMZ() + min_mz_peak_distance &&
619  tmp_peak_right.getMZ() > peak_right_border.getMZ() + min_mz_peak_distance)
620  << std::endl;
621 #endif
622  // all other features have to be either completely left or
623  // right of the current feature
624  if (!((tmp_peak_left.getMZ() < peak_left_border.getMZ() - min_mz_peak_distance &&
625  tmp_peak_right.getMZ() < peak_left_border.getMZ() - min_mz_peak_distance) ||
626  (tmp_peak_left.getMZ() > peak_right_border.getMZ() + min_mz_peak_distance &&
627  tmp_peak_right.getMZ() > peak_right_border.getMZ() + min_mz_peak_distance)))
628  {
629 #ifdef DEBUG_OPS
630  std::cout << "found overlapping peak" << std::endl;
631 #endif
632  overlapping_features = true;
633  break;
634  }
635  }
636  }
637  }
638  if (!overlapping_features)
639  {
640 #ifdef DEBUG_OPS
641  std::cout << "feature in spec ok" << mass_ranges[f][s_idx].second << " in spec "
642  << mass_ranges[f][s_idx].first << std::endl;
643 #endif
644  checked_mass_ranges_f.insert(checked_mass_ranges_f.end(),
645  mass_ranges[f].begin() + s_idx,
646  mass_ranges[f].begin() + s_idx + 2);
647  }
648  }
649  checked_mass_ranges.push_back(checked_mass_ranges_f);
650  }
651  mass_ranges.swap(checked_mass_ranges);
652  }
653 
654 
655 }
656 
657 #endif // OPENMS_ANALYSIS_ID_OFFLINEPRECURSORIONSELECTION_H
CoordinateType maxY() const
Accessor for max_ coordinate maximum.
Definition: DIntervalBase.h:259
void setMetaValue(const String &name, const DataValue &value)
sets the DataValue corresponding to a name
A more convenient string class.
Definition: String.h:57
Precursor meta information.
Definition: Precursor.h:56
void setMinY(CoordinateType const c)
Mutator for max_ coordinate of the smaller point.
Definition: DIntervalBase.h:272
Size size() const
Definition: MSExperiment.h:117
std::vector< double > DoubleList
Vector of double precision real types.
Definition: ListUtils.h:66
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:150
SOLVER
Definition: LPWrapper.h:129
void updateExclusionList_(ExclusionListType_ &exclusion_list) const
reduce scan count for each entry, and remove every entry which has reached 0 counts ...
Definition: PSLPFormulation.h:141
A container for features.
Definition: FeatureMap.h:93
#define OPENMS_PRECONDITION(condition, message)
Precondition macro.
Definition: openms/include/OpenMS/CONCEPT/Macros.h:107
void setMaxY(CoordinateType const c)
Mutator for max_ coordinate of the larger point.
Definition: DIntervalBase.h:286
LPWrapper::SOLVER getLPSolver()
Definition: OfflinePrecursorIonSelection.h:107
CoordinateType getMZ() const
Non-mutable access to m/z.
Definition: Peak1D.h:114
ContainerType::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSSpectrum.h:123
Iterator begin()
Definition: MSExperiment.h:147
bool enclosesBoundingBox(const Feature &f, typename MSExperiment< InputPeakType >::CoordinateType rt, typename MSExperiment< InputPeakType >::CoordinateType mz)
Definition: OfflinePrecursorIonSelection.h:140
Struct that holds the indices of the precursors in the feature map and the ilp formulation.
Definition: PSLPFormulation.h:69
std::vector< Int > IntList
Vector of signed integers.
Definition: ListUtils.h:59
void createAndSolveILPForKnownLCMSMapFeatureBased(const FeatureMap &features, const MSExperiment< InputPeakType > &experiment, std::vector< IndexTriple > &variable_indices, std::vector< std::vector< std::pair< Size, Size > > > &mass_ranges, std::set< Int > &charges_set, UInt ms2_spectra_per_rt_bin, std::vector< int > &solution_indices)
Encode ILP formulation for a given LC-MS map, but unknown protein sample.
Definition: PSLPFormulation.h:292
std::map< std::pair< double, double >, int, PairComparatorSecondElement< std::pair< double, double > > > ExclusionListType_
Definition: OfflinePrecursorIonSelection.h:114
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
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
ConstIterator RTEnd(CoordinateType rt) const
Fast search for spectrum range end (returns the past-the-end iterator)
Definition: MSExperiment.h:388
void makePrecursorSelectionForKnownLCMSMap(const FeatureMap &features, const MSExperiment< InputPeakType > &experiment, MSExperiment< InputPeakType > &ms2, std::set< Int > &charges_set, bool feature_based)
Makes the precursor selection for a given feature map, either feature or scan based.
Definition: OfflinePrecursorIonSelection.h:344
Iterator end()
Definition: MSExperiment.h:157
ConstIterator RTBegin(CoordinateType rt) const
Fast search for spectrum range begin.
Definition: MSExperiment.h:374
CoordinateType minY() const
Accessor for max_ coordinate minimum.
Definition: DIntervalBase.h:247
void getMassRanges(const FeatureMap &features, const MSExperiment< InputPeakType > &experiment, std::vector< std::vector< std::pair< Size, Size > > > &indices)
Calculates the mass ranges for each feature and stores them as indices of the raw data...
Definition: OfflinePrecursorIonSelection.h:156
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
void calculateXICs_(const FeatureMap &features, const std::vector< std::vector< std::pair< Size, Size > > > &mass_ranges, const MSExperiment< InputPeakType > &experiment, const std::set< Int > &charges_set, std::vector< std::vector< std::pair< Size, double > > > &xics)
Calculate the sum of intensities of relevant features for each scan separately.
Definition: OfflinePrecursorIonSelection.h:307
const std::vector< ConvexHull2D > & getConvexHulls() const
Non-mutable access to the convex hulls.
void sortByIntensity(bool reverse=false)
Lexicographically sorts the peaks by their intensity.
Definition: MSSpectrum.h:342
LPWrapper::SOLVER solver_
Definition: OfflinePrecursorIonSelection.h:136
double getRT() const
Definition: MSSpectrum.h:243
void setMSLevel(UInt ms_level)
Sets the MS level.
Definition: MSSpectrum.h:265
An LC-MS feature.
Definition: Feature.h:70
void setLPSolver(LPWrapper::SOLVER solver)
Definition: OfflinePrecursorIonSelection.h:101
void addSpectrum(const MSSpectrum< PeakT > &spectrum)
adds a spectra to the list
Definition: MSExperiment.h:758
void setRT(double rt)
Sets the absolute retention time (is seconds)
Definition: MSSpectrum.h:249
Invalid UInt exception.
Definition: Exception.h:304
void checkMassRanges_(std::vector< std::vector< std::pair< Size, Size > > > &mass_ranges, const MSExperiment< InputPeakType > &experiment)
Eliminates overlapping peaks.
Definition: OfflinePrecursorIonSelection.h:578
In-Memory representation of a mass spectrometry experiment.
Definition: MSExperiment.h:69
std::vector< SpectrumType >::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSExperiment.h:103
Implements different algorithms for precursor ion selection.
Definition: OfflinePrecursorIonSelection.h:62
void setPrecursors(const std::vector< Precursor > &precursors)
sets the precursors
Class for comparison of std::pair using second ONLY e.g. for use with std::sort.
Definition: ComparatorUtils.h:340
bool toBool() const
Conversion to bool.
bool empty() const
Definition: MSExperiment.h:127
void setCharge(Int charge)
Mutable access to the charge.
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:92
PSLPFormulation::IndexTriple IndexTriple
Definition: OfflinePrecursorIonSelection.h:66
int Int
Signed integer type.
Definition: Types.h:96
Implements ILP formulation of precursor selection problems.
Definition: PSLPFormulation.h:54

OpenMS / TOPP release 2.0.0 Documentation generated on Tue Aug 25 2015 05:53:50 using doxygen 1.8.9.1