Invert q-value

MetaMorpheus/EngineLayer/FdrAnalysis/FdrAnalysisEngine.cs

@@ -1,6 +1,8 @@
 using System;
 using System.Collections.Generic;
 using System.Linq;
+using EngineLayer;
+using EngineLayer.FdrAnalysis;
 
 namespace EngineLayer.FdrAnalysis
 {
@@ -10,13 +12,13 @@ public class FdrAnalysisEngine : MetaMorpheusEngine
         private readonly int MassDiffAcceptorNumNotches;
         private readonly double ScoreCutoff;
         private readonly string AnalysisType;
-        private readonly string OutputFolder; // used for storing PEP training models
+        private readonly string OutputFolder; // used for storing PEP training models  
         private readonly bool DoPEP;
 
         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)
         {
-            AllPsms = psms.OrderByDescending(p=>p).ToList();
+            AllPsms = psms.OrderByDescending(p => p).ToList();
             MassDiffAcceptorNumNotches = massDiffAcceptorNumNotches;
             ScoreCutoff = commonParameters.ScoreCutoff;
             AnalysisType = analysisType;
@@ -52,50 +54,14 @@ private void DoFalseDiscoveryRateAnalysis(FdrAnalysisResults myAnalysisResults)
                 QValueTraditional(psms);
                 if (psms.Count > 100)
                 {
-                    QValueInverted(psms);
-                }
-
-                // set q-value thresholds such that a lower scoring PSM can't have
-                // a higher confidence than a higher scoring PSM
-                //Populate min qValues
-                double qValueThreshold = 1.0;
-                double[] qValueNotchThreshold = new double[MassDiffAcceptorNumNotches + 1];
-                for (int i = 0; i < qValueNotchThreshold.Length; i++)
-                {
-                    qValueNotchThreshold[i] = 1.0;
-                }
-
-                for (int i = psms.Count - 1; i >= 0; i--)
-                {
-                    SpectralMatch psm = psms[i];
-
-                    // threshold q-values
-                    if (psm.FdrInfo.QValue > qValueThreshold)
+                    if (DoPEP)
                     {
-                        psm.FdrInfo.QValue = qValueThreshold;
-                    }
-                    else if (psm.FdrInfo.QValue < qValueThreshold)
-                    {
-                        qValueThreshold = psm.FdrInfo.QValue;
-                    }
-
-                    // threshold notch q-values
-                    int notch = psm.Notch ?? MassDiffAcceptorNumNotches;
-                    if (psm.FdrInfo.QValueNotch > qValueNotchThreshold[notch])
-                    {
-                        psm.FdrInfo.QValueNotch = qValueNotchThreshold[notch];
-                    }
-                    else if (psm.FdrInfo.QValueNotch < qValueNotchThreshold[notch])
-                    {
-                        qValueNotchThreshold[notch] = psm.FdrInfo.QValueNotch;
+                        Compute_PEPValue(myAnalysisResults);
                     }
+                    QValueInverted(psms);
                 }
+                CountPsm(psms);
             }
-            if (DoPEP)
-            {
-                Compute_PEPValue(myAnalysisResults);
-            }
-            CountPsm();
         }
 
         private static void QValueInverted(List<SpectralMatch> psms)
@@ -194,7 +160,6 @@ public void Compute_PEPValue(FdrAnalysisResults myAnalysisResults)
 
                     Compute_PEPValue_Based_QValue(AllPsms);
                 }
-                CountPsm(); // recounting Psm's after PEP based disambiguation
             }
 
             if (AnalysisType == "Peptide")
@@ -223,35 +188,31 @@ public static void Compute_PEPValue_Based_QValue(List<SpectralMatch> psms)
                 psms[psmsArrayIndicies[i]].FdrInfo.PEP_QValue = Math.Round(qValue, 6);
             }
         }
