IDFilter.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: Mathias Walzer $
// $Authors: Nico Pfeifer, Mathias Walzer, Hendrik Weisser $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CHEMISTRY/ProteaseDigestion.h>
#include <OpenMS/CONCEPT/LogStream.h>
#include <OpenMS/FORMAT/FASTAFile.h>
#include <OpenMS/KERNEL/ConsensusMap.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/METADATA/AnnotatedMSRun.h>
#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/METADATA/ID/IdentificationData.h>
#include <OpenMS/ANALYSIS/ID/IDScoreSwitcherAlgorithm.h>
#include <OpenMS/METADATA/PeptideEvidence.h>
#include <OpenMS/METADATA/PeptideIdentification.h>
#include <OpenMS/METADATA/ProteinIdentification.h>
#include <OpenMS/config.h>
#include <algorithm>
#include <climits>
#include <functional>
#include <map>
#include <set>
#include <unordered_set>
#include <vector>
 
namespace OpenMS
{
  template<typename T>
  concept IsPeptideOrProteinIdentification = 
    std::is_same_v<T, PeptideIdentification> || std::is_same_v<T, ProteinIdentification>;
 
  template<typename T>
  concept IsFeatureOrConsensusMap = 
    std::is_same_v<T, FeatureMap> || std::is_same_v<T, ConsensusMap>;
 
  class OPENMS_DLLAPI IDFilter
  {
  public:
    IDFilter() = default;
 
    virtual ~IDFilter() = default;
 
    typedef std::map<Int, PeptideHit*> ChargeToPepHitP;
    typedef std::unordered_map<std::string, ChargeToPepHitP> SequenceToChargeToPepHitP;
    typedef std::map<std::string, SequenceToChargeToPepHitP> RunToSequenceToChargeToPepHitP;
 
 
    template<class HitType>
    struct HasGoodScore {
      typedef HitType argument_type; // for use as a predicate
 
      double score;
      bool higher_score_better;
 
      HasGoodScore(double score_, bool higher_score_better_) : score(score_), higher_score_better(higher_score_better_)
      {
      }
 
      bool operator()(const HitType& hit) const
      {
        if (higher_score_better)
        {
          return hit.getScore() >= score;
        }
        return hit.getScore() <= score;
      }
    };
 
    template<class HitType>
    struct HasMetaValue {
      typedef HitType argument_type; // for use as a predicate
 
      String key;
      DataValue value;
 
      HasMetaValue(const String& key_, const DataValue& value_) : key(key_), value(value_)
      {
      }
 
      bool operator()(const HitType& hit) const
      {
        DataValue found = hit.getMetaValue(key);
        if (found.isEmpty())
          return false; // meta value "key" not set
        if (value.isEmpty())
          return true; // "key" is set, value doesn't matter
        return found == value;
      }
    };
 
    template<class HitType>
    struct HasMaxMetaValue {
      typedef HitType argument_type; // for use as a predicate
 
      String key;
      double value;
 
      HasMaxMetaValue(const String& key_, const double& value_) : key(key_), value(value_)
      {
      }
 
      bool operator()(const HitType& hit) const
      {
        DataValue found = hit.getMetaValue(key);
        if (found.isEmpty())
          return false; // meta value "key" not set
        return double(found) <= value;
      }
    };
 
        template<class HitType>
        struct HasMinMetaValue 
        {
          typedef HitType argument_type; // for use as a predicate
 
          String key;
          double value;
 
          HasMinMetaValue(const String& key_, const double& value_) : 
            key(key_), 
            value(value_)
          {
          }
 
          bool operator()(const HitType& hit) const
          {
            DataValue found = hit.getMetaValue(key);
            if (found.isEmpty())
            {
              return false; // meta value "key" not set
            }          
            return static_cast<double>(found) >= value;
          }
        };
 
 
    template<class HitType>
    struct HasDecoyAnnotation 
    {
      typedef HitType argument_type; // for use as a predicate
 
      struct HasMetaValue<HitType> target_decoy, is_decoy;
 
      HasDecoyAnnotation() : 
        target_decoy("target_decoy", "decoy"),
        is_decoy("isDecoy", "true")
      {
      }
 
      bool operator()(const HitType& hit) const
      {
        // @TODO: this could be done slightly more efficiently by returning
        // false if the "target_decoy" meta value is "target" or "target+decoy",
        // without checking for an "isDecoy" meta value in that case
        return target_decoy(hit) || is_decoy(hit);
      }
    };
 
