OpenSearchModificationAnalysis.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/METADATA/PeptideIdentificationList.h>
#include <OpenMS/DATASTRUCTURES/StringUtils.h>
 
#include <cmath>
#include <map>
#include <vector>
#include <unordered_map>
#include <unordered_set>
 
namespace OpenMS
{
  class OPENMS_DLLAPI OpenSearchModificationAnalysis
  {
  public:
 
    struct ModificationPattern
    {
      double count = 0.0;           
      std::vector<double> masses;   
      int num_charge_states = 0;    
    };
 
    struct ModificationSummary
    {
      int count;                    
      std::string name;                  
      int num_charge_states;        
      std::vector<double> masses;   
    };
 
    struct DeltaMassEntry
    {
      double delta_mass = 0.0;            
      int count = 0;                      
      int unique_peptides = 0;            
      int num_charge_states = 0;          
      double percentage = 0.0;            
      std::string mapped_modification;    
      bool is_known_modification = false; 
    };
 
    struct PTMEntry
    {
      std::string name;                        
      double theoretical_mass = 0.0;      
      double observed_mass = 0.0;         
      double mass_deviation = 0.0;        
      int count = 0;                      
      int unique_peptides = 0;            
      int num_charge_states = 0;          
      double percentage = 0.0;            
      std::map<char, int> residue_counts; 
      std::string target_residues;             
    };
 
    struct DeltaMassStatistics
    {
      std::vector<DeltaMassEntry> entries;  
      int total_psms = 0;                   
      int modified_psms = 0;                
      int unmodified_psms = 0;              
      int decoy_psms = 0;                   
      double mean_delta_mass = 0.0;         
      double median_delta_mass = 0.0;       
    };
 
    struct PTMStatistics
    {
      std::vector<PTMEntry> entries;        
      int total_modified_psms = 0;          
      int unknown_modification_psms = 0;    
      int num_unique_modifications = 0;     
    };
 
    struct OpenSearchAnalysisResult
    {
      DeltaMassStatistics delta_mass_stats; 
      PTMStatistics ptm_stats;              
      std::vector<ModificationSummary> summaries; 
    };
 
    struct FuzzyDoubleComparator
    {
      double epsilon;
      FuzzyDoubleComparator(double eps = 1e-9) : epsilon(eps) {}
      bool operator()(const double& a, const double& b) const 
      {
        return std::fabs(a - b) >= epsilon && a < b;
      }
    };
 
    using DeltaMassHistogram = std::map<double, double, FuzzyDoubleComparator>;
    using DeltaMassToChargeCount = std::map<double, int, FuzzyDoubleComparator>;
 
    OpenSearchModificationAnalysis() = default;
 
    ~OpenSearchModificationAnalysis() = default;
 
    std::pair<DeltaMassHistogram, DeltaMassToChargeCount> 
    analyzeDeltaMassPatterns(const PeptideIdentificationList& peptide_ids, 
                            bool use_smoothing = false, 
                            bool debug = false) const;
 
    std::vector<ModificationSummary>
    mapDeltaMassesToModifications(const DeltaMassHistogram& delta_mass_histogram,
                                 const DeltaMassToChargeCount& charge_histogram,
                                 PeptideIdentificationList& peptide_ids,
                                 double precursor_mass_tolerance = 5.0,
                                 bool precursor_mass_tolerance_unit_ppm = true,
                                 const std::string& output_file = "") const;
 
    std::vector<ModificationSummary>
    analyzeModifications(PeptideIdentificationList& peptide_ids,
                        double precursor_mass_tolerance = 5.0,
                        bool precursor_mass_tolerance_unit_ppm = true,
                        bool use_smoothing = false,
                        const std::string& output_file = "") const;
 
    OpenSearchAnalysisResult
    analyzeModificationsWithStatistics(PeptideIdentificationList& peptide_ids,
                                       double precursor_mass_tolerance = 5.0,
                                       bool precursor_mass_tolerance_unit_ppm = true,
                                       bool use_smoothing = false,
                                       const std::string& output_file = "") const;
 
    DeltaMassStatistics
    generateDeltaMassStatistics(const DeltaMassHistogram& histogram,
                                const DeltaMassToChargeCount& charge_histogram,
                                const PeptideIdentificationList& peptide_ids,
                                double precursor_mass_tolerance = 5.0,
                                bool precursor_mass_tolerance_unit_ppm = true) const;
 
    PTMStatistics
    generatePTMStatistics(const PeptideIdentificationList& peptide_ids,
                          double precursor_mass_tolerance = 5.0,
                          bool precursor_mass_tolerance_unit_ppm = true) const;
 
    std::map<char, int>
    analyzeResidueFrequency(const PeptideIdentificationList& peptide_ids,
                            double delta_mass,
                            double tolerance = 0.01) const;
 
    void writeDeltaMassStatistics(const DeltaMassStatistics& stats,
                                  const std::string& output_file) const;
 
    void writePTMStatistics(const PTMStatistics& stats,
                            const std::string& output_file) const;
 
    static constexpr double DELTA_MASS_ZERO_THRESHOLD_ = 0.05; // Da
 
    static bool isUnmodified(double delta_mass)
    {
      return std::abs(delta_mass) <= DELTA_MASS_ZERO_THRESHOLD_;
    }
 
  private:
 
    static constexpr double MAX_MOD_MAPPING_TOL_ = 0.02; // Da
 
    static constexpr double MOD_MAPPING_AMBIGUITY_TOL_ = 0.0005; // Da
 
    static constexpr double DELTA_MASS_BIN_WIDTH_ = 0.0005; // Da
 
    static Int64 deltaMassBin_(double delta_mass)
    {
      return static_cast<Int64>(std::llround(delta_mass / DELTA_MASS_BIN_WIDTH_));
    }
 
    static double gaussian_(double x, double sigma);
 
    static DeltaMassHistogram smoothDeltaMassHistogram_(const DeltaMassHistogram& histogram,
                                                       double sigma = 0.001);
 
    static DeltaMassHistogram findPeaksInHistogram_(const DeltaMassHistogram& histogram,
                                                   double count_threshold = 0.0,
                                                   double snr = 2.0);
 
    void writeModificationSummary_(const std::vector<ModificationSummary>& modifications,
                                  const std::string& output_file) const;
 
    std::map<double, std::string, FuzzyDoubleComparator>
    buildModificationMassLookup_() const;
 
    std::string getTargetResidues_(const std::string& mod_name) const;
 
    int countUniquePeptides_(const PeptideIdentificationList& peptide_ids,
                            double delta_mass,
                            double tolerance) const;
  };
 
} // namespace OpenMS