feat(zipper): add --restore-unconverted-bases flag for EM-seq

[nh13]

Summary

• Adds --restore-unconverted-bases flag to fgumi zipper for EM-seq post-alignment base restoration
• After bwameth alignment (C→T converted reads), restores bases to unconverted state at reference positions: T→C for top-strand (YD=f), A→G for bottom-strand (YD=r)
• Ensures consensus calling sees original bases so methylation is determined by the consensus caller
• Defines YD_FORWARD/YD_REVERSE constants (was inline byte literals)
• Uses reference_length() utility for ref span computation (was manual CIGAR walk)

Stack

3/4 — EM-seq support (depends on #169)

Test plan

• Tests for top-strand and bottom-strand base restoration
• Tests for indels, unmapped reads, missing YD tag, already-unconverted bases
• Existing zipper tests pass unchanged
• cargo ci-test && cargo ci-fmt && cargo ci-lint passes

[codecov]

Codecov Report

❌ Patch coverage is 90.83969% with 24 lines in your changes missing coverage. Please review.
✅ Project coverage is 89.23%. Comparing base (6226b8e) to head (6993095).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/lib/commands/zipper.rs	90.83%	24 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff            @@
##             main     #170    +/-   ##
========================================
  Coverage   89.23%   89.23%            
========================================
  Files         120      120            
  Lines       59278    59536   +258     
========================================
+ Hits        52894    53126   +232     
- Misses       6384     6410    +26     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[nh13]

@coderabbitai review

[nh13]

@coderabbitai review

[coderabbitai]

✅ Actions performed

Review triggered.

Note: CodeRabbit is an incremental review system and does not re-review already reviewed commits. This command is applicable only when automatic reviews are paused.

[coderabbitai]

Actionable comments posted: 3

Caution

Some comments are outside the diff and can’t be posted inline due to platform limitations.

⚠️ Outside diff range comments (2)

crates/fgumi-consensus/src/duplex_caller.rs (1)

973-1000: ⚠️ Potential issue | 🟠 Major

Don’t zero out real errors at conversion sites.

Line 974 turns the whole position into 0u16. That removes the expected C/T or G/A disagreement, but it also erases ordinary within-family errors already present in a.errors / b.errors or the underlying source reads. errors and cE will read artificially low anywhere a methylation site also has normal sequencing/PCR noise. Keep the conversion pair neutral, but still accumulate non-conversion errors here.

🤖 Prompt for AI Agents

Verify each finding against the current code and only fix it if needed.

In `@crates/fgumi-consensus/src/duplex_caller.rs` around lines 973 - 1000, The
current branch that returns 0u16 when is_conversion_artifact wipes out real
sequencing/PCR errors; instead, in the is_conversion_artifact case compute
error_count by accumulating non-conversion disagreements (just like the
source_reads branch or the approximate branch) but ignore mismatches that are
the expected conversion pair (e.g., C↔T or G↔A) — update the
is_conversion_artifact branch to: if source_reads is Some, iterate source_reads
and increment for Self::is_error(sr.bases[i], raw_base) only when the mismatch
is not a conversion-pair; otherwise use the same approximate logic with
a_err/b_err/a_dep/b_dep but subtract only the expected conversion-induced
disagreements (do not zero out a_err or b_err), then clamp and cast as before;
refer to variables is_conversion_artifact, source_reads, Self::is_error,
raw_base, a.errors, b.errors, a.depths, b.depths and index i when making the
change.

crates/fgumi-consensus/src/methylation.rs (1)

240-252: ⚠️ Potential issue | 🟠 Major

Combined MM/ML still drops opposite-strand sites.

This builder can emit only one base class per call. That works for per-strand annotations, but a combined duplex annotation can contain both top-strand C sites and bottom-strand G sites in the same read. In that case, one class is silently omitted from MM/ML while it still survives in cu/ct. Emit both segments in a single combined tag, or skip the combined MM/ML tag until multi-segment support exists.

Also applies to: 268-277

🤖 Prompt for AI Agents

Verify each finding against the current code and only fix it if needed.

In `@crates/fgumi-consensus/src/methylation.rs` around lines 240 - 252, The
current loop over annotation.evidence (using consensus_bases, ev.is_ref_c,
track_base, skips and skip_count) builds a single-base-class MM/ML segment and
thus drops opposite-strand sites when a combined duplex annotation contains both
C and G; change the logic to first scan the evidence for mixed-strand presence
(i.e., both ev.is_ref_c==true and ev.is_ref_c==false for bases matching
track_base/consensus_bases) and if both strands are present, do not emit the
single-class MM/ML segment here (skip creating/pushing MM/ML for this read) —
alternatively, implement multi-segment emission by aggregating both C and G
segments into one combined tag, but if that is not implemented yet, explicitly
skip emitting the combined MM/ML tag for reads with mixed-strand evidence (apply
the same guard in both places around the blocks using ev.is_ref_c at the two
ranges noted).

🧹 Nitpick comments (1)

src/lib/commands/zipper.rs (1)

3803-4080: Add one end-to-end test for the flag path.

These tests cover the walker itself, but they never hit execute/process_raw with restore_unconverted_bases: true, where the code loads ReferenceReader and switches from merge_raw to merge. One CLI-level case would lock down that wiring.

🤖 Prompt for AI Agents

Verify each finding against the current code and only fix it if needed.

In `@src/lib/commands/zipper.rs` around lines 3803 - 4080, Add an end-to-end test
that invokes the CLI path which sets restore_unconverted_bases: true so
execute/process_raw is exercised; specifically, create a test that calls the
function (or integration entrypoint) that runs execute/process_raw with
restore_unconverted_bases=true, ensuring the code path loads ReferenceReader and
picks merge (not merge_raw). The test should feed a small FASTA (use
create_test_fasta), a SAM/BAM-like input record covering converted bases, and
assert output sequence changes after process_raw completes; reference symbols to
locate code: execute/process_raw, restore_unconverted_bases flag, merge_raw,
merge, ReferenceReader.

🤖 Prompt for all review comments with AI agents

Verify each finding against the current code and only fix it if needed.

Inline comments:
In `@crates/fgumi-consensus/src/vanilla_caller.rs`:
- Around line 645-647: The code currently mutates SourceRead.bases with C/G
restorations and then calls create_consensus_from_source_reads to build the
emitted sequence, which causes restored bases to appear in the final consensus;
instead, preserve SourceRead.bases as the original converted-space bases and
only apply normalization/restoration on a separate temporary vector used for
scoring/error models (e.g., in the scorer path that currently overwrites bases
around the block noted). Update the scorer logic to work on a copy (e.g., make a
normalized_bases or scoring_reads copy) so create_consensus_from_source_reads
and any downstream consensus emission always consume the original
SourceRead.bases, while the error accounting uses the normalized copy; ensure
functions/methods referencing SourceRead.bases (and any callers of
create_consensus_from_source_reads) are adapted to accept/use the separate
scoring copy so no restored C/G ever gets emitted (consistent with methylation
handling).
- Around line 681-720: The anchor read (anchor / anchor_idx) and derived arrays
(ref_positions, ref_bases) are currently passed whole to
annotate_simplex_methylation which builds evidence sized from the anchor's full
length, but create_consensus_from_source_reads only emits consensus of length
consensus_len (e.g., when shorter reads determine consensus); to fix, clip the
anchor-based inputs to consensus_len before calling
annotate_simplex_methylation: compute consensus_len, then trim the anchor slice
used (source_reads[anchor_idx].bases), trim/refilter ref_positions and
corresponding ref_bases to the first consensus_len positions, and pass those
clipped slices to annotate_simplex_methylation so annotation length matches the
produced consensus.

In `@src/lib/commands/zipper.rs`:
- Around line 1125-1127: When you change the read sequence via
record.sequence_mut() (i.e., when modified and seq is applied to SEQ), also
remove or recompute mapper-derived mismatch/auxiliary tags (e.g., "NM", "MD",
"AS", and any other alignment-dependent tags) so the record stays consistent;
locate the block that sets *record.sequence_mut() = seq.into() and after that
call either remove those tags via
record.remove_aux(b"NM")/record.remove_aux(b"MD") (or invoke your existing
tag-recompute utility) or call the function that recalculates NM/MD from
CIGAR/sequence, ensuring the BAM/SAM record no longer contains stale NM/MD
values when SEQ changes.

---

Outside diff comments:
In `@crates/fgumi-consensus/src/duplex_caller.rs`:
- Around line 973-1000: The current branch that returns 0u16 when
is_conversion_artifact wipes out real sequencing/PCR errors; instead, in the
is_conversion_artifact case compute error_count by accumulating non-conversion
disagreements (just like the source_reads branch or the approximate branch) but
ignore mismatches that are the expected conversion pair (e.g., C↔T or G↔A) —
update the is_conversion_artifact branch to: if source_reads is Some, iterate
source_reads and increment for Self::is_error(sr.bases[i], raw_base) only when
the mismatch is not a conversion-pair; otherwise use the same approximate logic
with a_err/b_err/a_dep/b_dep but subtract only the expected conversion-induced
disagreements (do not zero out a_err or b_err), then clamp and cast as before;
refer to variables is_conversion_artifact, source_reads, Self::is_error,
raw_base, a.errors, b.errors, a.depths, b.depths and index i when making the
change.

In `@crates/fgumi-consensus/src/methylation.rs`:
- Around line 240-252: The current loop over annotation.evidence (using
consensus_bases, ev.is_ref_c, track_base, skips and skip_count) builds a
single-base-class MM/ML segment and thus drops opposite-strand sites when a
combined duplex annotation contains both C and G; change the logic to first scan
the evidence for mixed-strand presence (i.e., both ev.is_ref_c==true and
ev.is_ref_c==false for bases matching track_base/consensus_bases) and if both
strands are present, do not emit the single-class MM/ML segment here (skip
creating/pushing MM/ML for this read) — alternatively, implement multi-segment
emission by aggregating both C and G segments into one combined tag, but if that
is not implemented yet, explicitly skip emitting the combined MM/ML tag for
reads with mixed-strand evidence (apply the same guard in both places around the
blocks using ev.is_ref_c at the two ranges noted).

---

Nitpick comments:
In `@src/lib/commands/zipper.rs`:
- Around line 3803-4080: Add an end-to-end test that invokes the CLI path which
sets restore_unconverted_bases: true so execute/process_raw is exercised;
specifically, create a test that calls the function (or integration entrypoint)
that runs execute/process_raw with restore_unconverted_bases=true, ensuring the
code path loads ReferenceReader and picks merge (not merge_raw). The test should
feed a small FASTA (use create_test_fasta), a SAM/BAM-like input record covering
converted bases, and assert output sequence changes after process_raw completes;
reference symbols to locate code: execute/process_raw, restore_unconverted_bases
flag, merge_raw, merge, ReferenceReader.

🪄 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: Organization UI

Review profile: CHILL

Plan: Pro

Run ID: 7eca68ef-c545-406b-85c3-314c61f081e5

📥 Commits

Reviewing files that changed from the base of the PR and between 32755a9 and dc4ca4b.

📒 Files selected for processing (6)

• crates/fgumi-consensus/src/duplex_caller.rs
• crates/fgumi-consensus/src/filter.rs
• crates/fgumi-consensus/src/methylation.rs
• crates/fgumi-consensus/src/vanilla_caller.rs
• crates/fgumi-raw-bam/src/cigar.rs
• src/lib/commands/zipper.rs

✅ Files skipped from review due to trivial changes (1)

• crates/fgumi-consensus/src/filter.rs

🚧 Files skipped from review as they are similar to previous changes (1)

• crates/fgumi-raw-bam/src/cigar.rs

[coderabbitai]

♻️ Duplicate comments (1)

crates/fgumi-consensus/src/vanilla_caller.rs (1)

635-647: ⚠️ Potential issue | 🔴 Critical

Keep EM-seq restoration out of emitted consensus bases.

annotate_and_normalize() rewrites SourceRead.bases, and both call sites then feed those mutated reads into create_consensus_from_source_reads(). At converted EM-seq sites this makes the emitted consensus switch from aligned-space T/A to restored C/G, which changes the sequence itself instead of just the error model. Keep a separate normalized/scoring copy, and continue emitting bases plus stored source_reads from the original converted-space reads.

Possible direction

-        let (methylation, source_reads) = if self.options.em_seq {
-            let (annot, normalized) = self.annotate_and_normalize(source_reads);
-            (annot, normalized)
+        let original_source_reads = source_reads;
+        let (methylation, scoring_reads) = if self.options.em_seq {
+            let (annot, normalized) = self.annotate_and_normalize(original_source_reads.clone());
+            (annot, normalized)
         } else {
-            (None, source_reads)
+            (None, original_source_reads.clone())
         };
 
-        let (bases, quals, depths, errors) =
-            self.create_consensus_from_source_reads(&source_reads)?;
+        let (bases, _, _, _) =
+            self.create_consensus_from_source_reads(&original_source_reads)?;
+        let (_, quals, depths, errors) =
+            self.create_consensus_from_source_reads(&scoring_reads)?;

Apply the same split in process_subgroup, and keep source_reads: Some(original_source_reads).

Also applies to: 726-738, 1213-1224

🤖 Prompt for AI Agents

Verify each finding against the current code and only fix it if needed.

In `@crates/fgumi-consensus/src/vanilla_caller.rs` around lines 635 - 647,
annotate_and_normalize mutates SourceRead.bases causing EM-seq converted bases
to be restored in the emitted consensus; preserve original reads and use a
separate normalized copy for scoring. Before calling annotate_and_normalize (in
the shown block and the other two locations called out and in process_subgroup),
clone or otherwise save original_source_reads, call let (annot, normalized) =
self.annotate_and_normalize(original_source_reads.clone()) (or make
annotate_and_normalize return a new vec) and use normalized for
scoring/quality/depth/error calculations but pass original_source_reads into the
code path that builds/emits consensus bases (or adjust
create_consensus_from_source_reads to accept separate scoring_reads vs
emission_reads). Ensure you reference annotate_and_normalize,
create_consensus_from_source_reads, process_subgroup, and the source_reads
variable when making the change.

🧹 Nitpick comments (2)

src/lib/commands/zipper.rs (1)

191-192: CLI flag pattern inconsistency.

Other boolean flags in this struct use num_args = 0..=1, default_missing_value = "true", action = clap::ArgAction::Set, and value_parser = parse_bool. This flag lacks that configuration, so --restore-unconverted-bases alone will error—users must write --restore-unconverted-bases true.

Suggested fix

-    #[arg(long = "restore-unconverted-bases", default_value = "false")]
+    #[arg(long = "restore-unconverted-bases", default_value = "false", num_args = 0..=1, default_missing_value = "true", action = clap::ArgAction::Set, value_parser = parse_bool)]
     pub restore_unconverted_bases: bool,

🤖 Prompt for AI Agents

Verify each finding against the current code and only fix it if needed.

In `@src/lib/commands/zipper.rs` around lines 191 - 192, The CLI boolean flag
field restore_unconverted_bases uses only default_value which forces users to
pass an explicit value; change its attribute to match other booleans by adding
num_args = 0..=1, default_missing_value = "true", action = clap::ArgAction::Set,
and value_parser = parse_bool (keep the existing default_value = "false"), e.g.
replace the current #[arg(long = "restore-unconverted-bases", default_value =
"false")] with the same full arg configuration used by the other boolean flags
so --restore-unconverted-bases can be passed without an explicit value.