    template<class HitType>
    struct HasMatchingAccessionUnordered {
      typedef HitType argument_type; // for use as a predicate
 
      const std::unordered_set<String>& accessions;
 
      HasMatchingAccessionUnordered(const std::unordered_set<String>& accessions_) : 
        accessions(accessions_)
      {
      }
 
      bool operator()(const PeptideHit& hit) const
      {
        for (const auto& it : hit.extractProteinAccessionsSet())
        {
          if (accessions.count(it) > 0)
            return true;
        }
        return false;
      }
 
      bool operator()(const ProteinHit& hit) const
      {
        return (accessions.count(hit.getAccession()) > 0);
      }
 
      bool operator()(const PeptideEvidence& evidence) const
      {
        return (accessions.count(evidence.getProteinAccession()) > 0);
      }
    };
 
    template<class HitType>
    struct HasMatchingAccession {
      typedef HitType argument_type; // for use as a predicate
 
      const std::set<String>& accessions;
 
      HasMatchingAccession(const std::set<String>& accessions_) : accessions(accessions_)
      {
      }
 
      bool operator()(const PeptideHit& hit) const
      {
        for (const auto& it : hit.extractProteinAccessionsSet())
        {
          if (accessions.count(it) > 0)
            return true;
        }
        return false;
      }
 
      bool operator()(const ProteinHit& hit) const
      {
        return (accessions.count(hit.getAccession()) > 0);
      }
 
      bool operator()(const PeptideEvidence& evidence) const
      {
        return (accessions.count(evidence.getProteinAccession()) > 0);
      }
    };
 
    template<class HitType, class Entry>
    struct GetMatchingItems {
      typedef HitType argument_type;            // for use as a predicate
      typedef std::map<String, Entry*> ItemMap; // Store pointers to avoid copying data
      ItemMap items;
 
      GetMatchingItems(std::vector<Entry>& records)
      {
        for (typename std::vector<Entry>::iterator rec_it = records.begin(); rec_it != records.end(); ++rec_it)
        {
          items[getKey(*rec_it)] = &(*rec_it);
        }
      }
 
      GetMatchingItems()
      {
      }
 
      const String& getKey(const FASTAFile::FASTAEntry& entry) const
      {
        return entry.identifier;
      }
 
      bool exists(const HitType& hit) const
      {
        return items.count(getHitKey(hit)) > 0;
      }
 
      const String& getHitKey(const PeptideEvidence& p) const
      {
        return p.getProteinAccession();
      }
 
      const Entry& getValue(const PeptideEvidence& evidence) const
      {
        if (!exists(evidence))
        {
          throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Accession: '" + getHitKey(evidence) + "'. peptide evidence accession not in data");
        }
        return *(items.find(getHitKey(evidence))->second);
      }
    };
 
 
 
 
    struct HasMinPeptideLength;
 
    struct HasMinCharge;
 
    struct HasLowMZError;
 
    struct HasMatchingModification;
 
    struct HasMatchingSequence;
 
    struct HasNoEvidence;
 
 
    class PeptideDigestionFilter
    {
    private:
      EnzymaticDigestion& digestion_;
      Int min_cleavages_;
      Int max_cleavages_;
 
    public:
      typedef PeptideHit argument_type;
      PeptideDigestionFilter(EnzymaticDigestion& digestion, Int min, Int max) : digestion_(digestion), min_cleavages_(min), max_cleavages_(max)
      {
      }
 
      static inline Int disabledValue()
      {
        return -1;
      }
 
      bool operator()(PeptideHit& p) const
      {
        const auto& fun = [&](const Int missed_cleavages) {
          bool max_filter = max_cleavages_ != disabledValue() ? missed_cleavages > max_cleavages_ : false;
          bool min_filter = min_cleavages_ != disabledValue() ? missed_cleavages < min_cleavages_ : false;
          return max_filter || min_filter;
        };
        return digestion_.filterByMissedCleavages(p.getSequence().toUnmodifiedString(), fun);
      }
 
      void filterPeptideSequences(std::vector<PeptideHit>& hits)
      {
        hits.erase(std::remove_if(hits.begin(), hits.end(), (*this)), hits.end());
      }
    };
 
 
    struct DigestionFilter {
      typedef PeptideEvidence argument_type;
 
