Fd#1042
Open
trishorts wants to merge 89 commits into
Open
Conversation
Codecov Report❌ Patch coverage is Please upload reports for the commit 35e6eea to get more accurate results. Additional details and impacted files@@ Coverage Diff @@
## master #1042 +/- ##
==========================================
- Coverage 81.63% 81.54% -0.10%
==========================================
Files 369 347 -22
Lines 47605 45981 -1624
Branches 5649 5396 -253
==========================================
- Hits 38863 37495 -1368
+ Misses 7643 7424 -219
+ Partials 1099 1062 -37
🚀 New features to boost your workflow:
|
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)
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.
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.
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>
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 smith-chem-wisc#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>
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>
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>
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>
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>
…chem-wisc#1069 conflicts) Upstream independently added the same logical changes for two pieces of this PR: - subtilisin entry in proteases.tsv. Upstream named the protease "subtilisin|P" (uppercase P, consistent with proteomics convention for the proline restriction notation in the cleavage column); the PR had "subtilisin|p". Adopting upstream's uppercase form, including the matching ProteinDigestionTests assertions (upstream also adds LoadProteaseDictionary_SubtilisinP_HasProlineRestrictedMotifs). - Newer-TopFD Ms1Feature test data and tests. The PR named its data file Ms1Feature_TopFDvLatest_ms1.feature; upstream named the identical content Ms1Feature_TopFDv1.7.0_ms1.feature. Adopting upstream's version-specific naming and removing the redundant TopFDvLatest copy. Conflict resolution: took upstream's version for proteases.tsv, Test.csproj, Test/FileReadingTests/ExternalFileReading/TestMsFeature.cs, and Test/ProteomicsTests/ProteolyticDigestion/ProteinDigestionTests.cs. Deleted Test/FileReadingTests/ExternalFileTypes/Ms1Feature_TopFDvLatest_ms1.feature. Build clean. The 80 tests in the affected fixtures (TestMsFeature, TestFromFileDeconvolution, TestMs1FeatureFromMassFeature, TestConsensusToMs1FeatureEndToEnd, ProteinDigestionTests) all pass. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
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>
[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>
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>
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 smith-chem-wisc#1069. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
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 smith-chem-wisc#1069. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
… 229 -> bug is in our code DumpDensestScanPeaks writes the densest MS1 scan's 19,184 centroid peaks for the new flashdeconv_snip OpenMS extract (difftest/avg_snip/flashdeconv_cpp.cpp, 2nd CMake target linking the built OpenMS) which runs the real FLASHDeconvAlgorithm::performSpectrumDeconvolution exactly as the TOPP tool does for MS1. With OUR exact params (min_charge=1, max_charge=60, tol=10, min_mass=50, max_mass=100000, min_isotope_cosine=0.85) OpenMS produces 17 masses on scan 2301 -- the true masses (1549.2, 3741.3, 3894.2, 4389.6, 4698.7, 4788.8, 5330.9, 5574.5, 6585.4, 6970.7, 6971.7, 7907.9, 8040.4 ~= the GT msalign) -- while our pipeline produces 229 and recovers only 3 of them. Same input, same params, same (separately confirmed identical) averagine. So the over-production is DEFINITIVELY a bug in our code, not averagine/params/GT. OpenMS also emits single-charge groups (z3-3, z4-4, z8-8) so "single-charge" was a red herring; we form WRONG monos (miss 3894.2@z4 cos 0.970 that OpenMS finds) and pile up halos. fd_cpp_z1_60.txt is the 17-mass ground truth to debug against. Prime suspect: scoreAndFilterPeakGroups_ pre-loop offset seeding our loop omits. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
… truths, scoring drops them DenseScan_OursVsOpenMSTruth reads OpenMS's 17 ground-truth masses (fd_cpp_z1_60.txt) and, for each, reports our nearest CANDIDATE and our nearest FINAL group. Result: our candidate generation forms the true masses (3894.2@z4 EXACT, plus 1549.2, 4698.7, 4788.8, 3741.3, ...; a few off by ~1 Da / one charge), but our FINAL output drops most of them (ourFinal NONE for 1549.2, 3894.2, 5330.9, 5435.8, 6970.7, 7907.9, 8040.4, 12112.3) while keeping 229 halos. So the over-production bug is entirely in our SCORING pipeline (10-iter recruit->updateQscore loop + gates + removeChargeError + removeOverlapping): it discards correct candidates and retains wrong ones. Candidate gen is fine. Next: instrument the real OpenMS scoreAndFilterPeakGroups_ to record per-candidate intermediate state to emulate/debug against. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
… real OpenMS
Instrumented real OpenMS (env-gated FD_TRACE/FD_CANDTRACE/FD_BINTRACE/FD_PEAKTRACE
in difftest/avg_snip/flashdeconv_snip via OpenMSConfig linkage) to dump per-stage
intermediate state for scan 2301. With OUR exact params (min_charge=1, max_charge=60)
OpenMS produces 17 masses; ours was producing 229. Stage counts pinpointed the
divergence to getCandidatePeakGroups_ (post-fix 12157 emit vs 5840 -- still off but
this commit closes the easily-identified transcription gaps). Each fix is confirmed
verbatim against the OpenMS source:
1. Candidate MinAbsCharge = 1 (was crLo+1) -- OpenMS line 678:
PeakGroup pg(1, per_mass_abs_charge_ranges(1,bin)+1, is_positive_)
hardcodes minAbsCharge to 1. Using crLo+1 collapsed candidates to narrow/single-
charge (85% span0) where OpenMS produces mostly wide (63% span5+) -- flipping it
the right way. Scoring then narrows via updateChargeRange_.
2. Forward-final harmonic add = RAW h_max (was Math.Min(max_peak, h_max)) --
OpenMS line 842: total_harmonic_intensity[l] += h_max_isotope_intensity[l]
3. Backward-final harmonic add = Math.Min(max_peak, h_max) (was RAW) --
OpenMS line 925: total_harmonic_intensity[l] += std::min(max_peak_intensity, h_max_isotope_intensity[l])
(The forward/backward "final" adds were swapped vs OpenMS: forward uses RAW,
backward uses MIN; the in-loop adds use the opposite. Our port had them flipped.)
4. Added signal-peak bin-validity check `sbin < massBins.Length` in both forward
and backward signal branches (OpenMS lines 783, 864 bin < mass_bin_size).
Currently a no-op on this scan but faithful.
Effect on scan 2301 (z1-60): final masses 229 -> 165 (target 17); candidate
charge-span flipped from 85% single-charge to span5+=7652; harmonic gate rejection
13350 -> 13686 (target 19682). Remaining gap localized to an m/z-bin-boundary
peak-collection difference (we collect an extra iso-0 peak 602.1932 at bin 601
charge1 that OpenMS's walk excludes) -- to be debugged in a follow-up.
Note: 2 CandidateGen differential tests (CandidateGen_MatchesOpenMS,
CandidateGen_DenseRealScan_MatchesOpenMS) now fail because they diff against
cand_cpp (a hand-transcribed snip that still has the pre-fix logic). The real
OpenMS link snip (flashdeconv_snip) supersedes it; either update cand_cpp to
match these fixes or retire it. Full unit suite otherwise green (the 42 KoinaProsit/
PFly failures are unrelated network-dependent tests).
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…Groups -> 12478 -> 5822 candidates (matches OpenMS 5840)
THE candidate-gen bug. OpenMS FLASHDeconvAlgorithm::updateMembers_ (cpp:170-175)
narrows the tolerance fed into the finder bins + getCandidatePeakGroups_:
for (double& j : tolerance_) { j *= 1e-6; j /= tol_div_factor; }
So tolerance_[ms_level-1] = 10ppm / 2.5 = 4e-6 (NOT 1e-5). That tol drives the
mz_delta tolerance in the peak walk (cpp:759 mz_delta = min(0.16/abs_charge,
2*tol*mz)), so mz_delta = ~0.0048 for charge 1 (vs ~0.012 with un-narrowed tol).
Our FindCandidates was passing the UN-narrowed `p.DeconvolutionTolerancePpm * 1e-6`
(= 1e-5) to GetCandidatePeakGroups. That made our mz_delta 2.5x too wide, so
borderline peaks ~7 mDa off the isotope grid (e.g. 602.1932 at iso0 for the
candidate at 601) fell in our signal branch when OpenMS's tighter mz_delta puts
them in the harmonic branch. Result: our total_signal was inflated -> harmonic
gate too lenient -> 12478 candidates emitted vs OpenMS's 5840.
Fix: pass `p.DeconvolutionTolerancePpm * 1e-6 / p.TolDivFactor` (= 4e-6) to
GetCandidatePeakGroups, matching OpenMS exactly.
Verified on scan 2301 (z1-60):
candidate count: 12478 -> 5822 (OpenMS 5840; within 0.3%)
ours-within-10ppm-of-OpenMS: 4961
Stage counts (ours): 5822 -> 158 (scoring gates) -> 145 (removeCE) -> 100 (overlap)
vs OpenMS: 5840 -> 32 -> 29 -> 17.
Candidate generation now essentially matches OpenMS. Remaining ~5x overproduction
is in the SCORING pipeline (next target).
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…> 16 final masses on scan 2301
Same tol-narrowing bug as the candidate-gen fix (parent commit), now in the
scoring loop. ScoreCandidatesViaPeakGroups was passing the un-narrowed
p.DeconvolutionTolerancePpm * 1e-6 (= 1e-5) as `tolFraction` to
RecruitAllPeaksInSpectrum, while OpenMS scoreAndFilterPeakGroups_ (cpp:1133)
uses `tolerance_[ms_level-1]` which updateMembers_ has already narrowed by
tol_div_factor (= 4e-6).
Effect: our recruit was gathering 2.5x too many borderline-isotope peaks per
candidate, building mis-scored envelopes that survived the 0.85 cosine / SNR /
Qscore gates.
Fix: tolFraction = p.DeconvolutionTolerancePpm * 1e-6 / p.TolDivFactor.
Effect on scan 2301 (z1-60), OURS vs OpenMS ground truth (from instrumented
flashdeconv_snip):
Stage counts:
OURS OpenMS
candidates 5822 5840
afterScoringGates 37 32
afterRemoveChargeError 29 29 <- exact match
afterRemoveOverlap 16 17 <- within 1
Truth recovery (DenseScan_OursVsOpenMSTruth): 11/17 of OpenMS's truths now
recovered with the EXACT charge ranges, vs 3/17 before this session:
3741.3 z3, 3894.2 z4, 4389.6 z6, 4698.7 z4, 4788.8 z4, 5435.8 z8,
6585.4 z8, 7576.6 z7-8, 7907.9 z8, 12112.3 z12-15, 27889.7 z23-36.
Missed 6: 1549.2 z1, 5330.9 z5, 5574.5 z5, 6970.7 z6-8, 6971.7 z6-9, 8040.4 z7-10.
Per-spectrum over-production essentially solved (229 -> 16; OpenMS 17). Remaining
work: recover the 6 missed truths.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…nMS faithful)
OpenMS scoreAndFilterPeakGroups_ (cpp:1126) only drops a peak group if
`peak_group.empty()`. There is NO `>= MinIsotopicPeakCount` gate at this stage.
Our prior `if (pg.SignalPeaks.Count < p.MinIsotopicPeakCount /*=3*/) continue;`
was an unfaithful strictness that, at the now-tightened tolerance, was dropping
real low-mass / few-isotope candidates (e.g. 1549.2 z1) that OpenMS retains.
Result on scan 2301 (z1-60):
OURS OpenMS
candidates 5822 5840
afterScoringGates 39 32
afterRemoveChargeError 31 29
afterRemoveOverlap 17 17 <- EXACT MATCH
Truth recovery 12/17 (vs 11/17 before): now picks up 1549.2 z1-1 (cos 0.912).
The synthetic placeholder test MetaFlashDecon_ReturnsScoredEnvelopesOnRealSpectrum
(broken since tightening tol) is green again.
5 truths still missed (5330.9, 5574.5, 6970.7, 6971.7, 8040.4) — all have
candidates at the right mass region but with monos shifted 0.3-1.7 Da; suspect
the 10-iter recruit/updateQscore loop converging to a wrong isotope offset (no
pre-loop seed from getIsotopeCosineAndDetermineIsotopeIndex on the candidate's
per_isotope_intensities, FLASHDeconvAlgorithm.cpp:1091-1093).
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
… seed (OpenMS faithful)
OpenMS scoreAndFilterPeakGroups_ (cpp:1090-1101) does an initial cosine-vs-averagine
alignment BEFORE the 10-iter recruit->updateQscore loop, to pick a starting offset
that corrects ±isotope mis-anchors in the candidate's mono mass:
float prev_cos = peak_group.getIsotopeCosine();
float cos = getIsotopeCosineAndDetermineIsotopeIndex(monoMass, peak_group.getIsotopeIntensities(),
offset, avg, -minNegIsoIdx, -1, allowed_iso_error_, ...);
prev_mono_mass = monoMass + offset * iso_da_distance_;
peak_group.setIsotopeCosine(cos);
if (cos < min(0.5, min_isotope_cosine) - 0.1) continue;
Our ScoreCandidatesViaPeakGroups skipped this step (started at offset=0), so the
loop could stabilise at the wrong isotope for candidates whose getCandidatePeakGroups_
mono came out off by +/-1 Da.
Implementation:
- Added PerIsotopeIntensities to MetaFlashDeconAlgorithm.CandidateMass (shape =
MetaFlashDeconPeakGroup.PerIsotopeInt: indexed by `anchored_iso - MinNegativeIsotopeIndex`).
- MetaFlashDeconCandidateFinder.GetCandidatePeakGroups now computes & stores this
vector right after the apex-anchored mono walk (uses the same peaks list +
isoTol filter that builds the mono mass).
- ScoreCandidatesViaPeakGroups adds the OpenMS pre-loop block: calls
GetIsotopeCosineAndDetermineIsotopeIndex(perIso, b, apex, isoIntShift=1,
windowWidth=-1, minIsoSize=2 [matches OpenMS const min_iso_size_=2],
out int offset), applies the cos<min(0.5,minCos)-0.1 pre-filter, and sets
prevMono = seedMono + offset*isoDa (matches cpp:1093).
Result on scan 2301 (z1-60):
OURS OpenMS
candidates 5822 5840
afterScoringGates 40 32
afterRemoveChargeError 34 29
afterRemoveOverlap 22 17
Truth recovery 14/17 (up from 12/17): now recovers 5574.5 z5 (cos 0.864) and
8040.4 z7-10 (cos 0.877). Still missed (3): 5330.9 z5, 6970.7 z6-8, 6971.7 z6-9
(candidates exist with monos shifted 0.3-1.5 Da; suspect convergence loop
detail or removeChargeError/overlap subtlety).
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…C13-C12 (1.0033548) -> exact 5840=5840, 32=32
OpenMS Constants::ISOTOPE_MASSDIFF_55K_U = 1.002371 (Constants.h:100) -- the
average isotope mass difference for a 55 kDa peptide, used as `iso_da_distance_`
throughout FLASHDeconv (PeakGroup.h:315; FLASHDeconvAlgorithm::performSpectrumDeconvolution
sets it from this constant at cpp:76-78). Ours was using mzLib's
Chemistry.Constants.C13MinusC12 = 1.0033548 (the pure 13C-12C diff).
The ~0.001 Da/isotope drift was small but flipped borderline isotope assignments
in dense regions and pushed mono masses off by a fraction of a Da at higher
masses (compounding through the 10-iter loop). Centralised the OpenMS value as
`MetaFlashDeconAlgorithm.IsoDaDistance55K` and switched all 3 files
(MetaFlashDeconAlgorithm, MetaFlashDeconCandidateFinder, MetaFlashDeconPeakGroup)
to use it instead of Constants.C13MinusC12.
Result on scan 2301 (z1-60) -- now matches OpenMS at the early stages exactly:
OURS OpenMS
candidates 5840 5840 <- exact match
afterScoringGates 32 32 <- exact match
afterRemoveChargeError 24 29
afterRemoveOverlap 16 17
Truth recovery 15/17 (up from 14/17): now recovers 5330.9 z5 (cos 0.854).
Remaining 2 missed: 6970.7 z6-8 and 6971.7 z6-9 (a 1-Da-apart pair).
For DenseScan_OursVsOpenMSTruth, ourCand for 6585.4 now lands at 6585.4 z1-8
(was 6586.3 z1-9) and ourCand for 5330.9 at 5329.9 z1-5 (was 5332.6 z1-7) --
the more accurate iso_delta shifts candidate monos closer to truth.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…r-group CE overlap-ratio trace
Added to MetaFlashDeconVsRealPerSpectrum (both DenseScan_OursVsOpenMSTruth and
DenseScan_ClassifyExtrasVsReal):
- Dump of every entry in our `scored` (post-gates, n=32) and `afterCE` (n=24)
next to ourFinal (16), for direct comparison against the OpenMS FD_TRACE
AFTERCE dump (now also adds AFTERCE per-group lines in OpenMS).
- PER-GROUP table in the CE-instrument block listing (mono, z-range, rep,
overlap_intensity, GetIntensity, ratio, DROP/keep) so we can see exactly
which groups exceed the 0.5 threshold and by how much.
- TRACE block for one named target mass (currently 6970.7) listing every
shared-peak contribution to its overlap_intensity (pmz, pint, j, mass_j,
repz comparison, snr-comparison verdict).
Used these to localise the 2 missed truths on scan 2301 (6970.7 z6-8 and
6971.7 z6-9). Diagnosis: borderline overlap-ratio drops (0.524 / 0.532) from
accidental m/z collisions between mid-mass groups' early isotopes and the
27889.7 z23-36 protein's DEEP isotopes (iso-17 of 27889@z32 = m/z 872.84 = iso-2
of 6970@z8, etc.). 27889 has chargeSnr ~8.78 at z32 vs 6970's ~0.91 at z8 ->
6970 loses the snr-comparison and accumulates overlap that just barely exceeds
the 0.5 drop threshold. OpenMS does the same computation but keeps these groups
(snr borderline must fall the other way for it) -- gap is small (~5%) and the
remaining investigation requires instrumenting OpenMS's per-charge SNRs.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
… extra peak in 27889@z32 Instrumented OpenMS scoreAndFilterPeakGroups_ to dump per-charge SNRs alongside SCORED entries + the per-peak (mz, intensity, isotopeIndex) content of the 27889 z32 group. Mirrored on our side in DenseScan_ClassifyExtrasVsReal. Side-by-side comparison (scan 2301, mass 27889.7 z32): Ours: 17 z32 signal peaks OpenMS: 16 z32 signal peaks diff: ONE EXTRA peak at iso 25, m/z 873.3480, int 3.52e4 (only in ours) That single extra peak inflates our PerChargeInt[32] and PerChargeSumSignalSq[32] just enough to push chargeSnr_27889(z32) from OpenMS's 7.33 to our 8.78 (20% higher). Per-charge SNRs match at every other charge of 27889 (z23-z36) and the averagine for mass 27889 matches OpenMS exactly (apex 17/17, left 10/10, right 12/12, b-cos 1.0000). The extra peak triggers a cascade: the 27889@z32 SNR inflation slightly shifts the snr-comparison budget in removeChargeError, just enough to flip the overlap-intensity ratio for 6970.7's 4 dups + 6971.7 from ~0.49 to 0.52-0.58 (crossing the 0.5 drop threshold). Result: we drop the entire 6970-6971 cluster where OpenMS keeps it. Mechanism of the extra peak (signal-check threshold): pmz = 873.3480, cmz_z32 = 872.5604, iso_delta = 0.031323 iso_index = round((873.348 - 872.5604) / 0.031323) = 25 error = |0.7876 - 25*0.031323| = 0.00453 threshold = pmz * tol = 873.348 * 4e-6 = 0.00349 Nominal: 0.00453 > 0.00349 -> should NOT be signal (OpenMS correctly excludes). Our recruit nevertheless adds it. Suspect: our pg.MonoisotopicMass after the 10-iter loop differs from OpenMS's by ~0.04 Da at this position -- enough to shift cmz and bring 873.3480's iso-25 error just under threshold. This commit only extends the test diagnostics; the fix requires either matching the convergence behaviour of OpenMS's 10-iter loop more precisely (probably updateMonoMassAndIsotopeIntensities or updatePerChargeInformation precision) or accepting the borderline 1-mass delta on dense scans. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…-25 borderline -> 2 missed truths
Instrumented OpenMS PeakGroup::recruitAllPeaksInSpectrum (gated FD_RECRUITTRACE,
+include <fstream>/<cstdlib>) to dump the per-iteration signal-check computation
for peak m/z 873.348 at charge 32 of the 27889 group. Our C# test side dumps
the same. Side-by-side:
OpenMS 10-iter loop hits m/z 873.348 three times (3 different mono values):
iter A: mono=27885.7258 -> iso=29, err=0.00341 < threshold 0.00349 -> SIGNAL
iter B: mono=27890.7279 -> iso=24, err=0.00372 > threshold 0.00349 -> NOT signal
iter C: mono=27889.7321 -> iso=25, err=0.003516 > threshold 0.003493 -> NOT signal (just over by 2e-5)
Our final iter: mono=27889.7368 -> iso=25, err=0.003345 < threshold 0.003493 -> SIGNAL
So OpenMS's final-iter mono (27889.7321) lands the iso-25 check ~2e-5 above
threshold; ours (27889.7368) lands it ~1e-4 below. The 0.005 Da convergence
drift between the two pipelines' 10-iter loops is the entire cause of the 2
missed truths cascade.
Root of the convergence drift is NOT:
- ProtonMass: mzLib 1.007276466879 vs OpenMS PROTON_MASS_U 1.0072764667710
(diff 1e-10, negligible).
- iso_da_distance: both 1.002371 exactly (fixed earlier in 3badfcb).
- averagine: apex/left/right/b cosine = 1.0000 at mass 27889 (matches OpenMS).
- Pre-loop offset: matches OpenMS shape (cpp:1091-1093, min_iso_size_=2).
Suspect: subtle precision/order in updateMonoMassAndIsotopeIntensities (the
intensity-weighted mass formula or how it accumulates across iterations through
the recruit-then-updateQscore loop). Fixing it would close the last 2 missed
truths but requires per-iter trace dumps from BOTH pipelines (updateQscore
result + mono after each iter).
Best result remains: candidates 5840=5840, gates 32=32, removeCE 24 vs 29,
removeOL 16 vs 17, truth recovery 15/17. Investigation committed for the future.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…T scoring loop)
Instrumented BOTH 10-iter scoring loops to dump per-iteration recruit-mono,
returned offset, and mono-after-updateQscore for a target mass window (gated by
FD_ITERTRACE / FD_ITERTRACE_CS, target 27880-27900):
Cand A (start ~27885):
OpenMS iter 0: recruit=27885.72584 offset=0 mono_after=27885.73788
Ours iter 0: recruit=27885.72160 offset=0 mono_after=27885.73111
Cand B (start ~27890):
OpenMS iter 0: recruit=27890.72791 offset=-1 mono_after=27890.73444
Ours iter 0: recruit=27890.73700 offset=-1 mono_after=27890.73747
OpenMS iter 1: recruit=27889.73207 offset= 0 mono_after=27889.74032
Ours iter 1: recruit=27889.73509 offset= 0 mono_after=27889.73683
Both pipelines: (a) same number of iterations, (b) same offset sequence
(0 for A; -1, 0 for B), (c) loop converges in the same step.
But the iter 0 recruit_mono ALREADY differs by 0.004 (cand A) / 0.009 (cand B).
That's `candidate.Mass + offset_preloop * isoDa`. With offset_preloop=0 for cand A,
our candidate.Mass = 27885.722 vs OpenMS 27885.726 (CAND line). With
offset_preloop=1 for cand B, our candidate.Mass = 27889.735 vs OpenMS 27889.726.
So the entire 0.005 Da convergence drift is INHERITED from candidate generation
-- the candidate's intensity-weighted mono in `getCandidatePeakGroups` is what
differs. The 10-iter scoring loop just propagates it. Looking at the algorithm
(nominator += pi * (um - iso*isoDa); intensity += pi; mono = nominator/intensity)
the candidates have identical PEAK COUNT (verified 5840=5840) and same averagine,
threshold, and isoDa -- yet the per-candidate mono diverges 4-9 mDa. Suspect a
precision issue in either:
(a) `um = (pk.mz - ProtonMass) * pk.charge` -- OpenMS uses (float)PROTON_MASS_U
vs our double Constants.ProtonMass; per-peak diff ~1.8e-6 * z, but for
z23-36 of a 150-peak group this could cumulate.
(b) max-intensity peak (tMass) selection -- OpenMS Peak1D::intensity_ is float
while our peaks list stores double intensity; ties could break differently
and pick a different reference peak -> different anchored iso indices ->
different mono.
Reverting the float-intensity casts in updateMonoMassAndIsotopeIntensities +
updatePerChargeInformation (tried; no effect). The drift is at candidate-gen
stage, not scoring stage. Next investigation point.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…sses -> fixes +2 mass-bin shift
Found the root cause of the 134-bin discrepancy and the entire 2-missed-truths
cascade: our `MetaFlashDeconAveragine.GetEntry` computed `avgDelta = (weighted
mean K) * iso_da_distance`, treating the +k peak as mono + k*1.002371. But the
ACTUAL C-13 isotope shift is ~1.00335 (and varies per element), so for low-mass
peptides where the +1 peak is dominated by C-13, our `avgDelta(50)` was 0.026 vs
OpenMS's 0.025 (0.0004 Da bias).
That 0.0004 Da bias propagates into `mass_bin_min_value` (cpp:1074 / our
FindCandidates line 77: `log(max(1, min_mass - avgDelta(min_mass)))`), which
shifts massBinMin by ~8e-6 in log space => +2 in every mass bin index across the
entire candidate-gen pipeline. Same accumulated mass_intensities VALUES but
indexed 2 bins higher than OpenMS. This cascaded into filterMassBins top-3
running-max selecting different charges per m/z bin for the 27909 mass bin, then
fewer peaks in candidate pg (146 vs 151), then 0.005 Da drift in candidate mono
(27889.737 vs 27889.726), then borderline iso-25 inclusion at z32 of 27889
flipping, then 27889@z32 SNR 8.78 vs 7.33, then 6970/6971 borderline overlap.
Fix: track `intens[k] * masses[k]` per re-binned iso slot, apply the same Trim+
Normalize, then compute `avgDelta = sum_kept(intens*mass) / sum_kept(intens) -
mono` matching OpenMS `iso.averageMass() - iso[0].getMZ()` (cpp:107). This uses
the ACTUAL coarse-isotope masses (per-element averaged), not k * iso_da_distance.
Result at mass 50: ours 0.025 vs OpenMS 0.025 (was 0.026, now matches);
mass 1500: ours 0.952 vs 0.953 (was 0.979).
Scan 2301 stage counts:
OURS OpenMS
candidates 5796 5840
afterScoringGates 35 32
afterRemoveCE 27 29
afterRemoveOverlap 18 17
Truth recovery 15/17 (same count, different masses recovered): now picks up
7576.6, 8040.4, 12112.3, 27889.7 (high-mass) that the +2 bin shift had been
preserving accidentally. Still missing 6970.7, 6971.7 (the original +1 isotope
overlap-borderline pair; their fix likely involves a different downstream
threshold).
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…ctual masses -> fixes +2 mass-bin shift" This reverts commit 1ec5383.
…9-bin region (uncovered the +2 bin-shift -> avgDelta(50) chain)
…actual masses -> fixes +2 mass-bin shift" This reverts commit 548a242.
…iso-bin actual masses -> fixes +2 mass-bin shift"" This reverts commit ff93c10.
OpenMS PeakGroup::updateMonoMassAndIsotopeIntensities (PeakGroup.cpp:706-736) sorts peaks by logMz then sums the intensity-weighted mono with a FLOAT denominator (float pi, float intensity_). Our candidate-gen and scoring mono both summed in double and in recruit order, drifting the 27889 protein mono ~5-9 mDa. That tipped one borderline z32 signal peak (m/z 873.348, iso 25) into the envelope; it is shared with 6970.7/6971.7, and the extra shared intensity pushed their removeChargeError overlap ratio just over 0.5, dropping both true masses. Match OpenMS exactly: sort by logMz then float intensity, float pi, float denominator, exclude iso<0 from the mono. Densest real scan 2301 now matches OpenMS at every stage: candidates 5840=5840, gates 32=32, removeChargeError 29=29, final 17=17, truth recovery 17/17 (was 15/17). Relax the RealData peak-count invariant from MinIsotopicPeakCount(3) to 2: OpenMS drops only empty groups (cpp:1126) and a candidate needs >=2 isotopes, so a 2-peak envelope is faithful (commit 8af4e4d). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
ComparePerSpectrumToRealFlashDeconv now samples every 20th MS1 scan (~160 across the run) instead of every 80th (~40), for a tighter whole-run fidelity readout against the real _ms1.msalign. Every-scan (~3200) is impractical (hours through the full candidate-gen + scoring pipeline). Whole-run result at this density: matched 6605/6767 = 97.6% (15 ppm), our/real 1.04x, low-mass recovery complete. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add "do not change" notes at the fragile, differential-tested details that an innocent refactor would silently break - distilled from the STATUS investigation log: - IsoDaDistance55K = 1.002371 (ISOTOPE_MASSDIFF_55K_U), NOT Constants.C13MinusC12 - intensity-weighted mono: FLOAT pi + FLOAT denominator + logMz-sort + iso<0 excluded, identical in candidate-gen and scoring (the float/order detail behind 17/17) - GetBinNumber = floor(x+0.5), not Math.Round (banker's); binOffsets need AwayFromZero - scoring tolerance narrowed by tol_div_factor; empty-only post-gate (not >=3) - removeChargeError 0.5 overlap threshold + per-charge-SNR >2x comparison - getNoisePeakPower_ rounding must be AwayFromZero (caught 2 real bugs) - mass_bin_min/avgDelta +2-bin shift is a known, deliberate accidental-clean state - parameter defaults are literal OpenMS values, differential-tested, not free knobs - LogMzPeak.Intensity is float-origin; constant-averagine ctor is test-only Comments only; full suite 4041 pass / 0 fail / 5 skip. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
OpenMS MassFeatureTrace::findFeatures (MassFeatureTrace.cpp:84-143) drops every
mass trace whose summed per-isotope pattern has averagine cosine < 0.75. Our
independently-written native tracer never had this ("Phase 3b", deferred), so it
over-produced features ~4.7x (2845 vs the real FLASHDeconv 603 on rep2 fract10).
Port it verbatim:
- plumb each envelope's per-isotope vector (PeakGroup::getIsotopeIntensities,
min-negative-isotope based) through IsotopicEnvelope to the tracer
- per trace, accumulate member per-isotope intensities in FLOAT (std::vector<float>),
shifted by iso_off = (int)(0.5 + (mono-centroid)/isoDa) [truncate-toward-zero == C++ int()]
- cosine via the differential-tested getIsotopeCosineAndDetermineIsotopeIndex with
the exact MassFeatureTrace args (iso_int_shift=0, window_width=0, min_iso_size=2),
threshold 0.75 (MassFeatureTrace.cpp:32 default, NOT 0.85)
Differential-tested against the real OpenMS static function (massfeature_snip,
links built OpenMS): on 38 real traces, offsets identical, masses within 5e-6,
cosine within 2e-3 (averagine boundary), and zero traces flip the 0.75 gate.
Whole-file (rep2 fract10): 2845 -> 608 features vs real 603 (4.7x -> 1.008x).
Filter applies whenever per-isotope data is present (always true for real
MetaFlashDecon envelopes, as in OpenMS); skipped only for synthetic envelopes
that carry none. Full suite 4041 pass / 0 fail / 5 skip.
NOTE: feature-SET correspondence is still ~79% recall / ~57% precision vs OpenMS
because the trace DETECTION (our Kenar-style extender) is not yet the faithful
OpenMS MassTraceDetection (next: "piece C").
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Replace the independently-written Kenar-style trace extender with a verbatim port of OpenMS MassTraceDetection::run_ (MassTraceDetection.cpp:449-663) — the algorithm MassFeatureTrace actually runs (MassFeatureTrace.cpp:68-75): - apices = every collapsed neutral-mass peak, stable-sorted by intensity, seeded most-intense-first - bidirectional findNearest extension within +/-3*ftl_sd (ftl_sd = centroid/1e6 * mass_error_ppm), fixed (reestimate_mt_sd=false) - exact multiplicative iterative weighted-mean centroid (not naive Sum(wv)/Sum(w)) - outlier termination at >trace_termination_outliers consecutive non-empty misses; empty scans advance the index without incrementing the miss counter - peaks marked visited ONLY after isTraceValid_ (min_trace_length + min_sample_rate), so failed traces leave their peaks for other seeds Also correct the params to the real values MassFeatureTrace uses (our independent version had all three wrong): mass_error_ppm 20 (was effectively 10), outliers 5 (was 20), min_sample_rate 0.1 (was 0.05). Differential-tested vs the REAL OpenMS MassTraceDetection (masstrace_snip, links built OpenMS) on identical serialized per-scan peaks: 24/24 traces match both ways (centroid 10 ppm + RT overlap + identical size). So the port is faithful; the remaining whole-file feature-set divergence vs OpenMS is the upstream per-spectrum decon residual (97.6%), not the tracer. Whole-file (rep2 fract10): 488 features vs real 603. Full suite 4041 pass / 0 fail / 5 skip; tracer unit tests 7/7. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…eal data OpenMS MassFeatureTrace keys one Peak1D per PeakGroup with no ppm merging. Our per-scan collapse has a within-ppm merge, but MetaFlashDecon emits one envelope per neutral mass per scan and RemoveOverlapping keeps survivors apart, so the merge is dormant on the real path. Confirmed empirically (Collapse_MergeStats): 0/912 envelopes merged over 60 real scans (0.00%) -> already 1:1 and faithful. The merge only fires for non-faithful multi-same-mass-per-charge inputs (synthetic tests). Documented; left as-is (forcing strict 1:1 would break the synthetic charge-collapse test for zero real-path gain). Remaining whole-file feature-set gap vs OpenMS is the upstream per-spectrum decon residual (97.6%), not the tracer. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…gate/pg traces Adds [Explicit] harnesses to localize per-spectrum divergence from real OpenMS at identical params (the production _ms1.msalign uses different params, so it's not the fidelity reference): - DumpScansForSnip / CompareScansVsSnip: dump a sparse/medium/dense scan's peaks and compare our decon to flashdeconv_snip (real OpenMS). Result: dense/medium near-exact (scan 3551 135/135, 2625 26/27), but SPARSE scan 1332 = 15/20 (miss 5, extra 0). - DiagnoseScan1332Misses: stage-by-stage (cand/scored/afterCE/final) per OpenMS truth. - env-gated FD_GATETRACE_CS / FD_PGDUMP_CS in ScoreCandidatesViaPeakGroups (dormant unless set) dump per-candidate gate decisions + per-charge recruited peaks. ROOT CAUSE of the sparse-scan misses (confirmed, not guessed): for masses like 8166.8, OpenMS builds a clean z9-10 envelope (cos 0.893); ours forms fragmented, wrong-charge-range groups (z9-11 with a malformed z9 [iso 6-7 only] + spurious z11), cos 0.81 < the 0.85 gate -> UpdateQscore returns 0 (rep/snr/q zeroed) -> dropped. This per-spectrum loss disproportionately hurts FEATURE yield: a mass scoring just under 0.85 on a feature's weaker scans breaks the trace. Next: differential-test the recruit / updateChargeRange_ on scan 1332 vs OpenMS and copy the divergent step (rather than reasoning about why ours fragments). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…arse-scan masses The bin-axis origin mass_bin_min = log(min_mass - avgDelta(min_mass)) (OpenMS FLASHDeconvAlgorithm.cpp:1074) was using mzLib's Averagine value (~0.0251-0.0255) instead of OpenMS's CoarseIsotopePatternGenerator value (0.02514582). That ~mDa error shifts the whole mass-bin axis ~2 bins vs OpenMS. On DENSE scans it was a harmless near-integer shift, but on SPARSE scans it lands candidates in the wrong filterMassBins charge-range slot -> a spurious extra charge -> envelope cosine < 0.85 -> dropped. That silently cut per-spectrum recall, which disproportionately kills FEATURE yield (a mass scoring just under 0.85 on a feature's weaker scans breaks its trace). Copy-verified the full chain to this root via the OpenMS-linked snips: the dropped mass 8166.8 on scan 1332 -> spurious z11 -> candidate charge range z1-11 vs OpenMS z1-10 -> the 8170 bin sits ~4 mDa off -> the avgDelta axis. Recruit (z9/z10 peaks identical to OpenMS) and UpdateChargeRange are verified-faithful; the only divergence was the axis. Fix: use the exact-OpenMS avgDelta at the FLASHDeconv min_mass (50) for mass_bin_min (GetMassBinMinAvgDelta -> OpenMsAvgDeltaAtMass50 = 0.025145816488311823, from the OpenMS-linked avg_snip / FD_SCALARTRACE; mzLib cannot reproduce it without OpenMS's coarse isotope generator + abundances). avgDelta is used ONLY here. Result (vs real OpenMS at our params, flashdeconv_snip): sparse scan 1332 15/20 -> 20/20 (all 5 misses recovered, including 2 that previously formed NO candidate); dense scan 2301 preserved 17/17 (exact axis = zero shift, unlike a partial avgDelta fix which gave a bad 0.3-bin shift and regressed to 13/17). Full suite 4041 pass. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…our guesses Faithfulness audit (vs topp/FLASHDeconv.cpp, the tool that made the ground-truth _ms1.feature/_ms1.msalign) found four params where we'd diverged: - decon max_charge: 60 -> 100. FLASHDeconv TOPP default is 100 (FLASHDeconv.cpp:180); the production _ms1.msalign contains charges up to 100. At 60 we silently dropped every species above z60. - feature-tracer min_isotope_cosine: 0.75 -> 0.85. We used the MassFeatureTrace CLASS default (0.75); the TOPP tool sets it to -1 and resolves it to the MS1 min_isotope_cosine = 0.85 (FLASHDeconv.cpp:198, 573-575). - feature-tracer mass_error_ppm: 20 -> 10. We used the MassTraceDetection class default (20); the TOPP tool sets it to -1 and resolves it to the MS1 tol = 10 (FLASHDeconv.cpp:209, 563-565). - feature-tracer min_sample_rate: 0.1 -> 0.05. The TOPP tool overrides the MassFeatureTrace class default 0.1 to 0.05 (FLASHDeconv.cpp:210). min_charge=1, tol=10, min_mass=50, max_mass=100000, decon min_isotope_cosine=0.85, trace_termination_outliers=5 already matched the TOPP-resolved values. Goal is a faithful FLASHDeconv, not an imitation: match the params the real tool runs with. Updated the stale default-charge-range test (60 -> 100). Full suite green except that one assertion (now fixed). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The native feature tracer used DefaultMaxOutlierScans=5 (the MassTraceDetection class default), with a comment wrongly claiming the TOPP tool leaves it there. In fact FLASHDeconv.cpp:570 hard-sets trace_termination_outliers to 20 at runtime (lines 204-205 only remove it from the user-facing FeatureTracing param tree; the later setValue(...,20) is operative). At 5, a mass trace terminates after 5 consecutive non-empty scans with no in-window peak, so long features whose mass deconvolves elsewhere across a sparse elution tail collapse below the 10 s min_trace_length and are dropped. Worked example 22550.5 (RT 2529-2575 s): present at the apex (scans 1790-1798) and again ~36 s / 18 MS1 scans later (scan 1839); at 5 the trace dies at 9 s and is dropped, at 20 it bridges the gap and survives (45 s). Verified faithful, not guessed: our per-spectrum already matches FLASHDeconv's unfiltered _ms1.tsv mass-for-mass on the affected scans (the SNR/500-cap msalign is the wrong per-scan reference); this param was the only divergence. Whole-file z1-100 vs the real FLASHDeconv _ms1.feature: features 233 -> 374 (FD 366), recall 44% -> 68%, unique-mass precision 95.5%. Full suite green (4041 pass). Also adds null-gated tracer trace-dump diagnostics and [Explicit] whole-file / per-scan-window harnesses used to localize this; ComparePerSpectrum now runs at the faithful z1-100 defaults (the production msalign is z1-100). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ful) Ms1FeatureFile.FromMassFeatures emitted one row per contiguous charge run (ToMs1Features), so a gapped charge set became several rows. That was unfaithful: the real FLASHDeconv _ms1.feature writes exactly one row per mass feature with an aggregated Min/Max charge range (verified -- every multi-row mass in the ground-truth file differs in retention time, never only in charge), and it inflated intensity by writing the full feature.SummedIntensity on every split row (a 3-run feature reported 3x its real intensity). Now emits one row per feature via ToMs1Feature. Re-expanding Min..Max on reload may surface an unobserved intermediate charge, but that matches FLASHDeconv and the charge set is metadata, not used for mass/RT pairing. Updated the gapped-set test to pin the single-row, non-inflated-intensity contract. ToMs1Features kept as a public API. Shared by the consensus tracer (smith-chem-wisc#1069) -- both now emit one row per feature. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Revives this PR as MetaFlashDecon — a faithful C# reimplementation of OpenMS FLASHDeconv for top-down (and bottom-up) deconvolution in mzLib, plus a pluggable mass-feature tracer and an
_ms1.featurewriter that feeds MetaMorpheus FromFile search.What's here
Faithful FLASHDeconv port (
MetaFlashDecon*) — ported function-by-function from the OpenMS source and differential-tested against it:FLASHDeconvAlgorithm(setFilters_,updateMzBins_,updateCandidateMassBins_,filterMassBins_,getCandidatePeakGroups_) — continuous-charge support, the low-charge isotope-presence path, two-level harmonic rejection, and thetol_div_factorvoting/dedup tolerances.PeakGroup/PeakGroupScoring— the per-candidate recruit→Qscore offset loop, pre-loop isotope-cosine offset seed, charge-fit / isotope-cosine / charge-SNR gates,removeChargeErrorPeakGroups_, andremoveOverlappingPeakGroups_.IsotopicEnvelope.Scorecarries the logistic-regression Qscore (FLASHDeconv convention), not the raw averagine cosine.Pluggable feature tracing
IMassFeatureTracerwith (a) a native neutral-mass tracer (OpenMSMassFeatureTrace-style) and (b) a consensus-pipeline adapter (Promote consensus mass-tracing types + add Ms1Feature writer for MassFeature #1069).Ms1FeatureGenerator: whole mzML → per-scan deconvolution → tracer →_ms1.feature(TopFD-superset schema, RT in seconds), parallelized per scan.Supporting utilities
MsAlignWriter(MS1/MS2 msAlign export);SpectralSimilarity.CosineOfAlignedVectors.Validation — reproduces real OpenMS FLASHDeconv
Fidelity was established by differential testing against the real OpenMS library (MSVC extracts that link the built
OpenMS.dll) and against the production_ms1.msalignground truth on a Jurkat top-down fraction.Densest real MS1 scan (scan 2301, 19,184 peaks) — our pipeline vs real OpenMS, identical parameters:
Whole-run per-spectrum fidelity — 157 MS1 scans spread evenly across the full run, compared to the real FLASHDeconv
_ms1.msalign(RT-aligned, identical spectra):Reaching exact stage-level and per-spectrum agreement required matching OpenMS down to numeric detail — e.g. the intensity-weighted monoisotopic-mass summation is float-typed and logMz-sorted to match
PeakGroup::updateMonoMassAndIsotopeIntensities, which is what recovers the last borderline isotope-peak assignments.Tests
[Explicit](they require local Jurkat data and the built-OpenMS difftest harness), so they don't run in CI.