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-2018.
//
// 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/METADATA/PeptideIdentification.h>
#include <OpenMS/METADATA/PeptideEvidence.h>
#include <OpenMS/METADATA/ProteinIdentification.h>
#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/FORMAT/FASTAFile.h>

#include <algorithm>
#include <climits>
#include <vector>
#include <set>
#include <map>

namespace OpenMS
{
  class OPENMS_DLLAPI IDFilter
  {
public:

    IDFilter();

    virtual ~IDFilter();



    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 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
      {
        std::set<String> present_accessions = hit.extractProteinAccessionsSet();
        for (std::set<String>::iterator it = present_accessions.begin();
             it != present_accessions.end(); ++it)
        {
          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())
        {
          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())
          {
            LOG_WARN << "Peptide accession not available! Skipping Evidence." << std::endl;
          }
          else
          {
            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 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 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 filterHitsBySignificance(std::vector<IdentificationType>& ids,
                                         double threshold_fraction = 1.0)
    {
      for (typename std::vector<IdentificationType>::iterator id_it =
             ids.begin(); id_it != ids.end(); ++id_it)
      {
        double threshold_score = (threshold_fraction *
                                  id_it->getSignificanceThreshold());
        struct HasGoodScore<typename IdentificationType::HitType> score_filter(
          threshold_score, id_it->isHigherScoreBetter());
        keepMatchingItems(id_it->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 (typename std::vector<IdentificationType>::iterator id_it =
             ids.begin(); id_it != ids.end(); ++id_it)
      {
        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 (typename std::vector<IdentificationType>::iterator id_it =
             ids.begin(); id_it != ids.end(); ++id_it)
      {
        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 filterHitsBySignificance(PeakMap& experiment,
                                         double peptide_threshold_fraction,
                                         double protein_threshold_fraction)
    {
      // filter protein hits:
      filterHitsBySignificance(experiment.getProteinIdentifications(),
                               protein_threshold_fraction);
      // 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)
      {
        filterHitsBySignificance(exp_it->getPeptideIdentifications(),
                                 peptide_threshold_fraction);
        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);
    }

    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());
        }
      }
    }



  };

} // namespace OpenMS