fix(duplex): restore single-strand consensus via --min-reads 1,1,0 - #683
Conversation
|
Note Reviews pausedIt looks like this branch is under active development. To avoid overwhelming you with review comments due to an influx of new commits, CodeRabbit has automatically paused this review. You can configure this behavior by changing the Use 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 selected for processing (4)
WalkthroughThe duplex command accepts zero-valued thresholds, emits eligible single-strand consensus reads, and preprocesses overlapping mates in both execution modes. BA-only output preserves same-end R1/R2 mapping. Post-filter thresholds are enforced before emission. ChangesSingle-strand duplex consensus
Estimated code review effort: 4 (Complex) | ~45 minutes Sequence Diagram(s)sequenceDiagram
participant DuplexCommand
participant OverlappingConsensusPreprocessing
participant DuplexConsensusCaller
participant ConsensusOutput
DuplexCommand->>DuplexConsensusCaller: validate zero-valued --min-reads
DuplexCommand->>OverlappingConsensusPreprocessing: process eligible single-strand groups
OverlappingConsensusPreprocessing->>DuplexConsensusCaller: provide filtered groups
DuplexConsensusCaller->>DuplexConsensusCaller: recheck post-filter thresholds
DuplexConsensusCaller->>ConsensusOutput: emit mapped R1/R2 consensus reads
Possibly related PRs
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Comment |
|
@coderabbitai pause |
✅ Action performedReviews paused. |
01b6083 to
a748da5
Compare
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #683 +/- ##
==========================================
- Coverage 93.95% 93.95% -0.01%
==========================================
Files 178 178
Lines 108059 108301 +242
==========================================
+ Hits 101524 101750 +226
- Misses 6535 6551 +16 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
a748da5 to
a839a0d
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 2
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/lib/commands/duplex.rs (1)
1203-1225: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick winTest can no longer distinguish valid from invalid
min_reads.
validate()(lines 587-602) dropped every min-reads check; ordering is enforced only inDuplexConsensusCaller::new. Every case intest_validate_min_reads_accepts_zeroexpectsOk, so the test now passes even ifvalidate()ignoredmin_readsentirely — it cannot catch a regression that turnsvalidate()into a no-op with respect to this field. Ordering violations are covered separately byduplex_caller.rs::test_min_reads_validation_error. Consider renaming or annotating this test to state explicitly that it only guards against reintroducing the removed lower-bound rejection.Based on path instructions: "Flag assertions weaker than the stated contract."
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/lib/commands/duplex.rs` around lines 1203 - 1225, Rename or annotate test_validate_min_reads_accepts_zero to explicitly indicate it only guards that validate() does not reject zero-valued min_reads. Remove the expect_ok parameter and redundant always-true assertions, or otherwise make the test clearly assert the narrow lower-bound behavior; leave ordering validation to DuplexConsensusCaller::new and its existing tests.Source: Path instructions
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@src/lib/commands/duplex.rs`:
- Around line 604-611: Remove the duplicated padding check from
DuplexConsensusCaller’s allows_single_strand_consensus and reuse
DuplexConsensusCaller’s canonical min_yx_reads-based eligibility logic instead.
Expose a public instance method or static helper in DuplexConsensusCaller, then
update the CLI and threaded paths to call it so single-strand gating always
follows the caller’s padding semantics.
In `@tests/integration/test_duplex_command.rs`:
- Around line 438-484: Strengthen
test_duplex_single_strand_mode_emits_one_strand_molecules beyond
reader.records().count() by collecting the output records and validating both
consensus records’ identities: expected R1/R2 flags, matching molecule
identifier (MI), and non-empty sequences. Preserve the existing 2-record
expectation and the separate rejection assertion for requiring both strands.
---
Outside diff comments:
In `@src/lib/commands/duplex.rs`:
- Around line 1203-1225: Rename or annotate test_validate_min_reads_accepts_zero
to explicitly indicate it only guards that validate() does not reject
zero-valued min_reads. Remove the expect_ok parameter and redundant always-true
assertions, or otherwise make the test clearly assert the narrow lower-bound
behavior; leave ordering validation to DuplexConsensusCaller::new and its
existing tests.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: 50def66f-fc4b-46eb-a049-a58bbaf39b2a
📒 Files selected for processing (4)
crates/fgumi-consensus/src/duplex_caller.rsdocs/src/guide/duplex-consensus-calling.mdsrc/lib/commands/duplex.rstests/integration/test_duplex_command.rs
a839a0d to
7d695ac
Compare
|
Addressed the outside-diff comment on Two further fixes from a self-review pass:
|
7d695ac to
d0475cb
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
There was a problem hiding this comment.
Actionable comments posted: 1
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
crates/fgumi-consensus/src/duplex_caller.rs (1)
2169-2276: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick winSingle-strand arms skip the post-filtering minimum-reads check.
The both-strand arm (Lines 2110-2122) rebuilds
dr1/dr2and then callsduplex_consensus_has_minimum_readsbefore emitting, becausefilter_by_alignment(CIGAR filtering) can drop reads afterhas_minimum_number_of_readsalready passed on the raw pre-filter counts. The AB-only arm (Lines 2169-2216) and the BA-only arm (Lines 2217-2276) skip this recheck entirely: as soon asmin_yx_reads == 0and bothduplex_r1/duplex_r2build, the read is emitted regardless ofmin_total_readsormin_xy_reads.Before this PR,
min_yx_reads == 0was unreachable becauseDuplex::validate()rejected a zero YX value, so this gap was dormant. This PR removes that CLI-level rejection, so the gap is now reachable: with, for example,--min-reads 5,3,0, a single-strand group whose alignment-filtered read count drops below 5 (or 3) still emits a consensus, silently violating the configured threshold.Mirror the both-strand arm's check in both single-strand arms.
🐛 Proposed fix for both single-strand arms
(Some(ref r1_a), Some(ref r2_a), None, None) => { // Only AB strand - check if single-strand consensus is allowed (min_yx_reads == 0) if min_yx_reads == 0 { // Use duplex_consensus with only AB (no BA) let duplex_r1 = Self::duplex_consensus(Some(r1_a), None, None); let duplex_r2 = Self::duplex_consensus(Some(r2_a), None, None); - if let (Some(dr1), Some(dr2)) = (duplex_r1, duplex_r2) { + if let (Some(dr1), Some(dr2)) = (duplex_r1, duplex_r2) + && Self::duplex_consensus_has_minimum_reads(&dr1, min_total_reads, min_xy_reads, min_yx_reads) + && Self::duplex_consensus_has_minimum_reads(&dr2, min_total_reads, min_xy_reads, min_yx_reads) + { let empty: &[&RawRecord] = &[]; ... } } } (None, None, Some(ref r1_b), Some(ref r2_b)) => { ... if min_yx_reads == 0 { let duplex_r1 = Self::duplex_consensus(None, Some(r1_b), None); let duplex_r2 = Self::duplex_consensus(None, Some(r2_b), None); - if let (Some(dr1), Some(dr2)) = (duplex_r1, duplex_r2) { + if let (Some(dr1), Some(dr2)) = (duplex_r1, duplex_r2) + && Self::duplex_consensus_has_minimum_reads(&dr1, min_total_reads, min_xy_reads, min_yx_reads) + && Self::duplex_consensus_has_minimum_reads(&dr2, min_total_reads, min_xy_reads, min_yx_reads) + { let empty: &[&RawRecord] = &[]; ... } } }Add a test with a non-trivial threshold (e.g.
--min-reads 5,3,0) plus a dissimilar-CIGAR read that trims the surviving single-strand count below the threshold, and assert rejection.As per path instructions: "Flag silent behavior changes on the edge cases that differ between tools: ... min/max family-size thresholds."
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@crates/fgumi-consensus/src/duplex_caller.rs` around lines 2169 - 2276, Mirror the both-strand arm’s post-filter validation in the AB-only and BA-only branches of the consensus-building function: after constructing dr1 and dr2, require duplex_consensus_has_minimum_reads to pass using the filtered strand counts and configured min_total_reads/min_xy_reads/min_yx_reads before emitting or recording the pair. Add a regression test using a threshold such as 5,3,0 and a dissimilar-CIGAR read that lowers the surviving single-strand count below the threshold, asserting no consensus is emitted.Source: Path instructions
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@docs/src/guide/duplex-consensus-calling.md`:
- Around line 87-112: In the “Single-strand molecules” section, update the
sentence to use “whether it makes duplex” and change “afterwards” to
“afterward,” preserving all other wording and behavior descriptions.
---
Outside diff comments:
In `@crates/fgumi-consensus/src/duplex_caller.rs`:
- Around line 2169-2276: Mirror the both-strand arm’s post-filter validation in
the AB-only and BA-only branches of the consensus-building function: after
constructing dr1 and dr2, require duplex_consensus_has_minimum_reads to pass
using the filtered strand counts and configured
min_total_reads/min_xy_reads/min_yx_reads before emitting or recording the pair.
Add a regression test using a threshold such as 5,3,0 and a dissimilar-CIGAR
read that lowers the surviving single-strand count below the threshold,
asserting no consensus is emitted.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: c6f6f50c-7d20-4369-b88b-c0abfb8a762d
📒 Files selected for processing (4)
crates/fgumi-consensus/src/duplex_caller.rsdocs/src/guide/duplex-consensus-calling.mdsrc/lib/commands/duplex.rstests/integration/test_duplex_command.rs
d0475cb to
d149794
Compare
|
Addressed the outside-diff finding on The AB-only and BA-only arms now mirror the both-strand arm and re-check Regression test Also applied the two wording nits in |
|
@coderabbitai review |
✅ Action performedReview finished.
|
There was a problem hiding this comment.
Actionable comments posted: 1
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@src/lib/commands/duplex.rs`:
- Around line 593-612: Update validate() to enforce the --min-reads ordering
constraint (total >= XY >= YX) before execution begins, matching
DuplexConsensusCaller::new. Preserve support for unbounded nonnegative values
and valid single-strand configurations, while ensuring execute_threads_mode
rejects invalid input before creating output or starting workers.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: a7da703f-759e-49e8-9ae0-bc75d1773114
📒 Files selected for processing (4)
crates/fgumi-consensus/src/duplex_caller.rsdocs/src/guide/duplex-consensus-calling.mdsrc/lib/commands/duplex.rstests/integration/test_duplex_command.rs
d149794 to
ec29fcd
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
There was a problem hiding this comment.
Actionable comments posted: 3
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@docs/src/guide/duplex-consensus-calling.md`:
- Around line 83-86: Update the fgumi `--min-reads` example in the duplex
consensus calling documentation from space-separated values to comma-separated
`7,3,1`, while preserving the space-separated `7 3 1` form only in the
explicitly labeled fgbio command example.
In `@src/lib/commands/duplex.rs`:
- Around line 501-504: In both overlapping-consensus gates in
src/lib/commands/duplex.rs at lines 501-504 and 754-758, reorder the condition
so single_strand_allowed is evaluated before the expensive has_both_strands_raw
scan; update both sites to use short-circuit ordering while preserving the
existing consensus behavior.
In `@tests/integration/test_duplex_command.rs`:
- Around line 692-701: Update create_fully_overlapping_pair to assert that
r1_sequence and r2_sequence have equal lengths before deriving read_len or
constructing either record, matching the validation behavior of
create_duplex_read_pair_with_sequences.
🪄 Autofix
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: b1f8a47e-3f58-4485-9289-90bf24bb35bd
📒 Files selected for processing (4)
crates/fgumi-consensus/src/duplex_caller.rsdocs/src/guide/duplex-consensus-calling.mdsrc/lib/commands/duplex.rstests/integration/test_duplex_command.rs
ec29fcd to
f9d610f
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
Note GitHub couldn't provide a complete incremental comparison for this pull request, so CodeRabbit is performing a full review instead. This review may take a little longer. |
There was a problem hiding this comment.
Actionable comments posted: 1
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
crates/fgumi-consensus/src/duplex_caller.rs (1)
6576-6826: 🗄️ Data Integrity & Integration | 🟠 Major | 🏗️ Heavy liftPin single-strand output to an independent fgbio oracle.
The current tests derive expectations from fgumi inputs or compare two fgumi modes. They can pass when BA read-end mapping, post-filter threshold handling, or overlap masking differs from fgbio.
crates/fgumi-consensus/src/duplex_caller.rs#L6576-L6826: Add programmatically generated AB-only and BA-only1,1,0fixtures with complete captured fgbio 4.1.0 R1/R2 output identity, including flags, bases, MI/RX, depth/error tags, and applicable methylation tags.tests/integration/test_duplex_command.rs#L725-L786: Assert complete command output identity against the captured overlap-correction oracle in both execution modes.Keep fixture regeneration ignored. Do not require fgbio or a JVM in default CI.
Based on learnings, fgbio-backed tests must use fixed expected outputs and gated regeneration. As per path instructions, “New correctness-critical behavior needs an INDEPENDENT oracle … not just a self-consistency check.”
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@crates/fgumi-consensus/src/duplex_caller.rs` around lines 6576 - 6826, Add programmatically generated AB-only and BA-only 1,1,0 fixtures in DuplexConsensusCaller tests, using captured fgbio 4.1.0 outputs as fixed independent expectations for complete R1/R2 identity: flags, bases, MI/RX, depth/error tags, and applicable methylation tags; keep regeneration gated and ignored so default CI needs neither fgbio nor a JVM. Also update tests/integration/test_duplex_command.rs lines 725-786 to assert complete command-output identity against the captured overlap-correction oracle in both execution modes, while preserving the existing fixture-based execution path.Sources: Path instructions, Learnings
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@tests/integration/test_duplex_command.rs`:
- Around line 436-534: Extend
test_duplex_single_strand_mode_emits_one_strand_molecules to run the same
one-template input with --min-reads 2,2,0 for both thread configurations and
assert that no consensus records are emitted. Ensure the assertion verifies the
paired-end template counts as one template rather than two individual read
records, while preserving the existing successful 1,1,0 and rejected 1 cases.
---
Outside diff comments:
In `@crates/fgumi-consensus/src/duplex_caller.rs`:
- Around line 6576-6826: Add programmatically generated AB-only and BA-only
1,1,0 fixtures in DuplexConsensusCaller tests, using captured fgbio 4.1.0
outputs as fixed independent expectations for complete R1/R2 identity: flags,
bases, MI/RX, depth/error tags, and applicable methylation tags; keep
regeneration gated and ignored so default CI needs neither fgbio nor a JVM. Also
update tests/integration/test_duplex_command.rs lines 725-786 to assert complete
command-output identity against the captured overlap-correction oracle in both
execution modes, while preserving the existing fixture-based execution path.
🪄 Autofix
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: bbe0689d-15c9-4ef7-9207-5b4524717014
📒 Files selected for processing (4)
crates/fgumi-consensus/src/duplex_caller.rsdocs/src/guide/duplex-consensus-calling.mdsrc/lib/commands/duplex.rstests/integration/test_duplex_command.rs
f9d610f to
de7b817
Compare
|
Addressed both findings from the latest review. Template-count boundary (inline thread, resolved). The literal suggestion — Both accepting cases now assert full pair identity (two records, exactly one R1, molecule id, bases) through a shared fgbio oracle for single-strand output (outside-diff, no thread to resolve). Added Pinned per record: flags, bases, qualities, The Overlap-correction oracle.
fgumi's output is byte-identical to fgbio's on both settings. The asymmetry is fgbio's: the mates disagree at the last reference base, which masks to N at qual 2 in both, and in sequencing orientation that is the end of R1 (trimmed) and the start of R2 (kept as Three fixes from a self-review pass over the diff:
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
Note GitHub couldn't provide a complete incremental comparison for this pull request, so CodeRabbit is performing a full review instead. This review may take a little longer. |
There was a problem hiding this comment.
Actionable comments posted: 1
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@crates/fgumi-consensus/src/duplex_caller.rs`:
- Around line 7368-7378: Update the tag validation loop in the duplex record
assertion to use get_i16_array_tag for SamTag::BD_BASES and SamTag::BE_BASES,
while retaining get_string_tag checks for SamTag::BC_BASES and SamTag::BQ.
Ensure the absence assertion detects the emitted B:s arrays for bd and be
without changing the existing error behavior.
🪄 Autofix
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: f364ff97-8a94-472d-a0b9-bc0eac8abff4
📒 Files selected for processing (4)
crates/fgumi-consensus/src/duplex_caller.rsdocs/src/guide/duplex-consensus-calling.mdsrc/lib/commands/duplex.rstests/integration/test_duplex_command.rs
0.4.0 accepted a third `--min-reads` value of 0 -- fgbio's single-strand mode, where a molecule observed on only one strand still yields a consensus -- and `DuplexConsensusCaller` has always implemented it. #601 then added a lower bound to `Duplex::validate()` for consistency with simplex and codec, but its blanket `contains(&0)` also rejected the per-strand slots, taking the mode away in 0.5.0. Drop the bound entirely rather than narrowing it, matching fgbio, which validates only that the values run from least to most stringent. For simplex and codec a 0 really does make the minimum-family-size filter a silent no-op, which is what #601 was about; duplex is different in both slots. A 0 per strand is the single-strand mode, and a total of 0 is equivalent to 1 rather than degenerate, because the check it feeds (`min_total <= num_xy + num_yx`) is only reached for a group that already has at least one read. Restoring the flag alone was not enough to reproduce fgbio, because two paths that are unreachable while the mode is gated were wrong: The overlapping-bases consensus was skipped for groups without both strands, on the grounds that no duplex consensus can be built from them. That stops being true when the third value is 0, so single-strand molecules with overlapping mates kept uncorrected bases. Gate on whether the group can actually produce a consensus instead. The BA-only arm mapped output R1 from BA-R2 and R2 from BA-R1 (#298), to keep output R1 on the physical strand AB-R1 would have sequenced. That formula holds only while both strand groups are present: fgbio takes a lone group as its "AB" side, so output R1 comes from that group's R1s whether it is /A or /B. The previous mapping put the opposite end of the molecule in R1 for /B-only molecules than for every other molecule class. Build each consensus read from the input reads of the same end. #298's other half is preserved -- BA still goes in the BA slot, so is_ba_only stays true and per-strand methylation tags remain bm/bu/bt. On 20,000 simulated duplex molecules, output is now identical to fgbio 4.1.0 for `--min-reads 1,1,0` (40,000 records) and `--min-reads 0` (40,000 records), single- and multi-threaded; `--min-reads 1` is unchanged and still identical to fgbio. Closes #678
de7b817 to
fd78a27
Compare
|
Good catch — the finding is correct and the assertion was genuinely vacuous.
Swept the class rather than fixing only the cited site. Every tag-absence assertion in the changed files, checked against how the tag is actually written:
No other site had the mismatch. The presence assertions in the same test ( Added a comment at the site recording why the accessor has to match the written type, so the next edit does not reintroduce it. |
|
@coderabbitai review |
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
Note GitHub couldn't provide a complete incremental comparison for this pull request, so CodeRabbit is performing a full review instead. This review may take a little longer. |
This branch was previously deployed
1 participant
Closes #678.
fgumi duplex --min-reads 1,1,0— fgbio's single-strand mode, where a molecule observed on only one strand still yields a consensus — worked in 0.4.0 and errors out in 0.5.0.DuplexConsensusCallerhas always implemented the mode; what changed is the CLI. #601 added a lower bound toDuplex::validate()for consistency withsimplexandcodec, but its blanketcontains(&0)also rejected the per-strand slots, where a 0 is meaningful.The bound is dropped entirely rather than narrowed, matching fgbio, which validates only that the values run from least to most stringent.
simplexandcodeckeep their check, because there a 0 really does make the minimum-family-size filter a silent no-op — that part of #601 stands.duplexdiffers in both slots: a 0 per strand is the single-strand mode, and a total of 0 is equivalent to 1 rather than degenerate, because the check it feeds (min_total <= num_xy + num_yx) is only reached for a group that already has at least one read.Restoring the flag alone was not enough to reproduce fgbio, because two code paths that are unreachable while the mode is gated were wrong.
The overlapping-bases consensus was skipped for single-strand groups. It is gated on the group having both strands, "no duplex possible anyway" — which stops being true once the third value is 0, so single-strand molecules with overlapping mates kept uncorrected bases. It is now gated on whether the group can actually produce a consensus.
BA-only molecules came out with R1/R2 swapped relative to fgbio. #298 mapped output R1 from BA-R2 and R2 from BA-R1, to keep output R1 on the physical strand that AB-R1 would have sequenced. That formula holds only while both strand groups are present; fgbio does not extend it to a lone group, taking whichever group is present as its "AB" side, so output R1 comes from that group's R1s whether it is
/Aor/B. Verified in the data — within a molecule, AB-R1 and BA-R1 are on opposite strands at opposite ends:so the previous mapping put the opposite end of the molecule in R1 for
/B-only molecules than for every other molecule class. Each consensus read is now built from the input reads of the same end. #298's other half is preserved: BA still goes in the BA slot, sois_ba_onlystays true and per-strand methylation tags are still emitted asbm/bu/bt— that flag records which strand group the bases came from, which the end does not change.test_duplex_ba_only_methylation_tags_use_bottom_strandis untouched and passes;test_duplex_ba_only_mate_pair_mappingis updated to the new expectation.Verification
On 20,000 simulated duplex molecules (81,368 raw reads), against fgbio 4.1.0, comparing with
fgumi compare bams --command duplex:--min-reads11,1,01,1,0 --threads 40Before this change, 14,014 of those 40,000 records differed from fgbio, in two ways: R1/R2 swaps on
/B-only molecules, and missing overlap corrections on single-strand molecules whose mates overlap. Disabling the overlapping consensus on both sides isolates the first defect, and it accounts for exactly 9,048 differing records — the R1/R2 pairs of all 4,524/B-only molecules, with the 4,760/A-only and 10,716 two-strand molecules matching fgbio exactly.--min-reads 1is unchanged.New tests, each watched failing first:
--min-readslower-bound cases; single-strand molecules are emitted with1,1,0and still rejected with1; consensus reads follow input read ends for/A-only and/B-only molecules; the overlapping consensus is applied to single-strand groups. All are parameterized over single-threaded and--threads 2.cargo ci-test(6,903 tests),cargo ci-fmt, andcargo ci-lintpass.Note on fgbio's docs
fgbio implements this mode but does not document it — neither
CallDuplexConsensusReadsnorFilterConsensusReadsmentions a 0 value, andDuplexConsensusCaller.duplexConsensus's own scaladoc ("If either of the incoming reads are undefined, the duplex read will be undefined") contradicts itscase (Some(a), None)arm. This PR documents the mode on the fgumi side, in the duplex consensus guide and in--min-reads' help text.Summary by CodeRabbit
New Features
--min-reads 1,1,0.bD:i:0.Bug Fixes
--min-readsconfigurations are rejected before output is created.Documentation