OpenMS  3.0.0
SpectraMerger.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-2022.
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: Chris Bielow $
32 // $Authors: Chris Bielow, Andreas Bertsch, Lars Nilse $
33 // --------------------------------------------------------------------------
34 //
35 #pragma once
36 
49 #include <vector>
50 
51 namespace OpenMS
52 {
53 
62  class OPENMS_DLLAPI SpectraMerger :
63  public DefaultParamHandler, public ProgressLogger
64  {
65 
66 protected:
67 
68  /* Determine distance between two spectra
69 
70  Distance is determined as
71 
72  (d_rt/rt_max_ + d_mz/mz_max_) / 2
73 
74  */
76  public DefaultParamHandler
77  {
78 public:
80  DefaultParamHandler("SpectraDistance")
81  {
82  defaults_.setValue("rt_tolerance", 10.0, "Maximal RT distance (in [s]) for two spectra's precursors.");
83  defaults_.setValue("mz_tolerance", 1.0, "Maximal m/z distance (in Da) for two spectra's precursors.");
84  defaultsToParam_(); // calls updateMembers_
85  }
86 
87  void updateMembers_() override
88  {
89  rt_max_ = (double) param_.getValue("rt_tolerance");
90  mz_max_ = (double) param_.getValue("mz_tolerance");
91  }
92 
93  double getSimilarity(const double d_rt, const double d_mz) const
94  {
95  // 1 - distance
96  return 1 - ((d_rt / rt_max_ + d_mz / mz_max_) / 2);
97  }
98 
99  // measure of SIMILARITY (not distance, i.e. 1-distance)!!
100  double operator()(const BaseFeature& first, const BaseFeature& second) const
101  {
102  // get RT distance:
103  double d_rt = fabs(first.getRT() - second.getRT());
104  double d_mz = fabs(first.getMZ() - second.getMZ());
105 
106  if (d_rt > rt_max_ || d_mz > mz_max_)
107  {
108  return 0;
109  }
110 
111  // calculate similarity (0-1):
112  double sim = getSimilarity(d_rt, d_mz);
113 
114  return sim;
115  }
116 
117 protected:
118  double rt_max_;
119  double mz_max_;
120 
121  }; // end of SpectraDistance
122 
123 public:
124 
126  typedef std::map<Size, std::vector<Size> > MergeBlocks;
127 
129  typedef std::map<Size, std::vector<std::pair<Size, double> > > AverageBlocks;
130 
131  // @name Constructors and Destructors
132  // @{
134  SpectraMerger();
135 
137  SpectraMerger(const SpectraMerger& source);
138 
140  SpectraMerger(SpectraMerger&& source) = default;
141 
143  ~SpectraMerger() override;
144  // @}
145 
146  // @name Operators
147  // @{
149  SpectraMerger& operator=(const SpectraMerger& source);
150 
152  SpectraMerger& operator=(SpectraMerger&& source) = default;
153  // @}
154 
155  // @name Merging functions
156  // @{
158  template <typename MapType>
160  {
161  IntList ms_levels = param_.getValue("block_method:ms_levels");
162  Int rt_block_size(param_.getValue("block_method:rt_block_size"));
163  double rt_max_length = (param_.getValue("block_method:rt_max_length"));
164 
165  if (rt_max_length == 0) // no rt restriction set?
166  {
167  rt_max_length = (std::numeric_limits<double>::max)(); // set max rt span to very large value
168  }
169 
170  for (IntList::iterator it_mslevel = ms_levels.begin(); it_mslevel < ms_levels.end(); ++it_mslevel)
171  {
172  MergeBlocks spectra_to_merge;
173  Size idx_block(0);
174  SignedSize block_size_count(rt_block_size + 1);
175  Size idx_spectrum(0);
176  for (typename MapType::const_iterator it1 = exp.begin(); it1 != exp.end(); ++it1)
177  {
178  if (Int(it1->getMSLevel()) == *it_mslevel)
179  {
180  // block full if it contains a maximum number of scans or if maximum rt length spanned
181  if (++block_size_count >= rt_block_size ||
182  exp[idx_spectrum].getRT() - exp[idx_block].getRT() > rt_max_length)
183  {
184  block_size_count = 0;
185  idx_block = idx_spectrum;
186  }
187  else
188  {
189  spectra_to_merge[idx_block].push_back(idx_spectrum);
190  }
191  }
192 
193  ++idx_spectrum;
194  }
195  // check if last block had sacrifice spectra
196  if (block_size_count == 0) //block just got initialized
197  {
198  spectra_to_merge[idx_block] = std::vector<Size>();
199  }
200 
201  // merge spectra, remove all old MS spectra and add new consensus spectra
202  mergeSpectra_(exp, spectra_to_merge, *it_mslevel);
203  }
204 
205  exp.sortSpectra();
206  }
207 
209  template <typename MapType>
211  {
212 
213  // convert spectra's precursors to clusterizable data
214  Size data_size;
215  std::vector<BinaryTreeNode> tree;
216  std::map<Size, Size> index_mapping;
217  // local scope to save memory - we do not need the clustering stuff later
218  {
219  std::vector<BaseFeature> data;
220 
221  for (Size i = 0; i < exp.size(); ++i)
222  {
223  if (exp[i].getMSLevel() != 2)
224  {
225  continue;
226  }
227 
228  // remember which index in distance data ==> experiment index
229  index_mapping[data.size()] = i;
230 
231  // make cluster element
232  BaseFeature bf;
233  bf.setRT(exp[i].getRT());
234  const auto& pcs = exp[i].getPrecursors();
235  // keep the first Precursor
236  if (pcs.empty())
237  {
238  throw Exception::MissingInformation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, String("Scan #") + String(i) + " does not contain any precursor information! Unable to cluster!");
239  }
240  if (pcs.size() > 1)
241  {
242  OPENMS_LOG_WARN << "More than one precursor found. Using first one!" << std::endl;
243  }
244  bf.setMZ(pcs[0].getMZ());
245  data.push_back(bf);
246  }
247  data_size = data.size();
248 
249  SpectraDistance_ llc;
250  llc.setParameters(param_.copy("precursor_method:", true));
251  SingleLinkage sl;
252  DistanceMatrix<float> dist; // will be filled
254 
255  //ch.setThreshold(0.99);
256  // clustering ; threshold is implicitly at 1.0, i.e. distances of 1.0 (== similarity 0) will not be clustered
257  ch.cluster<BaseFeature, SpectraDistance_>(data, llc, sl, tree, dist);
258  }
259 
260  // extract the clusters
261  ClusterAnalyzer ca;
262  std::vector<std::vector<Size> > clusters;
263  // count number of real tree nodes (not the -1 ones):
264  Size node_count = 0;
265  for (Size ii = 0; ii < tree.size(); ++ii)
266  {
267  if (tree[ii].distance >= 1)
268  {
269  tree[ii].distance = -1; // manually set to disconnect, as SingleLinkage does not support it
270  }
271  if (tree[ii].distance != -1)
272  {
273  ++node_count;
274  }
275  }
276  ca.cut(data_size - node_count, tree, clusters);
277 
278  //std::cerr << "Treesize: " << (tree.size()+1) << " #clusters: " << clusters.size() << std::endl;
279  //std::cerr << "tree:\n" << ca.newickTree(tree, true) << "\n";
280 
281  // convert to blocks
282  MergeBlocks spectra_to_merge;
283 
284  for (Size i_outer = 0; i_outer < clusters.size(); ++i_outer)
285  {
286  if (clusters[i_outer].size() <= 1)
287  {
288  continue;
289  }
290  // init block with first cluster element
291  Size cl_index0 = clusters[i_outer][0];
292  spectra_to_merge[index_mapping[cl_index0]] = std::vector<Size>();
293  // add all other elements
294  for (Size i_inner = 1; i_inner < clusters[i_outer].size(); ++i_inner)
295  {
296  spectra_to_merge[index_mapping[cl_index0]].push_back(index_mapping[clusters[i_outer][i_inner]]);
297  }
298  }
299 
300  // do it
301  mergeSpectra_(exp, spectra_to_merge, 2);
302 
303  exp.sortSpectra();
304  }
305 
312  template <typename MapType>
313  void average(MapType& exp, const String& average_type)
314  {
315  // MS level to be averaged
316  int ms_level = param_.getValue("average_gaussian:ms_level");
317  if (average_type == "tophat")
318  {
319  ms_level = param_.getValue("average_tophat:ms_level");
320  }
321 
322  // spectrum type (profile, centroid or automatic)
323  std::string spectrum_type = param_.getValue("average_gaussian:spectrum_type");
324  if (average_type == "tophat")
325  {
326  spectrum_type = std::string(param_.getValue("average_tophat:spectrum_type"));
327  }
328 
329  // parameters for Gaussian averaging
330  double fwhm(param_.getValue("average_gaussian:rt_FWHM"));
331  double factor = -4 * log(2.0) / (fwhm * fwhm); // numerical factor within Gaussian
332  double cutoff(param_.getValue("average_gaussian:cutoff"));
333 
334  // parameters for Top-Hat averaging
335  bool unit(param_.getValue("average_tophat:rt_unit") == "scans"); // true if RT unit is 'scans', false if RT unit is 'seconds'
336  double range(param_.getValue("average_tophat:rt_range")); // range of spectra to be averaged over
337  double range_seconds = range / 2; // max. +/- <range_seconds> seconds from master spectrum
338  int range_scans = static_cast<int>(range); // in case of unit scans, the param is used as integer
339  if ((range_scans % 2) == 0)
340  {
341  ++range_scans;
342  }
343  range_scans = (range_scans - 1) / 2; // max. +/- <range_scans> scans from master spectrum
344 
345  AverageBlocks spectra_to_average_over;
346 
347  // loop over RT
348  int n(0); // spectrum index
349  for (typename MapType::const_iterator it_rt = exp.begin(); it_rt != exp.end(); ++it_rt)
350  {
351  if (Int(it_rt->getMSLevel()) == ms_level)
352  {
353  int m; // spectrum index
354  int steps;
355  bool terminate_now;
356  typename MapType::const_iterator it_rt_2;
357 
358  // go forward (start at next downstream spectrum; the current spectrum will be covered when looking backwards)
359  steps = 0;
360  m = n + 1;
361  it_rt_2 = it_rt + 1;
362  terminate_now = false;
363  while (it_rt_2 != exp.end() && !terminate_now)
364  {
365  if (Int(it_rt_2->getMSLevel()) == ms_level)
366  {
367  double weight = 1;
368  if (average_type == "gaussian")
369  {
370  //factor * (rt_2 -rt)^2
371  double base = it_rt_2->getRT() - it_rt->getRT();
372  weight = std::exp(factor * base * base);
373  }
374  std::pair<Size, double> p(m, weight);
375  spectra_to_average_over[n].push_back(p);
376  ++steps;
377  }
378  if (average_type == "gaussian")
379  {
380  // Gaussian
381  double base = it_rt_2->getRT() - it_rt->getRT();
382  terminate_now = std::exp(factor * base * base) < cutoff;
383  }
384  else if (unit)
385  {
386  // Top-Hat with RT unit = scans
387  terminate_now = (steps > range_scans);
388  }
389  else
390  {
391  // Top-Hat with RT unit = seconds
392  terminate_now = (std::abs(it_rt_2->getRT() - it_rt->getRT()) > range_seconds);
393  }
394  ++m;
395  ++it_rt_2;
396  }
397 
398  // go backward
399  steps = 0;
400  m = n;
401  it_rt_2 = it_rt;
402  terminate_now = false;
403  while (it_rt_2 != exp.begin() && !terminate_now)
404  {
405  if (Int(it_rt_2->getMSLevel()) == ms_level)
406  {
407  double weight = 1;
408  if (average_type == "gaussian")
409  {
410  double base = it_rt_2->getRT() - it_rt->getRT();
411  weight = std::exp(factor * base * base);
412  }
413  std::pair<Size, double> p (m, weight);
414  spectra_to_average_over[n].push_back(p);
415  ++steps;
416  }
417  if (average_type == "gaussian")
418  {
419  // Gaussian
420  double base = it_rt_2->getRT() - it_rt->getRT();
421  terminate_now = std::exp(factor * base * base) < cutoff;
422  }
423  else if (unit)
424  {
425  // Top-Hat with RT unit = scans
426  terminate_now = (steps > range_scans);
427  }
428  else
429  {
430  // Top-Hat with RT unit = seconds
431  terminate_now = (std::abs(it_rt_2->getRT() - it_rt->getRT()) > range_seconds);
432  }
433  --m;
434  --it_rt_2;
435  }
436 
437  }
438  ++n;
439  }
440 
441  // normalize weights
442  for (AverageBlocks::iterator it = spectra_to_average_over.begin(); it != spectra_to_average_over.end(); ++it)
443  {
444  double sum(0.0);
445  for (const auto& weight: it->second)
446  {
447  sum += weight.second;
448  }
449 
450  for (auto& weight: it->second)
451  {
452  weight.second /= sum;
453  }
454  }
455 
456  // determine type of spectral data (profile or centroided)
458  if (spectrum_type == "automatic")
459  {
460  Size idx = spectra_to_average_over.begin()->first; // index of first spectrum to be averaged
461  type = exp[idx].getType(true);
462  }
463  else if (spectrum_type == "profile")
464  {
466  }
467  else if (spectrum_type == "centroid")
468  {
470  }
471  else
472  {
473  throw Exception::InvalidParameter(__FILE__,__LINE__,OPENMS_PRETTY_FUNCTION, "Spectrum type has to be one of automatic, profile or centroid.");
474  }
475 
476  // generate new spectra
477  if (type == SpectrumSettings::CENTROID)
478  {
479  averageCentroidSpectra_(exp, spectra_to_average_over, ms_level);
480  }
481  else
482  {
483  averageProfileSpectra_(exp, spectra_to_average_over, ms_level);
484  }
485 
486  exp.sortSpectra();
487  }
488 
489  // @}
490 
491 protected:
492 
503  template <typename MapType>
504  void mergeSpectra_(MapType& exp, const MergeBlocks& spectra_to_merge, const UInt ms_level)
505  {
506  double mz_binning_width(param_.getValue("mz_binning_width"));
507  std::string mz_binning_unit(param_.getValue("mz_binning_width_unit"));
508 
509  // merge spectra
510  MapType merged_spectra;
511 
512  std::map<Size, Size> cluster_sizes;
513  std::set<Size> merged_indices;
514 
515  // set up alignment
516  SpectrumAlignment sas;
517  Param p;
518  p.setValue("tolerance", mz_binning_width);
519  if (!(mz_binning_unit == "Da" || mz_binning_unit == "ppm"))
520  {
521  throw Exception::IllegalSelfOperation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION); // sanity check
522  }
523 
524  p.setValue("is_relative_tolerance", mz_binning_unit == "Da" ? "false" : "true");
525  sas.setParameters(p);
526  std::vector<std::pair<Size, Size> > alignment;
527 
528  Size count_peaks_aligned(0);
529  Size count_peaks_overall(0);
530 
531  // each BLOCK
532  for (auto it = spectra_to_merge.begin(); it != spectra_to_merge.end(); ++it)
533  {
534  ++cluster_sizes[it->second.size() + 1]; // for stats
535 
536  typename MapType::SpectrumType consensus_spec = exp[it->first];
537  consensus_spec.setMSLevel(ms_level);
538 
539  //consensus_spec.unify(exp[it->first]); // append meta info
540  merged_indices.insert(it->first);
541 
542  //typename MapType::SpectrumType all_peaks = exp[it->first];
543  double rt_average = consensus_spec.getRT();
544  double precursor_mz_average = 0.0;
545  Size precursor_count(0);
546  if (!consensus_spec.getPrecursors().empty())
547  {
548  precursor_mz_average = consensus_spec.getPrecursors()[0].getMZ();
549  ++precursor_count;
550  }
551 
552  count_peaks_overall += consensus_spec.size();
553 
554  // block elements
555  for (auto sit = it->second.begin(); sit != it->second.end(); ++sit)
556  {
557  consensus_spec.unify(exp[*sit]); // append meta info
558  merged_indices.insert(*sit);
559 
560  rt_average += exp[*sit].getRT();
561  if (ms_level >= 2 && exp[*sit].getPrecursors().size() > 0)
562  {
563  precursor_mz_average += exp[*sit].getPrecursors()[0].getMZ();
564  ++precursor_count;
565  }
566 
567  // merge data points
568  sas.getSpectrumAlignment(alignment, consensus_spec, exp[*sit]);
569  //std::cerr << "alignment of " << it->first << " with " << *sit << " yielded " << alignment.size() << " common peaks!\n";
570  count_peaks_aligned += alignment.size();
571  count_peaks_overall += exp[*sit].size();
572 
573  Size align_index(0);
574  Size spec_b_index(0);
575 
576  // sanity check for number of peaks
577  Size spec_a = consensus_spec.size(), spec_b = exp[*sit].size(), align_size = alignment.size();
578  for (auto pit = exp[*sit].begin(); pit != exp[*sit].end(); ++pit)
579  {
580  if (alignment.empty() || alignment[align_index].second != spec_b_index)
581  // ... add unaligned peak
582  {
583  consensus_spec.push_back(*pit);
584  }
585  // or add aligned peak height to ALL corresponding existing peaks
586  else
587  {
588  Size counter(0);
589  Size copy_of_align_index(align_index);
590 
591  while (!alignment.empty() &&
592  copy_of_align_index < alignment.size() &&
593  alignment[copy_of_align_index].second == spec_b_index)
594  {
595  ++copy_of_align_index;
596  ++counter;
597  } // Count the number of peaks in a which correspond to a single b peak.
598 
599  while (!alignment.empty() &&
600  align_index < alignment.size() &&
601  alignment[align_index].second == spec_b_index)
602  {
603  consensus_spec[alignment[align_index].first].setIntensity(consensus_spec[alignment[align_index].first].getIntensity() +
604  (pit->getIntensity() / (double)counter)); // add the intensity divided by the number of peaks
605  ++align_index; // this aligned peak was explained, wait for next aligned peak ...
606  if (align_index == alignment.size())
607  {
608  alignment.clear(); // end reached -> avoid going into this block again
609  }
610  }
611  align_size = align_size + 1 - counter; //Decrease align_size by number of
612  }
613  ++spec_b_index;
614  }
615  consensus_spec.sortByPosition(); // sort, otherwise next alignment will fail
616  if (spec_a + spec_b - align_size != consensus_spec.size())
617  {
618  OPENMS_LOG_WARN << "wrong number of features after merge. Expected: " << spec_a + spec_b - align_size << " got: " << consensus_spec.size() << "\n";
619  }
620  }
621  rt_average /= it->second.size() + 1;
622  consensus_spec.setRT(rt_average);
623 
624  if (ms_level >= 2)
625  {
626  if (precursor_count)
627  {
628  precursor_mz_average /= precursor_count;
629  }
630  auto& pcs = consensus_spec.getPrecursors();
631  //if (pcs.size()>1) OPENMS_LOG_WARN << "Removing excessive precursors - leaving only one per MS2 spectrum.\n";
632  pcs.resize(1);
633  pcs[0].setMZ(precursor_mz_average);
634  }
635 
636  if (consensus_spec.empty())
637  {
638  continue;
639  }
640  else
641  {
642  merged_spectra.addSpectrum(std::move(consensus_spec));
643  }
644  }
645 
646  OPENMS_LOG_INFO << "Cluster sizes:\n";
647  for (const auto& cl_size : cluster_sizes)
648  {
649  OPENMS_LOG_INFO << " size " << cl_size.first << ": " << cl_size.second << "x\n";
650  }
651 
652  char buffer[200];
653  sprintf(buffer, "%d/%d (%.2f %%) of blocked spectra", (int)count_peaks_aligned,
654  (int)count_peaks_overall, float(count_peaks_aligned) / float(count_peaks_overall) * 100.);
655  OPENMS_LOG_INFO << "Number of merged peaks: " << String(buffer) << "\n";
656 
657  // remove all spectra that were within a cluster
658  typename MapType::SpectrumType empty_spec;
659  MapType exp_tmp;
660  for (Size i = 0; i < exp.size(); ++i)
661  {
662  if (merged_indices.count(i) == 0) // save unclustered ones
663  {
664  exp_tmp.addSpectrum(exp[i]);
665  exp[i] = empty_spec;
666  }
667  }
668 
669  //typedef std::vector<typename MapType::SpectrumType> Base;
670  //exp.Base::operator=(exp_tmp);
671  //Meta_Data will not be cleared
672  exp.clear(false);
673  exp.getSpectra().insert(exp.end(), std::make_move_iterator(exp_tmp.begin()),
674  std::make_move_iterator(exp_tmp.end()));
675 
676  // exp.erase(remove_if(exp.begin(), exp.end(), InMSLevelRange<typename MapType::SpectrumType>(ListUtils::create<int>(String(ms_level)), false)), exp.end());
677 
678  // ... and add consensus spectra
679  exp.getSpectra().insert(exp.end(), std::make_move_iterator(merged_spectra.begin()),
680  std::make_move_iterator(merged_spectra.end()));
681 
682  }
683 
704  template <typename MapType>
705  void averageProfileSpectra_(MapType& exp, const AverageBlocks& spectra_to_average_over, const UInt ms_level)
706  {
707  MapType exp_tmp; // temporary experiment for averaged spectra
708 
709  double mz_binning_width(param_.getValue("mz_binning_width"));
710  std::string mz_binning_unit(param_.getValue("mz_binning_width_unit"));
711 
712  unsigned progress = 0;
713  std::stringstream progress_message;
714  progress_message << "averaging profile spectra of MS level " << ms_level;
715  startProgress(0, spectra_to_average_over.size(), progress_message.str());
716 
717  // loop over blocks
718  for (AverageBlocks::const_iterator it = spectra_to_average_over.begin(); it != spectra_to_average_over.end(); ++it)
719  {
720  setProgress(++progress);
721 
722  // loop over spectra in blocks
723  std::vector<double> mz_positions_all; // m/z positions from all spectra
724  for (const auto& spec : it->second)
725  {
726  // loop over m/z positions
727  for (typename MapType::SpectrumType::ConstIterator it_mz = exp[spec.first].begin(); it_mz < exp[spec.first].end(); ++it_mz)
728  {
729  mz_positions_all.push_back(it_mz->getMZ());
730  }
731  }
732 
733  sort(mz_positions_all.begin(), mz_positions_all.end());
734 
735  std::vector<double> mz_positions; // positions at which the averaged spectrum should be evaluated
736  std::vector<double> intensities;
737  double last_mz = std::numeric_limits<double>::min(); // last m/z position pushed through from mz_position to mz_position_2
738  double delta_mz(mz_binning_width); // for m/z unit Da
739  for (const auto mz_pos : mz_positions_all)
740  {
741  if (mz_binning_unit == "ppm")
742  {
743  delta_mz = mz_binning_width * mz_pos / 1000000;
744  }
745 
746  if ((mz_pos - last_mz) > delta_mz)
747  {
748  mz_positions.push_back(mz_pos);
749  intensities.push_back(0.0);
750  last_mz = mz_pos;
751  }
752  }
753 
754  // loop over spectra in blocks
755  for (const auto& spec : it->second)
756  {
757  SplineInterpolatedPeaks spline(exp[spec.first]);
759 
760  // loop over m/z positions
761  for (Size i = spline.getPosMin(); i < mz_positions.size(); ++i)
762  {
763  if ((spline.getPosMin() < mz_positions[i]) && (mz_positions[i] < spline.getPosMax()))
764  {
765  intensities[i] += nav.eval(mz_positions[i]) * (spec.second); // spline-interpolated intensity * weight
766  }
767  }
768  }
769 
770  // update spectrum
771  typename MapType::SpectrumType average_spec = exp[it->first];
772  average_spec.clear(false); // Precursors are part of the meta data, which are not deleted.
773  //average_spec.setMSLevel(ms_level);
774 
775  // refill spectrum
776  for (Size i = 0; i < mz_positions.size(); ++i)
777  {
778  typename MapType::PeakType peak;
779  peak.setMZ(mz_positions[i]);
780  peak.setIntensity(intensities[i]);
781  average_spec.push_back(peak);
782  }
783 
784  // store spectrum temporarily
785  exp_tmp.addSpectrum(std::move(average_spec));
786  }
787 
788  endProgress();
789 
790  // loop over blocks
791  int n(0);
792  //typename MapType::SpectrumType empty_spec;
793  for (AverageBlocks::const_iterator it = spectra_to_average_over.begin(); it != spectra_to_average_over.end(); ++it)
794  {
795  exp[it->first] = exp_tmp[n];
796  //exp_tmp[n] = empty_spec;
797  ++n;
798  }
799  }
800 
816  template <typename MapType>
817  void averageCentroidSpectra_(MapType& exp, const AverageBlocks& spectra_to_average_over, const UInt ms_level)
818  {
819  MapType exp_tmp; // temporary experiment for averaged spectra
820 
821  double mz_binning_width(param_.getValue("mz_binning_width"));
822  std::string mz_binning_unit(param_.getValue("mz_binning_width_unit"));
823 
824  unsigned progress = 0;
825  ProgressLogger logger;
826  std::stringstream progress_message;
827  progress_message << "averaging centroid spectra of MS level " << ms_level;
828  logger.startProgress(0, spectra_to_average_over.size(), progress_message.str());
829 
830  // loop over blocks
831  for (AverageBlocks::const_iterator it = spectra_to_average_over.begin(); it != spectra_to_average_over.end(); ++it)
832  {
833  logger.setProgress(++progress);
834 
835  // collect peaks from all spectra
836  // loop over spectra in blocks
837  std::vector<std::pair<double, double> > mz_intensity_all; // m/z positions and peak intensities from all spectra
838  for (const auto& weightedMZ: it->second)
839  {
840  // loop over m/z positions
841  for (typename MapType::SpectrumType::ConstIterator it_mz = exp[weightedMZ.first].begin(); it_mz < exp[weightedMZ.first].end(); ++it_mz)
842  {
843  std::pair<double, double> mz_intensity(it_mz->getMZ(), (it_mz->getIntensity() * weightedMZ.second)); // m/z, intensity * weight
844  mz_intensity_all.push_back(mz_intensity);
845  }
846  }
847 
848  sort(mz_intensity_all.begin(), mz_intensity_all.end());
849 
850  // generate new spectrum
851  std::vector<double> mz_new;
852  std::vector<double> intensity_new;
853  double last_mz = std::numeric_limits<double>::min();
854  double delta_mz = mz_binning_width;
855  double sum_mz(0);
856  double sum_intensity(0);
857  Size count(0);
858  for (const auto& mz_pos : mz_intensity_all)
859  {
860  if (mz_binning_unit == "ppm")
861  {
862  delta_mz = mz_binning_width * (mz_pos.first) / 1000000;
863  }
864 
865  if (((mz_pos.first - last_mz) > delta_mz) && (count > 0))
866  {
867  mz_new.push_back(sum_mz / count);
868  intensity_new.push_back(sum_intensity); // intensities already weighted
869 
870  sum_mz = 0;
871  sum_intensity = 0;
872 
873  last_mz = mz_pos.first;
874  count = 0;
875  }
876 
877  sum_mz += mz_pos.first;
878  sum_intensity += mz_pos.second;
879  ++count;
880  }
881  if (count > 0)
882  {
883  mz_new.push_back(sum_mz / count);
884  intensity_new.push_back(sum_intensity); // intensities already weighted
885  }
886 
887  // update spectrum
888  typename MapType::SpectrumType average_spec = exp[it->first];
889  average_spec.clear(false); // Precursors are part of the meta data, which are not deleted.
890  //average_spec.setMSLevel(ms_level);
891 
892  // refill spectrum
893  for (Size i = 0; i < mz_new.size(); ++i)
894  {
895  typename MapType::PeakType peak;
896  peak.setMZ(mz_new[i]);
897  peak.setIntensity(intensity_new[i]);
898  average_spec.push_back(peak);
899  }
900 
901  // store spectrum temporarily
902  exp_tmp.addSpectrum(std::move(average_spec));
903  }
904 
905  logger.endProgress();
906 
907  // loop over blocks
908  int n(0);
909  for (const auto& spectral_index : spectra_to_average_over)
910  {
911  exp[spectral_index.first] = std::move(exp_tmp[n]);
912  ++n;
913  }
914  }
915  };
916 }
LogStream.h
DefaultParamHandler.h
OpenMS::SpectrumAlignment::getSpectrumAlignment
void getSpectrumAlignment(std::vector< std::pair< Size, Size > > &alignment, const SpectrumType1 &s1, const SpectrumType2 &s2) const
Definition: SpectrumAlignment.h:88
OpenMS::SpectrumSettings::CENTROID
centroid data or stick data
Definition: SpectrumSettings.h:73
OpenMS::ProgressLogger::setProgress
void setProgress(SignedSize value) const
Sets the current progress.
OpenMS::TOPPBase
Base class for TOPP applications.
Definition: TOPPBase.h:147
FileHandler.h
OpenMS::Math::sum
static double sum(IteratorType begin, IteratorType end)
Calculates the sum of a range of values.
Definition: StatisticFunctions.h:107
double
OpenMS::Param::setValue
void setValue(const std::string &key, const ParamValue &value, const std::string &description="", const std::vector< std::string > &tags=std::vector< std::string >())
Sets a value.
OpenMS::SpectraMerger::SpectraDistance_::getSimilarity
double getSimilarity(const double d_rt, const double d_mz) const
Definition: SpectraMerger.h:93
OpenMS::SpectraMerger::SpectraDistance_::rt_max_
double rt_max_
Definition: SpectraMerger.h:118
OpenMS::SpectraMerger
Merges blocks of MS or MS2 spectra.
Definition: SpectraMerger.h:62
OpenMS::SpectraMerger::average
void average(MapType &exp, const String &average_type)
average over neighbouring spectra
Definition: SpectraMerger.h:313
OpenMS::MSExperiment::sortSpectra
void sortSpectra(bool sort_mz=true)
Sorts the data points by retention time.
OpenMS::Peak1D::setMZ
void setMZ(CoordinateType mz)
Mutable access to m/z.
Definition: Peak1D.h:119
OpenMS::ClusterHierarchical::cluster
void cluster(std::vector< Data > &data, const SimilarityComparator &comparator, const ClusterFunctor &clusterer, std::vector< BinaryTreeNode > &cluster_tree, DistanceMatrix< float > &original_distance)
Clustering function.
Definition: ClusterHierarchical.h:112
OpenMS::String
A more convenient string class.
Definition: String.h:58
OpenMS::MSExperiment::begin
Iterator begin()
Definition: MSExperiment.h:150
OpenMS::BaseFeature
A basic LC-MS feature.
Definition: BaseFeature.h:58
MzMLFile.h
OpenMS::MSExperiment
In-Memory representation of a mass spectrometry run.
Definition: MSExperiment.h:70
OpenMS::SpectrumSettings::getPrecursors
const std::vector< Precursor > & getPrecursors() const
returns a const reference to the precursors
OpenMS::Size
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
SpectraMerger.h
OpenMS::MSSpectrum::getRT
double getRT() const
ClusterAnalyzer.h
OpenMS::Peak1D::setIntensity
void setIntensity(IntensityType intensity)
Mutable access to the data point intensity (height)
Definition: Peak1D.h:110
OpenMS::ProgressLogger::startProgress
void startProgress(SignedSize begin, SignedSize end, const String &label) const
Initializes the progress display.
OpenMS::SplineInterpolatedPeaks::getNavigator
SplineInterpolatedPeaks::Navigator getNavigator(double scaling=0.7)
returns an iterator for access of spline packages
OpenMS::IntList
std::vector< Int > IntList
Vector of signed integers.
Definition: ListUtils.h:55
OpenMS::SpectraMerger::SpectraDistance_::mz_max_
double mz_max_
Definition: SpectraMerger.h:119
OpenMS::MSExperiment::const_iterator
Base::const_iterator const_iterator
Definition: MSExperiment.h:118
OPENMS_LOG_WARN
#define OPENMS_LOG_WARN
Macro if a warning, a piece of information which should be read by the user, should be logged.
Definition: LogStream.h:460
OpenMS::SpectrumSettings::PROFILE
profile data
Definition: SpectrumSettings.h:74
OpenMS::ProgressLogger::endProgress
void endProgress() const
Ends the progress display.
OpenMS::MSExperiment::size
Size size() const
Definition: MSExperiment.h:120
OpenMS::SpectraMerger::mergeSpectra_
void mergeSpectra_(MapType &exp, const MergeBlocks &spectra_to_merge, const UInt ms_level)
merges blocks of spectra of a certain level
Definition: SpectraMerger.h:504
OpenMS::DefaultParamHandler
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:92
OpenMS::Exception::InvalidParameter
Exception indicating that an invalid parameter was handed over to an algorithm.
Definition: Exception.h:339
OpenMS::SpectrumAlignment
Aligns the peaks of two sorted spectra Method 1: Using a banded (width via 'tolerance' parameter) ali...
Definition: SpectrumAlignment.h:67
OpenMS::SpectraMerger::mergeSpectraBlockWise
void mergeSpectraBlockWise(MapType &exp)
Definition: SpectraMerger.h:159
OpenMS::Int
int Int
Signed integer type.
Definition: Types.h:102
OpenMS::Peak2D::setRT
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:214
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
OpenMS::MSExperiment::addSpectrum
void addSpectrum(const MSSpectrum &spectrum)
adds a spectrum to the list
OpenMS::SpectraMerger::SpectraDistance_
Definition: SpectraMerger.h:75
OpenMS::SplineInterpolatedPeaks::getPosMax
double getPosMax() const
returns the maximum m/z (or RT) of the spectrum
OpenMS::ProgressLogger
Base class for all classes that want to report their progress.
Definition: ProgressLogger.h:52
ProgressLogger.h
OpenMS::SpectrumSettings::SpectrumType
SpectrumType
Spectrum peak type.
Definition: SpectrumSettings.h:70
OpenMS::ClusterAnalyzer
Bundles analyzing tools for a clustering (given as sequence of BinaryTreeNode's)
Definition: ClusterAnalyzer.h:51
OpenMS::FileHandler
Facilitates file handling by file type recognition.
Definition: FileHandler.h:66
BaseFeature.h
OpenMS::FileHandler::loadExperiment
bool loadExperiment(const String &filename, MSExperiment &exp, FileTypes::Type force_type=FileTypes::UNKNOWN, ProgressLogger::LogType log=ProgressLogger::NONE, const bool rewrite_source_file=true, const bool compute_hash=true)
Loads a file into an MSExperiment.
RangeUtils.h
OpenMS::FileTypes::Type
Type
Actual file types enum.
Definition: FileTypes.h:56
OpenMS::FileHandler::getType
static FileTypes::Type getType(const String &filename)
Tries to determine the file type (by name or content)
OpenMS::SplineInterpolatedPeaks
Data structure for spline interpolation of MS1 spectra and chromatograms.
Definition: SplineInterpolatedPeaks.h:59
OpenMS::ClusterAnalyzer::cut
void cut(const Size cluster_quantity, const std::vector< BinaryTreeNode > &tree, std::vector< std::vector< Size > > &clusters)
Method to calculate a partition resulting from a certain step in clustering given by the number of cl...
OpenMS::Peak2D::setMZ
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:202
ClusterHierarchical.h
OpenMS::SpectraMerger::averageProfileSpectra_
void averageProfileSpectra_(MapType &exp, const AverageBlocks &spectra_to_average_over, const UInt ms_level)
average spectra (profile mode)
Definition: SpectraMerger.h:705
OpenMS::MSSpectrum::clear
void clear(bool clear_meta_data)
Clears all data and meta data.
OpenMS::DefaultParamHandler::setParameters
void setParameters(const Param &param)
Sets the parameters.
OpenMS::DistanceMatrix
A two-dimensional distance matrix, similar to OpenMS::Matrix.
Definition: DistanceMatrix.h:67
OpenMS::MSSpectrum::setMSLevel
void setMSLevel(UInt ms_level)
Sets the MS level.
OpenMS::SplineInterpolatedPeaks::Navigator::eval
double eval(double pos)
returns spline interpolated intensity at this position (fast access since we can start search from la...
OpenMS::DefaultParamHandler::getParameters
const Param & getParameters() const
Non-mutable access to the parameters.
CompleteLinkage.h
OpenMS::SpectraMerger::mergeSpectraPrecursors
void mergeSpectraPrecursors(MapType &exp)
merges spectra with similar precursors (must have MS2 level)
Definition: SpectraMerger.h:210
OpenMS::Peak1D
A 1-dimensional raw data point or peak.
Definition: Peak1D.h:53
SingleLinkage.h
OpenMS::Peak2D::getRT
CoordinateType getRT() const
Returns the RT coordinate (index 0)
Definition: Peak2D.h:208
OpenMS::MSSpectrum::sortByPosition
void sortByPosition()
Lexicographically sorts the peaks by their position.
OpenMS::UInt
unsigned int UInt
Unsigned integer type.
Definition: Types.h:94
main
int main(int argc, const char **argv)
Definition: INIFileEditor.cpp:71
OpenMS::SpectraMerger::AverageBlocks
std::map< Size, std::vector< std::pair< Size, double > > > AverageBlocks
blocks of spectra (master-spectrum index to update to spectra to average over)
Definition: SpectraMerger.h:129
MSExperiment.h
OpenMS::SpectraMerger::SpectraDistance_::operator()
double operator()(const BaseFeature &first, const BaseFeature &second) const
Definition: SpectraMerger.h:100
OpenMS::SignedSize
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
OpenMS::SpectrumSettings::unify
void unify(const SpectrumSettings &rhs)
merge another spectrum setting into this one (data is usually appended, except for spectrum type whic...
OpenMS::MSExperiment::end
Iterator end()
Definition: MSExperiment.h:160
OpenMS::MSExperiment::getSpectra
const std::vector< MSSpectrum > & getSpectra() const
returns the spectrum list
OpenMS::SplineInterpolatedPeaks::Navigator
iterator class for access of spline packages
Definition: SplineInterpolatedPeaks.h:109
OpenMS::SplineInterpolatedPeaks::getPosMin
double getPosMin() const
returns the minimum m/z (or RT) of the spectrum
OpenMS::SingleLinkage
SingleLinkage ClusterMethod.
Definition: SingleLinkage.h:55
SplineInterpolatedPeaks.h
OpenMS::ClusterHierarchical
Hierarchical clustering with generic clustering functions.
Definition: ClusterHierarchical.h:63
OpenMS::Param
Management and storage of parameters / INI files.
Definition: Param.h:69
OpenMS::MSExperiment::clear
void clear(bool clear_meta_data)
Clears all data and meta data.
OpenMS::FileHandler::storeExperiment
void storeExperiment(const String &filename, const MSExperiment &exp, ProgressLogger::LogType log=ProgressLogger::NONE)
Stores an MSExperiment to a file.
OpenMS::MSExperiment::ConstIterator
std::vector< SpectrumType >::const_iterator ConstIterator
Non-mutable iterator.
Definition: MSExperiment.h:106
OpenMS::Exception::MissingInformation
Not all required information provided.
Definition: Exception.h:186
OpenMS::Peak2D::getMZ
CoordinateType getMZ() const
Returns the m/z coordinate (index 1)
Definition: Peak2D.h:196
OPENMS_LOG_INFO
#define OPENMS_LOG_INFO
Macro if a information, e.g. a status should be reported.
Definition: LogStream.h:465
OpenMS::SpectraMerger::MergeBlocks
std::map< Size, std::vector< Size > > MergeBlocks
blocks of spectra (master-spectrum index to sacrifice-spectra(the ones being merged into the master-s...
Definition: SpectraMerger.h:126
SpectrumAlignment.h
OpenMS::SpectraMerger::SpectraDistance_::updateMembers_
void updateMembers_() override
This method is used to update extra member variables at the end of the setParameters() method.
Definition: SpectraMerger.h:87
OpenMS::MSSpectrum
The representation of a 1D spectrum.
Definition: MSSpectrum.h:66
OpenMS::MSSpectrum::setRT
void setRT(double rt)
Sets the absolute retention time (in seconds)
OpenMS::SpectraMerger::averageCentroidSpectra_
void averageCentroidSpectra_(MapType &exp, const AverageBlocks &spectra_to_average_over, const UInt ms_level)
average spectra (centroid mode)
Definition: SpectraMerger.h:817
StandardTypes.h
OpenMS::Exception::IllegalSelfOperation
Illegal self operation exception.
Definition: Exception.h:370
OpenMS::ProgressLogger::setLogType
void setLogType(LogType type) const
Sets the progress log that should be used. The default type is NONE!
TOPPBase.h
OpenMS::SpectraMerger::SpectraDistance_::SpectraDistance_
SpectraDistance_()
Definition: SpectraMerger.h:79