FragmentIndex.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:  $
// $Authors:  $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CHEMISTRY/AASequence.h>
#include <OpenMS/CHEMISTRY/EnzymaticDigestion.h>
#include <OpenMS/CHEMISTRY/ResidueModification.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/DATASTRUCTURES/RegularExpression.h>
#include <OpenMS/FORMAT/FASTAFile.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/KERNEL/Peak1D.h>
 
 
#include <array>
#include <mutex>
#include <vector>
#include <functional>
#include <algorithm>   // std::max (used by inline static isOpenSearchMode)
 
namespace OpenMS
{
  class OPENMS_DLLAPI FragmentIndex : public DefaultParamHandler
  {
  public:
 
 
    struct Peptide {
 
      Peptide() {}
      // We need a constructor in order to emplace back
      Peptide(UInt32 protein_idx, uint32_t mod_bitmask, std::pair<uint16_t , uint16_t> sequence, float precursor_mz):
          protein_idx(protein_idx),
        mod_bitmask_(mod_bitmask),
        sequence_(sequence),
        precursor_mz_(precursor_mz)
        {}
 
        struct Span
        {
          Span() = default;
          Span(std::pair<uint16_t, uint16_t> span) : first(span.first), second(span.second) {}
          operator std::pair<uint16_t, uint16_t>() const { return {first, second}; }
          bool operator==(const Span&) const = default;
          uint16_t first;
          uint16_t second;
        };
 
        UInt32 protein_idx;            
        uint32_t mod_bitmask_;         
        Span sequence_;                
        float precursor_mz_;           
    };
 
    struct SpectrumMatch
    {
      uint32_t num_matched_{};       
      uint32_t subset_bitmask_{};    
      float    sigma_delta_{};       
      uint16_t precursor_charge_{};  
      int16_t  isotope_error_{};     
      uint16_t realized_length_{};   
      size_t   peptide_idx_{};       
    };
 
 
    struct SpectrumMatchesTopN
    {
      std::vector<SpectrumMatch> hits_;     
 
      SpectrumMatchesTopN() = default;
 
      SpectrumMatchesTopN& operator+=(const SpectrumMatchesTopN& other)
      {
        this->hits_.insert(this->hits_.end(), other.hits_.begin(), other.hits_.end());
        return *this;
      }
 
      void clear()
      {
        hits_.clear();
      }
    };
    FragmentIndex();
 
    ~FragmentIndex() override = default;
 
    FragmentIndex(const FragmentIndex&) = default;
    FragmentIndex& operator=(const FragmentIndex&) = default;
    FragmentIndex(FragmentIndex&&) = default;
    FragmentIndex& operator=(FragmentIndex&&) = default;
 
    bool isBuild() const;
 
    const std::vector<Peptide>& getPeptides() const;
 
    struct RemovedOccurrence
    {
      UInt32 peptide_idx;   
      UInt32 protein_idx;   
      uint16_t start;       
    };
 
    const std::vector<RemovedOccurrence>& getRemovedOccurrences() const noexcept { return removed_occurrences_; }
 
    bool hasProteinOccurrences(const std::vector<FASTAFile::FASTAEntry>& fasta_entries) const;
 
    void getProteinOccurrences(Size peptide_idx, const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                               std::vector<std::pair<UInt32, UInt32>>& occurrences) const;
 
    Size getNumFragments() const noexcept { return fi_fragments_.size() + electron_fragments_.size(); }
 
#ifdef DEBUG_FRAGMENT_INDEX
    void addSpecialPeptide(AASequence& peptide, Size source_idx);
#endif
 
    void build(const std::vector<FASTAFile::FASTAEntry> & fasta_entries);
 
    void build(const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
               const std::function<const MSExperiment*(Size)>& searched_spectra);
 
    void clear();
 
 
    std::pair<size_t, size_t> getPeptidesInMassWindow(float precursor_mass,
                                                      const std::pair<float, float>& window) const;
 
    static bool isOpenSearchMode(double lower_magnitude,
                                 double upper_magnitude,
                                 bool unit_ppm) noexcept
    {
      const double threshold = unit_ppm ? 1000.0 : 1.0;
      return std::max(lower_magnitude, upper_magnitude) > threshold;
    }
 
