feat(parallelsearch): per-ID PEP for transient databases

MetaMorpheus/EngineLayer/FdrAnalysis/PEPAnalysisEngine.cs

@@ -1,4 +1,4 @@
-using Chemistry;
+using Chemistry;
 using Chromatography.RetentionTimePrediction;
 using Chromatography.RetentionTimePrediction.SSRCalc;
 using EngineLayer.CrosslinkSearch;
@@ -26,7 +26,7 @@ namespace EngineLayer
 {
     public class PepAnalysisEngine
     {
-        private int _randomSeed = 42;
+        protected int _randomSeed = 42;
 
         /// <summary>
         /// This method contains the hyper-parameters that will be used when training the machine learning model

MetaMorpheus/EngineLayer/ParallelSearch/TransientPepAnalysisEngine.cs

@@ -0,0 +1,171 @@
+#nullable enable
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using Chromatography.RetentionTimePrediction;
+using EngineLayer.FdrAnalysis;
+using Microsoft.ML;
+using Microsoft.ML.Data;
+
+// The concrete trained-model type returned by the FastTree pipeline; aliased so the train-once / apply
+// methods (this engine reuses one base-trained model across thousands of transient databases) stay readable.
+using PepModel = Microsoft.ML.Data.TransformerChain<Microsoft.ML.Data.BinaryPredictionTransformer<Microsoft.ML.Calibrators.CalibratedModelParametersBase<Microsoft.ML.Trainers.FastTree.FastTreeBinaryModelParameters, Microsoft.ML.Calibrators.PlattCalibrator>>>;
+
+namespace EngineLayer.ParallelSearch;
+
+/// <summary>
+/// ParallelSearch-specific PEP engine. The generic <see cref="PepAnalysisEngine"/> trains a fresh
+/// cross-validated model per dataset; the transient databases searched by ParallelSearchTask are far too
+/// small to support that. This derived engine trains ONE model on the large, decoy-rich base (human) search
+/// and REUSES it — plus a snapshotted background PEP -> PEP_QValue curve — to score every transient
+/// database's hits out-of-sample, without retraining. Keeping this behavior here leaves the base engine
+/// focused on the standard workflow (mirrors <see cref="TransientProteinParsimonyEngine"/> and
+/// <see cref="TransientProteinScoringAndFdrEngine"/>).
+/// </summary>
+public sealed class TransientPepAnalysisEngine : PepAnalysisEngine
+{
+    private MLContext? _trainedContext;
+    private PepModel? _trainedModel;
+
+    // Background PEP -> PEP_QValue relationship, snapshotted from the base (human) search, which has
+    // enough decoys to compute a real PEP-based q-value. Sorted by PEP ascending; PEP_QValue is monotone
+    // non-decreasing along it. Transient peptides (far too small for their own PEP target/decoy) get their
+    // PEP_QValue by looking up their MODEL PEP on this curve. This is deliberately PEP-based, NOT borrowed
+    // from the score-based QValue: q-value and PEP_QValue rank matches differently and are not interchangeable.
+    private double[] _bgPepAscPsm = Array.Empty<double>(), _bgQByPepPsm = Array.Empty<double>();
+    private double[] _bgPepAscPep = Array.Empty<double>(), _bgQByPepPep = Array.Empty<double>();
+
+    public TransientPepAnalysisEngine(List<SpectralMatch> psms, string searchType,
+        List<(string fileName, CommonParameters fileSpecificParameters)> fileSpecificParameters,
+        string outputFolder, IRetentionTimePredictor? rtPredictor = null)
+        : base(psms, searchType, fileSpecificParameters, outputFolder, rtPredictor)
+    {
+    }
+
+    /// <summary>True once <see cref="TrainSingleModelAndAssignBasePep"/> has produced a reusable model.</summary>
+    public bool HasTrainedModel => _trainedModel != null;
+
+    /// <summary>
+    /// Snapshots the (PEP ascending, PEP_QValue) curve over <paramref name="matches"/>, which must already
+    /// carry a PEP in GetFdrInfo(<paramref name="peptideLevel"/>).PEP and contain decoys. Computes the
+    /// PEP-based q-value on this set, then captures the two monotone arrays. Mutates the matches' PEP_QValue
+    /// and cumulative target/decoy (harmless — the score-based QValue lives in a different field).
+    /// </summary>
+    internal static (double[] pepAsc, double[] qByPep) BuildPepQValueCurve(List<SpectralMatch> matches, bool peptideLevel)
+    {
+        var ordered = matches.Where(m => m?.GetFdrInfo(peptideLevel) != null)
+                             .OrderBy(m => m.GetFdrInfo(peptideLevel).PEP).ToList();
+        if (ordered.Count == 0)
+            return (Array.Empty<double>(), Array.Empty<double>());
+        FdrAnalysisEngine.CalculateQValue(ordered, peptideLevelCalculation: peptideLevel, pepCalculation: true);
+        var pepAsc = new double[ordered.Count];
+        var qByPep = new double[ordered.Count];
+        for (int i = 0; i < ordered.Count; i++)
+        {
+            pepAsc[i] = ordered[i].GetFdrInfo(peptideLevel).PEP;
+            qByPep[i] = ordered[i].GetFdrInfo(peptideLevel).PEP_QValue;
+        }
+        return (pepAsc, qByPep);
+    }
+
+    /// <summary>
+    /// Maps a single PEP onto the background curve: returns the PEP_QValue of the background match whose PEP
+    /// is the smallest value >= <paramref name="pep"/> (lower-bound). Because PEP_QValue is monotone in PEP,
+    /// this is the q-value at that confidence level. Returns 2 (the FdrInfo sentinel) when no curve exists.
+    /// </summary>
+    internal static double LookupBackgroundPepQValue(double pep, double[] pepAsc, double[] qByPep)
+    {
+        if (pepAsc.Length == 0)
+            return 2.0;
+        int l = 0, r = pepAsc.Length;
+        while (l < r)
+        {
+            int mid = (l + r) >> 1;
+            if (pepAsc[mid] >= pep) r = mid; else l = mid + 1;
+        }
+        int idx = l < pepAsc.Length ? l : pepAsc.Length - 1;
+        return qByPep[idx];
+    }
+
+    /// <summary>
+    /// Trains ONE PEP model on all of this engine's PSMs (the base search) and assigns their PEP, then keeps
+    /// the model + ML context so it can be REUSED — via <see cref="AssignPepFromTrainedModel"/> — to score
+    /// out-of-sample PSMs (e.g. the per-database transient hits) WITHOUT retraining. The transient databases
+    /// are far too small to train their own model and are out-of-sample relative to this one, so the
+    /// cross-validation used by <see cref="PepAnalysisEngine.ComputePEPValuesForAllPSMs"/> is unnecessary for
+    /// them. The feature dictionaries were already built from these PSMs in the constructor, so transient
+    /// feature vectors are computed against the same (base-run) calibration. Returns false when the base PSMs
+    /// lack both target and decoy training examples.
+    /// </summary>
+    public bool TrainSingleModelAndAssignBasePep()
+    {
+        List<SpectralMatchGroup> peptideGroups = UsePeptideLevelQValueForTraining
+            ? SpectralMatchGroup.GroupByBaseSequence(AllPsms)
+            : SpectralMatchGroup.GroupByIndividualPsm(AllPsms);
+        var allIndices = Enumerable.Range(0, peptideGroups.Count).ToList();
+
+        var psmData = CreatePsmData(SearchType, peptideGroups, allIndices).ToList();
+        if (!psmData.Any(p => p.Label) || !psmData.Any(p => !p.Label))
+            return false; // need both positive (target) and negative (decoy) examples
+
+        var mlContext = new MLContext(seed: _randomSeed);
+        var trainer = mlContext.BinaryClassification.Trainers.FastTree(BGDTreeOptions);
+        var pipeline = mlContext.Transforms.Concatenate("Features", TrainingVariables).Append(trainer);
+        PepModel model = pipeline.Fit(mlContext.Data.LoadFromEnumerable(psmData));
+        _trainedContext = mlContext;
+        _trainedModel = model;
+
+        // Assign PEP to the base PSMs as well (the base assignment is for the base-search output only;
+        // transient PSMs scored later are genuinely out-of-sample, so no out-of-fold concern applies).
+        Compute_PSM_PEP(peptideGroups, allIndices, mlContext, model, SearchType, OutputFolder);
+
+        // Snapshot the background PEP -> PEP_QValue curves used to assign transient PEP_QValue by lookup.
+        // PSM-level: over all base PSMs. Peptide-level: over one representative (best-scoring) PSM per base
+        // sequence, so the q-value reflects unique peptides. Both read PEP (identical in both FdrInfos).
+        (_bgPepAscPsm, _bgQByPepPsm) = BuildPepQValueCurve(AllPsms, peptideLevel: false);
+        var basePeptideReps = SpectralMatchGroup.GroupByBaseSequence(AllPsms)
+            .Select(g => g.OrderByDescending(p => p.Score).First())
+            .ToList();
+        (_bgPepAscPep, _bgQByPepPep) = BuildPepQValueCurve(basePeptideReps, peptideLevel: false);
+        return true;
+    }
+
+    /// <summary>
+    /// Assigns PEP to a NEW set of PSMs (a transient database's hits) using the model trained on the base
+    /// search by <see cref="TrainSingleModelAndAssignBasePep"/>. Reuses the base-trained model and the base
+    /// feature dictionaries — no retraining. Safe to call from many databases concurrently. No-op until a
+    /// model has been trained.
+    /// </summary>
+    public void AssignPepFromTrainedModel(List<SpectralMatch> psms, bool peptideLevel = false)
+    {
+        if (_trainedModel == null || psms == null)
+            return;
+        var scorable = psms.Where(p => p != null).ToList();
+        if (scorable.Count == 0)
+            return;
+        // Compute_PSM_PEP writes BOTH PsmFdrInfo.PEP and PeptideFdrInfo.PEP. Transient PSMs may not have a
+        // PeptideFdrInfo yet (peptide-level FDR isn't always run per database), so create whichever is
+        // missing — otherwise PEP is silently never assigned (the old guard skipped these PSMs entirely).
+        foreach (var p in scorable)
+        {
+            p.PsmFdrInfo ??= new FdrInfo();
+            p.PeptideFdrInfo ??= new FdrInfo();
+        }
+        List<SpectralMatchGroup> groups = UsePeptideLevelQValueForTraining
+            ? SpectralMatchGroup.GroupByBaseSequence(scorable)
+            : SpectralMatchGroup.GroupByIndividualPsm(scorable);
+        Compute_PSM_PEP(groups, Enumerable.Range(0, groups.Count).ToList(), _trainedContext!, _trainedModel, SearchType, OutputFolder);
+
+        // Assign PEP_QValue by mapping each match's MODEL PEP through the BACKGROUND PEP->PEP_QValue curve.
+        // The transient database is far too small for its own PEP target/decoy (few/no decoys), so we borrow
+        // the relationship from the base search — but strictly on PEP, not on the score-based QValue, since
+        // q-value and PEP_QValue rank matches differently and are not interchangeable.
+        double[] pepAsc = peptideLevel ? _bgPepAscPep : _bgPepAscPsm;
+        double[] qByPep = peptideLevel ? _bgQByPepPep : _bgQByPepPsm;
+        foreach (var p in scorable)
+        {
+            var info = p.GetFdrInfo(peptideLevel);
+            info.PEP_QValue = LookupBackgroundPepQValue(info.PEP, pepAsc, qByPep);
+        }
+    }
+}

MetaMorpheus/TaskLayer/ParallelSearch/Analysis/Collectors/BasicMetricCollector.cs

@@ -1,4 +1,4 @@
-#nullable enable
+#nullable enable
 using System;
 using System.Collections.Generic;
 using System.Linq;
@@ -22,6 +22,14 @@ public class BasicMetricCollector : IMetricCollector
     public const string TargetPeptidesFromTransientDb = "TargetPeptidesFromTransientDb";
     public const string TargetPeptidesFromTransientDbAtQValueThreshold = "TargetPeptidesFromTransientDbAtQValueThreshold";
 
+    // PEP-based confident counts, reported at BOTH 1% and 5% PEP_QValue. PEP_QValue is the PEP-ranked q-value
+    // assigned by mapping each match's model PEP onto the background curve (see PepAnalysisEngine); it is a
+    // different confidence axis from the score-based QValue above, so both are reported side by side.
+    public const string TargetPsmsFromTransientDbAtPepQ01 = "TargetPsmsFromTransientDbAtPepQ01";
+    public const string TargetPsmsFromTransientDbAtPepQ05 = "TargetPsmsFromTransientDbAtPepQ05";
+    public const string TargetPeptidesFromTransientDbAtPepQ01 = "TargetPeptidesFromTransientDbAtPepQ01";
+    public const string TargetPeptidesFromTransientDbAtPepQ05 = "TargetPeptidesFromTransientDbAtPepQ05";
+
     public string CollectorName => "ResultCount";
 
     public IEnumerable<string> GetOutputColumns()
@@ -35,6 +43,10 @@ public IEnumerable<string> GetOutputColumns()
         yield return TargetPeptidesAtQValueThreshold;
         yield return TargetPeptidesFromTransientDb;
         yield return TargetPeptidesFromTransientDbAtQValueThreshold;
+        yield return TargetPsmsFromTransientDbAtPepQ01;
+        yield return TargetPsmsFromTransientDbAtPepQ05;
+        yield return TargetPeptidesFromTransientDbAtPepQ01;
+        yield return TargetPeptidesFromTransientDbAtPepQ05;
     }
 
     public bool CanCollectData(TransientDatabaseContext context)
@@ -61,7 +73,14 @@ public Dictionary<string, object> CollectData(TransientDatabaseContext context)
 
             [TargetPeptidesAtQValueThreshold] = context.AllPeptides.Count(p => !p.IsDecoy && p.PeptideFdrInfo?.QValue <= qValueThreshold),
             [TargetPeptidesFromTransientDb] = context.TransientPeptides.Count(p => !p.IsDecoy),
-            [TargetPeptidesFromTransientDbAtQValueThreshold] = context.TransientPeptides.Count(p => !p.IsDecoy && p.PeptideFdrInfo?.QValue <= qValueThreshold)
+            [TargetPeptidesFromTransientDbAtQValueThreshold] = context.TransientPeptides.Count(p => !p.IsDecoy && p.PeptideFdrInfo?.QValue <= qValueThreshold),
+
+            // PEP-based confident counts at 1% and 5%. PEP_QValue is only meaningful once a PEP model has been
+            // trained (transient matches default to the sentinel 2 otherwise), so these read 0 when PEP is off.
+            [TargetPsmsFromTransientDbAtPepQ01] = context.TransientPsms.Count(p => !p.IsDecoy && p.GetFdrInfo(false)?.PEP_QValue < 0.01),
+            [TargetPsmsFromTransientDbAtPepQ05] = context.TransientPsms.Count(p => !p.IsDecoy && p.GetFdrInfo(false)?.PEP_QValue < 0.05),
+            [TargetPeptidesFromTransientDbAtPepQ01] = context.TransientPeptides.Count(p => !p.IsDecoy && p.PeptideFdrInfo?.PEP_QValue < 0.01),
+            [TargetPeptidesFromTransientDbAtPepQ05] = context.TransientPeptides.Count(p => !p.IsDecoy && p.PeptideFdrInfo?.PEP_QValue < 0.05)
         };
     }
 }

