Merge branch 'develop' of https://github.com/OpenMS/OpenMS into toppa…

…s_outdir

CHANGELOG

@@ -26,11 +26,11 @@ Library:
 - better documentation for all SpectraFilter... tools (#7183)
 - TOPPView: offer Ion mobility view from 2D spectra view (#7423)
 - TOPPView: view ion mobility frames, irrespective of its MS level (formerly only MS1 was supported) (#7427)
+- OpenSwath: Add output on peak shape metrics to .osw file (#7222)
 
 New Tools:
 - IonMobilityBinning - Merges spectra with similar IM values and creates @p N output mzML's by discretizing the IM range (#7459)
 
-
 Fixes:
 - FileConverter: more robust (#7176)
 - MSFragger: allow relative path to database (#7155)

src/openms/include/OpenMS/ANALYSIS/OPENSWATH/DIAScoring.h

@@ -30,19 +30,22 @@ namespace OpenMS
     compute a number of scores based on the full MS2 spectrum available. The scores are the following:
 
     - isotope scores:
-      -- isotope_corr: computes the correlation of each fragment ion with the
+      - isotope_corr: computes the correlation of each fragment ion with the
          theoretical isotope distribution. This is the pearson correlation to
          the theoretical isotope pattern weighted by the relative intensity of
          the transition (more is better).
-      -- isotope_overlap: checks whether a signal at position (mz - 1) / charge
+      - isotope_overlap: checks whether a signal at position (mz - 1) / charge
          exists and how strong it is. This would be an indication that the current
          peak is an isotopic signal of another peak. This simply counts how
          often a peak was observed that is higher than the current peak, thus
          number is then weighted by the relative intensity of the transition
          (thus less is better here).
 
-    - massdiff score: computes the difference in ppm of the experimental signal
-         to the expected signal (thus less is better)
+    - massdiff score: computes the difference in ppm of the experimental signal to the expected signal (thus less is better). 
+      - Equation: sum(ppm_difference) / # transitions
+      - Notes: 
+        - Divide by the total number of transitions and is thus quite punishing if a transition is missing
+        - Also outputs a list of all the ppm differences, if signal is not found output -1.0
 
     - b/y ion score: checks for the presence of b/y ions of the peptide in question
 

src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathOSWWriter.h

@@ -92,17 +92,13 @@ namespace OpenMS
     String input_filename_;
     OpenMS::UInt64 run_id_;
     bool doWrite_;
-    bool use_ms1_traces_;
-    bool sonar_;
     bool enable_uis_scoring_;
 
   public:
 
     OpenSwathOSWWriter(const String& output_filename,
                        const UInt64 run_id,
                        const String& input_filename = "inputfile",
-                       bool ms1_scores = false,
-                       bool sonar = false,
                        bool uis_scores = false);
 
     bool isActive() const;

src/openms/include/OpenMS/ANALYSIS/OPENSWATH/OpenSwathScores.h

@@ -43,6 +43,7 @@ namespace OpenMS
     bool use_uis_scores = true;
     bool use_ionseries_scores = true;
     bool use_ms2_isotope_scores = true;
+    bool use_peak_shape_metrics = false;
   };
 
   /** @brief A structure to hold the different scores computed by OpenSWATH
@@ -182,11 +183,28 @@ namespace OpenMS
     std::vector<double> ind_total_area_intensity;
     std::vector<double> ind_intensity_score;
     std::vector<double> ind_apex_intensity;
+    std::vector<double> ind_apex_position;
+    std::vector<double> ind_fwhm;
     std::vector<double> ind_total_mi;
     std::vector<double> ind_log_intensity;
     std::vector<double> ind_intensity_ratio;
     std::vector<double> ind_mi_ratio;
     std::vector<double> ind_mi_score;
+
+    // peak shape metrics
+    std::vector<double> ind_start_position_at_5;
+    std::vector<double> ind_end_position_at_5;
+    std::vector<double> ind_start_position_at_10;
+    std::vector<double> ind_end_position_at_10;
+    std::vector<double> ind_start_position_at_50;
+    std::vector<double> ind_end_position_at_50;
+    std::vector<double> ind_total_width;
+    std::vector<double> ind_tailing_factor;
+    std::vector<double> ind_asymmetry_factor;
+    std::vector<double> ind_slope_of_baseline;
+    std::vector<double> ind_baseline_delta_2_height;
+    std::vector<double> ind_points_across_baseline;
+    std::vector<double> ind_points_across_half_height;
 
     OpenSwath_Ind_Scores() = default;
 

src/openms/include/OpenMS/ANALYSIS/OPENSWATH/PeakIntegrator.h

@@ -772,11 +772,22 @@ namespace OpenMS
       psm.width_at_50 = psm.end_position_at_50 - psm.start_position_at_50;
       psm.total_width = (p.PosEnd(right) - 1)->getPos() - p.PosBegin(left)->getPos();
       psm.slope_of_baseline = (p.PosEnd(right) - 1)->getIntensity() - p.PosBegin(left)->getIntensity();
-      psm.baseline_delta_2_height = psm.slope_of_baseline / peak_height;
+      if (peak_height != 0.0) // avoid division by zero
+      {
+        psm.baseline_delta_2_height = psm.slope_of_baseline / peak_height;
+      }
       // Source of tailing_factor and asymmetry_factor formulas:
       // USP 40 - NF 35 The United States Pharmacopeia and National Formulary - Supplementary
-      psm.tailing_factor = psm.width_at_5 / (2*(peak_apex_pos - psm.start_position_at_5));
-      psm.asymmetry_factor = (psm.end_position_at_10 - peak_apex_pos) / (peak_apex_pos - psm.start_position_at_10);
+
+      // Can only compute if start and peak apex are different 
+      if (psm.start_position_at_5 != peak_apex_pos)
+      {
+        psm.tailing_factor = psm.width_at_5 / (2*(peak_apex_pos - psm.start_position_at_5));
+      }
+      if (psm.start_position_at_10 != peak_apex_pos)
+      {
+        psm.asymmetry_factor = (psm.end_position_at_10 - peak_apex_pos) / (peak_apex_pos - psm.start_position_at_10);
+      }
       return psm;
     }
 

src/openms/source/ANALYSIS/OPENSWATH/DIAScoring.cpp

@@ -131,11 +131,11 @@ namespace OpenMS
       // about the mass difference if no signal is present.
       if (!signalFound)
       {
+        diff_ppm.push_back(-1); // if no signal is found than we set the ppm to -1
         continue;
       }
 
       double ppm = Math::getPPM(mz, transition.getProductMZ());
-      diff_ppm.push_back(transition.getProductMZ());
       diff_ppm.push_back(ppm);
       ppm_score += std::fabs(ppm);
       ppm_score_weighted += std::fabs(ppm) * normalized_library_intensity[k];

src/openms/source/ANALYSIS/OPENSWATH/MRMFeatureFinderScoring.cpp

@@ -130,6 +130,7 @@ namespace OpenMS
     scores_to_use.setValidStrings("use_ms1_mi", {"true","false"});
     scores_to_use.setValue("use_uis_scores", "false", "Use UIS scores for peptidoform identification", {"advanced"});
     scores_to_use.setValidStrings("use_uis_scores", {"true","false"});
+    scores_to_use.setValue("use_peak_shape_metrics", "false", "Use peak shape metrics for scoring", {"advanced"});
     scores_to_use.setValue("use_ionseries_scores", "true", "Use MS2-level b/y ion-series scores for peptidoform identification", {"advanced"});
     scores_to_use.setValidStrings("use_ionseries_scores", {"true","false"});
     scores_to_use.setValue("use_ms2_isotope_scores", "true", "Use MS2-level isotope scores (pearson & manhattan) across product transitions (based on ID if annotated or averagine)", {"advanced"});
@@ -344,16 +345,6 @@ namespace OpenMS
                                               signal_noise_estimators_identification,
                                               idscores);
 
-      std::vector<String> ind_transition_names;
-      std::vector<double> ind_area_intensity;
-      std::vector<double> ind_total_area_intensity;
-      std::vector<double> ind_intensity_score;
-      std::vector<double> ind_apex_intensity;
-      std::vector<double> ind_total_mi;
-      std::vector<double> ind_log_intensity;
-      std::vector<double> ind_intensity_ratio;
-      std::vector<double> ind_mi_ratio;
-
       std::vector<double> ind_mi_score;
       if (su_.use_mi_score_)
       {
@@ -362,7 +353,7 @@ namespace OpenMS
 
       for (size_t i = 0; i < native_ids_identification.size(); i++)
       {
-        ind_transition_names.emplace_back(native_ids_identification[i]);
+        idscores.ind_transition_names.emplace_back(native_ids_identification[i]);
         if (idmrmfeature.getFeature(native_ids_identification[i]).getIntensity() > 0)
         {
           double intensity_score = double(idmrmfeature.getFeature(native_ids_identification[i]).getIntensity()) / double(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("total_xic"));
@@ -384,37 +375,67 @@ namespace OpenMS
             if (mi_ratio > 1) { mi_ratio = 1 / mi_ratio; }
           }
 
-          ind_area_intensity.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getIntensity());
-          ind_total_area_intensity.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("total_xic"));
-          ind_intensity_score.push_back(intensity_score);
-          ind_apex_intensity.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("peak_apex_int"));
-          ind_total_mi .push_back(total_mi);
-          ind_log_intensity.push_back(std::log(idmrmfeature.getFeature(native_ids_identification[i]).getIntensity()));
-          ind_intensity_ratio.push_back(intensity_ratio);
-          ind_mi_ratio.push_back(mi_ratio);
+          idscores.ind_area_intensity.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getIntensity());
+          idscores.ind_total_area_intensity.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("total_xic"));
+          idscores.ind_intensity_score.push_back(intensity_score);
+          idscores.ind_apex_intensity.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("peak_apex_int"));
+          idscores.ind_apex_position.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("peak_apex_position"));
+          idscores.ind_fwhm.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("width_at_50"));
+          idscores.ind_total_mi .push_back(total_mi);
+          idscores.ind_log_intensity.push_back(std::log(idmrmfeature.getFeature(native_ids_identification[i]).getIntensity()));
+          idscores.ind_intensity_ratio.push_back(intensity_ratio);
+          idscores.ind_mi_ratio.push_back(mi_ratio);
+
+          if (su_.use_peak_shape_metrics)
+          {
+            idscores.ind_start_position_at_5.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("start_position_at_5"));
+            idscores.ind_end_position_at_5.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("end_position_at_5"));
+            idscores.ind_start_position_at_10.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("start_position_at_10"));
+            idscores.ind_end_position_at_10.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("end_position_at_10"));
+            idscores.ind_start_position_at_50.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("start_position_at_50"));
+            idscores.ind_end_position_at_50.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("end_position_at_50"));
+            idscores.ind_total_width.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("total_width"));
+            idscores.ind_tailing_factor.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("tailing_factor"));
+            idscores.ind_asymmetry_factor.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("asymmetry_factor"));
+            idscores.ind_slope_of_baseline.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("slope_of_baseline"));
+            idscores.ind_baseline_delta_2_height.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("baseline_delta_2_height"));
+            idscores.ind_points_across_baseline.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("points_across_baseline"));
+            idscores.ind_points_across_half_height.push_back(idmrmfeature.getFeature(native_ids_identification[i]).getMetaValue("points_across_half_height"));
+          }
         }
         else
         {
-          ind_area_intensity.push_back(0);
-          ind_total_area_intensity.push_back(0);
-          ind_intensity_score.push_back(0);
-          ind_apex_intensity.push_back(0);
-          ind_total_mi.push_back(0);
-          ind_log_intensity.push_back(0);
-          ind_intensity_ratio.push_back(0);
-          ind_mi_ratio.push_back(0);
+          idscores.ind_area_intensity.push_back(0);
+          idscores.ind_total_area_intensity.push_back(0);
+          idscores.ind_intensity_score.push_back(0);
+          idscores.ind_apex_intensity.push_back(0);
+          idscores.ind_apex_position.push_back(0);
+          idscores.ind_fwhm.push_back(0);
+          idscores.ind_total_mi.push_back(0);
+          idscores.ind_log_intensity.push_back(0);
+          idscores.ind_intensity_ratio.push_back(0);
+          idscores.ind_mi_ratio.push_back(0);
+
+          if (su_.use_peak_shape_metrics)
+          {
+            idscores.ind_start_position_at_5.push_back(0);
+            idscores.ind_end_position_at_5.push_back(0);
+            idscores.ind_start_position_at_10.push_back(0);
+            idscores.ind_end_position_at_10.push_back(0);
+            idscores.ind_start_position_at_50.push_back(0);
+            idscores.ind_end_position_at_50.push_back(0);
+            idscores.ind_total_width.push_back(0);
+            idscores.ind_tailing_factor.push_back(0);
+            idscores.ind_asymmetry_factor.push_back(0);
+            idscores.ind_slope_of_baseline.push_back(0);
+            idscores.ind_baseline_delta_2_height.push_back(0);
+            idscores.ind_points_across_baseline.push_back(0);
+            idscores.ind_points_across_half_height.push_back(0);
+          }
         }
-      }
-      idscores.ind_transition_names = ind_transition_names;
-      idscores.ind_area_intensity = ind_area_intensity;
-      idscores.ind_total_area_intensity = ind_total_area_intensity;
-      idscores.ind_intensity_score = ind_intensity_score;
-      idscores.ind_apex_intensity = ind_apex_intensity;
-      idscores.ind_total_mi = ind_total_mi;
-      idscores.ind_log_intensity = ind_log_intensity;
-      idscores.ind_intensity_ratio = ind_intensity_ratio;
-      idscores.ind_mi_ratio = ind_mi_ratio;
       idscores.ind_num_transitions = native_ids_identification.size();
+
+      }
     }
 
     // Compute DIA scores only on the identification transitions
@@ -1062,6 +1083,7 @@ namespace OpenMS
     su_.use_uis_scores           = param_.getValue("Scores:use_uis_scores").toBool();
     su_.use_ionseries_scores     = param_.getValue("Scores:use_ionseries_scores").toBool();
     su_.use_ms2_isotope_scores   = param_.getValue("Scores:use_ms2_isotope_scores").toBool();
+    su_.use_peak_shape_metrics   = param_.getValue("Scores:use_peak_shape_metrics").toBool();
   }
 
   void MRMFeatureFinderScoring::mapExperimentToTransitionList(const OpenSwath::SpectrumAccessPtr& input,