crates/fgumi-consensus/src/duplex_caller.rs (1)

5811-5845: Assert the zero-error behavior too.

This test now proves the bottom-strand artifact is recognized, but it still doesn’t lock in the actual regression fix from Lines 973-1000. Add an assertion on duplex.errors[0] == 0 so cE inflation can’t slip back in unnoticed.

Suggested assertion

         // Should call G (unconverted), not A
         assert_eq!(duplex.bases[0], b'G');
         // Quality summed (conversion artifact path)
         assert_eq!(duplex.quals[0], 55);
+        // Conversion artifacts should not be counted as duplex errors
+        assert_eq!(duplex.errors[0], 0);

🤖 Prompt for AI Agents

Verify each finding against the current code and only fix it if needed.

In `@crates/fgumi-consensus/src/duplex_caller.rs` around lines 5811 - 5845, Add an
assertion to lock in zero-error behavior in
test_duplex_em_seq_bottom_strand_conversion_artifact: after calling
DuplexConsensusCaller::duplex_consensus and extracting duplex, assert that
duplex.errors[0] == 0 to ensure cE inflation doesn't regress; update the test
function test_duplex_em_seq_bottom_strand_conversion_artifact to include this
check alongside the existing base and quality assertions.

🤖 Prompt for all review comments with AI agents