MetaMorpheus/TaskLayer/ParallelSearch/Analysis/TransientDatabaseMetrics.cs

@@ -61,7 +61,13 @@ public TransientDatabaseMetrics() { }
     public int TargetPeptidesAtQValueThreshold { get; set; }
     public int TargetPeptidesFromTransientDb { get; set; }
     public int TargetPeptidesFromTransientDbAtQValueThreshold { get; set; }
-
+
+    // PEP-based confident counts, reported at both 1% and 5% PEP_QValue (distinct confidence axis from QValue).
+    [Optional] public int TargetPsmsFromTransientDbAtPepQ01 { get; set; }
+    [Optional] public int TargetPsmsFromTransientDbAtPepQ05 { get; set; }
+    [Optional] public int TargetPeptidesFromTransientDbAtPepQ01 { get; set; }
+    [Optional] public int TargetPeptidesFromTransientDbAtPepQ05 { get; set; }
+
     // Protein group metrics (0 if parsimony not run)
     public int TargetProteinGroupsAtQValueThreshold { get; set; }
     public int TargetProteinGroupsFromTransientDb { get; set; }
@@ -365,6 +371,10 @@ public void PopulateResultsFromProperties()
         Results[BasicMetricCollector.TargetPeptidesAtQValueThreshold] = TargetPeptidesAtQValueThreshold;
         Results[BasicMetricCollector.TargetPeptidesFromTransientDb] = TargetPeptidesFromTransientDb;
         Results[BasicMetricCollector.TargetPeptidesFromTransientDbAtQValueThreshold] = TargetPeptidesFromTransientDbAtQValueThreshold;
