Skip to content

Commit e5ab254

Browse files
trishortsclaudeMICHAEL SHORTREED
authored
Promote consensus mass-tracing types + add Ms1Feature writer for MassFeature (#1069)
* add new protease subtilisin * feat: add subtilisin|p protease and unit test Add subtilisin|p to the embedded proteases.tsv with full cleavage specificity and proline-inhibition motifs (N[P]|, S[P]|, L[P]|, K[P]|, I[P]|, D[P]|, Y[P]|, V[P]|, G[P]|, F[P]|, T[P]|, E[P]|, Q[P]|, A[P]|, R[P]|). Add TestSubtilisinP_DigestsCorrectlyAndRespectsProlineRestriction to ProteinDigestionTests to verify: - subtilisin|p is present in the embedded protease dictionary with CleavageSpecificity.Full - All expected cleavage sites fire on a proline-free sequence (ANKTIDE) - The K[P]| proline-inhibition rule is respected (AKPIDE keeps KP intact) * Accept newer-TopFD column names in Ms1Feature PR 1064 ships an Ms1Feature reader that recognises only the older FlashDeconv / TopFD-v1.6.2 _ms1.feature schema (Sample_ID, ID, Time_begin, Time_end, Minimum_charge_state, Maximum_charge_state, Minimum_fraction_id, Maximum_fraction_id). Newer TopFD output keeps the same _ms1.feature extension but uses different column names (File_name, Fraction_ID, Feature_ID, Min_time, Max_time, Min_charge, Max_charge), plus extras like Envelope_num and EC_score. Format detection picks Ms1FeatureFile by extension; CsvHelper then throws because none of the expected [Name(...)] columns are present. Discovered while integrating PR 1064's FromFileDeconvolutionParameters into MetaMorpheus and pointing it at a real TopFD .ms1.feature from a top-down yeast run -- the file parsed by hand looks identical in shape to the old schema, just relabelled. Fix is column-name aliases on the existing record, plus [Optional] on fields that the newer schema omits entirely. Downstream GetSingleChargeFeatures() reads only Mass, ChargeStateMin/Max, RetentionTimeBegin/End, and IntensityApex -- all aliased to a column present in both schemas, so the join algorithm behaves identically regardless of which producer wrote the file. Per-field summary: Sample_ID -> [Optional] (newer TopFD has File_name instead; type-incompatible -- int vs string path -- and not used downstream) ID -> alias "Feature_ID" + [Optional] Time_begin -> alias "Min_time" Time_end -> alias "Max_time" Minimum_charge_state -> alias "Min_charge" Maximum_charge_state -> alias "Max_charge" Minimum_fraction_id -> [Optional] (newer TopFD has a single Fraction_ID column, not a min/max pair; not used downstream so alias would add no value) Maximum_fraction_id -> [Optional] No tests added in this commit -- a follow-up should drop a newer-TopFD _ms1.feature sample into Test/FileReadingTests/ExternalFileTypes/ and extend the existing Ms1FeatureFile read-roundtrip tests to cover both schemas. The current Ms1Feature_FlashDeconvOpenMs3.0.0_ms1.feature Ms1Feature_TopFDv1.6.2_ms1.feature fixtures keep passing because the existing [Name(...)] heads remain the first entry in every alias list. * Drops a 4-row fixture (Ms1Feature_TopFDvLatest_ms1.feature) captured from real TopFD output that uses the File_name / Fraction_ID / Feature_ID / Min_time / Max_time / Min_charge / Max_charge schema, and wires it into the existing TestMsFeature parameterised tests: * TestFeaturesLoadAndCountIsCorrect gains a TestCase asserting the fixture loads four features end-to-end via FileReader. * TestTopFDLatestMs1FeatureFirstAndLastAreCorrect locks the per- field mapping for both the aliased columns (Time_begin/Min_time, Minimum_charge_state/Min_charge, etc.) and the [Optional] fields absent from the newer schema (SampleId, FractionIdMin/Max default to 0). Covers the single-charge edge case (Min_charge == Max_charge == 1) in the last record. * TestTopFDLatestMs1GetSingleChargeFeatureFunctions confirms charge- range expansion is identical across the two TopFD schemas: a 6-14 range yields 9 envelopes; a 1-1 range yields exactly one; GetMs1Features() flattens to 9 + 12 + 10 + 1 = 32 across the four fixture features. * TestMs1FeatureReadWrite gains the new fixture as a TestCase. The writer emits the older-schema headers (newer columns aren't on the record), so the round-trip converts schema-newer -> schema- older + Optional defaults. Comment in the test explains why that is correct: every field downstream consumers actually read survives the round-trip; the columns that don't are exactly the ones marked [Optional] and unused. All 17 TestMsFeature tests pass, as do the 148 tests covering the related Ms1Feature / FromFile / SupportedFileExtensions / DinosaurTsv surface area. * Promote consensus mass-tracing types out of Development Move the post-decon consensus pipeline (per-charge trace grouping, median off-by-one correction, cross-charge feature stitching) from Development.Deconvolution.ConsensusTracing into a new production namespace MassSpectrometry.Deconvolution.Consensus. The types had been parked in Development as the research arc explored their design (see NOTES.md, phases 1-13); promotion was flagged as the next move once the design stabilised. New files under mzLib/MassSpectrometry/Deconvolution/Consensus/: * MassTrace -- per-scan envelope list at one charge, anchored to the first envelope's mass. * MassTraceBuilder -- greedy charge-locked grouper, gap-tolerant within a configurable mass tolerance. * CorrectedEnvelope -- one envelope inside a corrected trace, carrying both original and post-correction mass plus a WasCorrected flag. * CorrectedTrace -- a MassTrace after weighted-median off-by-one correction. * TraceCorrector -- the corrector itself, with uniform/intensity/ scorer-weighted variants via an EnvelopeWeight delegate. * MassFeature -- a cross-charge feature: a group of CorrectedTrace entries whose consensus masses agree within a ppm tolerance and whose RT windows overlap. * MassFeatureBuilder -- the union-find sweep that produces them. No behaviour change. The classes, methods, and signatures are identical to the research-namespace versions; only their namespace and file locations changed (each type now lives in its own file). Doc comments updated to drop references to phase numbers and to the "Development" location. The research scaffolding in Development.Deconvolution.ConsensusTracing (Phase*.cs fixtures, MakeConsensusSnips, LargeTestDataLocator) stays on the consensus-tracing branch; it will be updated to consume the new namespace when that branch is rebased on top of this PR. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Add Ms1FeatureFile.FromMassFeatures factory Lets the consensus-tracing pipeline persist its cross-charge feature list as a FLASHDeconv-style _ms1.feature file -- the existing wire format any Ms1FeatureFile consumer (notably MetaMorpheus's FromFileDeconvolutionParameters) already understands. The consensus-to-search bridge then needs no new file format and no new MetaMorpheus surface area. Two pieces: * MassFeature.ToMs1Feature(sequentialId, sampleId = 0, fractionId = 0) extension method (under Readers, because the return type is the Readers-namespace Ms1Feature record). Maps: Mass <- MassFeature.ConsensusMass Intensity <- MassFeature.SummedIntensity RetentionTimeBegin <- MassFeature.RTStart RetentionTimeEnd <- MassFeature.RTEnd RetentionTimeApex <- RT of the highest-intensity envelope on the highest-summed-intensity constituent trace ("apex of the dominant charge state at its most intense scan") IntensityApex <- intensity of that same envelope ChargeStateMin/Max <- min/max of MassFeature.Charges SampleId / FractionId <- callers' choice (defaults 0) * Static Ms1FeatureFile.FromMassFeatures(IEnumerable<MassFeature>, sampleId, fractionId, software) factory. Builds the records via ToMs1Feature, populates a fresh Ms1FeatureFile instance, and returns it. The caller invokes .WriteResults(path) on the returned file to actually emit bytes -- no new IO path, the existing reader's writer is reused. Software label defaults to FLASHDeconv because that's the canonical schema the writer emits; the [Name(...)] alias machinery from PR #1064 + schema-aliases PR means newer-TopFD readers also see the output as valid. No filter applied: every MassFeature becomes one row. The original question surfaced four filtering options (all, multi-charge only, score-thresholded, all + side artifact); per the design pin "all features, no filter" is the right default. Downstream consumers (MetaMorpheus's precursor HashSet + FDR machinery) handle the noise. The writer is decoupled from the producer: callers pass a sequence of finalised MassFeature objects; the factory doesn't care where they came from. Future producers (mass-trace results from non- Classic algorithms, etc.) can use the same path. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Cover Ms1FeatureFile.FromMassFeatures with a round-trip test Five NUnit cases under Test\FileReadingTests\ExternalFileReading\TestMs1FeatureFromMassFeature.cs: * ToMs1Feature_SingletonFeature_MapsAllFields -- one trace, one envelope. Locks the mapping: Mass, Intensity, RT begin/end/apex, apex intensity, charge bounds, fraction id all come from Finalise()-derived MassFeature fields. * ToMs1Feature_MultiChargeFeature_ApexIsDominantTraceMaxEnvelope -- two-charge feature with five envelopes. Confirms apex is "highest- intensity envelope on the highest-summed-intensity constituent trace", and that cross-charge consensus mass is the intensity- weighted mean of per-trace consensus masses. * ToMs1Feature_SampleIdAndFractionIdHonoured -- callers can override Sample_ID and Fraction_ID for multi-file/multi-fraction outputs. * FromMassFeatures_WriteThenRead_AllFieldsSurvive -- three features built in-memory, written via WriteResults, re-read with FileReader.ReadFile<Ms1FeatureFile>, every downstream-consumed field compared to the original ToMs1Feature output. Catches any serialisation drift between writer and reader on the new path. * FromMassFeatures_AssignsSequentialIdsFromZero -- producer's internal MassFeature.Id is ignored; written rows get 0..N-1. Keeps the output file's IDs dense and stable even when the upstream pipeline has filtered features and left holes. Fix discovered while writing the second-to-last test: the original FromMassFeatures factory used `file.Results.Add(...)` in a loop, which goes through ResultFile<T>.Results' lazy-load getter. The getter calls LoadResults() whenever the backing list is empty; for an in-memory- constructed Ms1FeatureFile with no FilePath, LoadResults throws ArgumentException on the empty path. Reshaped the factory to build the records up front, then assign Results once via the setter (which bypasses the lazy-load). Inline comment explains why. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Add [Explicit] end-to-end driver from mzML to _ms1.feature Demonstrates the consensus-to-search loop the rest of this PR was built for: load an mzML, run Classic decon per MS1 scan with top-down parameters (charge 1..60), build mass traces with Pass-B grouping (loose 1.5 Da to cohort off-by-one twins), apply weighted-median correction, stitch cross-charge features at 10 ppm, write the resulting MassFeature list as a FLASHDeconv-style _ms1.feature via the new Ms1FeatureFile.FromMassFeatures factory, and round-trip read the written file as a sanity check. Lives as an NUnit [Explicit] test because it consumes a 165 MB local mzML at a hard-coded path that's not part of the repo or CI: E:\TestData\MetaMorpheus\05-26-17_B7A_yeast_td_fract7_rep1.mzML Assume.That on File.Exists keeps it silently skipped on machines without the data. Run with: dotnet test --filter FullyQualifiedName~TestConsensusToMs1FeatureEndToEnd On the yeast top-down test file (4837 scans total, ~2700 MS1) the driver writes a 24 MB _ms1.feature with 266,373 features. That file, fed to MetaMorpheus's FromFile precursor path (#2650), yields a Cal -> GPTMD -> Search pipeline that converges to a tighter precursor mass tolerance (1.6 ppm) than the equivalent classic-only baseline (2.4 ppm) and finds 907 PSMs at 1% FDR vs 785 -- a 15.5% increase in identifications from consensus precursors alone, before any additive combination with the classic source. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Parameterise consensus driver via TestCaseSource Lets a single dotnet test invocation run the consensus pipeline over multiple raw/mzML files in succession. Each case captures its own source path, output path, and decon charge ceiling -- bottom-up uses MaxCharge=12 (matching MetaMorpheus's standard BU default), top-down would use 60. Currently shipping two bottom-up HEK293 Velos raws as the parameterised test cases: * 20100609_Velos1_TaGe_SA_293_3.raw (5163 MS1 scans, 76878 features after correction, 7756 multi-charge -- 8.9 s wall time) * 20100609_Velos1_TaGe_SA_293_4.raw (5111 MS1 scans, 74173 features, 7318 multi-charge -- 5.7 s) Each test case still gates on File.Exists via Assume.That; cases for machines without the data skip silently. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Normalise feature-file RT to minutes when loading FlashDeconv and TopFD both canonically emit _ms1.feature files with RetentionTime columns in seconds (e.g. Time_begin = 2787 for a 46-min LC run). MetaMorpheus and mzLib's MsDataScan.RetentionTime are always in minutes (mzML / Thermo convention). FromFileDeconvolutionAlgorithm compared the file's RT directly against the scan's RT, so on any real FlashDeconv or TopFD output the seconds-vs-minutes mismatch produced zero overlapping windows -- and zero PSMs at search time. Caught while running a Cal -> GPTMD -> Search pipeline over a 20-mzML Jurkat top-down dataset using its FlashDeconv companion _ms1.feature files. Classic-only produced 23,636 PSMs at 1% FDR; FromFile-only produced 0 with no error or warning -- silent failure mode. Fix: in the file-path FromFileDeconvolutionParameters constructor, detect seconds-as-loaded via a max-RetentionTimeEnd > 500 heuristic (no realistic LC run exceeds 8 hours = 500 min) and divide all RT fields by 60. Files already in minutes pass through unchanged. The internal pre-loaded-features constructor (test seam) is unaffected so callers supplying explicit synthetic units still get raw values. Verification: re-running the same Jurkat pipeline after the fix yields 23,976 PSMs in FromFile-only mode and 39,700 PSMs in additive Both mode (vs Classic-only's 23,636) -- +68% PSMs / +68% proteoforms over the Classic baseline on this dataset. Tests: * Existing EndToEnd_FilePathCtor_ResolvesExpectedChargeState updated to use MS2 RT in minutes (39.85) instead of seconds (2390.0) -- the fixtures are FlashDeconv / TopFD with seconds-RT and the fix now normalises them at load, so the test had to match the new units. * New FromFileDeconvolutionParameters_FileWithRtInSeconds_NormalisesToMinutes writes a synthetic seconds-RT file, loads via the file-path ctor, asserts every per-charge feature's RT lands in the expected minutes range (~39.83 / ~40.17), and pairs against an in-minutes MS2 to confirm overlap is found. * New FromFileDeconvolutionParameters_FileWithRtInMinutes_LeavesUnchanged sister test: file already in minutes (max RT < 500) is NOT double-converted. 31/31 tests in TestFromFileDeconvolution + TestMs1FeatureFromMassFeature pass after the fix. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * feat(decon): split resolvable co-grouped species in trace correction TraceCorrector now derives its off-by-one window from per-trace mass scatter (3*sigma via MAD, capped below half the isotope spacing) instead of a fixed +/-0.05 Da. An envelope at a resolvable, non-isotope offset (e.g. a deamidated form at +0.98402 Da, which the old fixed window snapped to the unmodified mass) is split into its own CorrectedTrace and surfaced as a distinct feature; when scatter is too large to resolve it from the 1.00335 Da isotope spacing it conservatively merges. Correct() now returns List<CorrectedTrace>. Bottom-up search needs this -- deamidation detection is a requirement there. Also in this hardening pass: - MassFeatureBuilder sorts components by (mass, charge, RT) before assigning IDs, so feature IDs and written _ms1.feature row order are deterministic instead of depending on Dictionary enumeration. - Guard clauses with clear messages on Finalise (empty traces), ToMs1Feature (un-finalised feature), and BuildTraces (scan/envelope length mismatch). - FromFileDeconvolutionParameters exposes RetentionTimeNormalizedFromSeconds and warns when the seconds->minutes heuristic fires, rather than silently mutating RT. - Unit + performance tests for the consensus engines (0 -> 100% line coverage on the promoted types) covering the behaviours above. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * test(decon): add TD/BU feature generators + per-feature decon scoring [Explicit], local-data drivers used to benchmark the consensus pipeline: JurkatMs1FeatureGenerator and VelosMs1FeatureGenerator regenerate the _ms1.feature files for the Jurkat top-down (charge 1..60) and HEK293 Velos bottom-up (charge 1..12) datasets via the new deamidation-split pipeline. ConsensusFeatureScoring computes the generic per-envelope deconvolution score (DeconvolutionScorer.ScoreEnvelope: averagine cosine, ppm error, peak completeness, intensity-ratio consistency) for every Classic envelope and aggregates it per feature into a sidecar TSV aligned to the _ms1.feature row order. This gives a uniform, search-independent feature-quality value for a computed noise-filter cutoff -- scored at generation time because the consensus pipeline keeps only mass+intensity per envelope and drops the peaks the scorer needs. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * fix(readers): preserve gapped charge sets on _ms1.feature round-trip Hardening fixes to the consensus MassFeature -> _ms1.feature path surfaced by code review. None were observed in normal use, but each is a real edge-case correctness or crash issue: - Write one row per contiguous charge run so a gapped set (e.g. {10,12,15}, when intermediate charges fall below the score cutoff) no longer reads back as the fabricated full range 10..15. Adds ToMs1Features. - FromMassFeatures: empty input writes a valid header-only file instead of crashing; also fixes a LoadResults infinite recursion when reading a zero-row file (the Results getter re-entered LoadResults on an empty set). - Ignore non-finite RetentionTimeEnd when sniffing seconds-vs-minutes, so a single NaN row can't flip an in-minutes file into a /60 conversion. - Default the produced file's Software to TopFD to match how it is re-detected on reload (the writer always emits Apex_intensity); fix the misleading "FLASHDeconv canonical" doc. - Document that per-charge Intensity comes from Apex_intensity, not the summed Intensity column. Adds NUnit coverage for each; tidies now-stale FLASHDeconv comments. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * refactor(consensus): use Constants.C13MinusC12 for isotope spacing TraceCorrector defined its own IsotopeSpacingDa = 1.00335, a less-precise copy of the canonical Chemistry.Constants.C13MinusC12 (1.00335483810). Alias the shared constant so the C12/C13 spacing has one source of truth; the 4.8e-6 Da change is far inside the off-by-one detection window, so results are unchanged. Addresses review feedback on #1069. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * fix(readers): only lazy-load Results when the file exists ResultFile<T>.Results lazy-loaded whenever _results was empty, which crashed for a factory-built file with an empty FilePath (e.g. a zero-feature Ms1FeatureFile.FromMassFeatures). Guarding the load with File.Exists(FilePath) lets an in-memory file return its set records directly, removing the need for Ms1FeatureFile's parallel _factoryRecords store so the records have a single source of truth. Addresses review feedback on #1069. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * test(consensus): add XIC-vs-MassTrace grouping comparison harness Two [Explicit] drivers (local Jurkat top-down data, not CI) that back the decision to keep MassTraceBuilder over reusing the existing GetAllXics: TestXicVsMassTraceComparison times both groupers head-to-head on the same deconvoluted envelopes (wall-clock, allocations, trace counts, grouping agreement); JurkatMs1FeatureGeneratorViaXic emits XIC-grouped _ms1.feature files through the identical downstream so end-to-end proteoform yield can be compared (GetAllXics measured -18% vs MassTraceBuilder). Refs #1069. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com> * Add consensus pipeline hierarchy overview to MassFeature doc comment Defines how IsotopicEnvelope, MassTrace, CorrectedTrace, and MassFeature nest and which builder produces each level, per review request. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> * test(rt): compare duplicate-peptide predictions within tolerance, not bit-exactly PredictRetentionTimeEquivalents_DuplicatePeptides_ReturnsSameValueForEach asserted exact equality via Distinct().Count() == 1. The Chronologer batch path runs the duplicate peptides as identical rows in a single libtorch forward pass, whose CPU conv/matmul kernels are not bit-reproducible across thread counts / CI runners, so identical rows can differ in the last float32 bits (intermittent CI failure: distinct == 2). The values are correct to float32 precision — the rest of this suite already pins batch agreement with .Within(1e-4). Assert every duplicate prediction is non-null and equal to the first within 1e-4, matching the suite's existing tolerance. Production code unchanged; the real heap-corruption non-determinism was already fixed in #1075. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011QVxqemPqC9BjvYtVg6aKY --------- Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com> Co-authored-by: MICHAEL SHORTREED <mrshortreed@wisc.edu>
1 parent 125c995 commit e5ab254

23 files changed

Lines changed: 2823 additions & 9 deletions
Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
namespace MassSpectrometry.Deconvolution.Consensus
2+
{
3+
/// <summary>
4+
/// One envelope inside a <see cref="CorrectedTrace"/>, carrying both its
5+
/// original (algorithm-reported) mass and the post-correction mass.
6+
/// <see cref="WasCorrected"/> is true iff <see cref="TraceCorrector"/>
7+
/// flagged this envelope as an off-by-one outlier and snapped its
8+
/// mass to the trace consensus.
9+
/// </summary>
10+
public sealed class CorrectedEnvelope
11+
{
12+
public int ScanIndex;
13+
public int ScanNumber;
14+
public double RT;
15+
public double OriginalMass;
16+
public double CorrectedMass;
17+
public int Charge;
18+
public double Intensity;
19+
public bool WasCorrected;
20+
}
21+
}
Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
using System.Collections.Generic;
2+
using System.Linq;
3+
4+
namespace MassSpectrometry.Deconvolution.Consensus
5+
{
6+
/// <summary>
7+
/// A <see cref="MassTrace"/> after off-by-one correction by
8+
/// <see cref="TraceCorrector"/>. <see cref="ConsensusMass"/> is the
9+
/// per-trace (weighted) median; envelopes whose original mass differed
10+
/// by approximately +/-1.00335 Da from it have been snapped to the
11+
/// consensus and have <see cref="CorrectedEnvelope.WasCorrected"/> set.
12+
/// Original masses are preserved on each envelope for diagnostic
13+
/// inspection.
14+
/// </summary>
15+
public sealed class CorrectedTrace
16+
{
17+
public int Id;
18+
public int Charge;
19+
public double ConsensusMass;
20+
public List<CorrectedEnvelope> Envelopes = new();
21+
public double OriginalSpread;
22+
public double CorrectedSpread;
23+
public int CorrectionCount => Envelopes.Count(e => e.WasCorrected);
24+
25+
public double FirstRT => Envelopes[0].RT;
26+
public double LastRT => Envelopes[^1].RT;
27+
public double TotalIntensity => Envelopes.Sum(e => e.Intensity);
28+
}
29+
}
Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
using System.Collections.Generic;
2+
using System.Linq;
3+
4+
namespace MassSpectrometry.Deconvolution.Consensus
5+
{
6+
/// <summary>
7+
/// CONSENSUS PIPELINE OVERVIEW (the type hierarchy, narrowest to widest).
8+
///
9+
/// The consensus mass-tracing pipeline turns raw per-scan deconvolution
10+
/// output into cross-charge features through four nested levels. Each level
11+
/// groups the one below it:
12+
///
13+
/// IsotopicEnvelope one deconvolved species (mass, intensity, charge)
14+
/// | reported in a SINGLE MS1 scan. Raw input; not a
15+
/// | type in this namespace.
16+
/// | grouped by <see cref="MassTraceBuilder"/> (charge-locked,
17+
/// | anchor-mass + scan-adjacency)
18+
/// v
19+
/// <see cref="MassTrace"/> the same species followed across adjacent
20+
/// | scans at ONE charge state. Its envelopes are the
21+
/// | raw per-scan tuples above.
22+
/// | corrected by <see cref="TraceCorrector"/> (off-by-one rescue +
23+
/// | resolution-aware splitting)
24+
/// v
25+
/// <see cref="CorrectedTrace"/> a MassTrace after correction, holding
26+
/// | <see cref="CorrectedEnvelope"/> entries (original
27+
/// | + corrected mass per envelope) and a per-trace
28+
/// | <see cref="CorrectedTrace.ConsensusMass"/>. Still
29+
/// | one charge state.
30+
/// | stitched by <see cref="MassFeatureBuilder"/> (cross-charge,
31+
/// | ppm mass agreement + RT overlap)
32+
/// v
33+
/// MassFeature THIS type: one species across ALL of its charge
34+
/// states. The widest grouping, and the unit the
35+
/// writer turns into an Ms1Feature row.
36+
///
37+
/// So: envelopes nest inside traces, traces (after correction) nest inside
38+
/// features. Charge is fixed within a trace and varies across a feature.
39+
///
40+
/// A cross-charge-state consensus feature: a group of
41+
/// <see cref="CorrectedTrace"/> entries (each at one charge state) whose
42+
/// consensus masses agree within a ppm tolerance and whose RT windows
43+
/// overlap.
44+
///
45+
/// A real proteoform or peptide produces envelopes at several charge
46+
/// states (BU peptides commonly +2/+3; TD proteoforms across +10..+15
47+
/// or wider). The trace builder is charge-locked, so each charge gets
48+
/// its own trace. <see cref="MassFeatureBuilder"/> stitches those
49+
/// per-charge traces back together. Charge multiplicity then becomes
50+
/// a confidence signal: a feature seen at multiple charges is
51+
/// corroborated by independent charge calculations; a single-charge,
52+
/// single-envelope feature is more likely to be noise.
53+
///
54+
/// Mutation policy: callers append to <see cref="Traces"/> during
55+
/// construction, then call <see cref="Finalise"/> exactly once to
56+
/// derive the aggregate fields. After Finalise, the feature is
57+
/// treated as read-only by downstream consumers.
58+
/// </summary>
59+
public sealed class MassFeature
60+
{
61+
public int Id;
62+
public List<CorrectedTrace> Traces = new();
63+
public double ConsensusMass;
64+
public HashSet<int> Charges = new();
65+
public int ChargeCount => Charges.Count;
66+
public int MaxTraceLength;
67+
public double RTStart;
68+
public double RTEnd;
69+
public double SummedIntensity;
70+
71+
/// <summary>
72+
/// Populate derived fields from the current <see cref="Traces"/>
73+
/// list. Idempotent; safe to re-run after the trace list changes.
74+
/// <see cref="ConsensusMass"/> is the intensity-weighted mean of
75+
/// per-trace consensus masses (heavier traces contribute more).
76+
/// </summary>
77+
public void Finalise()
78+
{
79+
if (Traces.Count == 0)
80+
throw new System.InvalidOperationException("Cannot Finalise a MassFeature with no traces.");
81+
82+
Charges = new HashSet<int>(Traces.Select(t => t.Charge));
83+
SummedIntensity = Traces.Sum(t => t.TotalIntensity);
84+
RTStart = Traces.Min(t => t.FirstRT);
85+
RTEnd = Traces.Max(t => t.LastRT);
86+
MaxTraceLength = Traces.Max(t => t.Envelopes.Count);
87+
88+
// Intensity-weighted mean of per-trace consensus masses. Heavier
89+
// traces (more intense, more confident) contribute more.
90+
double w = Traces.Sum(t => t.TotalIntensity);
91+
ConsensusMass = w == 0
92+
? Traces.Average(t => t.ConsensusMass)
93+
: Traces.Sum(t => t.ConsensusMass * t.TotalIntensity) / w;
94+
}
95+
}
96+
}
Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
using System.Collections.Generic;
2+
using System.Linq;
3+
4+
namespace MassSpectrometry.Deconvolution.Consensus
5+
{
6+
/// <summary>
7+
/// Stitches per-charge <see cref="CorrectedTrace"/> entries back into
8+
/// cross-charge <see cref="MassFeature"/> objects.
9+
///
10+
/// Algorithm:
11+
/// 1. Sort the corrected traces by ConsensusMass.
12+
/// 2. Sweep with a ppm-wide window; for each pair of traces inside
13+
/// the window, if (a) their charges differ and (b) their RT
14+
/// ranges overlap, union-find them.
15+
/// 3. Each connected component becomes a <see cref="MassFeature"/>.
16+
///
17+
/// Same-charge co-mass traces are intentionally NOT merged here: if
18+
/// two traces share a charge and mass, the trace builder would already
19+
/// have merged them when their scan ranges overlapped. If they didn't
20+
/// overlap (co-eluting twins or a Phase-2 split), they remain separate
21+
/// features.
22+
/// </summary>
23+
public static class MassFeatureBuilder
24+
{
25+
public static List<MassFeature> BuildFeatures(
26+
IReadOnlyList<CorrectedTrace> corrected,
27+
double massPpm)
28+
{
29+
// Sort traces by consensus mass to enable a single sliding-window pass.
30+
var sorted = corrected.OrderBy(t => t.ConsensusMass).ToArray();
31+
int n = sorted.Length;
32+
33+
// Union-find over trace indices.
34+
int[] parent = new int[n];
35+
for (int i = 0; i < n; i++) parent[i] = i;
36+
int Find(int x)
37+
{
38+
while (parent[x] != x) { parent[x] = parent[parent[x]]; x = parent[x]; }
39+
return x;
40+
}
41+
void Union(int x, int y)
42+
{
43+
int rx = Find(x), ry = Find(y);
44+
if (rx != ry) parent[rx] = ry;
45+
}
46+
47+
for (int i = 0; i < n; i++)
48+
{
49+
double anchorMass = sorted[i].ConsensusMass;
50+
double window = anchorMass * massPpm * 1e-6;
51+
for (int j = i + 1; j < n; j++)
52+
{
53+
double dm = sorted[j].ConsensusMass - anchorMass;
54+
if (dm > window) break;
55+
if (sorted[i].Charge == sorted[j].Charge) continue;
56+
if (!RtOverlap(sorted[i], sorted[j])) continue;
57+
Union(i, j);
58+
}
59+
}
60+
61+
// Collect connected components.
62+
var byRoot = new Dictionary<int, List<CorrectedTrace>>();
63+
for (int i = 0; i < n; i++)
64+
{
65+
int r = Find(i);
66+
if (!byRoot.TryGetValue(r, out var bucket))
67+
{
68+
bucket = new List<CorrectedTrace>();
69+
byRoot[r] = bucket;
70+
}
71+
bucket.Add(sorted[i]);
72+
}
73+
74+
// Finalise each component, then assign IDs in a deterministic order. Dictionary
75+
// enumeration is not order-stable, so without this sort the feature IDs -- and the
76+
// row order of the written _ms1.feature file -- would vary run-to-run on identical
77+
// input. Order by (consensus mass, then min charge, then RT start).
78+
var features = new List<MassFeature>(byRoot.Count);
79+
foreach (var (_, traces) in byRoot)
80+
{
81+
var f = new MassFeature { Traces = traces };
82+
f.Finalise();
83+
features.Add(f);
84+
}
85+
features.Sort((a, b) =>
86+
{
87+
int c = a.ConsensusMass.CompareTo(b.ConsensusMass);
88+
if (c != 0) return c;
89+
c = a.Charges.Min().CompareTo(b.Charges.Min());
90+
if (c != 0) return c;
91+
return a.RTStart.CompareTo(b.RTStart);
92+
});
93+
for (int i = 0; i < features.Count; i++)
94+
features[i].Id = i + 1;
95+
return features;
96+
}
97+
98+
private static bool RtOverlap(CorrectedTrace a, CorrectedTrace b)
99+
=> a.FirstRT <= b.LastRT && b.FirstRT <= a.LastRT;
100+
}
101+
}
Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
using System.Collections.Generic;
2+
using System.Linq;
3+
4+
namespace MassSpectrometry.Deconvolution.Consensus
5+
{
6+
/// <summary>
7+
/// A single mass trace: a sequence of per-scan envelopes that the grouper
8+
/// believes belong to the same species. Same charge across all entries,
9+
/// scan-adjacency &lt;= MaxGap, mass within tolerance of <see cref="AnchorMass"/>.
10+
///
11+
/// <see cref="MassTraceBuilder"/> populates these in scan-order. Off-by-one
12+
/// correction is then handled by <see cref="TraceCorrector"/>, which wraps
13+
/// each trace in a <see cref="CorrectedTrace"/> rather than mutating it.
14+
/// Mutation policy on the trace itself: contents are appended during
15+
/// construction and read-only thereafter.
16+
/// </summary>
17+
public sealed class MassTrace
18+
{
19+
public int Id;
20+
public int Charge;
21+
22+
/// <summary>
23+
/// First envelope's mass; never updated after trace creation. Keeping
24+
/// the anchor fixed bounds the trace to AnchorMass +/- tolerance and
25+
/// prevents drift over long traces.
26+
/// </summary>
27+
public double AnchorMass;
28+
29+
public List<(int ScanIndex, int ScanNumber, double RT, double Mass, double Intensity)> Envelopes
30+
= new();
31+
32+
public int LastScanIndex => Envelopes[^1].ScanIndex;
33+
public double MinMass => Envelopes.Min(e => e.Mass);
34+
public double MaxMass => Envelopes.Max(e => e.Mass);
35+
public double MassSpread => MaxMass - MinMass;
36+
}
37+
}
Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
using System.Collections.Generic;
2+
3+
namespace MassSpectrometry.Deconvolution.Consensus
4+
{
5+
/// <summary>
6+
/// Greedy mass-trace builder, charge-locked.
7+
///
8+
/// For each MS1 scan in input order, every envelope is offered to the
9+
/// open traces. An envelope joins a trace if (a) charges match, (b) the
10+
/// trace has been touched within MaxGap+1 scans of the current one, and
11+
/// (c) the envelope's mass is within tolerance of the trace's anchor.
12+
/// Best match (closest mass) wins ties. No drift: the anchor is the
13+
/// first envelope's mass, fixed for the life of the trace.
14+
///
15+
/// Inputs: per-scan envelope lists in scan order, aligned to a parallel
16+
/// list of <see cref="MsDataScan"/>. Output: every trace (open or closed
17+
/// at end-of-input) in a flat list. Length-1 singletons are included --
18+
/// downstream code applies length filters per its own taste.
19+
/// </summary>
20+
public static class MassTraceBuilder
21+
{
22+
public static List<MassTrace> BuildTraces(
23+
IReadOnlyList<MsDataScan> ms1Scans,
24+
IReadOnlyList<IReadOnlyList<IsotopicEnvelope>> perScanEnvelopes,
25+
double toleranceDa,
26+
int maxGap)
27+
{
28+
if (perScanEnvelopes.Count != ms1Scans.Count)
29+
throw new System.ArgumentException(
30+
$"perScanEnvelopes ({perScanEnvelopes.Count}) must align 1:1 with ms1Scans ({ms1Scans.Count}).",
31+
nameof(perScanEnvelopes));
32+
33+
var open = new List<MassTrace>();
34+
var closed = new List<MassTrace>();
35+
int nextId = 1;
36+
37+
for (int scanIdx = 0; scanIdx < ms1Scans.Count; scanIdx++)
38+
{
39+
// Retire open traces whose last-touched scan is too old.
40+
for (int i = open.Count - 1; i >= 0; i--)
41+
{
42+
int gap = scanIdx - open[i].LastScanIndex - 1;
43+
if (gap > maxGap)
44+
{
45+
closed.Add(open[i]);
46+
open.RemoveAt(i);
47+
}
48+
}
49+
50+
foreach (var env in perScanEnvelopes[scanIdx])
51+
{
52+
MassTrace best = null!;
53+
double bestDelta = double.MaxValue;
54+
foreach (var t in open)
55+
{
56+
if (t.Charge != env.Charge) continue;
57+
// Don't add a second envelope from the same scan to the same trace.
58+
if (t.LastScanIndex == scanIdx) continue;
59+
double d = System.Math.Abs(env.MonoisotopicMass - t.AnchorMass);
60+
if (d <= toleranceDa && d < bestDelta)
61+
{
62+
best = t;
63+
bestDelta = d;
64+
}
65+
}
66+
67+
var entry = (scanIdx,
68+
ms1Scans[scanIdx].OneBasedScanNumber,
69+
ms1Scans[scanIdx].RetentionTime,
70+
env.MonoisotopicMass,
71+
env.TotalIntensity);
72+
73+
if (best != null)
74+
{
75+
best.Envelopes.Add(entry);
76+
}
77+
else
78+
{
79+
var nt = new MassTrace
80+
{
81+
Id = nextId++,
82+
Charge = env.Charge,
83+
AnchorMass = env.MonoisotopicMass,
84+
};
85+
nt.Envelopes.Add(entry);
86+
open.Add(nt);
87+
}
88+
}
89+
}
90+
91+
closed.AddRange(open);
92+
return closed;
93+
}
94+
}
95+
}

0 commit comments

Comments
 (0)