    static void restrictSearchParameters(Param& param);
 
    static void checkSearchParameters(const Param& param, const std::string& name_prefix = "");
 
    static void precursorCharges(const Precursor& precursor, int min_charge, int max_charge, std::vector<uint16_t>& charges);
 
    static void checkFixedModifications(const StringList& fixed_modifications, const std::string& name_prefix = "");
 
    static StringList shadowedVariableTerminalModifications(const StringList& fixed_modifications,
                                                            const StringList& variable_modifications);
 
    static StringList shadowedVariableResidueModifications(const StringList& fixed_modifications,
                                                           const StringList& variable_modifications);
 
 
    bool isPeptideMassMode() const noexcept { return is_peptide_mass_mode_; }
 
 
    struct Hit
    {
      Hit(UInt32 peptide_idx, float fragment_mz) :
        peptide_idx(peptide_idx),
        fragment_mz(fragment_mz)
      {}
      UInt32 peptide_idx; // index in database
      float fragment_mz;
    };
 
    std::vector<Hit> query(const Peak1D& peak,
                           const std::pair<size_t,size_t>& peptide_idx_range,
                           uint16_t peak_charge);
 
    void querySpectrum(const MSSpectrum& spectrum,
                       SpectrumMatchesTopN& sms);
 
    void querySpectrum(const MSSpectrum& spectrum,
                       const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                       SpectrumMatchesTopN& sms);
 
    void querySpectrum(const MSSpectrum& spectrum,
                       const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                       SpectrumMatchesTopN& sms,
                       bool with_electron_ions);
 
    AASequence reconstructModifiedSequence(const Peptide& peptide,
                                           const std::vector<FASTAFile::FASTAEntry>& fasta_entries) const;
 
    int realizePrefixLength(const Peptide& mother,
                            const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                            double target_mh_plus,
                            double tolerance_lower_magnitude,
                            double tolerance_upper_magnitude,
                            bool tolerance_ppm) const;
 
    AASequence reconstructRealizedSubSequence(const Peptide& mother,
                                              const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                                              size_t realized_length,
                                              uint32_t subset_bitmask = 0) const;
 
protected:
 
 
  struct Fragment
  {
      Fragment() {}
      Fragment(UInt32 peptide_idx, float fragment_mz):
          peptide_idx_(peptide_idx),
          fragment_mz_(fragment_mz)
      {}
      UInt32 peptide_idx_; // 32 bit in sage
      float fragment_mz_;
  };
 
    bool is_build_{false};              
 
    void updateMembers_() override;
 
    void generatePeptides(const std::vector<FASTAFile::FASTAEntry>& fasta_entries);
 
    void generateMothers_(const std::vector<FASTAFile::FASTAEntry>& fasta_entries);
 
    static void sortPeptides_(std::vector<Peptide>& peptides, size_t min_task_size = 8192);
 
    struct PeptidoformRendering_
    {
      std::vector<std::pair<const ResidueModification*, uint32_t>> mod_class;
      bool runs_hold_peptidoforms{true};
      bool unique_variable_renderings{true};
      bool slots_depend_on_context{false};
 
      uint32_t classOf(const ResidueModification* mod) const
      {
        for (const auto& [ptr, cls] : mod_class)
        {
          if (ptr == mod) return cls;
        }
        return 0;
      }
    };
 
    PeptidoformRendering_ peptidoformRendering_() const;
 
    void renderPeptidoform_(const Peptide& peptide, const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                            const PeptidoformRendering_& rendering, std::vector<uint32_t>& tokens) const;
 
    bool samePeptidoform_(const Peptide& a, const Peptide& b, const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                          const PeptidoformRendering_& rendering, std::vector<uint32_t>& tokens_a,
                          std::vector<uint32_t>& tokens_b) const;
 
    Size deduplicateByString_(const std::vector<FASTAFile::FASTAEntry>& fasta_entries);
 
    struct VarModEntry
    {
      double delta_mass;                    
      const ResidueModification* mod_ptr;   
      ResidueModification::TermSpecificity term_spec; 
    };
 
    struct ModSlot
    {
      uint16_t position;                    
      double delta_mass;                    
      const ResidueModification* mod_ptr;   
 