-
-        public void CountPsm()
+        /// <summary>
+        /// This method gets the PSM count for each psm with q-value < 0.01 and adds that information to the individual psm.
+        /// </summary>
+        public void CountPsm(List<SpectralMatch> proteasePsms)
         {
-            var psmsGroupedByProtease = AllPsms.GroupBy(p => p.DigestionParams.Protease);
+            // exclude ambiguous psms and has a fdr cutoff = 0.01
+            var allUnambiguousPsms = proteasePsms.Where(psm => psm.FullSequence != null).ToList();
 
-            foreach (var proteasePsms in psmsGroupedByProtease)
-            {
-                // exclude ambiguous psms and has a fdr cutoff = 0.01
-                var allUnambiguousProteasePsms = proteasePsms.Where(p => p.FullSequence != null).ToList();
+            var unambiguousPsmsLessThanOnePercentFdr = allUnambiguousPsms.Where(psm =>
+                    psm.FdrInfo.QValue <= 0.01
+                    && psm.FdrInfo.QValueNotch <= 0.01)
+                .GroupBy(p => p.FullSequence);
 
-                var fullSequenceGroups = allUnambiguousProteasePsms.Where(p => p.FdrInfo.QValue < 0.01 && p.FdrInfo.QValueNotch < 0.01)
-                    .Select(p => p.FullSequence).GroupBy(s => s);
+            Dictionary<string, int> sequenceToPsmCount = new Dictionary<string, int>();
 
-                Dictionary<string, int> sequenceToPsmCount = new Dictionary<string, int>();
-                foreach (var fullSequence in fullSequenceGroups)
-                {
-                    sequenceToPsmCount.Add(fullSequence.Key, fullSequence.Count());
-                }
+            foreach (var sequenceGroup in unambiguousPsmsLessThanOnePercentFdr)
+            {
+                sequenceToPsmCount.TryAdd(sequenceGroup.First().FullSequence, sequenceGroup.Count());
+            }
 
-                foreach (SpectralMatch psm in allUnambiguousProteasePsms)
+            foreach (SpectralMatch psm in allUnambiguousPsms)
+            {
+                if (sequenceToPsmCount.TryGetValue(psm.FullSequence, out int count))
                 {
-                    if (sequenceToPsmCount.TryGetValue(psm.FullSequence, out int count))
-                    {
-                        psm.PsmCount = count;
-                    }
-                    else
-                    {
-                        psm.PsmCount = 0;
-                    }
+                    psm.PsmCount = count;
                 }
             }
         }

MetaMorpheus/EngineLayer/PsmTsv/PsmTsvWriter.cs

@@ -10,6 +10,7 @@
 using Omics;
 using Omics.Modifications;
 
+
 namespace EngineLayer
 {
     public static class PsmTsvWriter
@@ -161,6 +162,7 @@ internal static (string ResolvedString, string ResolvedValue) Resolve(IEnumerabl
             }
         }
 
+        //The null-coalescing operator ?? returns the value of its left-hand operand if it isn't null; otherwise, it evaluates the right-hand operand and returns its result.
         internal static void AddBasicMatchData(Dictionary<string, string> s, SpectralMatch psm)
         {
             s[PsmTsvHeader.FileName] = psm == null ? " " : Path.GetFileNameWithoutExtension(psm.FullFilePath);
@@ -177,6 +179,7 @@ internal static void AddBasicMatchData(Dictionary<string, string> s, SpectralMat
             s[PsmTsvHeader.Notch] = psm == null ? " " : Resolve(psm.BestMatchingBioPolymersWithSetMods.Select(p => p.Notch)).ResolvedString;
         }
 
+        //The null-coalescing operator ?? returns the value of its left-hand operand if it isn't null; otherwise, it evaluates the right-hand operand and returns its result.
         internal static void AddPeptideSequenceData(Dictionary<string, string> s, SpectralMatch sm, IReadOnlyDictionary<string, int> ModsToWritePruned)
         {
             bool pepWithModsIsNull = sm == null || sm.BestMatchingBioPolymersWithSetMods == null || !sm.BestMatchingBioPolymersWithSetMods.Any();
@@ -216,7 +219,7 @@ internal static void AddPeptideSequenceData(Dictionary<string, string> s, Spectr
                         .Where(d => Includes(b, d))
                         .Select(d => $"{d.OneBasedBeginPosition.ToString()}-{d.OneBasedEndPosition.ToString()}")))).ResolvedString;
             }
-           
+
             s[PsmTsvHeader.Contaminant] = pepWithModsIsNull ? " " : Resolve(pepsWithMods.Select(b => b.Parent.IsContaminant ? "Y" : "N")).ResolvedString;
             s[PsmTsvHeader.Decoy] = pepWithModsIsNull ? " " : Resolve(pepsWithMods.Select(b => b.Parent.IsDecoy ? "Y" : "N")).ResolvedString;
             s[PsmTsvHeader.PeptideDesicription] = pepWithModsIsNull ? " " : Resolve(pepsWithMods.Select(b => b.Description)).ResolvedString;

