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-2023.
//
// 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/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/KERNEL/StandardTypes.h>
#include <OpenMS/METADATA/ID/IdentificationData.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
{
  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 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) 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<class 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);
        }
      }
    }
 
 
 
 
    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(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<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);
      }
    }
 
    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)
    {
      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 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(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 metadata 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());
        // TODO super-duper inefficient.
        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);
    }
 
    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(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 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