MSExperiment.h

Go to the documentation of this file.

// Copyright (c) 2002-present, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: Marc Sturm, Tom Waschischeck $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/KERNEL/StandardDeclarations.h>
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/KERNEL/AreaIterator.h>
#include <OpenMS/KERNEL/MSChromatogram.h>
#include <OpenMS/KERNEL/MSSpectrum.h>
#include <OpenMS/METADATA/ExperimentalSettings.h>
 
#include <vector>
 
 
namespace OpenMS
{
  class Peak1D;
  class ChromatogramPeak;
 
  class OPENMS_DLLAPI MSExperiment final : public RangeManagerContainer<RangeRT, RangeMZ, RangeIntensity, RangeMobility>,
    public ExperimentalSettings
  {
 
public:
    typedef Peak1D PeakT;
    typedef ChromatogramPeak ChromatogramPeakT;
 
 
    typedef PeakT PeakType;
    typedef ChromatogramPeakT ChromatogramPeakType;
    typedef PeakType::CoordinateType CoordinateType;
    typedef PeakType::IntensityType IntensityType;
    typedef RangeManager<RangeRT, RangeMZ, RangeIntensity, RangeMobility> RangeManagerType;
    typedef RangeManagerContainer<RangeRT, RangeMZ, RangeIntensity, RangeMobility> RangeManagerContainerType;
    typedef MSSpectrum SpectrumType;
    typedef MSChromatogram ChromatogramType;
    typedef std::vector<SpectrumType> Base;
 
 
    typedef std::vector<SpectrumType>::iterator Iterator;
    typedef std::vector<SpectrumType>::const_iterator ConstIterator;
    typedef Internal::AreaIterator<PeakT, PeakT&, PeakT*, Iterator, SpectrumType::Iterator> AreaIterator;
    typedef Internal::AreaIterator<const PeakT, const PeakT&, const PeakT*, ConstIterator, SpectrumType::ConstIterator> ConstAreaIterator;
 
    // Attention: these refer to the spectra vector only!
    typedef Base::value_type value_type;
    typedef Base::iterator iterator;
    typedef Base::const_iterator const_iterator;
 
    MSExperiment();
 
    MSExperiment(const MSExperiment & source);
 
    MSExperiment(MSExperiment&&) = default;
 
    MSExperiment & operator=(const MSExperiment & source);
 
    MSExperiment& operator=(MSExperiment&&) & = default;
 
    MSExperiment & operator=(const ExperimentalSettings & source);
 
    ~MSExperiment() override;
 
    bool operator==(const MSExperiment & rhs) const;
 
    bool operator!=(const MSExperiment & rhs) const;
    
    inline Size size() const noexcept
    {
      return spectra_.size();
    }
 
    inline void resize(Size n)
    {
      spectra_.resize(n);
    }
 
    inline bool empty() const noexcept
    {
      return spectra_.empty();
    }
    
    inline void reserve(Size n)
    {
      spectra_.reserve(n);
    }
 
    inline SpectrumType& operator[](Size n)
    {
      return spectra_[n];
    }
 
    inline const SpectrumType& operator[](Size n) const
    {
      return spectra_[n];
    }
 
    inline Iterator begin() noexcept
    {
      return spectra_.begin();
    }
 
    inline ConstIterator begin() const noexcept
    {
      return spectra_.cbegin();
    }
 
    inline ConstIterator cbegin() const noexcept
    {
      return spectra_.cbegin();
    }
 
    inline Iterator end()
    {
      return spectra_.end();
    }
 
    inline ConstIterator end() const noexcept
    {
      return spectra_.cend();
    }
    
    inline ConstIterator cend() const noexcept
    {
      return spectra_.cend();
    }
 
    // Aliases / chromatograms
    void reserveSpaceSpectra(Size s);
    void reserveSpaceChromatograms(Size s);
 
 
 
    template <class Container>
    void get2DData(Container& cont) const
    {
      for (typename Base::const_iterator spec = spectra_.begin(); spec != spectra_.end(); ++spec)
      {
        if (spec->getMSLevel() != 1)
        {
          continue;
        }
        typename Container::value_type s; // explicit object here, since instantiation within push_back() fails on VS<12
        for (typename SpectrumType::const_iterator it = spec->begin(); it != spec->end(); ++it)
        {
          cont.push_back(s);
          cont.back().setRT(spec->getRT());
          cont.back().setMZ(it->getMZ());
          cont.back().setIntensity(it->getIntensity());
        }
      }
    }
 
    template <class Container>
    void set2DData(const Container& container)
    {
      set2DData<false, Container>(container);
    }
 
    template <class Container>
    void set2DData(const Container& container, const StringList& store_metadata_names)
    {
      // clean up the container first
      clear(true);
      SpectrumType* spectrum = nullptr;
      typename PeakType::CoordinateType current_rt = -std::numeric_limits<typename PeakType::CoordinateType>::max();
      for (typename Container::const_iterator iter = container.begin(); iter != container.end(); ++iter)
      {
        // check if the retention time has changed
        if (current_rt != iter->getRT() || spectrum == nullptr)
        {
          // append new spectrum
          if (current_rt > iter->getRT())
          {
            throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Input container is not sorted!");
          }
          current_rt =  iter->getRT();
          spectrum = createSpec_(current_rt, store_metadata_names);
        }
 
        // add either data point or mass traces (depending on template argument value)
        ContainerAdd_<typename Container::value_type, false>::addData_(spectrum, &(*iter), store_metadata_names);
      }
    }
 
    template <bool add_mass_traces, class Container>
    void set2DData(const Container& container)
    {
      // clean up the container first
      clear(true);
      SpectrumType* spectrum = nullptr;
      typename PeakType::CoordinateType current_rt = -std::numeric_limits<typename PeakType::CoordinateType>::max();
      for (typename Container::const_iterator iter = container.begin(); iter != container.end(); ++iter)
      {
        // check if the retention time has changed
        if (current_rt != iter->getRT() || spectrum == nullptr)
        {
          // append new spectrum
          if (current_rt > iter->getRT())
          {
            throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Input container is not sorted!");
          }
          current_rt =  iter->getRT();
          spectrum = createSpec_(current_rt);
        }
 
        // add either data point or mass traces (depending on template argument value)
        ContainerAdd_<typename Container::value_type, add_mass_traces>::addData_(spectrum, &(*iter));
      }
    }
 
 
 
 
    AreaIterator areaBegin(CoordinateType min_rt, CoordinateType max_rt, CoordinateType min_mz, CoordinateType max_mz, UInt  ms_level = 1);
 
    AreaIterator areaBegin(const RangeManagerType& range, UInt ms_level = 1);
 
    AreaIterator areaEnd();
 
    ConstAreaIterator areaBeginConst(CoordinateType min_rt, CoordinateType max_rt, CoordinateType min_mz, CoordinateType max_mz, UInt ms_level = 1) const;
 
    ConstAreaIterator areaBeginConst(const RangeManagerType& range, UInt ms_level = 1) const;
 
    ConstAreaIterator areaEndConst() const;
 
    // for fast pyOpenMS access to MS1 peak data in format: [rt, [mz, intensity]]
    void get2DPeakData(CoordinateType min_rt, CoordinateType max_rt, CoordinateType min_mz, CoordinateType max_mz, 
      std::vector<float>& rt, 
      std::vector<std::vector<float>>& mz, 
      std::vector<std::vector<float>>& intensity) const
    {
      float t = -1.0;
      for (auto it = areaBeginConst(min_rt, max_rt, min_mz, max_mz); it != areaEndConst(); ++it)
      {
        if (it.getRT() != t) 
        {
          t = (float)it.getRT();
          rt.push_back(t);
        }
        mz.back().push_back((float)it->getMZ());
        intensity.back().push_back(it->getIntensity());
      }
    }
 
    // for fast pyOpenMS access to MS1 peak data in format: [rt, [mz, intensity, ion mobility]]
    void get2DPeakData(CoordinateType min_rt, CoordinateType max_rt, CoordinateType min_mz, CoordinateType max_mz, 
      std::vector<float>& rt, 
      std::vector<std::vector<float>>& mz,
      std::vector<std::vector<float>>& intensity, 
      std::vector<std::vector<float>>& ion_mobility) const
    {
      float t = -1.0;
      for (auto it = areaBeginConst(min_rt, max_rt, min_mz, max_mz); it != areaEndConst(); ++it)
      {
        if (it.getRT() != t)
        {
          t = (float)it.getRT();
          rt.push_back(t);
        }
        
        const MSSpectrum& spectrum = it.getSpectrum();
        bool has_IM = spectrum.containsIMData();
        float peak_IM{-1.0f};
        if (has_IM)
        {
          const auto& im_data = spectrum.getIMData();
          const Size peak_index = it.getPeakIndex().peak;
          if (spectrum.getFloatDataArrays()[im_data.first].size() == spectrum.size())
          {
            peak_IM = spectrum.getFloatDataArrays()[im_data.first][peak_index];
          }          
        }
        ion_mobility.back().push_back(peak_IM);
        mz.back().push_back((float)it->getMZ());
        intensity.back().push_back(it->getIntensity());
      }
    }
 
    // for fast pyOpenMS access to MS1 peak data in format: [rt, mz, intensity]
    void get2DPeakData(
      CoordinateType min_rt,
      CoordinateType max_rt,
      CoordinateType min_mz,
      CoordinateType max_mz,
      std::vector<float>& rt,
      std::vector<float>& mz,
      std::vector<float>& intensity) const
    {
      for (auto it = areaBeginConst(min_rt, max_rt, min_mz, max_mz); it != areaEndConst(); ++it)
      {
        rt.push_back((float)it.getRT());
        mz.push_back((float)it->getMZ());
        intensity.push_back(it->getIntensity());
      }
    }
 
    // for fast pyOpenMS access to MS1 peak data in format: [rt, mz, intensity, ion mobility]
    void get2DPeakDataIon(
      CoordinateType min_rt,
      CoordinateType max_rt,
      CoordinateType min_mz,
      CoordinateType max_mz,
      std::vector<float>& rt,
      std::vector<float>& mz,
      std::vector<float>& intensity,
      std::vector<float>& ion_mobility) const
    {
      for (auto it = areaBeginConst(min_rt, max_rt, min_mz, max_mz); it != areaEndConst(); ++it)
      {
        rt.push_back((float)it.getRT());
        mz.push_back((float)it->getMZ());
        intensity.push_back(it->getIntensity());
 
        const MSSpectrum& spectrum = it.getSpectrum();
        bool has_IM = spectrum.containsIMData();
        float peak_IM = -1.0;
        if (has_IM)
        {
          const auto& im_data = spectrum.getIMData();
          const Size& peak_index = it.getPeakIndex().peak;
          if (spectrum.getFloatDataArrays()[im_data.first].size() == spectrum.size())
          {
            peak_IM = spectrum.getFloatDataArrays()[im_data.first][peak_index];
          }          
        }        
        ion_mobility.push_back(peak_IM);
      }
    }
 
    ConstIterator RTBegin(CoordinateType rt) const;
 
    ConstIterator RTEnd(CoordinateType rt) const;
 
    Iterator RTBegin(CoordinateType rt);
 
    Iterator RTEnd(CoordinateType rt);
 
 
    ConstIterator IMBegin(CoordinateType im) const;
 
    ConstIterator IMEnd(CoordinateType im) const;
 
    // Docu in base class
    void updateRanges() override;
 
    void updateRanges(Int ms_level);
 
    UInt64 getSize() const;
 
    const std::vector<UInt>& getMSLevels() const;
 
 
    UInt64 getSqlRunID() const;
 
    void setSqlRunID(UInt64 id);
 
 
    void sortSpectra(bool sort_mz = true);
 
    void sortChromatograms(bool sort_rt = true);
 
    bool isSorted(bool check_mz = true) const;
 
 
    void reset();
 
    bool clearMetaDataArrays();
 
    const ExperimentalSettings& getExperimentalSettings() const;
 
    ExperimentalSettings& getExperimentalSettings();
 
    void getPrimaryMSRunPath(StringList& toFill) const;
 
    ConstIterator getPrecursorSpectrum(ConstIterator iterator) const;
 
    int getPrecursorSpectrum(int zero_based_index) const;
 
    void swap(MSExperiment& from);
 
    void setSpectra(const std::vector<MSSpectrum>& spectra);
    void setSpectra(std::vector<MSSpectrum>&& spectra);
 
    void addSpectrum(const MSSpectrum& spectrum);
    void addSpectrum(MSSpectrum&& spectrum);
 
    const std::vector<MSSpectrum>& getSpectra() const;
 
    std::vector<MSSpectrum>& getSpectra();
 
    ConstIterator getClosestSpectrumInRT(const double RT) const;
    Iterator getClosestSpectrumInRT(const double RT);
 
    ConstIterator getClosestSpectrumInRT(const double RT, UInt ms_level) const;
    Iterator getClosestSpectrumInRT(const double RT, UInt ms_level);
 
    void setChromatograms(const std::vector<MSChromatogram>& chromatograms);
    void setChromatograms(std::vector<MSChromatogram>&& chromatograms);
 
    void addChromatogram(const MSChromatogram& chromatogram);
    void addChromatogram(MSChromatogram&& chrom);
 
    const std::vector<MSChromatogram>& getChromatograms() const;
 
    std::vector<MSChromatogram>& getChromatograms();
 
 
    MSChromatogram& getChromatogram(Size id);
 
    MSSpectrum& getSpectrum(Size id);
 
    Size getNrSpectra() const;
 
    Size getNrChromatograms() const;
 
    const MSChromatogram calculateTIC(float rt_bin_size = 0, UInt ms_level = 1) const;
 
    void clear(bool clear_meta_data);
 
    bool containsScanOfLevel(size_t ms_level) const;
 
    bool hasZeroIntensities(size_t ms_level) const;
 
    bool hasPeptideIdentifications() const;
 
    bool isIMFrame() const;
 
  protected:
    std::vector<UInt> ms_levels_;
    UInt64 total_size_;
    std::vector<MSChromatogram > chromatograms_;
    std::vector<SpectrumType> spectra_;
 
private:
 
    template<typename ContainerValueType, bool addMassTraces>
    struct ContainerAdd_
    {
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item);
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item, const StringList& store_metadata_names);
    };
 
    template<typename ContainerValueType>
    struct ContainerAdd_<ContainerValueType, false>
    {
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item)
      {
        // create temporary peak and insert it into spectrum
        spectrum->insert(spectrum->end(), PeakType());
        spectrum->back().setIntensity(item->getIntensity());
        spectrum->back().setPosition(item->getMZ());
      }
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item, const StringList& store_metadata_names)
      {
        addData_(spectrum, item);
        for (StringList::const_iterator itm = store_metadata_names.begin(); itm != store_metadata_names.end(); ++itm)
        {
          float val = std::numeric_limits<float>::quiet_NaN();
          if (item->metaValueExists(*itm)) val = item->getMetaValue(*itm);
          spectrum->getFloatDataArrays()[itm - store_metadata_names.begin()].push_back(val);
        }
      }
    };
 
    template<typename ContainerValueType>
    struct ContainerAdd_<ContainerValueType, true>
    {
      static void addData_(SpectrumType* spectrum, const ContainerValueType* item)
      {
        if (item->metaValueExists("num_of_masstraces"))
        {
          Size mts = item->getMetaValue("num_of_masstraces");
          int charge = (item->getCharge()==0 ? 1 : item->getCharge()); // set to 1 if charge is 0, otherwise div/0 below
          for (Size i = 0; i < mts; ++i)
          {
            String meta_name = String("masstrace_intensity_") + i;
            if (!item->metaValueExists(meta_name))
            {
              throw Exception::Precondition(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, String("Meta value '") + meta_name + "' expected but not found in container.");
            }
            ContainerValueType p;
            p.setIntensity(item->getMetaValue(meta_name));
            p.setPosition(item->getMZ() + Constants::C13C12_MASSDIFF_U / charge * i);
            ContainerAdd_<ContainerValueType, false>::addData_(spectrum, &p);
          }
        }
        else ContainerAdd_<ContainerValueType, false>::addData_(spectrum, item);
      }
    };
 
    /*
      @brief Append a spectrum to current MSExperiment
 
      @param rt RT of new spectrum
      @return Pointer to newly created spectrum
    */
    SpectrumType* createSpec_(PeakType::CoordinateType rt);
 
    /*
      @brief Append a spectrum including floatdata arrays to current MSExperiment
 
      @param rt RT of new spectrum
      @param metadata_names Names of floatdata arrays attached to this spectrum
      @return Pointer to newly created spectrum
    */
    SpectrumType* createSpec_(PeakType::CoordinateType rt, const StringList& metadata_names);
 
  };
 
  OPENMS_DLLAPI std::ostream& operator<<(std::ostream& os, const MSExperiment& exp);
 
} // namespace OpenMS
 
#include <OpenMS/KERNEL/StandardTypes.h>