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);
 
@@ -161,6 +169,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)
                 {
@@ -199,7 +208,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);
@@ -217,39 +228,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/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();

MetaMorpheus/Test/MetaDraw/MetaDrawTest.cs

@@ -20,6 +20,7 @@
 using Proteomics.ProteolyticDigestion;
 using Readers;
 using TaskLayer;
+using EngineLayer.FdrAnalysis;
 
 namespace Test.MetaDraw
 {
@@ -229,6 +230,7 @@ public static void MetaDraw_TestStationarySequencePositioning()
         [Test]
         public static void MetaDraw_SearchTaskTest()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false; 
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"MetaDraw_SearchTaskTest");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\smalldb.fasta");
             string spectraFile = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SmallCalibratible_Yeast.mzML");
@@ -693,6 +695,7 @@ public static void TestMetaDrawXlSpectralLibrary()
         [Test]
         public static void MetaDraw_GlycoSearchTaskWithChildScansTest()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"MetaDraw_GlycoSearchTaskTest");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"GlycoTestData\leukosialin.fasta");
             string spectraFile = Path.Combine(TestContext.CurrentContext.TestDirectory, @"GlycoTestData\sliced_glyco_hcd_ethcd.raw");
@@ -854,6 +857,7 @@ public static void MetaDraw_GlycoSearchTaskWithChildScansTest()
         [Test]
         public static void MetaDraw_TestChimeraScanSpectrumDisplay()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"MetaDraw_SearchTaskTest");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\hela_snip_for_unitTest.fasta");
             string spectraFile = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\TaGe_SA_HeLa_04_subset_longestSeq.mzML");
@@ -1081,9 +1085,9 @@ public static void TestMetaDrawErrors()
         }
 
         [Test]
-        [NonParallelizable]
         public static void TestMetaDrawLoadingWithWeirdFileNames()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             // test loading when the file has a periods, commas, spaces in the name
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestMetaDrawLoadingWithWeirdFileNames");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\smalldb.fasta");
@@ -1133,6 +1137,7 @@ public static void TestMetaDrawLoadingWithWeirdFileNames()
         [Test]
         public static void TestMetaDrawWithSpectralLibrary()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestMetaDrawWithSpectraLibrary");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\P16858.fasta");
             string library1 = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\P16858_target.msp");
@@ -1210,6 +1215,7 @@ public static void TestMetaDrawWithSpectralLibrary()
         [Test]
         public static void TestPsmFromTsvIonParsing()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestPsmFromTsvIonParsing");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\P16858.fasta");
             string spectraFile = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\slicedMouse.raw");
@@ -1318,6 +1324,7 @@ public static void SequenceCoverageMapTest()
         [Test]
         public static void TestMetaDrawLogicCleanUp()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestMetaDrawWithSpectraLibrary");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\P16858.fasta");
             string library1 = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SpectralLibrarySearch\P16858_target.msp");
@@ -1400,6 +1407,7 @@ public static void TestMetaDrawLogicCleanUp()
         [Test]
         public static void TestMetaDrawOutputFormats()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"MetaDraw_SearchTaskTest");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\smalldb.fasta");
             string spectraFile = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SmallCalibratible_Yeast.mzML");
@@ -1509,6 +1517,7 @@ public static void EnsureNoCrashesWithNGlyco()
         [Test]
         public static void TestMetaDrawSequenceDisplayOutputs()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"MetaDraw_SearchTaskTest");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\smalldb.fasta");
             string spectraFile = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SmallCalibratible_Yeast.mzML");

MetaMorpheus/Test/MetaDraw/SpectrumMatchPlotTests.cs

@@ -14,6 +14,7 @@
 using Omics.Fragmentation;
 using Readers;
 using TaskLayer;
+using EngineLayer.FdrAnalysis;
 
 namespace Test.MetaDraw
 {
@@ -31,6 +32,7 @@ public class SpectrumMatchPlotTests
         [OneTimeSetUp]
         public void Setup()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"MetaDraw_PeakAnnotaitonTest");
             string proteinDatabase = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\smalldb.fasta");
             string spectraFile = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestData\SmallCalibratible_Yeast.mzML");

MetaMorpheus/Test/PostSearchAnalysisTaskTests.cs

@@ -2,7 +2,9 @@
 using System.Diagnostics.CodeAnalysis;
 using System.IO;
 using System.Linq;