+        Results[BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ01] = TargetPsmsFromTransientDbAtPepQ01;
+        Results[BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ05] = TargetPsmsFromTransientDbAtPepQ05;
+        Results[BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ01] = TargetPeptidesFromTransientDbAtPepQ01;
+        Results[BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ05] = TargetPeptidesFromTransientDbAtPepQ05;
         Results[ProteinGroupCollector.TargetProteinGroupsAtQValueThreshold] = TargetProteinGroupsAtQValueThreshold;
         Results[ProteinGroupCollector.TargetProteinGroupsFromTransientDb] = TargetProteinGroupsFromTransientDb;
         Results[ProteinGroupCollector.TargetProteinGroupsFromTransientDbAtQValueThreshold] = TargetProteinGroupsFromTransientDbAtQValueThreshold;
@@ -486,6 +496,10 @@ public void PopulatePropertiesFromResults()
         TargetPeptidesAtQValueThreshold = GetValue<int>(BasicMetricCollector.TargetPeptidesAtQValueThreshold);
         TargetPeptidesFromTransientDb = GetValue<int>(BasicMetricCollector.TargetPeptidesFromTransientDb);
         TargetPeptidesFromTransientDbAtQValueThreshold = GetValue<int>(BasicMetricCollector.TargetPeptidesFromTransientDbAtQValueThreshold);
