IDConflictResolverAlgorithm.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Hendrik Weisser $
// $Authors: Hendrik Weisser, Lucia Espona, Moritz Freidank $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/KERNEL/ConsensusMap.h>
#include <OpenMS/KERNEL/FeatureMap.h>
#include <OpenMS/METADATA/PeptideIdentification.h>
#include <OpenMS/METADATA/PeptideIdentificationList.h>
 
//-------------------------------------------------------------
// Doxygen docu
//-------------------------------------------------------------
 
 
 
namespace OpenMS
{
 
class OPENMS_DLLAPI IDConflictResolverAlgorithm
{
public:
  static void resolve(FeatureMap& features, bool keep_matching = false);
 
  static void resolve(ConsensusMap& features, bool keep_matching = false);
 
  static void resolveAllHitRankAggregation(FeatureMap& features);
 
  static void resolveAllHitRankAggregation(ConsensusMap& features);
 
  static void resolveBetweenFeatures(FeatureMap& features);
  
  static void resolveBetweenFeatures(ConsensusMap& features);
 
  struct OPENMS_DLLAPI UnresolvedIdentifications
  {
    Size removed = 0;
    Size multiply_identified_spectra = 0;
    Size without_spectrum_reference = 0;
    Size inconsistent_score_direction = 0;
    std::string example;
  };
 
  static UnresolvedIdentifications reduceToOnePerSpectrum(PeptideIdentificationList& ids);
 
protected:
 
  template<class T>
  static void resolveConflict_(T& map, bool keep_matching)
  {
    // annotate as not part of the resolution
    for (PeptideIdentification& p : map.getUnassignedPeptideIdentifications())
    {
      p.setMetaValue("feature_id", "not mapped"); // not mapped to a feature
    }
 
    for (auto& c : map)
    {
      c.setMetaValue("feature_id",StringUtils::toStr(c.getUniqueId()));
      if (!keep_matching)
      {
        resolveConflict_(c.getPeptideIdentifications(),
                         map.getUnassignedPeptideIdentifications(),
                         c.getUniqueId());
      }
      else
      {
        resolveConflictKeepMatching_(c.getPeptideIdentifications(),
                         map.getUnassignedPeptideIdentifications(),
                         c.getUniqueId());
      }
    }
  }
  
  // compare peptide IDs by score of best hit (hits must be sorted first!)
  // (note to self: the "static" is necessary to avoid cryptic "no matching
  // function" errors from gcc when the comparator is used below)
  static bool compareIDsSmallerScores_(const PeptideIdentification & left,
                          const PeptideIdentification & right);
 
  static void resolveConflict_(
    PeptideIdentificationList & peptides,
    PeptideIdentificationList & removed,
    UInt64 uid);
 
  static void resolveConflictKeepMatching_(
      PeptideIdentificationList & peptides,
      PeptideIdentificationList & removed,
      UInt64 uid);
 
  static void resolveAggregateConflict_(
      PeptideIdentificationList & peptides,
      PeptideIdentificationList & removed,
      UInt64 uid);
 
  template<class T>
  static void rankAggregation_(T& map)
  {
    // annotate unassigned IDs as not part of the resolution
    for (PeptideIdentification& p : map.getUnassignedPeptideIdentifications())
    {
      p.setMetaValue("feature_id", "not mapped"); // not mapped to a feature
    }
 
    for (auto& c : map)
    {
      c.setMetaValue("feature_id",StringUtils::toStr(c.getUniqueId()));
      resolveAggregateConflict_(c.getPeptideIdentifications(),
                                map.getUnassignedPeptideIdentifications(),
                                c.getUniqueId());
    }
  }
 
  template<class T>
  static void resolveBetweenFeatures_(T & map)
  {
    // unassigned peptide identifications in this map
    PeptideIdentificationList& unassigned = map.getUnassignedPeptideIdentifications();
    
    // A std::map tracking the set of unique features.
    // Uniqueness criterion/key is a pair <charge, sequence> for each feature. The peptide sequence may be modified, i.e. is not stripped.
    typedef std::map<std::pair<Int, AASequence>, typename T::value_type*> FeatureSet;
    FeatureSet feature_set;
    
    // Create a std::map `feature_set` mapping pairs <charge, sequence> to a pointer to
    // the feature with the highest intensity for this sequence.
    for (typename T::value_type& element : map)
    {
      PeptideIdentificationList& pep_ids = element.getPeptideIdentifications();
      
      if (!pep_ids.empty())
      {
        if (pep_ids.size() != 1)
        {
          // Should never happen. In IDConflictResolverAlgorithm TOPP tool
          // IDConflictResolverAlgorithm::resolve() is called before IDConflictResolverAlgorithm::resolveBetweenFeatures().
          throw OpenMS::Exception::IllegalArgument(__FILE__, __LINE__, __FUNCTION__, "Feature does contain multiple identifications.");
        }
        
        // Make sure best hit is in front, i.e. sort hits first.
        pep_ids.front().sort();
        const std::vector<PeptideHit>& hits = pep_ids.front().getHits();
        
        if (!hits.empty())
        {
          const PeptideHit& highest_score_hit = hits.front();
          
          // Pair <charge, sequence> of charge of the new feature and the sequence of its highest scoring peptide hit.
          std::pair<Int, AASequence> pair = std::make_pair(element.getCharge(), highest_score_hit.getSequence());
          
          // If a <charge, sequence> pair is not yet in the FeatureSet or new feature `feature_in_set`
          // has higher intensity than its counterpart `feature_set[<charge, sequence>]`
          // store a pointer to `feature_in_set` in `feature_set`.
          typename FeatureSet::iterator feature_in_set = feature_set.find(pair);
          if (feature_in_set != feature_set.end())
          {
            // Identical (charge, sequence) key found. Remove annotations from either the old or new feature.
            
            if (feature_in_set->second->getIntensity() < element.getIntensity())
            {
              // Remove annotations from the old low-intensity feature. But only after moving these annotations to the unassigned list.
              PeptideIdentificationList& obsolete = feature_in_set->second->getPeptideIdentifications();
              unassigned.insert(unassigned.end(), obsolete.begin(), obsolete.end());
              PeptideIdentificationList pep_ids_empty;
              feature_in_set->second->setPeptideIdentifications(pep_ids_empty);
              
              // Replace feature in the set.
              feature_in_set->second = &(element);
            }
            else
            {
              // Remove annotations from the new low-intensity feature. But only after moving these annotations to the unassigned list.
              PeptideIdentificationList& obsolete = element.getPeptideIdentifications();
              unassigned.insert(unassigned.end(), obsolete.begin(), obsolete.end());
              PeptideIdentificationList pep_ids_empty;
              element.setPeptideIdentifications(pep_ids_empty);
            }
          }
          else
          {
            // Feature is not yet in our set -- add it.
            feature_set[pair] = &(element);
          }
        }
      }
    }
  }
  
};
 
}// namespace OpenMS