-using NUnit.Framework; using Assert = NUnit.Framework.Legacy.ClassicAssert;
+using EngineLayer.FdrAnalysis;
+using NUnit.Framework; 
+using Assert = NUnit.Framework.Legacy.ClassicAssert;
 
 namespace Test
 {
@@ -18,6 +20,7 @@ public static class PostSearchAnalysisTaskTests
         [Test]
         public static void AllResultsAndResultsTxtContainsCorrectValues_QValue_BottomUp()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = true;
             //First test that AllResults and Results display correct numbers of peptides and psms with q-value filter on
             EverythingRunnerEngineTestCase.TryGetTestCase(EverythingRunnerEngineTestCases.BottomUpQValue, out var testCase);
             string outputFolder = testCase.OutputDirectory;
@@ -68,11 +71,13 @@ public static void AllResultsAndResultsTxtContainsCorrectValues_QValue_BottomUp(
             Assert.AreEqual("All target PSMs with q-value <= 0.01: " + TaGe_SA_A549_3_snip_2ExpectedPsms, singleFileResults[5]);
             Assert.AreEqual("All target peptides with q-value <= 0.01: " + TaGe_SA_A549_3_snip_2ExpectedPeptides, singleFileResults[6]);
             Assert.AreEqual("All target protein groups with q-value <= 0.01 (1% FDR): 165", singleFileResults[7]);
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
         }
 
         [Test]
         public static void AllResultsAndResultsTxtContainsCorrectValues_PepQValue_BottomUp()
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = true;
             //First test that AllResults and Results display correct numbers of peptides and psms with pep q-value filter on
             EverythingRunnerEngineTestCase.TryGetTestCase(EverythingRunnerEngineTestCases.BottomUpPepQValue, out var testCase);
             string outputFolder = testCase.OutputDirectory;
@@ -99,6 +104,7 @@ public static void AllResultsAndResultsTxtContainsCorrectValues_PepQValue_Bottom
             Assert.AreEqual("TaGe_SA_A549_3_snip_2 - Target PSMs with pep q-value <= 0.01: 190", results[13]);
             Assert.AreEqual("TaGe_SA_A549_3_snip_2 - Target peptides with pep q-value <= 0.01: 153", results[14]);
             Assert.AreEqual("TaGe_SA_A549_3_snip_2 - Target protein groups within 1 % FDR: 140", results[15]);
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
         }
 
         /// <summary>
@@ -108,6 +114,7 @@ public static void AllResultsAndResultsTxtContainsCorrectValues_PepQValue_Bottom
         [TestCaseSource(nameof(GetTestCases))]
         public static void AllResultTxtContainsCorrectNumberOfResultLines(EverythingRunnerEngineTestCases testCaseIdentifier)
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false; 
             var testCase = EverythingRunnerEngineTestCase.GetTestCase(testCaseIdentifier);
 
             int expectedIndividualFileLines = testCase.DataFileList.Count == 1 || !testCase.WriteIndividualResults 
@@ -153,6 +160,7 @@ public static void AllResultTxtContainsCorrectNumberOfResultLines(EverythingRunn
         [TestCaseSource(nameof(GetTestCases))]
         public static void CorrectFilesAreWrittenWithCorrectName(EverythingRunnerEngineTestCases testCaseIdentifier)
         {
+            FdrAnalysisEngine.QvalueThresholdOverride = false; 
             var testCase = EverythingRunnerEngineTestCase.GetTestCase(testCaseIdentifier);
             var psmFiles = Directory.GetFiles(testCase.OutputDirectory, "*PSMs.psmtsv", SearchOption.AllDirectories);
             var pepXmlFiles = Directory.GetFiles(testCase.OutputDirectory, "*.pep.xml", SearchOption.AllDirectories);

MetaMorpheus/Test/ProteinGroupTest.cs

@@ -14,6 +14,7 @@
 using Omics.Digestion;
 using Omics.Modifications;
 using UsefulProteomicsDatabases;
+using EngineLayer.FdrAnalysis;
 
 namespace Test
 {
@@ -206,6 +207,8 @@ public static void TestModificationInfoListInProteinGroupsOutput()
             string mzmlName = @"TestData\PrunedDbSpectra.mzml";
             string fastaName = @"TestData\DbForPrunedDb.fasta";
             string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, @"TestPrunedGeneration");
+
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             var engine = new EverythingRunnerEngine(taskList, new List<string> { mzmlName }, new List<DbForTask> { new DbForTask(fastaName, false) }, outputFolder);
             engine.Run();
             string final = Path.Combine(MySetUpClass.outputFolder, "task2", "DbForPrunedDbGPTMDproteinPruned.xml");

MetaMorpheus/Test/QuantificationTest.cs

@@ -15,6 +15,7 @@
 using System.Text;
 using Omics.Modifications;
 using TaskLayer;
+using EngineLayer.FdrAnalysis;
 
 namespace Test
 {
@@ -224,6 +225,7 @@ public static void TestProteinQuantFileHeaders(bool hasDefinedExperimentalDesign
             }
 
             // run the search/quantification
+            FdrAnalysisEngine.QvalueThresholdOverride = false;
             SearchTask task = new SearchTask();
             task.RunTask(unitTestFolder, new List<DbForTask> { new DbForTask(dbName, false) }, fileInfos.Select(p => p.FullFilePathWithExtension).ToList(), "");