ALT-aware mapping via post-extension liftover groups - #7
Draft
nh13 wants to merge 25 commits into
Draft
Conversation
nh13
force-pushed
the
feat/alt-aware-liftgroup
branch
from
June 25, 2026 01:52
d6e59d4 to
fa30992
Compare
nh13
force-pushed
the
feat/alt-aware-liftgroup
branch
from
June 25, 2026 17:16
f4d71a9 to
ddc1ddf
Compare
Optional ALT-contig awareness. With a .alt file (--alt FILE, or auto-detected <idx>.alt), reads whose confident placement is deflated by a homologous ALT contig recover their MAPQ, while genuinely ambiguous reads (paralogs, segmental duplications) stay at MAPQ 0. Without a .alt file the aligner is byte-for-byte identical to baseline minibwa. Liftover-group design (not the bwa-mem .alt/postalt scheme): ALT copies are homologous, not identical, so each is scored independently and the primary representative is chosen post-extension. - l2bit: load .alt POS+CIGAR into an ALT->primary span-lift (holes + reverse strand); unmapped/decoy lines skipped. - mb_hit_place: per-hit lifted primary placement; grouping key = co-located lifted start, multi-interval (SV-breakpoint aware). - mb_select_sub: survival guard keeps diverged ALT twins through both prune passes. - mb_reconcile_alt: post-extension grouping by lifted placement, best-DP representative, group-scoped subsc/dp_max2/n_sub recompute with query-overlap + transitive guard (paralogs never merge). - pe.c: PE hooks (exclude non-rep twins from pairing; rescue regroup; chimeric demotion guard keyed on lifted placement). - seed.c: ALT-seed -> primary anchor projection for segdup recovery, scoped to primary-absent-from-subsample; per-read scratch on the kalloc arena. - CLI: --alt, --alt-records, --alt-lift-tol. - test/altlg: 13 end-to-end fixtures; baseline byte-identical without .alt.
An .alt CIGAR's structural indels can offset an ALT twin's lifted placement from the read's true primary copy beyond --alt-lift-tol, or drop it into an insertion hole (unliftable). The per-base lift then fails to co-locate the twin, so it survives mb_reconcile_alt as a co-equal representative -- taking the SAM-primary slot and pinning the primary's MAPQ to 0. Fold an ungrouped ALT hit into a non-ALT hit's group when their query spans overlap and the ALT contig's .alt-mapped primary region covers the non-ALT hit's locus (the correspondence the .alt record establishes). Not a categorical discard: an ALT hit with no overlapping primary stays independent, and two primary loci never merge (paralog-safe). Within a group the non-ALT member is the representative. Ground truth (wgsim from GRCh38 primary, ALT-homologous reads): placement accuracy 86.2% -> 91.9%; the confident+correct gap to bwa mem shrinks 7,162 -> 731. Byte-identical without .alt. Test: test/altlg/test-altalt.sh.
PE pairing already chooses the mate-consistent endpoint among near-equal paralog copies, but mb_set_sam_pri marked the SAM primary by per-read DP/hash order -- which can be the other, mate-inconsistent copy. (Masked in the non-ALT path by the chimeric demotion in pe.c; exposed where the ALT paralog-safety guard skips that demotion for distinct lifted groups.) --pe-pair-primary makes mb_set_sam_pri honor the pair-chosen endpoint. Off by default (plain `minibwa mem` stays byte-identical to baseline); one O(1) index check per pair, no measurable runtime change. Ground truth (--alt --pe-pair-primary vs --alt): placement accuracy 91.9% -> 94.1% (bwa mem 94.0%); confident+correct gap to bwa 731 -> 210; confident-but- wrong 11 -> 18 (FDR 0.0003) -- the cost of committing the pair-chosen copy in hard repeats where bwa stays cautious. Test: test/altlg/test-pe-primary.sh.
Replace the MB_NO_ALT_PROJECT / MB_PROJ_TRACE environment seams in seed.c with --dbg-no-alt-proj / --dbg-alt-proj, matching the existing --dbg-* family (kom_dbg_flag); no getenv remains in the aligner. Rename the lettered pe.c "hook" comments to describe what they enforce. Tests pass the flags instead of env vars. No behaviour change; all test/altlg fixtures green, baseline byte-identical.
PE-pair-primary recovers bwa-parity accuracy in ALT-homologous repeats but was opt-in (off by default). Without a .alt it is a provable no-op -- the chimeric demotion in pe.c masks it -- so defaulting it off there cost nothing AND kept plain `minibwa mem` byte-identical to baseline. With a .alt it is the whole accuracy win, so opt-in left it off in exactly the runs that want it. Make it default-on when (and only when) a .alt is loaded: - opt->pe_pair_primary is now tri-state (-1 auto / 0 off / 1 on), default auto. - --pe-pair-primary becomes an optional-argument flag matching --long/--adap: bare or =yes forces on, =no forces off, absent = auto. - Resolved to MB_F_PE_PAIR_PRI after index load: auto => on iff l2b->n_alt > 0 (mb_idx_has_alt), so "off unless .alt" stays the resting state. l2b records the loaded ALT-record count; this covers both --alt FILE and an auto-detected sibling <idx>.alt. Net: `minibwa mem --alt X` now delivers the recovery out of the box; `minibwa mem` is unchanged. Ground truth (--alt, wgsim from GRCh38 primary, ALT-homologous reads): placement accuracy 91.9% -> 94.1% (bwa mem 94.0%); confident+correct gap to bwa 731 -> 210; confident-but-wrong 11 -> 18 (FDR 0.0003). Byte-identical without .alt. Test: test/altlg/test-pe-primary.sh (locks =no -> copyB, auto/explicit -> copyA, chrM byte-identical).
With pair-primary on, the residual confident-but-wrong reads in ALT repeats are "coin-flip" pairs: the anchor end has >=2 co-optimal copies (a primary plus an ALT twin re-exposed as a distinct representative by the liftover grouping), the mate is placed only by rescue, and rescue fits it equally beside several of the anchor's copies. mb_matesw already detects this -- it sets skip[r] when it finds two equal-scoring rescued pairs -- but discards the signal, so the surviving single pair is reported confidently though it is an arbitrary pick of one copy. Surface that signal (rescue_tie[r]) and damp mapq_pe to ~0 when all hold: (a) one endpoint is rescue-only, (b) rescue_tie set for the anchor, and (c) the anchor has >=2 co-optimal representatives (DP second-best within one match score of best). (a)+(b) keep it off genuine recoveries, where the mate fits exactly one copy. Reuses already-computed fields; no new passes. Fires only when ALT lifting has re-exposed the second co-optimal copy as a representative, so it is a no-op without a .alt (byte-identical to baseline; verified 0 of 4,000,000 records change). Ground truth (--alt, on top of default pair-primary): confident-but-wrong 18 -> 8 (ALT FDR 0.00031 -> 0.00014); 16 of 4M correct reads drop below MAPQ 30 (10 of them bwa also places confidently), head-to-head gap to bwa 210 -> 220. No measurable runtime change. The byte-identical-without-.alt invariant is covered by the existing test/altlg suite.
The example programs link libminibwa.a, so they need whatever extra link flags that library was built with. $(LIBS) is empty unless something defines it, so this is a no-op on its own and supplies the flags when a build that does set it is in play -- for instance a parallel index build that pulls in the OpenMP runtime, where the omission fails at link time on Linux but not on macOS.
The coin-flip MAPQ damp and the duplicate-rescued-hit guard both ran on the baseline path. The damp carried an argument that it was a no-op without a .alt rather than a check, and the argument was wrong: two genuine paralogs surviving as representatives with dp_max within opt->a, plus a rescue tie, satisfies every condition with no ALT anywhere. The dedup argued it was 'safe' on the non-ALT path, which is true but is not the byte-identity the manifest promises. Measured on 100k HG002 WGS pairs against an index with no .alt: before this commit the SAM differed from stock r421, after it the two are byte-identical. The chrM fixture the distribution gates on has no ALT contigs and never reaches either site, so the gate passed the whole time the claim was false. Both now test mb_any_alt() on the hits actually in hand, which is the same predicate pe.c already uses for twin exclusion.
nh13
force-pushed
the
feat/alt-aware-liftgroup
branch
from
August 7, 2026 20:03
25ff6e0 to
6e4c3d2
Compare
…thing Three hygiene fixes to the ALT liftover-group feature. No change to mapping behaviour: the no-.alt path is byte-identical to stock, now demonstrably so. Keep the lift bookkeeping off the no-.alt path. mb_select_sub allocated a placement array and called mb_hit_place() for every kept hit whenever `l2b` was non-null -- but l2b is the sequence index and is always present, so every read paid for it. The survival guard it feeds can only fire on an ALT hit, and hit->is_alt is copied from l2b->ctg[tid].is_alt, so with no .alt loaded no hit is ever ALT and the array is written but never read. Count ALT contigs at load time and gate on that instead; skipping the block is provably a no-op. Give the suite a way to run. test/altlg holds 15 scripts that no target and no CI job invoked, and several drive the api-test helpers, which `all` does not build -- so on a fresh checkout they failed for want of a binary. `make test` now builds the helpers and runs the suite. Replace seven vacuous byte-identity checks. Each ran the SAME command twice and compared the output, then reported "hooks gated => inert without .alt". That only shows the aligner is deterministic; it cannot see a hook firing. The real claim needs an unpatched binary, so chrm_baseline() compares against $MB_STOCK and reports itself skipped when unset rather than passing vacuously. It also indexes inside its own tmpdir, which stops the suite leaving .mbw/.l2b beside the shared fixture. Unrelated pre-existing failures: with the helpers built, 11 of 15 scripts fail on this r421 base while all 15 pass on the pre-rebase branch (ddc1ddf, r383) with the same test files. Several failure messages anticipate exactly this ("SA ordering or max_occ changed"), so it is not yet clear whether the rebase regressed the feature or upstream's seeding changes invalidated the fixtures. Untouched here; the four that pass include the new stock comparison, which confirms the no-.alt guarantee.
Eleven of the fifteen ALT scripts failed on r421 with ALT records simply absent from the output. The cause is not the feature: `--outn` and `--alt-records` are wired into the `map` option chain only, and the bwa-compatible `mem` chain silently ignores unknown long options -- no diagnostic, exit status 0. Every script asked for secondary records via `mem --outn=N`, so on r421 they asked for nothing and then failed on the records that were never requested. `minibwa map --outn=5` emits the twin (secondary 0x100, MAPQ 0) alongside the MAPQ 60 primary, which is what the assertions describe. `mem` honoured `--outn` at the branch's original base (r383) and no longer does, which is why this surfaced only after the rebase onto r421. Switch the 42 mapping invocations to `map`. The two subcommands differ only in that `mem` also emits MD:Z, and nothing here asserts on MD. With the api-test helpers built and MB_STOCK pointing at a stock r421 build, the suite goes from 4/15 to 14/15, and the new no-.alt stock comparison passes on every script that runs it -- the PE hooks are byte-identical to stock without a .alt file, as the manifest claims. test-pe-primary is the one remaining failure and is a stale fixture, not a code defect: on r421 copyB comes out as a supplementary rather than a competing full-length hit, and the mate-consistent copy already wins the primary slot without the flag the test exists to exercise.
`make test` builds the api-test helper programs from tracked sources and, in two scripts that index the shared chrM fixture in place, leaves .l2b/.mbw beside it. Neither is source, and both showed up as untracked noise in every worktree that ran the suite.
The flag made the SAM primary follow the PE-pair-chosen endpoint instead of the first-by-DP/hash representative. That was a real fix when it was written: on the pre-r417 base it moved R1's primary from the mate-inconsistent paralog copy to the mate-consistent one, and it was verified to be a complete no-op without a .alt (0 of 400k primary records changed), because the competing copy only survives as a representative when Hook C's r_any_alt-gated paralog-safety skip protects it. Upstream then fixed the same thing itself. 54fb387 ("r417: fixed inconsistent 0x2 flag") sets `h[r]->sam_pri = 1` in the paired branch -- h[r] IS the pair-chosen endpoint, &hit[r][paux.i[r]] -- and sets reset_sam_pri = 0 so the per-read pass cannot override it. That is this flag's behaviour, unconditional and flagless. It has also been inert since the rebase onto r421: pri_idx is assigned only inside the paired branch, and the mb_set_sam_pri call that reads it now sits behind upstream's `if (reset_sam_pri)`, which is false on exactly that path. Both empirical checks agree -- --pe-pair-primary=yes vs =no differs by 0 records on its own fixture and on 100k HG002 pairs against hg38 with the real 3,177-record .alt (ALT machinery confirmed live there: 10,467 records differ with vs without the .alt). So this removes redundancy, not behaviour: SAM output on that same 100k/.alt run is byte-identical before and after. mb_set_sam_pri returns to its upstream 3-argument signature, mb_idx_has_alt goes with it (the auto default was its only caller), and test-pe-primary drops -- it asserted an on/off difference that upstream made impossible. The suite is now 14/14.
minibwa.h included l2bit.h on the line directly above MB_VERSION, which upstream bumps on every release. The two edits therefore landed in one conflict hunk, so assembling this feature onto current master always conflicted in minibwa.h: <<<<<<< master #define MB_VERSION "0.7-r424-dirty" ======= #include "l2bit.h" #define MB_VERSION "0.7-r421" >>>>>>> feat/alt-aware-liftgroup git rerere keys its cache on the conflict pre-image, and that pre-image contained the version string, so a resolution recorded against one release could never match the next one. The cache shows the result: resolved at r416 and r421, unresolved at r419 and r424, one dead entry per release. The feature dropped out of the distribution assembly every time upstream tagged. l2b_t is needed only by the mb_hit_place() prototype 200 lines below, so the include moves down beside the other forward declarations. No semantic change -- the header still compiles standalone and the aligner is untouched.
`make test` failed on any fresh checkout: test-output.sh: Abort trap: 6 minibwa map --outn=5 test/chrM-human.fa.gz ... test-output.sh used the shared chrM fixture without indexing it, and the index is a gitignored build artifact, so on a clean tree there was none. Stock minibwa aborts in that case (`[E::main_map] failed to load the index. ABORT!`), which is what the trap was -- not a fault in this feature. It passed in practice only because test-projrev.sh and test-segdup.sh indexed the shared fixture *in place* as a side effect, leaving untracked .mbw/.l2b in test/. So the suite was order-dependent: it passed on a second run, or whenever another script happened to go first, and failed on a clean tree when test-output.sh led. All three now copy chrM into TMPD and index it there, which is what lib-baseline.sh already documents and does. The suite is hermetic and order-independent, and test/ is left clean.
mb_idx_load() auto-detected <prefix>.alt; mb_idx_load_mmap() did not. The same index and the same command line therefore produced different MAPQs depending on whether --mmap was passed -- on the ALT fixture, 60 without it and 0 with it, because the ALT twin looked like an ordinary repeat when no lift blocks were loaded. Nothing warned; the only symptom was quietly worse mapping quality. Move the resolution out of the loaders and into main_map(), via a new mb_idx_set_alt_auto(). One place decides, so the two loaders cannot disagree, and the precedence between --alt FILE and the adjacent file is visible where it is chosen instead of being an override applied after the fact. That also drops a redundant parse: with --alt FILE, mb_idx_load() used to parse the adjacent .alt only for mb_idx_set_alt() to reset and re-parse it immediately after. The cost is that mb_idx_load() no longer loads .alt for library callers. The two api-test examples that need it already re-loaded it explicitly and now call mb_idx_set_alt_auto(); ex-one and ex-batch inspect no ALT fields. Covered by test/altlg/test-noalt.sh in the following commit, which asserts --mmap parity with ALT both on and off.
ALT-aware mapping turns itself on whenever a .alt sits beside the index, which matches bwa's convention but leaves no way to turn it back off. That matters more here than it does in bwa: bwa only marks ALT hits and defers the real work to bwa-postalt.js, which you run or you don't, whereas this recomputes MAPQ inside the aligner. Anyone on a stock hs38DH index -- which ships a .alt -- got different MAPQs from stock minibwa with nothing in the command line or the log to say why, and no way to check the difference short of moving the file aside. --no-alt loads no .alt at all. is_alt is set only by l2b_set_alt() and is never stored in the index, so loading nothing leaves mb_any_alt(), n_alt_ctg and mb_hit_place() all on their stock paths; there is no new branch in the hot loop and no ALT parse to pay for. It beats an explicit --alt in either order, so it is usable as a late override on a command line assembled elsewhere. test-noalt.sh keys on the fixture's 60-vs-0 MAPQ split and checks that --no-alt does something, that it is byte-identical to deleting the .alt outright, that it wins over --alt in both orders, and that --mmap agrees with the normal loader with ALT on and off -- the regression fixed in the preceding commit.
… --alt Two gaps in the preceding two commits, both found by self-review. Moving .alt resolution out of mb_idx_load() and into the caller updated one caller and missed the other: main_mem() loads the index at map-main.c:709 and never resolved .alt, so the bwa-compat `mem` subcommand silently stopped being ALT-aware. On the ALT fixture it went from MAPQ 60 with two supplementaries to a single MAPQ 0 record. The suite did not catch it because 6ca77f2 moved the ALT tests off `mem` and onto `map` on purpose, leaving that entry point with no coverage at all. Of the seven mb_idx_load*() call sites, this was the only remaining behaviour-bearing one: fastmap is SMEM-only and ex-one/ex-batch read no ALT fields. Separately, mb_idx_set_alt() returned void and dropped l2b_set_alt()'s -1, so a mistyped `--alt /path/typo.alt` loaded nothing and aligned on as if ALT-awareness were off -- no diagnostic even at -v 3, and output indistinguishable from a correct non-ALT run. It now returns the record count like l2b_set_alt(), and an explicitly named file that will not load is an error. The auto-detect branch still ignores its status: an absent adjacent .alt is the normal case, not a mistake, and that asymmetry is the point. test-noalt.sh grows a `mem`-vs-`map` parity check (both with and without the .alt, so it cannot pass for an unrelated reason) and a check that an unreadable --alt reports the failure and emits no alignments. Both were confirmed to fail when their fix is reverted. The bad-path check asserts on stderr rather than exit status because this branch predates r422 (lh3#67), where main() still discarded every subcommand's return value. Also corrects l2b_set_alt()'s header doc, which said it returns the number of ALT contigs; l2bit.c counts records there and keeps the contig count separately in n_alt_ctg.
README told the reader to use a reference without ALT contigs and listed "does not recognize alternate haplotypes" under Limitations, both of which this branch makes false. On the distribution that matters more than it looks: `dist` is the fork's default branch, so this file is the front page every clone lands on, and update_readme() inserts the generated feature banner while deliberately preserving upstream's prose around it. The result was a page advertising alt-liftgroup fifty lines from a paragraph denying it. Corrected here rather than in the renderer so the claim tracks the feature. If alt-liftgroup ever drops from an assembly, this commit drops with it and upstream's original wording -- which is accurate for a build without the feature -- comes back on its own.
Three defects in loading a .alt, all silent. The CIGAR walk treated H as consuming nothing, so a supplementary record whose leading hard clip IS its offset into the ALT contig lifted every block to the wrong place. H removes bases from SEQ but they still occupy the contig, so it advances the ALT axis like S and I. Reverse records walked the reverse complement and stored RC coordinates, which l2b_lift()'s forward binary search then read as forward ones; they are now remapped to ALT-forward. And blocks were appended in file order, which is non-monotonic once reverse records or multi-record contigs exist, so the binary search was operating on unsorted data: each contig's blocks are now sorted, with a total order on (alt_st, alt_en, pri_tid, pri_st, rev) because qsort is unstable and ties would otherwise let libc decide the output. Input from a .alt is untrusted. POS < 1 would underflow to a garbage primary coordinate; an unknown CIGAR operator or a length with no operator silently truncated the walk; and the reverse remap subtracts from ctg->len, so a record whose query span runs past the contig it names wrapped on uint64 to a near-UINT64_MAX coordinate that sorted to the end of lift[] and read back as a valid block. That last one raised MAPQ from 0 to 60 on a 300M record against a 150bp ALT contig -- silent, and in the confidence-inflating direction. All four are now rejected per-record, rolling back that record's partial blocks. Separately, kom_grow() uses realloc(), which does not zero new slots, so the ALT fields of a contig were uninitialized on the plain `minibwa index` path where l2b_set_alt() never runs. l2b_destroy() then freed a garbage ctg->lift, crashing `minibwa index` on any multi-contig reference. test-badalt.sh gates the four rejections. It first proves a well-formed reverse record loads (MAPQ 60) and that a no-.alt run differs (MAPQ 0): without those controls the rejection assertions would pass just as loudly if ALT support were broken outright. The suite's shared helpers (lib.sh, fixlib.sh) land here rather than with the later dedup commit, because this is the first test that sources them and a commit whose own test cannot run is not worth bisecting to.
mb_hit_place/mb_place_t/mb_subpl_t/MB_MAX_SUBPL lived in the public minibwa.h, which dragged the entire private l2bit interface into the public API (via an `#include "l2bit.h"`) and exposed a function that cannot be called through the public surface (mb_idx_t is opaque; the api-test probe already reaches into mbpriv.h for the argument). Move them to mbpriv.h, which already includes both minibwa.h and l2bit.h, and drop the l2bit.h include from the public header. MB_LIFT_TOL stays public (opt->lift_tol references it). Also add the MB_DBG_NO_ALT_SURVIVE debug bit used by the survival-guard test seam.
…on seam The doc described a `pref` parameter the function does not take, and claimed the pairing pass disambiguates the primary here. It does not: selection in mb_set_sam_pri() is per-read, and pe.c sets sam_pri on the pair-chosen endpoint itself (`h[r]->sam_pri = 1`). A doc describing a parameter that does not exist is worse than none. --dbg-no-alt-survive ablates just the mb_select_sub survival guard, so a fixture can assert the RED state -- the co-located ALT twin dropped -- as well as the GREEN one. A guard whose test only ever sees it working cannot tell "the guard holds" from "the fixture never exercises it".
`mem` auto-loaded an adjacent .alt with no way to point at another one or turn it off, so the bwa-compat entry point silently had a different ALT surface from `map`. Both now resolve through one mb_load_alt(), and `mem` accepts --alt, --no-alt, --alt-records and --alt-lift-tol. The resolution is also no longer silent about a user error. An explicitly named --alt that will not load is fatal rather than a quiet fall back to stock behaviour; one that loads zero usable records warns; and --alt together with --no-alt says which one won. An adjacent .alt that is simply absent stays quiet, because that is the normal case rather than a mistake. `mem` carries its own long-option table rather than sharing map's. Sharing it would make `mem` PARSE every knob map has while implementing only these, so one it has no arm for -- --meth, --mmap, --outn, --eqx, --hic -- was accepted and dropped on the floor: `mem --meth` ran without methylation and said nothing. It now rejects unknown options and a missing argument the way map does, and handles --version and --help itself. --alt-records no longer consumes the non-ALT secondary budget: it is emitted as an extra that bypasses the out_n cap, so counting it against n_sec let it crowd out genuine secondaries the user asked for with -h. usage_mem now documents the ALT knobs, and annotates bwa's `-j` as the primitive minibwa replaces rather than leaving it an unexplained ignored flag.
mb_distinct_lifted_group recomputed mb_hit_place(q) for q == h[r] (the pair-chosen endpoint, loop-invariant) on every iteration of the chimeric-hit loop. Add an mb_distinct_lifted_group_pb() variant that takes q's precomputed placement and hoist it out of the loop. Only runs when the mate has ALT hits, so the no-.alt path is unchanged.
…e guard The suite had no fixture for a reverse .alt record or for an ALT contig carrying more than one, which is exactly where the span-lift walked the reverse complement and where blocks arrived out of order -- the two defects the lift fix addresses. test-revmulti.sh covers both. test-survive.sh asserted only that a co-located ALT twin survives sub-hit selection. A test that never sees the guard off cannot distinguish "the guard holds" from "the fixture never reaches it", so it now drives --dbg-no-alt-survive and asserts the twin is dropped without the guard. The helpers were copy-pasted across the suite, in places with incompatible per-file variants of the same name. lib.sh and fixlib.sh hold one version of each, with signatures that are supersets of every prior call site so the optional trailing arguments default to "match anything". Also drops .gitignore patterns for build artifacts the suite no longer produces.
The step-2 comment claimed the transitive close-up is O(n_hit^2), but the all-cross-pairs merge guard makes it O(n_hit^4) in the worst case. Correct the comment: document the real bound, why it stays near-quadratic in practice (n_hit is small under default max_occ/out_n subsampling -- a 40-copy segdup measured n_hit=41; the grp[j]==grp[i] short-circuit plus the fact that co-locating hits of one read also query-overlap means they merge on first contact rather than re-triggering the guard), and why a union-find is not a drop-in replacement (it cannot express the all-cross-pairs guard). No behavior change.
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
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.
ALT-aware mapping via post-extension liftover groups
Draft for design feedback, opened against this fork's
master. It adds optional ALT-contigawareness: given a
.altfile, reads whose confident placement is deflated by a homologous ALTcontig recover their MAPQ, while genuinely ambiguous reads (paralogs, segmental duplications) stay at
MAPQ 0. Without a
.altfile the output is byte-for-byte identical to baseline minibwa — the wholepath is gated on
mb_any_alt.This is the liftover-group design — the same primitive
DRAGEN documents
("liftover groups… each scored taking properly paired alignments into account; the winning group
provides its primary representative as the primary alignment, with MAPQ from the score difference to the
second-best group") — not a port of the bwa-mem
.alt/bwa-postalt.jsscheme.The main question for reviewers is at the bottom: is post-extension lifted-placement grouping the right
primitive?
Problem
GRCh38 ships ALT contigs — alternate versions of regions too variable for one linear sequence (HLA/MHC,
many segmental duplications). A read from such a region scores about equally on the main chromosome and
its ALT copy, so a normal aligner drives it to MAPQ 0 or places it on the ALT contig — lost for
primary-assembly variant calling, in exactly the regions that matter most.
Design (liftover groups, post-extension)
Two principles: (1) ALT copies are homologous, not identical — score each on its own merits and defer
the pick to the end; (2) never lift into the chainer (it breaks
comput_sc's collinearity/monotonicitygates) — all grouping is post-extension, on real DP scores.
Span-lift (
l2bit): load each.altPOS+CIGAR into an ALT→primary coordinate lift (holes andreverse strand handled; unmapped/decoy lines skipped).
Per-hit lifted placement (
mb_hit_place): each hit's primary footprint. Grouping key is theco-located lifted start (same
pri_tid/strand,|Δlifted_st| ≤ tol), multi-interval so anSV-breakpoint-spanning ALT hit still groups via a matching sub-interval — not interval-overlap
fraction (conflates segdup paralogs) and not a scalar k-anchor.
Survival guard (
mb_select_sub): diverged ALT twins survive both prune passes to reachreconciliation.
Reconciliation (
mb_reconcile_alt): group co-located hits; the non-ALT member is therepresentative; recompute group-scoped
subsc/dp_max2/n_subover representatives only, gated onquery-span overlap with an all-cross-pairs transitive guard so genuine paralogs never merge. An
ungrouped ALT hit whose query span overlaps a non-ALT hit inside the ALT contig's own primary region
is folded in too (handles
.alt-indel displacement and unliftable insertion holes).PE integration (
pe.c): non-representative ALT twins are excluded from pair enumeration so thesuboptimal measures the gap to the second-best group-pair; plus a post-rescue regroup and a
cross-group demotion guard keyed on lifted placement.
(Updated: this step used to also add
--pe-pair-primary, making the SAM primary follow thepair-chosen endpoint. Upstream
54fb387"r417: fixed inconsistent 0x2 flag" now does exactly thatunconditionally — it sets
h[r]->sam_pri = 1for the pair-chosen endpoint in the paired branch — sothe flag was redundant and has been removed. SAM output is byte-identical with and without it on
100k HG002 pairs against hg38 with the real
.alt.)ALT-seed → primary projection (
seed.c): when a primary seed is dropped bymax_occsubsamplingin a segdup, inject a native primary anchor at the lifted coordinate so the primary candidate exists
(scoped to "primary absent from the subsample"; chains normally — not a chainer-merge).
Coin-flip MAPQ damp (
pe.c): the residual confident-but-wrong reads are rescue-manufacturedcoin-flips — the anchor end has ≥2 co-optimal copies (a primary plus an ALT twin re-exposed as a
distinct representative by the grouping), and the mate is placed only by rescue, fitting several of
them equally. Mate rescue already detects this (it sets
skip[r]on a rescued-pair tie) but discardedthe signal; surface it (
rescue_tie) and dampmapq_peto ~0 for that exact triple. A no-op withouta
.alt(the second co-optimal copy only survives under lifting) — byte-identical to baseline.CLI:
--alt FILE(auto-detected as<idx>.alt),--alt-lift-tol INT(default 10),--alt-records(emit ALT alignments).
--dbg-alt-proj/--dbg-no-alt-projare diagnostics in the existing--dbg-*family.Results
Ground truth is
wgsim/holodeckreads with known origin; arms arebwa memon hs38DH,minibwa --alton the GATK index +
hs38DH.fa.alt, andminibwano-.alt. Accuracy = primary placed within toleranceof the true locus.
(a) Reads from the primary assembly (
wgsim, ALT-homologous stratum, 66,082 reads):.alt--alt(pair-primary auto-on + coin-flip damp)The follow-up commits took
--altfrom 86.2% → 91.9% → 94.1%; the confident-and-correct gap to bwa(reads bwa places confidently+correctly that minibwa does not) collapsed 7,162 → 731 → 210. The
coin-flip damp then halves the confident-but-wrong tail (18 → 8; ALT FDR 0.00031 → 0.00014) at a
+10-read gap (210 → 220) and no runtime cost.
(b) Reads from the ALT haplotype — the direction the feature most directly targets, with exact,
liftover-noise-free truth (
holodeckfrom the 786 ALT/HLA contigs, 1.86M pairs; each read's originlifted to primary through the
.alt):--alt.altWithout the feature, ALT-haplotype reads are never confidently placed on primary; with it, minibwa
matches bwa. (Both arms show ~9% "wrong" here — a shared scoring artifact from single-point lift across
ALT↔primary structural indels; the minibwa−bwa delta is +0.27 pt.)
Corroboration. Against
bwa mem+bwa-postalt.json 5M PE WGS (hs38DH): 97.6% MAPQ-bucketagreement; the residual is bidirectional and structural, and the simulations above show minibwa is the
more conservative arm on genuine multi-mappers. Real-data DeepVariant concordance vs GIAB HG002 v4.2.1:
MHC F1 +0.39 SNP / +0.33 INDEL, segdup +0.03–0.04, control flat.
Performance & baseline. Byte-identical without
.alt; no-.altCPU equals baseline.--altcosts≈ +5.4% CPU / +12.7% wall (1kg HG00096 5M PE WGS); the pair-primary selection and the coin-flip damp
both reuse already-computed fields — no measurable runtime change (2M-pair WGS: 283s vs 284s CPU). The
ksw2 SW kernels are untouched.
Cost and honest limitations
mate-consistent endpoint, but
mb_set_sam_pripicked the primary by per-read DP/hash order — emittingthe wrong paralog copy, which this feature's paralog guard exposed. That was fixed here behind
--pe-pair-primary; upstream54fb387(r417) then fixed it unconditionally, so the flag has beenremoved and this is no longer a cost or a caveat of the feature.
where the truth copy scores ≤ a paralog and bwa stays at MAPQ 0; the coin-flip damp then claws that
back to +8 by catching the rescue-manufactured ties (16 of 4M correct reads also drop below MAPQ 30, 10
of which bwa places confidently). The base feature lowers FDR (0.00165 → 0.00027); pair-primary +
damp lands at 0.00014 with a +10-read recall cost.
not
max_occ(raising it 20× recovers ~14 for +8% runtime and worsens FDR). This is thesegmental-duplication mapping limit DRAGEN also handles cautiously at the mapper and resolves
downstream. A blunt
n_pp-aware MAPQ cap was prototyped and refuted (it destroys 1,300–2,900 correctrecoveries for ~6 fewer confident-wrong); the coin-flip damp succeeds where it failed because it gates
on the rescue-tie + co-optimal-anchor conjunction, which a genuine recovery does not satisfy.
design — this fixes MAPQ via grouping rather than rewriting coordinates.
--alt-lift-toldefault 10: sweeping 10→100 changed ~50 reads out of 10M and merged no paralogs.Validation suite
test/altlg/— 15 end-to-end fixtures: span-lift, per-hit placement, survival guard, reconciliation(incl. chimera/query-overlap and displaced/hole twins), PE happy/paralog/rescue, PE-pair-primary,
segdup fwd+rev, indel reject,
--alt-records,--alt-lift-tol, unliftable demotion, SV-breakpointgrouping. All green; baseline byte-identical.
For reviewers
for ALT-awareness in minibwa?
.altthe right default for pair-primary (off without.altkeeps baselinebyte-identical), and is the rescue-tie + co-optimal-anchor conjunction the right coin-flip
discriminator for the MAPQ damp?
Commits: (1) the feature; (2) fold displaced/unliftable ALT twins; (3) PE-pair-primary selection;
(4) chore (idiomatic
--dbg-*flags, comment clarity); (5) pair-primary default-on under.alt;(6) coin-flip MAPQ damp. Suggested reading order:
l2bit.c(span-lift) →mb_hit_place/mb_reconcile_altinmap-algo.c→pe.c→seed.c. Thetest/altlg/fixtures double as workedexamples.