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/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 {
 
      // 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)
        {}
 
        UInt32 protein_idx;            
        uint32_t mod_bitmask_;         
        std::pair<uint16_t , uint16_t> 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_{};     
      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;
 
    bool isBuild() const;
 
    const std::vector<Peptide>& getPeptides() 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 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 constexpr uint32_t SNES_KIND_BIT_MASK = 1u << 31;   
    static constexpr uint32_t SNES_SLOT_MASK = ~SNES_KIND_BIT_MASK; 
 
    enum class SnesAnchor
    {
      NONE,        
      PROT_NTERM,  
      PROT_CTERM   
    };
 
    static bool isSingleCMother(uint32_t mod_bitmask) noexcept
    {
      return (mod_bitmask & SNES_KIND_BIT_MASK) != 0;
    }
    static bool isSingleNMother(uint32_t mod_bitmask) noexcept
    {
      return (mod_bitmask & SNES_KIND_BIT_MASK) == 0;
    }
 
    bool isSnesMode() const noexcept { return is_snes_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 realizeSNESLength(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() = default;
      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 generateSNESMothers_(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> computeSnesSigmaDeltaSet_(bool include_prot_nterm_mods,
                                                   bool include_prot_cterm_mods) 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};
 
    bool is_snes_mode_{false};
 
    bool snes_enabled_{false};
 
    std::vector<double> snes_sigma_delta_set_;
    std::vector<double> snes_sigma_delta_set_with_prot_nterm_;
    std::vector<double> snes_sigma_delta_set_with_prot_cterm_;
 
    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_;
 
    void generateFragmentsLightweight_(
      std::vector<Fragment>& fragments,
      std::vector<Fragment>& 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_mod_masses) const;
 
    void generateFragmentsForSeries_(
      std::vector<Fragment>& fragments,
      const char* sequence,
      size_t seq_len,
      UInt32 peptide_idx,
      double n_term_mod_mass,
      double c_term_mod_mass,
      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<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 querySpectrumSNES_(const MSSpectrum& spectrum,
                            const std::vector<FASTAFile::FASTAEntry>& fasta_entries,
                            SpectrumMatchesTopN& sms);
 
    void queryPeaks(SpectrumMatchesTopN& candidates,
                   const MSSpectrum& spectrum,
                   const std::pair<size_t, size_t>& candidates_range,
                   const int16_t isotope_error,
                   const uint16_t precursor_charge,
                   const bool with_electron_ions);
    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 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;
 
 
  };
 
}