MetaMorpheus/Test/MultiProteaseParsimonyTest.cs

@@ -990,7 +990,7 @@ public static void MultiProteaseParsimony_TestingGreedyAlgorithm()
         }
 
         /// <summary>
-        /// This test ensures that FDR for each psm is calculated accoriding to its protease
+        /// This test ensures that FDR for each psm is calculated according to its protease
         /// </summary>
         [Test]
         public static void MultiProteaseParsimony_TestingProteaseSpecificFDRCalculations()
@@ -1032,16 +1032,15 @@ public static void MultiProteaseParsimony_TestingProteaseSpecificFDRCalculations
             new FdrAnalysisEngine(psms, 0, new CommonParameters(), fsp, new List<string>()).Run();
             psms = psms.OrderByDescending(p => p.Score).ToList();
 
-            Assert.AreEqual(0.00, Math.Round(psms[0].FdrInfo.QValue, 2));
             Assert.AreEqual(0.00, Math.Round(psms[1].FdrInfo.QValue, 2));
             Assert.AreEqual(0.00, Math.Round(psms[2].FdrInfo.QValue, 2));
             Assert.AreEqual(0.00, Math.Round(psms[3].FdrInfo.QValue, 2));
-            Assert.AreEqual(0.33, Math.Round(psms[4].FdrInfo.QValue, 2));
+            Assert.AreEqual(0.5, Math.Round(psms[4].FdrInfo.QValue, 2));
             Assert.AreEqual(0.33, Math.Round(psms[5].FdrInfo.QValue, 2));
             Assert.AreEqual(0.00, Math.Round(psms[6].FdrInfo.QValue, 2));
             Assert.AreEqual(0.33, Math.Round(psms[7].FdrInfo.QValue, 2));
-            Assert.AreEqual(0.50, Math.Round(psms[8].FdrInfo.QValue, 2));
-            Assert.AreEqual(0.50, Math.Round(psms[9].FdrInfo.QValue, 2));
+            Assert.AreEqual(0.67, Math.Round(psms[8].FdrInfo.QValue, 2));
+            Assert.AreEqual(0.5, Math.Round(psms[9].FdrInfo.QValue, 2));
         }
     }
 }

MetaMorpheus/Test/PeptideSpectralMatchTest.cs

