ChromatogramExtractor.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-2017.
//
// 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: Hannes Roest $
// $Authors: Hannes Roest $
// --------------------------------------------------------------------------

#pragma once

#include <OpenMS/ANALYSIS/OPENSWATH/ChromatogramExtractorAlgorithm.h>

#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/ANALYSIS/TARGETED/TargetedExperiment.h>
#include <OpenMS/ANALYSIS/MAPMATCHING/TransformationDescription.h>

#include <OpenMS/ANALYSIS/OPENSWATH/DATAACCESS/DataAccessHelper.h>

namespace OpenMS
{

  class OPENMS_DLLAPI ChromatogramExtractor :
    public ProgressLogger
  {

public:

    typedef ChromatogramExtractorAlgorithm::ExtractionCoordinates ExtractionCoordinates;

    template <typename ExperimentT>
    void extractChromatograms(const ExperimentT& input, ExperimentT& output, 
        OpenMS::TargetedExperiment& transition_exp, double mz_extraction_window, bool ppm,
        TransformationDescription trafo, double rt_extraction_window, const String& filter)
    {
      // invert the trafo because we want to transform nRT values to "real" RT values
      trafo.invert();

      Size input_size = input.size();
      if (input_size < 1)
      {
        return;
      }

      int used_filter = getFilterNr_(filter);
      populatePeptideRTMap_(transition_exp, rt_extraction_window);

      // sort the transition experiment by product mass
      // this is essential because the algorithm assumes sorted transitions!
      transition_exp.sortTransitionsByProductMZ();

      // prepare all the spectra (but leave them empty)
      SpectrumSettings settings = input[0];
      std::vector<typename ExperimentT::ChromatogramType> chromatograms;
      prepareSpectra_(settings, chromatograms, transition_exp);

      //go through all spectra
      startProgress(0, input_size, "Extracting chromatograms");
      for (Size scan_idx = 0; scan_idx < input_size; ++scan_idx)
      {
        setProgress(scan_idx);

        if (input[scan_idx].size() == 0)
          continue;

        Size peak_idx = 0;

        double mz;
        double integrated_intensity = 0;

        // go through all transitions / chromatograms which are sorted by
        // ProductMZ. We can use this to step through the spectrum and at the
        // same time step through the transitions. We increase the peak counter
        // until we hit the next transition and then extract the signal.
        for (Size k = 0; k < chromatograms.size(); ++k)
        {

          double current_rt = input[scan_idx].getRT();
          if (outsideExtractionWindow_(transition_exp.getTransitions()[k], current_rt, trafo, rt_extraction_window))
          {
            continue;
          }

          typename ExperimentT::ChromatogramType::PeakType p;
          mz = transition_exp.getTransitions()[k].getProductMZ();

          if (used_filter == 1)
          {
            extract_value_tophat(input[scan_idx], mz, peak_idx, integrated_intensity, mz_extraction_window, ppm);
          }
          else if (used_filter == 2)
          {
            extract_value_bartlett(input[scan_idx], mz, peak_idx, integrated_intensity, mz_extraction_window, ppm);
          }


          p.setRT(current_rt);
          p.setIntensity(integrated_intensity);
          chromatograms[k].push_back(p);
        }
      }
      endProgress();

      // add all the chromatograms to the output
      output.setChromatograms(chromatograms);
    }

    void extractChromatograms(const OpenSwath::SpectrumAccessPtr input, 
        std::vector< OpenSwath::ChromatogramPtr >& output, 
        const std::vector<ExtractionCoordinates>& extraction_coordinates,
        double mz_extraction_window, bool ppm, const String& filter)
    {
      ChromatogramExtractorAlgorithm().extractChromatograms(input, output, 
          extraction_coordinates, mz_extraction_window, ppm, -1, filter);
    }

    void extractChromatograms(const OpenSwath::SpectrumAccessPtr input, 
        std::vector< OpenSwath::ChromatogramPtr >& output, 
        const std::vector<ExtractionCoordinates>& extraction_coordinates,
        double mz_extraction_window, bool ppm, double im_extraction_window, const String& filter)
    {
      ChromatogramExtractorAlgorithm().extractChromatograms(input, output, 
          extraction_coordinates, mz_extraction_window, ppm, im_extraction_window, filter);
    }

public:

    void prepare_coordinates(std::vector< OpenSwath::ChromatogramPtr > & output_chromatograms,
      std::vector< ExtractionCoordinates > & coordinates,
      const OpenMS::TargetedExperiment & transition_exp,
      const double rt_extraction_window,
      const bool ms1) const;

    template <typename TransitionExpT>
    static void return_chromatogram(std::vector< OpenSwath::ChromatogramPtr > & chromatograms,
      std::vector< ChromatogramExtractor::ExtractionCoordinates > & coordinates,
      TransitionExpT& transition_exp_used,
      SpectrumSettings settings,
      std::vector<OpenMS::MSChromatogram > & output_chromatograms,
      bool ms1)
    {
      typedef std::map<String, const typename TransitionExpT::Transition* > TransitionMapType;
      TransitionMapType trans_map;
      for (Size i = 0; i < transition_exp_used.getTransitions().size(); i++)
      {
        trans_map[transition_exp_used.getTransitions()[i].getNativeID()] = &transition_exp_used.getTransitions()[i];
      }

      for (Size i = 0; i < chromatograms.size(); i++)
      {
        const OpenSwath::ChromatogramPtr & chromptr = chromatograms[i];
        const ChromatogramExtractor::ExtractionCoordinates & coord = coordinates[i];

        // copy data
        OpenMS::MSChromatogram chrom;
        OpenSwathDataAccessHelper::convertToOpenMSChromatogram(chromptr, chrom);
        chrom.setNativeID(coord.id);

        // Create precursor and set
        // 1) the target m/z
        // 2) the isolation window (upper/lower)
        // 3) the peptide sequence
        Precursor prec;
        if (ms1) 
        {
          prec.setMZ(coord.mz);
          chrom.setChromatogramType(ChromatogramSettings::BASEPEAK_CHROMATOGRAM);

          // extract compound / peptide id from transition and store in
          // more-or-less default field
          String r = extract_id_(transition_exp_used, coord.id);
          prec.setMetaValue("peptide_sequence", r);
        }
        else 
        {
          typename TransitionExpT::Transition transition = (*trans_map[coord.id]);

          prec.setMZ(transition.getPrecursorMZ());
          if (settings.getPrecursors().size() > 0)
          {
            prec.setIsolationWindowLowerOffset(settings.getPrecursors()[0].getIsolationWindowLowerOffset());
            prec.setIsolationWindowUpperOffset(settings.getPrecursors()[0].getIsolationWindowUpperOffset());
          }

          // Create product and set its m/z
          Product prod;
          prod.setMZ(transition.getProductMZ());
          chrom.setProduct(prod);
          chrom.setChromatogramType(ChromatogramSettings::SELECTED_REACTION_MONITORING_CHROMATOGRAM);

          // extract compound / peptide id from transition and store in
          // more-or-less default field
          if (!transition.getPeptideRef().empty())
          {
            String r = extract_id_(transition_exp_used, transition.getPeptideRef());
            prec.setMetaValue("peptide_sequence", r);
          }
          else
          {
            String r = extract_id_(transition_exp_used, transition.getCompoundRef());
            prec.setMetaValue("peptide_sequence", r);
          }
        }
        chrom.setPrecursor(prec);

        // Set the rest of the meta-data
        chrom.setInstrumentSettings(settings.getInstrumentSettings());
        chrom.setAcquisitionInfo(settings.getAcquisitionInfo());
        chrom.setSourceFile(settings.getSourceFile());

        for (Size j = 0; j < settings.getDataProcessing().size(); ++j)
        {
          settings.getDataProcessing()[j]->setMetaValue("performed_on_spectra", "true");
          chrom.getDataProcessing().push_back(settings.getDataProcessing()[j]);
        }
        output_chromatograms.push_back(chrom);
      }
    }

    template <typename SpectrumT>
    void extract_value_tophat(const SpectrumT& input, const double& mz, Size& peak_idx,
        double& integrated_intensity, const double& extract_window, const bool ppm)
    {
      integrated_intensity = 0;
      if (input.size() == 0)
      {
        return;
      }

      // calculate extraction window
      double left, right;
      if (ppm)
      {
        left  = mz - mz * extract_window / 2.0 * 1.0e-6;
        right = mz + mz * extract_window / 2.0 * 1.0e-6;
      }
      else
      {
        left  = mz - extract_window / 2.0;
        right = mz + extract_window / 2.0;
      }

      Size walker;

      // advance the peak_idx until we hit the m/z value of the next transition
      while (peak_idx < input.size() && input[peak_idx].getMZ() < mz)
      {
        peak_idx++;
      }

      // walk right and left and add to our intensity
      walker = peak_idx;
      // if we moved past the end of the spectrum, we need to try the last peak
      // of the spectrum (it could still be within the window)
      if (peak_idx >= input.size())
      {
        walker = input.size() - 1;
      }

      // add the current peak if it is between right and left
      if (input[walker].getMZ() > left && input[walker].getMZ() < right)
      {
        integrated_intensity += input[walker].getIntensity();
      }

      // (i) Walk to the left one step and then keep walking left until we go
      // outside the window. Note for the first step to the left we have to
      // check for the walker becoming zero.
      walker = peak_idx;
      if (walker > 0)
      {
        walker--;
        // special case: walker is now zero
        if (walker == 0 && input[walker].getMZ() > left && input[walker].getMZ() < right)
        {
          integrated_intensity += input[walker].getIntensity();
        }
      }
      while (walker > 0 && input[walker].getMZ() > left && input[walker].getMZ() < right)
      {
        integrated_intensity += input[walker].getIntensity(); walker--;
      }

      // (ii) Walk to the right one step and then keep walking right until we
      // are outside the window
      walker = peak_idx;
      if (walker < input.size() )
      {
        walker++;
      }
      while (walker < input.size() && input[walker].getMZ() > left &&  input[walker].getMZ() < right)
      {
        integrated_intensity += input[walker].getIntensity(); walker++;
      }
    }

    template <typename SpectrumT>
    void extract_value_bartlett(const SpectrumT& input, const double& mz, Size& peak_idx,
        double& integrated_intensity, const double& extract_window, const bool ppm)
    {
      integrated_intensity = 0;
      if (input.size() == 0)
      {
        return;
      }

      // calculate extraction window
      double left, right, half_window_size, weight;
      if (ppm)
      {
        half_window_size = mz * extract_window / 2.0 * 1.0e-6;
        left  = mz - mz * extract_window / 2.0 * 1.0e-6;
        right = mz + mz * extract_window / 2.0 * 1.0e-6;
      }
      else
      {
        half_window_size = extract_window / 2.0;
        left  = mz - extract_window / 2.0;
        right = mz + extract_window / 2.0;
      }

      Size walker;

      // advance the peak_idx until we hit the m/z value of the next transition
      while (peak_idx < input.size() && input[peak_idx].getMZ() < mz)
      {
        peak_idx++;
      }

      // walk right and left and add to our intensity
      walker = peak_idx;
      // if we moved past the end of the spectrum, we need to try the last peak
      // of the spectrum (it could still be within the window)
      if (peak_idx >= input.size())
      {
        walker = input.size() - 1;
      }

      // add the current peak if it is between right and left
      if (input[walker].getMZ() > left && input[walker].getMZ() < right)
      {
        weight =  1 - fabs(input[walker].getMZ() - mz) / half_window_size;
        integrated_intensity += input[walker].getIntensity() * weight;
      }

      // (i) Walk to the left one step and then keep walking left until we go
      // outside the window. Note for the first step to the left we have to
      // check for the walker becoming zero.
      walker = peak_idx;
      if (walker > 0)
      {
        walker--;
        // special case: walker is now zero
        if (walker == 0 && input[walker].getMZ() > left && input[walker].getMZ() < right)
        {
          integrated_intensity += input[walker].getIntensity();
        }
      }
      while (walker > 0 && input[walker].getMZ() > left && input[walker].getMZ() < right)
      {
        weight =  1 - fabs(input[walker].getMZ() - mz) / half_window_size;
        integrated_intensity += input[walker].getIntensity() * weight; walker--;
      }

      // (ii) Walk to the right one step and then keep walking right until we
      // are outside the window
      walker = peak_idx;
      if (walker < input.size() )
      {
        walker++;
      }
      while (walker<input.size() && input[walker].getMZ()> left &&  input[walker].getMZ() < right)
      {
        weight = 1 - fabs(input[walker].getMZ() - mz) / half_window_size;
        integrated_intensity += input[walker].getIntensity() * weight; walker++;
      }
    }

private:

    template <typename TransitionExpT>
    static String extract_id_(TransitionExpT& transition_exp_used, String id);
    
    template <class SpectrumSettingsT, class ChromatogramT>
    void prepareSpectra_(SpectrumSettingsT& settings, std::vector<ChromatogramT>& chromatograms, OpenMS::TargetedExperiment& transition_exp)
    {
      // first prepare all the spectra (but leave them empty)
      for (Size i = 0; i < transition_exp.getTransitions().size(); i++)
      {
        const ReactionMonitoringTransition* transition = &transition_exp.getTransitions()[i];

        // 1) and 2) Extract precursor m/z and isolation window
        ChromatogramT chrom;
        Precursor prec;
        prec.setMZ(transition->getPrecursorMZ());
        if (settings.getPrecursors().size() > 0)
        {
          prec.setIsolationWindowLowerOffset(settings.getPrecursors()[0].getIsolationWindowLowerOffset());
          prec.setIsolationWindowUpperOffset(settings.getPrecursors()[0].getIsolationWindowUpperOffset());
        }

        // 3) set precursor peptide sequence / compound id in more-or-less default field
        String pepref = transition->getPeptideRef();
        for (Size pep_idx = 0; pep_idx < transition_exp.getPeptides().size(); pep_idx++)
        {
          const OpenMS::TargetedExperiment::Peptide* pep = &transition_exp.getPeptides()[pep_idx];
          if (pep->id == pepref)
          {
            prec.setMetaValue("peptide_sequence", pep->sequence);
            break;
          }
        }
        String compref = transition->getCompoundRef();
        for (Size comp_idx = 0; comp_idx < transition_exp.getCompounds().size(); comp_idx++)
        {
          const OpenMS::TargetedExperiment::Compound* comp = &transition_exp.getCompounds()[comp_idx];
          if (comp->id == compref)
          {
            prec.setMetaValue("peptide_sequence", String(comp->id) );
            break;
          }
        }

        // add precursor to spectrum
        chrom.setPrecursor(prec);

        // 4) Create product and set its m/z
        Product prod;
        prod.setMZ(transition->getProductMZ());
        chrom.setProduct(prod);

        // 5) Set the rest of the meta-data
        chrom.setInstrumentSettings(settings.getInstrumentSettings());
        chrom.setAcquisitionInfo(settings.getAcquisitionInfo());
        chrom.setSourceFile(settings.getSourceFile());

        for (Size j = 0; j < settings.getDataProcessing().size(); ++j)
        {
          settings.getDataProcessing()[j]->setMetaValue("performed_on_spectra", "true");
          chrom.getDataProcessing().push_back(settings.getDataProcessing()[j]);
        }

        // Set the id of the chromatogram, using the id of the transition (this gives directly the mapping of the two)
        chrom.setNativeID(transition->getNativeID());
        chrom.setChromatogramType(ChromatogramSettings::SELECTED_REACTION_MONITORING_CHROMATOGRAM);
        chromatograms.push_back(chrom);
      }

    }

    bool outsideExtractionWindow_(const ReactionMonitoringTransition& transition, double current_rt,
                                   const TransformationDescription& trafo, double rt_extraction_window);

    int getFilterNr_(const String& filter);

    void populatePeptideRTMap_(OpenMS::TargetedExperiment& transition_exp, double rt_extraction_window);

    std::map<OpenMS::String, double> PeptideRTMap_;

  };
    
  // Specialization for template (LightTargetedExperiment)
  template<>
  inline String ChromatogramExtractor::extract_id_<OpenSwath::LightTargetedExperiment>(OpenSwath::LightTargetedExperiment& transition_exp_used, String id)
  {
    OpenSwath::LightCompound comp = transition_exp_used.getCompoundByRef(id);
    if (!comp.sequence.empty())
    {
      return comp.sequence;
    }
    else
    {
      return comp.compound_name;
    }
  }


  // Specialization for template (TargetedExperiment)
  template<>
  inline String ChromatogramExtractor::extract_id_<OpenMS::TargetedExperiment>(OpenMS::TargetedExperiment& transition_exp_used, String id)
  {
    if (transition_exp_used.hasPeptide(id))
    {
      TargetedExperiment::Peptide p = transition_exp_used.getPeptideByRef(id);
      return p.sequence;
    }
    else if (transition_exp_used.hasCompound(id))
    {
      TargetedExperiment::Compound c = transition_exp_used.getCompoundByRef(id);
      return c.id;
    }
    else
    {
      return "";
    }
  }
}