Repository navigation
fix(clip): unmap fully-clipped reads and match fgbio unmap/mate fidelity (CLIP-01) - #492
Conversation
WalkthroughFully clipped reads now become unmapped in typed and raw paths while preserving and normalizing sequence data. Paired clipping now synchronizes mate coordinates, flags, tags, and template lengths through ChangesClipper unmap behavior
Mate-pair info synchronization
Estimated code review effort: 4 (Complex) | ~60 minutes Possibly related PRs
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Comment |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #492 +/- ##
==========================================
+ Coverage 91.12% 91.17% +0.05%
==========================================
Files 78 78
Lines 51540 51916 +376
==========================================
+ Hits 46964 47334 +370
- Misses 4576 4582 +6 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…t-order (R2-CLIP-02/03/04) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
c7535af to
389e375
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@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 `@crates/fgumi-sam/src/clipper.rs`:
- Around line 5969-6057: The new clipper tests are parameterized with runtime
`for` loops, which should be converted to `rstest` cases for per-case reporting.
Refactor the `crosscheck_unmap_*` tests in `clipper.rs` to use `#[rstest]` with
`#[case(...)]` for `ClippingMode` (and any clip/bases variations) instead of
looping inside the test body, while keeping the existing assertions and coverage
intact.
🪄 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: d27c0889-7069-4598-8416-07f25da78ae2
📒 Files selected for processing (2)
crates/fgumi-sam/src/clipper.rssrc/lib/commands/clip.rs
…t-order (R2-CLIP-02/03/04) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
389e375 to
fb85484
Compare
|
@coderabbitai review |
✅ Action performedReview finished.
|
…t-order (R2-CLIP-02/03/04) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
|
@coderabbitai review |
✅ Action performedReview finished.
|
fgbio SamRecordClipper.clip{Start,End}OfAlignment unmaps a read (via
htsjdk SAMUtils.makeReadUnmapped) when the requested clip is >= the
number of clippable (alignment + insertion) bases, keeping the bases and
qualities intact and emptying the CIGAR. fgumi instead clipped every
aligned base, leaving the read mapped at a stale POS with an all-clip
CIGAR (e.g. 50H) and an empty sequence — silent corruption (CLIP-01).
Add the numberOfClippableBases guard to both SamRecordClipper and
RawRecordClipper (clip_start_of_alignment / clip_end_of_alignment): when
clippable <= requested, unmap the read (clear reference/pos/mapq/CIGAR/
TLEN and the duplicate/secondary/supplementary/proper-pair flags, set
UNMAPPED, reverse-complement bases + reverse quals for reverse-strand
reads) and return the clippable count. Two fgumi tests that pinned the
old all-soft-clip behavior are updated to assert the unmap.
make_read_unmapped[_raw] also reorients the strand-sensitive aux tags of
reverse-strand reads, matching htsjdk reverseComplement(true): the
base-oriented tags E2/SQ are reverse-complemented and the
quality-oriented tags OQ/U2 are reversed (SAMRecord.TAGS_TO_REVERSE_-
COMPLEMENT / TAGS_TO_REVERSE). Without this, a reverse-strand read
carrying e.g. OQ that clipping unmaps kept a forward-oriented OQ,
diverging from fgbio. Only the Z-string encoding is handled (the sole
form these tags take in practice); the typed and raw paths use the same
tag sets so they stay in lockstep. Adds E2/U2/SQ to the SamTag catalog.
Because fgbio ClipBam then calls SamPairUtil.setMateInfo on every pair,
replace clip.rs's partial MC/MQ-only mate update with set_mate_info_raw,
a faithful port of htsjdk setMateInfo (all three branches: both mapped,
both unmapped, one of each) covering mate ref/pos/strand, the
mate-unmapped flag, MQ/MC tags, and TLEN. This makes paired output match
fgbio when a mate unmaps, and also fixes stale mate-pos/TLEN in the
ordinary both-mapped case.
In the one-mapped setMateInfo branch, the mapped read's mate-reverse
flag now reflects the unmapped read's actual strand (SamPairUtil.java:267)
rather than an assumed false. The clipper's own unmap clears REVERSE, but
a read that arrived unmapped-on-input may still carry it, and
set_mate_info_raw runs on every pair, so such a read reaches this branch.
Parity verified end-to-end against fgbio 4.1.0 ClipBam
(reports/fgbio-parity-fixtures/clip-01-unmap-fully-clipped.sh):
before -> DIFFER, after -> IDENTICAL (fgumi compare bams --mode content)
for a fully-clipped fragment, a pair where both mates unmap, and normal
overlapping / fixed / past-mate paired clipping. New unit tests cross-check
the tag reorientation (typed vs raw, negative- and positive-strand) and
the one-mapped mate-reverse flag.
Refs: CLIP-01.
fb85484 to
6794f10
Compare
…t-order (R2-CLIP-02/03/04) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
…t-order (R2-CLIP-02/03/04) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
…t-order (R2-CLIP-02/03/04) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
…t-order (R2-CLIP-02/03/04) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
…t-order (R2-CLIP-02/03/04) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
…t-order (R2-CLIP-02/03/04) (#502) Three round-2 clip parity fixes, stacked on #492 (CLIP-01 / setMateInfo): - R2-CLIP-02: the clip loop was `if len==1 {frag} else if len==2 {pair}` with no else, so any template with a secondary/supplementary read (len > 2) was left COMPLETELY unclipped with no mate fixing. Now locate the primary R1/R2 by SAM flags (fgbio `Template.r1`/`r2`), clip them, and repair mate info on every supplementary alignment via a port of htsjdk 5.0.0 `SamPairUtil.setMateInformationOnSupplementalAlignment` (mate ref/pos/strand, mate-unmapped flag, TLEN = -matePrimary.TLEN, MC when mate mapped, and MQ unconditionally). Both the single-threaded and threaded paths are refactored. - R2-CLIP-03: `--clip-bases-past-mate` bounded the mate window with unclipped start/end (soft + HARD clips); fgbio uses `mate.unSoftClipped{Start,End}` (soft only, since hard-clipped bases are physically absent). Added soft-only `unsoftclipped_{start,end}[_raw]` helpers and used them in both clippers, so a hard-clipped mate no longer widens the window and under-clips. - R2-CLIP-04: `--sort-order` on clip only relabeled the SO header field (a no-op); removed the flag entirely (rely on a downstream `fgumi sort`). fgbio-oracle parity (fgbio 4.1.0 ClipBam, default Hard): - CHIM (len-3 chimeric template, --clip-overlapping-reads): before DIFFER (fgumi left 100M/100M unclipped, supp mate-info stale) -> after MATCH (75M25H/25H75M, supp gets MC/MQ/mate-pos/TLEN like fgbio). - HARDMATE (FR dovetail, reverse mate 50M50H, --clip-bases-past-mate): before DIFFER (fgumi 110M10H) -> after MATCH (60M60H, past mate's soft end). - SORTFLAG: --sort-order now rejected ("unexpected argument"). Fixture: reports/fgbio-parity-fixtures/r2-clip-02-03-04.sh (+ before/after evidence). Tests: soft-only helper unit tests (typed + raw), find_primary_pair_indices partition test; full suite 2213 passed; ci-fmt/ci-lint clean.
What & why
clip(and every internal caller ofSamRecordClipper/RawRecordClipper) diverged from fgbio when clipping removes all of a read's aligned bases. This PR closes the gap and, in doing so, ports the full fgbio/htsjdk fidelity for the unmap and mate-info paths.fgbio
SamRecordClipper.clip{Start,End}OfAlignmentcomputesnumberOfClippableBases(the sum of alignmentM/=/Xand insertionICIGAR op lengths) and, when the requested clip is>=that count, calls htsjdkSAMUtils.makeReadUnmappedand returns — the read becomes unmapped with an empty CIGAR, keeping its bases/qualities. fgumi instead clipped every aligned base, leaving the read mapped at a stalePOSwith an all-clip CIGAR (e.g.50H) and an empty sequence — silent corruption. This is tracker finding CLIP-01.Fix
1. Unmap when the whole alignment is clipped (CLIP-01).
numberOfClippableBasesguard to bothSamRecordClipperandRawRecordClipper, inclip_start_of_alignmentandclip_end_of_alignment: whenclippable <= requested, unmap the read and return the clippable count, rather than proceeding to clip.make_read_unmapped[_raw]mirrors htsjdkSAMUtils.makeReadUnmapped: clear reference / pos / mapq / CIGAR / TLEN and the duplicate, secondary, supplementary, and proper-pair flags; setUNMAPPED; reverse-complement bases + reverse quals (then clear the reverse flag) for reverse-strand reads; bases/quals are otherwise preserved.10S10Mclip 10,10Mclip 10) are updated to assert the unmap.2. Reorient strand-sensitive aux tags on unmap (htsjdk
reverseComplement(true)fidelity).htsjdk
makeReadUnmappedreverse-complements a reverse-strand read viareverseComplement(true), which does more thanSEQ/QUAL: it also reverse-complements the base-oriented tags inSAMRecord.TAGS_TO_REVERSE_COMPLEMENT(E2,SQ) and reverses the quality-oriented tags inTAGS_TO_REVERSE(OQ,U2). fgumi'smake_read_unmapped[_raw]previously left these forward-oriented, so a reverse-strand read carrying e.g.OQthat clipping unmapped diverged from fgbio — squarely in the CLIP-01 path.make_read_unmapped(typed) andmake_read_unmapped_rawnow reorient these tags for reverse-strand reads, driven by sharedTAGS_TO_REVERSE_COMPLEMENT/TAGS_TO_REVERSEconstants so the two paths stay in lockstep.Z-string encoding is handled (the sole form these tags take in practice); htsjdk additionally reverses array-encoded forms, but none of these tags is emitted as an array by real tooling. Documented inline.E2/U2/SQto theSamTagcatalog.3. Full mate-info parity, including the one-mapped mate-reverse flag.
Because fgbio
ClipBamthen callsSamPairUtil.setMateInfo(r1, r2, true)on every pair, the clipper unmap surfaced a paired-read gap: fgumi'scliponly refreshed the MC/MQ tags. This PR replaces that withset_mate_info_raw, a faithful port of htsjdksetMateInfo(all three branches — both mapped, both unmapped, one of each) covering mate ref/pos/strand, the mate-unmapped flag, MQ/MC tags, and TLEN. This makes paired output match fgbio when a mate unmaps, and also fixes stale mate-pos/TLEN in the ordinary both-mapped case (previously never updated after a clip shifted a start).SamPairUtil.java:267) rather than an assumedfalse. The clipper's own unmap clearsREVERSE, but a read that arrived unmapped-on-input may still carry it, andset_mate_info_rawruns on every pair, so such a read reaches this branch.Parity evidence (fgbio 4.1.0
ClipBamvsfgumi clip,compare bams --mode content)Reproducer:
reports/fgbio-parity-fixtures/clip-01-unmap-fully-clipped.sh.Scenario A — mapped fragment,
50M@chr1:100,--read-one-five-prime 60:FLAG=0 POS=150 CIGAR=50H(mapped, corrupt)FLAG=4 POS=0 CIGAR=*(unmapped)Scenario B — FR pair, both clipped past alignment (
--read-one-five-prime 101 --read-two-five-prime 101, ClipBamTest.scala:621):MATE_REVERSEset,RNEXT=chr1)FLAG=77/141,RNEXT=* PNEXT=0 TLEN=0Also verified IDENTICAL to fgbio for normal paired clipping that does not unmap (overlap clip, fixed 3′ clip, clip-past-mate), confirming the
setMateInfochange is safe on the common path.Tests / CI
crosscheck_unmap_*): start/end, boundary (clip == clippable), no-unmap regression guard, and negative-strand revcomp — on bothSamRecordClipperandRawRecordClipper.OQ/U2(reversed) andE2(reverse-complemented) identically on the typed and raw paths; positive-strand leaves them untouched.test_set_mate_info_raw_all_branches(all three htsjdk branches) plus a new one-mapped case asserting the mate-reverse flag follows the unmapped read's actual strand.cargo ci-fmt,cargo ci-lint,cargo ci-test(2217 tests) all green.Refs: CLIP-01.