From c88769123be3458b52e35ba85c836db4e4b05362 Mon Sep 17 00:00:00 2001 From: minibwa-distro-bot Date: Sun, 23 Aug 2026 00:08:01 +0000 Subject: [PATCH] sync: assemble 0.7-nh13.dev+3d10d42 on upstream f0e1174 --- .gitignore | 5 -- Makefile | 4 + api-test/ex-lift-check.c | 40 ++++++++-- dist-manifest.json | 2 +- l2bit.c | 117 +++++++++++++++++++++++++---- l2bit.h | 1 - map-algo.c | 39 ++++++---- map-main.c | 104 +++++++++++++++++++++---- mbpriv.h | 61 +++++++++++++++ minibwa.h | 79 +------------------ pe.c | 17 +++-- test/altlg/fixlib.sh | 16 ++++ test/altlg/lib.sh | 34 +++++++++ test/altlg/mkfixture-indel.sh | 6 +- test/altlg/mkfixture-paralog.sh | 6 +- test/altlg/mkfixture-place.sh | 2 +- test/altlg/mkfixture-projrev.sh | 11 +-- test/altlg/mkfixture-revmulti.sh | 36 +++++++++ test/altlg/mkfixture-segdup-rev.sh | 11 +-- test/altlg/mkfixture-segdup.sh | 6 +- test/altlg/mkfixture-svbreak.sh | 11 +-- test/altlg/test-altalt.sh | 10 +-- test/altlg/test-badalt.sh | 99 ++++++++++++++++++++++++ test/altlg/test-hitplace.sh | 2 +- test/altlg/test-hookc.sh | 9 +-- test/altlg/test-indel.sh | 3 +- test/altlg/test-lifttol.sh | 7 +- test/altlg/test-noalt.sh | 7 +- test/altlg/test-output.sh | 4 +- test/altlg/test-pe.sh | 7 +- test/altlg/test-projrev.sh | 4 +- test/altlg/test-reconcile.sh | 45 ++++++----- test/altlg/test-revmulti.sh | 55 ++++++++++++++ test/altlg/test-segdup-rev.sh | 3 +- test/altlg/test-segdup.sh | 4 +- test/altlg/test-survive.sh | 27 +++++-- test/altlg/test-svbreak.sh | 11 +-- 37 files changed, 637 insertions(+), 268 deletions(-) create mode 100644 test/altlg/fixlib.sh create mode 100644 test/altlg/lib.sh create mode 100644 test/altlg/mkfixture-revmulti.sh create mode 100755 test/altlg/test-badalt.sh create mode 100644 test/altlg/test-revmulti.sh diff --git a/.gitignore b/.gitignore index f8e05a4..ceb86cd 100644 --- a/.gitignore +++ b/.gitignore @@ -12,8 +12,3 @@ __pycache__/ api-test/ex-*-check api-test/mbmap-batch api-test/mbmap-one - -# Indexes the ALT suite builds beside the shared chrM fixture. test-projrev.sh and -# test-segdup.sh index it in place; the rest index inside their own tmpdir. -test/*.l2b -test/*.mbw diff --git a/Makefile b/Makefile index 1d96164..63b075f 100644 --- a/Makefile +++ b/Makefile @@ -85,6 +85,10 @@ minibwa:libminibwa.a $(MALLOC_O) $(AOBJS) main.o # minibwa to enable the no-.alt byte-identity comparisons; they are skipped without it. test:minibwa $(MAKE) -C api-test + @if [ -z "$$MB_STOCK" ]; then \ + echo "NOTE: MB_STOCK unset -> the no-.alt byte-identity-vs-stock comparisons are SKIPPED."; \ + echo " Set MB_STOCK= to exercise that invariant."; \ + fi @fail=0; for t in test/altlg/test-*.sh; do \ printf '%s: ' "$$t"; \ if sh "$$t" > /tmp/altlg-log.$$$$ 2>&1; then echo PASS; \ diff --git a/api-test/ex-lift-check.c b/api-test/ex-lift-check.c index c8b835b..ffe3141 100644 --- a/api-test/ex-lift-check.c +++ b/api-test/ex-lift-check.c @@ -1,18 +1,46 @@ #include #include +#include #include "minibwa.h" #include "l2bit.h" -/* Loads .l2b, sets .alt, prints lift of (chrP_alt, pos) -> (chrP, pos). - * For pos<300 expect identity; for pos>=305 expect pos-5 (the 5bp insertion). */ +/* Loads .l2b, sets .alt, then either: + * (default, no extra args) prints lift of (chrP_alt, pos) -> (chrP, pos): + * for pos<300 expect identity; for pos>=305 expect pos-5 (the 5bp insertion). + * (check mode) with trailing " " triples, asserts each + * lift's primary position equals (or expected<0 means "in a hole / + * not liftable"); prints one line per triple and exits nonzero on any + * mismatch. Also verifies the contig's lift[] is sorted by alt_st. */ +static int64_t find_tid(const l2b_t *l2b, const char *name){ + int64_t i; for(i=0;i<(int64_t)l2b->n_ctg;i++) if(!strcmp(l2b->ctg[i].name,name)) return i; return -1; +} int main(int argc, char **argv){ - if (argc<2){fprintf(stderr,"usage: ex-lift-check \n");return 2;} + if (argc<2){fprintf(stderr,"usage: ex-lift-check [ ...]\n");return 2;} char buf[1024]; snprintf(buf,sizeof buf,"%s.l2b",argv[1]); l2b_t *l2b=l2b_load(buf); if(!l2b){fprintf(stderr,"no l2b\n");return 2;} snprintf(buf,sizeof buf,"%s.alt",argv[1]); l2b_set_alt(l2b,buf); - /* find chrP_alt tid */ - int64_t at=-1,i; for(i=0;i<(int64_t)l2b->n_ctg;i++) if(!strcmp(l2b->ctg[i].name,"chrP_alt")) at=i; + if (argc>2){ /* check mode: (ctg,pos,expected) triples */ + int rc=0, a; + if ((argc-2)%3!=0){fprintf(stderr,"check mode needs triples\n");return 2;} + for(a=2;a+2ctg[at].is_alt){fprintf(stderr,"%s not flagged ALT\n",ctg);return 1;} + unsigned k; int sorted=1; + for(k=1;kctg[at].n_lift;k++) if(l2b->ctg[at].lift[k].alt_st < l2b->ctg[at].lift[k-1].alt_st) sorted=0; + int64_t pt; uint64_t pp=0; uint8_t rev; + int ok=l2b_lift(l2b,at,pos,&pt,&pp,&rev); + long got = ok? (long)pp : -1; + int pass = (got==exp) && sorted; + printf("lift(%s,%llu)=%ld rev=%d sorted=%d expected=%ld %s\n", + ctg,(unsigned long long)pos,got,ok?rev:-1,sorted,exp,pass?"OK":"MISMATCH"); + if(!pass) rc=1; + } + l2b_destroy(l2b); + return rc; + } + /* default (legacy) mode: chrP_alt at 100 and 400 */ + int64_t at=find_tid(l2b,"chrP_alt"); if(at<0||!l2b->ctg[at].is_alt){fprintf(stderr,"chrP_alt not flagged ALT\n");return 1;} - /* lift two positions via the block list */ int64_t pt1, pt2; uint64_t pp1, pp2; uint8_t rev1, rev2; int ok1=l2b_lift(l2b,at,100,&pt1,&pp1,&rev1); int ok2=l2b_lift(l2b,at,400,&pt2,&pp2,&rev2); diff --git a/dist-manifest.json b/dist-manifest.json index 1841d49..7be5ea3 100644 --- a/dist-manifest.json +++ b/dist-manifest.json @@ -1 +1 @@ -{"base": "f0e117436c28addc359b67123d2353f0d4a1f9e8", "head": "06ab16caf2925c7f0963218293f6877127af420a", "version": "0.7-nh13.dev+06ab16c", "merged": ["ops-distro", "ll-affine-reassoc", "ksw2-extension-kernels", "single-copy-parser", "index-threads-v2", "extd2-avx512", "inline-appenders", "parallel-encode", "meth-cleanups", "meth-sam-tags", "soft-clip-penalty", "submem-ablation", "alt-liftgroup"], "dropped": [], "regressed": []} +{"base": "f0e117436c28addc359b67123d2353f0d4a1f9e8", "head": "3d10d42ddc7c6c195b2a69be15be1ec082b5105a", "version": "0.7-nh13.dev+3d10d42", "merged": ["ops-distro", "ll-affine-reassoc", "ksw2-extension-kernels", "single-copy-parser", "index-threads-v2", "extd2-avx512", "inline-appenders", "parallel-encode", "meth-cleanups", "meth-sam-tags", "soft-clip-penalty", "submem-ablation", "alt-liftgroup"], "dropped": [], "regressed": []} diff --git a/l2bit.c b/l2bit.c index 3334ff2..efc8b59 100644 --- a/l2bit.c +++ b/l2bit.c @@ -1,6 +1,7 @@ #define _POSIX_C_SOURCE 200809L #include #include +#include #include #include "kommon.h" #include "l2bit.h" @@ -146,6 +147,13 @@ static void l2b_add_seq(l2b_t *l2b, uint64_t len, const char *seq, const char *n ctg->comm = comm? kom_strdup(comm) : 0; ctg->len = len; ctg->off = l2b->tot_len; + /* kom_grow() uses realloc(), which does NOT zero new slots, so the ALT fields + * must be initialized explicitly here. Without this, l2b_destroy() frees an + * uninitialized ctg->lift on the plain `minibwa index` path (l2b_set_alt is + * never called there), which crashes on any multi-contig reference. */ + ctg->is_alt = 0; + ctg->n_lift = 0; + ctg->lift = 0; l2b->tot_len += len; m_pac_old = l2b->m_pac; @@ -250,6 +258,25 @@ void l2b_destroy(l2b_t *l2b) * ALT liftover index * ****************************/ +/* Order lift blocks by ALT-forward start (then end) so l2b_lift()'s binary + * search over [alt_st, alt_en) is valid. Blocks are appended in file order, + * which is non-monotonic for reverse records (emitted descending after the + * RC->forward remap) and for ALT contigs carrying multiple .alt records. */ +static int l2b_lift_cmp(const void *a, const void *b) +{ + const l2b_lift_t *x = (const l2b_lift_t*)a, *y = (const l2b_lift_t*)b; + if (x->alt_st != y->alt_st) return x->alt_st < y->alt_st? -1 : 1; + if (x->alt_en != y->alt_en) return x->alt_en < y->alt_en? -1 : 1; + /* qsort is not stable, so equal keys would reorder arbitrarily between libc + * implementations and make the lift -- and therefore the SAM -- differ across + * platforms. Two records CAN tie on the ALT span and differ in target, which + * multi-record ALT contigs make reachable. Break every tie. */ + if (x->pri_tid != y->pri_tid) return x->pri_tid < y->pri_tid? -1 : 1; + if (x->pri_st != y->pri_st) return x->pri_st < y->pri_st? -1 : 1; + if (x->rev != y->rev) return x->rev < y->rev? -1 : 1; + return 0; +} + int l2b_set_alt(l2b_t *l2b, const char *fn) { FILE *fp; @@ -260,7 +287,6 @@ int l2b_set_alt(l2b_t *l2b, const char *fn) uint64_t i; /* Reset: free any existing lift blocks and clear flags. */ - l2b->n_alt = 0; for (i = 0; i < l2b->n_ctg; ++i) { free(l2b->ctg[i].lift); l2b->ctg[i].lift = 0; @@ -277,10 +303,12 @@ int l2b_set_alt(l2b_t *l2b, const char *fn) char *fields[12]; int nf; int64_t alt_tid, pri_tid; + long pos1; uint64_t pri_pos, alt_cursor, pri_cursor; uint32_t flag; uint8_t rev; - uint32_t m_lift; + uint32_t m_lift, n_lift0; + int bad; l2b_ctg_t *ctg; /* Skip SAM header lines. */ @@ -329,26 +357,41 @@ int l2b_set_alt(l2b_t *l2b, const char *fn) flag = (uint32_t)atol(fields[1]); rev = (flag & 0x10) ? 1 : 0; - /* POS is 1-based; convert to 0-based. */ - pri_pos = (uint64_t)(atol(fields[3]) - 1); + /* POS is 1-based; a valid mapped record has POS >= 1. A malformed POS=0 + * (or negative) would underflow to a garbage primary coordinate, so skip. */ + pos1 = atol(fields[3]); + if (pos1 < 1) { + if (kom_verbose >= 2) + fprintf(stderr, "[W::%s] QNAME '%s' has invalid POS '%s', skipping\n", __func__, fields[0], fields[3]); + continue; + } + pri_pos = (uint64_t)(pos1 - 1); ctg = &l2b->ctg[alt_tid]; - /* n_alt counts records (a contig may carry several); n_alt_ctg counts contigs, - * so only the first record for a contig advances it. */ + /* n_alt_ctg counts contigs (a contig may carry several records via + * supplementary lines), so only the first record for a contig advances it. */ if (!ctg->is_alt) ++l2b->n_alt_ctg; ctg->is_alt = 1; ++n_alt; - /* Walk CIGAR to build lift blocks. */ + /* Walk CIGAR to build lift blocks. alt_cursor tracks the position along + * the CIGAR's query axis: for a forward record this IS the ALT-forward + * coordinate; for a reverse record it is the reverse-complement coordinate, + * so each M block is remapped to ALT-forward below. A final sort over all + * of the contig's blocks restores the alt_st ordering l2b_lift() needs. */ alt_cursor = 0; pri_cursor = pri_pos; m_lift = ctg->n_lift; + n_lift0 = ctg->n_lift; /* rollback point if the CIGAR is malformed */ + bad = 0; p = fields[5]; while (*p) { uint64_t len = 0; int op; + if (*p < '0' || *p > '9') { bad = 1; break; } /* expected an op length */ while (*p >= '0' && *p <= '9') len = len * 10 + (*p++ - '0'); + if (*p == '\0') { bad = 1; break; } /* length with no operator (ran off the field) */ op = *p++; if (op == 'M' || op == '=' || op == 'X') { /* Emit one lift block. */ @@ -360,33 +403,75 @@ int l2b_set_alt(l2b_t *l2b, const char *fn) blk.pri_en = pri_cursor + len; blk.rev = rev; - if (ctg->n_lift >= m_lift) { - m_lift = ctg->n_lift + 1; - m_lift += m_lift >> 1; - ctg->lift = (l2b_lift_t*)realloc(ctg->lift, m_lift * sizeof(l2b_lift_t)); + /* Reverse record: the CIGAR walks the reverse-complement of the ALT + * contig, so [alt_st, alt_en) is in RC coordinates. Remap to + * ALT-forward ([len_alt-alt_en, len_alt-alt_st)) while leaving + * pri_st/pri_en unchanged; l2b_lift()'s reverse branch then maps the + * forward query position back with its decreasing correspondence. */ + if (rev) { + /* The remap subtracts from ctg->len, so a record whose query span + * runs past the contig it names would wrap to ~UINT64_MAX and sort + * to the end of lift[], where the binary search reads it as a valid + * block -- silently wrong lifts, no crash. Reject the record. */ + if (blk.alt_en > ctg->len) { bad = 1; break; } + uint64_t fst = ctg->len - blk.alt_en; + uint64_t fen = ctg->len - blk.alt_st; + blk.alt_st = fst; + blk.alt_en = fen; } + + kom_grow(l2b_lift_t, ctg->lift, ctg->n_lift, m_lift); ctg->lift[ctg->n_lift++] = blk; alt_cursor += len; pri_cursor += len; - } else if (op == 'I' || op == 'S') { - /* ALT-only insertion: advance only alt_cursor (hole in primary map). */ + } else if (op == 'I' || op == 'S' || op == 'H') { + /* Query/ALT-consuming ops with no primary footprint. A hard clip + * removes bases from SEQ but they still occupy the ALT contig, so + * (like S and I) it advances the ALT axis -- essential for a + * supplementary .alt record whose leading H is its offset into the + * ALT contig. */ alt_cursor += len; } else if (op == 'D' || op == 'N') { /* Deletion from ALT: advance only pri_cursor. */ pri_cursor += len; + } else if (op == 'P') { + /* Padding: consumes neither coordinate. */ + } else { + bad = 1; break; /* unknown CIGAR operator */ } - /* H (hard clip) and P (padding) consume nothing in either coordinate. */ + } + if (bad) { + /* Roll back this record's partial blocks. The buffer keeps whatever + * capacity the partial append grew it to; shrink so the next record's + * `m_lift = ctg->n_lift` still describes the real allocation, as it does + * on every non-rollback path. */ + ctg->n_lift = n_lift0; + if (ctg->n_lift > 0 && ctg->n_lift < m_lift) + ctg->lift = kom_realloc(l2b_lift_t, ctg->lift, ctg->n_lift); + if (kom_verbose >= 2) + fprintf(stderr, "[W::%s] QNAME '%s' has an unparsable CIGAR '%s', skipping record\n", __func__, fields[0], fields[5]); + continue; } /* Shrink the lift array to exact size. */ if (ctg->n_lift > 0 && ctg->n_lift < m_lift) - ctg->lift = (l2b_lift_t*)realloc(ctg->lift, ctg->n_lift * sizeof(l2b_lift_t)); + ctg->lift = kom_realloc(l2b_lift_t, ctg->lift, ctg->n_lift); } free(line); fclose(fp); - l2b->n_alt = n_alt; + + /* Blocks were appended in file order. Reverse records and multi-record ALT + * contigs make lift[] non-monotonic in alt_st, but l2b_lift()'s binary search + * requires it sorted. Sort each ALT contig's blocks now that every record is + * loaded. */ + for (i = 0; i < l2b->n_ctg; ++i) { + l2b_ctg_t *c = &l2b->ctg[i]; + if (c->n_lift > 1) + qsort(c->lift, c->n_lift, sizeof(l2b_lift_t), l2b_lift_cmp); + } + return n_alt; } diff --git a/l2bit.h b/l2bit.h index c32f4eb..05389db 100644 --- a/l2bit.h +++ b/l2bit.h @@ -38,7 +38,6 @@ typedef struct { uint64_t n_pac, m_pac; uint64_t n_ambi, m_ambi; uint64_t n_mask, m_mask; - uint64_t n_alt; // number of ALT lift records loaded from a .alt (0 => no ALT awareness) l2b_intv_t *ambi, *mask; uint64_t *pac; char *cat_name, *cat_comm; diff --git a/map-algo.c b/map-algo.c index 34d6b66..2260333 100644 --- a/map-algo.c +++ b/map-algo.c @@ -571,15 +571,10 @@ void mb_set_parent(void *km, float mask_level, int mask_len, int n, mb_hit_t *r, /* Mark the SAM primary among the group representatives (parent==id). * - * `pref` is the PE-pair-chosen endpoint (mb_pair sets it to paux.i[r] when a - * proper pair was applied; -1 otherwise / on the SE path). When present and it - * survived as a representative, it IS the read's primary placement: the pair - * score (DP + insert-size consistency) disambiguated near-equal paralog copies - * that this per-read pass cannot (the reps are sorted by DP score then hash, so - * among equal-scoring subtelomeric/segdup paralogs the first-by-index rep is - * effectively arbitrary -- and was emitting a different copy than the mate- - * consistent one the pairing chose). Otherwise fall back to the 5'-most - * (is_primary5) or the first representative. */ + * Selection is per-read: the 5'-most representative (smallest query start) when + * is_primary5 is set, otherwise the first representative by index. The pairing + * pass does not disambiguate the primary here -- upstream sets sam_pri on the + * pair-chosen endpoint separately (see the paired branch in pe.c). */ int32_t mb_set_sam_pri(int32_t n, mb_hit_t *r, int32_t is_primary5) { int32_t i, new_pri, n_pri = 0, min_i = -1, min_qs = -1, first_i = -1; @@ -652,7 +647,10 @@ void mb_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_, * intentionally generous: over-keeping is cheap; the authoritative grouping * happens later. Never make this tolerance tighter than MB_LIFT_TOL. */ if (use_lift) { - for (i = 0; i < n; ++i) { + /* --dbg-no-alt-survive ablates just this guard (testing), so a fixture + * can assert the RED state (twin dropped) as well as the GREEN one. */ + int ablate = (kom_dbg_flag & MB_DBG_NO_ALT_SURVIVE) != 0; + if (!ablate) for (i = 0; i < n; ++i) { if (!keep[i] && r[i].is_alt) { if (mb_place_matches_any(l2b, &r[i], kept_pl, n_kept, lift_tol)) keep[i] = 1; @@ -748,10 +746,23 @@ void mb_reconcile_alt(void *km, const l2b_t *l2b, int n_hit, mb_hit_t *hit, int for (i = 0; i < n_hit; ++i) pl[i] = mb_hit_place(l2b, &hit[i]); - /* 2. group by lifted placement. n_hit is small (a handful), so an O(n^2) - * transitive close-up is cheaper and clearer than a real union-find: assign - * each hit to the lowest-index hit it co-locates with that already has a - * group. Unliftable hits never match and stay singletons. */ + /* 2. group by lifted placement via transitive close-up: for each co-locating + * pair in different groups, unite the groups -- but only after an all-cross- + * pairs guard (the t/u loops below) confirms EVERY member pair co-locates AND + * query-overlaps, which is what stops chimeric segments and paralogs from + * drifting together. Unliftable hits never match and stay singletons. + * + * Complexity: the guard makes a single merge O(n_hit^2), so the loop is + * O(n_hit^4) in the worst case, NOT O(n_hit^2). In practice it stays near + * O(n_hit^2): (a) n_hit is small -- default max_occ/out_n subsampling caps the + * surviving hits, and even a 40-copy segdup measured n_hit=41; (b) the + * `grp[j]==grp[i]` short-circuit skips within-group pairs, and hits of one + * read that co-locate also query-overlap, so they merge on first contact and + * collapse to one group instead of re-triggering the guard. A true union- + * find would drop the exponent but cannot express the all-cross-pairs guard + * (it unites on a single edge), so it would merge groups this must keep apart. + * If a profile ever shows this hot, cap n_hit rather than change the grouping + * semantics. */ for (i = 0; i < n_hit; ++i) grp[i] = i; /* initially own group */ for (i = 0; i < n_hit; ++i) { if (!pl[i].liftable) continue; diff --git a/map-main.c b/map-main.c index 90223d5..79d574f 100644 --- a/map-main.c +++ b/map-main.c @@ -161,7 +161,11 @@ static void worker_for_format(void *data, long i, int tid) } } mb_format(km, out, idx->l2b, t, seg_en - seg_st, &s->n_hit[seg_st], &s->hit[seg_st], j, opt, k - seg_st, mate_qlen); - n_sec += (h->parent != h->id); + /* An --alt-records hit is emitted as an extra and bypasses + * the out_n cap, so it must not consume the non-ALT + * secondary budget either -- otherwise it can crowd out + * genuine secondaries the user asked for with -h/out_n. */ + n_sec += (h->parent != h->id) && !alt_rec; } } } else if (!(opt->flag & MB_F_NO_UNMAP)) { @@ -406,6 +410,7 @@ static ko_longopt_t long_options[] = { { "dbg-an-pos", ko_no_argument, 606 }, // anchor position { "dbg-no-alt-proj", ko_no_argument, 607 }, // ablate ALT->primary projection (testing) { "dbg-alt-proj", ko_no_argument, 608 }, // trace projected primary anchors (testing) + { "dbg-no-alt-survive", ko_no_argument, 609 }, // ablate the mb_select_sub survival guard (testing) { "version", ko_no_argument, 901 }, { "help", ko_no_argument, 902 }, { 0, 0, 0 } @@ -542,6 +547,25 @@ static int32_t parse_meth_tags(const char *spec) return is_excl? (MB_METH_TAG_ALL & ~mask) : mask; } +/* `mem`'s own long options: exactly the knobs it implements. Sharing map's + * table instead would make `mem` PARSE every option map has while handling only + * these, so one it has no arm for -- --meth, --mmap, --outn, --eqx, --hic -- + * would be accepted and silently dropped, and `mem --meth` would run without + * methylation and say nothing. Codes are shared with long_options[] so the two + * parsers cannot drift on what a number means. */ +static ko_longopt_t mem_long_options[] = { + { "alt", ko_required_argument, 316 }, + { "no-alt", ko_no_argument, 319 }, + { "alt-records", ko_no_argument, 317 }, + { "alt-lift-tol", ko_required_argument, 318 }, + { "dbg-no-alt-proj", ko_no_argument, 607 }, + { "dbg-alt-proj", ko_no_argument, 608 }, + { "dbg-no-alt-survive", ko_no_argument, 609 }, + { "version", ko_no_argument, 901 }, + { "help", ko_no_argument, 902 }, + { 0, 0, 0 } +}; + static void set_ins_size(mb_opt_t *opt, const char *arg) { char *q; @@ -553,6 +577,32 @@ static void set_ins_size(mb_opt_t *opt, const char *arg) opt->flag |= MB_F_PE_PREDEF; } +/* Resolve the .alt for a loaded index, shared by `map` and `mem` so their ALT + * surface stays consistent. Returns 0 on success, 1 on a hard error (an + * explicitly-named --alt that will not load). --no-alt loads nothing, leaving + * every ALT code path inert. Emits a warning when --alt is overridden by + * --no-alt, and when an explicitly-named --alt yields zero usable records. */ +static int mb_load_alt(mb_idx_t *idx, const char *prefix, const char *alt_fn, int no_alt) +{ + if (no_alt) { + if (alt_fn) + fprintf(stderr, "[W::%s] --alt '%s' is ignored because --no-alt was also given\n", __func__, alt_fn); + return 0; + } + if (alt_fn) { /* explicitly named: failure to load is a user error */ + int nr = mb_idx_set_alt(idx, alt_fn); + if (nr < 0) { + fprintf(stderr, "[ERROR] failed to load the ALT file '%s'\n", alt_fn); + return 1; + } + if (nr == 0) + fprintf(stderr, "[W::%s] --alt '%s' loaded 0 usable ALT records; running as if without a .alt\n", __func__, alt_fn); + } else { /* auto: an absent adjacent .alt is the normal case, not an error */ + mb_idx_set_alt_auto(idx, prefix); + } + return 0; +} + int main_map(int argc, char *argv[]) { const char *opt_str = "x:o:k:c:m:p:A:B:U:b:O:E:L:t:K:N:PyYR:H:aul:w:W:g:5s:fI:"; @@ -676,6 +726,8 @@ int main_map(int argc, char *argv[]) kom_dbg_flag |= MB_DBG_NO_ALT_PROJ; } else if (c == 608) { // --dbg-alt-proj kom_dbg_flag |= MB_DBG_ALT_PROJ; + } else if (c == 609) { // --dbg-no-alt-survive + kom_dbg_flag |= MB_DBG_NO_ALT_SURVIVE; } else if (c == 'K') { mo.mb_size = mo.max_mb_size = kom_parse_num(o.arg, &s); if (*s == ',') mo.max_mb_size = kom_parse_num(s + 1, &s); @@ -718,16 +770,8 @@ int main_map(int argc, char *argv[]) * normal path cannot disagree about whether this index is ALT-aware. --no-alt * loads nothing at all, which leaves every ALT code path inert -- is_alt is set * only by l2b_set_alt(), never stored in the index. */ - if (!no_alt) { - /* An absent adjacent .alt is the normal case, so the auto branch ignores its - * status. An explicitly named one that will not load is a user error. */ - if (alt_fn) { - if (mb_idx_set_alt(idx, alt_fn) < 0) { - fprintf(stderr, "[ERROR] failed to load the ALT file '%s'\n", alt_fn); - return 1; - } - } else mb_idx_set_alt_auto(idx, argv[o.ind]); - } + if (mb_load_alt(idx, argv[o.ind], alt_fn, no_alt) != 0) + return 1; if (kom_verbose >= 3) fprintf(stderr, "[M::%s::%.3f*%.2f] index loaded\n", __func__, kom_realtime(), kom_percent_cpu()); @@ -775,7 +819,11 @@ static int usage_mem(FILE *fp, const mb_opt_t *opt) fprintf(fp, " -R STR SAM read group line in a format like '@RG\\tID:foo\\tSM:bar' []\n"); fprintf(fp, " -H STR if STR starts with @, insert to header; or insert lines in file STR []\n"); fprintf(fp, " -o FILE output file name [stdout]\n"); - fprintf(fp, " *j treat ALT contigs as part of the primary assembly\n"); + fprintf(fp, " --no-alt ignore .alt; align as if no ALT file exists\n"); + fprintf(fp, " --alt FILE path to the .alt file (default: auto-detected .alt)\n"); + fprintf(fp, " --alt-records emit ALT-contig alignments with full SEQ\n"); + fprintf(fp, " --alt-lift-tol INT bp tolerance for grouping ALT twins by lifted locus [%d]\n", MB_LIFT_TOL); + fprintf(fp, " *j (bwa) treat ALT as primary; minibwa uses the liftover-group scheme above instead\n"); fprintf(fp, " -5 take the alignment with the smallest query position as primary\n"); fprintf(fp, " *q don't modify mapQ of supplementary alignments\n"); fprintf(fp, " *K NUM batch size []\n"); @@ -794,6 +842,8 @@ static int usage_mem(FILE *fp, const mb_opt_t *opt) fprintf(fp, "Notes:\n"); fprintf(fp, " - \"minibwa mem\" aims to match the \"bwa mem\" command-line interface\n"); fprintf(fp, " - '*' options are ignored as they are missing or incompatible with minibwa\n"); + fprintf(fp, " - an adjacent .alt is auto-loaded and applied (liftover-group ALT\n"); + fprintf(fp, " awareness); pass --no-alt to disable it\n"); fprintf(fp, " - minibwa and bwa-mem may output different alignments\n"); return fp == stdout? 0 : 1; } @@ -804,12 +854,19 @@ int main_mem(int argc, char *argv[]) ketopt_t o = KETOPT_INIT; mb_opt_t mo; char *fn_out = 0, *rg_line = 0; + const char *alt_fn = 0; + int32_t no_alt = 0; kstring_t hdr_ins = {0,0,0}, hdr = {0,0,0}; mb_idx_t *idx; mb_opt_init(&mo); mo.flag |= MB_F_WRITE_MD; // bwa-mem always writes MD - while ((c = ketopt(&o, argc, argv, 1, "t:k:w:d:r:y:c:D:W:m:SPA:B:O:E:L:U:x:pR:H:o:j5qK:v:T:h:z:aCVYMuI:", 0)) >= 0) { + /* mem_long_options carries minibwa's ALT knobs (--alt/--no-alt/--alt-records/ + * --alt-lift-tol and the --dbg-* diagnostics), so `mem` exposes the same ALT + * surface as `map` rather than silently auto-loading .alt with no way to + * override or disable it -- but only those, so an option `mem` does not + * implement is rejected below instead of parsed and ignored. */ + while ((c = ketopt(&o, argc, argv, 1, "t:k:w:d:r:y:c:D:W:m:SPA:B:O:E:L:U:x:pR:H:o:j5qK:v:T:h:z:aCVYMuI:", mem_long_options)) >= 0) { // algorithm if (c == 't') mo.n_thread = atoi(o.arg); else if (c == 'k') mo.min_len = atoi(o.arg); @@ -837,12 +894,31 @@ int main_mem(int argc, char *argv[]) else if (c == 'C') mo.flag |= MB_F_COPY_COMMENT; else if (c == 'Y') mo.flag |= MB_F_SUPP_SOFT; else if (c == 'I') set_ins_size(&mo, o.arg); + // ALT liftover-group knobs (shared with `map` via long_options) + else if (c == 316) alt_fn = o.arg; // --alt + else if (c == 317) mo.flag |= MB_F_ALT_RECORDS; // --alt-records + else if (c == 318) { mo.lift_tol = atoi(o.arg); if (mo.lift_tol < 0) mo.lift_tol = 0; } // --alt-lift-tol + else if (c == 319) no_alt = 1; // --no-alt + else if (c == 607) kom_dbg_flag |= MB_DBG_NO_ALT_PROJ; // --dbg-no-alt-proj + else if (c == 608) kom_dbg_flag |= MB_DBG_ALT_PROJ; // --dbg-alt-proj + else if (c == 609) kom_dbg_flag |= MB_DBG_NO_ALT_SURVIVE; // --dbg-no-alt-survive + else if (c == 901) { puts(MB_VERSION); return 0; } + else if (c == 902) return usage_mem(stdout, &mo); + /* Reject what map rejects, rather than swallowing it. */ + else if (c == ':') { + fprintf(stderr, "[ERROR] missing option argument\n"); + return 1; + } else if (c == '?') { + fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.erri]); + return 1; + } } if (argc - o.ind < 2) return usage_mem(stderr, &mo); idx = mb_idx_load(argv[o.ind], !!(mo.flag & MB_F_METH)); kom_assert(idx, "failed to load the index."); - mb_idx_set_alt_auto(idx, argv[o.ind]); /* as main_map does; the loader no longer resolves it */ + if (mb_load_alt(idx, argv[o.ind], alt_fn, no_alt) != 0) + return 1; if (kom_verbose >= 3) fprintf(stderr, "[M::%s::%.3f*%.2f] index loaded\n", __func__, kom_realtime(), kom_percent_cpu()); diff --git a/mbpriv.h b/mbpriv.h index 6887b90..b3aaf55 100644 --- a/mbpriv.h +++ b/mbpriv.h @@ -7,6 +7,59 @@ #include "kommon.h" #include "bseq.h" +/* Maximum number of per-block lifted sub-placements retained in mb_place_t. + * A read footprint (~150 bp pre/post DP) overlaps at most this many .alt lift + * blocks in practice; if MORE blocks overlap (a pathologically fragmented .alt + * CIGAR over the footprint) the FIRST MB_MAX_SUBPL are kept and the rest spill + * (documented in mb_hit_place). Spilling only ever DROPS candidate co-location + * intervals -- it can never invent a spurious match -- so it is conservatively + * safe for paralog isolation. */ +#define MB_MAX_SUBPL 8 + +/* One lifted sub-placement: where a single overlapping .alt lift block maps the + * footprint onto primary coordinates. + * st representative (min) primary coordinate of this block's lifted span + * pri_tid primary contig this block lands on + * rev strand of this block's footprint on primary (.alt block strand XOR h->rev) */ +typedef struct { + int64_t st; + int64_t pri_tid; + uint8_t rev; +} mb_subpl_t; + +/* The lifted PLACEMENT of one hit: the primary footprint it occupies over its + * liftable portion. Computed by mb_hit_place(). + * + * MULTI-INTERVAL placement (SV-breakpoint-aware grouping): instead of collapsing + * every overlapping .alt lift block into a single [lifted_st, lifted_en], the + * placement records ONE sub-placement per overlapping block in subpl[]. A + * breakpoint-spanning ALT hit whose footprint straddles an SV-scale indel then + * exposes BOTH the near-breakpoint primary position AND the far one as separate + * sub-placements, so it can still co-locate with its primary twin via the + * matching sub-interval (mb_places_colocate) instead of being dragged thousands + * of bp away by a min/max collapse. Co-location requires a SHARED sub-interval, + * so two distinct primary loci that happen to land in one inflated span are NOT + * merged (paralog safety). + * + * pri_tid REPRESENTATIVE primary contig (== subpl[0].pri_tid; == h->tid for + * non-ALT hits). Kept for back-compat readers. + * lifted_st REPRESENTATIVE primary coordinate (== subpl[0].st). Back-compat + * grouping key for any reader not yet on the multi-interval API. + * lifted_en max primary coordinate over all sub-placements (cosmetic: nothing + * reads it for grouping decisions). + * rev REPRESENTATIVE strand (== subpl[0].rev; .alt block strand XOR h->rev). + * liftable 1 iff n_subpl >= 1 (at least one block of the footprint lifts); + * 0 iff the ENTIRE footprint falls in holes (ALT-specific -> own group). + * n_subpl number of valid sub-placements (1..MB_MAX_SUBPL; 0 when !liftable). + * subpl the per-block lifted sub-placements (first n_subpl entries valid). */ +typedef struct { + int64_t pri_tid; + int64_t lifted_st, lifted_en; + uint8_t rev, liftable; + int n_subpl; + mb_subpl_t subpl[MB_MAX_SUBPL]; +} mb_place_t; + #define MB_DBG_ALN_SEQ (0x1LL) #define MB_DBG_ANCHOR (0x2LL) #define MB_DBG_SEED (0x4LL) @@ -15,6 +68,7 @@ #define MB_DBG_AN_POS (0x20LL) #define MB_DBG_NO_ALT_PROJ (0x40LL) // --dbg-no-alt-proj: ablate ALT-seed->primary projection (testing) #define MB_DBG_ALT_PROJ (0x80LL) // --dbg-alt-proj: trace projected primary anchors (testing) +#define MB_DBG_NO_ALT_SURVIVE (0x100LL) // --dbg-no-alt-survive: ablate the mb_select_sub survival guard (testing) #define MB_SEED_LONG_JOIN 0x1 #define MB_SEED_IGNORE 0x2 @@ -76,6 +130,13 @@ void mb_sync_hits(void *km, int n_regs, mb_hit_t *regs); void mb_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_, mb_hit_t *r, const l2b_t *l2b, int lift_tol); int mb_any_alt(int n_hit, const mb_hit_t *hit); +/* Compute the lifted PLACEMENT (primary footprint) of one hit. For a non-ALT + * hit this is the identity placement on its own contig. For an ALT hit it lifts + * the aligned footprint (the chain interval [ts,te) pre-DP, or the exact CIGAR + * span post-DP) through the .alt span-lift to primary coordinates. h->p may be + * NULL (pre-DP, coarse) or non-NULL (post-DP, exact); liftable==0 if the whole + * footprint is in holes. See mb_place_t. */ +mb_place_t mb_hit_place(const l2b_t *l2b, const mb_hit_t *h); int mb_places_colocate(const mb_place_t *a, const mb_place_t *b, int lift_tol); void mb_reconcile_alt(void *km, const l2b_t *l2b, int n_hit, mb_hit_t *hit, int sub_diff, int lift_tol); void mb_filter_hits(const mb_opt_t *opt, int qlen, int *n_regs, mb_hit_t *regs); diff --git a/minibwa.h b/minibwa.h index c0ba3d4..2f4abf0 100644 --- a/minibwa.h +++ b/minibwa.h @@ -3,7 +3,7 @@ #include -#define MB_VERSION "0.7-nh13.dev+06ab16c" +#define MB_VERSION "0.7-nh13.dev+3d10d42" #define MB_F_PAF (0x1LL) // output in the PAF format #define MB_F_NO_UNMAP (0x2LL) // output unmapped query sequences @@ -99,14 +99,6 @@ typedef struct { int32_t lift_tol; // ALT liftover-group co-location tolerance in bp (default MB_LIFT_TOL) } mb_opt_t; -/* Deliberately here rather than beside at the top: l2b_t is needed only - * by the mb_hit_place() prototype below, and an include on the line above - * MB_VERSION shares a conflict hunk with every upstream version bump. Because - * git rerere keys its cache on the conflict pre-image, and that pre-image would - * then contain the version string, no recorded resolution could ever be reused - * and this feature dropped out of the assembly once per release. Keep it here. */ -#include "l2bit.h" - struct mb_idx_s; typedef struct mb_idx_s mb_idx_t; @@ -149,59 +141,6 @@ typedef struct { * genuine paralogs even at a looser tolerance. */ #define MB_LIFT_TOL 10 -/* Maximum number of per-block lifted sub-placements retained in mb_place_t. - * A read footprint (~150 bp pre/post DP) overlaps at most this many .alt lift - * blocks in practice; if MORE blocks overlap (a pathologically fragmented .alt - * CIGAR over the footprint) the FIRST MB_MAX_SUBPL are kept and the rest spill - * (documented in mb_hit_place). Spilling only ever DROPS candidate co-location - * intervals -- it can never invent a spurious match -- so it is conservatively - * safe for paralog isolation. */ -#define MB_MAX_SUBPL 8 - -/* One lifted sub-placement: where a single overlapping .alt lift block maps the - * footprint onto primary coordinates. - * st representative (min) primary coordinate of this block's lifted span - * pri_tid primary contig this block lands on - * rev strand of this block's footprint on primary (.alt block strand XOR h->rev) */ -typedef struct { - int64_t st; - int64_t pri_tid; - uint8_t rev; -} mb_subpl_t; - -/* The lifted PLACEMENT of one hit: the primary footprint it occupies over its - * liftable portion. Computed by mb_hit_place(). - * - * MULTI-INTERVAL placement (SV-breakpoint-aware grouping): instead of collapsing - * every overlapping .alt lift block into a single [lifted_st, lifted_en], the - * placement records ONE sub-placement per overlapping block in subpl[]. A - * breakpoint-spanning ALT hit whose footprint straddles an SV-scale indel then - * exposes BOTH the near-breakpoint primary position AND the far one as separate - * sub-placements, so it can still co-locate with its primary twin via the - * matching sub-interval (mb_places_colocate) instead of being dragged thousands - * of bp away by a min/max collapse. Co-location requires a SHARED sub-interval, - * so two distinct primary loci that happen to land in one inflated span are NOT - * merged (paralog safety). - * - * pri_tid REPRESENTATIVE primary contig (== subpl[0].pri_tid; == h->tid for - * non-ALT hits). Kept for back-compat readers. - * lifted_st REPRESENTATIVE primary coordinate (== subpl[0].st). Back-compat - * grouping key for any reader not yet on the multi-interval API. - * lifted_en max primary coordinate over all sub-placements (cosmetic: nothing - * reads it for grouping decisions). - * rev REPRESENTATIVE strand (== subpl[0].rev; .alt block strand XOR h->rev). - * liftable 1 iff n_subpl >= 1 (at least one block of the footprint lifts); - * 0 iff the ENTIRE footprint falls in holes (ALT-specific -> own group). - * n_subpl number of valid sub-placements (1..MB_MAX_SUBPL; 0 when !liftable). - * subpl the per-block lifted sub-placements (first n_subpl entries valid). */ -typedef struct { - int64_t pri_tid; - int64_t lifted_st, lifted_en; - uint8_t rev, liftable; - int n_subpl; - mb_subpl_t subpl[MB_MAX_SUBPL]; -} mb_place_t; - struct mb_tbuf_s; typedef struct mb_tbuf_s mb_tbuf_t; @@ -229,22 +168,6 @@ int mb_idx_set_alt_auto(mb_idx_t *idx, const char *prefix); const char *mb_idx_ctg_name(const mb_idx_t *idx, int32_t tid); int64_t mb_idx_ctg_len(const mb_idx_t *idx, int32_t tid); -/** - * Compute the lifted PLACEMENT (primary footprint) of one hit. - * - * For a non-ALT hit this is the identity placement on its own contig. For an - * ALT hit it lifts the aligned footprint (the chain interval [ts,te) pre-DP, or - * the exact CIGAR span post-DP) through the .alt span-lift to primary - * coordinates, taking the min/max over the LIFTED (primary) outputs so a reverse - * .alt block folds correctly and a footprint end sitting in a hole walks inward - * to the first/last liftable base. See mb_place_t. - * - * @param l2b the span-lift index (must have .alt loaded for ALT hits) - * @param h the hit; h->p may be NULL (pre-DP, coarse) or non-NULL (post-DP, exact) - * @return the placement; liftable==0 if the whole footprint is in holes - */ -mb_place_t mb_hit_place(const l2b_t *l2b, const mb_hit_t *h); - void mb_opt_init(mb_opt_t *opt); int mb_opt_preset(mb_opt_t *opt, const char *preset); diff --git a/pe.c b/pe.c index eef9d2c..5909e80 100644 --- a/pe.c +++ b/pe.c @@ -25,16 +25,19 @@ static inline int mb_insert_dir(const mb_hit_t *h0, const mb_hit_t *h1, int64_t * NOT distinct (same locus). An UNLIFTABLE hit is NOT a protected distinct group: * it returns 0 here so it falls through to normal demotion (preserving baseline * behavior and not shielding spurious hits). */ -static inline int mb_distinct_lifted_group(const l2b_t *l2b, const mb_hit_t *a, const mb_hit_t *b, int lift_tol) +/* As mb_distinct_lifted_group() but with b's placement precomputed, so a caller + * that tests many `a` against a fixed `b` (the pair-chosen endpoint) does not + * recompute mb_hit_place(b) on every iteration. */ +static inline int mb_distinct_lifted_group_pb(const l2b_t *l2b, const mb_hit_t *a, const mb_place_t *pb, int lift_tol) { - mb_place_t pa = mb_hit_place(l2b, a), pb = mb_hit_place(l2b, b); - if (!pa.liftable || !pb.liftable) return 0; /* unliftable -> NOT a protected distinct group */ + mb_place_t pa = mb_hit_place(l2b, a); + if (!pa.liftable || !pb->liftable) return 0; /* unliftable -> NOT a protected distinct group */ /* Distinct (protected) iff they do NOT co-locate: no shared lifted * sub-interval (different pri_tid, opposite strand, or all sub-placements * > lift_tol apart). A breakpoint-spanning ALT hit that shares a * sub-interval with its primary twin is NOT distinct (same locus -> may be * demoted), while two genuine paralog loci stay protected. */ - return !mb_places_colocate(&pa, &pb, lift_tol); + return !mb_places_colocate(&pa, pb, lift_tol); } static inline double mb_pair_score(const mb_hit_t *h0, const mb_hit_t *h1, const mb_pestat_t pes[4], int32_t match_sc) @@ -659,6 +662,10 @@ void mb_pair(void *km, const mb_opt_t *opt, const l2b_t *l2b, int32_t n_hit[2], * this mate actually has ALT hits, so a reference with no .alt loaded is * byte-identical (mb_any_alt is false => no behavior change). */ int r_any_alt = mb_any_alt(n_hit[r], hit[r]); + /* q == h[r] is invariant across the loop, so lift its placement once + * rather than recomputing it (a CIGAR walk + lift searches) per hit. */ + mb_place_t pq; + if (r_any_alt) pq = mb_hit_place(l2b, h[r]); for (i = 0; i < n_hit[r]; ++i) { // handle other chimeric hits const mb_hit_t *q = h[r]; mb_hit_t *p = &hit[r][i]; @@ -673,7 +680,7 @@ void mb_pair(void *km, const mb_opt_t *opt, const l2b_t *l2b, int32_t n_hit[2], * inflate MAPQ, undoing the Task-4 paralog guard for PE. Skip the * demotion across distinct groups. Gated on r_any_alt so the * no-.alt baseline is unchanged. */ - if (r_any_alt && mb_distinct_lifted_group(l2b, p, q, opt->lift_tol)) + if (r_any_alt && mb_distinct_lifted_group_pb(l2b, p, &pq, opt->lift_tol)) continue; if (ol > opt->mask_level * (p->qe - p->qs)) { // if p overlaps with h[r] a lot, make it a secondary hit for (j = 0; j < n_hit[r]; ++j) // FIXME: quadratic time complexity, but almost never an issue on real data diff --git a/test/altlg/fixlib.sh b/test/altlg/fixlib.sh new file mode 100644 index 0000000..2742334 --- /dev/null +++ b/test/altlg/fixlib.sh @@ -0,0 +1,16 @@ +# Shared fixture-builder helpers for the ALT liftover-group suite (sourced by +# mkfixture-*.sh). Previously copy-pasted (gen into 7 scripts, rc into 4). + +# Deterministic random DNA of a given length: gen . +gen() { python3 -c " +import random +random.seed($1) +print(''.join(random.choice('ACGT') for _ in range($2))) +"; } + +# Reverse-complement. With an argument it complements that string; with none it +# acts as a stdin filter. Handles upper and lower case. +rc() { + if [ "$#" -ge 1 ]; then printf '%s' "$1" | rev | tr 'ACGTacgt' 'TGCAtgca' + else rev | tr 'ACGTacgt' 'TGCAtgca'; fi +} diff --git a/test/altlg/lib.sh b/test/altlg/lib.sh new file mode 100644 index 0000000..fbec691 --- /dev/null +++ b/test/altlg/lib.sh @@ -0,0 +1,34 @@ +# Shared helpers for the ALT liftover-group test suite (sourced by test-*.sh). +# +# Every helper here was previously copy-pasted across many test-*.sh scripts, +# sometimes with small and occasionally INCOMPATIBLE per-file variations. The +# signatures below are supersets that are backward-compatible with every prior +# call site: trailing arguments are optional and default to "match anything", so +# both e.g. `primary_mapq f` (any primary) and `primary_mapq f q` (that read's +# primary) work, and both `flag_of f q c` and `flag_of f q c pos` work. + +ok() { echo " ok: $1"; } +fail() { echo "FAIL: $1"; exit 1; } + +# SAM body: drop @ header lines. +sam_body() { mawk '$1!~/^@/{print}' "$1"; } + +# Is bit set in flag ? (mawk has no and(); test by arithmetic.) +has_bit() { mawk -v f="$1" -v b="$2" 'BEGIN{ f=int(f); b=int(b); printf "%d\n", (int(f/b)%2==1)?1:0 }'; } + +# FLAG / MAPQ of the record on contig , optionally at 1-based POS . +flag_of() { mawk -v q="$2" -v c="$3" -v p="${4:-}" '$1==q && $3==c && (p=="" || $4==p){print $2; exit}' "$1"; } +mapq_of() { mawk -v q="$2" -v c="$3" -v p="${4:-}" '$1==q && $3==c && (p=="" || $4==p){print $5; exit}' "$1"; } + +# AS:i: tag of the record on contig . +as_of() { mawk -v q="$2" -v c="$3" '$1==q && $3==c { + for(i=12;i<=NF;i++){ if(substr($i,1,5)=="AS:i:"){ print substr($i,6); exit } } }' "$1"; } + +# FLAG / MAPQ of the PRIMARY record (neither secondary 0x100 nor supplementary +# 0x800). With , restrict to that read; without it, the first primary. +primary_flag() { mawk -v q="${2:-}" '$1!~/^@/ && (q=="" || $1==q){f=int($2); if(int(f/256)%2==0 && int(f/2048)%2==0){print $2; exit}}' "$1"; } +primary_mapq() { mawk -v q="${2:-}" '$1!~/^@/ && (q=="" || $1==q){f=int($2); if(int(f/256)%2==0 && int(f/2048)%2==0){print $5; exit}}' "$1"; } + +# Is there a record on chrP at 1-based POS ? (chrP is the primary contig in +# every fixture that uses this.) +has_pos() { mawk -v p="$2" 'BEGIN{f=0} $1!~/^@/ && $3=="chrP" && $4==p{f=1} END{print f}' "$1"; } diff --git a/test/altlg/mkfixture-indel.sh b/test/altlg/mkfixture-indel.sh index 71482ea..d00fd6f 100755 --- a/test/altlg/mkfixture-indel.sh +++ b/test/altlg/mkfixture-indel.sh @@ -26,13 +26,9 @@ # # Usage: mkfixture-indel.sh set -eu +. "$(dirname "$0")/fixlib.sh" d="$1"; mkdir -p "$d" -gen() { python3 -c " -import random -random.seed($1) -print(''.join(random.choice('ACGT') for _ in range($2))) -"; } LEFT=$(gen 11 75) # 75bp block A (shared ALT<->primary) RIGHT=$(gen 22 75) # 75bp block B (shared ALT<->primary) diff --git a/test/altlg/mkfixture-paralog.sh b/test/altlg/mkfixture-paralog.sh index 6ca4e59..6f2e324 100755 --- a/test/altlg/mkfixture-paralog.sh +++ b/test/altlg/mkfixture-paralog.sh @@ -35,15 +35,11 @@ # distinct competitor, so the multi-mapper verdict holds regardless; the # assertion that proves the boundary is checked by the placement dump. set -eu +. "$(dirname "$0")/fixlib.sh" d="$1"; mkdir -p "$d" # Unique-ish building blocks (LCG PRNG) so each region is internally unique but # COPY1==COPY2 exactly (the paralog). -gen() { python3 -c " -import random -random.seed($1) -print(''.join(random.choice('ACGT') for _ in range($2))) -"; } PAD1=$(gen 1 200) COPY=$(gen 2 150) diff --git a/test/altlg/mkfixture-place.sh b/test/altlg/mkfixture-place.sh index af5fb53..1b8e95c 100755 --- a/test/altlg/mkfixture-place.sh +++ b/test/altlg/mkfixture-place.sh @@ -21,6 +21,7 @@ # taken from the first liftable base. # twin-indel chrP read with a tiny (<=tol) deletion vs chrP_alt: still groups. set -eu +. "$(dirname "$0")/fixlib.sh" d="$1"; mkdir -p "$d" S=$(printf '%s' 'GATCCTAGCATGCTAGGCTAACGTTAGCCGATCGTAGCTAGGCATCGATCGTAGCTAGCTAGGCATCGATTACGATCGGCTAATCGATCGTAGCTGATCGA'\ @@ -31,7 +32,6 @@ S=$(printf '%s' 'GATCCTAGCATGCTAGGCTAACGTTAGCCGATCGTAGCTAGGCATCGATCGTAGCTAGCTAGG L=${#S} # revcomp helper (reads stdin, writes revcomp to stdout) -rc() { rev | tr 'ACGTacgt' 'TGCAtgca'; } # --- contigs --- printf '>chrP\n%s\n' "$S" > "$d/ref.fa" diff --git a/test/altlg/mkfixture-projrev.sh b/test/altlg/mkfixture-projrev.sh index 1e711bc..00462a7 100755 --- a/test/altlg/mkfixture-projrev.sh +++ b/test/altlg/mkfixture-projrev.sh @@ -55,19 +55,10 @@ # # Usage: mkfixture-projrev.sh [N-copies] set -eu +. "$(dirname "$0")/fixlib.sh" d="$1"; mkdir -p "$d" N="${2:-20}" -gen() { python3 -c " -import random -random.seed($1) -print(''.join(random.choice('ACGT') for _ in range($2))) -"; } -rc() { python3 -c " -import sys -s=sys.argv[1] -print(s.translate(str.maketrans('ACGT','TGCA'))[::-1]) -" "$1"; } # Sequence seeds chosen so the SA-array lexicographic order places the # ALT-contig hit BEFORE the primary LIFT hit in the 150bp SMEM's SA interval, diff --git a/test/altlg/mkfixture-revmulti.sh b/test/altlg/mkfixture-revmulti.sh new file mode 100644 index 0000000..822730e --- /dev/null +++ b/test/altlg/mkfixture-revmulti.sh @@ -0,0 +1,36 @@ +#!/bin/sh +# Fixture for the reverse-multi-block and multi-record/hard-clip .alt lift paths +# (regression for the l2b_set_alt / l2b_lift coordinate bugs). Every OTHER +# fixture emits a single full-length M block per ALT contig, which is exactly the +# coincidentally-correct regime; these two contigs are the cases that regime hid: +# +# chrPrev_alt (150bp) : reverse record (FLAG 16) with an internal deletion, +# CIGAR 75M2000D75M @ chrP POS 201. A reverse record's +# CIGAR walks the reverse-complement of the ALT contig, so +# the block coordinates need an RC->forward remap; before +# the fix the two blocks were swapped. +# chrPhc_alt (100bp) : two records for one ALT contig — a primary 50M @ POS +# 1001 (FLAG 0) and a supplementary 50H50M @ POS 2001 +# (FLAG 2048). The leading hard clip is the offset into +# the ALT contig; before the fix H consumed nothing and +# the array was left unsorted, so the supplementary span +# lifted to the wrong place (or read as a hole). +# +# Emits ref.fa + ref.fa.alt only; the test drives ex-lift-check in check mode. +set -eu +. "$(dirname "$0")/fixlib.sh" +d="$1"; mkdir -p "$d" + +# Sequence content is irrelevant to liftover (pure coordinate arithmetic from the +# .alt CIGARs); gen comes from fixlib.sh. +{ + printf '>chrP\n%s\n' "$(gen 71 3000)" + printf '>chrPrev_alt\n%s\n' "$(gen 72 150)" + printf '>chrPhc_alt\n%s\n' "$(gen 73 100)" +} > "$d/ref.fa" + +{ + printf 'chrPrev_alt\t16\tchrP\t201\t60\t75M2000D75M\t*\t0\t0\t*\t*\n' + printf 'chrPhc_alt\t0\tchrP\t1001\t60\t50M\t*\t0\t0\t*\t*\n' + printf 'chrPhc_alt\t2048\tchrP\t2001\t60\t50H50M\t*\t0\t0\t*\t*\n' +} > "$d/ref.fa.alt" diff --git a/test/altlg/mkfixture-segdup-rev.sh b/test/altlg/mkfixture-segdup-rev.sh index 4b51b24..8f29a67 100755 --- a/test/altlg/mkfixture-segdup-rev.sh +++ b/test/altlg/mkfixture-segdup-rev.sh @@ -45,19 +45,10 @@ # # Usage: mkfixture-segdup-rev.sh [N-copies] set -eu +. "$(dirname "$0")/fixlib.sh" d="$1"; mkdir -p "$d" N="${2:-20}" -gen() { python3 -c " -import random -random.seed($1) -print(''.join(random.choice('ACGT') for _ in range($2))) -"; } -rc() { python3 -c " -import sys -s=sys.argv[1] -print(s.translate(str.maketrans('ACGT','TGCA'))[::-1]) -" "$1"; } CORE=$(gen 100 100) # 100bp shared seed core (segdup): on every copy + ALT ALTTAIL=$(gen 777 50) # 50bp tail UNIQUE to the ALT contig (makes ALT SMEM unique) diff --git a/test/altlg/mkfixture-segdup.sh b/test/altlg/mkfixture-segdup.sh index 8beff85..354d651 100755 --- a/test/altlg/mkfixture-segdup.sh +++ b/test/altlg/mkfixture-segdup.sh @@ -37,14 +37,10 @@ # # Usage: mkfixture-segdup.sh [N-copies] set -eu +. "$(dirname "$0")/fixlib.sh" d="$1"; mkdir -p "$d" N="${2:-20}" # number of primary copies of CORE (>= the max_occ used in tests) -gen() { python3 -c " -import random -random.seed($1) -print(''.join(random.choice('ACGT') for _ in range($2))) -"; } CORE=$(gen 100 100) # 100bp shared seed core (segdup): on every copy + ALT ALTTAIL=$(gen 777 50) # 50bp tail UNIQUE to the ALT contig (makes ALT SMEM unique) diff --git a/test/altlg/mkfixture-svbreak.sh b/test/altlg/mkfixture-svbreak.sh index a6ee4dc..46462d6 100755 --- a/test/altlg/mkfixture-svbreak.sh +++ b/test/altlg/mkfixture-svbreak.sh @@ -37,18 +37,9 @@ # # Usage: mkfixture-svbreak.sh set -eu +. "$(dirname "$0")/fixlib.sh" d="$1"; mkdir -p "$d" -gen() { python3 -c " -import random -random.seed($1) -print(''.join(random.choice('ACGT') for _ in range($2))) -"; } -rc() { python3 -c " -import sys -s=sys.argv[1] -print(s.translate(str.maketrans('ACGT','TGCA'))[::-1]) -" "$1"; } A=75 # LEFT block length B=75 # RIGHT block length diff --git a/test/altlg/test-altalt.sh b/test/altlg/test-altalt.sh index 16ab39d..d6f8024 100644 --- a/test/altlg/test-altalt.sh +++ b/test/altlg/test-altalt.sh @@ -16,6 +16,7 @@ # # Usage: test/altlg/test-altalt.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -28,18 +29,9 @@ MK="$MDIR/test/altlg/mkfixture-altalt.sh" TMPD=$(mktemp -d /tmp/altlg-altalt.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } -has_bit() { - mawk -v f="$1" -v b="$2" 'BEGIN{ f=int(f); b=int(b); - printf "%d\n", (int(f / b) % 2 == 1) ? 1 : 0; }' -} # primary record (no 0x100, no 0x800) fields for a qname -primary_flag() { mawk -v q="$2" '$1==q{f=int($2); if(int(f/256)%2==0 && int(f/2048)%2==0){print $2; exit}}' "$1"; } primary_rname() { mawk -v q="$2" '$1==q{f=int($2); if(int(f/256)%2==0 && int(f/2048)%2==0){print $3; exit}}' "$1"; } -primary_mapq() { mawk -v q="$2" '$1==q{f=int($2); if(int(f/256)%2==0 && int(f/2048)%2==0){print $5; exit}}' "$1"; } -flag_of() { mawk -v q="$2" -v c="$3" '$1==q && $3==c {print $2; exit}' "$1"; } # group-check field for the ALT contig hit of a qname (is_alt=1 row on ) grp_field() { diff --git a/test/altlg/test-badalt.sh b/test/altlg/test-badalt.sh new file mode 100755 index 0000000..5934263 --- /dev/null +++ b/test/altlg/test-badalt.sh @@ -0,0 +1,99 @@ +#!/bin/sh +# Malformed .alt records must be rejected, not half-loaded. +# +# l2b_set_alt() parses .alt lines out of a user-supplied SAM. Every coordinate it +# derives is attacker- or accident-controlled, and one of them is a SUBTRACTION: +# a reverse record's ALT span is remapped to forward as `ctg->len - alt_en`. If +# the record's CIGAR query span runs past the ALT contig it names -- a .alt built +# against a different reference revision is enough -- that wraps on uint64 to a +# near-UINT64_MAX coordinate, the block sorts to the end of lift[], and +# l2b_lift()'s binary search reads it as valid. The failure is silent and it +# inflates confidence: a garbage lift makes an ALT twin look co-located with its +# primary, and MAPQ goes UP. +# +# Verifies: +# 1. Over-long reverse record -> rejected (MAPQ not inflated by a garbage lift). +# 2. POS < 1 -> rejected. +# 3. Unknown CIGAR operator -> rejected. +# 4. A WELL-FORMED reverse record still loads -- the guards must reject only +# what is actually malformed, or they would silently disable ALT support. +# +# Usage: test/altlg/test-badalt.sh [] +set -eu +. "$(dirname "$0")/lib.sh" +. "$(dirname "$0")/fixlib.sh" + +MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" +MINIBWA="$MDIR/minibwa" + +TMPD=$(mktemp -d /tmp/altlg-badalt.XXXXXX) +trap 'rm -rf "$TMPD"' EXIT + +# chrP (400bp primary) + chrP_alt, the reverse complement of chrP[100,250). +# A reverse .alt record is the one that exercises the subtraction. +PRI=$(gen 7 400) +ALT=$(rc "$(printf '%s' "$PRI" | cut -c101-250)") +printf '>chrP\n%s\n>chrP_alt\n%s\n' "$PRI" "$ALT" > "$TMPD/ref.fa" +printf '@r1\n%s\n+\n%s\n' \ + "$(printf '%s' "$PRI" | cut -c121-220)" "$(printf '%100s' '' | tr ' ' 'I')" > "$TMPD/reads.fq" +"$MINIBWA" index "$TMPD/ref.fa" 2>/dev/null || fail "indexing the fixture failed" + +# Run with the given .alt content and report the primary's MAPQ. +mapq_with() { + printf '%s' "$1" > "$TMPD/ref.fa.alt" + "$MINIBWA" map --outn=5 "$TMPD/ref.fa" "$TMPD/reads.fq" 2>/dev/null > "$TMPD/out.sam" \ + || fail "minibwa map failed on .alt <$1>" + [ -s "$TMPD/out.sam" ] || fail "empty output for .alt <$1>" + primary_mapq "$TMPD/out.sam" +} + +# ============================================================ +# 4 first: establish that a well-formed reverse record DOES load. +# Without this the rejection checks below could all pass for the +# trivial reason that ALT support is broken outright. +# ============================================================ +echo "== well-formed reverse record (control) ==" +good=$(mapq_with 'chrP_alt 16 chrP 101 60 150M * 0 0 * * +') +[ "$good" = "60" ] \ + || fail "well-formed reverse .alt: primary MAPQ is '${good:-}', expected 60 -- ALT support is not working, so the rejection checks below would be vacuous" +ok "well-formed reverse record loads: primary MAPQ $good" + +# The same read with NO .alt at all: the baseline a rejected record must match. +rm -f "$TMPD/ref.fa.alt" +"$MINIBWA" map --outn=5 "$TMPD/ref.fa" "$TMPD/reads.fq" 2>/dev/null > "$TMPD/none.sam" +none=$(primary_mapq "$TMPD/none.sam") +[ "$none" != "$good" ] \ + || fail "no-.alt MAPQ ($none) equals the loaded one ($good); this fixture cannot tell the two apart" +ok "no .alt at all: primary MAPQ $none (the baseline a rejected record must match)" + +# ============================================================ +# 1. Over-long reverse record: 300M on a 150bp ALT contig. +# ============================================================ +echo "== malformed records are rejected ==" +span=$(mapq_with 'chrP_alt 16 chrP 101 60 300M * 0 0 * * +') +[ "$span" = "$none" ] \ + || fail "over-long reverse record: primary MAPQ is '$span', expected '$none' -- the span was accepted and a wrapped coordinate lifted it" +ok "reverse record longer than its ALT contig is rejected (MAPQ $span, not inflated)" + +# ============================================================ +# 2. POS < 1 (a mapped record must be 1-based >= 1). +# ============================================================ +pos=$(mapq_with 'chrP_alt 16 chrP 0 60 150M * 0 0 * * +') +[ "$pos" = "$none" ] \ + || fail "POS=0 record: primary MAPQ is '$pos', expected '$none'" +ok "POS < 1 is rejected (MAPQ $pos)" + +# ============================================================ +# 3. Unknown CIGAR operator. +# ============================================================ +cig=$(mapq_with 'chrP_alt 16 chrP 101 60 100M50Z * 0 0 * * +') +[ "$cig" = "$none" ] \ + || fail "unknown CIGAR op: primary MAPQ is '$cig', expected '$none'" +ok "unknown CIGAR operator is rejected (MAPQ $cig)" + +echo "[test-badalt] PASS" +exit 0 diff --git a/test/altlg/test-hitplace.sh b/test/altlg/test-hitplace.sh index 841f99b..ae9911b 100755 --- a/test/altlg/test-hitplace.sh +++ b/test/altlg/test-hitplace.sh @@ -15,6 +15,7 @@ # # Usage: test/altlg/test-hitplace.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MKFIXTURE="$MDIR/test/altlg/mkfixture-place.sh" @@ -40,7 +41,6 @@ cat "$TMPD/place.txt" echo "----- coarse placement -----" cat "$TMPD/place-coarse.txt" -fail() { echo "FAIL: $1"; exit 1; } # field extractor: print value of key= on the line matching [rev=] # usage: field (first matching hit) diff --git a/test/altlg/test-hookc.sh b/test/altlg/test-hookc.sh index 8888be6..51adfef 100755 --- a/test/altlg/test-hookc.sh +++ b/test/altlg/test-hookc.sh @@ -32,6 +32,7 @@ # # Usage: test/altlg/test-hookc.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -44,15 +45,7 @@ EX_PLACE="$MDIR/api-test/ex-place-check" TMPD=$(mktemp -d /tmp/altlg-hookc.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } -has_bit() { mawk -v f="$1" -v b="$2" 'BEGIN{ f=int(f); b=int(b); - printf "%d\n", (int(f/b)%2==1)?1:0; }'; } -flag_of() { mawk -v q="$2" -v c="$3" -v p="${4:-}" \ - '$1==q && $3==c && (p=="" || $4==p){print $2; exit}' "$1"; } -mapq_of() { mawk -v q="$2" -v c="$3" -v p="${4:-}" \ - '$1==q && $3==c && (p=="" || $4==p){print $5; exit}' "$1"; } # ========================================================================= echo "[test-hookc] building fixture ..." diff --git a/test/altlg/test-indel.sh b/test/altlg/test-indel.sh index c68ba8f..4200531 100755 --- a/test/altlg/test-indel.sh +++ b/test/altlg/test-indel.sh @@ -20,6 +20,7 @@ # # Usage: test/altlg/test-indel.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -28,8 +29,6 @@ MK="$MDIR/test/altlg/mkfixture-indel.sh" TMPD=$(mktemp -d /tmp/altlg-indel.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } echo "[test-indel] building indel fixture (.alt 75M5D75M) ..." /bin/sh "$MK" "$TMPD" 2>/dev/null diff --git a/test/altlg/test-lifttol.sh b/test/altlg/test-lifttol.sh index f6aba21..e4c785c 100755 --- a/test/altlg/test-lifttol.sh +++ b/test/altlg/test-lifttol.sh @@ -10,6 +10,7 @@ # # Usage: test/altlg/test-lifttol.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -17,14 +18,8 @@ MK="$MDIR/test/altlg/mkfixture-lifttol.sh" TMPD=$(mktemp -d /tmp/altlg-lifttol.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # primary MAPQ of : the record with neither 0x100 nor 0x800 set. -primary_mapq() { - mawk -v q="$2" '$1==q { f=int($2); - if (int(f/256)%2==0 && int(f/2048)%2==0) { print $5; exit } }' "$1" -} echo "[test-lifttol] building fixture ..." /bin/sh "$MK" "$TMPD" 2>/dev/null diff --git a/test/altlg/test-noalt.sh b/test/altlg/test-noalt.sh index 49237fe..b305dc6 100755 --- a/test/altlg/test-noalt.sh +++ b/test/altlg/test-noalt.sh @@ -26,6 +26,7 @@ # # Usage: test/altlg/test-noalt.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -33,18 +34,12 @@ MINIBWA="$MDIR/minibwa" TMPD=$(mktemp -d /tmp/altlg-noalt.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # Alignment records only: @PG carries the command line, which legitimately # differs between two runs we expect to align identically. -sam_body() { mawk '$1!~/^@/{print}' "$1"; } # MAPQ of the primary (neither secondary 0x100 nor supplementary 0x800). # mawk has no and(), so test the two bits by arithmetic. -primary_mapq() { - mawk '$1!~/^@/ && int($2/256)%2==0 && int($2/2048)%2==0 {print $5; exit}' "$1" -} # Run `map` with the fixture and the given extra flags, into $TMPD/$1.sam. run() { diff --git a/test/altlg/test-output.sh b/test/altlg/test-output.sh index 2dded10..889353d 100755 --- a/test/altlg/test-output.sh +++ b/test/altlg/test-output.sh @@ -11,6 +11,7 @@ # # Usage: test/altlg/test-output.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -18,8 +19,6 @@ MINIBWA="$MDIR/minibwa" TMPD=$(mktemp -d /tmp/altlg-output.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # ============================================================ # Build a minimal fixture: @@ -167,7 +166,6 @@ CHRM_FA="$MDIR/test/chrM-human.fa.gz" CHRM_R1="$MDIR/test/chrM-read_1.fa.gz" CHRM_R2="$MDIR/test/chrM-read_2.fa.gz" # Filter helper: strip header lines for body-only comparison -sam_body() { mawk '$1!~/^@/{print}' "$1"; } if [ -f "$CHRM_FA" ] && [ -f "$CHRM_R1" ]; then # Copy and index inside TMPD rather than beside the shared fixture (same reason as diff --git a/test/altlg/test-pe.sh b/test/altlg/test-pe.sh index 61413c9..98cff32 100755 --- a/test/altlg/test-pe.sh +++ b/test/altlg/test-pe.sh @@ -34,6 +34,7 @@ # # Usage: test/altlg/test-pe.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -47,15 +48,9 @@ MK_RESC="$MDIR/test/altlg/mkfixture-pe-rescue.sh" TMPD=$(mktemp -d /tmp/altlg-pe.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # has_bit -> "1" if (flag & bit) else "0" (single power-of-two bit). -has_bit() { mawk -v f="$1" -v b="$2" 'BEGIN{ f=int(f); b=int(b); - printf "%d\n", (int(f/b)%2==1)?1:0; }'; } # flag_of -> FLAG of the matching record (pos optional). -flag_of() { mawk -v q="$2" -v c="$3" -v p="${4:-}" '$1==q && $3==c && (p=="" || $4==p){print $2; exit}' "$1"; } -mapq_of() { mawk -v q="$2" -v c="$3" -v p="${4:-}" '$1==q && $3==c && (p=="" || $4==p){print $5; exit}' "$1"; } # proper-pair primary FLAG for a given mate-bit (0x40 first / 0x80 last). prim_flag_mate() { mawk -v q="$2" -v mb="$3" '$1==q { f=int($2); if (int(f/256)%2==0 && int(f/2048)%2==0 && int(f/mb)%2==1) { print $2; exit } }' "$1"; } diff --git a/test/altlg/test-projrev.sh b/test/altlg/test-projrev.sh index 2ff8a63..539a2d5 100755 --- a/test/altlg/test-projrev.sh +++ b/test/altlg/test-projrev.sh @@ -37,6 +37,7 @@ # # Usage: test/altlg/test-projrev.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -51,11 +52,8 @@ N=20 TMPD=$(mktemp -d /tmp/altlg-projrev.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # chrP POS lookup helpers (same style as test-segdup.sh). -has_pos() { mawk -v p="$2" 'BEGIN{f=0} $1!~/^@/ && $3=="chrP" && $4==p{f=1} END{print f}' "$1"; } mapq_at_pos() { mawk -v p="$2" '$1!~/^@/ && $3=="chrP" && $4==p{print $5; exit}' "$1"; } flag_at_pos() { mawk -v p="$2" '$1!~/^@/ && $3=="chrP" && $4==p{print $2; exit}' "$1"; } diff --git a/test/altlg/test-reconcile.sh b/test/altlg/test-reconcile.sh index 3ab25ea..8479757 100755 --- a/test/altlg/test-reconcile.sh +++ b/test/altlg/test-reconcile.sh @@ -33,6 +33,7 @@ # # Usage: test/altlg/test-reconcile.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -48,28 +49,12 @@ MK_CHIM="$MDIR/test/altlg/mkfixture-chimera.sh" TMPD=$(mktemp -d /tmp/altlg-reconcile.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # flag_of -> FLAG of the first matching record -flag_of() { mawk -v q="$2" -v c="$3" '$1==q && $3==c {print $2; exit}' "$1"; } -mapq_of() { mawk -v q="$2" -v c="$3" '$1==q && $3==c {print $5; exit}' "$1"; } # has_bit -> "1" if (flag & bit) else "0". must be a single # power-of-two (FLAG bits via integer arithmetic; mawk lacks a portable &). -has_bit() { - mawk -v f="$1" -v b="$2" 'BEGIN{ f=int(f); b=int(b); - printf "%d\n", (int(f / b) % 2 == 1) ? 1 : 0; }' -} # primary_flag -> FLAG of the record without 0x100 (secondary) # and without 0x800 (supplementary) — i.e. the SAM primary line. -primary_flag() { - mawk -v q="$2" '$1==q { f=int($2); - if (int(f/256)%2==0 && int(f/2048)%2==0) { print $2; exit } }' "$1" -} -primary_mapq() { - mawk -v q="$2" '$1==q { f=int($2); - if (int(f/256)%2==0 && int(f/2048)%2==0) { print $5; exit } }' "$1" -} # ========================================================================= echo "[test-reconcile] building reconcile fixture ..." @@ -109,8 +94,6 @@ fi # --- case (c): homologous-divergent twin --- echo "== case (c): divergent twin does not zero the MAPQ ==" # as_of -> integer of the AS:i: tag of the first match -as_of() { mawk -v q="$2" -v c="$3" '$1==q && $3==c { - for(i=12;i<=NF;i++){ if(substr($i,1,5)=="AS:i:"){ print substr($i,6); exit } } }' "$1"; } cc_flag=$(flag_of "$TMPD/rec.sam" r-diverge chrP) cc_mapq=$(mapq_of "$TMPD/rec.sam" r-diverge chrP) [ -n "$cc_flag" ] || fail "(c) no chrP record for r-diverge" @@ -167,6 +150,32 @@ pb_mapq=$(primary_mapq "$PARD/par.sam" r-para) [ "$pb_mapq" = "0" ] || fail "(b) r-para primary MAPQ=$pb_mapq; expected 0 (paralogs must not merge)" ok "(b) r-para primary MAPQ=$pb_mapq (low-MAPQ multi-mapper; groups not merged)" +# The MAPQ==0 check above is necessary but NOT sufficient: COPY1 and COPY2 are two +# distinct non-ALT primary loci, so they force MAPQ 0 regardless of what +# reconciliation does with altB — a wrongful merge of altB into COPY1's group +# would still leave two co-scoring groups and MAPQ 0, so it passes either way. +# Probe the grouping DIRECTLY: altB is 2*tol from COPY1 and MUST remain its own +# group representative (parent == id); a wrongful cross-tolerance merge would set +# its parent to COPY1's id. This is the assertion that actually exercises the +# paralog guard under test. +[ -x "$EX_GROUP" ] || fail "(b) ex-group-check probe not built ($EX_GROUP); run 'make -C api-test'" +echo "[test-reconcile] paralog grouping dump ..." +"$EX_GROUP" "$PARD/ref.fa" "$PARD/reads.fq" 2>/dev/null > "$PARD/grp.txt" +[ -s "$PARD/grp.txt" ] || fail "(b) ex-group-check emitted no grouping rows" +cat "$PARD/grp.txt" +grp_field() { # -> value of key= on the first row for that ctg + mawk -v c="$2" -v k="$3" ' + { ctg=""; want=""; for(i=1;i<=NF;i++){ n=index($i,"="); + if(n>0){ key=substr($i,1,n-1); val=substr($i,n+1); + if(key=="ctg") ctg=val; if(key==k) want=val; } } + if(ctg==c){ print want; exit } }' "$1" +} +b_id=$(grp_field "$PARD/grp.txt" chrP_altB id) +b_parent=$(grp_field "$PARD/grp.txt" chrP_altB parent) +[ -n "$b_id" ] && [ -n "$b_parent" ] || fail "(b) chrP_altB not present in grouping dump" +[ "$b_parent" = "$b_id" ] || fail "(b) chrP_altB merged into another group (parent=$b_parent != id=$b_id); paralog twin 2*tol away must NOT merge" +ok "(b) chrP_altB is its own group rep (parent=$b_parent == id=$b_id): paralog twin not merged across tolerance" + # ========================================================================= # --- case (d): chimera — disjoint query spans must NOT merge --- # This case LOCKS the query-overlap requirement of mb_reconcile_alt's merge diff --git a/test/altlg/test-revmulti.sh b/test/altlg/test-revmulti.sh new file mode 100644 index 0000000..00967b7 --- /dev/null +++ b/test/altlg/test-revmulti.sh @@ -0,0 +1,55 @@ +#!/bin/sh +# Regression test for reverse-multi-block and multi-record/hard-clip .alt lifting. +# Drives ex-lift-check in check mode: each " " triple asserts +# l2b_lift(ctg, pos) maps to primary , and that the contig's lift[] is +# sorted by alt_st (the invariant l2b_lift's binary search relies on). +# +# Expected primary coordinates (0-based), derived by hand from the .alt CIGARs: +# chrPrev_alt (reverse, 75M2000D75M @ chrP POS 201 -> 0-based 200): +# forward-ALT base 0 = RC base 149 -> primary 2349 +# forward-ALT base 74 = RC base 75 -> primary 2275 +# forward-ALT base 75 = RC base 74 -> primary 274 +# forward-ALT base 149 = RC base 0 -> primary 200 +# chrPhc_alt (50M @1001 + 50H50M @2001, both forward; POS 0-based 1000/2000): +# base 0 -> 1000 ; base 49 -> 1049 ; base 50 -> 2000 ; base 99 -> 2049 +# Usage: test/altlg/test-revmulti.sh [] +set -eu + +MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" +MKFIXTURE="$MDIR/test/altlg/mkfixture-revmulti.sh" +MINIBWA="$MDIR/minibwa" +EX_LIFT="$MDIR/api-test/ex-lift-check" + +TMPD=$(mktemp -d /tmp/altlg-test.XXXXXX) +trap 'rm -rf "$TMPD"' EXIT + +echo "[test-revmulti] building fixture in $TMPD ..." +/bin/sh "$MKFIXTURE" "$TMPD" + +echo "[test-revmulti] indexing with minibwa ..." +"$MINIBWA" index "$TMPD/ref.fa" 2>/dev/null + +echo "[test-revmulti] running ex-lift-check (check mode) ..." +# ex-lift-check exits non-zero on any MISMATCH or unsorted lift[]. +OUT=$("$EX_LIFT" "$TMPD/ref.fa" \ + chrPrev_alt 0 2349 \ + chrPrev_alt 74 2275 \ + chrPrev_alt 75 274 \ + chrPrev_alt 149 200 \ + chrPhc_alt 0 1000 \ + chrPhc_alt 49 1049 \ + chrPhc_alt 50 2000 \ + chrPhc_alt 99 2049) +echo "$OUT" | sed 's/^/ /' + +if echo "$OUT" | grep -q MISMATCH; then + echo "FAIL: at least one lift MISMATCH (see above)" + exit 1 +fi +if echo "$OUT" | grep -q 'sorted=0'; then + echo "FAIL: a contig's lift[] is not sorted by alt_st" + exit 1 +fi + +echo "[test-revmulti] PASS: reverse-multi-block and hard-clip multi-record .alt lift correctly" +exit 0 diff --git a/test/altlg/test-segdup-rev.sh b/test/altlg/test-segdup-rev.sh index 7f349c5..4866ac5 100755 --- a/test/altlg/test-segdup-rev.sh +++ b/test/altlg/test-segdup-rev.sh @@ -24,6 +24,7 @@ # # Usage: test/altlg/test-segdup-rev.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -35,8 +36,6 @@ N=20 TMPD=$(mktemp -d /tmp/altlg-segdup-rev.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # proj_line -> the MB_PROJ POS for chrP on that strand, # for the FULL-read seed (len=150); empty if none. diff --git a/test/altlg/test-segdup.sh b/test/altlg/test-segdup.sh index bf41276..9b68e33 100755 --- a/test/altlg/test-segdup.sh +++ b/test/altlg/test-segdup.sh @@ -15,6 +15,7 @@ # # Usage: test/altlg/test-segdup.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -26,13 +27,10 @@ N=20 # primary copies of CORE TMPD=$(mktemp -d /tmp/altlg-segdup.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # chrP_pos_list -> sorted space-separated list of chrP POS values chrP_pos_list() { mawk '$1!~/^@/ && $3=="chrP"{print $4}' "$1" | sort -n | tr '\n' ' '; } # has_pos -> "1" if any chrP record has that POS, else "0" -has_pos() { mawk -v p="$2" 'BEGIN{f=0} $1!~/^@/ && $3=="chrP" && $4==p{f=1} END{print f}' "$1"; } echo "[test-segdup] building segdup fixture (N=$N copies) ..." /bin/sh "$MK_SEG" "$TMPD" "$N" 2>/dev/null diff --git a/test/altlg/test-survive.sh b/test/altlg/test-survive.sh index 05090dc..a0c8746 100755 --- a/test/altlg/test-survive.sh +++ b/test/altlg/test-survive.sh @@ -11,16 +11,18 @@ # With the survival guard the twin is force-kept because its lifted placement # co-locates with the kept chrP primary (same pri_tid=chrP, same rev, |Δlst|=0). # -# RED assertion (guard absent): chrP_alt NOT in SAM output (-p 0.9 --outn=50). # GREEN assertion (guard present): chrP_alt IS in SAM output. -# Note: the RED state can be verified by temporarily passing NULL for l2b in -# the two mb_select_sub driver calls in mb_map_sai. +# RED assertion (guard ablated): with --dbg-no-alt-survive the guard is skipped +# and chrP_alt must be ABSENT. Asserting both states proves the fixture actually +# exercises the guard's drop path (not that chrP_alt survives for some unrelated +# reason), so the test cannot silently stop stressing the guard. # # Baseline (chrM, no .alt): guard never fires; output byte-identical without # guard because l2b has no ALT contigs (guard gated on l2b ALT presence). # # Usage: test/altlg/test-survive.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MKFIXTURE="$MDIR/test/altlg/mkfixture-survive.sh" @@ -32,8 +34,6 @@ MINIBWA="$MDIR/minibwa" TMPD=$(mktemp -d /tmp/altlg-survive.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # --------------------------------------------------------------------------- # Part 1: survival fixture (unique chrP + chrP_alt with 1bp SNP at read pos 55) @@ -83,6 +83,23 @@ if [ -n "$pri_pos" ] && [ -n "$alt_pos" ]; then ok "chrP_alt POS=$alt_pos co-locates with chrP POS=$pri_pos (|diff|=$d <= 10)" fi +# --- RED state: ablate ONLY the survival guard and confirm chrP_alt is dropped. +# This proves the fixture geometry genuinely drives the twin below the +# score-ratio/min-diff thresholds, so the GREEN assertion above is meaningful. --- +echo "[test-survive] mapping with --dbg-no-alt-survive (guard ablated) ..." +"$MINIBWA" map -p 0.9 --outn=50 --dbg-no-alt-survive "$TMPD/ref.fa" "$TMPD/reads.fq" 2>/dev/null \ + > "$TMPD/out.red.sam" +echo "----- SAM output (guard ablated) -----" +grep -v "^@" "$TMPD/out.red.sam" +echo "----- end -----" +red_pri=$(mawk '$1=="r-survive" && $3=="chrP" && ($2+0==0) {found=1} END{print found+0}' "$TMPD/out.red.sam") +[ "$red_pri" = "1" ] || fail "(RED) primary chrP alignment should still be present with the guard ablated" +red_alt=$(mawk '$1=="r-survive" && $3=="chrP_alt" {found=1} END{print found+0}' "$TMPD/out.red.sam") +if [ "$red_alt" != "0" ]; then + fail "(RED) chrP_alt should be DROPPED with the guard ablated, but it survived; the fixture is not exercising the guard's drop path" +fi +ok "chrP_alt dropped when guard ablated (RED): fixture provably stresses the survival guard" + # --------------------------------------------------------------------------- # Part 2: chrM baseline — guard never fires (no .alt -> l2b has no ALT ctgs) # --------------------------------------------------------------------------- diff --git a/test/altlg/test-svbreak.sh b/test/altlg/test-svbreak.sh index 7dab7bd..c9c2b54 100755 --- a/test/altlg/test-svbreak.sh +++ b/test/altlg/test-svbreak.sh @@ -26,6 +26,7 @@ # # Usage: test/altlg/test-svbreak.sh [] set -eu +. "$(dirname "$0")/lib.sh" MDIR="${1:-$(cd "$(dirname "$0")/../.." && pwd)}" MINIBWA="$MDIR/minibwa" @@ -37,18 +38,8 @@ MK_SV="$MDIR/test/altlg/mkfixture-svbreak.sh" TMPD=$(mktemp -d /tmp/altlg-svbreak.XXXXXX) trap 'rm -rf "$TMPD"' EXIT -fail() { echo "FAIL: $1"; exit 1; } -ok() { echo " ok: $1"; } # primary_flag/_mapq -> the SAM-primary record (no 0x100, no 0x800) -primary_flag() { - mawk -v q="$2" '$1==q { f=int($2); - if (int(f/256)%2==0 && int(f/2048)%2==0) { print $2; exit } }' "$1" -} -primary_mapq() { - mawk -v q="$2" '$1==q { f=int($2); - if (int(f/256)%2==0 && int(f/2048)%2==0) { print $5; exit } }' "$1" -} # secondary present on a given contig (0x100 set)? prints the FLAG of the first. alt_secondary_flag() { mawk -v q="$2" -v c="$3" '$1==q && $3==c { f=int($2);