      static constexpr uint16_t NTERM_SLOT = UINT16_MAX - 1; 
      static constexpr uint16_t CTERM_SLOT = UINT16_MAX;      
    };
 
    static constexpr size_t MAX_MOD_SLOTS = 32; 
 
    void initModificationTables_();
 
    size_t buildModSlots_(const char* sequence, size_t seq_len, ModSlot* out_slots,
                          bool is_protein_nterm = false, bool is_protein_cterm = false) const;
 
    std::vector<double> computeSigmaDeltaSet_(bool include_prot_nterm_mods,
                                              bool include_prot_cterm_mods) const;
 
    std::vector<std::vector<double>> modShiftLevels_(bool include_prot_nterm_mods,
                                                     bool include_prot_cterm_mods) const;
 
    bool nonemptyZeroSigmaSubsetExists_() const;
 
    template <typename FragmentSink>
    void generatePeptidoformFragments_(FragmentSink& fragments, FragmentSink& electron_fragments, const char* sequence,
                                       size_t seq_len, UInt32 peptide_idx, uint32_t slot_bits, const ModSlot* slots,
                                       size_t n_slots, std::vector<double>& mod_masses) const;
 
    template <typename F>
    void forEachPrefixMass_(const char* sequence, size_t len, F&& f) const;
 
    std::array<double, 128> fixed_mod_deltas_{};
    std::array<const ResidueModification*, 128> fixed_mod_ptrs_{};
    double fixed_nterm_delta_{0.0};   
    double fixed_cterm_delta_{0.0};   
    const ResidueModification* fixed_nterm_mod_ptr_{nullptr};
    const ResidueModification* fixed_cterm_mod_ptr_{nullptr};
 
    std::array<std::vector<VarModEntry>, 128> variable_mod_table_{};
    std::vector<VarModEntry> variable_nterm_mods_;
    std::vector<VarModEntry> variable_cterm_mods_;
 
    bool mod_tables_initialized_{false};
 
    std::array<double, 128> fixed_residue_masses_{};
 
    bool is_peptide_mass_mode_{false};
 
    std::vector<double> sigma_delta_set_;
    std::vector<double> sigma_delta_set_with_prot_nterm_;
    std::vector<double> sigma_delta_set_with_prot_cterm_;
    bool zero_sigma_subsets_{false};
 
    struct ModShift_
    {
      double sigma;            
      bool protein_terminal;   
      size_t min_mods;         
    };
    std::vector<ModShift_> mod_shifts_;
    double min_mod_shift_{0.0};
    double max_mod_shift_{0.0};
 
    static constexpr uint32_t ENTRY_SITES_SHIFT_ = 28;
    static constexpr uint32_t ENTRY_MOTHER_MASK_ = (1u << ENTRY_SITES_SHIFT_) - 1;
    static constexpr uint32_t ENTRY_MAX_SITES_ = (1u << (32 - ENTRY_SITES_SHIFT_)) - 1;
    bool entry_sites_{false}; 
 
    static std::array<double, 128> residue_mass_table_;
    static std::once_flag mass_table_once_flag_;
    static void initResidueMassTable_();
 
    struct IonOffsets
    {
      double b_offset{0.0};
      double y_offset{0.0};
      double a_offset{0.0};
      double c_offset{0.0};
      double x_offset{0.0};
      double z_offset{0.0};
      double zp1_offset{0.0}; 
    };
    static IonOffsets ion_offsets_;
 
    template <typename FragmentSink>
    void generateFragmentsLightweight_(
      FragmentSink& fragments,
      FragmentSink& electron_fragments,
      const char* sequence,
      size_t seq_len,
      UInt32 peptide_idx,
      double n_term_mod_mass,
      double c_term_mod_mass,
      const double* residue_masses,
      const double* residue_mod_masses) const;
 
    template <typename FragmentSink>
    void generateFragmentsForSeries_(
      FragmentSink& fragments,
      const char* sequence,
      size_t seq_len,
      UInt32 peptide_idx,
      double n_term_mod_mass,
      double c_term_mod_mass,
      const double* residue_masses,
      const double* residue_mod_masses,
      bool add_b,
      bool add_a,
      bool add_c,
      bool add_y,
      bool add_x,
      bool add_z,
      bool add_zp1) const;
 