+        TargetPsmsFromTransientDbAtPepQ01 = GetValue<int>(BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ01);
+        TargetPsmsFromTransientDbAtPepQ05 = GetValue<int>(BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ05);
+        TargetPeptidesFromTransientDbAtPepQ01 = GetValue<int>(BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ01);
+        TargetPeptidesFromTransientDbAtPepQ05 = GetValue<int>(BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ05);
         TargetProteinGroupsAtQValueThreshold = GetValue<int>(ProteinGroupCollector.TargetProteinGroupsAtQValueThreshold);
         TargetProteinGroupsFromTransientDb = GetValue<int>(ProteinGroupCollector.TargetProteinGroupsFromTransientDb);
         TargetProteinGroupsFromTransientDbAtQValueThreshold = GetValue<int>(ProteinGroupCollector.TargetProteinGroupsFromTransientDbAtQValueThreshold);

MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs

@@ -121,6 +121,9 @@ public override SearchParameters SearchParameters
 
     #endregion
 
+    // PEP model TRAINED ONCE on the base (human) search and reused to assign a PEP to every transient
+    // database's PSMs (they're too small to train their own model, and are out-of-sample vs the base).
+    [TomlIgnore] private TransientPepAnalysisEngine? _pepEngine;
     protected override MyTaskResults RunSpecific(string outputFolder,
         List<DbForTask> dbFilenameList, List<string> currentRawFileList,
         string taskId, FileSpecificParameters[] fileSettingsList)
