Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
SimpleExtender.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: Clemens Groepl $
32 // $Authors: $
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_SIMPLEEXTENDER_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_SIMPLEEXTENDER_H
37 
40 
41 #include <queue>
42 #include <iostream>
43 #include <fstream>
44 
45 namespace OpenMS
46 {
47 
78  template <class PeakType, class FeatureType>
80  public FeaFiModule<PeakType, FeatureType>,
81  public FeatureFinderDefs
82  {
83 public:
85 
92 
95  Base(map, features, ff),
97  {
98  this->setName("SimpleExtender");
99 
100  this->defaults_.setValue("dist_mz_up", 6.0, "Maximum high m/z distance of peak in the region/boundary from the seed.");
101  this->defaults_.setMinFloat("dist_mz_up", 0.0);
102  this->defaults_.setValue("dist_mz_down", 2.0, "Maximum low m/z distance of peak in the region/boundary from the seed.");
103  this->defaults_.setMinFloat("dist_mz_down", 0.0);
104  this->defaults_.setValue("dist_rt_up", 5.0, "Maximum high RT distance of peak in the region/boundary from the seed.");
105  this->defaults_.setMinFloat("dist_rt_up", 0.0);
106  this->defaults_.setValue("dist_rt_down", 5.0, "Maximum low RT distance of peak in the region/boundary from the seed.");
107  this->defaults_.setMinFloat("dist_rt_down", 0.0);
108 
109  // priority check is per default switched off
110  // these values were used for the Myoglobin quantification project
111  // DON'T REMOVE THIS
112  this->defaults_.setValue("priority_thr", -0.1, "Minimum priority for data points to be included into the boundary of the feature (default 0.0). The priority of a data point is a function of its intensity and its distance to the last point included into the feature region. Setting this threshold to zero or a very small value is usually a good idea.", StringList::create("advanced"));
113 
114  this->defaults_.setValue("intensity_factor", 0.03, "Influences for intensity (ion count) threshold in the feature extension. We include only raw data points into this region if their intensity is larger than [intensity_factor * (intensity of the seed)].");
115  this->defaults_.setMinFloat("intensity_factor", 0.0);
116  this->defaults_.setMaxFloat("intensity_factor", 1.0);
117 
118  this->defaultsToParam_();
119  }
120 
122  virtual ~SimpleExtender()
123  {
124  }
125 
127  void extend(const ChargedIndexSet & seed_region, ChargedIndexSet & result_region)
128  {
129  // empty region and boundary datastructures
130  result_region.clear();
131  priorities_.clear();
133  boundary_ = std::priority_queue<IndexWithPriority, std::vector<IndexWithPriority>, typename IndexWithPriority::PriorityLess>();
134 
135 #ifdef DEBUG_FEATUREFINDER
136  std::vector<IndexPair> debug_vector;
137 #endif
138  // find maximum of region (seed)
139  CoordinateType max_intensity = 0.0;
140  IndexPair seed;
141 
142  for (IndexSet::const_iterator citer = seed_region.begin(); citer != seed_region.end(); ++citer)
143  {
144  if (this->getPeakIntensity(*citer) > max_intensity)
145  {
146  seed = *citer;
147  max_intensity = this->getPeakIntensity(seed);
148  }
149  }
150 
151  // remember last extracted point (in this case the seed !)
154 
155  // Add peaks received from seeder directly to boundary
156  for (IndexSet::const_iterator citer = seed_region.begin(); citer != seed_region.end(); ++citer)
157  {
158  ProbabilityType priority = computePeakPriority_(*citer);
159  priorities_[*citer] = priority;
160  boundary_.push(IndexWithPriority(*citer, priority));
161  }
162  // pass on charge information
163  result_region.charge = seed_region.charge;
164 
165  // re-compute intensity threshold
166  intensity_threshold_ = (DoubleReal)(this->param_).getValue("intensity_factor") * this->getPeakIntensity(seed);
167 
168 #ifdef DEBUG_FEATUREFINDER
169  std::cout << "\n";
170  std::cout << "Extending from " << this->getPeakRt(seed) << "/" << this->getPeakMz(seed) << std::endl;
171  std::cout << "Intensity of seed " << this->getPeakIntensity(seed);
172  std::cout << " (" << seed.first << "/" << seed.second << ")" << std::endl;
173  std::cout << "Intensity_threshold: " << intensity_threshold_ << std::endl;
174 #endif
175 
176  while (!boundary_.empty())
177  {
178  // remove peak with highest priority
179  const IndexPair current_index = boundary_.top().index;
180  boundary_.pop();
181 
182  // check for corrupt index
183  OPENMS_PRECONDITION(current_index.first < (*this->map_).size(), "Scan index outside of map!");
184  OPENMS_PRECONDITION(current_index.second < (*this->map_)[current_index.first].size(), "Peak index outside of scan!");
185 
186  // remember last extracted peak
187  last_pos_extracted_[Peak2D::RT] = this->getPeakRt(current_index);
188  last_pos_extracted_[Peak2D::MZ] = this->getPeakMz(current_index);
189 
190  // Now we explore the neighbourhood of the current peak. Points in this area are included
191  // into the boundary if their intensity is not too low and they are not too
192  // far away from the seed.
193  // Add position to the current average of positions weighted by intensity
194  running_avg_.add(last_pos_extracted_, this->getPeakIntensity(current_index));
195 
196  // explore neighbourhood of current peak
197  moveMzUp_(current_index);
198  moveMzDown_(current_index);
199  moveRtUp_(current_index);
200  moveRtDown_(current_index);
201 
202  // set peak flags and add to boundary
203  this->ff_->getPeakFlag(current_index) = USED;
204 #ifdef DEBUG_FEATUREFINDER
205  debug_vector.push_back(current_index);
206 #endif
207  result_region.insert(current_index);
208 
209  } // end of while ( !boundary_.empty() )
210 
211 #ifdef DEBUG_FEATUREFINDER
212  std::cout << "Feature region size: " << result_region.size() << std::endl;
213 #endif
214 
215 #ifdef DEBUG_FEATUREFINDER
216  static UInt number = 1;
217  writeDebugFile_(debug_vector, number++);
218  debug_vector.clear();
219 #endif
220 
221  return;
222  } // end of extend
223 
237  {
239  index(i),
240  priority(p)
241  {
242  }
243 
246 
249  {
250  inline bool operator()(const IndexWithPriority & x, const IndexWithPriority & y) const
251  {
252  return x.priority < y.priority;
253  }
254 
255  };
256  };
257 
258 protected:
259 
260  virtual void updateMembers_()
261  {
262  dist_mz_up_ = this->param_.getValue("dist_mz_up");
263  dist_mz_down_ = this->param_.getValue("dist_mz_down");
264  dist_rt_up_ = this->param_.getValue("dist_rt_up");
265  dist_rt_down_ = this->param_.getValue("dist_rt_down");
266  priority_threshold_ = this->param_.getValue("priority_thr");
267  }
268 
270  void writeDebugFile_(const std::vector<IndexPair> & peaks, UInt nr_feat)
271  {
272  String filename = String(nr_feat).fillLeft('0', 4) + "_Extension.dta2d";
273  std::ofstream file(filename.c_str());
274  for (Size i = 0; i < peaks.size(); ++i)
275  {
276  file << this->getPeakRt(peaks[i]) << " " << this->getPeakMz(peaks[i]) << " " << peaks.size() - i << std::endl;
277  }
278  file.close();
279  }
280 
282  bool isTooFarFromCentroid_(const IndexPair & index)
283  {
284  //Corrupt index
285  OPENMS_PRECONDITION(index.first < (*this->map_).size(), "Scan index outside of map!");
286  OPENMS_PRECONDITION(index.second < (*this->map_)[index.first].size(), "Peak index outside of scan!");
287 
288  const DPosition<2> & curr_mean = running_avg_.getPosition();
289 
290  if (this->getPeakMz(index) > curr_mean[Peak2D::MZ] + dist_mz_up_ ||
291  this->getPeakMz(index) < curr_mean[Peak2D::MZ] - dist_mz_down_ ||
292  this->getPeakRt(index) > curr_mean[Peak2D::RT] + dist_rt_up_ ||
293  this->getPeakRt(index) < curr_mean[Peak2D::RT] - dist_rt_down_)
294  {
295  //too far
296  return true;
297  }
298 
299  //close enough
300  return false;
301  }
302 
304  void moveMzUp_(const IndexPair & index)
305  {
306  try
307  {
308  IndexPair tmp = index;
309  while (true)
310  {
311  this->getNextMz(tmp);
312  if (isTooFarFromCentroid_(tmp)) break;
313  checkNeighbour_(tmp);
314  }
315  }
316  catch (NoSuccessor)
317  {
318  }
319  }
320 
322  void moveMzDown_(const IndexPair & index)
323  {
324  try
325  {
326  IndexPair tmp = index;
327  while (true)
328  {
329  this->getPrevMz(tmp);
330  if (isTooFarFromCentroid_(tmp)) break;
331  checkNeighbour_(tmp);
332  }
333  }
334  catch (NoSuccessor)
335  {
336  }
337  }
338 
340  void moveRtUp_(const IndexPair & index)
341  {
342  try
343  {
344  IndexPair tmp = index;
345 
346  while (true)
347  {
348  this->getNextRt(tmp);
349  if (isTooFarFromCentroid_(tmp)) break;
350  checkNeighbour_(tmp);
351  }
352  }
353  catch (NoSuccessor)
354  {
355  }
356  }
357 
359  void moveRtDown_(const IndexPair & index)
360  {
361  try
362  {
363  IndexPair tmp = index;
364  while (true)
365  {
366  this->getPrevRt(tmp);
367  if (isTooFarFromCentroid_(tmp)) break;
368  checkNeighbour_(tmp);
369  }
370  }
371  catch (NoSuccessor)
372  {
373  }
374  }
375 
378  {
379  return (*this->map_)[index.first][index.second].getIntensity();
380  }
381 
383  void checkNeighbour_(const IndexPair & index)
384  {
385  //Corrupt index
386  OPENMS_PRECONDITION(index.first < (*this->map_).size(), "Scan index outside of map!");
387  OPENMS_PRECONDITION(index.second < (*this->map_)[index.first].size(), "Peak index outside of scan!");
388 
389  // skip this point if its intensity is too low
390  if (this->getPeakIntensity(index) <= intensity_threshold_)
391  {
392  return;
393  }
394  if (this->ff_->getPeakFlag(index) == UNUSED)
395  {
396  DoubleReal pr_new = computePeakPriority_(index);
397 
398  if (pr_new > priority_threshold_)
399  {
400  //std::map<IndexPair, DoubleReal>::iterator piter = priorities_.find(index);
401  this->ff_->getPeakFlag(index) = USED;
402  priorities_[index] = pr_new;
403  boundary_.push(IndexWithPriority(index, pr_new));
404  }
405  }
406  }
407 
410 
412  std::map<IndexPair, ProbabilityType> priorities_;
413 
416 
418  std::priority_queue<IndexWithPriority, std::vector<IndexWithPriority>, typename IndexWithPriority::PriorityLess> boundary_;
419 
422 
431 
434 
437 
438 private:
440  SimpleExtender();
445 
446  };
447 }
448 #endif // OPENMS_TRANSFORMATIONS_FEATUREFINDER_SIMPLEEXTENDER_H
void extend(const ChargedIndexSet &seed_region, ChargedIndexSet &result_region)
return next seed
Definition: SimpleExtender.h:127
Implements a module of the FeatureFinder algorithm.
Definition: FeaFiModule.h:157
IsotopeCluster::IndexPair IndexPair
Index to peak consisting of two UInts (scan index / peak index)
Definition: FeatureFinderDefs.h:54
Definition: FeatureFinderDefs.h:63
Param defaults_
Container for default parameters. This member should be filled in the constructor of derived classes!...
Definition: DefaultParamHandler.h:155
void setValue(const String &key, const DataValue &value, const String &description="", const StringList &tags=StringList())
Sets a value.
A more convenient string class.
Definition: String.h:56
Base::IntensityType IntensityType
Intensity of a data point.
Definition: SimpleExtender.h:87
void getNextMz(FeatureFinderDefs::IndexPair &index) const
fills index with the index of next peak in m/z dimension
Definition: FeaFiModule.h:225
CoordinateType getPeakRt(const FeatureFinderDefs::IndexPair &index) const
Returns the retention time of a peak.
Definition: FeaFiModule.h:210
void clear()
Reset everything. (Note that update() will cause a division by zero after that.)
Definition: AveragePosition.h:104
A helper structure to sort indizes by their priority.
Definition: SimpleExtender.h:236
Base::CoordinateType CoordinateType
Coordinates of a point (m/z and rt)
Definition: SimpleExtender.h:89
IntensityType getPeakIntensity(const FeatureFinderDefs::IndexPair &index) const
Returns the intensity of a peak.
Definition: FeaFiModule.h:190
void moveRtUp_(const IndexPair &index)
Extension into positive rt dimension.
Definition: SimpleExtender.h:340
std::map< IndexPair, ProbabilityType > priorities_
Keeps track of peaks already included in the boundary (value is priority of peak) ...
Definition: SimpleExtender.h:412
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:148
Definition: FeatureFinderDefs.h:63
Compares two indizes by priority.
Definition: SimpleExtender.h:248
Retention time dimension id (0 if used as a const int)
Definition: Peak2D.h:76
void getNextRt(FeatureFinderDefs::IndexPair &index)
fills index with the index of the nearest peak in the next scan
Definition: FeaFiModule.h:267
PeakType::IntensityType IntensityType
Input intensity type.
Definition: FeaFiModule.h:168
A container for features.
Definition: FeatureMap.h:111
#define OPENMS_PRECONDITION(condition, message)
Precondition macro.
Definition: Macros.h:107
SimpleExtender()
Not implemented.
Mass-to-charge dimension id (1 if used as a const int)
Definition: Peak2D.h:77
void getPrevMz(FeatureFinderDefs::IndexPair &index) const
fills index with the index of previous peak in m/z dimension
Definition: FeaFiModule.h:246
FeatureFinder * ff_
Pointer to the calling FeatureFinder that is used to access the feature flags and report progress...
Definition: FeaFiModule.h:415
The purpose of this struct is to provide definitions of classes and typedefs which are used throughou...
Definition: FeatureFinderDefs.h:51
const MapType * map_
Input data pointer.
Definition: FeaFiModule.h:411
DPosition< 2 > last_pos_extracted_
Position of last peak extracted from the boundary (used to compute the priority of neighbouring peaks...
Definition: SimpleExtender.h:415
std::priority_queue< IndexWithPriority, std::vector< IndexWithPriority >, typename IndexWithPriority::PriorityLess > boundary_
Represents the boundary of a feature.
Definition: SimpleExtender.h:418
ProbabilityType priority
Definition: SimpleExtender.h:245
SimpleExtender & operator=(const SimpleExtender &)
Not implemented.
CoordinateType getPeakMz(const FeatureFinderDefs::IndexPair &index) const
Returns the m/z of a peak.
Definition: FeaFiModule.h:200
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
ChargedIndexSet region_
charged index set
Definition: SimpleExtender.h:436
void getPrevRt(FeatureFinderDefs::IndexPair &index)
fills index with the index of the nearest peak in the previous scan
Definition: FeaFiModule.h:323
IndexPair index
Definition: SimpleExtender.h:244
PositionType const & getPosition() const
Returns the current average position.
Definition: AveragePosition.h:92
FeaFiModule< PeakType, FeatureType > Base
Definition: SimpleExtender.h:84
void writeDebugFile_(const std::vector< IndexPair > &peaks, UInt nr_feat)
write DTA2D debug file for the feature with index nr_feat
Definition: SimpleExtender.h:270
CoordinateType dist_rt_down_
Maximum distance to seed in negative retention time.
Definition: SimpleExtender.h:430
ProbabilityType computePeakPriority_(const IndexPair &index)
Computes the priority of a peak as function of intensity and distance from seed.
Definition: SimpleExtender.h:377
Simple feature extension algorithm.
Definition: SimpleExtender.h:79
static StringList create(const String &list, const char splitter= ',')
Returns a list that is created by splitting the given (comma-separated) string (String are not trimme...
CoordinateType dist_rt_up_
Maximum distance to seed in positive retention time.
Definition: SimpleExtender.h:428
Exception that is thrown if a method an invalid IndexPair is given.
Definition: FeatureFinderDefs.h:66
String & fillLeft(char c, UInt size)
Adds c on the left side until the size of the string is size.
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
CoordinateType dist_mz_down_
Maximum distance to seed in negative m/z.
Definition: SimpleExtender.h:426
void add(PositionType position, CoordinateType const weight=1)
Add a position.
Definition: AveragePosition.h:113
IndexWithPriority(const FeatureFinderDefs::IndexPair &i, DoubleReal p)
Definition: SimpleExtender.h:238
virtual void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
Definition: SimpleExtender.h:260
SimpleExtender(const MSExperiment< PeakType > *map, FeatureMap< FeatureType > *features, FeatureFinder *ff)
Constructor.
Definition: SimpleExtender.h:94
bool operator()(const IndexWithPriority &x, const IndexWithPriority &y) const
Definition: SimpleExtender.h:250
index set with associated charge estimate
Definition: IsotopeCluster.h:53
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
void moveMzUp_(const IndexPair &index)
Extends the seed into positive m/z direction.
Definition: SimpleExtender.h:304
void moveMzDown_(const IndexPair &index)
Extends the seed into negative m/z direction.
Definition: SimpleExtender.h:322
void checkNeighbour_(const IndexPair &index)
Checks the neighbours of the current for insertion into the boundary.
Definition: SimpleExtender.h:383
bool isTooFarFromCentroid_(const IndexPair &index)
Checks if the current peak is too far from the centroid.
Definition: SimpleExtender.h:282
Int charge
charge estimate (convention: zero means &quot;no charge estimate&quot;)
Definition: IsotopeCluster.h:62
CoordinateType dist_mz_up_
Maximum distance to seed in positive m/z.
Definition: SimpleExtender.h:424
The main feature finder class.
Definition: FeatureFinder.h:57
void moveRtDown_(const IndexPair &index)
Extends the seed into negative retention time direction.
Definition: SimpleExtender.h:359
void setName(const String &name)
Mutable access to the name.
DoubleReal ProbabilityType
Priority of a point (see below)
Definition: SimpleExtender.h:91
const Flag & getPeakFlag(const IndexPair &index) const
Returns a non-mutable reference to a peak flag.
Definition: FeatureFinder.h:91
Math::AveragePosition< 2 > running_avg_
keeps an running average of the peak coordinates weighted by the intensities
Definition: SimpleExtender.h:409
void setMinFloat(const String &key, DoubleReal min)
Sets the minimum value for the floating point or floating point list parameter key.
ProbabilityType priority_threshold_
Minium priority for points in the feature region (priority is function of intensity and distance to s...
Definition: SimpleExtender.h:433
virtual ~SimpleExtender()
destructor
Definition: SimpleExtender.h:122
void defaultsToParam_()
Updates the parameters after the defaults have been set in the constructor.
void setMaxFloat(const String &key, DoubleReal max)
Sets the maximum value for the floating point or floating point list parameter key.
IntensityType intensity_threshold_
Mininum intensity of a boundary point. Calculated from &#39;intensity_factor&#39; and the seed intensity...
Definition: SimpleExtender.h:421
double DoubleReal
Double-precision real type.
Definition: Types.h:118

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