      // Build an accession index to avoid the linear search cost
      GetMatchingItems<PeptideEvidence, FASTAFile::FASTAEntry> accession_resolver_;
      ProteaseDigestion& digestion_;
      bool ignore_missed_cleavages_;
      bool methionine_cleavage_;
 
      DigestionFilter(std::vector<FASTAFile::FASTAEntry>& entries, ProteaseDigestion& digestion, bool ignore_missed_cleavages, bool methionine_cleavage) :
          accession_resolver_(entries), digestion_(digestion), ignore_missed_cleavages_(ignore_missed_cleavages), methionine_cleavage_(methionine_cleavage)
      {
      }
 
      bool operator()(const PeptideEvidence& evidence) const
      {
        if (!evidence.hasValidLimits())
        {
          OPENMS_LOG_WARN << "Invalid limits! Peptide '" << evidence.getProteinAccession() << "' not filtered" << std::endl;
          return true;
        }
 
        if (accession_resolver_.exists(evidence))
        {
          return digestion_.isValidProduct(AASequence::fromString(accession_resolver_.getValue(evidence).sequence), evidence.getStart(), evidence.getEnd() - evidence.getStart(),
                                           ignore_missed_cleavages_, methionine_cleavage_);
        }
        else
        {
          if (evidence.getProteinAccession().empty())
          {
            OPENMS_LOG_WARN << "Peptide accession not available! Skipping Evidence." << std::endl;
          }
          else
          {
            OPENMS_LOG_WARN << "Peptide accession '" << evidence.getProteinAccession() << "' not found in fasta file!" << std::endl;
          }
          return true;
        }
      }
 
      void filterPeptideEvidences(std::vector<PeptideIdentification>& peptides)
      {
        IDFilter::FilterPeptideEvidences<IDFilter::DigestionFilter>(*this, peptides);
      }
    };
 
 
 
 
    template<class IdentificationType>
    struct HasNoHits {
      typedef IdentificationType argument_type; // for use as a predicate
 
      bool operator()(const IdentificationType& id) const
      {
        return id.getHits().empty();
      }
    };
 
 
 
 
    struct HasRTInRange;
 
    struct HasMZInRange;
 
 
 
 
    template<class Container, class Predicate>
    static void removeMatchingItems(Container& items, const Predicate& pred)
    {
      items.erase(std::remove_if(items.begin(), items.end(), pred), items.end());
    }
 
    template<class Container, class Predicate>
    static void keepMatchingItems(Container& items, const Predicate& pred)
    {
      items.erase(std::remove_if(items.begin(), items.end(), std::not_fn(pred)), items.end());
    }
 
    template<class Container, class Predicate>
    static void moveMatchingItems(Container& items, const Predicate& pred, Container& target)
    {
      auto part = std::partition(items.begin(), items.end(), std::not_fn(pred));
      std::move(part, items.end(), std::back_inserter(target));
      items.erase(part, items.end());
    }
 
    template<class IDContainer, class Predicate>
    static void removeMatchingItemsUnroll(IDContainer& items, const Predicate& pred)
    {
      for (auto& item : items)
      {
        removeMatchingItems(item.getHits(), pred);
      }
    }
 
    template<class IDContainer, class Predicate>
    static void keepMatchingItemsUnroll(IDContainer& items, const Predicate& pred)
    {
      for (auto& item : items)
      {
        keepMatchingItems(item.getHits(), pred);
      }
    }
 
    template<class MapType, class Predicate>
    static void keepMatchingPeptideHits(MapType& prot_and_pep_ids, Predicate& pred)
    {
      for (auto& feat : prot_and_pep_ids)
      {
        keepMatchingItemsUnroll(feat.getPeptideIdentifications(), pred);
      }
      keepMatchingItemsUnroll(prot_and_pep_ids.getUnassignedPeptideIdentifications(), pred);
    }
 
    template<class MapType, class Predicate>
    static void removeMatchingPeptideHits(MapType& prot_and_pep_ids, Predicate& pred)
    {
      for (auto& feat : prot_and_pep_ids)
      {
        removeMatchingItemsUnroll(feat.getPeptideIdentifications(), pred);
      }
      removeMatchingItemsUnroll(prot_and_pep_ids.getUnassignedPeptideIdentifications(), pred);
    }
 
