Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
MRMTransitionGroupPicker.h
Go to the documentation of this file.
1 // --------------------------------------------------------------------------
2 // OpenMS -- Open-Source Mass Spectrometry
3 // --------------------------------------------------------------------------
4 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
5 // ETH Zurich, and Freie Universitaet Berlin 2002-2013.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Hannes Roest $
32 // $Authors: Hannes Roest $
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_ANALYSIS_OPENSWATH_MRMTRANSITIONGROUPPICKER_H
36 #define OPENMS_ANALYSIS_OPENSWATH_MRMTRANSITIONGROUPPICKER_H
37 
43 
47 
50 
52 
53 //#define DEBUG_TRANSITIONGROUPPICKER
54 
55 namespace OpenMS
56 {
57 
76  class OPENMS_DLLAPI MRMTransitionGroupPicker :
77  public DefaultParamHandler
78  {
79 
80 public:
81 
82  // this is the type in which we store the chromatograms for this analysis
84 
85 protected:
89  bool use_gauss_;
90 
92 
95 
98 
101 
102  std::vector<RichPeakChromatogram> picked_chroms_;
103  std::vector<RichPeakChromatogram> smoothed_chroms_;
104 
106 
107  template <typename SpectrumT, typename TransitionT>
110  SpectrumT & master_peak_container, int chr_idx, double best_left, double best_right)
111  {
112  const SpectrumT & ref_chromatogram = transition_group.getChromatograms()[chr_idx];
113 
114  // search for begin / end of the reference chromatogram (and add one more point)
115  typename SpectrumT::const_iterator begin = ref_chromatogram.begin();
116  while (begin != ref_chromatogram.end() && begin->getMZ() < best_left) {begin++; }
117  if (begin != ref_chromatogram.begin()) {begin--; }
118 
119  typename SpectrumT::const_iterator end = begin;
120  while (end != ref_chromatogram.end() && end->getMZ() < best_right) {end++; }
121  if (end != ref_chromatogram.end()) {end++; }
122 
123  // resize the master container and set the m/z values to the ones of the master container
124  master_peak_container.resize(distance(begin, end));
125  typename SpectrumT::iterator it = master_peak_container.begin();
126  for (typename SpectrumT::const_iterator chrom_it = begin; chrom_it != end; chrom_it++, it++)
127  {
128  it->setMZ(chrom_it->getMZ());
129  }
130  }
131 
133  template <typename SpectrumT>
134  SpectrumT resampleChromatogram_(const SpectrumT & chromatogram,
135  SpectrumT & master_peak_container, double best_left, double best_right)
136  {
137  // get the start / end point of this chromatogram => go one past
138  // best_left / best_right to make the resampling accurate also at the
139  // edge.
140  typename SpectrumT::const_iterator begin = chromatogram.begin();
141  while (begin != chromatogram.end() && begin->getMZ() < best_left) {begin++; }
142  if (begin != chromatogram.begin()) {begin--; }
143 
144  typename SpectrumT::const_iterator end = begin;
145  while (end != chromatogram.end() && end->getMZ() < best_right) {end++; }
146  if (end != chromatogram.end()) {end++; }
147 
148  SpectrumT resampled_peak_container = master_peak_container; // copy the master container, which contains the RT values
149  LinearResamplerAlign lresampler;
150  lresampler.raster(begin, end, resampled_peak_container.begin(), resampled_peak_container.end());
151 
152 #if DEBUG_TRANSITIONGROUPPICKER
153  {
154  std::cout << "=========================================================================== " << std::endl;
155  double tot;
156  tot = 0;
157  for (typename SpectrumT::const_iterator it = begin; it != end; it++)
158  {
159  std::cout << " before resampl " << *it << std::endl;
160  tot += it->getIntensity();
161  }
162  std::cout << " total " << tot << std::endl;
163 
164  tot = 0;
165  for (typename SpectrumT::iterator it = resampled_peak_container.begin(); it != resampled_peak_container.end(); it++)
166  {
167  std::cout << " resampl " << *it << std::endl;
168  tot += it->getIntensity();
169  }
170  std::cout << " total " << tot << std::endl;
171  std::cout << " resampled size " << resampled_peak_container.size() << std::endl;
172  }
173 #endif
174  return resampled_peak_container;
175  }
176 
178 
180  double calculateBgEstimation_(const RichPeakChromatogram& smoothed_chromat, double best_left, double best_right);
181 
183  void updateMembers_();
184 
186  MRMTransitionGroupPicker& operator=(const MRMTransitionGroupPicker& rhs);
187 
188 public:
189 
193 
197 
201  template <typename SpectrumT, typename TransitionT>
203  {
204  // Pick chromatograms
205  picked_chroms_.clear();
206  for (Size k = 0; k < transition_group.getChromatograms().size(); k++)
207  {
208  RichPeakChromatogram& chromatogram = transition_group.getChromatograms()[k];
209  if (!chromatogram.isSorted()) { chromatogram.sortByPosition(); }
210 
211  // pickChromatogram will return the picked and the smoothed chromatogram
212  RichPeakChromatogram picked_chrom, smoothed_chrom;
213  pickChromatogram(chromatogram, smoothed_chrom, picked_chrom);
214  picked_chrom.sortByIntensity(); // we could do without that
215  picked_chroms_.push_back(picked_chrom);
216  smoothed_chroms_.push_back(smoothed_chrom);
217  }
218 
219  // Find features (peak groups) in this group of transitions.
220  // While there are still peaks left, one will be picked and used to create
221  // a feature. Whenever we run out of peaks, we will get -1 back as index
222  // and terminate.
223  int chr_idx, peak_idx, cnt = 0;
224  while (true)
225  {
226  chr_idx = -1; peak_idx = -1;
227  findLargestPeak(picked_chroms_, chr_idx, peak_idx);
228  if (chr_idx == -1 && peak_idx == -1) break;
229 
230  /*
231  // FEATURE (hroest) check that this left/right do not collide with any already present features -- if so, re-set the left/right
232  double best_left = picked_chroms[chr_idx].getFloatDataArrays()[1][peak_idx];
233  double best_right = picked_chroms[chr_idx].getFloatDataArrays()[2][peak_idx];
234  */
235 
236  // get feature, prevent non-extended zero features to be added
237  MRMFeature mrm_feature = createMRMFeature(transition_group, picked_chroms_, chr_idx, peak_idx);
238  if (mrm_feature.getIntensity() > 0)
239  {
240  transition_group.addFeature(mrm_feature);
241  }
242 
243  cnt++;
244  if ((stop_after_feature_ > 0 && cnt > stop_after_feature_) && mrm_feature.getIntensity() / (double)mrm_feature.getMetaValue("total_xic") < stop_after_intensity_ratio_)
245  {
246  break;
247  }
248  }
249  }
250 
252  // This function will return a smoothed chromatogram and a picked chromatogram
253  void pickChromatogram(const RichPeakChromatogram& chromatogram, RichPeakChromatogram& smoothed_chrom, RichPeakChromatogram& picked_chrom);
254 
256  template <typename SpectrumT, typename TransitionT>
258  std::vector<SpectrumT> & picked_chroms, int & chr_idx, int & peak_idx)
259  {
260  MRMFeature mrmFeature;
261  const double best_left = picked_chroms[chr_idx].getFloatDataArrays()[1][peak_idx];
262  const double best_right = picked_chroms[chr_idx].getFloatDataArrays()[2][peak_idx];
263  const double peak_apex = picked_chroms[chr_idx][peak_idx].getRT();
264 
265  // Remove other, overlapping, picked peaks (in this and other
266  // chromatograms) and then ensure that at least one peak is set to zero
267  // (the currently best peak).
268  remove_overlapping_features(picked_chroms, best_left, best_right);
269  picked_chroms[chr_idx][peak_idx].setIntensity(0.0);
270 
271  // Prepare linear resampling of all the chromatograms, here creating the
272  // empty master_peak_container with the same RT (m/z) values as the reference
273  // chromatogram.
274  SpectrumT master_peak_container;
275  prepareMasterContainer_(transition_group, master_peak_container, chr_idx, best_left, best_right);
276 
277  double total_intensity = 0; double total_peak_apices = 0; double total_xic = 0;
278  for (Size k = 0; k < transition_group.getChromatograms().size(); k++)
279  {
280  const SpectrumT& chromatogram = transition_group.getChromatograms()[k];
281  for (typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
282  {
283  total_xic += it->getIntensity();
284  }
285 
286  // resample the current chromatogram
287  const SpectrumT used_chromatogram = resampleChromatogram_(chromatogram, master_peak_container, best_left, best_right);
288  // const SpectrumT& used_chromatogram = chromatogram; // instead of resampling
289 
290  Feature f;
291  double quality = 0;
292  f.setQuality(0, quality);
293  f.setOverallQuality(quality);
294 
295  ConvexHull2D::PointArrayType hull_points;
296  DoubleReal intensity_sum(0.0), rt_sum(0.0);
297  double peak_apex_int = -1;
298  double peak_apex_dist = std::fabs(used_chromatogram.begin()->getMZ() - peak_apex);
299  // FEATURE : use RTBegin / MZBegin -> for this we need to know whether the template param is a real chromatogram or a spectrum!
300  for (typename SpectrumT::const_iterator it = used_chromatogram.begin(); it != used_chromatogram.end(); it++)
301  {
302  if (it->getMZ() > best_left && it->getMZ() < best_right)
303  {
304  DPosition<2> p;
305  p[0] = it->getMZ();
306  p[1] = it->getIntensity();
307  hull_points.push_back(p);
308  if (std::fabs(it->getMZ() - peak_apex) <= peak_apex_dist)
309  {
310  peak_apex_int = p[1];
311  peak_apex_dist = std::fabs(it->getMZ() - peak_apex);
312  }
313  rt_sum += it->getMZ();
314  intensity_sum += it->getIntensity();
315  }
316  }
317 
318  if (background_subtraction_ != "none")
319  {
320  double background = 0;
321  // we use the smoothed chromatogram here to have a more accurate estimatation of the noise at the flanks of the peak
322  if (background_subtraction_ == "smoothed")
323  {
324  if (smoothed_chroms_.size() <= k)
325  {
326  std::cerr << "Tried to calculate background estimation without any smoothed chromatograms" << std::endl;
327  background = 0;
328  }
329  else
330  {
331  background = calculateBgEstimation_(smoothed_chroms_[k], best_left, best_right);
332  }
333  }
334  else if (background_subtraction_ == "original")
335  {
336  background = calculateBgEstimation_(used_chromatogram, best_left, best_right);
337  }
338  intensity_sum -= background;
339  if (intensity_sum < 0)
340  {
341  std::cerr << "Warning: Intensity was below 0 after background subtraction: " << intensity_sum << ". Setting it to 0." << std::endl;
342  intensity_sum = 0;
343  }
344  }
345 
346  f.setRT(picked_chroms[chr_idx][peak_idx].getMZ());
347  f.setMZ(chromatogram.getMetaValue("product_mz"));
348  f.setIntensity(intensity_sum);
349  ConvexHull2D hull;
350  hull.setHullPoints(hull_points);
351  f.getConvexHulls().push_back(hull);
352  f.setMetaValue("MZ", chromatogram.getMetaValue("product_mz"));
353  f.setMetaValue("native_id", chromatogram.getNativeID());
354  f.setMetaValue("peak_apex_int", peak_apex_int);
355  //f.setMetaValue("leftWidth", best_left);
356  //f.setMetaValue("rightWidth", best_right);
357 
358  total_intensity += intensity_sum;
359  total_peak_apices += peak_apex_int;
360  mrmFeature.addFeature(f, chromatogram.getNativeID()); //map index and feature
361  }
362  mrmFeature.setRT(picked_chroms[chr_idx][peak_idx].getMZ());
363  mrmFeature.setIntensity(total_intensity);
364  mrmFeature.setMetaValue("PeptideRef", transition_group.getTransitionGroupID());
365  mrmFeature.setMetaValue("leftWidth", best_left);
366  mrmFeature.setMetaValue("rightWidth", best_right);
367  mrmFeature.setMetaValue("total_xic", total_xic);
368  mrmFeature.setMetaValue("peak_apices_sum", total_peak_apices);
369 
370  return mrmFeature;
371  }
372 
373  // maybe private, but we have tests
374 
376  template <typename SpectrumT>
377  void remove_overlapping_features(std::vector<SpectrumT> & picked_chroms, double best_left, double best_right)
378  {
379  // delete all seeds that lie within the current seed
380  for (Size k = 0; k < picked_chroms.size(); k++)
381  {
382  for (Size i = 0; i < picked_chroms[k].size(); i++)
383  {
384  if (picked_chroms[k][i].getMZ() >= best_left && picked_chroms[k][i].getMZ() <= best_right)
385  {
386  picked_chroms[k][i].setIntensity(0.0);
387  }
388  }
389  }
390 
391  // delete all seeds that overlap within the current seed
392  for (Size k = 0; k < picked_chroms.size(); k++)
393  {
394  for (Size i = 0; i < picked_chroms[k].size(); i++)
395  {
396  double left = picked_chroms[k].getFloatDataArrays()[1][i];
397  double right = picked_chroms[k].getFloatDataArrays()[2][i];
398  if ((left >= best_left && left <= best_right)
399  || (right >= best_left && right <= best_right))
400  {
401  picked_chroms[k][i].setIntensity(0.0);
402  }
403  }
404  }
405  }
406 
408  void findLargestPeak(std::vector<RichPeakChromatogram>& picked_chroms, int& chr_idx, int& peak_idx);
409  };
410 }
411 
412 #endif
const double k
const DataValue & getMetaValue(const String &name) const
returns the value corresponding to a string
void raster(SpecT< PeakType > &spectrum)
Applies the resampling algorithm to an MSSpectrum.
Definition: LinearResamplerAlign.h:69
void setMetaValue(const String &name, const DataValue &value)
sets the DataValue corresponding to a name
std::vector< RichPeakChromatogram > smoothed_chroms_
Definition: MRMTransitionGroupPicker.h:103
String background_subtraction_
Definition: MRMTransitionGroupPicker.h:91
A more convenient string class.
Definition: String.h:56
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:197
void sortByPosition()
Lexicographically sorts the peaks by their position.
Definition: MSSpectrum.h:391
IntensityType getIntensity() const
Definition: Peak2D.h:161
const std::vector< SpectrumType > & getChromatograms() const
Definition: MRMTransitionGroup.h:148
std::vector< PointType > PointArrayType
Definition: ConvexHull2D.h:79
bool isSorted() const
Checks if all peaks are sorted with respect to ascending m/z.
Definition: MSSpectrum.h:453
std::vector< RichPeakChromatogram > picked_chroms_
Definition: MRMTransitionGroupPicker.h:102
int stop_after_feature_
Definition: MRMTransitionGroupPicker.h:99
A 2-dimensional hull representation in [counter]clockwise direction - depending on axis labelling...
Definition: ConvexHull2D.h:75
DoubleReal gauss_width_
Definition: MRMTransitionGroupPicker.h:88
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:167
UInt sgolay_polynomial_order_
Definition: MRMTransitionGroupPicker.h:87
void addFeature(Feature &feature, const String &key)
Adds an feature from a single chromatogram into the feature.
bool use_gauss_
Definition: MRMTransitionGroupPicker.h:89
DoubleReal sn_win_len_
Definition: MRMTransitionGroupPicker.h:96
The representation of a 1D spectrum.
Definition: MSSpectrum.h:67
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:209
The MRMTransitionGroupPicker finds peaks in chromatograms that belong to the same precursors...
Definition: MRMTransitionGroupPicker.h:76
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 setQuality(Size index, QualityType q)
Set the quality in dimension c.
UInt sgolay_frame_length_
Definition: MRMTransitionGroupPicker.h:86
Linear Resampling of raw data with alignment.
Definition: LinearResamplerAlign.h:58
void setHullPoints(const PointArrayType &points)
accessor for the outer(!) points (no checking is performed if this is actually a convex hull) ...
An LC-MS feature.
Definition: Feature.h:66
void addFeature(MRMFeature &feature)
Definition: MRMTransitionGroup.h:184
void prepareMasterContainer_(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group, SpectrumT &master_peak_container, int chr_idx, double best_left, double best_right)
create an empty master peak container that has the correct mz / RT values set
Definition: MRMTransitionGroupPicker.h:109
void setOverallQuality(QualityType q)
Set the overall quality.
SpectrumT resampleChromatogram_(const SpectrumT &chromatogram, SpectrumT &master_peak_container, double best_left, double best_right)
use the master container from above to resample a chromatogram at those points stored in the master c...
Definition: MRMTransitionGroupPicker.h:134
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
void pickTransitionGroup(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group)
Definition: MRMTransitionGroupPicker.h:202
UInt sn_bin_count_
Definition: MRMTransitionGroupPicker.h:97
DoubleReal stop_after_intensity_ratio_
Definition: MRMTransitionGroupPicker.h:100
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:90
MRMFeature createMRMFeature(MRMTransitionGroup< SpectrumT, TransitionT > &transition_group, std::vector< SpectrumT > &picked_chroms, int &chr_idx, int &peak_idx)
Create feature from a vector of chromatograms and a specified peak.
Definition: MRMTransitionGroupPicker.h:257
void remove_overlapping_features(std::vector< SpectrumT > &picked_chroms, double best_left, double best_right)
Remove overlaping features that are within the current seed feature or overlap with it...
Definition: MRMTransitionGroupPicker.h:377
A multi-chromatogram MRM feature.
Definition: MRMFeature.h:50
MSSpectrum< ChromatogramPeak > RichPeakChromatogram
Definition: MRMTransitionGroupPicker.h:83
const String & getTransitionGroupID() const
Definition: MRMTransitionGroup.h:112
DoubleReal peak_width_
Definition: MRMTransitionGroupPicker.h:93
DoubleReal signal_to_noise_
Definition: MRMTransitionGroupPicker.h:94

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