MSExperiment.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-2012.
//
// 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: Stephan Aiche$
// $Authors: Marc Sturm $
// --------------------------------------------------------------------------

#ifndef OPENMS_KERNEL_MSEXPERIMENT_H
#define OPENMS_KERNEL_MSEXPERIMENT_H

#include <OpenMS/KERNEL/MSSpectrum.h>
#include <OpenMS/KERNEL/MSChromatogram.h>
#include <OpenMS/METADATA/ExperimentalSettings.h>
#include <OpenMS/DATASTRUCTURES/DRange.h>
#include <OpenMS/FORMAT/DB/PersistentObject.h>
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/KERNEL/AreaIterator.h>

#include <vector>
#include <algorithm>
#include <limits>

namespace OpenMS
{
  class Peak1D;

  template <typename PeakT = Peak1D, typename ChromatogramPeakT = ChromatogramPeak>
  class MSExperiment :
    public std::vector<MSSpectrum<PeakT> >,
    public RangeManager<2>,
    public ExperimentalSettings,
    public PersistentObject
  {
public:


    typedef PeakT PeakType;
    typedef ChromatogramPeakT ChromatogramPeakType;
    typedef DRange<2> AreaType;
    typedef typename PeakType::CoordinateType CoordinateType;
    typedef typename PeakType::IntensityType IntensityType;
    typedef RangeManager<2> RangeManagerType;
    typedef MSSpectrum<PeakType> SpectrumType;
    typedef MSChromatogram<ChromatogramPeakType> ChromatogramType;
    typedef std::vector<SpectrumType> Base;



    typedef typename std::vector<SpectrumType>::iterator Iterator;
    typedef typename std::vector<SpectrumType>::const_iterator ConstIterator;
    typedef Internal::AreaIterator<PeakT, PeakT &, PeakT *, Iterator, typename SpectrumType::Iterator> AreaIterator;
    typedef Internal::AreaIterator<const PeakT, const PeakT &, const PeakT *, ConstIterator, typename SpectrumType::ConstIterator> ConstAreaIterator;

    MSExperiment() :
      Base(),
      RangeManagerType(),
      ExperimentalSettings(),
      PersistentObject(),
      ms_levels_(),
      total_size_(0)
    {}

    MSExperiment(const MSExperiment & source) :
      std::vector<MSSpectrum<PeakT> >(source),
      RangeManagerType(source),
      ExperimentalSettings(source),
      PersistentObject(source),
      ms_levels_(source.ms_levels_),
      total_size_(source.total_size_),
      chromatograms_(source.chromatograms_)
    {}

    MSExperiment & operator=(const MSExperiment & source)
    {
      if (&source == this) return *this;

      Base::operator=(source);
      RangeManagerType::operator=(source);
      ExperimentalSettings::operator=(source);
      PersistentObject::operator=(source);

      ms_levels_     = source.ms_levels_;
      total_size_    = source.total_size_;
      chromatograms_ = source.chromatograms_;

      //no need to copy the alloc?!
      //alloc_

      return *this;
    }

    MSExperiment & operator=(const ExperimentalSettings & source)
    {
      ExperimentalSettings::operator=(source);
      return *this;
    }

    bool operator==(const MSExperiment & rhs) const
    {
      return ExperimentalSettings::operator==(rhs) && std::operator==(rhs, *this) && chromatograms_ == rhs.chromatograms_;
    }

    bool operator!=(const MSExperiment & rhs) const
    {
      return !(operator==(rhs));
    }



    template <class Container>
    void get2DData(Container & cont) const
    {
      for (typename Base::const_iterator spec = Base::begin(); spec != Base::end(); ++spec)
      {
        if (spec->getMSLevel() != 1)
        {
          continue;
        }
        for (typename SpectrumType::const_iterator it = spec->begin(); it != spec->end(); ++it)
        {
          cont.push_back(typename Container::value_type());
          cont.back().setRT(spec->getRT());
          cont.back().setMZ(it->getMZ());
          cont.back().setIntensity(it->getIntensity());
        }
      }
    }

    template <class Container>
    void set2DData(const Container & cont)
    {
      SpectrumType * spectrum = 0;
      // If the container is empty, nothing will happen
      if (cont.empty()) return;

      typename PeakType::CoordinateType current_rt = -(std::numeric_limits<typename PeakType::CoordinateType>::max)();

      for (typename Container::const_iterator iter = cont.begin(); iter != cont.end(); ++iter)
      {
        // check if the retention time time has changed
        if (current_rt != iter->getRT() || spectrum == 0)
        {
          if (current_rt > iter->getRT())
          {
            throw Exception::Precondition(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Input container is not sorted!");
          }
          current_rt =  iter->getRT();
          Base::insert(Base::end(), SpectrumType());
          spectrum = &(Base::back());
          spectrum->setRT(current_rt);
          spectrum->setMSLevel(1);
        }

        // create temporary peak and insert it into spectrum
        spectrum->insert(spectrum->end(), PeakType());
        spectrum->back().setIntensity(iter->getIntensity());
        spectrum->back().setPosition(iter->getMZ());
      }
    }




    AreaIterator areaBegin(CoordinateType min_rt, CoordinateType max_rt, CoordinateType min_mz, CoordinateType max_mz)
    {
      OPENMS_PRECONDITION(min_rt <= max_rt, "Swapped RT range boundaries!")
      OPENMS_PRECONDITION(min_mz <= max_mz, "Swapped MZ range boundaries!")
      //std::cout << "areaBegin: " << min_rt << " " << max_rt << " " << min_mz << " " << max_mz << std::endl;
      return AreaIterator(this->begin(), RTBegin(min_rt), RTEnd(max_rt), min_mz, max_mz);
    }

    AreaIterator areaEnd()
    {
      return AreaIterator();
    }

    ConstAreaIterator areaBeginConst(CoordinateType min_rt, CoordinateType max_rt, CoordinateType min_mz, CoordinateType max_mz) const
    {
      OPENMS_PRECONDITION(min_rt <= max_rt, "Swapped RT range boundaries!")
      OPENMS_PRECONDITION(min_mz <= max_mz, "Swapped MZ range boundaries!")
      //std::cout << "areaBeginConst: " << min_rt << " " << max_rt << " " << min_mz << " " << max_mz << std::endl;
      return ConstAreaIterator(this->begin(), RTBegin(min_rt), RTEnd(max_rt), min_mz, max_mz);
    }

    ConstAreaIterator areaEndConst() const
    {
      return ConstAreaIterator();
    }

    ConstIterator RTBegin(CoordinateType rt) const
    {
      SpectrumType s;
      s.setRT(rt);
      return lower_bound(Base::begin(), Base::end(), s, typename SpectrumType::RTLess());
    }

    ConstIterator RTEnd(CoordinateType rt) const
    {
      SpectrumType s;
      s.setRT(rt);
      return upper_bound(Base::begin(), Base::end(), s, typename SpectrumType::RTLess());
    }

    Iterator RTBegin(CoordinateType rt)
    {
      SpectrumType s;
      s.setRT(rt);
      return lower_bound(Base::begin(), Base::end(), s, typename SpectrumType::RTLess());
    }

    Iterator RTEnd(CoordinateType rt)
    {
      SpectrumType s;
      s.setRT(rt);
      return upper_bound(Base::begin(), Base::end(), s, typename SpectrumType::RTLess());
    }



    // Docu in base class
    virtual void updateRanges()
    {
      updateRanges(-1);
    }

    void updateRanges(Int ms_level)
    {
      //clear MS levels
      ms_levels_.clear();

      //reset mz/rt/int range
      this->clearRanges();
      //reset point count
      total_size_ = 0;

      //empty
      if (this->size() == 0 && chromatograms_.empty())
      {
        return;
      }

      //update
      for (typename Base::iterator it = this->begin(); it != this->end(); ++it)
      {
        if (ms_level < Int(0) || Int(it->getMSLevel()) == ms_level)
        {
          //ms levels
          if (std::find(ms_levels_.begin(), ms_levels_.end(), it->getMSLevel()) == ms_levels_.end())
          {
            ms_levels_.push_back(it->getMSLevel());
          }

          // calculate size
          total_size_ += it->size();

          //rt
          if (it->getRT() < RangeManagerType::pos_range_.minX()) RangeManagerType::pos_range_.setMinX(it->getRT());
          if (it->getRT() > RangeManagerType::pos_range_.maxX()) RangeManagerType::pos_range_.setMaxX(it->getRT());

          //do not update mz and int when the spectrum is empty
          if (it->size() == 0) continue;

          it->updateRanges();

          //mz
          if (it->getMin()[0] < RangeManagerType::pos_range_.minY()) RangeManagerType::pos_range_.setMinY(it->getMin()[0]);
          if (it->getMax()[0] > RangeManagerType::pos_range_.maxY()) RangeManagerType::pos_range_.setMaxY(it->getMax()[0]);

          //int
          if (it->getMinInt() < RangeManagerType::int_range_.minX()) RangeManagerType::int_range_.setMinX(it->getMinInt());
          if (it->getMaxInt() > RangeManagerType::int_range_.maxX()) RangeManagerType::int_range_.setMaxX(it->getMaxInt());

        }
        // for MS level = 1 we extend the range for all the MS2 precursors
        if (ms_level == 1 && it->getMSLevel() == 2)
        {
          if (!it->getPrecursors().empty())
          {
            DoubleReal pc_rt = it->getRT();
            if (pc_rt < RangeManagerType::pos_range_.minX()) RangeManagerType::pos_range_.setMinX(pc_rt);
            if (pc_rt > RangeManagerType::pos_range_.maxX()) RangeManagerType::pos_range_.setMaxX(pc_rt);
            DoubleReal pc_mz = it->getPrecursors()[0].getMZ();
            if (pc_mz < RangeManagerType::pos_range_.minY()) RangeManagerType::pos_range_.setMinY(pc_mz);
            if (pc_mz > RangeManagerType::pos_range_.maxY()) RangeManagerType::pos_range_.setMaxY(pc_mz);
          }

        }

      }
      std::sort(ms_levels_.begin(), ms_levels_.end());




      if (this->chromatograms_.empty())
      {
        return;
      }

      //TODO CHROM update intensity, m/z and RT according to chromatograms as well! (done????)

      for (typename std::vector<ChromatogramType>::iterator it = chromatograms_.begin(); it != chromatograms_.end(); ++it)
      {

        // ignore TICs and ECs (as these are usually positioned at 0 and therefor lead to a large white margin in plots if included)
        if (it->getChromatogramType() == ChromatogramSettings::TOTAL_ION_CURRENT_CHROMATOGRAM ||
            it->getChromatogramType() == ChromatogramSettings::EMISSION_CHROMATOGRAM)
        {
          continue;
        }

        // update MZ
        if (it->getMZ() < RangeManagerType::pos_range_.minY()) RangeManagerType::pos_range_.setMinY(it->getMZ());
        if (it->getMZ() > RangeManagerType::pos_range_.maxY()) RangeManagerType::pos_range_.setMaxY(it->getMZ());

        // do not update RT and in if the specturm is empty
        if (it->size() == 0) continue;

        total_size_ += it->size();

        it->updateRanges();

        // RT
        if (it->getMin()[0] < RangeManagerType::pos_range_.minX()) RangeManagerType::pos_range_.setMinX(it->getMin()[0]);
        if (it->getMax()[0] > RangeManagerType::pos_range_.maxX()) RangeManagerType::pos_range_.setMaxX(it->getMax()[0]);

        // int
        if (it->getMinInt() < RangeManagerType::int_range_.minX()) RangeManagerType::int_range_.setMinX(it->getMinInt());
        if (it->getMaxInt() > RangeManagerType::int_range_.maxX()) RangeManagerType::int_range_.setMaxX(it->getMaxInt());
      }
    }

    CoordinateType getMinMZ() const
    {
      return RangeManagerType::pos_range_.minPosition()[1];
    }

    CoordinateType getMaxMZ() const
    {
      return RangeManagerType::pos_range_.maxPosition()[1];
    }

    CoordinateType getMinRT() const
    {
      return RangeManagerType::pos_range_.minPosition()[0];
    }

    CoordinateType getMaxRT() const
    {
      return RangeManagerType::pos_range_.maxPosition()[0];
    }

    const AreaType & getDataRange() const
    {
      return RangeManagerType::pos_range_;
    }

    UInt64 getSize() const
    {
      return total_size_;
    }

    const std::vector<UInt> & getMSLevels() const
    {
      return ms_levels_;
    }



    void sortSpectra(bool sort_mz = true)
    {
      std::sort(this->begin(), this->end(), typename SpectrumType::RTLess());

      if (sort_mz)
      {
        // sort each spectrum by m/z
        for (Iterator iter = this->begin(); iter != this->end(); ++iter)
        {
          iter->sortByPosition();
        }
      }
    }

    void sortChromatograms(bool sort_rt = true)
    {
      // sort the chromatograms according to their product m/z
      std::sort(chromatograms_.begin(), chromatograms_.end(), typename ChromatogramType::MZLess());

      if (sort_rt)
      {
        for (typename std::vector<ChromatogramType>::iterator it = chromatograms_.begin(); it != chromatograms_.end(); ++it)
        {
          it->sortByPosition();
        }
      }
    }

    bool isSorted(bool check_mz = true) const
    {
      //check RT positions
      for (Size i = 1; i < this->size(); ++i)
      {
        if (this->operator[](i - 1).getRT() > this->operator[](i).getRT()) return false;
      }
      //check spectra
      if (check_mz)
      {
        for (Size i = 0; i < this->size(); ++i)
        {
          if (!this->operator[](i).isSorted()) return false;
        }
      }
      // TODO CHROM
      return true;
    }


    void reset()
    {
      Base::clear();           //remove data
      RangeManagerType::clearRanges();           //reset range manager
      ExperimentalSettings::operator=(ExperimentalSettings());           //reset meta info
    }

    bool clearMetaDataArrays()
    {
      bool meta_present = false;
      for (Size i = 0; i < this->size(); ++i)
      {
        if (this->operator[](i).getFloatDataArrays().size() != 0 || this->operator[](i).getIntegerDataArrays().size() != 0 || this->operator[](i).getStringDataArrays().size() != 0)
        {
          meta_present = true;
        }
        this->operator[](i).getStringDataArrays().clear();
        this->operator[](i).getIntegerDataArrays().clear();
        this->operator[](i).getFloatDataArrays().clear();
      }
      return meta_present;
    }

    const ExperimentalSettings & getExperimentalSettings() const
    {
      return *this;
    }

    ExperimentalSettings & getExperimentalSettings()
    {
      return *this;
    }

    ConstIterator getPrecursorSpectrum(ConstIterator iterator) const
    {
      if (iterator == this->end() || iterator == this->begin())
      {
        return this->end();
      }
      UInt ms_level = iterator->getMSLevel();
      do
      {
        --iterator;
        if (iterator->getMSLevel() < ms_level)
        {
          return iterator;
        }
      }
      while (iterator != this->begin());

      return this->end();
    }

    void swap(MSExperiment & from)
    {
      MSExperiment tmp;

      //swap range information
      tmp.RangeManagerType::operator=(* this);
      this->RangeManagerType::operator=(from);
      from.RangeManagerType::operator=(tmp);

      //swap experimental settings
      tmp.ExperimentalSettings::operator=(* this);
      this->ExperimentalSettings::operator=(from);
      from.ExperimentalSettings::operator=(tmp);

      //swap persistent object
      tmp.PersistentObject::operator=(* this);
      this->PersistentObject::operator=(from);
      from.PersistentObject::operator=(tmp);

      // swap chromatograms
      std::swap(chromatograms_, from.chromatograms_);

      //swap peaks
      Base::swap(from);

      //swap remaining members
      ms_levels_.swap(from.ms_levels_);
      std::swap(total_size_, from.total_size_);
    }

    void setChromatograms(const std::vector<MSChromatogram<ChromatogramPeakType> > & chromatograms)
    {
      chromatograms_ = chromatograms;
    }

    void addChromatogram(const MSChromatogram<ChromatogramPeakType> & chromatogram)
    {
      chromatograms_.push_back(chromatogram);
    }

    const std::vector<MSChromatogram<ChromatogramPeakType> > & getChromatograms() const
    {
      return chromatograms_;
    }

    const MSChromatogram<ChromatogramPeakType> getTIC() const
    {
      // The TIC is (re)calculated from the MS1 spectra. Even if MSExperiment does not contain a TIC chromatogram explicitly, it can be reported.
      MSChromatogram<ChromatogramPeakType> TIC;
      for (typename Base::const_iterator spec_it = this->begin(); spec_it != this->end(); ++spec_it)
      {
        if (spec_it->getMSLevel() == 1)
        {
          DoubleReal totalIntensity = 0;
          // sum intensities of a spectrum
          for (typename SpectrumType::const_iterator peak_it = spec_it->begin(); peak_it != spec_it->end(); ++peak_it)
          {
            totalIntensity += peak_it->getIntensity();
          }
          // fill chromatogram
          ChromatogramPeakType peak;
          peak.setRT(spec_it->getRT());
          peak.setIntensity(totalIntensity);
          TIC.push_back(peak);
        }
      }
      return TIC;
    }

    void clear(bool clear_meta_data)
    {
      Base::clear();

      if (clear_meta_data)
      {
        clearRanges();
        clearId();
        this->ExperimentalSettings::operator=(ExperimentalSettings());             // no "clear" method
        chromatograms_.clear();
        ms_levels_.clear();
        total_size_ = 0;
      }
    }

protected:

    // Docu in base class
    virtual void clearChildIds_()
    {
      for (Size i = 0; i < this->size(); ++i)
      {
        this->operator[](i).clearId(true);
      }
    }

    std::vector<UInt> ms_levels_;
    UInt64 total_size_;

    std::vector<MSChromatogram<ChromatogramPeakType> > chromatograms_;
  };

  template <typename PeakT, typename ChromatogramPeakT>
  std::ostream & operator<<(std::ostream & os, const MSExperiment<PeakT, ChromatogramPeakT> & exp)
  {
    os << "-- MSEXPERIMENT BEGIN --" << std::endl;

    //experimental settings
    os << static_cast<const ExperimentalSettings &>(exp);

    //spectra
    for (typename MSExperiment<PeakT>::const_iterator it = exp.begin(); it != exp.end(); ++it)
    {
      os << *it;
    }

    //chromatograms
    for (typename std::vector<MSChromatogram<ChromatogramPeakT> >::const_iterator it = exp.getChromatograms().begin(); it != exp.getChromatograms().end(); ++it)
    {
      os << *it;
    }

    os << "-- MSEXPERIMENT END --" << std::endl;

    return os;
  }

} // namespace OpenMS

#endif // OPENMS_KERNEL_MSEXPERIMENT_H