    template<IsFeatureOrConsensusMap MapType, class Predicate>
    static void removeMatchingPeptideIdentifications(MapType& prot_and_pep_ids, Predicate& pred)
    {
      for (auto& feat : prot_and_pep_ids)
      {
        removeMatchingItems(feat.getPeptideIdentifications(), pred);
      }
      removeMatchingItems(prot_and_pep_ids.getUnassignedPeptideIdentifications(), pred);
    }
 
 
 
 
    template<class IdentificationType>
    static Size countHits(const std::vector<IdentificationType>& ids)
    {
      Size counter = 0;
      for (typename std::vector<IdentificationType>::const_iterator id_it = ids.begin(); id_it != ids.end(); ++id_it)
      {
        counter += id_it->getHits().size();
      }
      return counter;
    }
 
    template<class IdentificationType>
    static bool getBestHit(const std::vector<IdentificationType>& identifications, bool assume_sorted, typename IdentificationType::HitType& best_hit)
    {
      if (identifications.empty())
        return false;
 
      typename std::vector<IdentificationType>::const_iterator best_id_it = identifications.end();
      typename std::vector<typename IdentificationType::HitType>::const_iterator best_hit_it;
 
      for (typename std::vector<IdentificationType>::const_iterator id_it = identifications.begin(); id_it != identifications.end(); ++id_it)
      {
        if (id_it->getHits().empty())
          continue;
 
        if (best_id_it == identifications.end()) // no previous "best" hit
        {
          best_id_it = id_it;
          best_hit_it = id_it->getHits().begin();
        }
        else if (best_id_it->getScoreType() != id_it->getScoreType())
        {
          throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Can't compare scores of different types", best_id_it->getScoreType() + "/" + id_it->getScoreType());
        }
 
        bool higher_better = best_id_it->isHigherScoreBetter();
        for (typename std::vector<typename IdentificationType::HitType>::const_iterator hit_it = id_it->getHits().begin(); hit_it != id_it->getHits().end(); ++hit_it)
        {
          if ((higher_better && (hit_it->getScore() > best_hit_it->getScore())) || (!higher_better && (hit_it->getScore() < best_hit_it->getScore())))
          {
            best_hit_it = hit_it;
          }
          if (assume_sorted)
            break; // only consider the first hit
        }
      }
 
      if (best_id_it == identifications.end())
      {
        return false; // no hits in any IDs
      }
 
      best_hit = *best_hit_it;
      return true;
    }
 
    static void extractPeptideSequences(const std::vector<PeptideIdentification>& peptides, std::set<String>& sequences, bool ignore_mods = false);
 
    static std::map<String, std::vector<ProteinHit>> extractUnassignedProteins(ConsensusMap& cmap);
 
    template<class EvidenceFilter>
    static void FilterPeptideEvidences(EvidenceFilter& filter, std::vector<PeptideIdentification>& peptides)
    {
      for (std::vector<PeptideIdentification>::iterator pep_it = peptides.begin(); pep_it != peptides.end(); ++pep_it)
      {
        for (std::vector<PeptideHit>::iterator hit_it = pep_it->getHits().begin(); hit_it != pep_it->getHits().end(); ++hit_it)
        {
          std::vector<PeptideEvidence> evidences;
          remove_copy_if(hit_it->getPeptideEvidences().begin(), hit_it->getPeptideEvidences().end(), back_inserter(evidences), std::not_fn(filter));
          hit_it->setPeptideEvidences(evidences);
        }
      }
    }
 
 
 
    static void removeUnreferencedProteins(ConsensusMap& cmap, bool include_unassigned);
 
    static void removeUnreferencedProteins(std::vector<ProteinIdentification>& proteins, const std::vector<PeptideIdentification>& peptides);
    static void removeUnreferencedProteins(ProteinIdentification& proteins, const std::vector<PeptideIdentification>& peptides);
 
    static void updateProteinReferences(std::vector<PeptideIdentification>& peptides, const std::vector<ProteinIdentification>& proteins, bool remove_peptides_without_reference = false);
 
    static void updateProteinReferences(ConsensusMap& cmap, bool remove_peptides_without_reference = false);
 
    static void updateProteinReferences(ConsensusMap& cmap, const ProteinIdentification& ref_run, bool remove_peptides_without_reference = false);
 
    static bool updateProteinGroups(std::vector<ProteinIdentification::ProteinGroup>& groups, const std::vector<ProteinHit>& hits);
 
