FeatureFinderIdentificationAlgorithm.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: Hendrik Weisser $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/ANALYSIS/TARGETED/TargetedExperiment.h>
#include <OpenMS/ANALYSIS/MAPMATCHING/TransformationDescription.h>
#include <OpenMS/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.h>
#include <OpenMS/CONCEPT/ProgressLogger.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/FEATUREFINDER/FeatureFinderAlgorithmPickedHelperStructs.h>
#include <OpenMS/METADATA/PeptideIdentificationList.h>
 
#include <vector>
#include <fstream>
#include <map>
 
namespace OpenMS {
class OPENMS_DLLAPI FeatureFinderIdentificationAlgorithm :
  public DefaultParamHandler
{
public:
  FeatureFinderIdentificationAlgorithm();
 
  void run(
    PeptideIdentificationList peptides,
    const std::vector<ProteinIdentification>& proteins,
    FeatureMap& features,
    const FeatureMap& seeds = FeatureMap(),
    const std::string& spectra_file = ""
    );
 
  void runOnCandidates(FeatureMap& features);
 
  PeakMap& getMSData();
  const PeakMap& getMSData() const;
 
  void setMSData(const PeakMap& ms_data); // for pyOpenMS
  void setMSData(PeakMap&& ms_data); // moves peak data and saves the copy. Note that getMSData() will give back a processed/modified version.
 
  PeakMap& getChromatograms();
  const PeakMap& getChromatograms() const;
 
  ProgressLogger& getProgressLogger();
  const ProgressLogger& getProgressLogger() const;
 
  TargetedExperiment& getLibrary();
  const TargetedExperiment& getLibrary() const;
 
protected:
 
  typedef FeatureFinderAlgorithmPickedHelperStructs::MassTrace MassTrace;
  typedef FeatureFinderAlgorithmPickedHelperStructs::MassTraces MassTraces;
 
  typedef std::multimap<double, PeptideIdentification*> RTMap;
  typedef std::map<Int, RTMap> ChargeMap;
  typedef std::map<AASequence, ChargeMap> PeptideMap;
  typedef std::map<std::string, RTMap> PeptideRefRTMap;
 
  PeptideMap peptide_map_;
 
  Size n_peptides_ = 0; 
 
  Size batch_size_; 
  double rt_window_; 
  double mz_window_; 
  bool mz_window_ppm_; 
 
  double mapping_tolerance_; 
 
  double isotope_pmin_; 
  Size n_isotopes_; 
 
  double peak_width_;
  double min_peak_width_;
  double signal_to_noise_;
 
  std::string elution_model_;
 
  // output file (before filtering)
  std::string candidates_out_;
 
  FeatureMap* candidates_sink_ = nullptr;
 
  Size debug_level_;
 
  void updateMembers_() override;
 
  void storeCandidates_(const FeatureMap& candidates) const;
 
  struct RTRegion
  {
    double start, end;
    ChargeMap ids; 
  };
 
  struct IMStats
  {
    double median = -1.0; 
    double min = -1.0;    
    double max = -1.0;    
  };
 
  struct FeatureFilterQuality
  {
    bool operator()(const Feature& feature)
    {
      return feature.getOverallQuality() == 0.0;
    }
  } feature_filter_quality_;
 
  struct FeatureFilterPeptides
  {
    bool operator()(const Feature& feature)
    {
      return feature.getPeptideIdentifications().empty();
    }
  } feature_filter_peptides_;
 
  struct PeptideCompare
  {
    bool operator()(const PeptideIdentification& p1,
                    const PeptideIdentification& p2)
    {
      const std::string& seq1 = p1.getHits()[0].getSequence().toString();
      const std::string& seq2 = p2.getHits()[0].getSequence().toString();
      if (seq1 == seq2)
      {
        Int charge1 = p1.getHits()[0].getCharge();
        Int charge2 = p2.getHits()[0].getCharge();
        if (charge1 == charge2)
        {
          return p1.getRT() < p2.getRT();
        }
        return charge1 < charge2;
      }
      return seq1 < seq2;
    }
  } peptide_compare_;
 
  struct FeatureCompare
  {
    bool operator()(const Feature& f1, const Feature& f2)
    {
      const std::string ref1 = StringUtils::toStr(f1.getMetaValue("PeptideRef"));
      const std::string ref2 = StringUtils::toStr(f2.getMetaValue("PeptideRef"));
      if (ref1 == ref2)
      {
        return f1.getRT() < f2.getRT();
      }
      return ref1 < ref2;
    }
  } feature_compare_;
 
  PeakMap ms_data_; 
  PeakMap chrom_data_; 
  TargetedExperiment library_; 
  TargetedExperiment output_library_; 
 
  bool quantify_decoys_;
  double add_mass_offset_peptides_{0.0}; 
  double seed_apex_rt_tolerance_{5.0}; 
  bool use_psm_cutoff_;
  double psm_score_cutoff_;
  PeptideIdentificationList unassignedIDs_;
 
  const double seed_rt_window_ = 60.0; 
 
  std::map<std::string, double> isotope_probs_; 
  std::map<std::string, IMStats> im_stats_;
 
  IMStats global_im_stats_;
 
  MRMFeatureFinderScoring feat_finder_; 
 
  ProgressLogger prog_log_;
 
  void generateTransitions_(const std::string& peptide_id, double mz, Int charge,
                            const IsotopeDistribution& iso_dist);
 
  void addPeptideRT_(TargetedExperiment::Peptide& peptide, double rt) const;
 
  void getRTRegions_(ChargeMap& peptide_data, std::vector<RTRegion>& rt_regions, bool clear_IDs = true) const;
 
  IMStats getRTRegionIMStats_(const RTRegion& r);
 
  void calculateGlobalIMStats_();
 
  void annotateFeaturesFinalizeAssay_(
    FeatureMap& features,
    std::map<Size, std::vector<PeptideIdentification*> >& feat_ids,
    RTMap& rt_internal);
 
  void annotateFeatures_(FeatureMap& features, PeptideRefRTMap& ref_rt_map);
 
  void ensureConvexHulls_(Feature& feature) const;
 
  void postProcess_(FeatureMap& features, bool store_candidates = true);
 
  void initializeFeatureFinder_();
  double calculateRTWindow_() const;
  void removeSeedPseudoIDs_(FeatureMap& features);
 
  static bool isSeedPseudoHit_(const PeptideHit& hit);
  static bool isPseudoID_(const PeptideIdentification& pid);
 
  static bool isEmptyPseudoID_(const PeptideIdentification& pid);
 
  std::pair<double, double> calculateRTBounds_(double rt_min, double rt_max) const;
 
  void statistics_(const FeatureMap& features) const;
 
  void createAssayLibrary_(const PeptideMap::iterator& begin, const PeptideMap::iterator& end, PeptideRefRTMap& ref_rt_map, bool clear_IDs = true);
 
  void addPeptideToMap_(PeptideIdentification& peptide,
    PeptideMap& peptide_map);
 
  void filterFeatures_(FeatureMap& features);
 
  void runSingleGroup_(
    PeptideIdentificationList peptides,
    const std::vector<ProteinIdentification>& proteins,
    FeatureMap& features,
    const FeatureMap& seeds,
    const std::string& spectra_file);
 
  // seeds for untargeted extraction
  Size addSeeds_(PeptideIdentificationList& peptides, const FeatureMap& seeds);
 
  // quant. decoys
  Size addOffsetPeptides_(PeptideIdentificationList& peptides, double offset);
 
  template <typename It>
  std::vector<std::pair<It,It>>
  chunk_(It range_from, It range_to, const std::ptrdiff_t batch_size)
  {
    /* Aliases, to make the rest of the code more readable. */
    using std::vector;
    using std::pair;
    using std::make_pair;
    using std::distance;
    using diff_t = std::ptrdiff_t;
 
    /* Total item number and batch_size size. */
    const diff_t total {distance(range_from, range_to)};
    const diff_t num {total / batch_size};
 
    vector<pair<It,It>> chunks(num);
 
    It batch_end {range_from};
 
    /* Use the 'generate' algorithm to create batches. */
    std::generate(begin(chunks), end(chunks), [&batch_end, batch_size]()
    {
      It batch_start {batch_end };
 
      std::advance(batch_end, batch_size);
      return make_pair(batch_start, batch_end);
    });
 
    /* The last batch_size's end must always be 'range_to'. */
    if (chunks.empty())
    {
      chunks.emplace_back(range_from, range_to);
    }
    else
    {
      chunks.back().second = range_to;
    }
 
    return chunks;
  }
}; // namespace OpenMS
} // namespace OpenMS