fix(consensus)!: make consensus read downsampling deterministic - #727
Merged
Merged
Conversation
|
Note Reviews pausedUse the following commands to manage reviews:
Use the checkboxes below for quick actions:
No actionable comments were generated in the recent review. 🎉 ℹ️ Recent review info⚙️ Run configurationConfiguration used: Path: .coderabbit.yaml Review profile: ASSERTIVE Plan: Pro Run ID: ⛔ Files ignored due to path filters (2)
📒 Files selected for processing (14)
💤 Files with no reviewable changes (2)
WalkthroughConsensus downsampling now uses deterministic fgbio-compatible read-name hash ranks. Vanilla, CODEC, simplex, and duplex workflows remove seeded RNG behavior, preserve pair semantics, validate zero caps, and add unit and integration coverage. ChangesDeterministic consensus downsampling
Estimated code review effort: 4 (Complex) | ~60 minutes Sequence Diagram(s)sequenceDiagram
participant Command
participant ConsensusCaller
participant fgbio_read_name_rank
participant select_lowest_ranking
participant BAMOutput
Command->>ConsensusCaller: run with read cap
ConsensusCaller->>fgbio_read_name_rank: hash read names
fgbio_read_name_rank-->>ConsensusCaller: signed ranks
ConsensusCaller->>select_lowest_ranking: ranks and cap
select_lowest_ranking-->>ConsensusCaller: retained indices
ConsensusCaller->>BAMOutput: write consensus records
Possibly related issues
Possibly related PRs
Suggested labels: 🚥 Pre-merge checks | ✅ 3✅ Passed checks (3 passed)
Comment |
Member
Author
|
@coderabbitai pause |
✅ Action performedReviews paused. |
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #727 +/- ##
==========================================
+ Coverage 94.03% 94.12% +0.08%
==========================================
Files 178 181 +3
Lines 108899 109415 +516
==========================================
+ Hits 102408 102982 +574
+ Misses 6491 6433 -58 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Downsampling selected which reads to keep by shuffling with a seeded StdRng held on the consensus caller. The RNG is stateful, so the surviving subset for a tag family depended on how many families that caller instance had already processed. fgumi did not have fgbio's thread-count bug: the unified pipeline builds a fresh caller per batch at a fixed batch size, so --threads 1 and --threads 8 agreed. But the no---threads path holds a single caller for the entire run, putting its RNG at a different position for every family, so the two execution modes disagreed under a cap. On a 77k-molecule library, duplex output at --max-reads-per-strand 2 differed on 566 of 77,804 records (0.73%). Uncapped runs always agreed. Rank reads by a Murmur3 hash of the read name and keep the lowest-ranking ones. The retained subset becomes a pure function of the family: independent of caller history, execution mode, and the order reads arrive in. Because both ends of a template share a read name they share a rank, so a template is retained or discarded as a unit -- which also recovers reads the old shuffle lost by splitting templates across ends (simplex --max-reads 2 on that library goes from 726,854 to 734,536 consensus reads, matching fgbio's count). The rank is computed once per read and stored, never inside a sort key closure: Rust's sort_by_key may evaluate its closure more than once per element, so hashing there would re-hash on every comparison. Ranks compare as signed i32, matching fgbio's Scala Int, and the sort is stable so tied ranks keep input order. Selection now lives in one place, caller::select_lowest_ranking, shared by the raw-record and codec paths. With no RNG left on the caller there is no per-instance state for downsampling to depend on, so the seed option and the rand dependency are removed and the determinism becomes structural rather than incidental. Which reads are retained differs from previous releases for any run that sets a cap. There is no way to fix the reproducibility bug while preserving the old selection, since the old selection was the bug. Runs without a cap are byte-identical, verified against a pre-change binary. Ports fulcrumgenomics/fgbio#1166. Once that merges, this also restores byte parity with fgbio under a cap; a parity oracle test is deliberately deferred until then rather than pinned to the 4.1.0 shuffle it would disagree with. Closes #725
A cap of zero empties every strand, so no molecule can produce a consensus, and fgumi exited 0 having written an empty BAM. simplex and codec already validate their equivalents. Only a lower bound is checked, matching fgbio: duplex's --min-reads is a vector with its own ordering rule, so the max >= min comparison the sibling tools make does not apply, and a cap below it is already handled gracefully -- consensus_call returns no consensus rather than panicking. A negative value is rejected by clap, since the option is an Option<usize>, so validate() only has to cover zero; fgbio needs the wider check because its option is an Option[Int]. Ports the second commit of fulcrumgenomics/fgbio#1166. Refs #725
nh13
force-pushed
the
725/nhomer/fix-deterministic-consensus-downsampling
branch
from
August 9, 2026 05:12
cf25828 to
00bba7b
Compare
nh13
marked this pull request as ready for review
August 9, 2026 17:44
Member
Author
|
@coderabbitai review |
✅ Action performedReview finished.
|
This was referenced Aug 9, 2026
Merged
Merged
nh13
added a commit
that referenced
this pull request
Aug 18, 2026
…785) Reads discarded by the `--max-reads` consensus-downsampling cap were counted as used and never reached `--rejects`. Because `raw_reads_used = total_input_reads - filtered_reads` and only a recorded rejection increments `filtered_reads`, a family of 10 capped to 3 reported `raw_reads_used = 10` against a consensus of depth 3, and the 7 discarded reads were dropped silently. Add a `Downsampled` rejection reason and record the discarded reads at each drop site so `raw_reads_used` excludes them: - simplex (`VanillaUmiConsensusCaller::process_group`): `downsample_reads` now returns the discarded records alongside the survivors; they are counted and, when rejects tracking is enabled, routed to the `--rejects` output byte-for-byte in input order. - codec (`CodecConsensusCaller`): the per-strand cap counts the discarded reads, mirroring how the codec caller already counts its other per-read rejections (e.g. MinorityAlignment). Per-read routing to `--rejects` for the codec waits on its reject-mask machinery (#751), so this only fixes the count there. The duplex/codec per-strand `consensus_call` path is unaffected: since #727 it retains the full uncapped source reads (the cap shapes only the consensus bases/quals/depths), so no reads are dropped there. `Downsampled` is fgumi-specific: fgbio has no equivalent metric yet (fulcrumgenomics/fgbio#1166, #1167), so `raw_reads_used` and the new `raw_reads_rejected_for_downsampled` row diverge from fgbio under a cap. The new row is emitted only when non-zero.
This branch was previously deployed
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.
Problem
Consensus downsampling selected which reads to keep by shuffling with a seeded
StdRngheld on the consensus caller:The RNG is stateful, so which reads survive a tag family depends on how many families that caller instance has already downsampled.
fgumi does not have fgbio's thread-count bug. The unified pipeline builds a fresh caller per batch at a fixed batch size, so
--threads 1and--threads 8already agreed — I verified this byte-for-byte. But the no---threadspath holds a single caller for the entire run (simplex.rs:378,duplex.rs:420), putting its RNG at a different position for every family. The two execution modes therefore disagreed under a cap:--max-reads-per-strand 2: 566 of 77,804 consensus records differed (0.73%) between--threads Nand no---threads--max-reads 2: differedTwo further consequences of RNG-based selection:
consensus_callonce per strand per end (duplex_caller.rs:1992-1995), so R1 and R2 sampled different templates. Beyond reproducibility, this silently lost reads: families whose split left one end empty produced no consensus pair at all.StdRng(ChaCha); fgbio usesjava.util.Random(an LCG). Against fgbio 4.1.0 at--max-reads-per-strand 2, 487 of 77,804 duplex records (0.63%) differed in SEQ/QUAL. Depth distributions matched exactly, so selection was the only difference. This was already documented as a known limitation indownsample_source_reads.Fix
Rank reads by a Murmur3 hash of the read name and keep the lowest-ranking ones, porting fulcrumgenomics/fgbio#1166. The retained subset becomes a pure function of the family — independent of caller history, execution mode, and arrival order — and because mates share a read name they share a rank, so a template is retained or discarded as a unit.
fgumi already contained a byte-exact port of htsjdk's
Murmur3.hashUnencodedChars, validated against reference vectors captured from htsjdk 3.1.2. It moves from the root crate down intofgumi-raw-bam, where the consensus callers can reach it, alongside the samtools-parity read-name comparator already living there.Three details are load-bearing for parity, and each has a dedicated test:
SourceRead/ClippedRecordInfo, never inside a sort key closure. Rust'ssort_by_keymay evaluate its closure more than once per element, so hashing there would re-hash on every comparison — the same defect the review of fgbio#1166 flagged in its Scala (25–30× the element count at n=100k).i32, matching fgbio's ScalaInt. About half of all read names hash negative, so an unsigned comparison silently reorders them.sortInPlaceBy, so tied ranks keep input order.Selection now lives in one place,
caller::select_lowest_ranking, shared by the raw-record and codec paths. With no RNG left on the caller there is no per-instance state for downsampling to depend on, so theseedoption and theranddependency are removed — the determinism is structural rather than incidental.Behaviour changes
1. Which reads are retained differs from previous releases for any run that sets
--max-reads,--max-reads-per-strand, or codec--max-reads. There is no way to fix the reproducibility bug while preserving the old selection, since the old selection was the bug. Runs without a cap are byte-identical, verified against a pre-change binary.2. Both ends of a template are now retained or discarded together, which also recovers reads the old shuffle lost. On a 77k-molecule library:
cDdistribution1: 726,8541: 734,5361: 702,047 ·2: 32,489The count now matches fgbio. The depth does not, because simplex still applies the cap to the whole tag family rather than to each end — tracked separately as #723, and deliberately not fixed here.
3.
duplex --max-reads-per-strand 0is now rejected (second commit). It previously exited 0 having written an empty BAM.Tests
Every selection test was verified red under a mutant that breaks the property it claims to guard, not merely green afterward:
wrapping_abs(the abs form that exists elsewhere in this repo) fails itsort_unstable_by_keyfails it, using mixed tied/distinct keys at n=30, since all-equal keys can never separate stable from unstablehash(name + end)mutant fails exactly this test and nothing elseMurmur3(42)by running real htsjdk. Without it, a deterministic-but-biased rule passes everything above — sorting by read name, which on Illumina names sorts by tile/x/y, would systematically favour one region of the flowcell.Plus the regression test for the bug that is actually live in fgumi: capped output must be identical across
--threads 1,--threads 8, and no---threads, for simplex and duplex. Verified failing ate2ad27bf— both capped cases fail, both uncapped cases pass.Full suite: 7182 passed, 30 skipped.
Deliberately not included
No fgbio byte-parity oracle. fgbio#1166 is still open, so there is nothing stable to pin against; a fixture generated from released 4.1.0 would disagree by design, since it still uses the shuffle. Once that PR merges, this change also restores byte parity under a cap, and the oracle can follow. This is a known gap, not an oversight.
Also out of scope, both filed: #723 (simplex caps the whole MI group instead of each end) and #724 / fulcrumgenomics/fgbio#1167 (reads dropped by downsampling are counted as used in
raw_reads_usedand never reach--rejects— pre-existing and shared with fgbio).Reading order
Two commits, mirroring fgbio#1166's structure. Both build and test standalone.
fix(consensus)!: make consensus read downsampling deterministic— start atcrates/fgumi-raw-bam/src/hash.rs(the ranking primitive), thencaller::select_lowest_ranking, then the three call sites invanilla_caller.rsandcodec_caller.rs.fix(duplex): reject --max-reads-per-strand 0— reviewable independently.Closes #725
Risk: output changes for capped consensus runs, pinned by signed fgbio-compatible read-name hashes and stable pair selection;
unsafechanges: none, and CLAUDE.md allowlist changes: none; memory bounds, queue capacity, and thread/backpressure policy changes: none.seedoption.--max-reads-per-strand 0.