Biosaur2Algorithm.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: Timo Sachsenberg $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/KERNEL/FeatureMap.h>
#include <OpenMS/KERNEL/MSExperiment.h>
 
#include <map>
#include <vector>
 
namespace OpenMS
{
 
class OPENMS_DLLAPI Biosaur2Algorithm :
  public DefaultParamHandler
{
public:
  struct Hill
  {
    std::vector<Size> scan_indices;        
    std::vector<Size> peak_indices;        
    std::vector<double> mz_values;         
    std::vector<double> intensities;       
    std::vector<double> rt_values;         
    std::vector<double> drift_times;       
    std::vector<double> ion_mobilities;    
    double mz_weighted_mean = 0.0;         
    double rt_start = 0.0;                 
    double rt_end = 0.0;                   
    double rt_apex = 0.0;                  
    double intensity_apex = 0.0;           
    double intensity_sum = 0.0;            
    double drift_time_median = -1.0;       
    double ion_mobility_median = -1.0;     
    Size length = 0;                       
    Size hill_idx = 0;                     
  };
 
  struct IsotopeCandidate
  {
    Size hill_idx = 0;              
    Size isotope_number = 0;        
    double mass_diff_ppm = 0.0;     
    double cos_corr = 0.0;          
  };
 
  struct PeptideFeature
  {
    double mz = 0.0;                              
    double rt_start = 0.0;                        
    double rt_end = 0.0;                          
    double rt_apex = 0.0;                         
    double intensity_apex = 0.0;                  
    double intensity_sum = 0.0;                   
    int charge = 0;                               
    Size n_isotopes = 0;                          
    Size n_scans = 0;                             
    double mass_calib = 0.0;                      
    double drift_time = -1.0;                     
    double ion_mobility = -1.0;                   
    std::vector<IsotopeCandidate> isotopes;       
    Size mono_hill_idx = 0;                       
  };
 
  Biosaur2Algorithm();
 
  void setMSData(const MSExperiment& ms_data);
 
  void setMSData(MSExperiment&& ms_data);
 
  MSExperiment& getMSData();
 
  const MSExperiment& getMSData() const;
 
  void run(FeatureMap& feature_map);
 
  void run(FeatureMap& feature_map,
           std::vector<Hill>& hills,
           std::vector<PeptideFeature>& peptide_features);
 
  void writeTSV(const std::vector<PeptideFeature>& features, const std::string& filename) const;
 
  void writeHills(const std::vector<Hill>& hills, const std::string& filename) const;
 
protected:
  void updateMembers_() override;
 
private:
 
 
  struct FastHillEntry
  {
    Size hill_index = 0;
    Size first_scan = 0;
    Size last_scan = 0;
  };
 
  struct PatternCandidate
  {
    Size mono_index = 0;                          
    double mono_mz = 0.0;                         
    int charge = 0;                               
    double cos_cor_isotopes = 0.0;                
    std::vector<IsotopeCandidate> isotopes;       
    Size n_scans = 0;                             
  };
 
 
 
 
  double calculatePPM_(double mz1, double mz2) const;
 
  double calculateMedian_(const std::vector<double>& values) const;
 
  double cosineCorrelation1D_(const std::vector<double>& v1,
                             const std::vector<double>& v2) const;
 
  std::pair<double, Size> checkingCosCorrelationForCarbon_(const std::vector<double>& theor_full,
                                                           const std::vector<double>& exp_full,
                                                           double thresh) const;
 
  std::pair<std::vector<double>, Size> computeAveragine_(double neutral_mass,
                                                         double apex_intensity) const;
 
  std::vector<double> meanFilter_(const std::vector<double>& data, Size window) const;
 
  std::pair<double, double> calibrateMass_(const std::vector<double>& mass_errors, double bin_width = 0.05) const;
 
  double computeHillMzStep_(const MSExperiment& exp,
                            double htol_ppm,
                            double min_intensity,
                            double min_mz,
                            double max_mz) const;
 
  void processTOF_(MSExperiment& exp) const;
 
  void centroidProfileSpectra_(MSExperiment& exp) const;
 
  void centroidPASEFData_(MSExperiment& exp, double mz_step, double pasef_tolerance) const;
 
  std::vector<Hill> detectHills_(const MSExperiment& exp, double htol_ppm, double min_intensity, double min_mz, double max_mz, bool use_im, Size min_length, std::vector<double>* hill_mass_diffs = nullptr) const;
 
  struct ScanPrep
  {
    std::vector<int> ids;                       
    std::vector<int> bins;                      
    std::vector<std::pair<int, int>> by_bin;    
    int bin_lo = 0;                             
    std::vector<int> bin_start;                 
    std::vector<int> im_bin_per_peak;           
    bool use_im = false;
  };
 
  struct PeakRef { Size hill; Size scan; Size peak; };
 
  ScanPrep prepareScan_(const MSSpectrum& spectrum, double min_intensity, double min_mz, double max_mz,
                        double mz_step, bool use_im_global) const;
 
  void linkScanToHills_(const MSSpectrum& spectrum,
                        Size scan_idx,
                        double htol_ppm,
                        const ScanPrep& cur,
                        const ScanPrep* prev,
                        std::vector<PeakRef>& refs,
                        Size& hill_idx_counter,
                        std::vector<Size>& prev_peak_to_hill,
                        const MSSpectrum*& prev_spectrum_ptr,
                        std::vector<double>* hill_mass_diffs) const;
 
  std::vector<Hill> processHills_(std::vector<Hill> hills, Size min_length) const;
 
  std::vector<Hill> splitHills_(std::vector<Hill> hills, double hvf, Size min_length) const;
 
  Size checkIsotopeValleySplit_(const std::vector<IsotopeCandidate>& isotopes, const std::vector<Hill>& hills, double ivf) const;
 
  double buildFastMzLookup_(const std::vector<Hill>& hills,
                            bool use_im,
                            std::map<int, std::vector<FastHillEntry>>& hills_mz_fast,
                            std::vector<int>& hill_im_bins) const;
 
  std::vector<PatternCandidate> generateIsotopeCandidates_(const std::vector<Hill>& hills,
                                                           double itol_ppm,
                                                           int min_charge,
                                                           int max_charge,
                                                           double ivf,
                                                           double mz_step,
                                                           const std::map<int, std::vector<FastHillEntry>>& hills_mz_fast,
                                                           const std::vector<Size>& hill_idx_to_index,
                                                           const std::vector<int>& hill_im_bins,
                                                           bool use_im) const;
 
  std::vector<PatternCandidate> applyRtFiltering_(const std::vector<PatternCandidate>& candidates,
                                                  const std::vector<Hill>& hills,
                                                  const std::vector<Size>& hill_idx_to_index) const;
 
  std::map<int, std::pair<double, double>> refineIsotopeCalibration_(const std::vector<PatternCandidate>& candidates,
                                                                      double itol_ppm,
                                                                      bool enable_isotope_calib) const;
 
  std::vector<PatternCandidate> filterByCalibration_(const std::vector<PatternCandidate>& candidates,
                                                     const std::vector<Hill>& hills,
                                                     const std::vector<Size>& hill_idx_to_index,
                                                     const std::map<int, std::pair<double, double>>& isotope_calib_map_ready,
                                                     bool enable_isotope_calib) const;
 
  std::vector<PeptideFeature> selectNonOverlappingPatterns_(const std::vector<PatternCandidate>& filtered_ready,
                                                            const std::vector<Hill>& hills,
                                                            bool negative_mode,
                                                            int iuse,
                                                            double itol_ppm) const;
 
  std::vector<PeptideFeature> detectIsotopePatterns_(std::vector<Hill>& hills, double itol_ppm, int min_charge, int max_charge, bool negative_mode, double ivf, int iuse, bool enable_isotope_calib, bool use_im) const;
 
  FeatureMap convertToFeatureMap_(const std::vector<PeptideFeature>& features,
                                  const std::vector<Hill>& hills) const;
 
  void debugCheckIsotopeConsistency_(const char* stage_label,
                                     double mono_mz_center,
                                     double mono_rt_apex,
                                     Size mono_hill_idx,
                                     int charge,
                                     double itol_ppm,
                                     const Hill& iso_hill,
                                     Size isotope_number) const;
 
  double cosineCorrelation_(const std::vector<double>& intensities1, const std::vector<Size>& scans1,
                            const std::vector<double>& intensities2, const std::vector<Size>& scans2) const;
 
  bool shouldThrowForMissingIM_(const MSSpectrum& spectrum) const;
 
  void processFAIMSGroup_(double faims_cv,
                          MSExperiment& group_exp,
                          double original_paseftol,
                          std::vector<Hill>& hills_out,
                          std::vector<PeptideFeature>& features_out);
 
 
  MSExperiment ms_data_;      
 
  // Algorithm parameters (cached from Param for performance)
  double mini_;               
  double minmz_;              
  double maxmz_;              
  double htol_;               
  double itol_;               
  double hvf_;                
  double ivf_;                
  Size minlh_;                
  int cmin_;                  
  int cmax_;                  
  double pasefmini_;          
  Size pasefminlh_;           
  int iuse_;                  
  bool negative_mode_;        
  bool tof_mode_;             
  bool profile_mode_;         
  bool use_hill_calib_;       
  bool ignore_iso_calib_;     
  double paseftol_;           
  double hrttol_;             
  std::string convex_hull_mode_;   
  bool faims_merge_features_; 
};
 
} // namespace OpenMS