    static void removeUngroupedProteins(const std::vector<ProteinIdentification::ProteinGroup>& groups, std::vector<ProteinHit>& hits);
 
 
 
    template<IsPeptideOrProteinIdentification IdentificationType>
    static void removeEmptyIdentifications(std::vector<IdentificationType>& ids)
    {
      struct HasNoHits<IdentificationType> empty_filter;
      removeMatchingItems(ids, empty_filter);
    }
 
    template<class IdentificationType>
    static void filterHitsByScore(std::vector<IdentificationType>& ids, double threshold_score)
    {
      for (typename std::vector<IdentificationType>::iterator id_it = ids.begin(); id_it != ids.end(); ++id_it)
      {
        struct HasGoodScore<typename IdentificationType::HitType> score_filter(threshold_score, id_it->isHigherScoreBetter());
        keepMatchingItems(id_it->getHits(), score_filter);
      }
    }
 
    template<class IdentificationType>
    static void filterHitsByScore(std::vector<IdentificationType>& ids, double threshold_score, IDScoreSwitcherAlgorithm::ScoreType score_type)
    {
      IDScoreSwitcherAlgorithm switcher;
      bool at_least_one_found = false;
      for (IdentificationType& id : ids)
      {
        if (switcher.isScoreType(id.getScoreType(), score_type))
        {
          struct HasGoodScore<typename IdentificationType::HitType> score_filter(threshold_score, id.isHigherScoreBetter());
          keepMatchingItems(id.getHits(), score_filter);
        }
        else
        {
          // If one assumes they are all the same in the vector, this could be done in the beginning.
          String metaval = switcher.findScoreType(id, score_type);
          if (!metaval.empty())
          {
            if (switcher.isScoreTypeHigherBetter(score_type))
            {
              struct HasMinMetaValue<typename IdentificationType::HitType> score_filter(metaval, threshold_score);
              keepMatchingItems(id.getHits(), score_filter);
            }
            else
            {
              struct HasMaxMetaValue<typename IdentificationType::HitType> score_filter(metaval, threshold_score);
              keepMatchingItems(id.getHits(), score_filter);
            }
            at_least_one_found = true;
          }
        }
      }
      if (!at_least_one_found) OPENMS_LOG_WARN << String("Warning: No hit with the given score_type found. All hits removed.") << std::endl;
    }
 
    static void filterGroupsByScore(std::vector<ProteinIdentification::ProteinGroup>& grps, double threshold_score, bool higher_better);
 
    template<class IdentificationType>
    static void filterHitsByScore(IdentificationType& id, double threshold_score)
    {
      struct HasGoodScore<typename IdentificationType::HitType> score_filter(threshold_score, id.isHigherScoreBetter());
      keepMatchingItems(id.getHits(), score_filter);
    }
 
    template<class IdentificationType>
    static void keepNBestHits(std::vector<IdentificationType>& ids, Size n)
    {
      for (typename std::vector<IdentificationType>::iterator id_it = ids.begin(); id_it != ids.end(); ++id_it)
      {
        id_it->sort();
        if (n < id_it->getHits().size())
          id_it->getHits().resize(n);
      }
    }
 
    template<class IdentificationType>
    static void filterHitsByRank(std::vector<IdentificationType>& ids, Size min_rank, Size max_rank)
    {
      for (auto& id : ids)
      {
        auto& hits = id.getHits();
        if (hits.empty()) continue;
 
        id.sort(); // Ensure hits are properly sorted
 
        // ignore max_rank?
        if (max_rank < min_rank) max_rank = hits.size();
 
        Size rank = 1;
        double last_score = hits.front().getScore();
 
        // Remove hits not within [min_rank, max_rank], while computing rank on the fly
        hits.erase(
          std::remove_if(hits.begin(), hits.end(),
            [&](const auto& hit) {
              if (hit.getScore() != last_score)
              {
                ++rank;
                last_score = hit.getScore();
              }
              return rank < min_rank || rank > max_rank;
            }),
          hits.end()
        );
      }
    }
    
    template<class IdentificationType>
    static void removeDecoyHits(std::vector<IdentificationType>& ids)
    {
      struct HasDecoyAnnotation<typename IdentificationType::HitType> decoy_filter;
      for (typename std::vector<IdentificationType>::iterator id_it = ids.begin(); id_it != ids.end(); ++id_it)
      {
        removeMatchingItems(id_it->getHits(), decoy_filter);
      }
    }
 
