ProSEAlgorithm.h

Go to the documentation of this file.

// Copyright (c) 2002-present, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer:  $
// $Authors: Raphael Förster $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/ProgressLogger.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
 
#include <OpenMS/ANALYSIS/ID/FragmentIndex.h>
#include <OpenMS/ANALYSIS/ID/FragmentIonLikelihoodModel.h>
#include <OpenMS/ANALYSIS/ID/OpenSearchModificationAnalysis.h>
#include <OpenMS/CHEMISTRY/EnzymaticDigestion.h>
#include <OpenMS/CHEMISTRY/ModifiedPeptideGenerator.h>
#include <OpenMS/FORMAT/FASTAFile.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/METADATA/PeptideIdentificationList.h>
 
#include <algorithm>   // std::min (used by inline computeModMatchTolerance_)
#include <functional>  // std::hash (ProteinMapping_)
#include <iosfwd>      // std::ostream (renderRunSummary / renderModificationSummary)
#include <map>
#include <string>      // std::string (renderRunSummaryJson return / manifest)
#include <string_view>
#include <unordered_map>
#include <unordered_set>
#include <utility>     // std::pair (renderRunSummaryJson manifest)
#include <vector>
 
namespace OpenMS
{
 
class TheoreticalSpectrumGenerator;
 
class OPENMS_DLLAPI ProSEAlgorithm :
  public DefaultParamHandler,
  public ProgressLogger
{
  public:
    ProSEAlgorithm(); 
 
    enum class ExitCodes
    {
      EXECUTION_OK,
      INPUT_FILE_EMPTY,
      UNEXPECTED_RESULT,
      UNKNOWN_ERROR,
      ILLEGAL_PARAMETERS
    };
 
    struct RunStatistics
    {
      std::string input_file;                  
      Size ms2_spectra = 0;                    
      Size matched_spectra = 0;                
      Size target_psms = 0;                    
      Size decoy_psms = 0;                     
      bool fdr_applied = false;                
      double achieved_psm_fdr = -1.0;          
      Size unique_peptides = 0;                
      Size unique_proteins = 0;                
      std::map<Int, Size> charge_histogram;    
      std::map<Size, Size> missed_cleavage_histogram; 
      bool score_stats_valid = false;          
      double hyperscore_min = 0.0;
      double hyperscore_median = 0.0;
      double hyperscore_max = 0.0;
      bool prec_tol_valid = false;             
      double prec_err_median = 0.0, prec_err_mad = 0.0, prec_err_recommended = 0.0;
      bool frag_tol_valid = false;             
      double frag_err_median = 0.0, frag_err_mad = 0.0, frag_err_recommended = 0.0;
      double seconds_search = 0.0;             
      double seconds_calibration = 0.0;        
      double seconds_fdr = 0.0;                
    };
 
    struct SharedSearchStats
    {
      std::string database_file;               
      std::string enzyme;
      double precursor_tol_lower = 0.0, precursor_tol_upper = 0.0;
      std::string precursor_tol_unit;
      double fragment_tol = 0.0;
      std::string fragment_tol_unit;
      Int min_charge = 0, max_charge = 0;
      Size missed_cleavages = 0;
      std::vector<std::string> fixed_mods, variable_mods, ion_series;
      bool open_search = false;
      bool calibration_enabled = false;
      bool peptide_mass_mode = false;
      bool chunked = false;
      std::string decoy_mode;                  
      double psm_fdr_threshold = 0.0, protein_fdr_threshold = 0.0;
      Size db_target_proteins = 0;             
      Size db_decoy_proteins = 0;              
      Size indexed_peptides = 0;               
      Size indexed_fragments = 0;              
      double seconds_index_build = 0.0;        
      double seconds_total = 0.0;              
    };
 
    struct SearchResult
    {
      ExitCodes exit_code = ExitCodes::EXECUTION_OK;
      std::vector<ProteinIdentification> protein_ids;
      PeptideIdentificationList peptide_ids;
      OpenSearchModificationAnalysis::OpenSearchAnalysisResult modification_analysis;
      bool is_open_search = false;
      RunStatistics stats;
      PeptideIdentificationList pre_fdr_peptide_ids; 
    };
 
    struct MultiFileSearchResult
    {
      std::vector<SearchResult> per_file;
      SearchResult aggregate;
 
      std::string decoy_string;
      bool decoy_is_prefix = true;
      bool have_decoys = false;
 
      SharedSearchStats shared;
    };
 
    struct SearchContext
    {
      std::vector<FASTAFile::FASTAEntry> db;
      FragmentIndex fragment_index;
 
      bool release_fragment_index_after_scoring = false;
 
      std::string decoy_string;
      bool decoy_is_prefix = true;
      bool have_decoys = false;
      bool electron_ions = false;
    };
 
    ExitCodes search(const std::string& in_spectra,
      const std::string& in_db,
      std::vector<ProteinIdentification>& prot_ids,
      PeptideIdentificationList& pep_ids) const;
 
    static void applyCompleteSetProteinFDR(std::vector<ProteinIdentification>& protein_ids,
                                           PeptideIdentificationList& peptide_ids,
                                           const std::string& decoy_string,
                                           bool decoy_is_prefix,
                                           double protein_fdr);
 
    void annotatePsmQValues(PeptideIdentificationList& peptide_ids) const;
 
    SearchResult searchWithModificationAnalysis(const std::string& in_spectra,
                                                const std::string& in_db,
                                                const std::string& output_base_name = "") const;
 
    ExitCodes search(PeakMap& spectra,
                     const std::vector<FASTAFile::FASTAEntry>& fasta_db,
                     std::vector<ProteinIdentification>& prot_ids,
                     PeptideIdentificationList& pep_ids) const;
 
    SearchContext prepareContext(const std::vector<FASTAFile::FASTAEntry>& fasta_db) const;
 
    SearchContext prepareContext(const std::vector<FASTAFile::FASTAEntry>& fasta_db, bool electron_ions) const;
 
    ExitCodes search(PeakMap& spectra,
                     SearchContext& ctx,
                     std::vector<ProteinIdentification>& prot_ids,
                     PeptideIdentificationList& pep_ids) const;
 
    SearchResult searchWithModificationAnalysis(PeakMap& spectra,
                                                const std::vector<FASTAFile::FASTAEntry>& fasta_db,
                                                const std::string& output_base_name = "") const;
 
    MultiFileSearchResult searchWithModificationAnalysis(
        const std::vector<std::string>& in_spectra_files,
        const std::vector<FASTAFile::FASTAEntry>& fasta_db,
        const std::vector<std::string>& output_base_names = {},
        const std::string& aggregate_base_name = "",
        bool build_pooled_aggregate = true) const;
 
    MultiFileSearchResult searchWithModificationAnalysis(
        const std::vector<std::string>& in_spectra_files,
        const std::string& in_db,
        const std::vector<std::string>& output_base_names = {},
        const std::string& aggregate_base_name = "",
        bool build_pooled_aggregate = true) const;
 
    void setKeepPreFdrPsms(bool keep)
    {
      keep_pre_fdr_psms_ = keep;
      if (!keep) { pre_fdr_peptide_ids_ = PeptideIdentificationList(); }
    }
 
  protected:
    void updateMembers_() override;
 
    struct AnnotatedHit_
    {
      AASequence sequence;
      /*
      std::string_view sequence;
      SignedSize peptide_mod_index; ///< enumeration index of the non-RNA peptide modification
      */
      // Layout: doubles first, then floats, then int, then uint16_t — minimizes padding (40 bytes excluding AASequence)
      double score = 0; 
      double delta_mass = 0.0; 
      float prefix_fraction = 0; 
      float suffix_fraction = 0; 
      float mean_error = 0.0f; 
      int isotope_error = 0; 
      uint16_t applied_charge = 0; 
      uint16_t matched_prefix_ions = 0; 
      uint16_t matched_suffix_ions = 0; 
      UInt32 index_peptide = 0; 
 
      static bool hasBetterScore(const AnnotatedHit_& a, const AnnotatedHit_& b)
      {
        if (a.score != b.score) return a.score > b.score;
        return a.sequence < b.sequence;
      }
    };
 
    struct CandidatePoolStats_
    {
      double sum = 0.0;         
      double sumsq = 0.0;       
      double best = 0.0;        
      double second_best = 0.0; 
      Size count = 0;           
      std::unordered_set<std::string> seen_candidates;
 
      void add(double score)
      {
        if (score > best) { second_best = best; best = score; }
        else if (score > second_best) { second_best = score; }
        sum += score;
        sumsq += score * score;
        ++count;
      }
 
      double deltaScore() const { return best - second_best; }
 
      double zScore() const;
    };
 
    struct DeisotopingSettings_
    {
      unsigned int min_peaks{3};          
      bool charge_cap_precursor{false};   
      bool sum_intensity{false};          
    };
 
    static void preprocessSpectra_(PeakMap& exp,
                                   double fragment_mass_tolerance,
                                   bool fragment_mass_tolerance_unit_ppm,
                                   bool deisotope_requested,
                                   Size peaks_keep_n,
                                   Int peaks_window_top,
                                   const std::string& window_type,
                                   const DeisotopingSettings_& deisotoping,
                                   FragmentIonLikelihoodModel::PeakLists* ion_evidence = nullptr,
                                   bool ion_evidence_scored_peaks = false,
                                   PeakMap* evidence_spectra = nullptr,
                                   PeakMap* query_spectra = nullptr);
 
    static void preprocessSpectra_(PeakMap& exp,
                                   double fragment_mass_tolerance,
                                   bool fragment_mass_tolerance_unit_ppm,
                                   bool deisotope_requested,
                                   Size peaks_keep_n,
                                   Int peaks_window_top,
                                   const std::string& window_type = "auto");
 
    void annotateIonPriors_(const PeakMap& spectra,
                            const FragmentIonLikelihoodModel::PeakLists& evidence,
                            std::vector<ProteinIdentification>& protein_ids,
                            PeptideIdentificationList& peptide_ids) const;
 
    static AASequence reversedNoiseSequence_(const AASequence& sequence);
 
    static void filterLocalPeaks_(MSSpectrum& spectrum, Size peaks_per_window);
 
    enum class DecoyMode_
    {
      AUTO,     
      GENERATE, 
      IGNORE    
    };
 
    struct DecoyStrategy_
    {
      bool generate{false};       
      bool strip_existing{false}; 
      bool have_decoys{false};    
      std::string decoy_string;   
      bool is_prefix{true};       
      std::string strip_string;   
      bool strip_is_prefix{true}; 
    };
 
    DecoyStrategy_ resolveDecoyStrategy_(
        const std::vector<FASTAFile::FASTAEntry>& db) const;
 
    Param fragmentIndexParameters_(bool electron_ions = false) const;
 
    struct SpectrumGenerators_;
 
    SpectrumGenerators_ spectrumGenerators_() const;
 
    struct LocalFragmentEvidence_
    {
      double chance_match_surprise = 0.0;
      double mass_competition_evidence = 0.0;
    };
 
    static std::vector<double> localPeakDensities_(const MSSpectrum& spectrum);
 
    static LocalFragmentEvidence_ localFragmentEvidence_(const MSSpectrum& spectrum,
                                                        const MSSpectrum& theoretical,
                                                        const std::vector<double>& densities,
                                                        double tolerance, bool ppm);
 
    static LocalFragmentEvidence_ localFragmentEvidence_(const MSSpectrum& spectrum,
                                                        const MSSpectrum& theoretical,
                                                        int max_charge,
                                                        const std::vector<double>& densities,
                                                        double tolerance, bool ppm,
                                                        std::vector<double>& alternatives);
 
    static bool isElectronActivated_(const MSSpectrum& spectrum);
 
    Size countElectronActivated_(const PeakMap& spectra) const;
 
    std::vector<FASTAFile::FASTAEntry> buildDecoyAugmentedDB_(
        const std::vector<FASTAFile::FASTAEntry>& fasta_db,
        const DecoyStrategy_& strategy) const;
 
    std::vector<FASTAFile::FASTAEntry> buildDecoyAugmentedDB_(
        std::vector<FASTAFile::FASTAEntry>&& fasta_db,
        const DecoyStrategy_& strategy) const;
 
    SearchContext prepareContext_(std::vector<FASTAFile::FASTAEntry>&& fasta_db, bool electron_ions,
                                  const std::function<const PeakMap*(Size)>& searched_spectra = {}) const;
 
    bool restrictIndexToSpectra_() const { return !calibration_enabled_; }
 
    std::vector<FASTAFile::FASTAEntry> buildCalibrationSample_(
        const std::vector<FASTAFile::FASTAEntry>& full_db) const;
 
    ExitCodes searchChunked_(
        PeakMap& spectra,
        std::vector<FASTAFile::FASTAEntry>& full_db,
        const DecoyStrategy_& strategy,
        std::vector<ProteinIdentification>& protein_ids,
        PeptideIdentificationList& peptide_ids) const;
 
    void scoreSpectraAgainstIndex_(
        const PeakMap& spectra,
        FragmentIndex& fi,
        const std::vector<FASTAFile::FASTAEntry>& db,
        const SpectrumGenerators_& generators,
        double effective_fragment_tol,
        bool fragment_mass_tolerance_unit_ppm,
        bool open_search_mode,
        std::vector<std::vector<AnnotatedHit_>>& annotated_hits,
        std::vector<CandidatePoolStats_>& pool_stats,
        const std::string& progress_label,
        const PeakMap* query_spectra = nullptr) const;
 
    void postProcessHits_(const PeakMap& exp,
      std::vector<std::vector<ProSEAlgorithm::AnnotatedHit_> >& annotated_hits,
      const std::vector<CandidatePoolStats_>& pool_stats,
      std::vector<ProteinIdentification>& protein_ids,
      PeptideIdentificationList& peptide_ids,
      Size top_hits,
      const StringList& modifications_fixed,
      const StringList& modifications_variable,
      Int peptide_missed_cleavages,
      double precursor_mass_tolerance,
      double fragment_mass_tolerance,
      const std::string& precursor_mass_tolerance_unit_ppm,
      const std::string& fragment_mass_tolerance_unit_ppm,
      const Int precursor_min_charge,
      const Int precursor_max_charge,
      const std::string& enzyme,
      const std::string& database_name,
      const PeakMap* evidence_spectra = nullptr) const;
 
    void annotatePeptDeepFeatures_(const PeakMap& spectra,
                                   std::vector<ProteinIdentification>& protein_ids,
                                   PeptideIdentificationList& peptide_ids) const;
 
    bool peptdeep_enable_{false};
    std::string peptdeep_ms2_model_;
    std::string peptdeep_rt_model_;
    std::string peptdeep_instrument_{"QE"};
    double peptdeep_nce_{-1.0};
    std::string peptdeep_rt_model_type_{"b_spline"};
    int peptdeep_batch_size_{500};
 
    mutable double precursor_mass_tolerance_lower_{10.0};   
    mutable double precursor_mass_tolerance_upper_{10.0};   
    std::string precursor_mass_tolerance_unit_{"ppm"};
 
    Size precursor_min_charge_;
    Size precursor_max_charge_;
 
    double fragment_mass_tolerance_;
 
    std::string fragment_mass_tolerance_unit_;
 
    bool deisotope_requested_{true};
    bool scoring_multiple_charges_{false}; 
    int scoring_max_charge_{2};            
    bool mass_accuracy_score_{false};      
    double mass_error_sd_ppm_{7.0};        
 
    int scoringMaxCharge_(int precursor_charge) const
    {
      return scoring_multiple_charges_ ? std::max(1, std::min(precursor_charge - 1, scoring_max_charge_)) : 1;
    }
    Size peaks_keep_n_{0};     
    Int peaks_window_top_{20}; 
    std::string peaks_window_type_ {"auto"}; 
    DeisotopingSettings_ deisotoping_;        
    bool query_raw_spectrum_{false}; 
 
    StringList modifications_fixed_;
 
    StringList modifications_variable_;
 
    Size modifications_max_variable_mods_per_peptide_;
 
    std::string enzyme_;
 
    DecoyMode_ decoy_mode_{DecoyMode_::AUTO};
    std::string decoy_prefix_;
 
    double fdr_psm_{0.0};
    bool fdr_psm_by_scored_charges_{true}; 
    double fdr_protein_{0.0};
 
    StringList annotate_psm_;
 
    bool self_trained_ion_priors_{false};       
    FragmentIonLikelihoodModel::ContextSet ion_prior_model_{FragmentIonLikelihoodModel::ContextSet::RICH}; 
    bool ion_prior_scored_peaks_{false};        
    int ion_prior_max_fragment_charge_{2};      
    double ion_prior_train_fdr_{0.01};          
    Size ion_prior_min_psms_{100};              
    bool ion_prior_feature_llr_{true};          
    bool ion_prior_feature_explained_{true};    
    bool ion_prior_feature_topk_{false};        
 
    Size peptide_min_size_;
    Size peptide_max_size_;
    Size peptide_missed_cleavages_;
    EnzymaticDigestion::Specificity peptide_enzyme_specificity_{EnzymaticDigestion::SPEC_FULL};
 
    std::string peptide_motif_;
 
    Size report_top_hits_;
    bool isotope_error_observed_minus_theoretical_{true};
 
    bool add_a_ions_{false};
    bool add_b_ions_{true};
    bool add_c_ions_{false};
    bool add_x_ions_{false};
    bool add_y_ions_{true};
    bool add_z_ions_{false};
    bool add_zp1_ions_{false};
    bool ions_by_activation_{true}; 
    bool protein_mapping_from_index_{true}; 
 
    Size database_chunk_size_{0};  
 
    bool calibration_enabled_{false};
    double calibration_subset_ratio_{0.1};
    Size calibration_min_psms_{50};
 
    struct CalibrationResult_
    {
      double precursor_shift{0};     
      double precursor_spread{0};    
      double cal_lower{0};           
      double cal_upper{0};           
      double fragment_tolerance{0};  
      double fragment_shift{0};      
      bool extreme_bias{false};      
      bool success{false};           
    };
 
    mutable CalibrationResult_ last_calibration_result_;
 
    mutable RunStatistics last_run_stats_;
 
    bool keep_pre_fdr_psms_ = false; 
    mutable PeptideIdentificationList pre_fdr_peptide_ids_;
 
    mutable double last_mod_match_tolerance_used_{-1.0};
 
    double computeModMatchTolerance_() const
    {
      if (precursor_mass_tolerance_lower_ <= 0.0) return precursor_mass_tolerance_upper_;
      if (precursor_mass_tolerance_upper_ <= 0.0) return precursor_mass_tolerance_lower_;
      return std::min(precursor_mass_tolerance_lower_, precursor_mass_tolerance_upper_);
    }
 
    CalibrationResult_ runCalibrationPass_(PeakMap& spectra,
                                           FragmentIndex& fragment_index,
                                           const std::vector<FASTAFile::FASTAEntry>& db,
                                           const PeakMap* query_spectra = nullptr) const;
 
    struct ProteinMapping_
    {
      static constexpr Size SHARDS = 64;
      std::vector<std::string_view> sequences;
      std::vector<std::vector<std::pair<UInt32, UInt32>>> occurrences;
      std::vector<std::unordered_map<std::string_view, Size>> entries;
      std::string fallback_reason;
 
      Size find(std::string_view sequence) const
      {
        if (entries.empty()) return sequences.size();
        const auto& table = entries[std::hash<std::string_view>{}(sequence) % SHARDS];
        const auto entry = table.find(sequence);
        return entry == table.end() ? sequences.size() : entry->second;
      }
    };
 
    static ProteinMapping_ buildProteinMapping_(const FragmentIndex& index,
                                                const std::vector<FASTAFile::FASTAEntry>& db,
                                                std::vector<UInt32> candidates,
                                                const Param& indexer_parameters);
 
    static bool applyProteinMapping_(const ProteinMapping_& mapping,
                                     const std::vector<FASTAFile::FASTAEntry>& db,
                                     const Param& indexer_parameters,
                                     std::vector<ProteinIdentification>& protein_ids,
                                     PeptideIdentificationList& peptide_ids,
                                     std::string& fallback_reason);
 
    static bool accessionHasDecoyMarker_(const std::string& accession,
                                         const std::string& marker, bool is_prefix);
 
    static void capturePreFdrStats_(const PeptideIdentificationList& peptide_ids,
                                    RunStatistics& stats);
 
    static double maxRetainedScore_(const PeptideIdentificationList& peptide_ids);
 
    void collectRunStatistics_(const PeakMap& spectra,
                               const std::vector<ProteinIdentification>& protein_ids,
                               const PeptideIdentificationList& peptide_ids,
                               RunStatistics& stats) const;
 
  public:
    static void updateFinalStats(RunStatistics& stats,
                                 const PeptideIdentificationList& peptide_ids,
                                 const std::string& enzyme,
                                 bool fdr_applied);
 
    static void renderRunSummary(const RunStatistics& stats,
                                 const SharedSearchStats& shared,
                                 const OpenSearchModificationAnalysis::OpenSearchAnalysisResult& mod_analysis,
                                 bool is_open_search,
                                 std::ostream& os);
 
    static void renderModificationSummary(const OpenSearchModificationAnalysis::OpenSearchAnalysisResult& mod_analysis,
                                          std::ostream& os);
 
    static std::string renderRunSummaryYaml(
        const MultiFileSearchResult& mfres,
        const std::vector<std::pair<std::string, std::vector<std::string>>>& manifest,
        Size files_failed,
        Size files_total);
 
  private:
 
    bool isOpenSearchMode_() const
    {
      return FragmentIndex::isOpenSearchMode(precursor_mass_tolerance_lower_,
                                             precursor_mass_tolerance_upper_,
                                             precursor_mass_tolerance_unit_ == "ppm");
    }
};
 
} // namespace