Verify each finding against the current code and only fix it if needed.

Duplicate comments:
In `@crates/fgumi-consensus/src/vanilla_caller.rs`:
- Around line 635-647: annotate_and_normalize mutates SourceRead.bases causing
EM-seq converted bases to be restored in the emitted consensus; preserve
original reads and use a separate normalized copy for scoring. Before calling
annotate_and_normalize (in the shown block and the other two locations called
out and in process_subgroup), clone or otherwise save original_source_reads,
call let (annot, normalized) =
self.annotate_and_normalize(original_source_reads.clone()) (or make
annotate_and_normalize return a new vec) and use normalized for
scoring/quality/depth/error calculations but pass original_source_reads into the
code path that builds/emits consensus bases (or adjust
create_consensus_from_source_reads to accept separate scoring_reads vs
emission_reads). Ensure you reference annotate_and_normalize,
create_consensus_from_source_reads, process_subgroup, and the source_reads
variable when making the change.

---

Nitpick comments:
In `@crates/fgumi-consensus/src/duplex_caller.rs`:
- Around line 5811-5845: Add an assertion to lock in zero-error behavior in
test_duplex_em_seq_bottom_strand_conversion_artifact: after calling
DuplexConsensusCaller::duplex_consensus and extracting duplex, assert that
duplex.errors[0] == 0 to ensure cE inflation doesn't regress; update the test
function test_duplex_em_seq_bottom_strand_conversion_artifact to include this
check alongside the existing base and quality assertions.