    template<class IdentificationType>
    static void removeHitsMatchingProteins(std::vector<IdentificationType>& ids, const std::set<String> accessions)
    {
      struct HasMatchingAccession<typename IdentificationType::HitType> acc_filter(accessions);
      for (auto& id_it : ids)
      {
        removeMatchingItems(id_it.getHits(), acc_filter);
      }
    }
 
    template<IsPeptideOrProteinIdentification IdentificationType>
    static void keepHitsMatchingProteins(IdentificationType& id, const std::set<String>& accessions)
    {
      struct HasMatchingAccession<typename IdentificationType::HitType> acc_filter(accessions);
      keepMatchingItems(id.getHits(), acc_filter);      
    }
 
    template<class IdentificationType>
    static void keepHitsMatchingProteins(std::vector<IdentificationType>& ids, const std::set<String>& accessions)
    {
      for (auto& id_it : ids) keepHitsMatchingProteins(id_it, accessions);
    }
 
 
 
 
    static void keepBestPeptideHits(std::vector<PeptideIdentification>& peptides, bool strict = false);
 
    static void filterPeptidesByLength(std::vector<PeptideIdentification>& peptides, Size min_length, Size max_length = UINT_MAX);
 
    static void filterPeptidesByCharge(std::vector<PeptideIdentification>& peptides, Int min_charge, Int max_charge);
 
    static void filterPeptidesByRT(std::vector<PeptideIdentification>& peptides, double min_rt, double max_rt);
 
    static void filterPeptidesByMZ(std::vector<PeptideIdentification>& peptides, double min_mz, double max_mz);
 
    static void filterPeptidesByMZError(std::vector<PeptideIdentification>& peptides, double mass_error, bool unit_ppm);
 
 
    template<class Filter>
    static void filterPeptideEvidences(Filter& filter, std::vector<PeptideIdentification>& peptides);
 
    static void filterPeptidesByRTPredictPValue(std::vector<PeptideIdentification>& peptides, const String& metavalue_key, double threshold = 0.05);
 
    static void removePeptidesWithMatchingModifications(std::vector<PeptideIdentification>& peptides, const std::set<String>& modifications);
 
    static void removePeptidesWithMatchingRegEx(std::vector<PeptideIdentification>& peptides, const String& regex);
 
    static void keepPeptidesWithMatchingModifications(std::vector<PeptideIdentification>& peptides, const std::set<String>& modifications);
 
    static void removePeptidesWithMatchingSequences(std::vector<PeptideIdentification>& peptides, const std::vector<PeptideIdentification>& bad_peptides, bool ignore_mods = false);
 
    static void keepPeptidesWithMatchingSequences(std::vector<PeptideIdentification>& peptides, const std::vector<PeptideIdentification>& good_peptides, bool ignore_mods = false);
 
    static void keepUniquePeptidesPerProtein(std::vector<PeptideIdentification>& peptides);
 
    static void removeDuplicatePeptideHits(std::vector<PeptideIdentification>& peptides, bool seq_only = false);
 
 
 
 
    static void filterHitsByScore(AnnotatedMSRun& annotated_data,
                                  double peptide_threshold_score,
                                  double protein_threshold_score)
    {
      // filter protein hits:
      filterHitsByScore(annotated_data.getProteinIdentifications(),
                        protein_threshold_score);
      // don't remove empty protein IDs - they contain search meta data and may
      // be referenced by peptide IDs (via run ID)
 
      // filter peptide hits:
      for (PeptideIdentification& peptide_id : annotated_data.getPeptideIdentifications())
      {
        filterHitsByScore(peptide_id, peptide_threshold_score);
      }
      updateProteinReferences(annotated_data.getPeptideIdentifications(), annotated_data.getProteinIdentifications());
    }
 
    static void keepNBestHits(AnnotatedMSRun& annotated_data, Size n)
    {
      // don't filter the protein hits by "N best" here - filter the peptides
      // and update the protein hits!
      std::vector<PeptideIdentification> all_peptides; // IDs from all spectra
      // filter peptide hits:
      for (PeptideIdentification& peptide_id : annotated_data.getPeptideIdentifications())
      {
        // Create a temporary vector with a single PeptideIdentification
        std::vector<PeptideIdentification> temp_vec = {peptide_id};
        keepNBestHits(temp_vec, n);
        // Copy back the filtered hits
        if (!temp_vec.empty())
        {
          peptide_id = temp_vec[0];
        }
        else
        {
          peptide_id.getHits().clear();
        }
        
        // Since we're working with individual PeptideIdentifications, we don't need to remove empty ones
        // but we still need to update protein references
        temp_vec = {peptide_id};
        updateProteinReferences(temp_vec, annotated_data.getProteinIdentifications());
        all_peptides.push_back(peptide_id);
      }
      // update protein hits:
      removeUnreferencedProteins(annotated_data.getProteinIdentifications(), all_peptides);
    }
 
