Q value threshold adjustment and bug fix

MetaMorpheus/EngineLayer/FdrAnalysis/FdrAnalysisEngine.cs

@@ -1,9 +1,7 @@
 using EngineLayer.CrosslinkSearch;
 using System;
 using System.Collections.Generic;
-using System.IO;
 using System.Linq;
-using System.Text.RegularExpressions;
 
 namespace EngineLayer.FdrAnalysis
 {
@@ -14,7 +12,17 @@ public class FdrAnalysisEngine : MetaMorpheusEngine
         private readonly string AnalysisType;
         private readonly string OutputFolder; // used for storing PEP training models  
         private readonly bool DoPEP;
-
+        private static bool qvalueThresholdOverride;
+        private readonly int PsmCountThresholdForInvertedQvalue = 1000;
+        /// <summary>
+        /// This is to be used only for unit testing. Threshold for q-value calculation is set to 1000
+        /// However, many unit tests don't generate that many PSMs. Therefore, this property is used to override the threshold
+        /// to enable PEP calculation in unit tests with lower number of PSMs
+        /// </summary>
+        public static bool QvalueThresholdOverride   // property
+        {
+            set { qvalueThresholdOverride = value; }  // set method
+        }
         public FdrAnalysisEngine(List<SpectralMatch> psms, int massDiffAcceptorNumNotches, CommonParameters commonParameters,
             List<(string fileName, CommonParameters fileSpecificParameters)> fileSpecificParameters, List<string> nestedIds, string analysisType = "PSM", 
             bool doPEP = true, string outputFolder = null) : base(commonParameters, fileSpecificParameters, nestedIds)
@@ -71,10 +79,10 @@ private void DoFalseDiscoveryRateAnalysis(FdrAnalysisResults myAnalysisResults)
                         .Select(b => b.FirstOrDefault())
                         .ToList();
 
-                if (psms.Count > 100 & DoPEP)
+                if ((psms.Count > PsmCountThresholdForInvertedQvalue || qvalueThresholdOverride) & DoPEP)
                 {
                     CalculateQValue(psms, peptideLevelCalculation: false, pepCalculation: false);
-                    if (peptides.Count > 100 )
+                    if (peptides.Count > PsmCountThresholdForInvertedQvalue || qvalueThresholdOverride)
                     {
                         CalculateQValue(peptides, peptideLevelCalculation: true, pepCalculation: false);
 
@@ -160,6 +168,7 @@ public void CalculateQValue(List<SpectralMatch> psms, bool peptideLevelCalculati
                 // Stop if canceled
                 if (GlobalVariables.StopLoops) { break; }
 
+                // we have to keep track of q-values separately for each notch
                 int notch = psm.Notch ?? MassDiffAcceptorNumNotches;
                 if (psm.IsDecoy)
                 {
@@ -198,7 +207,9 @@ public void CalculateQValue(List<SpectralMatch> psms, bool peptideLevelCalculati
             }
             else
             {
-                if(psms.Count < 100)
+                //the QValueThreshodOverride condition here can be problematic in unit tests. If tests fail unexpectedly,
+                //try setting FdrAnalysisEngine.QvalueThresholdOverride = false; in the test method
+                if (psms.Count < PsmCountThresholdForInvertedQvalue && !qvalueThresholdOverride)
                 {
 
                    QValueTraditional(psms, peptideLevelAnalysis: peptideLevelCalculation);
@@ -216,39 +227,49 @@ public void CalculateQValue(List<SpectralMatch> psms, bool peptideLevelCalculati
         private void QValueTraditional(List<SpectralMatch> psms, bool peptideLevelAnalysis)
         {
             double qValue = 0;
-            double qValueNotch = 0;
+            double[] qValueNotch = new double[MassDiffAcceptorNumNotches + 1];
+
             for (int i = 0; i < psms.Count; i++)
             {
                 // Stop if canceled
                 if (GlobalVariables.StopLoops) { break; }
-
+                int notch = psms[i].Notch ?? MassDiffAcceptorNumNotches;
                 qValue = Math.Max(qValue, psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoy / Math.Max(psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeTarget, 1));
-                qValueNotch = Math.Max(qValueNotch, psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoyNotch / Math.Max(psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeTargetNotch, 1));
+                qValueNotch[notch] = Math.Max(qValueNotch[notch], psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoyNotch / Math.Max(psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeTargetNotch, 1));
 
                 psms[i].GetFdrInfo(peptideLevelAnalysis).QValue = Math.Min(qValue, 1);
-                psms[i].GetFdrInfo(peptideLevelAnalysis).QValueNotch = Math.Min(qValueNotch, 1);
+                psms[i].GetFdrInfo(peptideLevelAnalysis).QValueNotch = Math.Min(qValueNotch[notch], 1);
             }
         }
 
-        private static void QValueInverted(List<SpectralMatch> psms, bool peptideLevelAnalysis)
+        private void QValueInverted(List<SpectralMatch> psms, bool peptideLevelAnalysis)
         {
+            double[] qValueNotch = new double[MassDiffAcceptorNumNotches + 1];
+            bool[] qValueNotchCalculated = new bool[MassDiffAcceptorNumNotches + 1];
             psms.Reverse();
             //this calculation is performed from bottom up. So, we begin the loop by computing qValue
             //and qValueNotch for the last/lowest scoring psm in the bunch
             double qValue = (psms[0].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoy + 1) / psms[0].GetFdrInfo(peptideLevelAnalysis).CumulativeTarget;
-            double qValueNotch = (psms[0].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoyNotch + 1) / psms[0].GetFdrInfo(peptideLevelAnalysis).CumulativeTargetNotch;
 
             //Assign FDR values to PSMs
             for (int i = 0; i < psms.Count; i++)
             {
                 // Stop if canceled
                 if (GlobalVariables.StopLoops) { break; }
+                int notch = psms[i].Notch ?? MassDiffAcceptorNumNotches;
+
+                // populate the highest q-Value for each notch 
+                if (!qValueNotchCalculated[notch])
+                {
+                    qValueNotch[notch] = (psms[0].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoyNotch + 1) / psms[0].GetFdrInfo(peptideLevelAnalysis).CumulativeTargetNotch;
+                    qValueNotchCalculated[notch] = true;
+                }
 
                 qValue = Math.Min(qValue, (psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoy + 1) / Math.Max(psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeTarget, 1));
-                qValueNotch = Math.Min(qValueNotch, (psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoyNotch + 1) / Math.Max(psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeTargetNotch, 1));
+                qValueNotch[notch] = Math.Min(qValueNotch[notch], (psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeDecoyNotch + 1) / Math.Max(psms[i].GetFdrInfo(peptideLevelAnalysis).CumulativeTargetNotch, 1));
 
                 psms[i].GetFdrInfo(peptideLevelAnalysis).QValue = Math.Min(qValue, 1);
-                psms[i].GetFdrInfo(peptideLevelAnalysis).QValueNotch = Math.Min(qValueNotch, 1);
+                psms[i].GetFdrInfo(peptideLevelAnalysis).QValueNotch = Math.Min(qValueNotch[notch], 1);
             }
             psms.Reverse(); //we inverted the psms for this calculation. now we need to put them back into the original order
         }

MetaMorpheus/TaskLayer/CalibrationTask/CalibrationTask.cs

@@ -82,56 +82,86 @@
                 string originalUncalibratedFilePath = currentRawFileList[spectraFileIndex];
                 string originalUncalibratedFilenameWithoutExtension = Path.GetFileNameWithoutExtension(originalUncalibratedFilePath);
                 string calibratedFilePath = Path.Combine(OutputFolder, originalUncalibratedFilenameWithoutExtension + CalibSuffix + ".mzML");
-
+                bool calibrated = false;
                 // mark the file as in-progress
                 StartingDataFile(originalUncalibratedFilePath, new List<string> { taskId, "Individual Spectra Files", originalUncalibratedFilePath });
-
                 CommonParameters combinedParams = SetAllFileSpecificCommonParams(CommonParameters, fileSettingsList[spectraFileIndex]);
 
                 // load the file
                 Status("Loading spectra file...", new List<string> { taskId, "Individual Spectra Files" });
-
                 MsDataFile myMsDataFile = myFileManager.LoadFile(originalUncalibratedFilePath, combinedParams);
 
                 // get datapoints to fit calibration function to
                 Status("Acquiring calibration data points...", new List<string> { taskId, "Individual Spectra Files" });
                 DataPointAquisitionResults acquisitionResults = null;
-
                 acquisitionResults = GetDataAcquisitionResults(myMsDataFile, originalUncalibratedFilePath, variableModifications, fixedModifications, proteinList, taskId, combinedParams, new PpmTolerance(MaxPrecursorTolerance), new PpmTolerance(MaxProductTolerance));
-                if(acquisitionResults.Psms.Count >= NumRequiredPsms && acquisitionResults.Ms1List.Count >= NumRequiredMs1Datapoints && acquisitionResults.Ms2List.Count >= NumRequiredMs2Datapoints)
-                {
+                int numPsms = acquisitionResults.Psms.Count;
 
+                // check if we have enough data points to calibrate and then calibrate
+                if (acquisitionResults.Psms.Count >= NumRequiredPsms && acquisitionResults.Ms1List.Count >= NumRequiredMs1Datapoints && acquisitionResults.Ms2List.Count >= NumRequiredMs2Datapoints)
+                {
                     // generate calibration function and shift data points
                     Status("Calibrating...", new List<string> { taskId, "Individual Spectra Files" });
                     CalibrationEngine engine = new(myMsDataFile, acquisitionResults, combinedParams, FileSpecificParameters, new List<string> { taskId, "Individual Spectra Files", originalUncalibratedFilenameWithoutExtension });
                     _ = engine.Run();
 
-                    //update file
-                    myMsDataFile = engine.CalibratedDataFile;
-
-                    // write the calibrated mzML file
-                    myMsDataFile.ExportAsMzML(calibratedFilePath, CalibrationParameters.WriteIndexedMzml);
-                    myFileManager.DoneWithFile(originalUncalibratedFilePath);
+                    // get the calibrated data points again to see if there was an increase
+                    acquisitionResults = GetDataAcquisitionResults(engine.CalibratedDataFile, originalUncalibratedFilePath, variableModifications, fixedModifications, proteinList, taskId, combinedParams, new PpmTolerance(MaxPrecursorTolerance), new PpmTolerance(MaxProductTolerance));
 
                     // write toml settings for the calibrated file
                     string newTomlFileName = Path.Combine(OutputFolder, originalUncalibratedFilenameWithoutExtension + CalibSuffix + ".toml");
-
                     FileSpecificParameters fileSpecificParams = new();
 
                     // carry over file-specific parameters from the uncalibrated file to the calibrated one
                     if (fileSettingsList[spectraFileIndex] != null)
                     {
                         fileSpecificParams = fileSettingsList[spectraFileIndex].Clone();
                     }
-
-                    //suggest 4 * interquartile range as the ppm tolerance
                     fileSpecificParams.PrecursorMassTolerance = new PpmTolerance((4.0 * acquisitionResults.PsmPrecursorIqrPpmError) + Math.Abs(acquisitionResults.PsmPrecursorMedianPpmError));
                     fileSpecificParams.ProductMassTolerance = new PpmTolerance((4.0 * acquisitionResults.PsmProductIqrPpmError) + Math.Abs(acquisitionResults.PsmProductMedianPpmError));
 
-                    Toml.WriteFile(fileSpecificParams, newTomlFileName, tomlConfig);
 
-                    FinishedWritingFile(newTomlFileName, new List<string> { taskId, "Individual Spectra Files", originalUncalibratedFilenameWithoutExtension });
+                    // if we have more PSMs after calibration, we should re-calibrate
+                    if (acquisitionResults.Psms.Count > numPsms)
+                    {
+                        calibrated = true;
+                        numPsms = acquisitionResults.Psms.Count;
+                        //update file
+                        myMsDataFile = engine.CalibratedDataFile;
 
+                        //suggest 4 * interquartile range as the ppm tolerance
+                        fileSpecificParams.PrecursorMassTolerance = new PpmTolerance((4.0 * acquisitionResults.PsmPrecursorIqrPpmError) + Math.Abs(acquisitionResults.PsmPrecursorMedianPpmError));
+                        fileSpecificParams.ProductMassTolerance = new PpmTolerance((4.0 * acquisitionResults.PsmProductIqrPpmError) + Math.Abs(acquisitionResults.PsmProductMedianPpmError));
+
+                        // generate calibration function and shift data points
+                        Status("Calibrating...", new List<string> { taskId, "Individual Spectra Files" });
+                        engine = new(myMsDataFile, acquisitionResults, combinedParams, FileSpecificParameters, new List<string> { taskId, "Individual Spectra Files", originalUncalibratedFilenameWithoutExtension });
+                        _ = engine.Run();
+
+                        // get the calibrated data points again to see if there was an increase
+                        acquisitionResults = GetDataAcquisitionResults(engine.CalibratedDataFile, originalUncalibratedFilePath, variableModifications, fixedModifications, proteinList, taskId, combinedParams, new PpmTolerance((4.0 * acquisitionResults.PsmPrecursorIqrPpmError) + Math.Abs(acquisitionResults.PsmPrecursorMedianPpmError)), new PpmTolerance((4.0 * acquisitionResults.PsmProductIqrPpmError) + Math.Abs(acquisitionResults.PsmProductMedianPpmError)));
+
+                        if (acquisitionResults.Psms.Count > numPsms)
+                        {
+                            //update file
+                            myMsDataFile = engine.CalibratedDataFile;
+                            //suggest 4 * interquartile range as the ppm tolerance
+                            fileSpecificParams.PrecursorMassTolerance = new PpmTolerance((4.0 * acquisitionResults.PsmPrecursorIqrPpmError) + Math.Abs(acquisitionResults.PsmPrecursorMedianPpmError));
+                            fileSpecificParams.ProductMassTolerance = new PpmTolerance((4.0 * acquisitionResults.PsmProductIqrPpmError) + Math.Abs(acquisitionResults.PsmProductMedianPpmError));
+
+                        }
+
+                        // write the calibrated mzML file
+                        myMsDataFile.ExportAsMzML(calibratedFilePath, CalibrationParameters.WriteIndexedMzml);
+
+                    }
+                    myFileManager.DoneWithFile(originalUncalibratedFilePath);
+                    if(!calibrated)
+                    {
+                        myMsDataFile.ExportAsMzML(calibratedFilePath, CalibrationParameters.WriteIndexedMzml);
+                    }
+                    Toml.WriteFile(fileSpecificParams, newTomlFileName, tomlConfig);
+                    FinishedWritingFile(newTomlFileName, new List<string> { taskId, "Individual Spectra Files", originalUncalibratedFilenameWithoutExtension });
                     // finished calibrating this file
                     FinishedWritingFile(calibratedFilePath, new List<string> { taskId, "Individual Spectra Files", originalUncalibratedFilenameWithoutExtension });
                     MyTaskResults.NewSpectra.Add(calibratedFilePath);

MetaMorpheus/Test/MatchIonsOfAllCharges.cs

@@ -20,6 +20,7 @@
 using Omics.Modifications;
 using Omics.SpectrumMatch;
 using static System.Net.WebRequestMethods;
+using EngineLayer.FdrAnalysis;
 
 namespace Test
 {
@@ -408,6 +409,7 @@ public static void TestLibraryGeneration()
 
         public static void TestLibraryUpdate()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             string thisTaskOutputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\UpdateLibrary");
             _ = Directory.CreateDirectory(thisTaskOutputFolder);
             SearchTask task = Toml.ReadFile<SearchTask>(Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\SpectralSearchTask.toml"), MetaMorpheusTask.tomlConfig);
@@ -425,7 +427,6 @@ public static void TestLibraryUpdate()
 
             string rawCopy = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\UpdateLibrary\rawCopy.mzML");
             System.IO.File.Copy(raw1, rawCopy);
-
             EverythingRunnerEngine UpdateLibrary = new(new List<(string, MetaMorpheusTask)> { ("UpdateSpectraFileOutput", task) }, new List<string> { raw1, raw2 }, new List<DbForTask> { new DbForTask(lib, false), new DbForTask( db1,false), new DbForTask(db2, false) }, thisTaskOutputFolder);
 
             UpdateLibrary.Run();