MapAlignmentAlgorithmIdentification.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-2022.
//
// 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: 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