SamRecordClipper still clips past-mate overlap by reference distance, not query distance

[nh13]

Summary

fgumi-raw-bam::overlap was moved to query-space overlap-clip arithmetic in #759 (fgumi #752, the fgumi side of fgbio#1090). Three sibling implementations in fgumi-sam still carry the superseded reference-space formula:

• SamRecordClipper::num_bases_extending_past_mate (crates/fgumi-sam/src/clipper.rs)
• SamRecordClipper::num_bases_extending_past_mate_raw (same file)
• record_utils::num_bases_extending_past_mate (crates/fgumi-sam/src/record_utils.rs)

All three build a reference coordinate — the mate's soft-only unclipped end — by adding the mate's trailing soft clip, a query distance, to the mate's alignment end, then look up the read position there. With an indel near the read's 3' end that coordinate lands somewhere the read never sequenced: inside a deletion it has no read position at all, read_pos_at_ref_pos(..., false) returns 0, and read_length - 0 clips the entire read. An insertion moves it the other way and under-clips by the insertion's length.

Who is affected

clipper.rs's pair is reached from SamRecordClipper::clip_extending_past_mate_ends, which fgumi clip calls (src/lib/commands/clip.rs). So fgumi clip --clip-overlapping-reads-style usage still mis-clips read-through pairs with a terminal indel.

record_utils::num_bases_extending_past_mate is currently only reached from a #[cfg(test)] bridge in vanilla_caller.rs, so it has no production impact today, but it is public API and will drift further from the raw-byte path.

The consensus callers (simplex, duplex, codec) are not affected — they all route through fgumi-raw-bam::overlap, fixed in #759.

Why it was not fixed in #759

The correct algorithm needs the mate's CIGAR, not two derived reference boundaries. clipper.rs's two functions take mate_unclipped_start / mate_unclipped_end as parameters, so fixing them means changing a public signature and its call sites — a wider change than the consensus fix warranted, and one worth doing on its own so its own before/after clipping counts can be measured.

Suggested fix

Port the algorithm from bases_extending_past_mate in crates/fgumi-raw-bam/src/overlap.rs: take the last reference position both alignments cover (the first, for a negative-strand read), count the query bases each read has past it, and clip the read down to the mate's count.

Two disjoint-alignment cases come with it, and they are not the same case — see the comment below for why this distinction was worth pinning down:

• The read's 3' end faces the mate but its alignment stops short. Keep the existing extrapolation: project the read's soft-clipped tail one query base per reference base against the mate's soft-only unclipped boundary. Do not answer 0 — read-through is entirely possible here, and answering 0 leaves adapter on the read. The estimate is capped by the read's own soft clipping, so it never takes an aligned base.
• The read's 3' end faces away from the mate, i.e. the mate lies entirely on the far side. Clip nothing. That is an outward-facing template rather than read-through, and measuring it against the anchor like any other read counts all of the read's query bases as past the mate and takes its aligned ones — the destructive answer overlap clipping reproduces fgbio #1090's reference/query distance conflation #752 exists to remove.

Note that reading a deleted base as the last one before it (htsjdk's returnLastBaseIfDeleted = true) is not a fix — it replaces one reference-space answer with another and still carries the mate's soft clip in the wrong space.

overlap.rs's test suite is a ready-made oracle: the deletion-at-boundary pair, fgbio#1090's own insertion example, an ungapped control that must not move, test_disjoint_alignments_still_clip_their_read_through for the first bullet above, and test_num_bases_extending_past_mate_raw_read_entirely_past_mate for the second.

[nh13]

The "Suggested fix" section has been updated, but not because it was wrong — the part it now spells out was genuinely contested, and it is worth recording how it was settled before this gets ported.

fgbio's rework of the same defect (fulcrumgenomics/fgbio#1172) initially took the opposite reading of the disjoint case: with no shared reference position it answered 0, on the argument that without an overlap "the reads cannot have read through into each other's adapter anyway". That is the intuitive reading and it is false. Both alignments lie within the insert, so a pair that read through can still leave them disjoint — it happens whenever each alignment stops short of the other's, which soft clipping at a variant or a mis-called base is enough to produce.

Three things settled it:

1. Continuity. Two 100-base reads over a 39-base insert, 20M80S at 1000 against 80S20M at 1019, share exactly one reference position and each give up 61 bases by the query-space count. Move the mate one base right and they share none. Answering 0 steps from 61 to 0 for a one-base shift; the extrapolation gives 60, which is the correct read-through for a one-base-longer insert. The extrapolation is the analytic continuation of the query-space count, not a separate heuristic.
2. fgbio#842. That branch was added deliberately, because without it "the numBasesExtendingPastMate function always returned 0 [...] generally resulting in reads now failing when they used to pass" in FilterConsensusReads. Answering 0 reintroduces the subset of that bug where the alignments are disjoint.
3. Measurement. On a shared set of 20,000 randomly generated FR pairs spanning overlapping, touching and disjoint geometries, the 0-answering revision differed from the pre-fix formula on 28 indel-free evaluations — all of them this geometry, all of them a loss of clipping. A fix for a reference/query conflation must be inert on indel-free input, and that revision was not. The conflation provably cannot arise in this branch anyway: neither alignment places a base between the two, so neither can place an indel between them.

fgbio#1172 has been amended accordingly. The two implementations now agree on all 40,000 evaluations of that shared set.

The second bullet is the piece this issue did not have. A read whose 3' end points away from its mate is not the same case, and clipping it by the same rule takes its aligned bases — on the forward strand, 20M at 301-320 against a reverse mate 10M5S at 101-110 gave up 15 of its 20 aligned bases. That is reachable per strand by a different inconsistency (is_fr_pair_raw's forward arm reads TLEN; its reverse arm reads the mate position recorded on the record, which the mate record in hand can contradict), and #759 now returns 0 for it. A port that carries the first bullet without the second reintroduces a whole-read clip through a different door.

Worth noting separately from this issue: is_fr_pair_raw mirrors htsjdk's per-record orientation test, whose forward arm reads TLEN. htsjdk 5 changed that to prefer the mate CIGAR when the MC tag is present, which is why fgbio does not admit the outward-facing geometry at all. Bringing is_fr_pair_raw in line would close the forward-strand door properly rather than guarding behind it, but it changes FR classification everywhere and belongs in its own change.