MapAlignmentAlgorithmIdentification.h

Go to the documentation of this file.

// Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Hendrik Weisser $
// $Authors: Eva Lange, Clemens Groepl, Hendrik Weisser $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/ANALYSIS/MAPMATCHING/TransformationDescription.h>
#include <OpenMS/CONCEPT/ProgressLogger.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/DATASTRUCTURES/ListUtils.h>
#include <OpenMS/KERNEL/ConsensusMap.h>
#include <OpenMS/KERNEL/FeatureMap.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/METADATA/PeptideIdentification.h>
#include <OpenMS/METADATA/ID/IdentificationData.h>
 
#include <cmath> // for "abs"
#include <limits> // for "max"
#include <map>
 
namespace OpenMS
{
  class OPENMS_DLLAPI MapAlignmentAlgorithmIdentification :
    public DefaultParamHandler,
    public ProgressLogger
  {
public:
    MapAlignmentAlgorithmIdentification();
 
    ~MapAlignmentAlgorithmIdentification() override;
 
    // Set a reference for the alignment
    template <typename DataType> void setReference(DataType& data)
    {
      reference_.clear();
      if (data.empty()) return; // empty input resets the reference
      SeqToList rt_data;
      // set these here because "checkParameters_" may not have been called yet:
      use_feature_rt_ = param_.getValue("use_feature_rt").toBool();
      score_cutoff_ = param_.getValue("score_cutoff").toBool();
      score_type_ = (std::string)param_.getValue("score_type");
      bool sorted = getRetentionTimes_(data, rt_data);
      computeMedians_(rt_data, reference_, sorted);
 
      if (reference_.empty())
      {
        throw Exception::MissingInformation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Could not extract retention time information from the reference file");
      }
    }
 
    template <typename DataType>
    void align(std::vector<DataType>& data,
               std::vector<TransformationDescription>& transformations,
               Int reference_index = -1)
    {
      checkParameters_(data.size());
      startProgress(0, 3, "aligning maps");
 
      reference_index_ = reference_index;
      // is reference one of the input files?
      bool use_internal_reference = (reference_index >= 0);
      if (use_internal_reference)
      {
        if (reference_index >= Int(data.size()))
        {
          throw Exception::IndexOverflow(__FILE__, __LINE__,
                                         OPENMS_PRETTY_FUNCTION,
                                         reference_index, data.size());
        }
        setReference(data[reference_index]);
      }
 
      // one set of RT data for each input map, except reference (if any):
      std::vector<SeqToList> rt_data(data.size() - use_internal_reference);
      bool all_sorted = true;
      for (Size i = 0, j = 0; i < data.size(); ++i)
      {
        if ((reference_index >= 0) && (i == Size(reference_index)))
        {
          continue; // skip reference map, if any
        }
        all_sorted &= getRetentionTimes_(data[i], rt_data[j++]);
      }
      setProgress(1);
 
      computeTransformations_(rt_data, transformations, all_sorted);
      setProgress(2);
 
      setProgress(3);
      endProgress();
    }
 
protected:
 
    typedef std::map<String, DoubleList> SeqToList;
 
    typedef std::map<String, double> SeqToValue;
 
    Int reference_index_;
 
    SeqToValue reference_;
 
    Size min_run_occur_;
 
    bool use_feature_rt_{};
 
    bool use_adducts_{};
 
    double min_score_;
 
    bool score_cutoff_{};
 
    String score_type_;
 
    bool (*better_) (double, double) = [](double, double) {return true;};
 
    void computeMedians_(SeqToList& rt_data, SeqToValue& medians,
                         bool sorted = false);
 
    bool getRetentionTimes_(std::vector<PeptideIdentification>& peptides,
                            SeqToList& rt_data);
 
    // "id_data" can't be "const" here or template resolution will fail
    bool getRetentionTimes_(IdentificationData& id_data, SeqToList& rt_data);
 
    bool getRetentionTimes_(PeakMap& experiment, SeqToList& rt_data);
 
    template <typename MapType>
    bool getRetentionTimes_(MapType& features, SeqToList& rt_data)
    {
      if (!score_cutoff_)
      {
        better_ = [](double, double)
        {return true;};
      }
      else if (features[0].getPeptideIdentifications()[0].isHigherScoreBetter())
      {
        better_ = [](double a, double b)
        { return a >= b; };
      }
      else
      {
        better_ = [](double a, double b)
        { return a <= b; };
      }
 
      for (typename MapType::Iterator feat_it = features.begin();
           feat_it != features.end(); ++feat_it)
      {
        if (use_feature_rt_)
        {
          // find the peptide ID closest in RT to the feature centroid:
          String sequence;
          double rt_distance = std::numeric_limits<double>::max();
          bool any_hit = false;
          for (std::vector<PeptideIdentification>::iterator pep_it =
                 feat_it->getPeptideIdentifications().begin(); pep_it !=
                 feat_it->getPeptideIdentifications().end(); ++pep_it)
          {
            if (!pep_it->getHits().empty())
            {
              any_hit = true;
              double current_distance = fabs(pep_it->getRT() -
                                             feat_it->getRT());
              if (current_distance < rt_distance)
              {
                pep_it->sort();
                if (better_(pep_it->getHits()[0].getScore(), min_score_))
                {
                  sequence = pep_it->getHits()[0].getSequence().toString();
                  rt_distance = current_distance;
                }
              }
            }
          }
 
          if (any_hit) rt_data[sequence].push_back(feat_it->getRT());
        }
        else
        {
          getRetentionTimes_(feat_it->getPeptideIdentifications(), rt_data);
        }
      }
 
      if (!use_feature_rt_ &&
          param_.getValue("use_unassigned_peptides").toBool())
      {
        getRetentionTimes_(features.getUnassignedPeptideIdentifications(),
                           rt_data);
      }
 
      // remove duplicates (can occur if a peptide ID was assigned to several
      // features due to overlap or annotation tolerance):
      for (SeqToList::iterator rt_it = rt_data.begin(); rt_it != rt_data.end();
           ++rt_it)
      {
        DoubleList& rt_values = rt_it->second;
        sort(rt_values.begin(), rt_values.end());
        DoubleList::iterator it = unique(rt_values.begin(), rt_values.end());
        rt_values.resize(it - rt_values.begin());
      }
      return true; // RTs were already sorted for duplicate detection
    }
 
    void computeTransformations_(std::vector<SeqToList>& rt_data,
                                 std::vector<TransformationDescription>&
                                 transforms, bool sorted = false);
 
    void checkParameters_(const Size runs);
 
    void getReference_();
 
    IdentificationData::ScoreTypeRef handleIdDataScoreType_(const IdentificationData& id_data);
 
private:
 
    MapAlignmentAlgorithmIdentification(const MapAlignmentAlgorithmIdentification&);
 
    MapAlignmentAlgorithmIdentification& operator=(const MapAlignmentAlgorithmIdentification&);
 
  };
 
} // namespace OpenMS