sync: upstream f0e1174

.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

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=<path to a stock minibwa> 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; \

api-test/ex-lift-check.c

@@ -1,18 +1,46 @@
 #include <stdio.h>
 #include <string.h>
+#include <stdlib.h>
 #include "minibwa.h"
 #include "l2bit.h"
-/* Loads <idx>.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 <idx>.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 "<ctg> <pos> <expected>" triples, asserts each
+ *     lift's primary position equals <expected> (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 <idxprefix>\n");return 2;}
+    if (argc<2){fprintf(stderr,"usage: ex-lift-check <idxprefix> [<ctg> <pos> <expected> ...]\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+2<argc;a+=3){
+            const char *ctg=argv[a]; uint64_t pos=strtoull(argv[a+1],0,10); long exp=strtol(argv[a+2],0,10);
+            int64_t at=find_tid(l2b,ctg);
+            if(at<0||!l2b->ctg[at].is_alt){fprintf(stderr,"%s not flagged ALT\n",ctg);return 1;}
+            unsigned k; int sorted=1;
+            for(k=1;k<l2b->ctg[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);

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": []}

l2bit.c

@@ -1,6 +1,7 @@
 #define _POSIX_C_SOURCE 200809L
 #include <zlib.h>
 #include <stdio.h>
+#include <stdlib.h>
 #include <assert.h>
 #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;
 }
 

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;

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;