PeakTypeEstimator.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2018.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Chris Bielow $
// $Authors: Chris Bielow $
// --------------------------------------------------------------------------

#pragma once

#include <OpenMS/METADATA/SpectrumSettings.h>

#include <cmath>
#include <numeric>

namespace OpenMS
{
  class OPENMS_DLLAPI PeakTypeEstimator
  {
public:
    template <typename PeakConstIterator>
    static SpectrumSettings::SpectrumType estimateType(const PeakConstIterator& begin, const PeakConstIterator& end)
    {
      typedef typename PeakConstIterator::value_type PeakT;
      // abort if there are less than 5 peak in the iterator range
      if (end - begin < 5)
      {
        return SpectrumSettings::UNKNOWN;
      }

      const int max_peaks = 5; // maximal number of peaks we are looing at
      int profile_evidence = 0; // number of peaks found to be profile
      int centroid_evidence = 0; // number of peaks found to be centroided

      // copy data, since we need to modify
      std::vector<PeakT> data(begin, end);
      // total intensity of spectrum
      double total_int = std::accumulate(begin, end, 0.0, [](double int_, const PeakT& p) { return int_ + p.getIntensity(); } );
      double explained_int = 0;
      // get the 5 highest peaks
      for (int i = 0; i < max_peaks; ++i)
      {
        // stop if we explained +50% of all intensity 
        // (due to danger of interpreting noise - usually wrongly classified as centroided data)
        if (explained_int > 0.5 * total_int) break;
        
        double int_max = 0;
        Size idx = std::numeric_limits<Size>::max();
        // find highest peak position
        for (Size i = 0; i < data.size(); ++i)
        {
          if (data[i].getIntensity() > int_max)
          {
            int_max = data[i].getIntensity();
            idx = i;
          }
        } 
        // no more peaks
        if (idx == std::numeric_limits<Size>::max()) break;

        // check left and right peak shoulders and count number of sample points
        typedef typename std::vector<PeakT>::iterator PeakIterator; // non-const version, since we need to modify the peaks
        PeakIterator it_max = data.begin() + idx;
        PeakIterator it = it_max;
        double int_last = int_max;
        while (it != data.begin() 
               && it->getIntensity() <= int_last        // at most 100% of last sample point
               && it->getIntensity() > 0 
               && (it->getIntensity() / int_last) > 0.1 // at least 10% of last sample point
               && it->getMZ() + 1 > it_max->getMZ())    // at most 1 Th away
        {
          int_last = it->getIntensity();
          explained_int += int_last;
          it->setIntensity(0); // remove peak from future consideration
          --it;
        }
        // if the current point is rising again, restore the intensity of the 
        // previous 'sink' point (because it could belong to a neighbour peak)
        // e.g. imagine intensities: 1-2-3-2-1-2-4-2-1. We do not want to destroy the middle '1'
        if (it->getIntensity() > int_last) (it+1)->setIntensity(int_last);

        //std::cerr << "  Peak candidate: " << it_max->getMZ() << " ...";
        bool break_left = false;
        if (it_max - it < 2+1)  // 'it' does not fulfill the conditions, i.e. does not count
        { // fewer than two sampling points on left shoulder
          //std::cerr << " break left " << it_max - it << " points\n";
          break_left = true;
          // do not end loop here.. we still need to clean up the right side
        }
        it_max->setIntensity(int_max); // restore center intensity
        explained_int -= int_max;
        it = it_max;
        int_last = int_max;
        while (it != data.end() 
               && it->getIntensity() <= int_last        // at most 100% of last sample point
               && it->getIntensity() > 0 
               && (it->getIntensity() / int_last) > 0.1 // at least 10% of last sample point
               && it->getMZ() - 1 < it_max->getMZ())    // at most 1 Th away
        {
          int_last = it->getIntensity();
          explained_int += int_last;
          it->setIntensity(0); // remove peak from future consideration
          ++it;
        }
        // if the current point is rising again, restore the intensity of the 
        // previous 'sink' point (because it could belong to a neighbour peak)
        // e.g. imagine intensities: 1-2-4-2-1-2-3-2-1. We do not want to destroy the middle '1'
        // (note: the sequence is not identical to the one of the left shoulder)
        if (it != data.end() && it->getIntensity() > int_last) (it-1)->setIntensity(int_last);

        if (break_left || it - it_max < 2+1)  // 'it' does not fulfill the conditions, i.e. does not count
        { // fewer than two sampling points on right shoulder
          //std::cerr << " break right " << it - it_max << " points\n";
          ++centroid_evidence;
          continue;
        }
        // peak has at least two sampling points on either side within 1 Th
        ++profile_evidence;
        //std::cerr << " PROFILE " << it - it_max << " points right\n";

      }
      
      float evidence_ratio = profile_evidence / float(profile_evidence + centroid_evidence);
      //std::cerr << "--> Evidence ratio: " << evidence_ratio; 

      if (evidence_ratio > 0.75) // 80% are profile
      {
        //std::cerr << "  PROFILE\n";
        return SpectrumSettings::PROFILE;
      }
      else
      { 
        //std::cerr << "  CENTROID\n";
        return SpectrumSettings::CENTROID;
      }
    }
  };

} // namespace OpenMS