Created factory method for MbrScorer (#887)

* revised the RevisedModifiedPeptides function

* revised  the run error checking

* add the tester for merge

* Edit the tester for the new detectionType

* fix the bug and set intensity to 0 for ambiguous peak

* slove the ambiguous peak intensity problem

* Fix the bug on MBR score when there are all ambiguous peak inside

* change some comment

* refactored mbr scorer constructoin

* Edited mbrScorer test

* Added check to avoid invalid distributions in MbrScorer

* Broke some stuff

* Fixed remaining issues

* removed NaNs in ppm error calculations

* Changes to scorer changed the real data MBR test. Test assert statements have been amended

---------

Co-authored-by: RayMSMS <150720362+RayMSMS@users.noreply.github.com>

Author: Alexander-Sol

mzLib/FlashLFQ/FlashLFQResults.cs

@@ -810,6 +810,7 @@ internal void RevisedModifiedPeptides()
                     .DistinctBy(p=>p.ModifiedSequence)
                     .Select(p=>p.ModifiedSequence)
                     .ToList();
+
                 foreach (var modSeq in allIDs)
                 {
                     if (PeptideModifiedSequences.ContainsKey(modSeq))

mzLib/FlashLFQ/FlashLfqEngine.cs

@@ -796,55 +796,6 @@ internal RtInfo PredictRetentionTime(
             return new RtInfo(predictedRt: donorPeak.Apex.IndexedPeak.RetentionTime - medianRtDiff, width: rtRange);
         }
 
-        /// <summary>
-        /// Constructs a MbrScorer object that is used to score all MBR peaks for a given acceptor file
-        /// </summary>
-        /// <param name="acceptorFileIdentifiedPeaks"> All MSMS identified peaks in the acceptor file </param>
-        /// <param name="fileSpecificMbrTolerance">A ppm tolerance specific to the given file</param>
-        /// <returns> A MbrScorer object </returns>
-        private MbrScorer BuildMbrScorer(List<ChromatographicPeak> acceptorFileIdentifiedPeaks, out PpmTolerance fileSpecificMbrTolerance)
-        {
-            // Construct a distribution of ppm errors for all MSMS peaks in the acceptor file
-            var apexToAcceptorFilePeakDict = new Dictionary<IIndexedPeak, ChromatographicPeak>();
-            List<double> ppmErrors = new List<double>();
-            foreach (var peak in acceptorFileIdentifiedPeaks.Where(p => p.Apex != null
-                && PeptideModifiedSequencesToQuantify.Contains(p.Identifications.First().ModifiedSequence)
-                && p.Identifications.First().QValue < FlashParams.DonorQValueThreshold)) 
-            {
-                if (!apexToAcceptorFilePeakDict.ContainsKey(peak.Apex.IndexedPeak))
-                {
-                    apexToAcceptorFilePeakDict.Add(peak.Apex.IndexedPeak, peak);
-                }
-
-                ppmErrors.Add(peak.MassError);
-            }
-            if (ppmErrors.Count < 3)
-            {
-                fileSpecificMbrTolerance = null;
-                return null;
-            }
-            double ppmSpread = ppmErrors.Count > 30 ? ppmErrors.InterquartileRange() / 1.36 : ppmErrors.StandardDeviation();
-            Normal ppmDistribution = new Normal(ppmErrors.Median(), ppmSpread);
-            double fileSpecificMbrPpmTolerance = Math.Min(Math.Abs(ppmErrors.Median()) + ppmSpread * 4, FlashParams.MbrPpmTolerance);
-            fileSpecificMbrTolerance = new PpmTolerance(fileSpecificMbrPpmTolerance); // match between runs PPM tolerance
-
-            // Construct a distribution of peak log intensities for all MSMS peaks in the acceptor file
-            var acceptorFileLogIntensities = acceptorFileIdentifiedPeaks
-                .Where(p => p.Intensity > 0)
-                .Select(p => Math.Log(p.Intensity, 2))
-                .ToList();
-            double medianAcceptorLogIntensity = acceptorFileLogIntensities.Median();
-            Normal logIntensityDistribution = new Normal(acceptorFileLogIntensities.Median(), acceptorFileLogIntensities.InterquartileRange() / 1.36);
-            try // if the constructor fails, we don't want to crash the program
-            {
-                return new MbrScorer(apexToAcceptorFilePeakDict, acceptorFileIdentifiedPeaks, ppmDistribution, logIntensityDistribution);
-            }
-            catch
-            {
-                return null;
-            }
-        }
-
         /// <summary>
         /// Returns a pseudo-randomly selected peak that does not have the same mass as the donor
         /// </summary>
@@ -911,14 +862,16 @@ private void QuantifyMatchBetweenRunsPeaks(SpectraFileInfo acceptorFile)
                 .Count() > 1;
 
             // acceptor file known peaks
-            var acceptorFileIdentifiedPeaks = _results.Peaks[acceptorFile];
+            var acceptorFileIdentifiedPeaks = _results.Peaks[acceptorFile]
+                .Where(p => p.Identifications.Any(id => PeptideModifiedSequencesToQuantify.Contains(id.ModifiedSequence)))
+                .ToList();
 
             // these are the analytes already identified in this run. we don't need to try to match them from other runs
             var acceptorFileIdentifiedSequences = new HashSet<string>(acceptorFileIdentifiedPeaks
                 .Where(peak => peak.IsotopicEnvelopes.Any() && peak.Identifications.Min(id => id.QValue) < 0.01)
                 .SelectMany(p => p.Identifications.Select(d => d.ModifiedSequence)));
 
-            MbrScorer scorer = BuildMbrScorer(acceptorFileIdentifiedPeaks, out var mbrTol);
+            MbrScorer scorer = MbrScorerFactory.BuildMbrScorer(acceptorFileIdentifiedPeaks, FlashParams, out var mbrTol);
             if (scorer == null)
                 return;
 

mzLib/FlashLFQ/MBR/MbrScorer.cs

@@ -1,12 +1,13 @@
 using Easy.Common.EasyComparer;
+using MassSpectrometry;
 using MathNet.Numerics.Distributions;
 using MathNet.Numerics.Statistics;
+using MzLibUtil;
 using System;
 using System.Collections.Generic;
 using System.Data;
 using System.Data.Entity.ModelConfiguration.Conventions;
 using System.Linq;
-using MassSpectrometry;
 
 namespace FlashLFQ
 {
@@ -17,10 +18,11 @@ namespace FlashLFQ
     internal class MbrScorer
     {
         // Intensity and ppm distributions are specific to each acceptor file
-        private readonly Normal _logIntensityDistribution;
-        private readonly Normal _ppmDistribution;
-        private readonly Normal _scanCountDistribution;
-        private readonly Gamma _isotopicCorrelationDistribution;
+        private Normal _logIntensityDistribution;
+        private Normal _ppmDistribution;
+        private Normal _scanCountDistribution;
+        private  Gamma _isotopicCorrelationDistribution;
+
         // The logFcDistributions and rtDifference distributions are unique to each donor file - acceptor file pair
         private Dictionary<SpectraFileInfo, Normal> _logFcDistributionDictionary;
         private Dictionary<SpectraFileInfo, Normal> _rtPredictionErrorDistributionDictionary;
@@ -35,23 +37,69 @@ internal class MbrScorer
         /// </summary>
         internal MbrScorer(
             Dictionary<IIndexedPeak, ChromatographicPeak> apexToAcceptorFilePeakDict,
-            List<ChromatographicPeak> acceptorFileMsmsPeaks,
-            Normal ppmDistribution, 
-            Normal logIntensityDistribution)
+            List<ChromatographicPeak> unambiguousAcceptorFilePeaks)
         {
             ApexToAcceptorFilePeakDict = apexToAcceptorFilePeakDict;
-            UnambiguousMsMsAcceptorPeaks = acceptorFileMsmsPeaks.Where(p => p.Apex != null && p.DetectionType != DetectionType.MBR && p.NumIdentificationsByFullSeq == 1).ToList();
-            MaxNumberOfScansObserved = acceptorFileMsmsPeaks.Max(peak => peak.ScanCount);
-            _logIntensityDistribution = logIntensityDistribution;
-            _ppmDistribution = ppmDistribution;
-            _isotopicCorrelationDistribution = GetIsotopicEnvelopeCorrDistribution();
+            UnambiguousMsMsAcceptorPeaks = unambiguousAcceptorFilePeaks;
+            MaxNumberOfScansObserved = unambiguousAcceptorFilePeaks.Max(peak => peak.ScanCount);
+
+            // Initialize the dictionaries that will hold the log fold change and RT prediction error distributions
+            // for each donor file
             _logFcDistributionDictionary = new();
             _rtPredictionErrorDistributionDictionary = new();
+        }
+
+        /// <summary>
+        /// Constructs the distributions that are used to score MBR matches
+        /// </summary>
+        /// <returns>Returns true if the scorer was initialized successfully, false otherwise</returns>
+        internal bool InitializeScorer()
+        {
+            if (UnambiguousMsMsAcceptorPeaks.Count < 3)
+                return false;
+
+            // Populate distributions for scoring MBR matches
+            _logIntensityDistribution = GetLogIntensityDistribution();
+            _ppmDistribution = GetPpmErrorDistribution();
+            _isotopicCorrelationDistribution = GetIsotopicEnvelopeCorrDistribution();
+            _scanCountDistribution = GetScanCountDistribution();
+
+            return IsValid();
+        }
+
+        private Normal GetPpmErrorDistribution()
+        {
+            // Construct a distribution of ppm errors for all MSMS peaks in the acceptor file
+            List<double> ppmErrors = UnambiguousMsMsAcceptorPeaks.Select(p => p.MassError).Where(e => !double.IsNaN(e)).ToList();
+            if (ppmErrors.Count < 2)
+                return null;
+            double ppmSpread = ppmErrors.Count > 30 ? ppmErrors.InterquartileRange() / 1.36 : ppmErrors.StandardDeviation();
+            Normal ppmDistribution = new Normal(ppmErrors.Median(), ppmSpread);
+            return ppmDistribution;
+        }
+
+        private Normal GetLogIntensityDistribution()
+        {
+            var logIntensities = UnambiguousMsMsAcceptorPeaks
+                .Where(p => p.Intensity > 0)
+                .Select(p => Math.Log(p.Intensity, 2))
+                .ToList();
+
+            if (logIntensities.Count < 2)
+                return null;
+            
+            double mean = logIntensities.Median();
+            double stdDev = logIntensities.InterquartileRange() / 1.36;
+            return new Normal(mean, stdDev);
+        }
 
-            // This is kludgey, because scan counts are discrete
+        // This is kludgey, because scan counts are discrete
+        private Normal GetScanCountDistribution()
+        {
             List<double> scanList = UnambiguousMsMsAcceptorPeaks.Select(peak => (double)peak.ScanCount).ToList();
+
             // build a normal distribution for the scan list of the acceptor peaks
-            _scanCountDistribution = new Normal(scanList.Average(), scanList.Count > 30 ? scanList.StandardDeviation() : scanList.InterquartileRange() / 1.36);
+            return new Normal(scanList.Average(), scanList.Count > 30 ? scanList.StandardDeviation() : scanList.InterquartileRange() / 1.36);
         }
 
         /// <summary>
@@ -61,11 +109,13 @@ internal MbrScorer(
         private Gamma GetIsotopicEnvelopeCorrDistribution()
         {
             var pearsonCorrs = UnambiguousMsMsAcceptorPeaks.Select(p => 1 - p.IsotopicPearsonCorrelation).Where(p => p > 0).ToList();
-            if (pearsonCorrs.Count <= 1) return null;
+            if (pearsonCorrs.Count < 2) return null;
             double mean = pearsonCorrs.Mean();
             double variance = pearsonCorrs.Variance();
             var alpha = Math.Pow(mean, 2) / variance;
             var beta = mean / variance;
+            if (!Gamma.IsValidParameterSet(alpha, beta))
+                return null;
             return new Gamma(alpha, beta);
         }
 
@@ -98,20 +148,17 @@ internal void AddRtPredErrorDistribution(SpectraFileInfo donorFile, List<double>
                 rtPredictionErrors.Add(avgDiff - anchorPeptideRtDiffs[i]);
             }
 
-            Normal rtPredictionErrorDist = new Normal(0, 0); 
             // Default distribution. Effectively assigns a RT Score of zero if no alignment can be performed
             // between the donor and acceptor based on shared MS/MS IDs
-
-            if(rtPredictionErrors.Any())
+            Normal rtPredictionErrorDist = new Normal(0, 0);
+            if (rtPredictionErrors.Count >= 2)
             {
                 double medianRtError = rtPredictionErrors.Median();
                 double stdDevRtError = rtPredictionErrors.StandardDeviation();
-                if(stdDevRtError >= 0.0 && !double.IsNaN(medianRtError))
-                {
+                if (Normal.IsValidParameterSet(medianRtError, stdDevRtError))
                     rtPredictionErrorDist = new Normal(medianRtError, 1);
-                }
             }
-            
+
             _rtPredictionErrorDistributionDictionary.Add(donorFile, rtPredictionErrorDist);
         }
 
@@ -139,6 +186,14 @@ internal double ScoreMbr(MbrChromatographicPeak acceptorPeak, ChromatographicPea
                 * acceptorPeak.IsotopicDistributionScore, 0.20);
         }
 
+        /// <summary>
+        /// Returns the standard deviation of the Ppm error distribution + the median of the Ppm error distribution
+        /// </summary>
+        internal double GetPpmErrorTolerance()
+        {
+            return _ppmDistribution.StdDev * 4 + Math.Abs(_ppmDistribution.Median);
+        }
+
         // Setting a minimum score prevents the MBR score from going to zero if one component of that score is 0
         // 3e-7 is the fraction of a normal distribution that lies at least 5 stdDev away from the mean
         private double _minScore = 3e-7;
@@ -189,10 +244,6 @@ internal double CalculateIntensityScore(double acceptorIntensity, Chromatographi
         /// <param name="idDonorPeaks"> List of peaks in the donoro file. </param>
         internal void CalculateFoldChangeBetweenFiles(List<ChromatographicPeak> idDonorPeaks)
         {
-
-            var donorFileLogIntensities = idDonorPeaks.Where(p => p.Intensity > 0).Select(p => Math.Log(p.Intensity, 2)).ToList();
-            double medianDonorLogIntensity = donorFileLogIntensities.Median();
-
             // Find the difference in peptide intensities between donor and acceptor files
             // this intensity score creates a conservative bias in MBR
             List<double> listOfFoldChangesBetweenTheFiles = new List<double>();
@@ -244,7 +295,19 @@ internal bool IsValid(SpectraFileInfo donorFile)
         {
             return _rtPredictionErrorDistributionDictionary.TryGetValue(donorFile, out var rtDist)
                 && rtDist != null
-                && _ppmDistribution != null
+                && IsValid();
+        }
+
+        /// <summary>
+        /// This method checks whether the scorer is validly parameterized and capable of scoring MBR transfers
+        /// Notably, it is indifferent to the isotopic correlation distribution being null, as a null isotopic distribution correlation
+        /// results in all MBR transfers receiving the minimum score for the isotopic distribution component.
+        /// This could be changed in the future, but currently multiple tests results in null isotopic distributions, and will break if they can't do MBR
+        /// </summary>
+        /// <returns></returns>
+        internal bool IsValid()
+        {
+            return _ppmDistribution != null
                 && _scanCountDistribution != null
                 && _logIntensityDistribution != null;
         }

mzLib/FlashLFQ/MBR/MbrScorerFactory.cs

@@ -0,0 +1,45 @@
+using MassSpectrometry;
+using MathNet.Numerics.Distributions;
+using MathNet.Numerics.Statistics;
+using MzLibUtil;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace FlashLFQ
+{
+    internal static class MbrScorerFactory
+    {
+        /// <summary>
+        /// Constructs a MbrScorer object that is used to score all MBR peaks for a given acceptor file
+        /// </summary>
+        /// <param name="acceptorFileMsmsPeaks"> All MSMS identified peaks in the acceptor file that contain quantifiable peptides </param>
+        /// <param name="fileSpecificMbrTolerance">A ppm tolerance specific to the given file</param>
+        /// <returns> A MbrScorer object </returns>
+        public static MbrScorer BuildMbrScorer(List<ChromatographicPeak> acceptorFileMsmsPeaks,
+            FlashLfqParameters flashParams, out PpmTolerance fileSpecificMbrTolerance)
+        {
+            // Construct a dictionary linking each MSMS peaks to the indexed peak of its apex.
+            // This is to ensure MBR doesn't assign a peptide to a peak that is already claimed by one or more other peptides.
+            var apexToAcceptorFilePeakDict = acceptorFileMsmsPeaks
+                .Where(p => p.Apex != null)
+                .DistinctBy(p => p.Apex.IndexedPeak)
+                .ToDictionary(p => p.Apex.IndexedPeak, p => p);
+
+            MbrScorer scorer = new MbrScorer(apexToAcceptorFilePeakDict, acceptorFileMsmsPeaks);
+            // Try to initialize the scorer, which will fail if the acceptor file is empty or has no MSMS peaks.
+            if (!scorer.InitializeScorer())
+            {
+                fileSpecificMbrTolerance = null;
+                return null;
+            }
+
+            double mbrPpmTolerance = Math.Min(scorer.GetPpmErrorTolerance(), flashParams.MbrPpmTolerance);
+            fileSpecificMbrTolerance = new PpmTolerance(mbrPpmTolerance); // match between runs PPM tolerance
+            return scorer;
+        }
+
+    }
+}

mzLib/TestFlashLFQ/TestFlashLFQ.cs

@@ -1512,8 +1512,8 @@ public static void RealDataMbrTest()
             // Any change to ML.NET or the PEP Analysis engine will cause these to change.
             Console.WriteLine("r1 PIP event count: " + f1r1MbrResults.Count);
             Console.WriteLine("r2 PIP event count: " + f1r2MbrResults.Count);
-            Assert.AreEqual(141, f1r1MbrResults.Count);
-            Assert.AreEqual(77, f1r2MbrResults.Count);
+            Assert.AreEqual(140, f1r1MbrResults.Count);
+            Assert.AreEqual(78, f1r2MbrResults.Count);
 
             // Check that MS/MS identified peaks and MBR identified peaks have similar intensities 
             List<(double, double)> peptideIntensities = f1r1MbrResults.Select(pep => (Math.Log(pep.Value.GetIntensity(f1r1)), Math.Log(pep.Value.GetIntensity(f1r2)))).ToList();