@@ -0,0 +1,102 @@
+using EngineLayer;
+using MassSpectrometry;
+using NUnit.Framework;
+using Proteomics.ProteolyticDigestion;
+using Proteomics;
+using System;
+using System.Collections.Generic;
+using Omics.Digestion;
+using Omics.Fragmentation;
+using Omics.Modifications;
+
+namespace Test
+{
+    [TestFixture]
+    internal class PeptideSpectralMatchTest
+    {
+        [Test]
+        public static void GetAminoAcidCoverageTest()
+        {
+            CommonParameters commonParams = new CommonParameters(digestionParams: new DigestionParams(protease: "trypsin", minPeptideLength: 1));
+
+            MsDataScan scanNumberOne = new MsDataScan(new MzSpectrum(new double[] { 10 }, new double[] { 1 }, false), 1, 2, true, Polarity.Positive, double.NaN, null, null, MZAnalyzerType.Orbitrap, double.NaN, null, null, "scan=1", 10, 2, 100, double.NaN, null, DissociationType.AnyActivationType, 0, null);
+
+            Ms2ScanWithSpecificMass ms2ScanOneMzTen = new Ms2ScanWithSpecificMass(scanNumberOne, 10, 2, "File", new CommonParameters());
+
+            string sequence = "";
+            Dictionary<string, Modification> allKnownMods = new();
+            int numFixedMods = 0;
+            DigestionParams digestionParams = new DigestionParams();
+            Protein myProtein = new Protein(sequence, "ACCESSION");
+            int oneBasedStartResidueInProtein = 0;
+            int oneBasedEndResidueInProtein = Math.Max(sequence.Length, 0);
+            int missedCleavages = 0;
+            CleavageSpecificity cleavageSpecificity = CleavageSpecificity.Full;
+            string peptideDescription = null;
+            int? pairedTargetDecoyHash = null;
+
+            PeptideWithSetModifications pwsmNoBaseSequence = new(sequence, allKnownMods, numFixedMods, digestionParams, myProtein,
+                oneBasedStartResidueInProtein, oneBasedEndResidueInProtein, missedCleavages, cleavageSpecificity,
+                peptideDescription, pairedTargetDecoyHash);
+            PeptideSpectralMatch psmNoBaseSequenceNoMFI = new(pwsmNoBaseSequence, 0, 10, 0, ms2ScanOneMzTen, commonParams,
+                new List<MatchedFragmentIon>());
+            psmNoBaseSequenceNoMFI.ResolveAllAmbiguities();
+
+            //PSM has neither sequence nor matched fragment ions
+            var b = psmNoBaseSequenceNoMFI.BaseSequence;
+            Assert.AreEqual("", b);
+            var m = psmNoBaseSequenceNoMFI.MatchedFragmentIons;
+            Assert.AreEqual(0, m.Count);
+            psmNoBaseSequenceNoMFI.GetAminoAcidCoverage();
+
+            sequence = "PEPTIDE";
+            oneBasedEndResidueInProtein = Math.Max(sequence.Length, 0);
+            myProtein = new Protein(sequence, "ACCESSION");
+            PeptideWithSetModifications pwsmBaseSequence = new(sequence, allKnownMods, numFixedMods, digestionParams, myProtein,
+                oneBasedStartResidueInProtein, oneBasedEndResidueInProtein, missedCleavages, cleavageSpecificity,
+                peptideDescription, pairedTargetDecoyHash);
+            PeptideSpectralMatch psmBaseSequenceNoMFI = new(pwsmBaseSequence, 0, 10, 0, ms2ScanOneMzTen, commonParams,
+                new List<MatchedFragmentIon>());
+
+            //PSM has sequence but does not have matched fragment ions
+            psmBaseSequenceNoMFI.ResolveAllAmbiguities();
+            b = psmBaseSequenceNoMFI.BaseSequence;
+            Assert.AreEqual(sequence, b);
+            m = psmBaseSequenceNoMFI.MatchedFragmentIons;
+            Assert.AreEqual(0, m.Count);
+            psmBaseSequenceNoMFI.GetAminoAcidCoverage();
+
+            //PSM has no sequence but does have matched fragment ions
+            Product ntProduct = new Product(ProductType.y, FragmentationTerminus.N, 1, 1, 1, 0);
+            MatchedFragmentIon mfi = new(ntProduct, 1, 1, 1);
+            List<MatchedFragmentIon> mfiList = new List<MatchedFragmentIon>() { mfi };
+            PeptideSpectralMatch psmNoBaseSequenceMFI = new(pwsmNoBaseSequence, 0, 10, 0, ms2ScanOneMzTen, commonParams,
+                mfiList);
+            psmNoBaseSequenceMFI.ResolveAllAmbiguities();
+
+            b = psmNoBaseSequenceMFI.BaseSequence;
+            Assert.AreEqual("", b);
+            m = psmNoBaseSequenceMFI.MatchedFragmentIons;
+            Assert.AreEqual(1, m.Count);
+            psmNoBaseSequenceMFI.GetAminoAcidCoverage();
+
+            //PSM has sequence and matched fragment ions
+            PeptideSpectralMatch psmBaseSequenceMFI = new(pwsmBaseSequence, 0, 10, 0, ms2ScanOneMzTen, commonParams,
+                mfiList);
+            psmBaseSequenceMFI.ResolveAllAmbiguities();
+
+            b = psmBaseSequenceMFI.BaseSequence;
+            Assert.AreEqual("PEPTIDE", b);
+            m = psmBaseSequenceMFI.MatchedFragmentIons;
+            Assert.AreEqual(1, m.Count);
+            psmBaseSequenceMFI.GetAminoAcidCoverage();
+
+            List<PeptideSpectralMatch> psms = new List<PeptideSpectralMatch>() { psmNoBaseSequenceNoMFI, psmBaseSequenceNoMFI, psmNoBaseSequenceMFI, psmBaseSequenceMFI };
+
+            foreach (var psm in psms)
+            {
+                psm.ToString();
+            }
+        }
+    }
+}