    static void keepNBestSpectra(std::vector<PeptideIdentification>& peptides, Size n);
 
    template<class MapType>
    static void keepNBestPeptideHits(MapType& map, Size n)
    {
      // The rank predicate needs annotated ranks, not sure if they are always updated. Use the following instead,
      // which sorts Hits first.
      for (auto& feat : map)
      {
        keepNBestHits(feat.getPeptideIdentifications(), n);
      }
      keepNBestHits(map.getUnassignedPeptideIdentifications(), n);
    }
 
    template<class MapType>
    static void removeEmptyIdentifications(MapType& prot_and_pep_ids)
    {
      const auto pred = HasNoHits<PeptideIdentification>();
      removeMatchingPeptideIdentifications(prot_and_pep_ids, pred);
    }
 
    static void keepBestPerPeptide(std::vector<PeptideIdentification>& pep_ids, bool ignore_mods, bool ignore_charges, Size nr_best_spectrum)
    {
      annotateBestPerPeptide(pep_ids, ignore_mods, ignore_charges, nr_best_spectrum);
      HasMetaValue<PeptideHit> best_per_peptide {"best_per_peptide", 1};
      keepMatchingItemsUnroll(pep_ids, best_per_peptide);
    }
 
    static void keepBestPerPeptidePerRun(std::vector<ProteinIdentification>& prot_ids, std::vector<PeptideIdentification>& pep_ids, bool ignore_mods, bool ignore_charges, Size nr_best_spectrum)
    {
      annotateBestPerPeptidePerRun(prot_ids, pep_ids, ignore_mods, ignore_charges, nr_best_spectrum);
      HasMetaValue<PeptideHit> best_per_peptide {"best_per_peptide", 1};
      keepMatchingItemsUnroll(pep_ids, best_per_peptide);
    }
 
    // TODO allow skipping unassigned?
    template<class MapType>
    static void annotateBestPerPeptidePerRun(MapType& prot_and_pep_ids, bool ignore_mods, bool ignore_charges, Size nr_best_spectrum)
    {
      const auto& prot_ids = prot_and_pep_ids.getProteinIdentifications();
 
      RunToSequenceToChargeToPepHitP best_peps_per_run;
      for (const auto& idrun : prot_ids)
      {
        best_peps_per_run[idrun.getIdentifier()] = SequenceToChargeToPepHitP();
      }
 
      for (auto& feat : prot_and_pep_ids)
      {
        annotateBestPerPeptidePerRunWithData(best_peps_per_run, feat.getPeptideIdentifications(), ignore_mods, ignore_charges, nr_best_spectrum);
      }
 
      annotateBestPerPeptidePerRunWithData(best_peps_per_run, prot_and_pep_ids.getUnassignedPeptideIdentifications(), ignore_mods, ignore_charges, nr_best_spectrum);
    }
 
    template<class MapType>
    static void keepBestPerPeptidePerRun(MapType& prot_and_pep_ids, bool ignore_mods, bool ignore_charges, Size nr_best_spectrum)
    {
      annotateBestPerPeptidePerRun(prot_and_pep_ids, ignore_mods, ignore_charges, nr_best_spectrum);
      HasMetaValue<PeptideHit> best_per_peptide {"best_per_peptide", 1};
      keepMatchingPeptideHits(prot_and_pep_ids, best_per_peptide);
    }
 
    static void annotateBestPerPeptidePerRun(const std::vector<ProteinIdentification>& prot_ids, std::vector<PeptideIdentification>& pep_ids, bool ignore_mods, bool ignore_charges,
                                             Size nr_best_spectrum)
    {
      RunToSequenceToChargeToPepHitP best_peps_per_run;
      for (const auto& id : prot_ids)
      {
        best_peps_per_run[id.getIdentifier()] = SequenceToChargeToPepHitP();
      }
      annotateBestPerPeptidePerRunWithData(best_peps_per_run, pep_ids, ignore_mods, ignore_charges, nr_best_spectrum);
    }
 
