IDFilter.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-2020.
//
// 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: Mathias Walzer $
// $Authors: Nico Pfeifer, Mathias Walzer, Hendrik Weisser $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/config.h>
#include <OpenMS/CHEMISTRY/ProteaseDigestion.h>
#include <OpenMS/CONCEPT/LogStream.h>
#include <OpenMS/FORMAT/FASTAFile.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/METADATA/PeptideIdentification.h>
#include <OpenMS/METADATA/PeptideEvidence.h>
#include <OpenMS/METADATA/ProteinIdentification.h>
#include <OpenMS/METADATA/ID/IdentificationData.h>
 
#include <algorithm>
#include <climits>
#include <vector>
#include <set>
#include <map>
#include <unordered_set>
 
namespace OpenMS
{
  class OPENMS_DLLAPI IDFilter
  {
public:
 
    IDFilter();
 
    virtual ~IDFilter();
 
    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 HasMaxRank
    {
      typedef HitType argument_type; // for use as a predicate
 
      Size rank;
 
      HasMaxRank(Size rank_):
        rank(rank_)
      {
        if (rank_ == 0)
        {
          throw Exception::IllegalArgument(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "The cut-off value for rank filtering must not be zero!");
        }
      }
 
      bool operator()(const HitType& hit) const
      {
        Size hit_rank = hit.getRank();
        if (hit_rank == 0)
        {
          throw Exception::MissingInformation(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "No rank assigned to peptide or protein hit");
        }
        return hit_rank <= rank;
      }
    };
 
    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 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)
      {
        return digestion_.filterByMissedCleavages(
          p.getSequence().toUnmodifiedString(),
          [&](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;
          });
      }
 
      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::not1(pred)),
                  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 <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);
 
    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::not1(filter));
          hit_it->setPeptideEvidences(evidences);
        }
      }
    }
 
 
 
 
    template <class IdentificationType>
    static void updateHitRanks(std::vector<IdentificationType>& ids)
    {
      for (typename std::vector<IdentificationType>::iterator it = ids.begin();
           it != ids.end(); ++it)
      {
        it->assignRanks();
      }
    }
 
    static void removeUnreferencedProteins(
      std::vector<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 bool updateProteinGroups(
      std::vector<ProteinIdentification::ProteinGroup>& groups,
      const std::vector<ProteinHit>& hits);
 
 
 
 
    template <class 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(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)
    {
      updateHitRanks(ids);
      if (min_rank > 1)
      {
        struct HasMaxRank<typename IdentificationType::HitType>
          rank_filter(min_rank - 1);
        for (typename std::vector<IdentificationType>::iterator id_it =
               ids.begin(); id_it != ids.end(); ++id_it)
        {
          removeMatchingItems(id_it->getHits(), rank_filter);
        }
      }
      if (max_rank >= min_rank)
      {
        struct HasMaxRank<typename IdentificationType::HitType>
          rank_filter(max_rank);
        for (typename std::vector<IdentificationType>::iterator id_it =
               ids.begin(); id_it != ids.end(); ++id_it)
        {
          keepMatchingItems(id_it->getHits(), rank_filter);
        }
      }
    }
 
    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 <class IdentificationType>
    static void keepHitsMatchingProteins(std::vector<IdentificationType>& ids,
                                         const std::set<String>& accessions)
    {
      struct HasMatchingAccession<typename IdentificationType::HitType> acc_filter(accessions);
      for (auto& id_it : ids)
      {
        keepMatchingItems(id_it.getHits(), acc_filter);
      }
    }
 
 
 
 
    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 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(PeakMap& experiment,
                                  double peptide_threshold_score,
                                  double protein_threshold_score)
    {
      // filter protein hits:
      filterHitsByScore(experiment.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 (PeakMap::Iterator exp_it = experiment.begin();
           exp_it != experiment.end(); ++exp_it)
      {
        filterHitsByScore(exp_it->getPeptideIdentifications(),
                          peptide_threshold_score);
        removeEmptyIdentifications(exp_it->getPeptideIdentifications());
        updateProteinReferences(exp_it->getPeptideIdentifications(),
                                experiment.getProteinIdentifications());
      }
      // @TODO: remove proteins that aren't referenced by peptides any more?
    }
 
    static void keepNBestHits(PeakMap& experiment, 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 (PeakMap::Iterator exp_it = experiment.begin();
           exp_it != experiment.end(); ++exp_it)
      {
        std::vector<PeptideIdentification>& peptides =
          exp_it->getPeptideIdentifications();
        keepNBestHits(peptides, n);
        removeEmptyIdentifications(peptides);
        updateProteinReferences(peptides,
                                experiment.getProteinIdentifications());
        all_peptides.insert(all_peptides.end(), peptides.begin(),
                            peptides.end());
      }
      // update protein hits:
      removeUnreferencedProteins(experiment.getProteinIdentifications(),
                                 all_peptides);
    }
 
    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)
    {
      removeMatchingPeptideHits(prot_and_pep_ids, HasNoHits<PeptideHit>());
    }
 
    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);
    }
 
    template <class MapType>
    static void keepBestPerPeptidePerRun(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);
 
      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(
      PeakMap& experiment,
      const std::vector<FASTAFile::FASTAEntry>& proteins)
    {
      std::set<String> accessions;
      for (std::vector<FASTAFile::FASTAEntry>::const_iterator it =
             proteins.begin(); it != proteins.end(); ++it)
      {
        accessions.insert(it->identifier);
      }
 
      // filter protein hits:
      keepHitsMatchingProteins(experiment.getProteinIdentifications(),
                               accessions);
      updateHitRanks(experiment.getProteinIdentifications());
 
      // filter peptide hits:
      for (PeakMap::Iterator exp_it = experiment.begin();
           exp_it != experiment.end(); ++exp_it)
      {
        if (exp_it->getMSLevel() == 2)
        {
          keepHitsMatchingProteins(exp_it->getPeptideIdentifications(),
                                   accessions);
          removeEmptyIdentifications(exp_it->getPeptideIdentifications());
          updateHitRanks(exp_it->getPeptideIdentifications());
        }
      }
    }
 
 
 
    static void keepBestMatchPerQuery(
      IdentificationData& id_data,
      IdentificationData::ScoreTypeRef score_ref);
 
    static void filterQueryMatchesByScore(
      IdentificationData& id_data,
      IdentificationData::ScoreTypeRef score_ref, double cutoff);
 
    static void removeDecoys(IdentificationData& id_data);
 
  };
 
} // namespace OpenMS