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

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