    std::vector<Peptide> fi_peptides_;   
    std::vector<RemovedOccurrence> removed_occurrences_; 
    std::vector<Fragment> fi_fragments_; 
    std::vector<Fragment> electron_fragments_;
    std::vector<float> electron_bucket_min_mz_; 
 
    std::vector<uint32_t> protein_lengths_;
 
    float fragment_min_mz_;  
    float fragment_max_mz_;  
    size_t min_ion_index_{0}; 
    size_t bucketsize_;       
    std::vector<float> bucket_min_mz_;  
    double precursor_mass_tolerance_lower_{20.0};   
    double precursor_mass_tolerance_upper_{20.0};   
    bool precursor_mass_tolerance_unit_ppm_{true};
    float fragment_mz_tolerance_;
    bool fragment_mz_tolerance_unit_ppm_{true};    
private:
 
 
    void querySpectrumPeptideMasses_(const MSSpectrum& spectrum,
                                     const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                                     SpectrumMatchesTopN& sms,
                                     bool with_electron_ions);
 
    struct CandidateBlock_
    {
      size_t first;           
      size_t second;
      size_t cell_offset;     
      int16_t isotope_error;  
    };
 
    void queryPeaks(SpectrumMatchesTopN& candidates,
                   const MSSpectrum& spectrum,
                   const std::vector<CandidateBlock_>& blocks,
                   const uint16_t precursor_charge,
                   const bool with_electron_ions);
 
    static constexpr size_t SKIP_STRIDE_ = 64;
    size_t skip_per_bucket_{0};
    std::vector<UInt32> bucket_skip_;
    std::vector<UInt32> electron_bucket_skip_;
    void buildSkipTables_();
 
    void searchDifferentPrecursorRanges(const MSSpectrum& spectrum,
                                        float precursor_mass,
                                        SpectrumMatchesTopN& sms,
                                        uint16_t charge,
                                        bool with_electron_ions);
 
    void sortAndBucketFragments_(std::vector<Fragment>& fragments,
                                 std::vector<float>& bucket_min_mz,
                                 int num_threads);
 
    void trimHits(SpectrumMatchesTopN& init_hits) const;
 
    //since we work with TheoreticalSpectrumGenerator, we must transfer some of those member variables
    bool add_b_ions_;
    bool add_y_ions_;
    bool add_a_ions_;
    bool add_c_ions_;
    bool add_x_ions_;
    bool add_z_ions_;
    bool add_zp1_ions_;
    bool electron_ions_{false}; 
 
    // SpectrumGenerator independend member variables
    std::string digestion_enzyme_;
    EnzymaticDigestion::Specificity enzyme_specificity_{EnzymaticDigestion::SPEC_FULL}; 
 
    bool isProteinNTerminal_(const std::string& protein, Size start) const;
 
    bool clip_nterm_methionine_ {false}; 
    size_t missed_cleavages_; 
    float peptide_min_mass_;
    float peptide_max_mass_;
    size_t peptide_min_length_;
    size_t peptide_max_length_;
  
    StringList modifications_fixed_;    
    StringList modifications_variable_; 
    size_t max_variable_mods_per_peptide_;
 
    // Search Related member variables
 
    uint16_t min_matched_peaks_;  
    int16_t min_isotope_error_;   
    int16_t max_isotope_error_;   
    uint16_t min_precursor_charge_; 
    uint16_t max_precursor_charge_; 
    uint16_t max_fragment_charge_;  
    uint32_t max_processed_hits_;   
    bool deduplicate_{false};       
    RegularExpression peptide_motif_; 
 
    bool matchesMotif_(const char* sequence, size_t length) const;
    
    bool isOpenSearchMode_() const noexcept
    {
      return isOpenSearchMode(precursor_mass_tolerance_lower_,
                              precursor_mass_tolerance_upper_,
                              precursor_mass_tolerance_unit_ppm_);
    }
 
    std::pair<float, float> computeMassWindow_(float precursor_mass) const;
 
    void keepPeptidesInPrecursorWindows_(const MSExperiment& spectra);
 
  };
 
}