@@ -393,6 +396,19 @@ private void Initialize(string taskId, List<DbForTask> dbFilenameList,
         ProseCreatedWhileRunning.Append(
             $"Searching {ParallelSearchParameters.TransientDatabases.Count} transient databases against {currentRawFileList.Count} spectra files. ");
 
+        // PEP: train ONE model on the base search PSMs (decoy-rich, large) before writing the base output
+        // files, so the base PSM/peptide/protein files include PEP and PEP_QValue.
+        var trainingPsms = baselinePsms
+            .Where(p => !p.IsDecoy && p.PsmFdrInfo?.QValue <= CommonParameters.QValueThreshold)
+            .ToList();
+        if (trainingPsms.Count >= 100)
+        {
+            var engine = new TransientPepAnalysisEngine(
+                trainingPsms, "standard", FileSpecificParameters, outputFolder, null);
+            if (engine.TrainSingleModelAndAssignBasePep())
+                _pepEngine = engine;
+        }
+
         Task.WhenAll([psmTask, peptideTask, proteinTask]).Wait();
     }
 
@@ -707,6 +723,16 @@ private void PerformSearch(List<IBioPolymer> proteinsToSearch, SpectralMatch[] s
                 .ToList();
         }
 
+        // PEP: assign each transient PSM/peptide a posterior error probability + PEP_QValue (mapped onto the
+        // background curve) BEFORE the metric collectors and statistics run, so confident counts and family
+        // tests can use PEP_QValue. Runs after the FDR alignment so the borrowed score-based QValue is in place.
+        if (_pepEngine != null)
+        {
+            _pepEngine.AssignPepFromTrainedModel(transientPsms, peptideLevel: false);
+            if (transientPeptides.Count > 0)
+                _pepEngine.AssignPepFromTrainedModel(transientPeptides, peptideLevel: true);
+        }
+
         #region Result Caching and Analysis
 
         // Create analysis context

MetaMorpheus/TaskLayer/Properties/AssemblyInfo.cs

@@ -0,0 +1,2 @@
+// Exposes internal members (e.g. extracted pure decision helpers) to the Test assembly for unit testing.
+[assembly: System.Runtime.CompilerServices.InternalsVisibleTo("Test")]