    static void annotateBestPerPeptidePerRunWithData(RunToSequenceToChargeToPepHitP& best_peps_per_run, std::vector<PeptideIdentification>& pep_ids, bool ignore_mods, bool ignore_charges,
                                                     Size nr_best_spectrum)
    {
      for (auto& pep : pep_ids)
      {
        SequenceToChargeToPepHitP& best_pep = best_peps_per_run[pep.getIdentifier()];
        annotateBestPerPeptideWithData(best_pep, pep, ignore_mods, ignore_charges, nr_best_spectrum);
      }
    }
 
    static void annotateBestPerPeptide(std::vector<PeptideIdentification>& pep_ids, bool ignore_mods, bool ignore_charges, Size nr_best_spectrum)
    {
      SequenceToChargeToPepHitP best_pep;
      for (auto& pep : pep_ids)
      {
        annotateBestPerPeptideWithData(best_pep, pep, ignore_mods, ignore_charges, nr_best_spectrum);
      }
    }
 
    static void annotateBestPerPeptideWithData(SequenceToChargeToPepHitP& best_pep, PeptideIdentification& pep, bool ignore_mods, bool ignore_charges, Size nr_best_spectrum)
    {
      bool higher_score_better = pep.isHigherScoreBetter();
      // make sure that first = best hit
      pep.sort();
 
      auto pepIt = pep.getHits().begin();
      auto pepItEnd = nr_best_spectrum == 0 || pep.getHits().size() <= nr_best_spectrum ? pep.getHits().end() : pep.getHits().begin() + nr_best_spectrum;
      for (; pepIt != pepItEnd; ++pepIt)
      {
        PeptideHit& hit = *pepIt;
 
        String lookup_seq;
        if (ignore_mods)
        {
          lookup_seq = hit.getSequence().toUnmodifiedString();
        }
        else
        {
          lookup_seq = hit.getSequence().toString();
        }
 
        int lookup_charge = 0;
        if (!ignore_charges)
        {
          lookup_charge = hit.getCharge();
        }
 
        // try to insert
        auto it_inserted = best_pep.emplace(std::move(lookup_seq), ChargeToPepHitP());
        auto it_inserted_chg = it_inserted.first->second.emplace(lookup_charge, &hit);
 
        PeptideHit*& p = it_inserted_chg.first->second; // now this gets either the old one if already present, or this
        if (!it_inserted_chg.second)                    // was already present -> possibly update
        {
          if ((higher_score_better && (hit.getScore() > p->getScore())) || (!higher_score_better && (hit.getScore() < p->getScore())))
          {
            p->setMetaValue("best_per_peptide", 0);
            hit.setMetaValue("best_per_peptide", 1);
            p = &hit;
          }
          else // note that this was def. not the best
          {
            // TODO if it is only about filtering, we can omit writing this metavalue (absence = false)
            hit.setMetaValue("best_per_peptide", 0);
          }
        }
        else // newly inserted -> first for that sequence (and optionally charge)
        {
          hit.setMetaValue("best_per_peptide", 1);
        }
      }
    }
 
    static void keepHitsMatchingProteins(
      AnnotatedMSRun& experiment,
      const std::vector<FASTAFile::FASTAEntry>& proteins)
    {
      std::set<String> accessions;
      for (auto it = proteins.begin(); it != proteins.end(); ++it)
      {
        accessions.insert(it->identifier);
      }
 
      // filter protein hits:
      keepHitsMatchingProteins(experiment.getProteinIdentifications(), accessions);
 
      // filter peptide hits:
      // std::pair<OpenMS::MSSpectrum&, OpenMS::PeptideIdentification&>
      for (auto [spectrum, peptide_id] : experiment)
      {
        if (spectrum.getMSLevel() == 2)
        {
          keepHitsMatchingProteins(peptide_id, accessions);
        }
      }
      removeEmptyIdentifications(experiment.getPeptideIdentifications());
    }
 
 
 
 
    static void keepBestMatchPerObservation(IdentificationData& id_data, IdentificationData::ScoreTypeRef score_ref);
 
    static void filterObservationMatchesByScore(IdentificationData& id_data, IdentificationData::ScoreTypeRef score_ref, double cutoff);
 
    static void removeDecoys(IdentificationData& id_data);
  };
 
} // namespace OpenMS