Changes to MBR within FlashLFQ

mzLib/FlashLFQ/ChromatographicPeak.cs

@@ -1,9 +1,9 @@
-using MzLibUtil;
-using System;
+using System;
 using System.Collections.Generic;
 using System.Linq;
 using System.Text;
 using ClassExtensions = Chemistry.ClassExtensions;
+using FlashLFQ.PEP;
 
 namespace FlashLFQ
 {
@@ -16,20 +16,23 @@ public class ChromatographicPeak
         public int ScanCount => IsotopicEnvelopes.Count;
         public double SplitRT;
         public readonly bool IsMbrPeak;
-        public double PredictedRetentionTime { get; init; }
-
+        public double MbrScore;
+        public double PpmScore { get; set; }
+        public double IntensityScore { get; set; }
+        public double RtScore { get; set; }
+        public double ScanCountScore { get; set; }
+        public double IsotopicDistributionScore { get; set; }
         /// <summary>
-        /// A score bounded by 100 and 0, with more confident MBR-detections receiving higher scores
+        /// Stores the pearson correlation between the apex isotopic envelope and the theoretical isotopic distribution
         /// </summary>
-        public double MbrScore { get; private set; }
-
-        /// The four scores below are bounded by 0 and 1, with higher scores being better
-        public double PpmScore { get; private set; }
-        public double IntensityScore { get; private set; }
-        public double RtScore { get; private set; }
-        public double ScanCountScore { get; private set; }
-
-        public ChromatographicPeak(Identification id, bool isMbrPeak, SpectraFileInfo fileInfo)
+        public double IsotopicPearsonCorrelation => Apex?.PearsonCorrelation ?? -1;
+        public double RtPredictionError { get; set; }
+        public List<int> ChargeList { get; set; }
+        internal double MbrQValue { get; set; }
+        public ChromatographicPeakData PepPeakData { get; set; }
+        public double? MbrPep { get; set; }
+
+        public ChromatographicPeak(Identification id, bool isMbrPeak, SpectraFileInfo fileInfo, bool randomRt = false)
         {
             SplitRT = 0;
             NumChargeStatesObserved = 0;
@@ -40,12 +43,14 @@ public ChromatographicPeak(Identification id, bool isMbrPeak, SpectraFileInfo fi
             IsotopicEnvelopes = new List<IsotopicEnvelope>();
             IsMbrPeak = isMbrPeak;
             SpectraFileInfo = fileInfo;
+            RandomRt = randomRt;
         }
 
-        public ChromatographicPeak(Identification id, bool isMbrPeak, SpectraFileInfo fileInfo, double predictedRetentionTime) :
-            this(id, isMbrPeak, fileInfo)
+        public bool Equals(ChromatographicPeak peak)
         {
-            PredictedRetentionTime = predictedRetentionTime;
+            return SpectraFileInfo.Equals(peak.SpectraFileInfo) 
+                && Identifications.First().ModifiedSequence.Equals(peak.Identifications.First().ModifiedSequence)
+                && ApexRetentionTime == peak.ApexRetentionTime;
         }
 
         public IsotopicEnvelope Apex { get; private set; }
@@ -54,13 +59,18 @@ public ChromatographicPeak(Identification id, bool isMbrPeak, SpectraFileInfo fi
         public int NumIdentificationsByBaseSeq { get; private set; }
         public int NumIdentificationsByFullSeq { get; private set; }
         public double MassError { get; private set; }
+        /// <summary>
+        /// Bool that describes whether the retention time of this peak was randomized
+        /// If true, implies that this peak is a decoy peak identified by the MBR algorithm
+        /// </summary>
+        public bool RandomRt { get; }
+        public bool DecoyPeptide => Identifications.First().IsDecoy;
 
         public void CalculateIntensityForThisFeature(bool integrate)
         {
             if (IsotopicEnvelopes.Any())
             {
-                double maxIntensity = IsotopicEnvelopes.Max(p => p.Intensity);
-                Apex = IsotopicEnvelopes.First(p => p.Intensity == maxIntensity);
+                Apex = IsotopicEnvelopes.MaxBy(p => p.Intensity);
 
                 if (integrate)
                 {
@@ -123,25 +133,6 @@ public void ResolveIdentifications()
             this.NumIdentificationsByBaseSeq = Identifications.Select(v => v.BaseSequence).Distinct().Count();
             this.NumIdentificationsByFullSeq = Identifications.Select(v => v.ModifiedSequence).Distinct().Count();
         }
-
-        /// <summary>
-        /// Calculates four component scores and one overarching Mbr score for an MBR peak.
-        /// MBR Score is equal to 100 * the geometric mean of the four component scores.
-        /// </summary>
-        /// <param name="scorer"> An MbrScorer specific to the file where this peak was found </param>
-        /// <param name="donorPeak"> The donor peak used as the basis for the MBR identification. </param>
-        internal void CalculateMbrScore(MbrScorer scorer, ChromatographicPeak donorPeak)
-        {
-            if (SpectraFileInfo != scorer.AcceptorFile) throw new MzLibException("Error when performing match-between-runs: Mismatch between scorer and peak.");
-
-            IntensityScore = scorer.CalculateIntensityScore(this, donorPeak);
-            RtScore = scorer.CalculateRetentionTimeScore(this, donorPeak);
-            PpmScore = scorer.CalculatePpmErrorScore(this);
-            ScanCountScore = scorer.CalculateScanCountScore(this);
-
-            MbrScore = 100 * Math.Pow(IntensityScore * RtScore * PpmScore * ScanCountScore, 0.25);
-        }
-
         public static string TabSeparatedHeader
         {
             get
@@ -164,20 +155,18 @@ public static string TabSeparatedHeader
                 sb.Append("Peak Charge" + "\t");
                 sb.Append("Num Charge States Observed" + "\t");
                 sb.Append("Peak Detection Type" + "\t");
-                sb.Append("MBR Score" + "\t");
-                sb.Append("Ppm Score" + "\t");
-                sb.Append("Intensity Score" + "\t");
-                sb.Append("Rt Score" + "\t");
-                sb.Append("Scan Count Score" + "\t");
+                sb.Append("PIP Q-Value" + "\t");
+                sb.Append("PIP PEP" + "\t");
                 sb.Append("PSMs Mapped" + "\t");
                 sb.Append("Base Sequences Mapped" + "\t");
                 sb.Append("Full Sequences Mapped" + "\t");
                 sb.Append("Peak Split Valley RT" + "\t");
-                sb.Append("Peak Apex Mass Error (ppm)");
+                sb.Append("Peak Apex Mass Error (ppm)" + "\t");
+                sb.Append("Decoy Peptide" + "\t");
+                sb.Append("Random RT");
                 return sb.ToString();
             }
         }
-
         public override string ToString()
         {
             StringBuilder sb = new StringBuilder();
@@ -260,17 +249,16 @@ public override string ToString()
                 sb.Append("" + "MSMS" + "\t");
             }
 
-            sb.Append("" + (IsMbrPeak ? MbrScore.ToString() : "") + "\t");
-            sb.Append("" + (IsMbrPeak ? PpmScore.ToString() : "") + "\t");
-            sb.Append("" + (IsMbrPeak ? IntensityScore.ToString() : "") + "\t");
-            sb.Append("" + (IsMbrPeak ? RtScore.ToString() : "") + "\t");
-            sb.Append("" + (IsMbrPeak ? ScanCountScore.ToString() : "") + "\t");
+            sb.Append("" + (IsMbrPeak ? MbrQValue.ToString() : "") + "\t");
+            sb.Append("" + (IsMbrPeak ? MbrPep.ToString() : "") + "\t");
 
             sb.Append("" + Identifications.Count + "\t");
             sb.Append("" + NumIdentificationsByBaseSeq + "\t");
             sb.Append("" + NumIdentificationsByFullSeq + "\t");
             sb.Append("" + SplitRT + "\t");
             sb.Append("" + MassError);
+            sb.Append("\t" + DecoyPeptide);
+            sb.Append("\t" + RandomRt);
 
             return sb.ToString();
         }

mzLib/FlashLFQ/FlashLFQ.csproj

@@ -13,6 +13,8 @@
   <ItemGroup>
     <PackageReference Include="CsvHelper" Version="32.0.3" />
     <PackageReference Include="MathNet.Numerics" Version="5.0.0" />
+    <PackageReference Include="Microsoft.ML" Version="3.0.1" />
+    <PackageReference Include="Microsoft.ML.FastTree" Version="3.0.1" />
     <PackageReference Include="NetSerializer" Version="4.1.2" />
     <PackageReference Include="SharpLearning.Optimization" Version="[0.28.0]" />
   </ItemGroup>

mzLib/FlashLFQ/FlashLFQResults.cs

@@ -15,27 +15,26 @@ public class FlashLfqResults
         public readonly Dictionary<string, ProteinGroup> ProteinGroups;
         public readonly Dictionary<SpectraFileInfo, List<ChromatographicPeak>> Peaks;
         private readonly HashSet<string> _peptideModifiedSequencesToQuantify;
+        public string PepResultString { get; set; }
+        public double MbrQValueThreshold { get; set; }
 
-        public FlashLfqResults(List<SpectraFileInfo> spectraFiles, List<Identification> identifications, HashSet<string> peptides = null)
+        public FlashLfqResults(List<SpectraFileInfo> spectraFiles, List<Identification> identifications, double mbrQValueThreshold = 0.05,
+            HashSet<string> peptideModifiedSequencesToQuantify = null)
         {
             SpectraFiles = spectraFiles;
             PeptideModifiedSequences = new Dictionary<string, Peptide>();
             ProteinGroups = new Dictionary<string, ProteinGroup>();
             Peaks = new Dictionary<SpectraFileInfo, List<ChromatographicPeak>>();
-            if(peptides == null || !peptides.Any())
-            {
-                peptides = identifications.Select(id => id.ModifiedSequence).ToHashSet();
-            }
-            _peptideModifiedSequencesToQuantify = peptides;
+            MbrQValueThreshold = mbrQValueThreshold;
+            _peptideModifiedSequencesToQuantify = peptideModifiedSequencesToQuantify ?? identifications.Where(id => !id.IsDecoy).Select(id => id.ModifiedSequence).ToHashSet();
 
             foreach (SpectraFileInfo file in spectraFiles)
             {
                 Peaks.Add(file, new List<ChromatographicPeak>());
             }
 
-
             // Only quantify peptides within the set of valid peptide modified (full) sequences. This is done to enable pepitde-level FDR control of reported results
-            foreach (Identification id in identifications.Where(id => peptides.Contains(id.ModifiedSequence)))
+            foreach (Identification id in identifications.Where(id => !id.IsDecoy & _peptideModifiedSequencesToQuantify.Contains(id.ModifiedSequence)))
             {
                 if (!PeptideModifiedSequences.TryGetValue(id.ModifiedSequence, out Peptide peptide))
                 {
@@ -59,6 +58,17 @@ public FlashLfqResults(List<SpectraFileInfo> spectraFiles, List<Identification>
             }
         }
 
+        public void ReNormalizeResults(bool integrate = false, int maxThreads = 10, bool useSharedPeptides = false)
+        {
+            foreach(var peak in Peaks.SelectMany(p => p.Value))
+            {
+                peak.CalculateIntensityForThisFeature(integrate);
+            }
+            new IntensityNormalizationEngine(this, integrate, silent: true, maxThreads).NormalizeResults();
+            CalculatePeptideResults(quantifyAmbiguousPeptides: false);
+            CalculateProteinResultsMedianPolish(useSharedPeptides: useSharedPeptides);
+        }
+
         public void MergeResultsWith(FlashLfqResults mergeFrom)
         {
             this.SpectraFiles.AddRange(mergeFrom.SpectraFiles);
@@ -128,6 +138,8 @@ public void CalculatePeptideResults(bool quantifyAmbiguousPeptides)
             {
                 var groupedPeaks = filePeaks.Value
                     .Where(p => p.NumIdentificationsByFullSeq == 1)
+                    .Where(p => !p.Identifications.First().IsDecoy)
+                    .Where(p => !p.IsMbrPeak || (p.MbrQValue < MbrQValueThreshold && !p.RandomRt))
                     .GroupBy(p => p.Identifications.First().ModifiedSequence)
                     .Where(group => _peptideModifiedSequencesToQuantify.Contains(group.Key))
                     .ToList();
@@ -163,11 +175,15 @@ public void CalculatePeptideResults(bool quantifyAmbiguousPeptides)
                 // report ambiguous quantification
                 var ambiguousPeaks = filePeaks.Value
                     .Where(p => p.NumIdentificationsByFullSeq > 1)
+                    .Where(p => !p.Identifications.First().IsDecoy)
+                    .Where(p => !p.IsMbrPeak || (p.MbrQValue < MbrQValueThreshold && !p.RandomRt))
                     .ToList();
                 foreach (ChromatographicPeak ambiguousPeak in ambiguousPeaks)
                 {
-                    foreach (Identification id in ambiguousPeak.Identifications)
+                    foreach (Identification id in ambiguousPeak.Identifications.Where(id => !id.IsDecoy))
                     {
+                        if (!_peptideModifiedSequencesToQuantify.Contains(id.ModifiedSequence)) continue; // Ignore the ids/sequences we don't want to quantify
+
                         string sequence = id.ModifiedSequence;
 
                         double alreadyRecordedIntensity = PeptideModifiedSequences[sequence].GetIntensity(filePeaks.Key);
@@ -224,7 +240,7 @@ private void HandleAmbiguityInFractions()
 
                     foreach (SpectraFileInfo file in sample)
                     {
-                        foreach (ChromatographicPeak peak in Peaks[file])
+                        foreach (ChromatographicPeak peak in Peaks[file].Where(p => !p.IsMbrPeak || p.MbrQValue < MbrQValueThreshold))
                         {
                             foreach (Identification id in peak.Identifications)
                             {
@@ -549,7 +565,7 @@ public void CalculateProteinResultsMedianPolish(bool useSharedPeptides)
             }
         }
 
-        public void WriteResults(string peaksOutputPath, string modPeptideOutputPath, string proteinOutputPath, string bayesianProteinQuantOutput, bool silent = true)
+        public void WriteResults(string peaksOutputPath, string modPeptideOutputPath, string proteinOutputPath, string bayesianProteinQuantOutput, bool silent)
         {
             if (!silent)
             {