In `@src/lib/commands/zipper.rs`:
- Around line 191-192: The CLI boolean flag field restore_unconverted_bases uses
only default_value which forces users to pass an explicit value; change its
attribute to match other booleans by adding num_args = 0..=1,
default_missing_value = "true", action = clap::ArgAction::Set, and value_parser
= parse_bool (keep the existing default_value = "false"), e.g. replace the
current #[arg(long = "restore-unconverted-bases", default_value = "false")] with
the same full arg configuration used by the other boolean flags so
--restore-unconverted-bases can be passed without an explicit value.

---

ℹ️ Review info

⚙️ Run configuration

Configuration used: Organization UI

Review profile: CHILL

Plan: Pro

Run ID: 6e7c37ac-b984-4a40-b749-433eaf4bd2e4

📥 Commits

Reviewing files that changed from the base of the PR and between dc4ca4b and 28552ac.

📒 Files selected for processing (6)

• crates/fgumi-consensus/src/duplex_caller.rs
• crates/fgumi-consensus/src/filter.rs
• crates/fgumi-consensus/src/methylation.rs
• crates/fgumi-consensus/src/vanilla_caller.rs
• crates/fgumi-raw-bam/src/cigar.rs
• src/lib/commands/zipper.rs

🚧 Files skipped from review as they are similar to previous changes (1)

• crates/fgumi-raw-bam/src/cigar.rs