diff --git a/MetaMorpheus/CMD/CMD.csproj b/MetaMorpheus/CMD/CMD.csproj
index a60989660..78cdfd82b 100644
--- a/MetaMorpheus/CMD/CMD.csproj
+++ b/MetaMorpheus/CMD/CMD.csproj
@@ -8,6 +8,12 @@
full
true
+
+ true
@@ -24,7 +30,7 @@
-
+
diff --git a/MetaMorpheus/EngineLayer/EngineLayer.csproj b/MetaMorpheus/EngineLayer/EngineLayer.csproj
index 06e043759..987ce2d8f 100644
--- a/MetaMorpheus/EngineLayer/EngineLayer.csproj
+++ b/MetaMorpheus/EngineLayer/EngineLayer.csproj
@@ -29,7 +29,7 @@
-
+
diff --git a/MetaMorpheus/EngineLayer/FdrAnalysis/PEPAnalysisEngine.cs b/MetaMorpheus/EngineLayer/FdrAnalysis/PEPAnalysisEngine.cs
index 0dbf45b94..70ee8ded3 100644
--- a/MetaMorpheus/EngineLayer/FdrAnalysis/PEPAnalysisEngine.cs
+++ b/MetaMorpheus/EngineLayer/FdrAnalysis/PEPAnalysisEngine.cs
@@ -26,7 +26,7 @@ namespace EngineLayer
{
public class PepAnalysisEngine
{
- private int _randomSeed = 42;
+ protected int _randomSeed = 42;
///
/// This method contains the hyper-parameters that will be used when training the machine learning model
diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/MslPeptideReader.cs b/MetaMorpheus/EngineLayer/ParallelSearch/MslPeptideReader.cs
new file mode 100644
index 000000000..01e2f83e0
--- /dev/null
+++ b/MetaMorpheus/EngineLayer/ParallelSearch/MslPeptideReader.cs
@@ -0,0 +1,324 @@
+using System;
+using System.Collections.Generic;
+using System.Text;
+using Chemistry;
+using Omics;
+using Omics.Fragmentation;
+using Omics.Fragmentation.Peptide;
+using Omics.SpectralMatch.MslSpectralLibrary;
+using Proteomics;
+using Proteomics.ProteolyticDigestion;
+using Readers.SpectralLibrary;
+
+namespace EngineLayer.ParallelSearch
+{
+ ///
+ /// Reads a MetaMorpheus .msl spectral library as a PEPTIDE SOURCE for the parallel search,
+ /// using the library the way the binary format was designed for:
+ /// keeps fragments on disk; the precursor index (mass / charge / iRT metadata, NO fragment I/O) is
+ /// filtered on TWO orthogonal axes — precursor mass AND retention time — and only the surviving
+ /// CANDIDATES have their fragments fetched on demand.
+ ///
+ /// Retention time is the axis that makes the filter selective: precursor mass alone is degenerate when
+ /// tens of thousands of dense experimental precursors saturate the mass window, but a peptide's predicted
+ /// elution (stored as Chronologer iRT in the .msl, calibrated to this run from the base search) collides
+ /// with far fewer scans. Entries with iRT==0 (predictor could not place them) fall back to mass-only so
+ /// they are never lost.
+ ///
+ /// Each candidate's stored float fragments are turned into s directly (no
+ /// re-fragmentation). Accession handling: one shared per accession (concatenation
+ /// of its candidate peptides with per-peptide offsets) so parsimony/coverage stay correct and in-bounds.
+ ///
+ public static class MslPeptideReader
+ {
+ private const double MassPreFilterMarginDa = 0.01;
+
+ /// RT↔iRT calibration and learned precursor tolerance applied to the candidate pre-filter.
+ public readonly struct CandidatePriors
+ {
+ public readonly double[] SortedScanMasses; // experimental precursor masses, ASCENDING
+ public readonly double[] ScanRetentionTimes; // RT of each scan, SAME order as SortedScanMasses
+ public readonly double PrecursorTolPpm; // learned precursor tolerance (ppm)
+ public readonly double RtSlope, RtIntercept; // observedRT = RtSlope*iRT + RtIntercept
+ public readonly double RtWindowMin; // +/- window in observed-RT minutes (k * residualSD)
+
+ public CandidatePriors(double[] sortedScanMasses, double[] scanRetentionTimes,
+ double precursorTolPpm, double rtSlope, double rtIntercept, double rtWindowMin)
+ {
+ SortedScanMasses = sortedScanMasses;
+ ScanRetentionTimes = scanRetentionTimes;
+ PrecursorTolPpm = precursorTolPpm;
+ RtSlope = rtSlope;
+ RtIntercept = rtIntercept;
+ RtWindowMin = rtWindowMin;
+ }
+ }
+
+ private readonly struct PendingPeptide
+ {
+ public readonly string FullSequence;
+ public readonly int Start;
+ public readonly int Length;
+ public readonly List Fragments;
+ public PendingPeptide(string fullSequence, int start, int length, List fragments)
+ { FullSequence = fullSequence; Start = start; Length = length; Fragments = fragments; }
+ }
+
+ private sealed class AccessionGroup
+ {
+ public readonly StringBuilder Sequence = new();
+ public readonly List Peptides = new();
+ public bool IsDecoy;
+ }
+
+ public static List<(IBioPolymerWithSetMods Peptide, List Fragments)> ReadPeptides(
+ string mslPath, string databaseName, CandidatePriors priors)
+ {
+ var mods = GlobalVariables.AllModsKnownDictionary;
+ var digestionParams = new DigestionParams("trypsin", maxMissedCleavages: 0);
+ string fallbackAccession = $"{databaseName}_UNKNOWN";
+ var groups = new Dictionary();
+
+ int totalEntries = 0, massCandidates = 0, massRtCandidates = 0;
+ var sw = System.Diagnostics.Stopwatch.StartNew();
+
+ long loadMs, getEntryTicks = 0, buildTicks = 0;
+ var swPhase = System.Diagnostics.Stopwatch.StartNew();
+ var candidateIdx = new List();
+ using (var library = MslLibrary.LoadIndexOnly(mslPath))
+ {
+ loadMs = swPhase.ElapsedMilliseconds; // cost of LoadIndexOnly (metadata read + sort)
+
+ // PASS 1 — filter on mass + RT using ONLY the in-memory index (no fragment I/O).
+ foreach (var idx in library.QueryMzWindow(float.MinValue, float.MaxValue))
+ {
+ totalEntries++;
+ if (IsCandidate(idx.PrecursorMz, idx.Charge, idx.Irt, priors, out bool massMatched))
+ {
+ massRtCandidates++;
+ candidateIdx.Add(idx.PrecursorIdx);
+ }
+ if (massMatched)
+ massCandidates++;
+ }
+
+ // PASS 2 — fetch fragments in ON-DISK (PrecursorIdx/write) order, not the m/z order the
+ // window query yielded. Fragment blocks are stored in PrecursorIdx order, so a sorted walk
+ // reads the file front-to-back (sequential, OS read-ahead) instead of one random seek per
+ // candidate — the dominant reader cost. Result order doesn't affect IDs/scores; the shared
+ // protein concatenation just follows this order.
+ candidateIdx.Sort();
+ foreach (int pid in candidateIdx)
+ {
+ long t0 = swPhase.ElapsedTicks;
+ MslLibraryEntry entry = library.GetEntry(pid);
+ getEntryTicks += swPhase.ElapsedTicks - t0;
+ if (entry is null) continue;
+
+ string fullSequence = entry.FullSequence;
+ string baseSequence = IBioPolymerWithSetMods.GetBaseSequenceFromFullSequence(fullSequence);
+ string accession = string.IsNullOrEmpty(entry.ProteinAccession) ? fallbackAccession : entry.ProteinAccession;
+
+ if (!groups.TryGetValue(accession, out var group))
+ {
+ group = new AccessionGroup { IsDecoy = entry.IsDecoy };
+ groups[accession] = group;
+ }
+
+ int start = group.Sequence.Length + 1;
+ group.Sequence.Append(baseSequence);
+ long t1 = swPhase.ElapsedTicks;
+ // Lean libraries store no fragments → return null so the engine fragments the peptide on
+ // the fly (cheap, and matches the search's exact dissociation/precision). Full libraries
+ // build Products from the stored float ions.
+ List built = entry.MatchedFragmentIons is { Count: > 0 }
+ ? BuildProducts(entry.MatchedFragmentIons)
+ : null;
+ buildTicks += swPhase.ElapsedTicks - t1;
+ group.Peptides.Add(new PendingPeptide(fullSequence, start, baseSequence.Length, built));
+ }
+ }
+ double getEntryMs = getEntryTicks * 1000.0 / System.Diagnostics.Stopwatch.Frequency;
+ double buildMs = buildTicks * 1000.0 / System.Diagnostics.Stopwatch.Frequency;
+
+ // Per-database candidate-filter diagnostics. Opt-in (MM_PARALLELSEARCH_DIAG=1) so it is
+ // available when profiling the .msl reader but silent in normal runs.
+ if (Environment.GetEnvironmentVariable("MM_PARALLELSEARCH_DIAG") == "1")
+ Console.WriteLine($"MSL_RT {databaseName}: entries={totalEntries} massCand={massCandidates} " +
+ $"({Pct(massCandidates, totalEntries)}%) massRtCand={massRtCandidates} ({Pct(massRtCandidates, totalEntries)}%) " +
+ $"load={loadMs}ms getEntry={getEntryMs:F0}ms build={buildMs:F0}ms " +
+ $"tolPpm={priors.PrecursorTolPpm:F1} rtWin=+/-{priors.RtWindowMin:F1}min ms={sw.ElapsedMilliseconds}");
+
+ var result = new List<(IBioPolymerWithSetMods, List)>();
+ foreach (var kvp in groups)
+ {
+ var protein = new Protein(kvp.Value.Sequence.ToString(), kvp.Key, isDecoy: kvp.Value.IsDecoy);
+ foreach (var pending in kvp.Value.Peptides)
+ {
+ var peptide = new PeptideWithSetModifications(
+ pending.FullSequence, mods, p: protein, digestionParams: digestionParams,
+ oneBasedStartResidueInProtein: pending.Start,
+ oneBasedEndResidueInProtein: pending.Start + pending.Length - 1,
+ missedCleavages: 0);
+ result.Add((peptide, pending.Fragments));
+ }
+ }
+ return result;
+ }
+
+ ///
+ /// Reads a MERGED index (many databases in one .msl, each entry's accession stamped "db|accession")
+ /// with ONE open, runs the mass+RT candidate filter ONCE over all entries, and returns the surviving
+ /// candidate peptides GROUPED BY their source database. The per-database lists are then searched
+ /// INDEPENDENTLY (each against the shared base PSMs in its own copy) exactly as in the per-file path —
+ /// databases never compete; this only collapses 1000s of file opens into one shared in-memory index.
+ ///
+ public static Dictionary Fragments)>>
+ ReadCandidatesGroupedByDatabase(string mergedMslPath, CandidatePriors priors)
+ {
+ var mods = GlobalVariables.AllModsKnownDictionary;
+ var digestionParams = new DigestionParams("trypsin", maxMissedCleavages: 0);
+
+ // dbTag -> (accession-within-db -> shared-protein layout)
+ var byDb = new Dictionary>();
+ var candidateIdx = new List();
+ int totalEntries = 0, massRtCandidates = 0;
+ var swPhase = System.Diagnostics.Stopwatch.StartNew();
+ long loadMs;
+
+ using (var library = MslLibrary.LoadIndexOnly(mergedMslPath))
+ {
+ loadMs = swPhase.ElapsedMilliseconds;
+ // PASS 1 — mass+RT filter over the whole merged index (no fragment I/O).
+ foreach (var idx in library.QueryMzWindow(float.MinValue, float.MaxValue))
+ {
+ totalEntries++;
+ if (IsCandidate(idx.PrecursorMz, idx.Charge, idx.Irt, priors, out _))
+ {
+ massRtCandidates++;
+ candidateIdx.Add(idx.PrecursorIdx);
+ }
+ }
+
+ // PASS 2 — fetch candidates in on-disk order (sequential), group by source database.
+ candidateIdx.Sort();
+ foreach (int pid in candidateIdx)
+ {
+ MslLibraryEntry entry = library.GetEntry(pid);
+ if (entry is null) continue;
+
+ var (dbTag, accession) = ParseDbTagAndAccession(entry.ProteinAccession);
+
+ string fullSequence = entry.FullSequence;
+ string baseSequence = IBioPolymerWithSetMods.GetBaseSequenceFromFullSequence(fullSequence);
+
+ if (!byDb.TryGetValue(dbTag, out var accGroups))
+ {
+ accGroups = new Dictionary();
+ byDb[dbTag] = accGroups;
+ }
+ if (!accGroups.TryGetValue(accession, out var group))
+ {
+ group = new AccessionGroup { IsDecoy = entry.IsDecoy };
+ accGroups[accession] = group;
+ }
+ int start = group.Sequence.Length + 1;
+ group.Sequence.Append(baseSequence);
+ List built = entry.MatchedFragmentIons is { Count: > 0 } ? BuildProducts(entry.MatchedFragmentIons) : null;
+ group.Peptides.Add(new PendingPeptide(fullSequence, start, baseSequence.Length, built));
+ }
+ }
+
+ if (Environment.GetEnvironmentVariable("MM_PARALLELSEARCH_DIAG") == "1")
+ Console.WriteLine($"MSL_MERGED: entries={totalEntries} candidates={massRtCandidates} " +
+ $"({Pct(massRtCandidates, totalEntries)}%) databases={byDb.Count} load={loadMs}ms total={swPhase.ElapsedMilliseconds}ms");
+
+ // Materialize per-database peptide lists (each database is searched independently downstream).
+ var result = new Dictionary)>>(byDb.Count);
+ foreach (var dbKvp in byDb)
+ {
+ var list = new List<(IBioPolymerWithSetMods, List)>();
+ foreach (var accKvp in dbKvp.Value)
+ {
+ var protein = new Protein(accKvp.Value.Sequence.ToString(), accKvp.Key, isDecoy: accKvp.Value.IsDecoy);
+ foreach (var pending in accKvp.Value.Peptides)
+ {
+ var peptide = new PeptideWithSetModifications(
+ pending.FullSequence, mods, p: protein, digestionParams: digestionParams,
+ oneBasedStartResidueInProtein: pending.Start,
+ oneBasedEndResidueInProtein: pending.Start + pending.Length - 1,
+ missedCleavages: 0);
+ list.Add((peptide, pending.Fragments));
+ }
+ }
+ result[dbKvp.Key] = list;
+ }
+ return result;
+ }
+
+ private static double Pct(int n, int d) => d == 0 ? 0 : Math.Round(100.0 * n / d, 1);
+
+ /// First index i where sorted[i] >= value (sorted ascending).
+ internal static int LowerBound(double[] sorted, double value)
+ {
+ int i = Array.BinarySearch(sorted, value);
+ return i < 0 ? ~i : i;
+ }
+
+ ///
+ /// The .msl candidate pre-filter (extracted for testing): does any experimental scan match this library
+ /// entry on BOTH precursor mass (within the learned ppm window) AND retention time (within the calibrated
+ /// window of the entry's predicted RT)? An iRT of 0 (unpredicted) keeps the entry on mass alone.
+ /// reports whether the mass axis alone matched (for diagnostics).
+ ///
+ internal static bool IsCandidate(double precursorMz, int charge, float irt, in CandidatePriors priors, out bool massMatched)
+ {
+ massMatched = false;
+ double neutralMass = precursorMz.ToMass(charge);
+ double tolDa = neutralMass * priors.PrecursorTolPpm * 1e-6 + MassPreFilterMarginDa;
+ int a = LowerBound(priors.SortedScanMasses, neutralMass - tolDa);
+ double hi = neutralMass + tolDa;
+ if (a >= priors.SortedScanMasses.Length || priors.SortedScanMasses[a] > hi)
+ return false; // no mass match at all
+ massMatched = true;
+
+ if (irt == 0f)
+ return true; // unpredicted RT -> keep on mass alone
+ double predRt = priors.RtSlope * irt + priors.RtIntercept;
+ for (int i = a; i < priors.SortedScanMasses.Length && priors.SortedScanMasses[i] <= hi; i++)
+ if (Math.Abs(priors.ScanRetentionTimes[i] - predRt) <= priors.RtWindowMin)
+ return true;
+ return false;
+ }
+
+ ///
+ /// Splits a merged-index accession stamped "db|accession" into its (dbTag, accession) parts. An entry
+ /// with no bar falls back to dbTag "UNKNOWN"; an empty accession becomes "<dbTag>_UNKNOWN".
+ ///
+ internal static (string dbTag, string accession) ParseDbTagAndAccession(string tagged)
+ {
+ tagged ??= "";
+ int bar = tagged.IndexOf('|');
+ string dbTag = bar > 0 ? tagged.Substring(0, bar) : "UNKNOWN";
+ string accession = bar >= 0 ? tagged.Substring(bar + 1) : tagged;
+ if (string.IsNullOrEmpty(accession))
+ accession = dbTag + "_UNKNOWN";
+ return (dbTag, accession);
+ }
+
+ private static List BuildProducts(List ions)
+ {
+ var products = new List(ions.Count);
+ foreach (var ion in ions)
+ {
+ FragmentationTerminus terminus =
+ TerminusSpecificProductTypes.ProductTypeToFragmentationTerminus.TryGetValue(
+ ion.ProductType, out var t) ? t : FragmentationTerminus.None;
+ products.Add(new Product(
+ ion.ProductType, terminus, ion.Mz.ToMass(ion.Charge), ion.FragmentNumber,
+ ion.ResiduePosition, ion.NeutralLoss, ion.SecondaryProductType, ion.SecondaryFragmentNumber));
+ }
+ return products;
+ }
+ }
+}
diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs b/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs
index 7e8a6bb1a..e57a78116 100644
--- a/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs
+++ b/MetaMorpheus/EngineLayer/ParallelSearch/TransientClassicSearchEngine.cs
@@ -4,11 +4,9 @@
using System.Linq;
using System.Threading;
using System.Threading.Tasks;
-using Chemistry;
using EngineLayer.ClassicSearch;
using EngineLayer.FdrAnalysis;
using EngineLayer.Util;
-using MzLibUtil;
using Omics;
using Omics.Fragmentation;
using Omics.Modifications;
@@ -21,7 +19,7 @@ namespace EngineLayer.ParallelSearch
/// The base PSM list is shared across threads, and only the PSMs that are updated by the transient search are cloned and modified, while the rest of the PSMs remain shared and read-only.
/// Features that are removed from classic search engine for runtime efficiency and memory efficiency:
/// - Some detailed logging
- /// - Mass Differance acceptor is always exact.
+ /// - Mass Differance acceptor is always exact.
///
public class TransientClassicSearchEngine : ClassicSearchEngine
{
@@ -29,15 +27,20 @@ public class TransientClassicSearchEngine : ClassicSearchEngine
private bool _singleThreadMode;
private readonly ConcurrentDictionary UpdatedIndexes = new();
private readonly bool _copyOnWriteEnabled;
+ // When supplied (e.g. from a .msl spectral library), the search iterates these peptides
+ // directly instead of digesting Proteins, and matches each peptide's STORED (float) fragments
+ // instead of re-fragmenting — skipping both digestion and fragmentation.
+ private readonly List<(IBioPolymerWithSetMods Peptide, List Fragments)> _precomputedPeptides;
public TransientClassicSearchEngine(SpectralMatch[] globalPsms, Ms2ScanWithSpecificMass[] arrayOfSortedMS2Scans,
- List variableModifications, List fixedModifications,
+ List variableModifications, List fixedModifications,
List proteinList, MassDiffAcceptor searchMode, CommonParameters commonParameters, List<(string FileName, CommonParameters Parameters)> fileSpecificParameters, List nestedIds,
- bool copyOnWriteEnabled = false)
+ bool copyOnWriteEnabled = false, List<(IBioPolymerWithSetMods Peptide, List Fragments)> precomputedPeptides = null)
: base(globalPsms, arrayOfSortedMS2Scans, variableModifications, fixedModifications, null, null, null, proteinList, searchMode, commonParameters, fileSpecificParameters, null, nestedIds, false)
{
UpdatedIndexes = new ConcurrentDictionary();
_singleThreadMode = CommonParameters.MaxThreadsToUsePerFile <= 1;
_copyOnWriteEnabled = copyOnWriteEnabled;
+ _precomputedPeptides = precomputedPeptides;
if (!_singleThreadMode)
{
@@ -58,79 +61,56 @@ protected override MetaMorpheusEngineResults RunSpecific()
Status("Performing classic search...");
- if (Proteins.Any())
+ bool usePrecomputedPeptides = _precomputedPeptides != null && _precomputedPeptides.Count > 0;
+ if (Proteins.Any() || usePrecomputedPeptides)
{
+ // Match one peptide against every candidate scan and update PSMs directly, using the shared
+ // MatchFragmentIons / CalculatePeptideScore so MatchedFragmentIon stays the type at the engine
+ // boundary. The two memory/runtime wins over ClassicSearchEngine are kept: the base PSM list is
+ // shared (only updated PSMs are cloned — see AddPeptideCandidateToPsm), and a .msl library
+ // supplies precomputed fragments so the Fragment() call is skipped entirely.
+ void ProcessOnePeptide(IBioPolymerWithSetMods peptide, List peptideTheorProducts,
+ List precomputedFragments = null)
+ {
+ Interlocked.Increment(ref peptideCounter);
+
+ List products;
+ if (precomputedFragments != null)
+ {
+ products = precomputedFragments;
+ }
+ else
+ {
+ peptideTheorProducts.Clear();
+ peptide.Fragment(CommonParameters.DissociationType, CommonParameters.DigestionParams.FragmentationTerminus, peptideTheorProducts, CommonParameters.FragmentationParameters);
+ products = peptideTheorProducts;
+ }
+
+ foreach (ScanWithIndexAndNotchInfo scan in GetAcceptableScans(peptide.MonoisotopicMass, SearchMode))
+ {
+ Ms2ScanWithSpecificMass theScan = ArrayOfSortedMS2Scans[scan.ScanIndex];
+ List matchedIons = MatchFragmentIons(theScan, products, CommonParameters);
+ double thisScore = CalculatePeptideScore(theScan.TheScan, matchedIons);
+ AddPeptideCandidateToPsm(scan, thisScore, peptide, matchedIons);
+ }
+ }
+ // Protein (FASTA) peptide source: digest each protein and search each peptide.
Action processProteinRange = (start, end) =>
{
List peptideTheorProducts = new();
- HashSet matchedFragmentIons = new();
- Tolerance productTolerance = CommonParameters.ProductMassTolerance;
+
for (int i = start; i < end; i++)
{
- // Stop loop if canceled
- if (GlobalVariables.StopLoops) { return; }
+ if (GlobalVariables.StopLoops) return;
- // digest each protein into peptides and search for each peptide in all spectra within precursor mass tolerance
+ // digest each protein into peptides and search each peptide in all spectra within precursor mass tolerance
foreach (var specificBioPolymer in Proteins[i].Digest(CommonParameters.DigestionParams, FixedModifications, VariableModifications))
- {
- Interlocked.Increment(ref peptideCounter);
-
- peptideTheorProducts.Clear();
- specificBioPolymer.Fragment(CommonParameters.DissociationType, CommonParameters.DigestionParams.FragmentationTerminus, peptideTheorProducts, CommonParameters.FragmentationParameters);
-
- // score each scan that has an acceptable precursor mass
- foreach (ScanWithIndexAndNotchInfo scan in GetAcceptableScans(specificBioPolymer.MonoisotopicMass, SearchMode))
- {
- matchedFragmentIons.Clear();
- Ms2ScanWithSpecificMass theScan = ArrayOfSortedMS2Scans[scan.ScanIndex];
- int precursorCharge = theScan.PrecursorCharge;
-
- // Match Fragment Ions
- foreach (var product in peptideTheorProducts)
- {
- // unknown fragment mass; this only happens rarely for sequences with unknown amino acids
- if (double.IsNaN(product.NeutralMass))
- {
- continue;
- }
-
- // get the closest peak in the spectrum to the theoretical peak
- var closestExperimentalMass = theScan.GetClosestExperimentalIsotopicEnvelope(product.NeutralMass);
-
- // is the mass error acceptable?
- if (closestExperimentalMass != null
- && productTolerance.Within(closestExperimentalMass.MonoisotopicMass, product.NeutralMass)
- && Math.Abs(closestExperimentalMass.Charge) <= Math.Abs(precursorCharge))//TODO apply this filter before picking the envelope
- {
- matchedFragmentIons.Add(new MatchedFragmentIon(product, closestExperimentalMass.MonoisotopicMass.ToMz(closestExperimentalMass.Charge),
- closestExperimentalMass.Peaks.First().intensity, closestExperimentalMass.Charge));
- }
- }
-
- if (matchedFragmentIons.Count < CommonParameters.ScoreCutoff)
- continue;
-
- // Score the peptide-spectrum match
- double tic = theScan.TotalIonCurrent;
- double score = 0;
- foreach (var ion in matchedFragmentIons)
- {
- if (ion.NeutralTheoreticalProduct.ProductType != ProductType.D)
- {
- score += 1 + ion.Intensity / tic;
- }
- }
-
- var matchedIons = matchedFragmentIons.ToList(); // materialize before passing to another thread
- AddPeptideCandidateToPsm(scan, score, specificBioPolymer, matchedIons);
- }
- }
+ ProcessOnePeptide(specificBioPolymer, peptideTheorProducts);
// report search progress (proteins searched so far out of total proteins in database)
proteinsSearched++;
var percentProgress = (int)((proteinsSearched / Proteins.Count) * 100);
-
if (percentProgress > oldPercentProgress)
{
oldPercentProgress = percentProgress;
@@ -139,19 +119,35 @@ protected override MetaMorpheusEngineResults RunSpecific()
}
};
+ // Library (.msl) peptide source: iterate precomputed peptides directly, no digestion.
+ Action processPeptideRange = (start, end) =>
+ {
+ List peptideTheorProducts = new();
+
+ for (int i = start; i < end; i++)
+ {
+ if (GlobalVariables.StopLoops) return;
+ var (peptide, fragments) = _precomputedPeptides[i];
+ ProcessOnePeptide(peptide, peptideTheorProducts, fragments);
+ }
+ };
+
+ int itemCount = usePrecomputedPeptides ? _precomputedPeptides.Count : Proteins.Count;
+ Action processRange = usePrecomputedPeptides ? processPeptideRange : processProteinRange;
+
if (_singleThreadMode)
{
- processProteinRange(0, Proteins.Count);
+ processRange(0, itemCount);
}
else
{
- var proteinPartioner = Partitioner.Create(0, Proteins.Count);
+ var partitioner = Partitioner.Create(0, itemCount);
Parallel.ForEach(
- proteinPartioner,
+ partitioner,
new ParallelOptions { MaxDegreeOfParallelism = CommonParameters.MaxThreadsToUsePerFile },
(range, loopState) =>
{
- processProteinRange(range.Item1, range.Item2);
+ processRange(range.Item1, range.Item2);
});
}
}
@@ -167,6 +163,10 @@ protected override MetaMorpheusEngineResults RunSpecific()
private void AddPeptideCandidateToPsm(ScanWithIndexAndNotchInfo scan, double thisScore, IBioPolymerWithSetMods peptide, List matchedIons)
{
+ // matched ion count below the cutoff is not a candidate (mirrors ClassicSearchEngine's ScoreCutoff gate)
+ if (matchedIons.Count < CommonParameters.ScoreCutoff)
+ return;
+
int scanIndex = scan.ScanIndex;
if (_singleThreadMode)
@@ -265,8 +265,8 @@ private SpectralMatch EnsureWritablePsm(int scanIndex, SpectralMatch existingPsm
if (!_copyOnWriteEnabled || existingPsm == null || UpdatedIndexes.ContainsKey(scanIndex))
{
return existingPsm;
- }
-
+ }
+
SpectralMatches[scanIndex] = ClonePsmForWrite(existingPsm);
return SpectralMatches[scanIndex];
}
diff --git a/MetaMorpheus/EngineLayer/ParallelSearch/TransientPepAnalysisEngine.cs b/MetaMorpheus/EngineLayer/ParallelSearch/TransientPepAnalysisEngine.cs
new file mode 100644
index 000000000..79e39af76
--- /dev/null
+++ b/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>>;
+
+namespace EngineLayer.ParallelSearch;
+
+///
+/// ParallelSearch-specific PEP engine. The generic 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 and
+/// ).
+///
+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(), _bgQByPepPsm = Array.Empty();
+ private double[] _bgPepAscPep = Array.Empty(), _bgQByPepPep = Array.Empty();
+
+ public TransientPepAnalysisEngine(List psms, string searchType,
+ List<(string fileName, CommonParameters fileSpecificParameters)> fileSpecificParameters,
+ string outputFolder, IRetentionTimePredictor? rtPredictor = null)
+ : base(psms, searchType, fileSpecificParameters, outputFolder, rtPredictor)
+ {
+ }
+
+ /// True once has produced a reusable model.
+ public bool HasTrainedModel => _trainedModel != null;
+
+ ///
+ /// Snapshots the (PEP ascending, PEP_QValue) curve over , which must already
+ /// carry a PEP in GetFdrInfo().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).
+ ///
+ internal static (double[] pepAsc, double[] qByPep) BuildPepQValueCurve(List 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(), Array.Empty());
+ 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);
+ }
+
+ ///
+ /// Maps a single PEP onto the background curve: returns the PEP_QValue of the background match whose PEP
+ /// is the smallest value >= (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.
+ ///
+ 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];
+ }
+
+ ///
+ /// 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 — 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 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.
+ ///
+ public bool TrainSingleModelAndAssignBasePep()
+ {
+ List 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;
+ }
+
+ ///
+ /// Assigns PEP to a NEW set of PSMs (a transient database's hits) using the model trained on the base
+ /// search by . 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.
+ ///
+ public void AssignPepFromTrainedModel(List 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 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);
+ }
+ }
+}
diff --git a/MetaMorpheus/GUI/GUI.csproj b/MetaMorpheus/GUI/GUI.csproj
index f1b354242..c9da4341c 100644
--- a/MetaMorpheus/GUI/GUI.csproj
+++ b/MetaMorpheus/GUI/GUI.csproj
@@ -63,7 +63,7 @@
-
+
diff --git a/MetaMorpheus/GuiFunctions/GuiFunctions.csproj b/MetaMorpheus/GuiFunctions/GuiFunctions.csproj
index 9c2aaa847..ad079f318 100644
--- a/MetaMorpheus/GuiFunctions/GuiFunctions.csproj
+++ b/MetaMorpheus/GuiFunctions/GuiFunctions.csproj
@@ -16,7 +16,7 @@
-
+
diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/Analysis/Collectors/BasicMetricCollector.cs b/MetaMorpheus/TaskLayer/ParallelSearch/Analysis/Collectors/BasicMetricCollector.cs
index d6951fd05..3d90653d2 100644
--- a/MetaMorpheus/TaskLayer/ParallelSearch/Analysis/Collectors/BasicMetricCollector.cs
+++ b/MetaMorpheus/TaskLayer/ParallelSearch/Analysis/Collectors/BasicMetricCollector.cs
@@ -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 GetOutputColumns()
@@ -35,6 +43,10 @@ public IEnumerable 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 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)
};
}
}
\ No newline at end of file
diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/Analysis/TransientDatabaseMetrics.cs b/MetaMorpheus/TaskLayer/ParallelSearch/Analysis/TransientDatabaseMetrics.cs
index f837d1cd5..0db0a0ffe 100644
--- a/MetaMorpheus/TaskLayer/ParallelSearch/Analysis/TransientDatabaseMetrics.cs
+++ b/MetaMorpheus/TaskLayer/ParallelSearch/Analysis/TransientDatabaseMetrics.cs
@@ -61,6 +61,12 @@ 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; }
@@ -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(BasicMetricCollector.TargetPeptidesAtQValueThreshold);
TargetPeptidesFromTransientDb = GetValue(BasicMetricCollector.TargetPeptidesFromTransientDb);
TargetPeptidesFromTransientDbAtQValueThreshold = GetValue(BasicMetricCollector.TargetPeptidesFromTransientDbAtQValueThreshold);
+ TargetPsmsFromTransientDbAtPepQ01 = GetValue(BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ01);
+ TargetPsmsFromTransientDbAtPepQ05 = GetValue(BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ05);
+ TargetPeptidesFromTransientDbAtPepQ01 = GetValue(BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ01);
+ TargetPeptidesFromTransientDbAtPepQ05 = GetValue(BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ05);
TargetProteinGroupsAtQValueThreshold = GetValue(ProteinGroupCollector.TargetProteinGroupsAtQValueThreshold);
TargetProteinGroupsFromTransientDb = GetValue(ProteinGroupCollector.TargetProteinGroupsFromTransientDb);
TargetProteinGroupsFromTransientDbAtQValueThreshold = GetValue(ProteinGroupCollector.TargetProteinGroupsFromTransientDbAtQValueThreshold);
diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/IO/StatisticalTestResultFile.cs b/MetaMorpheus/TaskLayer/ParallelSearch/IO/StatisticalTestResultFile.cs
index c86caf07d..3731d0bcd 100644
--- a/MetaMorpheus/TaskLayer/ParallelSearch/IO/StatisticalTestResultFile.cs
+++ b/MetaMorpheus/TaskLayer/ParallelSearch/IO/StatisticalTestResultFile.cs
@@ -32,6 +32,16 @@ public StatisticalTestResultFile(double alpha = 0.05) : base()
public override void LoadResults()
{
+ // Writer-only usage (parameterless ctor with in-memory Results) has no file to read.
+ // When the in-memory result list is empty the base lazy getter still calls LoadResults;
+ // guard so that does not throw on an empty/missing FilePath (e.g. a run where no
+ // statistical results were produced). Set via the setter to avoid getter recursion.
+ if (string.IsNullOrEmpty(FilePath) || !File.Exists(FilePath))
+ {
+ Results = new List();
+ return;
+ }
+
using var reader = new StreamReader(FilePath);
using var csv = new CsvReader(reader, new CsvConfiguration(CultureInfo.InvariantCulture)
{
diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/IO/TestSummaryResultFile.cs b/MetaMorpheus/TaskLayer/ParallelSearch/IO/TestSummaryResultFile.cs
index f45121d2e..53b7c5736 100644
--- a/MetaMorpheus/TaskLayer/ParallelSearch/IO/TestSummaryResultFile.cs
+++ b/MetaMorpheus/TaskLayer/ParallelSearch/IO/TestSummaryResultFile.cs
@@ -1,4 +1,5 @@
#nullable enable
+using System.Collections.Generic;
using System.Globalization;
using System.IO;
using System.Linq;
@@ -27,6 +28,14 @@ public TestSummaryResultFile() : base() { }
public override void LoadResults()
{
+ // Writer-only usage (parameterless ctor with in-memory Results) has no file to read;
+ // guard against the base lazy getter calling this on an empty/missing FilePath.
+ if (string.IsNullOrEmpty(FilePath) || !File.Exists(FilePath))
+ {
+ Results = new List();
+ return;
+ }
+
using var csv = new CsvReader(new StreamReader(FilePath), CsvConfiguration);
Results = csv.GetRecords().ToList();
}
diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchParameters.cs b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchParameters.cs
index 7bd49d408..08e18f9df 100644
--- a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchParameters.cs
+++ b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchParameters.cs
@@ -16,6 +16,17 @@ public class ParallelSearchParameters : SearchParameters
public bool CompressTransientSearchOutputs { get; set; } = false;
public string? DeNovoMappingDataFilePath { get; set; } = null;
+ ///
+ /// When true, the entries in are merged-index .msl files: a single file
+ /// (or a set of shards) holds many databases, each entry tagged "db|accession". Merged mode loads each
+ /// file once, runs the candidate filter once, and emits one searchable database per source db-group, so
+ /// parallelism is driven by the database count INSIDE the file rather than the file-list count. Replaces
+ /// the former MM_PARALLELSEARCH_MERGED environment variable so the behavior is on the
+ /// parameter/TOML/CLI surface and reproducible. Only takes effect when every transient database ends in
+ /// .msl.
+ ///
+ public bool UseMergedTransientLibrary { get; set; } = false;
+
#region Follow-Up Search Parameters
///
diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs
index d9ef752aa..9b47ff42a 100644
--- a/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs
+++ b/MetaMorpheus/TaskLayer/ParallelSearch/ParallelSearchTask.cs
@@ -9,12 +9,15 @@
using EngineLayer.Util;
using FlashLFQ;
using Nett;
+using Chromatography.RetentionTimePrediction;
using Omics;
+using Omics.Fragmentation;
using Omics.Modifications;
using Omics.SpectrumMatch;
using System;
using System.Collections.Concurrent;
using System.Collections.Generic;
+using System.Diagnostics;
using System.IO;
using System.Linq;
using System.Threading;
@@ -38,6 +41,12 @@ public class ParallelSearchTask : SearchTask
private int _dashboardInitialFinishedCount;
private int _dashboardCompletedThisRun;
+ // Phase timing (P0 of GPU plan): reveals where wall-clock actually goes so GPU effort
+ // targets the real hot spot. Inner-loop split is accumulated across the parallel
+ // per-database loop via Interlocked; coarse phases are timed in RunSpecific.
+ private long _searchEngineTicks; // time inside TransientClassicSearchEngine.Run
+ private long _postAnalysisTicks; // time inside PerformPostSearchAnalysis
+
private const string DashboardPhaseInitializing = "Initializing";
private const string DashboardPhaseSearching = "Searching";
private const string DashboardPhaseStatisticalAnalysis = "Statistical analysis";
@@ -106,6 +115,20 @@ public override SearchParameters SearchParameters
[TomlIgnore] public Ms2ScanWithSpecificMass[] AllSortedMs2Scans { get; private set; } = [];
[TomlIgnore] private SpectralMatch[] BaseSearchPsms = null!; // PSMs from base database search
+ // Calibration LEARNED from the base search (S1), applied to the .msl candidate pre-filter (S3):
+ // precursor tolerance, and the iRT→observed-RT regression + window. Null until learned (and stays
+ // null when no transient database is a .msl library). See LearnMslCalibration.
+ [TomlIgnore] private MslCandidateCalibration? _mslCalibration;
+ // Scan precursor masses (ascending) + their RTs, cached once (identical across all .msl databases).
+ [TomlIgnore] private double[]? _sortedScanMasses;
+ [TomlIgnore] private double[]? _scanRetentionTimes;
+
+ // 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;
+ // The PEP model's RT-feature predictor (Chronologer); lives as long as _pepEngine, disposed at task end.
+ [TomlIgnore] private Chromatography.RetentionTimePrediction.IRetentionTimePredictor? _pepRtPredictor;
+
// Optimization caches for FDR alignment and parsimony
[TomlIgnore] private readonly PsmSpectralMatchFdrAlignmentService _psmFdrAlignmentService = new();
@@ -156,35 +179,166 @@ protected override MyTaskResults RunSpecific(string outputFolder,
goto Finalization;
}
+ // Phase timing (P0 of GPU plan) — coarse wall-clock per phase.
+ var swInit = Stopwatch.StartNew();
+
// Initialize all necessary data structures including base search
Initialize(taskId, dbFilenameList, currentRawFileList, fileSettingsList, outputFolder);
InitializeCompletedDatabaseWriter();
+ swInit.Stop();
Status($"Starting search of {TotalDatabases} transient databases...", taskId);
ReportTaskDashboard(taskId, ParallelSearchDashboardUpdateKind.TaskStatus, DashboardPhaseSearching,
$"Searching {TotalDatabases} transient databases...");
+ // MERGED-INDEX mode: one .msl file holds many databases (entries tagged "db|accession"). The number
+ // of databases to search comes from inside the file, not the (count==1) file list, so parallelism
+ // must not be capped by the file count.
+ // MERGED-INDEX mode also supports SHARDING: several merged .msl files, each holding many
+ // databases (entries tagged "db|accession"). A single merged file exceeds the format's 2^31
+ // precursor cap at production scale, so the library is split into shards with each database
+ // kept WHOLE in exactly one shard — the per-database union across shards is therefore clean.
+ bool mergedMode = ParallelSearchParameters.UseMergedTransientLibrary
+ && ParallelSearchParameters.TransientDatabases.Count >= 1
+ && ParallelSearchParameters.TransientDatabases.All(d => d.FilePath.EndsWith(".msl", StringComparison.OrdinalIgnoreCase));
+
// Determine optimal thread allocation
int totalAvailableThreads = Environment.ProcessorCount;
- int databaseParallelism = Math.Min(ParallelSearchParameters.MaxSearchesInParallel,
- ParallelSearchParameters.TransientDatabases.Count);
+ int databaseParallelism = mergedMode
+ ? ParallelSearchParameters.MaxSearchesInParallel
+ : Math.Min(ParallelSearchParameters.MaxSearchesInParallel,
+ ParallelSearchParameters.TransientDatabases.Count);
int threadsPerDatabase = Math.Max(1, totalAvailableThreads / databaseParallelism);
CommonParameters.MaxThreadsToUsePerFile = threadsPerDatabase;
- // Loop through each transient database
- try
+ // S4: LOAD-AHEAD PRODUCER/CONSUMER. Database loading (now an I/O-bound .msl index query + candidate
+ // fragment fetch) runs ahead on producer threads and hands PREPARED databases to searcher threads
+ // through a bounded queue, so loading overlaps with the CPU-bound search instead of blocking it.
+ // The result WRITE channel was already pipelined; this pipelines the read side too. The bounded
+ // capacity caps how many loaded databases are held in memory at once.
+ var swTransientLoop = Stopwatch.StartNew();
+ // Abort signal: if a CONSUMER (search) faults, cancel so the producer stops blocking on Add into a
+ // bounded queue that is no longer being drained, instead of being abandoned until queue disposal.
+ using (var loadAbort = new System.Threading.CancellationTokenSource())
+ using (var loadedQueue = new System.Collections.Concurrent.BlockingCollection(
+ boundedCapacity: Math.Max(2, databaseParallelism)))
{
- Parallel.ForEach(ParallelSearchParameters.TransientDatabases,
- new ParallelOptions { MaxDegreeOfParallelism = databaseParallelism },
- transientDbPath =>
+ // MERGED-INDEX mode (mergedMode computed above): a single .msl holds every database (each entry
+ // tagged "db|accession"). Load it ONCE and run the candidate filter ONCE, then emit one prepared
+ // database per source db-group. Each is searched INDEPENDENTLY by the consumer below (own
+ // base-PSM copy) — databases never compete; this only collapses 1000s of file opens into one
+ // shared in-memory index.
+
+ // Producers: load databases concurrently, blocking on the bounded queue when it is full.
+ var producer = Task.Run(() =>
+ {
+ try
{
- ProcessTransientDatabase(transientDbPath, outputFolder, taskId);
- });
- }
- finally
- {
- CompleteCompletedDatabaseWriter();
+ if (mergedMode)
+ {
+ // One or more merged shards. Load each in turn, run the candidate filter once per
+ // shard, and emit one prepared database per source db-group. Because every database
+ // is wholly contained in a single shard, db-groups never collide across shards.
+ var priors = BuildCandidatePriors();
+ // Parallelize the per-shard load+filter. It was serial (~72s/shard: LoadIndexOnly +
+ // single-threaded candidate filter over ~18M entries + GetEntry fetch), i.e. ~3.2h at
+ // 160 shards and the wall bottleneck (the consumer search parallelizes fine). Each shard
+ // has its OWN index/file handle, so K shards load concurrently; GetEntry within a shard
+ // stays single-threaded (one FileStream). Bounded by MM_PARALLELSEARCH_PRODUCERS (default
+ // 4) to not oversubscribe the consumer search or the disk; the loadedQueue still caps RAM.
+ int producerDop = int.TryParse(Environment.GetEnvironmentVariable("MM_PARALLELSEARCH_PRODUCERS"), out var pdv) && pdv > 0
+ ? pdv : 4;
+ producerDop = Math.Max(1, Math.Min(producerDop, ParallelSearchParameters.TransientDatabases.Count));
+ Parallel.ForEach(ParallelSearchParameters.TransientDatabases,
+ new ParallelOptions { MaxDegreeOfParallelism = producerDop },
+ mergedDb =>
+ {
+ if (GlobalVariables.StopLoops) return;
+ string shardName = Path.GetFileNameWithoutExtension(mergedDb.FilePath);
+ Status($"Loading merged transient index {shardName}...", taskId);
+ var grouped = MslPeptideReader.ReadCandidatesGroupedByDatabase(
+ mergedDb.FilePath, priors);
+ Status($"Merged index {shardName}: {grouped.Count} databases with candidates.", taskId);
+ // Emit this shard's prepared db-groups serially (light work; the outer loop is
+ // already K-way parallel, so an inner Parallel.ForEach would just oversubscribe).
+ foreach (var kvp in grouped)
+ {
+ if (GlobalVariables.StopLoops) return;
+ var loaded = BuildLoadedFromCandidates(kvp.Key, kvp.Value, outputFolder, taskId);
+ if (loaded != null) loadedQueue.Add(loaded, loadAbort.Token);
+ }
+ });
+ }
+ else
+ {
+ Parallel.ForEach(ParallelSearchParameters.TransientDatabases,
+ new ParallelOptions { MaxDegreeOfParallelism = databaseParallelism },
+ transientDbPath =>
+ {
+ if (GlobalVariables.StopLoops) return;
+ LoadedTransientDatabase? loaded;
+ try
+ {
+ loaded = LoadTransientDatabaseForPipeline(transientDbPath, outputFolder, taskId);
+ }
+ catch (Exception ex)
+ {
+ WriteTransientProcessError(transientDbPath, outputFolder, ex);
+ throw;
+ }
+ if (loaded != null)
+ loadedQueue.Add(loaded, loadAbort.Token);
+ });
+ }
+ }
+ finally
+ {
+ loadedQueue.CompleteAdding();
+ }
+ });
+
+ try
+ {
+ // Consumers: search each loaded database as it becomes available.
+ Parallel.ForEach(loadedQueue.GetConsumingEnumerable(),
+ new ParallelOptions { MaxDegreeOfParallelism = databaseParallelism },
+ loaded =>
+ {
+ if (GlobalVariables.StopLoops) return;
+ try
+ {
+ SearchLoadedTransientDatabase(loaded, taskId);
+ }
+ catch (Exception ex)
+ {
+ WriteTransientProcessError(loaded.TransientDb, outputFolder, ex);
+ throw;
+ }
+ });
+ }
+ catch
+ {
+ // A searcher faulted and Parallel.ForEach stopped draining the queue. Signal the producer so a
+ // load thread blocked on the bounded Add unwinds promptly instead of being abandoned; then rethrow.
+ loadAbort.Cancel();
+ throw;
+ }
+ finally
+ {
+ CompleteCompletedDatabaseWriter();
+ // ALWAYS observe the producer task: on the happy path this surfaces a load exception; on the
+ // failure path it prevents an unobserved task / abandoned loader. When WE triggered the abort
+ // above, the producer's resulting cancellation is a consequence, not the root cause — swallow it.
+ try { producer.GetAwaiter().GetResult(); }
+ catch (Exception) when (loadAbort.IsCancellationRequested) { }
+ }
}
+ swTransientLoop.Stop();
+ LogPhaseTimingBreakdown(taskId, outputFolder, swInit.Elapsed, swTransientLoop.Elapsed, databaseParallelism);
+
+ // PEP RT predictor (Chronologer/TorchSharp) is finished once the transient loop is done; release it.
+ _pepRtPredictor?.Dispose();
+ _pepRtPredictor = null;
Finalization:
@@ -237,7 +391,16 @@ protected override MyTaskResults RunSpecific(string outputFolder,
Status("Writing Final Results...", taskId);
ReportTaskDashboard(taskId, ParallelSearchDashboardUpdateKind.TaskStatus, DashboardPhaseWritingFinalResults,
"Writing final results...");
- WriteFinalOutputs(outputFolder, taskId, currentRawFileList.Count);
+ try
+ {
+ WriteFinalOutputs(outputFolder, taskId, currentRawFileList.Count);
+ }
+ catch (Exception ex)
+ {
+ // Dump the full stack to a labeled file — the task runner only surfaces the message.
+ try { File.WriteAllText(Path.Combine(outputFolder, "WriteFinalOutputs_ERROR.txt"), ex.ToString()); } catch { }
+ throw;
+ }
ReportTaskDashboard(taskId, ParallelSearchDashboardUpdateKind.TaskCompleted, DashboardPhaseCompleted,
"Many search task complete!");
@@ -246,6 +409,44 @@ protected override MyTaskResults RunSpecific(string outputFolder,
return MyTaskResults;
}
+ ///
+ /// Logs where wall-clock time went so GPU effort targets the real hot spot.
+ /// The per-database search and post-analysis times are summed across all parallel
+ /// workers (CPU-seconds); dividing by the database parallelism approximates the
+ /// wall-clock each contributed to the transient loop.
+ ///
+ private void LogPhaseTimingBreakdown(string taskId, string outputFolder, TimeSpan init, TimeSpan transientLoop, int databaseParallelism)
+ {
+ double searchCpuSec = _searchEngineTicks / (double)Stopwatch.Frequency;
+ double postCpuSec = _postAnalysisTicks / (double)Stopwatch.Frequency;
+ int par = Math.Max(1, databaseParallelism);
+
+ string summary =
+ "Phase timing — " +
+ $"Initialize (load + base search): {init.TotalSeconds:F1}s | " +
+ $"Transient loop wall-clock: {transientLoop.TotalSeconds:F1}s | " +
+ $"search engine: {searchCpuSec:F1} CPU-s (~{searchCpuSec / par:F1}s wall) | " +
+ $"post-search analysis: {postCpuSec:F1} CPU-s (~{postCpuSec / par:F1}s wall) | " +
+ $"db parallelism: {par}.";
+
+ Status(summary, taskId);
+
+ // Also persist to a flushed file so the profile survives any later-stage failure
+ // (the breakdown is the whole point of the run; don't let finalization bugs eat it).
+ try
+ {
+ string detail =
+ summary + Environment.NewLine + Environment.NewLine +
+ "Breakdown (transient loop only; Initialize is a one-time cost):" + Environment.NewLine +
+ $" search engine (TransientClassicSearchEngine): {searchCpuSec,10:F1} CPU-s" + Environment.NewLine +
+ $" post-search analysis (stats/collectors/FDR): {postCpuSec,10:F1} CPU-s" + Environment.NewLine +
+ $" search : post ratio = {searchCpuSec / Math.Max(1e-9, postCpuSec):F2} : 1" + Environment.NewLine +
+ $" db parallelism = {par}" + Environment.NewLine;
+ File.WriteAllText(Path.Combine(outputFolder, "PhaseTiming.txt"), detail);
+ }
+ catch { /* timing file is best-effort */ }
+ }
+
#region Initialization
private void Initialize(string taskId, List dbFilenameList,
@@ -317,6 +518,25 @@ private void Initialize(string taskId, List dbFilenameList,
_peptideFdrAlignmentService.BuildBaselineCache(baselinePsms);
BuildBaselinePeptideToScanIndexLookup();
+ // S1: when any transient database is a .msl library, learn the precursor tolerance and the
+ // RT (iRT→observed) calibration from the confident base-search PSMs. These become the .msl
+ // candidate pre-filter priors (S3). Cheap relative to the search and done once.
+ if (ParallelSearchParameters.TransientDatabases.Any(
+ d => d.FilePath.EndsWith(".msl", StringComparison.OrdinalIgnoreCase)))
+ {
+ _mslCalibration = LearnMslCalibration(baselinePsms, taskId);
+ // Cache the (mass-sorted) scan mass + RT arrays ONCE — they are identical for every .msl
+ // database, so rebuilding them per database (1000s of times) would be pure waste.
+ _sortedScanMasses = Array.ConvertAll(AllSortedMs2Scans, s => s.PrecursorMass);
+ _scanRetentionTimes = Array.ConvertAll(AllSortedMs2Scans, s => s.RetentionTime);
+ }
+
+ // PEP: train ONE model on the (large, high-quality) base human PSMs and reuse it to assign a PEP
+ // to every transient database's PSMs — those databases are far too small to train their own model.
+ // Independent of the .msl path, so FASTA transient searches get PEP too. Uses Chronologer for the
+ // RT-residual feature. Best-effort: any failure leaves transient PEP unset.
+ TrainBasePepModel(baselinePsms, outputFolder, taskId);
+
if (SearchParameters.DoParsimony)
{
Status("Preparing baseline parsimony cache...", taskId);
@@ -490,32 +710,99 @@ private TransientDatabaseResultsManager CreateResultsManager(string outputFolder
#region Search
- private void ProcessTransientDatabase(DbForTask transientDb, string outputFolder, string taskId)
+ ///
+ /// A transient database after the LOAD stage (producer): proteins/peptides ready, fragments fetched.
+ /// Carried through the bounded channel to a searcher (consumer) so loading overlaps with searching.
+ ///
+ private sealed class LoadedTransientDatabase
+ {
+ public DbForTask TransientDb = null!;
+ public string DbName = null!;
+ public string DbOutputFolder = null!;
+ public List NestedIds = null!;
+ public List TransientProteins = null!;
+ public List<(IBioPolymerWithSetMods Peptide, List Fragments)>? PrecomputedPeptides;
+ public HashSet TransientProteinAccessions = null!;
+ }
+
+ /// Builds the .msl candidate pre-filter priors (scan masses + RTs, learned precursor/RT
+ /// calibration). The scan arrays are cached once (identical for every database).
+ private MslPeptideReader.CandidatePriors BuildCandidatePriors()
+ {
+ var sortedScanMasses = _sortedScanMasses ?? Array.ConvertAll(AllSortedMs2Scans, s => s.PrecursorMass);
+ var scanRetentionTimes = _scanRetentionTimes ?? Array.ConvertAll(AllSortedMs2Scans, s => s.RetentionTime);
+ var cal = _mslCalibration ?? new MslCandidateCalibration(
+ CommonParameters.PrecursorMassTolerance.Value, 1, 0, double.PositiveInfinity);
+ return new MslPeptideReader.CandidatePriors(
+ sortedScanMasses, scanRetentionTimes,
+ precursorTolPpm: cal.PrecursorTolPpm, rtSlope: cal.RtSlope,
+ rtIntercept: cal.RtIntercept, rtWindowMin: cal.RtWindowMin);
+ }
+
+ ///
+ /// Builds a prepared database from a merged-index db-group's candidate peptides (no file load). The
+ /// returned database is searched INDEPENDENTLY by the consumer, identical to the per-file path.
+ ///
+ private LoadedTransientDatabase? BuildLoadedFromCandidates(string dbName,
+ List<(IBioPolymerWithSetMods Peptide, List Fragments)> candidates, string outputFolder, string taskId)
{
if (GlobalVariables.StopLoops)
- return;
+ return null;
+
+ List nestedIds = [taskId, dbName];
+ bool shouldProcess = !_resultsManager!.HasCachedResults(dbName) || ParallelSearchParameters.OverwriteTransientSearchOutputs;
+ if (!shouldProcess)
+ {
+ ReportProgress(new(100, $"Skipping {dbName} - results already exist in cache", nestedIds));
+ UpdateProgress(TotalDatabases, taskId);
+ return null;
+ }
+
+ var transientProteins = candidates.Select(p => p.Peptide.Parent).Distinct().ToList();
+ return new LoadedTransientDatabase
+ {
+ // Synthetic per-group identity (name only); no per-database file in merged mode.
+ TransientDb = new DbForTask(dbName, false),
+ DbName = dbName,
+ DbOutputFolder = Path.Combine(outputFolder, dbName),
+ NestedIds = nestedIds,
+ TransientProteins = transientProteins,
+ PrecomputedPeptides = candidates,
+ TransientProteinAccessions = new HashSet(transientProteins.Select(p => p.Accession)),
+ };
+ }
+
+ ///
+ /// PRODUCER stage: skip/overwrite handling + load the transient database (FASTA digest set, or the
+ /// .msl mass+RT-filtered candidate peptides with their fragments). Returns null when the database is
+ /// skipped (cached) or the run is stopping. I/O-bound; runs ahead of the searchers.
+ ///
+ private LoadedTransientDatabase? LoadTransientDatabaseForPipeline(DbForTask transientDb, string outputFolder, string taskId)
+ {
+ if (GlobalVariables.StopLoops)
+ return null;
- string dbName = Path.GetFileNameWithoutExtension(transientDb.FilePath);
- string dbOutputFolder = Path.Combine(outputFolder, dbName);
- List nestedIds = [taskId, dbName];
+ string dbName = Path.GetFileNameWithoutExtension(transientDb.FilePath);
+ string dbOutputFolder = Path.Combine(outputFolder, dbName);
+ List nestedIds = [taskId, dbName];
Status($"Processing {dbName}...", nestedIds);
- // Check if we should skip or overwrite this database
- bool shouldProcess = !_resultsManager!.HasCachedResults(dbName) ||
- ParallelSearchParameters.OverwriteTransientSearchOutputs;
+ // Check if we should skip or overwrite this database
+ bool shouldProcess = !_resultsManager!.HasCachedResults(dbName) ||
+ ParallelSearchParameters.OverwriteTransientSearchOutputs;
if (!shouldProcess)
{
ReportProgress(new(100, $"Skipping {dbName} - results already exist in cache", nestedIds));
UpdateProgress(TotalDatabases, taskId);
- return;
+ return null;
}
ReportDatabaseDashboard(taskId, ParallelSearchDashboardUpdateKind.DatabaseStarted, dbName,
$"Processing {dbName}...", DashboardDatabaseProcessingProgress);
- // Handle overwrite scenario
+ // Handle overwrite scenario
if (ParallelSearchParameters.OverwriteTransientSearchOutputs && _resultsManager.HasCachedResults(dbName))
{
Status($"Overwriting existing results for {dbName}...", nestedIds);
@@ -527,55 +814,100 @@ private void ProcessTransientDatabase(DbForTask transientDb, string outputFolder
}
}
- if (!Directory.Exists(dbOutputFolder))
- Directory.CreateDirectory(dbOutputFolder);
+ // NOTE: the output folder is created lazily by the writer, and ONLY for databases that produced
+ // transient PSMs — at 1000s of databases the vast majority match nothing, so creating a folder +
+ // header-only files for each was a large fraction of the runtime. See WriteCompletedDatabaseOutputsAsync.
Status($"Loading transient database {dbName}...", nestedIds);
ReportDatabaseDashboard(taskId, ParallelSearchDashboardUpdateKind.DatabaseProgress, dbName,
$"Loading transient database {dbName}...", DashboardDatabaseLoadingProgress);
- // Load transient database
- var transientProteins = LoadTransientDatabase(transientDb, nestedIds, taskId);
+ // Load transient database. FASTA/XML -> proteins (digested during search); a .msl spectral
+ // library -> precomputed peptides paired with their stored (float) fragments (the search
+ // iterates these and matches the stored fragments directly, skipping digestion+fragmentation).
+ List transientProteins;
+ List<(IBioPolymerWithSetMods Peptide, List Fragments)>? precomputedPeptides = null;
+ bool isMslLibrary = transientDb.FilePath.EndsWith(".msl", StringComparison.OrdinalIgnoreCase);
+ if (isMslLibrary)
+ {
+ // S3: index-only candidate pre-filter on precursor mass AND Chronologer-iRT (learned from the
+ // base search, S1). See BuildCandidatePriors.
+ precomputedPeptides = MslPeptideReader.ReadPeptides(transientDb.FilePath, dbName, BuildCandidatePriors());
+ // shared parent proteins (one per accession) back the accession filter and counts
+ transientProteins = precomputedPeptides.Select(p => p.Peptide.Parent).Distinct().ToList();
+ }
+ else
+ {
+ transientProteins = LoadTransientDatabase(transientDb, nestedIds, taskId);
+ }
- if (GlobalVariables.StopLoops)
- {
- return;
- }
+ if (GlobalVariables.StopLoops)
+ return null;
- // Create HashSet of transient bioPolymer accessions for later filtering
- var transientProteinAccessions = new HashSet(
- transientProteins.Select(p => p.Accession));
+ return new LoadedTransientDatabase
+ {
+ TransientDb = transientDb,
+ DbName = dbName,
+ DbOutputFolder = dbOutputFolder,
+ NestedIds = nestedIds,
+ TransientProteins = transientProteins,
+ PrecomputedPeptides = precomputedPeptides,
+ // Create HashSet of transient bioPolymer accessions for later filtering
+ TransientProteinAccessions = new HashSet(transientProteins.Select(p => p.Accession)),
+ };
+ }
+
+ ///
+ /// CONSUMER stage: search the loaded database against the shared spectra and run post-analysis, then
+ /// hand the results to the (already-pipelined) write channel. CPU-bound; overlaps with loaders.
+ ///
+ private void SearchLoadedTransientDatabase(LoadedTransientDatabase loaded, string taskId)
+ {
+ if (GlobalVariables.StopLoops)
+ return;
+
+ DbForTask transientDb = loaded.TransientDb;
+ string dbName = loaded.DbName;
+ string dbOutputFolder = loaded.DbOutputFolder;
+ List nestedIds = loaded.NestedIds;
+ List transientProteins = loaded.TransientProteins;
+ HashSet transientProteinAccessions = loaded.TransientProteinAccessions;
Status($"Searching {dbName} ({transientProteins.Count} transient proteins)...", nestedIds);
ReportDatabaseDashboard(taskId, ParallelSearchDashboardUpdateKind.DatabaseProgress, dbName,
$"Searching {dbName} ({transientProteins.Count} transient proteins)...", DashboardDatabaseSearchStartProgress,
DashboardDatabaseSearchStartProgress, DashboardDatabaseSearchEndProgress);
- // Reuse baseline PSMs with copy-on-write in peptide/proteoform mode.
- SpectralMatch[] psmArray = BaseSearchPsms.ToArray();
- PerformSearch(transientProteins, psmArray, nestedIds, out HashSet updatedPsmIndexes, out int transientPeptideCount, useCopyOnWrite: true);
+ // Reuse baseline PSMs with copy-on-write in peptide/proteoform mode.
+ SpectralMatch[] psmArray = BaseSearchPsms.ToArray();
+ long searchStart = Stopwatch.GetTimestamp();
+ PerformSearch(transientProteins, psmArray, nestedIds, out HashSet updatedPsmIndexes, out int transientPeptideCount, useCopyOnWrite: true, precomputedPeptides: loaded.PrecomputedPeptides);
+ Interlocked.Add(ref _searchEngineTicks, Stopwatch.GetTimestamp() - searchStart);
Status($"Performing post-search analysis for {dbName}...", nestedIds);
ReportDatabaseDashboard(taskId, ParallelSearchDashboardUpdateKind.DatabaseProgress, dbName,
$"Performing post-search analysis for {dbName}...", DashboardDatabasePostSearchProgress);
- int totalProteins = BaseBioPolymers.Count + transientProteins.Count;
-
- // Process database through unified manager (handles analysis + statistical caching)
- var (analysisContext, dbResults) = PerformPostSearchAnalysis(
- psmArray,
- dbOutputFolder,
- nestedIds,
- dbName,
- totalProteins,
- transientProteinAccessions,
- transientPeptideCount,
- transientProteins,
- transientDb,
- updatedPsmIndexes
- ).GetAwaiter().GetResult();
-
- _completedDatabaseWriteChannel!.Writer.WriteAsync((analysisContext, dbResults)).AsTask().GetAwaiter().GetResult();
+ int totalProteins = BaseBioPolymers.Count + transientProteins.Count;
+
+ // Process database through unified manager (handles analysis + statistical caching)
+ long postAnalysisStart = Stopwatch.GetTimestamp();
+ var (analysisContext, dbResults) = PerformPostSearchAnalysis(
+ psmArray,
+ dbOutputFolder,
+ nestedIds,
+ dbName,
+ totalProteins,
+ transientProteinAccessions,
+ transientPeptideCount,
+ transientProteins,
+ transientDb,
+ updatedPsmIndexes
+ ).GetAwaiter().GetResult();
+ Interlocked.Add(ref _postAnalysisTicks, Stopwatch.GetTimestamp() - postAnalysisStart);
+
+ // (PEP is now assigned inside PerformPostSearchAnalysis, before the metric collectors run.)
+ _completedDatabaseWriteChannel!.Writer.WriteAsync((analysisContext, dbResults)).AsTask().GetAwaiter().GetResult();
// Cleanup transient proteins to free memory
transientProteins.Clear();
@@ -583,21 +915,173 @@ private void ProcessTransientDatabase(DbForTask transientDb, string outputFolder
ReportProgress(new(100, $"Finished analysis for {dbName}; queued output writing", nestedIds));
}
+ /// Writes a per-database <db>_PROCESS_ERROR.txt diagnostic; never throws.
+ private static void WriteTransientProcessError(DbForTask transientDb, string outputFolder, Exception ex)
+ {
+ try
+ {
+ string dbName = Path.GetFileNameWithoutExtension(transientDb.FilePath);
+ File.WriteAllText(Path.Combine(outputFolder, dbName + "_PROCESS_ERROR.txt"), ex.ToString());
+ }
+ catch { }
+ }
+
+ /// Calibration learned from the base search and applied to the .msl candidate pre-filter.
+ private readonly struct MslCandidateCalibration
+ {
+ public readonly double PrecursorTolPpm; // symmetric precursor tolerance (covers offset+spread)
+ public readonly double RtSlope, RtIntercept; // observedRT = RtSlope*iRT + RtIntercept
+ public readonly double RtWindowMin; // +/- observed-RT window (k * residual SD); +inf = no RT filter
+ public MslCandidateCalibration(double tolPpm, double slope, double intercept, double rtWindowMin)
+ { PrecursorTolPpm = tolPpm; RtSlope = slope; RtIntercept = intercept; RtWindowMin = rtWindowMin; }
+ }
+
+ ///
+ /// S1 — learn the .msl candidate pre-filter priors from confident base-search PSMs:
+ /// • precursor tolerance = |mean| + 3·SD of the precursor mass error (clamped to the configured tol),
+ /// • RT calibration = OLS of observed RT on Chronologer iRT, with a ±2·residualSD window.
+ /// On any failure (too few PSMs, predictor/model unavailable) returns a SAFE fallback: the configured
+ /// precursor tolerance and NO RT filter (RtWindowMin = +inf) so nothing is lost.
+ ///
+ ///
+ /// Trains ONE PEP model on the base (human) search PSMs and keeps it for assigning PEP to each transient
+ /// database's hits (the databases are far too small to train their own). Uses Chronologer for the
+ /// RT-residual feature. Best-effort: any failure simply leaves the transient PEP unassigned.
+ ///
+ private void TrainBasePepModel(List baselinePsms, string outputFolder, string taskId)
+ {
+ try
+ {
+ var trainingPsms = baselinePsms.Where(p => p != null).ToList();
+ if (trainingPsms.Count < 100)
+ {
+ Status($"PEP: only {trainingPsms.Count} base PSMs — skipping PEP model.", taskId);
+ return;
+ }
+ _pepRtPredictor = RetentionTimePredictorFactory.Create(PredictorType.Chronologer);
+ var engine = new TransientPepAnalysisEngine(trainingPsms, "standard", FileSpecificParameters, outputFolder, _pepRtPredictor);
+ if (engine.TrainSingleModelAndAssignBasePep())
+ {
+ _pepEngine = engine;
+ Status("PEP: trained model on base search; assigning PEP to transient PSMs.", taskId);
+ }
+ else
+ {
+ Status("PEP: base PSMs lacked target/decoy training examples — PEP disabled.", taskId);
+ _pepRtPredictor.Dispose();
+ _pepRtPredictor = null;
+ }
+ }
+ catch (Exception ex)
+ {
+ Status($"PEP training failed ({ex.GetType().Name}: {ex.Message}); PEP disabled.", taskId);
+ _pepRtPredictor?.Dispose();
+ _pepRtPredictor = null;
+ }
+ }
+
+ private MslCandidateCalibration LearnMslCalibration(List baselinePsms, string taskId)
+ {
+ double configuredTolPpm = CommonParameters.PrecursorMassTolerance.Value;
+ var fallback = new MslCandidateCalibration(configuredTolPpm, 1, 0, double.PositiveInfinity);
+ try
+ {
+ // Confident, unambiguous target PSMs.
+ var conf = baselinePsms.Where(p => p != null && !p.IsDecoy
+ && p.PsmFdrInfo != null && p.PsmFdrInfo.QValue < 0.01
+ && p.BestMatchingBioPolymersWithSetMods.Count() == 1)
+ .ToList();
+ if (conf.Count < 50)
+ {
+ Status($"S1: only {conf.Count} confident base PSMs — using safe fallback (no RT filter).", taskId);
+ return fallback;
+ }
+
+ // Precursor mass error (ppm) and the (peptide, observedRT) pairs for the RT regression.
+ var ppmErrors = new List(conf.Count);
+ var obsByPeptide = new Dictionary(ReferenceEqualityComparer.Instance);
+ var peptides = new List(conf.Count);
+ foreach (var psm in conf)
+ {
+ var pep = psm.BestMatchingBioPolymersWithSetMods.First().SpecificBioPolymer;
+ double mass = pep.MonoisotopicMass;
+ if (mass <= 0) continue;
+ ppmErrors.Add((psm.ScanPrecursorMass - mass) / mass * 1e6);
+ if (pep is IRetentionPredictable rp && !obsByPeptide.ContainsKey(rp))
+ {
+ obsByPeptide[rp] = psm.ScanRetentionTime;
+ peptides.Add(rp);
+ }
+ }
+
+ double pMean = ppmErrors.Average();
+ double pSd = Math.Sqrt(ppmErrors.Sum(e => (e - pMean) * (e - pMean)) / ppmErrors.Count);
+ // Trust the calibration: the real precursor accuracy (|mean|+3σ of the confident base PSMs) IS
+ // the tolerance, regardless of the looser configured/standard value. This is applied to BOTH the
+ // candidate pre-filter AND the transient search's acceptor (see PerformSearch), so they are
+ // consistent — peptides outside it are rejected as loose-tolerance false positives, not lost.
+ double tolPpm = Math.Max(2.0, Math.Abs(pMean) + 3.0 * pSd);
+
+ // Chronologer iRT for the confident peptides, then OLS observedRT vs iRT.
+ using var predictor = RetentionTimePredictorFactory.Create(PredictorType.Chronologer);
+ var predRt = new List(peptides.Count);
+ var obsRt = new List(peptides.Count);
+ foreach (var r in predictor.PredictRetentionTimeEquivalents(peptides,
+ maxThreads: Math.Max(1, Environment.ProcessorCount - 2)))
+ {
+ if (r.PredictedValue is null) continue;
+ if (obsByPeptide.TryGetValue(r.Peptide, out double rt)) { predRt.Add(r.PredictedValue.Value); obsRt.Add(rt); }
+ }
+ if (predRt.Count < 50)
+ {
+ Status($"S1: only {predRt.Count} RT predictions — precursor tol {tolPpm:F1} ppm, no RT filter.", taskId);
+ return new MslCandidateCalibration(tolPpm, 1, 0, double.PositiveInfinity);
+ }
+
+ int n = predRt.Count;
+ double mx = predRt.Average(), my = obsRt.Average();
+ double sxy = 0, sxx = 0;
+ for (int i = 0; i < n; i++) { double dx = predRt[i] - mx; sxy += dx * (obsRt[i] - my); sxx += dx * dx; }
+ double slope = sxx > 0 ? sxy / sxx : 1.0;
+ double intercept = my - slope * mx;
+ double residSs = 0;
+ for (int i = 0; i < n; i++) { double e = obsRt[i] - (slope * predRt[i] + intercept); residSs += e * e; }
+ double residSd = Math.Sqrt(residSs / n);
+ double rtWindow = Math.Max(1.0, 2.0 * residSd); // ±2σ; floor 1 min
+
+ Status($"S1 learned: precursor ±{tolPpm:F1} ppm (μ={pMean:F2},σ={pSd:F2}); " +
+ $"RT obs={slope:F3}·iRT+{intercept:F2} residSD={residSd:F2}min → window ±{rtWindow:F1}min (n={n}).", taskId);
+ return new MslCandidateCalibration(tolPpm, slope, intercept, rtWindow);
+ }
+ catch (Exception ex)
+ {
+ Status($"S1 calibration failed ({ex.GetType().Name}: {ex.Message}); safe fallback (no RT filter).", taskId);
+ return fallback;
+ }
+ }
+
///
/// Populates and returns the spectral match array using classic search engine
///
- private void PerformSearch(List proteinsToSearch, SpectralMatch[] spectralMatchArray, List nestedIds, out HashSet updatedPsmIndexes, out int peptidesSearched, bool useCopyOnWrite = false)
+ private void PerformSearch(List proteinsToSearch, SpectralMatch[] spectralMatchArray, List nestedIds, out HashSet updatedPsmIndexes, out int peptidesSearched, bool useCopyOnWrite = false, List<(IBioPolymerWithSetMods Peptide, List Fragments)> precomputedPeptides = null)
{
+ // For a .msl search, use the LEARNED precursor tolerance (S1) so the search engine and the
+ // candidate pre-filter agree — the calibration is trusted over the looser configured tolerance.
+ MzLibUtil.Tolerance precursorTolerance = CommonParameters.PrecursorMassTolerance;
+ if (precomputedPeptides != null && _mslCalibration.HasValue)
+ precursorTolerance = new MzLibUtil.PpmTolerance(_mslCalibration.Value.PrecursorTolPpm);
+
var massDiffAcceptor = GetMassDiffAcceptor(
- CommonParameters.PrecursorMassTolerance,
+ precursorTolerance,
SearchParameters.MassDiffAcceptorType,
SearchParameters.CustomMdac);
- // Run the classic search engine
+ // Run the classic search engine. When precomputedPeptides is supplied (from a .msl library),
+ // the engine iterates those peptides directly instead of digesting proteinsToSearch.
var searchEngine = new TransientClassicSearchEngine(
spectralMatchArray, AllSortedMs2Scans, VariableModifications,
FixedModifications, proteinsToSearch, massDiffAcceptor, CommonParameters,
- FileSpecificParameters, nestedIds, copyOnWriteEnabled: useCopyOnWrite);
+ FileSpecificParameters, nestedIds, copyOnWriteEnabled: useCopyOnWrite, precomputedPeptides: precomputedPeptides);
var results = searchEngine.Run();
updatedPsmIndexes = (results as TransientSearchEngineResults)!.UpdatedSpectralMatchIndexes;
@@ -643,6 +1127,16 @@ private void PerformSearch(List proteinsToSearch, SpectralMatch[] s
var transientPeptides = FilterToTransientDatabaseOnly(allPeptides, transientProteinAccessions).ToList();
_ = _peptideFdrAlignmentService.ApplyBaseline(transientPeptides);
+ // 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);
+ }
+
List? proteinGroups = null;
if (SearchParameters.DoParsimony && transientPsms.Count > 0)
{
@@ -750,14 +1244,26 @@ private void PerformSearch(List proteinsToSearch, SpectralMatch[] s
private void InitializeCompletedDatabaseWriter()
{
+ // The completed-database output (per-db folder + PSM/peptide/results files) was written by ONE
+ // consumer draining a capacity-2 channel, so at 1000s of databases the parallel searchers blocked
+ // on a single serial writer (low CPU utilization). Run several writer consumers — each writes a
+ // different database's folder, so file writes never collide — and widen the channel so searchers
+ // don't stall waiting for a write slot. The three pieces of shared bookkeeping the consumers touch
+ // are each synchronized: AppendCheckpoint -> ParallelSearchResultCache.AppendToFile (lock _writeLock
+ // + _cacheLock over a ConcurrentDictionary), MarkDatabaseCompleted (lock _dashboardLock), and
+ // UpdateProgress (lock _progressLock).
+ int writerCount = Math.Max(2, Environment.ProcessorCount / 4);
_completedDatabaseWriteChannel = Channel.CreateBounded<(TransientDatabaseContext Context, TransientDatabaseMetrics Metrics)>(
- new BoundedChannelOptions(2)
+ new BoundedChannelOptions(Math.Max(writerCount * 4, 32))
{
- SingleReader = true,
+ SingleReader = false, // multiple parallel writer consumers
SingleWriter = false,
FullMode = BoundedChannelFullMode.Wait
});
- _completedDatabaseWriterTask = Task.Run(RunCompletedDatabaseWriterLoopAsync);
+ var writers = new Task[writerCount];
+ for (int i = 0; i < writerCount; i++)
+ writers[i] = Task.Run(RunCompletedDatabaseWriterLoopAsync);
+ _completedDatabaseWriterTask = Task.WhenAll(writers);
}
private void CompleteCompletedDatabaseWriter()
@@ -789,6 +1295,50 @@ private async Task RunCompletedDatabaseWriterLoopAsync()
}
}
+ /// Confidence thresholds for what gets written to the per-database output files.
+ private const double OutputQValueThreshold = 0.05;
+ private const double OutputPepQValueThreshold = 0.05;
+
+ ///
+ /// A match is confident when its PEP-based q-value is below — the
+ /// PEP_QValue is assigned by mapping the match's model PEP onto the background curve (see
+ /// TransientPepAnalysisEngine), so it is meaningful even though the tiny transient databases can't compute
+ /// their own. When PEP is active a match also counts as confident if its score-based QValue clears the
+ /// threshold (confident by EITHER metric). When PEP is unavailable (no trained model), only the borrowed
+ /// score-based QValue is used.
+ ///
+ private bool IsConfident(SpectralMatch p, bool peptideLevel)
+ => IsConfidentMatch(p?.GetFdrInfo(peptideLevel), _pepEngine != null, OutputPepQValueThreshold, OutputQValueThreshold);
+
+ ///
+ /// Pure confidence decision (extracted for testing). When , a match is confident
+ /// if EITHER its PEP_QValue is below OR its score-based QValue is at or
+ /// below — the per-database files exist for the user to investigate any
+ /// reasonable hit, so an edge-case match that one metric flags and the other misses is still written. When
+ /// PEP is inactive, only the score-based QValue is used. A null FdrInfo is never confident.
+ ///
+ internal static bool IsConfidentMatch(EngineLayer.FdrAnalysis.FdrInfo info, bool pepActive, double pepQThreshold, double qThreshold)
+ {
+ if (info == null)
+ return false;
+ bool qConfident = info.QValue <= qThreshold;
+ return pepActive
+ ? info.PEP_QValue < pepQThreshold || qConfident
+ : qConfident;
+ }
+
+ ///
+ /// Row-level confidence filter (see ). Uses the peptide-level FdrInfo when
+ /// is set, otherwise the PSM-level one. Returns the input unchanged when
+ /// null/empty so writers behave as before for empty lists.
+ ///
+ private List FilterToConfident(List matches, bool peptideLevel)
+ {
+ if (matches == null || matches.Count == 0)
+ return matches;
+ return matches.Where(p => IsConfident(p, peptideLevel)).ToList();
+ }
+
private async Task WriteCompletedDatabaseOutputsAsync(TransientDatabaseContext context, TransientDatabaseMetrics metrics)
{
string dbName = context.DatabaseName;
@@ -796,28 +1346,54 @@ private async Task WriteCompletedDatabaseOutputsAsync(TransientDatabaseContext c
List nestedIds = context.NestedIds;
bool writeAllResults = !ParallelSearchParameters.WriteTransientResultsOnly;
+ // Write per-database output ONLY for databases with a transient result at <= 5% FDR. At 1000s of
+ // databases the vast majority match nothing (or only sub-threshold noise); writing a folder + files
+ // for each was a dominant cost. The database is still recorded in the cross-database summary/checkpoint
+ // below, so it counts as processed (and resume still works — HasCachedResults reads the checkpoint).
+ bool hasTransientOutput =
+ (context.TransientPeptides != null && context.TransientPeptides.Any(p => IsConfident(p, true)))
+ || (context.TransientPsms != null && context.TransientPsms.Any(p => IsConfident(p, false)));
+ if (!hasTransientOutput)
+ {
+ _resultsManager!.AppendCheckpoint(metrics);
+ MarkDatabaseCompleted();
+ UpdateProgress(TotalDatabases, nestedIds[0]);
+ ReportDatabaseDashboard(nestedIds[0], ParallelSearchDashboardUpdateKind.DatabaseFinished, dbName,
+ $"Finished {dbName} (no transient hits)", DashboardDatabaseFinishedProgress);
+ ReportProgress(new(100, $"Finished {dbName} (no transient hits)", nestedIds));
+ return;
+ }
+
+ // Create the output folder lazily — only databases with transient hits get one.
+ if (!Directory.Exists(outputFolder))
+ Directory.CreateDirectory(outputFolder);
+
Status($"Writing results for {dbName}...", nestedIds);
ReportDatabaseDashboard(nestedIds[0], ParallelSearchDashboardUpdateKind.DatabaseProgress, dbName,
$"Writing results for {dbName}...", DashboardDatabaseWritingProgress);
+ // Output is limited to confident matches (q-value <= OutputQValueThreshold). The filtering is row-level:
+ // every file below contains only matches at or below 5% FDR.
+ var confidentTransientPsms = FilterToConfident(context.TransientPsms, peptideLevel: false);
string transientPsmFile = Path.Combine(outputFolder,
$"{dbName}_All{GlobalVariables.AnalyteType.GetSpectralMatchLabel()}s.{GlobalVariables.AnalyteType.GetSpectralMatchExtension()}");
- await WritePsmsToTsvAsync(context.TransientPsms, transientPsmFile, SearchParameters.ModsToWriteSelection, false);
+ await WritePsmsToTsvAsync(confidentTransientPsms, transientPsmFile, SearchParameters.ModsToWriteSelection, false);
FinishedWritingFile(transientPsmFile, nestedIds);
if (writeAllResults)
{
string psmFile = Path.Combine(outputFolder,
$"All{GlobalVariables.AnalyteType.GetSpectralMatchLabel()}s.{GlobalVariables.AnalyteType.GetSpectralMatchExtension()}");
- await WritePsmsToTsvAsync(context.AllPsms, psmFile, SearchParameters.ModsToWriteSelection, false);
+ await WritePsmsToTsvAsync(FilterToConfident(context.AllPsms, peptideLevel: false), psmFile, SearchParameters.ModsToWriteSelection, false);
FinishedWritingFile(psmFile, nestedIds);
}
- if (context.TransientPeptides is { Count: > 0 })
+ var confidentTransientPeptides = FilterToConfident(context.TransientPeptides, peptideLevel: true);
+ if (confidentTransientPeptides is { Count: > 0 })
{
string transientPeptideFile = Path.Combine(outputFolder,
$"{dbName}_All{GlobalVariables.AnalyteType}s.{GlobalVariables.AnalyteType.GetSpectralMatchExtension()}");
- await WritePsmsToTsvAsync(context.TransientPeptides, transientPeptideFile, SearchParameters.ModsToWriteSelection, true);
+ await WritePsmsToTsvAsync(confidentTransientPeptides, transientPeptideFile, SearchParameters.ModsToWriteSelection, true);
FinishedWritingFile(transientPeptideFile, nestedIds);
}
@@ -825,7 +1401,7 @@ private async Task WriteCompletedDatabaseOutputsAsync(TransientDatabaseContext c
{
string peptideFile = Path.Combine(outputFolder,
$"All{GlobalVariables.AnalyteType}s.{GlobalVariables.AnalyteType.GetSpectralMatchExtension()}");
- await WritePsmsToTsvAsync(context.AllPeptides, peptideFile, SearchParameters.ModsToWriteSelection, true);
+ await WritePsmsToTsvAsync(FilterToConfident(context.AllPeptides, peptideLevel: true), peptideFile, SearchParameters.ModsToWriteSelection, true);
FinishedWritingFile(peptideFile, nestedIds);
}
@@ -898,17 +1474,17 @@ private void WriteFinalOutputs(string outputFolder, string taskId, int numFiles)
// Write global summary text file
WriteGlobalResultsText(_resultsManager.TransientDatabaseMetricsDictionary, outputFolder, taskId, numFiles);
- // Deal with custom reduced database writing
- // TODO: Revise this to only be the family one.
- bool useFamilyRanking = ParallelSearchParameters.UseFamilyAwareRanking;
- var statsByDatabase = useFamilyRanking
- ? SelectDatabasesForWritingByFamily()
- : SelectDatabasesForWritingByTestRatio();
+ // Deal with custom reduced database writing.
+ // Always use the family-aware gate (combined q-value + a minimum number of evidence families).
+ // The old test-ratio gate required passing >=50% of ALL tests, which even a perfect spike-in
+ // (SARS-CoV-2 passed 31/78) can't clear — many sub-tests structurally cannot fire for a small
+ // genome (NullEvidence / Undefined / BelowEligibilityThreshold), so nothing was ever written.
+ var statsByDatabase = SelectDatabasesForWritingByFamily();
Task[] dbWritingTasks = new Task[3];
if (statsByDatabase.Count > 0)
{
- Log($"Found {statsByDatabase.Count} significant databases passing cutoff (mode: {(useFamilyRanking ? "family-aware" : "test-ratio")})", [taskId]);
+ Log($"Found {statsByDatabase.Count} significant databases passing cutoff (family-aware, >={MinFamiliesForSignificance}/7 families)", [taskId]);
dbWritingTasks[0] = ParallelSearchParameters.DatabasesToWriteAndSearch[DatabaseToProduce.AllSignificantOrganisms].Write
? Task.Run(() => CreateCombinedDatabaseWithAllProteins(taskId, statsByDatabase.Select(p => p.Key), outputFolder))
@@ -941,30 +1517,56 @@ private Dictionary> SelectDatabasesForWri
p => p.OrderBy(t => t.ToString()).ToList());
}
+ /// Minimum number (of 7) of independent evidence families a database must pass to be written
+ /// as a confidently-detected organism.
+ /// NOTE on tuning: on the SARS virus spike-in, SARS-CoV-2 passes 7/7 and SARS-CoV 5/7, while a tail of
+ /// phage databases that floor out at the combined-q value only reach 4/7. So a bar of 4 admits that phage
+ /// tail (more sensitive, noisier) and a bar of 5 isolates the genuine detections. Left at 4 by request;
+ /// raise to 5 if the phage tail proves to be false positives.
+ private const int MinFamiliesForSignificance = 4;
+
private Dictionary> SelectDatabasesForWritingByFamily()
{
double qValueThreshold = CommonParameters.QValueThreshold;
- int minFamilyPasses = Math.Max(1, ParallelSearchParameters.TestRatioForWriting > 0
- ? (int)Math.Ceiling(ParallelSearchParameters.TestRatioForWriting * 7)
- : 1);
+ int minFamilyPasses = MinFamiliesForSignificance;
+ // Resolve the source DbForTask for a database tag. In merged-index mode TransientDatabases holds only the
+ // single merged .msl, so there is no per-organism DbForTask — synthesize one keyed by the db tag.
+ DbForTask ResolveDb(string dbTag) =>
+ ParallelSearchParameters.TransientDatabases.FirstOrDefault(db => Path.GetFileNameWithoutExtension(db.FileName) == dbTag)
+ ?? new DbForTask(dbTag, false);
return _resultsManager!.StatisticalTestResultList
.GroupBy(p => p.DatabaseName)
- .Where(p =>
- {
- if (!_resultsManager!.TransientDatabaseMetricsDictionary.TryGetValue(p.Key, out var metrics))
- return false;
-
- int passedFamilies = metrics.PassedFamilyCount;
- double combinedQ = metrics.CombinedQValue;
+ .Where(g => QualifiesAsDetectedOrganism(g, minFamilyPasses, qValueThreshold))
+ .ToDictionary(g => ResolveDb(g.Key), g => g.OrderBy(t => t.ToString()).ToList());
+ }
- return passedFamilies >= minFamilyPasses
- && !double.IsNaN(combinedQ)
- && combinedQ <= qValueThreshold;
- })
- .ToDictionary(
- p => ParallelSearchParameters.TransientDatabases.First(db => Path.GetFileNameWithoutExtension(db.FileName) == p.Key),
- p => p.OrderBy(t => t.ToString()).ToList());
+ ///
+ /// The family-aware detection predicate (extracted for testing): a database qualifies as a confidently
+ /// detected organism when at least independent evidence families are
+ /// significant AND the overall combined q-value is at or below . Both
+ /// quantities are computed DIRECTLY from the test results — the same source StatisticalAnalysis_Results.csv
+ /// uses — rather than from the metrics-summary fields, which are not reliably populated for this writer
+ /// (and were silently selecting nothing, even for SARS-CoV-2 at 7/7).
+ ///
+ internal static bool QualifiesAsDetectedOrganism(
+ IEnumerable dbResults, int minFamilyPasses, double qValueThreshold,
+ double significanceAlpha = 0.05, string overallCombinedMetric = "All")
+ {
+ var results = dbResults as ICollection ?? dbResults.ToList();
+ int passedFamilies = results
+ .Where(r => !r.IsCombinedResult && r.EvidenceFamily.HasValue && r.IsSignificant(significanceAlpha))
+ .Select(r => r.EvidenceFamily!.Value)
+ .Distinct()
+ .Count();
+ double combinedQ = results
+ .Where(r => r.IsCombinedResult && r.MetricName == overallCombinedMetric)
+ .Select(r => r.QValue)
+ .DefaultIfEmpty(double.NaN)
+ .First();
+ return passedFamilies >= minFamilyPasses
+ && !double.IsNaN(combinedQ)
+ && combinedQ <= qValueThreshold;
}
diff --git a/MetaMorpheus/TaskLayer/ParallelSearch/Statistics/IsolationForest/AnomalyDetectionFeatureBuilder.cs b/MetaMorpheus/TaskLayer/ParallelSearch/Statistics/IsolationForest/AnomalyDetectionFeatureBuilder.cs
index 5fcc31afd..85225e0f3 100644
--- a/MetaMorpheus/TaskLayer/ParallelSearch/Statistics/IsolationForest/AnomalyDetectionFeatureBuilder.cs
+++ b/MetaMorpheus/TaskLayer/ParallelSearch/Statistics/IsolationForest/AnomalyDetectionFeatureBuilder.cs
@@ -121,6 +121,9 @@ private static void SetFamilyFeature(double[] vec, int index, double pValue)
vec[index] = NegLog10WithSentinel(pValue, AnomalySentinels.MissingPValue);
}
+ /// Replaces NaN/Infinity with 0 so feature vectors are always valid IsolationForest input.
+ private static double SanitizeFinite(double v) => double.IsNaN(v) || double.IsInfinity(v) ? 0.0 : v;
+
private static double NegLog10WithSentinel(double value, double sentinel)
{
if (double.IsNaN(value) || double.IsInfinity(value) || value <= 0)
@@ -144,11 +147,15 @@ private static List ScaleAndBuild(
for (int f = 0; f < featureCount; f++)
{
- var values = rawFeatures.Select(r => r.RawFeatures[f]).ToList();
+ // Sanitize raw values first (a non-null Inf effect size, or a 0/0 ratio, can slip through the
+ // per-feature sentinels) so the robust-scaling statistics aren't corrupted by NaN/Infinity.
+ var values = rawFeatures.Select(r => SanitizeFinite(r.RawFeatures[f])).ToList();
var scaledValues = ScaleFeature(values).ToList();
+ // And sanitize the scaled output — IsolationForestInput rejects any NaN/Infinity, and with
+ // thousands of (mostly empty) databases a single stray value would abort the whole analysis.
for (int i = 0; i < rawFeatures.Count; i++)
- scaled[i][f] = scaledValues[i];
+ scaled[i][f] = SanitizeFinite(scaledValues[i]);
}
var result = new List(rawFeatures.Count);
diff --git a/MetaMorpheus/TaskLayer/Properties/AssemblyInfo.cs b/MetaMorpheus/TaskLayer/Properties/AssemblyInfo.cs
new file mode 100644
index 000000000..c1e4d7b8f
--- /dev/null
+++ b/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")]
diff --git a/MetaMorpheus/TaskLayer/TaskLayer.csproj b/MetaMorpheus/TaskLayer/TaskLayer.csproj
index e2b793186..a1fa9df37 100644
--- a/MetaMorpheus/TaskLayer/TaskLayer.csproj
+++ b/MetaMorpheus/TaskLayer/TaskLayer.csproj
@@ -39,7 +39,7 @@
-
+
diff --git a/MetaMorpheus/Test/ParallelSearchTask/Analysis/BasicMetricCollectorTests.cs b/MetaMorpheus/Test/ParallelSearchTask/Analysis/BasicMetricCollectorTests.cs
index 99941fa62..2374a7c0c 100644
--- a/MetaMorpheus/Test/ParallelSearchTask/Analysis/BasicMetricCollectorTests.cs
+++ b/MetaMorpheus/Test/ParallelSearchTask/Analysis/BasicMetricCollectorTests.cs
@@ -73,4 +73,44 @@ public void CollectData_UsesMinimumThresholdAcrossPsmAndPeptide()
Assert.That(result[BasicMetricCollector.TargetPeptidesFromTransientDbAtQValueThreshold], Is.EqualTo(1));
});
}
+
+ [Test]
+ public void CollectData_PepQValueCounts_AtOneAndFivePercent_ExcludeDecoys()
+ {
+ var cp = ParallelSearchTestContextFactory.CreateCommonParameters(qValueThreshold: 0.01, pepQValueThreshold: 0.05);
+
+ // peptideQValue is wired to PeptideFdrInfo.PEP_QValue; psmQValue to PsmFdrInfo.PEP_QValue.
+ var transientPeptides = new List
+ {
+ ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 30, psmQValue: 0.5, peptideQValue: 0.005, scanNumber: 1), // < 1% and < 5%
+ ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 25, psmQValue: 0.5, peptideQValue: 0.03, scanNumber: 2), // < 5% only
+ ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 20, psmQValue: 0.5, peptideQValue: 0.10, scanNumber: 3), // neither
+ ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: true, score: 18, psmQValue: 0.5, peptideQValue: 0.001, scanNumber: 4), // decoy -> excluded
+ };
+ var transientPsms = new List
+ {
+ ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 30, psmQValue: 0.005, peptideQValue: 0.5, scanNumber: 5),
+ ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 25, psmQValue: 0.03, peptideQValue: 0.5, scanNumber: 6),
+ ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: true, score: 18, psmQValue: 0.001, peptideQValue: 0.5, scanNumber: 7), // decoy -> excluded
+ };
+
+ var context = ParallelSearchTestContextFactory.CreateContext(
+ cp,
+ transientPsms,
+ transientPsms: transientPsms,
+ transientPeptides,
+ transientPeptides: transientPeptides,
+ totalProteins: 10,
+ transientPeptideCount: 4);
+
+ var result = new BasicMetricCollector().CollectData(context);
+
+ Assert.Multiple(() =>
+ {
+ Assert.That(result[BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ01], Is.EqualTo(1));
+ Assert.That(result[BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ05], Is.EqualTo(2));
+ Assert.That(result[BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ01], Is.EqualTo(1));
+ Assert.That(result[BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ05], Is.EqualTo(2));
+ });
+ }
}
diff --git a/MetaMorpheus/Test/ParallelSearchTask/Analysis/TransientDatabaseMetricsPepColumnsTests.cs b/MetaMorpheus/Test/ParallelSearchTask/Analysis/TransientDatabaseMetricsPepColumnsTests.cs
new file mode 100644
index 000000000..87d3fb4ab
--- /dev/null
+++ b/MetaMorpheus/Test/ParallelSearchTask/Analysis/TransientDatabaseMetricsPepColumnsTests.cs
@@ -0,0 +1,48 @@
+using NUnit.Framework;
+using TaskLayer.ParallelSearch.Analysis;
+using TaskLayer.ParallelSearch.Analysis.Collectors;
+
+namespace Test.ParallelSearchTask.Analysis;
+
+///
+/// Verifies the PEP-based 1% / 5% confident-count columns survive the property <-> Results-dictionary
+/// round-trip used to (de)serialize TransientDatabaseMetrics to/from CSV.
+///
+[TestFixture]
+public class TransientDatabaseMetricsPepColumnsTests
+{
+ [Test]
+ public void PepQValueColumns_RoundTripThroughResultsDictionary()
+ {
+ var metrics = new TransientDatabaseMetrics("UP000464024")
+ {
+ TargetPsmsFromTransientDbAtPepQ01 = 91,
+ TargetPsmsFromTransientDbAtPepQ05 = 96,
+ TargetPeptidesFromTransientDbAtPepQ01 = 30,
+ TargetPeptidesFromTransientDbAtPepQ05 = 31,
+ };
+
+ metrics.PopulateResultsFromProperties();
+
+ // The four new keys must be present in the serialized dictionary with their values.
+ Assert.Multiple(() =>
+ {
+ Assert.That(metrics.Results[BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ01], Is.EqualTo(91));
+ Assert.That(metrics.Results[BasicMetricCollector.TargetPsmsFromTransientDbAtPepQ05], Is.EqualTo(96));
+ Assert.That(metrics.Results[BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ01], Is.EqualTo(30));
+ Assert.That(metrics.Results[BasicMetricCollector.TargetPeptidesFromTransientDbAtPepQ05], Is.EqualTo(31));
+ });
+
+ // And they must read back into a fresh metrics object.
+ var restored = new TransientDatabaseMetrics("UP000464024") { Results = metrics.Results };
+ restored.PopulatePropertiesFromResults();
+
+ Assert.Multiple(() =>
+ {
+ Assert.That(restored.TargetPsmsFromTransientDbAtPepQ01, Is.EqualTo(91));
+ Assert.That(restored.TargetPsmsFromTransientDbAtPepQ05, Is.EqualTo(96));
+ Assert.That(restored.TargetPeptidesFromTransientDbAtPepQ01, Is.EqualTo(30));
+ Assert.That(restored.TargetPeptidesFromTransientDbAtPepQ05, Is.EqualTo(31));
+ });
+ }
+}
diff --git a/MetaMorpheus/Test/ParallelSearchTask/ConfidenceFilterTests.cs b/MetaMorpheus/Test/ParallelSearchTask/ConfidenceFilterTests.cs
new file mode 100644
index 000000000..abf2b9511
--- /dev/null
+++ b/MetaMorpheus/Test/ParallelSearchTask/ConfidenceFilterTests.cs
@@ -0,0 +1,49 @@
+using EngineLayer.FdrAnalysis;
+using NUnit.Framework;
+using PST = TaskLayer.ParallelSearch.ParallelSearchTask;
+
+namespace Test.ParallelSearchTask;
+
+///
+/// Tests the row-level output confidence gate. When a PEP model is active a match is confident if EITHER its
+/// PEP_QValue is strictly below 5% OR its (inclusive) score-based QValue clears 5% — the per-database files
+/// exist so the user can investigate any reasonable hit, so an edge-case match that only one metric flags is
+/// still written. When PEP is inactive, only the score-based QValue is used.
+///
+[TestFixture]
+public class ConfidenceFilterTests
+{
+ [Test]
+ public void NullInfo_IsNeverConfident()
+ {
+ Assert.That(PST.IsConfidentMatch(null, pepActive: true, pepQThreshold: 0.05, qThreshold: 0.05), Is.False);
+ Assert.That(PST.IsConfidentMatch(null, pepActive: false, pepQThreshold: 0.05, qThreshold: 0.05), Is.False);
+ }
+
+ [Test]
+ public void PepActive_ConfidentWhenPepQValueStrictlyBelowThreshold()
+ {
+ // Poor score QValue, but PEP clears -> confident on the PEP axis.
+ Assert.That(PST.IsConfidentMatch(new FdrInfo { PEP_QValue = 0.049, QValue = 0.9 }, true, 0.05, 0.05), Is.True);
+ // PEP at the boundary is NOT confident on the PEP axis (strict <), and the score QValue is poor.
+ Assert.That(PST.IsConfidentMatch(new FdrInfo { PEP_QValue = 0.05, QValue = 0.9 }, true, 0.05, 0.05), Is.False);
+ // Both metrics fail -> not confident.
+ Assert.That(PST.IsConfidentMatch(new FdrInfo { PEP_QValue = 0.06, QValue = 0.9 }, true, 0.05, 0.05), Is.False);
+ }
+
+ [Test]
+ public void PepActive_ConfidentByEitherMetric()
+ {
+ // Poor PEP_QValue but excellent score-based QValue -> confident when PEP is active (either metric).
+ Assert.That(PST.IsConfidentMatch(new FdrInfo { PEP_QValue = 0.5, QValue = 0.001 }, true, 0.05, 0.05), Is.True);
+ // Score QValue at the inclusive boundary also clears.
+ Assert.That(PST.IsConfidentMatch(new FdrInfo { PEP_QValue = 0.5, QValue = 0.05 }, true, 0.05, 0.05), Is.True);
+ }
+
+ [Test]
+ public void PepInactive_FallsBackToQValue_InclusiveThreshold()
+ {
+ Assert.That(PST.IsConfidentMatch(new FdrInfo { QValue = 0.05, PEP_QValue = 2.0 }, false, 0.05, 0.05), Is.True); // <= inclusive
+ Assert.That(PST.IsConfidentMatch(new FdrInfo { QValue = 0.051, PEP_QValue = 0.0 }, false, 0.05, 0.05), Is.False);
+ }
+}
diff --git a/MetaMorpheus/Test/ParallelSearchTask/FamilyDetectionTests.cs b/MetaMorpheus/Test/ParallelSearchTask/FamilyDetectionTests.cs
new file mode 100644
index 000000000..eb0d6ffd7
--- /dev/null
+++ b/MetaMorpheus/Test/ParallelSearchTask/FamilyDetectionTests.cs
@@ -0,0 +1,125 @@
+using System.Collections.Generic;
+using NUnit.Framework;
+using TaskLayer.ParallelSearch.Statistics;
+using PST = TaskLayer.ParallelSearch.ParallelSearchTask;
+
+namespace Test.ParallelSearchTask;
+
+///
+/// Tests the family-aware organism-detection predicate: a database is a confident detection when at least
+/// N independent evidence families are significant AND the overall combined q-value clears the threshold.
+/// This replaced the test-ratio gate that required passing >=50% of ALL tests — unreachable for a small
+/// genome (SARS-CoV-2 passed 31/78), so nothing was ever written.
+///
+[TestFixture]
+public class FamilyDetectionTests
+{
+ private const int MinFamilies = 4;
+ private const double QThreshold = 0.01;
+
+ private static StatisticalTestResult Family(StatisticalEvidenceFamily fam, double qValue) => new()
+ {
+ TestName = "GaussianTest", // not a combined name -> IsCombinedResult == false
+ MetricName = fam.ToString(),
+ EvidenceFamily = fam,
+ IsDefined = true,
+ QValue = qValue,
+ };
+
+ private static StatisticalTestResult OverallCombined(double qValue, string metric = "All") => new()
+ {
+ TestName = "Combined", // IsCombinedResult == true
+ MetricName = metric,
+ IsDefined = true,
+ QValue = qValue,
+ };
+
+ private static readonly StatisticalEvidenceFamily[] SevenFamilies =
+ {
+ StatisticalEvidenceFamily.CountEnrichment,
+ StatisticalEvidenceFamily.AmbiguityOrTargetDecoy,
+ StatisticalEvidenceFamily.ScoreDistribution,
+ StatisticalEvidenceFamily.Fragmentation,
+ StatisticalEvidenceFamily.RetentionTime,
+ StatisticalEvidenceFamily.ProteinGroup,
+ StatisticalEvidenceFamily.PrecursorDeconvolution,
+ };
+
+ [Test]
+ public void StrongDetection_AllFamiliesSignificant_Qualifies()
+ {
+ var results = new List();
+ foreach (var f in SevenFamilies) results.Add(Family(f, 1e-300));
+ results.Add(OverallCombined(4.7e-300));
+
+ Assert.That(PST.QualifiesAsDetectedOrganism(results, MinFamilies, QThreshold), Is.True);
+ }
+
+ [Test]
+ public void TooFewFamilies_DoesNotQualify()
+ {
+ var results = new List
+ {
+ Family(StatisticalEvidenceFamily.CountEnrichment, 1e-50),
+ Family(StatisticalEvidenceFamily.Fragmentation, 1e-50),
+ Family(StatisticalEvidenceFamily.RetentionTime, 1e-50), // only 3 distinct families
+ OverallCombined(1e-100),
+ };
+ Assert.That(PST.QualifiesAsDetectedOrganism(results, MinFamilies, QThreshold), Is.False);
+ }
+
+ [Test]
+ public void CombinedQOverThreshold_DoesNotQualify()
+ {
+ var results = new List();
+ for (int i = 0; i < 5; i++) results.Add(Family(SevenFamilies[i], 1e-50));
+ results.Add(OverallCombined(0.02)); // 5 families pass but combined q 0.02 > 0.01
+
+ Assert.That(PST.QualifiesAsDetectedOrganism(results, MinFamilies, QThreshold), Is.False);
+ }
+
+ [Test]
+ public void NoOverallCombinedResult_NaN_DoesNotQualify()
+ {
+ var results = new List();
+ for (int i = 0; i < 5; i++) results.Add(Family(SevenFamilies[i], 1e-50)); // families pass, but no "All" combined
+
+ Assert.That(PST.QualifiesAsDetectedOrganism(results, MinFamilies, QThreshold), Is.False);
+ }
+
+ [Test]
+ public void InsignificantFamilyResults_DoNotCountTowardFamilies()
+ {
+ var results = new List();
+ foreach (var f in SevenFamilies) results.Add(Family(f, 0.5)); // present but NOT significant
+ results.Add(OverallCombined(1e-300));
+
+ Assert.That(PST.QualifiesAsDetectedOrganism(results, MinFamilies, QThreshold), Is.False);
+ }
+
+ [Test]
+ public void PerFamilyCombined_NotTreatedAsOverall()
+ {
+ // A combined result for a single family (MetricName != "All") must NOT supply the overall combined q.
+ var results = new List();
+ for (int i = 0; i < 5; i++) results.Add(Family(SevenFamilies[i], 1e-50));
+ results.Add(OverallCombined(1e-300, metric: "CountEnrichment")); // per-family combined, not the overall
+
+ Assert.That(PST.QualifiesAsDetectedOrganism(results, MinFamilies, QThreshold), Is.False);
+ }
+
+ [Test]
+ public void DuplicateSignificantFamily_CountsOnce()
+ {
+ // Same family significant twice should count as ONE family, so 3 distinct families < 4 -> no.
+ var results = new List
+ {
+ Family(StatisticalEvidenceFamily.CountEnrichment, 1e-50),
+ Family(StatisticalEvidenceFamily.CountEnrichment, 1e-60), // duplicate family
+ Family(StatisticalEvidenceFamily.Fragmentation, 1e-50),
+ Family(StatisticalEvidenceFamily.RetentionTime, 1e-50),
+ OverallCombined(1e-100),
+ };
+ Assert.That(PST.QualifiesAsDetectedOrganism(results, MinFamilies, QThreshold), Is.False);
+ }
+}
diff --git a/MetaMorpheus/Test/ParallelSearchTask/MslCandidateFilterTests.cs b/MetaMorpheus/Test/ParallelSearchTask/MslCandidateFilterTests.cs
new file mode 100644
index 000000000..0401b8e27
--- /dev/null
+++ b/MetaMorpheus/Test/ParallelSearchTask/MslCandidateFilterTests.cs
@@ -0,0 +1,111 @@
+using Chemistry;
+using NUnit.Framework;
+using Priors = EngineLayer.ParallelSearch.MslPeptideReader.CandidatePriors;
+using Reader = EngineLayer.ParallelSearch.MslPeptideReader;
+
+namespace Test.ParallelSearchTask;
+
+///
+/// Tests the .msl candidate pre-filter — the learned mass + RT gate that decides which library entries are
+/// worth fragmenting/scoring — plus the merged-index "db|accession" parse and the lower-bound helper.
+///
+[TestFixture]
+public class MslCandidateFilterTests
+{
+ private const int Charge = 2;
+ private const double Mz = 500.0;
+
+ // Build priors whose scan list has exactly one mass-matching scan (the middle one) at a chosen RT.
+ private static Priors PriorsWithMatch(double matchingScanRt, double tolPpm = 10, double rtWindowMin = 5,
+ double rtSlope = 1, double rtIntercept = 0)
+ {
+ double neutral = Mz.ToMass(Charge);
+ var scanMasses = new[] { neutral - 100.0, neutral, neutral + 100.0 }; // only the middle one is within tol
+ var scanRts = new[] { 0.0, matchingScanRt, 0.0 };
+ return new Priors(scanMasses, scanRts, tolPpm, rtSlope, rtIntercept, rtWindowMin);
+ }
+
+ [Test]
+ public void IsCandidate_MassMatch_UnpredictedRt_KeepsOnMassAlone()
+ {
+ var priors = PriorsWithMatch(matchingScanRt: 20.0);
+ // irt == 0 -> RT filter disabled; mass matches -> candidate.
+ bool cand = Reader.IsCandidate(Mz, Charge, 0f, priors, out bool massMatched);
+ Assert.That(cand, Is.True);
+ Assert.That(massMatched, Is.True);
+ }
+
+ [Test]
+ public void IsCandidate_NoMassMatch_NotCandidate()
+ {
+ var priors = PriorsWithMatch(matchingScanRt: 20.0);
+ // m/z + 25 -> neutral + 50, which lands in the 100-Da gap between scan masses -> no mass match.
+ bool cand = Reader.IsCandidate(Mz + 25.0, Charge, 0f, priors, out bool massMatched);
+ Assert.That(cand, Is.False);
+ Assert.That(massMatched, Is.False);
+ }
+
+ [Test]
+ public void IsCandidate_MassMatch_RtInWindow_IsCandidate()
+ {
+ var priors = PriorsWithMatch(matchingScanRt: 20.0); // slope 1, intercept 0 -> predRt == iRT
+ bool cand = Reader.IsCandidate(Mz, Charge, 20.0f, priors, out bool massMatched); // predRt 20 vs scan RT 20
+ Assert.That(cand, Is.True);
+ Assert.That(massMatched, Is.True);
+ }
+
+ [Test]
+ public void IsCandidate_MassMatch_RtOutOfWindow_NotCandidate_ButMassMatched()
+ {
+ var priors = PriorsWithMatch(matchingScanRt: 20.0, rtWindowMin: 5);
+ // predRt 50 vs the only mass-matching scan's RT 20 -> |50-20|=30 > 5 -> RT fails.
+ bool cand = Reader.IsCandidate(Mz, Charge, 50.0f, priors, out bool massMatched);
+ Assert.That(cand, Is.False);
+ Assert.That(massMatched, Is.True, "mass matched even though RT did not");
+ }
+
+ [Test]
+ public void IsCandidate_MassTolerance_BoundaryBehaviour()
+ {
+ double neutral = Mz.ToMass(Charge);
+ // The effective window is (ppm * mass) + a 0.01 Da margin. A scan 0.03 Da away is outside a 10 ppm
+ // window (~0.02 Da total) but inside a 40 ppm one (~0.05 Da total).
+ double offset = 0.03;
+ var priorsNarrow = new Priors(new[] { neutral + offset }, new[] { 0.0 }, precursorTolPpm: 10, rtSlope: 1, rtIntercept: 0, rtWindowMin: 5);
+ Assert.That(Reader.IsCandidate(Mz, Charge, 0f, priorsNarrow, out _), Is.False);
+
+ var priorsWide = new Priors(new[] { neutral + offset }, new[] { 0.0 }, precursorTolPpm: 40, rtSlope: 1, rtIntercept: 0, rtWindowMin: 5);
+ Assert.That(Reader.IsCandidate(Mz, Charge, 0f, priorsWide, out _), Is.True);
+ }
+
+ [Test]
+ public void LowerBound_ReturnsFirstIndexAtOrAboveValue()
+ {
+ var sorted = new[] { 1.0, 3.0, 5.0, 9.0 };
+ Assert.Multiple(() =>
+ {
+ Assert.That(Reader.LowerBound(sorted, 0.0), Is.EqualTo(0)); // below all
+ Assert.That(Reader.LowerBound(sorted, 3.0), Is.EqualTo(1)); // exact match
+ Assert.That(Reader.LowerBound(sorted, 4.0), Is.EqualTo(2)); // between -> next-higher index
+ Assert.That(Reader.LowerBound(sorted, 9.0), Is.EqualTo(3)); // exact last
+ Assert.That(Reader.LowerBound(sorted, 10.0), Is.EqualTo(4)); // above all -> length
+ });
+ }
+
+ [Test]
+ public void ParseDbTagAndAccession_SplitsOnFirstBar()
+ {
+ Assert.Multiple(() =>
+ {
+ Assert.That(Reader.ParseDbTagAndAccession("UP000464024|P0DTC9"), Is.EqualTo(("UP000464024", "P0DTC9")));
+ // Only the FIRST bar splits; later bars stay in the accession.
+ Assert.That(Reader.ParseDbTagAndAccession("db|acc|extra"), Is.EqualTo(("db", "acc|extra")));
+ // No bar -> UNKNOWN db, whole string is the accession.
+ Assert.That(Reader.ParseDbTagAndAccession("loneAccession"), Is.EqualTo(("UNKNOWN", "loneAccession")));
+ // Empty accession after the bar -> "_UNKNOWN".
+ Assert.That(Reader.ParseDbTagAndAccession("db|"), Is.EqualTo(("db", "db_UNKNOWN")));
+ // Null/empty input.
+ Assert.That(Reader.ParseDbTagAndAccession(null), Is.EqualTo(("UNKNOWN", "UNKNOWN_UNKNOWN")));
+ });
+ }
+}
diff --git a/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchEndToEndTests.cs b/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchEndToEndTests.cs
new file mode 100644
index 000000000..f4b32d3ec
--- /dev/null
+++ b/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchEndToEndTests.cs
@@ -0,0 +1,71 @@
+using System.Collections.Generic;
+using System.IO;
+using System.Linq;
+using EngineLayer.DatabaseLoading;
+using NUnit.Framework;
+using TaskLayer;
+using PSTask = TaskLayer.ParallelSearch.ParallelSearchTask;
+
+namespace Test.ParallelSearchTask;
+
+///
+/// End-to-end smoke test that drives a real ParallelSearchTask once over small bundled data, so the
+/// orchestration that the focused unit tests can't reach — base search, PEP training, the per-database
+/// search + PEP assignment + confident-output writing, the statistics, and the family-aware writer — is
+/// exercised. Base db = smalldb.fasta (the search that feeds PEP training); transient db = gapdh.fasta
+/// (abundant yeast GAPDH, so the transient search produces hits).
+///
+[TestFixture]
+public class ParallelSearchEndToEndTests
+{
+ [Test]
+ public void RunTask_SmallYeast_ProducesSummaryWithPepColumns()
+ {
+ string testData = Path.Combine(TestContext.CurrentContext.TestDirectory, "TestData");
+ string outputFolder = Path.Combine(TestContext.CurrentContext.TestDirectory, "ParallelSearchE2E");
+ if (Directory.Exists(outputFolder))
+ Directory.Delete(outputFolder, true);
+ Directory.CreateDirectory(outputFolder);
+
+ // Work on copies of the spectra file (the search writes file-specific artifacts next to it). Two copies
+ // are pooled into the base search so it clears the >=100-PSM bar that gates PEP-model training.
+ string src = Path.Combine(testData, "SmallCalibratible_Yeast.mzML");
+ string mzml1 = Path.Combine(outputFolder, "Yeast_1.mzML");
+ string mzml2 = Path.Combine(outputFolder, "Yeast_2.mzML");
+ File.Copy(src, mzml1, true);
+ File.Copy(src, mzml2, true);
+
+ string baseDb = Path.Combine(testData, "smalldb.fasta"); // base (yeast) -> PEP training set
+ string transientDb = Path.Combine(testData, "smalldb.fasta"); // transient db (same yeast proteins) -> guaranteed hits
+
+ var task = new PSTask(new List { new DbForTask(transientDb, false) });
+ task.ParallelSearchParameters.MaxSearchesInParallel = 1;
+
+ Assert.DoesNotThrow(() =>
+ task.RunTask(outputFolder, new List { new DbForTask(baseDb, false) }, new List { mzml1, mzml2 }, "test"));
+
+ // The cross-database summary must have been written, and must carry the new PEP 1%/5% columns.
+ var summary = Directory.EnumerateFiles(outputFolder, "ManySearchSummary.csv", SearchOption.AllDirectories).FirstOrDefault();
+ Assert.That(summary, Is.Not.Null, "ManySearchSummary.csv was not written");
+
+ var lines = File.ReadAllLines(summary);
+ var header = lines[0].Split(',');
+ int Col(string name) => System.Array.IndexOf(header, name);
+ Assert.Multiple(() =>
+ {
+ Assert.That(Col("TargetPeptidesFromTransientDbAtPepQ01"), Is.GreaterThanOrEqualTo(0));
+ Assert.That(Col("TargetPeptidesFromTransientDbAtPepQ05"), Is.GreaterThanOrEqualTo(0));
+ Assert.That(Col("TargetPsmsFromTransientDbAtPepQ01"), Is.GreaterThanOrEqualTo(0));
+ Assert.That(Col("TargetPsmsFromTransientDbAtPepQ05"), Is.GreaterThanOrEqualTo(0));
+ });
+
+ var row = lines.Skip(1).First(l => l.Trim().Length > 0).Split(',');
+ int transientPep = int.Parse(row[Col("TargetPeptidesFromTransientDb")]);
+ int pepQ05 = int.Parse(row[Col("TargetPeptidesFromTransientDbAtPepQ05")]);
+
+ // The transient search produced hits, and the base-trained PEP model assigned a confident PEP_QValue to
+ // them (this is what fails if PEP training or the per-database PEP-assignment block stops running).
+ Assert.That(transientPep, Is.GreaterThanOrEqualTo(1), "transient search should produce hits");
+ Assert.That(pepQ05, Is.GreaterThanOrEqualTo(1), "PEP model should assign a confident PEP_QValue to transient peptides");
+ }
+}
diff --git a/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchParametersTests.cs b/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchParametersTests.cs
index 7205f30f9..274f55c2e 100644
--- a/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchParametersTests.cs
+++ b/MetaMorpheus/Test/ParallelSearchTask/ParallelSearchParametersTests.cs
@@ -20,6 +20,7 @@ public void Constructor_SetsParallelSearchDefaults()
Assert.That(parameters.MaxSearchesInParallel, Is.EqualTo(4));
Assert.That(parameters.WriteTransientResultsOnly, Is.True);
Assert.That(parameters.WriteTransientSpectralLibrary, Is.False);
+ Assert.That(parameters.UseMergedTransientLibrary, Is.False);
Assert.That(parameters.DoParsimony, Is.True);
Assert.That(parameters.NoOneHitWonders, Is.True);
Assert.That(parameters.MassDiffAcceptorType, Is.EqualTo(MassDiffAcceptorType.Exact));
diff --git a/MetaMorpheus/Test/ParallelSearchTask/PepBackgroundCurveTests.cs b/MetaMorpheus/Test/ParallelSearchTask/PepBackgroundCurveTests.cs
new file mode 100644
index 000000000..fc905eaaf
--- /dev/null
+++ b/MetaMorpheus/Test/ParallelSearchTask/PepBackgroundCurveTests.cs
@@ -0,0 +1,126 @@
+using System.Collections.Generic;
+using System.Linq;
+using EngineLayer;
+using EngineLayer.ParallelSearch;
+using NUnit.Framework;
+using Test.ParallelSearchTask.Utility;
+
+namespace Test.ParallelSearchTask;
+
+///
+/// Tests for the background PEP -> PEP_QValue curve used to assign a calibrated PEP q-value to transient
+/// database IDs (which are far too small to compute their own PEP target/decoy). The curve is snapshotted
+/// from the decoy-rich base search; each transient match's model PEP is mapped onto it by lower-bound lookup.
+///
+[TestFixture]
+public class PepBackgroundCurveTests
+{
+ // A monotone curve: PEP ascending, PEP_QValue non-decreasing (the shape BuildPepQValueCurve produces).
+ private static readonly double[] PepAsc = { 0.01, 0.05, 0.20, 0.60 };
+ private static readonly double[] QByPep = { 0.001, 0.01, 0.05, 0.30 };
+
+ [Test]
+ public void Lookup_EmptyCurve_ReturnsSentinelTwo()
+ {
+ double q = TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.01, System.Array.Empty(), System.Array.Empty());
+ Assert.That(q, Is.EqualTo(2.0));
+ }
+
+ [Test]
+ public void Lookup_PepBelowAll_ReturnsLowestQValue()
+ {
+ // Most-confident transient PEP, below every background PEP -> the smallest (best) background q-value.
+ double q = TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.005, PepAsc, QByPep);
+ Assert.That(q, Is.EqualTo(0.001));
+ }
+
+ [Test]
+ public void Lookup_PepAboveAll_ReturnsHighestQValue()
+ {
+ // Worse than every background PEP -> the largest (worst) background q-value (clamped to last).
+ double q = TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.95, PepAsc, QByPep);
+ Assert.That(q, Is.EqualTo(0.30));
+ }
+
+ [Test]
+ public void Lookup_ExactPep_ReturnsThatQValue()
+ {
+ Assert.That(TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.01, PepAsc, QByPep), Is.EqualTo(0.001));
+ Assert.That(TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.20, PepAsc, QByPep), Is.EqualTo(0.05));
+ Assert.That(TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.60, PepAsc, QByPep), Is.EqualTo(0.30));
+ }
+
+ [Test]
+ public void Lookup_PepBetweenPoints_ReturnsLowerBoundQValue()
+ {
+ // 0.03 falls between background PEPs 0.01 and 0.05; lower-bound is 0.05 -> its q-value 0.01.
+ Assert.That(TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.03, PepAsc, QByPep), Is.EqualTo(0.01));
+ // 0.50 falls between 0.20 and 0.60; lower-bound is 0.60 -> 0.30.
+ Assert.That(TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.50, PepAsc, QByPep), Is.EqualTo(0.30));
+ }
+
+ [Test]
+ public void Lookup_SingleElementCurve_AlwaysReturnsThatQValue()
+ {
+ var pepAsc = new[] { 0.10 };
+ var qByPep = new[] { 0.02 };
+ Assert.That(TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.001, pepAsc, qByPep), Is.EqualTo(0.02));
+ Assert.That(TransientPepAnalysisEngine.LookupBackgroundPepQValue(0.99, pepAsc, qByPep), Is.EqualTo(0.02));
+ }
+
+ [Test]
+ public void Lookup_IsMonotoneNonDecreasingInPep()
+ {
+ // Sweeping PEP upward never yields a smaller q-value than a lower PEP did.
+ double prev = -1;
+ for (double pep = 0.0; pep <= 1.0; pep += 0.02)
+ {
+ double q = TransientPepAnalysisEngine.LookupBackgroundPepQValue(pep, PepAsc, QByPep);
+ Assert.That(q, Is.GreaterThanOrEqualTo(prev));
+ prev = q;
+ }
+ }
+
+ [Test]
+ public void BuildCurve_EmptyInput_ReturnsEmptyArrays()
+ {
+ var (pepAsc, qByPep) = TransientPepAnalysisEngine.BuildPepQValueCurve(new List(), peptideLevel: false);
+ Assert.That(pepAsc, Is.Empty);
+ Assert.That(qByPep, Is.Empty);
+ }
+
+ [Test]
+ public void BuildCurve_TargetsAndDecoys_ProducesSortedMonotoneCurve()
+ {
+ var cp = ParallelSearchTestContextFactory.CreateCommonParameters();
+ // Targets with low PEP (confident), decoys with high PEP (poor) — the realistic separation.
+ var matches = new List();
+ double[] targetPeps = { 0.001, 0.01, 0.05, 0.10 };
+ double[] decoyPeps = { 0.85, 0.90, 0.95, 0.99 };
+ int scan = 1;
+ foreach (var p in targetPeps)
+ {
+ var m = ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 20, psmQValue: 0.001, peptideQValue: 0.001, scanNumber: scan++);
+ m.PsmFdrInfo.PEP = p;
+ matches.Add(m);
+ }
+ foreach (var p in decoyPeps)
+ {
+ var m = ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: true, score: 20, psmQValue: 0.5, peptideQValue: 0.5, scanNumber: scan++);
+ m.PsmFdrInfo.PEP = p;
+ matches.Add(m);
+ }
+
+ var (pepAsc, qByPep) = TransientPepAnalysisEngine.BuildPepQValueCurve(matches, peptideLevel: false);
+
+ Assert.That(pepAsc.Length, Is.EqualTo(matches.Count));
+ // PEP axis sorted ascending.
+ for (int i = 1; i < pepAsc.Length; i++)
+ Assert.That(pepAsc[i], Is.GreaterThanOrEqualTo(pepAsc[i - 1]), "PEP axis must be ascending");
+ // q-value axis monotone non-decreasing.
+ for (int i = 1; i < qByPep.Length; i++)
+ Assert.That(qByPep[i], Is.GreaterThanOrEqualTo(qByPep[i - 1]), "PEP_QValue must be non-decreasing in PEP");
+ // The most-confident (lowest-PEP) target has a better q-value than the worst decoy.
+ Assert.That(qByPep.First(), Is.LessThan(qByPep.Last()));
+ }
+}
diff --git a/MetaMorpheus/Test/ParallelSearchTask/PepTransientAssignmentTests.cs b/MetaMorpheus/Test/ParallelSearchTask/PepTransientAssignmentTests.cs
new file mode 100644
index 000000000..c72753af6
--- /dev/null
+++ b/MetaMorpheus/Test/ParallelSearchTask/PepTransientAssignmentTests.cs
@@ -0,0 +1,66 @@
+using System.Collections.Generic;
+using EngineLayer;
+using EngineLayer.ParallelSearch;
+using NUnit.Framework;
+using Test.ParallelSearchTask.Utility;
+
+namespace Test.ParallelSearchTask;
+
+///
+/// Tests assigning PEP / PEP_QValue to transient-database matches from the base-trained model. The transient
+/// databases are far too small to train their own model, so one model trained on the (decoy-rich) base search
+/// is reused and each transient match's PEP_QValue is mapped onto the background curve.
+///
+[TestFixture]
+public class PepTransientAssignmentTests
+{
+ private static List<(string, CommonParameters)> Fsp(CommonParameters cp) => new() { ("file.raw", cp) };
+
+ [Test]
+ public void AssignPepFromTrainedModel_NoModel_IsNoOp()
+ {
+ var cp = ParallelSearchTestContextFactory.CreateCommonParameters();
+ var basePsms = new List
+ {
+ ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 20, psmQValue: 0.001, peptideQValue: 0.001),
+ };
+ var engine = new TransientPepAnalysisEngine(basePsms, "standard", Fsp(cp), TestContext.CurrentContext.TestDirectory);
+ Assert.That(engine.HasTrainedModel, Is.False);
+
+ var transient = ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 20, psmQValue: 0.5, peptideQValue: 0.5, scanNumber: 2);
+ transient.PeptideFdrInfo.PEP_QValue = 1.23; // a sentinel we expect to remain untouched
+
+ engine.AssignPepFromTrainedModel(new List { transient }, peptideLevel: true);
+
+ Assert.That(transient.PeptideFdrInfo.PEP_QValue, Is.EqualTo(1.23), "no trained model -> no-op");
+ }
+
+ [Test]
+ public void TrainAndAssign_CreatesMissingFdrInfo_AndAssignsPepQValue()
+ {
+ var cp = ParallelSearchTestContextFactory.CreateCommonParameters();
+ var basePsms = new List();
+ for (int i = 0; i < 12; i++)
+ basePsms.Add(ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 20 + i, psmQValue: 0.001, peptideQValue: 0.001, scanNumber: i + 1));
+ for (int i = 0; i < 12; i++)
+ basePsms.Add(ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: true, score: 5 + i, psmQValue: 0.5, peptideQValue: 0.5, scanNumber: 100 + i));
+
+ var engine = new TransientPepAnalysisEngine(basePsms, "standard", Fsp(cp), TestContext.CurrentContext.TestDirectory);
+ bool trained = engine.TrainSingleModelAndAssignBasePep();
+ // Need both target and decoy training examples; if the minimal synthetic features are degenerate, skip.
+ Assume.That(trained, Is.True);
+ Assert.That(engine.HasTrainedModel, Is.True);
+
+ // A transient peptide with NO PeptideFdrInfo: AssignPepFromTrainedModel must create it (this was the
+ // exact bug that left every transient PEP at the sentinel) and assign a mapped PEP_QValue.
+ var transient = ParallelSearchTestContextFactory.CreateSpectralMatch(cp, isDecoy: false, score: 25, psmQValue: 0.5, peptideQValue: 0.5, scanNumber: 500);
+ transient.PeptideFdrInfo = null;
+ transient.PsmFdrInfo = null;
+
+ engine.AssignPepFromTrainedModel(new List { transient }, peptideLevel: true);
+
+ Assert.That(transient.PeptideFdrInfo, Is.Not.Null, "missing PeptideFdrInfo must be created");
+ Assert.That(transient.PsmFdrInfo, Is.Not.Null, "missing PsmFdrInfo must be created");
+ Assert.That(transient.PeptideFdrInfo.PEP_QValue, Is.Not.EqualTo(2.0), "PEP_QValue must be assigned from the background curve");
+ }
+}
diff --git a/MetaMorpheus/Test/Test.csproj b/MetaMorpheus/Test/Test.csproj
index 37d8f8a6d..db8a42080 100644
--- a/MetaMorpheus/Test/Test.csproj
+++ b/MetaMorpheus/Test/Test.csproj
@@ -24,7 +24,7 @@
-
+