ALT-aware mapping via post-extension liftover groups

.gitignore

@@ -4,3 +4,9 @@
 Makefile.bak
 minibwa
 tmp
+
+# api-test helper programs. `make test` builds these from tracked sources in
+# api-test/, so the binaries are artifacts.
+api-test/ex-*-check
+api-test/mbmap-batch
+api-test/mbmap-one

Makefile

@@ -39,7 +39,7 @@ ifeq ($(ARCH), x86_64)
 endif
 
 .SUFFIXES:.c .o
-.PHONY:all clean depend
+.PHONY:all clean depend test
 
 .c.o:
 		$(CC) -c $(CFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
@@ -55,6 +55,23 @@ libminibwa.a:$(LOBJS)
 minibwa:libminibwa.a $(MALLOC_O) $(AOBJS) main.o
 		$(CC) $(CFLAGS) $(LDFLAGS) $(MALLOC_O) $(AOBJS) main.o -o $@ -L. -lminibwa $(LIBS)
 
+# Run the ALT liftover-group suite. Each script builds its own fixtures, but several
+# drive the api-test helpers, which live in their own Makefile and are not part of `all`
+# -- without them those scripts fail on a fresh checkout. Set MB_STOCK to a stock
+# 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; \
+			else echo FAIL; cat /tmp/altlg-log.$$$$; fail=1; fi; \
+			rm -f /tmp/altlg-log.$$$$; \
+		done; exit $$fail
+
 clean:
 		rm -fr *.o a.out $(PROG) *~ *.a *.dSYM
 

README.md

@@ -44,8 +44,10 @@ structural changes anyway.
 
 Minibwa is designed for mapping short reads and accurate long reads. It does
 not support spliced alignment and has not been tuned for aligning long contigs.
-For now, minibwa does not properly work with alternate contigs in the reference
-genome. Please use a version of the reference without such contigs.
+Alternate contigs are supported: when a `<idx>.alt` file sits beside the index
+it is loaded automatically and ALT hits are reconciled with their primary
+counterparts, so an ALT twin no longer collapses the primary's mapping quality.
+Pass `--no-alt` to ignore it and align exactly as if no ALT file existed.
 
 ### Installation
 
@@ -130,7 +132,6 @@ include GPL'd source code.
 
 * Minibwa does not work with noisy long reads or spliced RNA-seq reads.
 * Minibwa does not support undirectional bisulfite sequencing data.
-* Minibwa does not recognize alternate haplotypes.
 
 [apache2]: https://github.com/lh3/minibwa/tree/Apache2
 [zlib]: https://zlib.net/

api-test/Makefile

@@ -1,6 +1,6 @@
 CC=			gcc
 CFLAGS=		-std=gnu99 -Wall -O3
-EXE=		mbmap-one mbmap-batch
+EXE=		mbmap-one mbmap-batch ex-lift-check ex-place-check ex-group-check
 
 .PHONY:all clean
 
@@ -10,10 +10,19 @@ all:$(EXE)
 	(cd ..; make libminibwa.a)
 
 mbmap-one:ex-one.c ../libminibwa.a
-	$(CC) $(CFLAGS) -o $@ -I.. $^ -lz -lm
+	$(CC) $(CFLAGS) -o $@ -I.. $^ -lz -lm $(LIBS)
 
 mbmap-batch:ex-batch.c ../libminibwa.a
-	$(CC) $(CFLAGS) -o $@ -I.. $^ -lz -lm
+	$(CC) $(CFLAGS) -o $@ -I.. $^ -lz -lm $(LIBS)
+
+ex-lift-check:ex-lift-check.c ../libminibwa.a
+	$(CC) $(CFLAGS) -o $@ -I.. $^ -lz -lm $(LIBS)
+
+ex-place-check:ex-place-check.c ../libminibwa.a
+	$(CC) $(CFLAGS) -o $@ -I.. $^ -lz -lm $(LIBS)
+
+ex-group-check:ex-group-check.c ../libminibwa.a
+	$(CC) $(CFLAGS) -o $@ -I.. $^ -lz -lm $(LIBS)
 
 clean:
 	rm -f $(EXE)

api-test/ex-group-check.c

@@ -0,0 +1,59 @@
+#include <stdlib.h>
+#include <assert.h>
+#include <stdio.h>
+#include <zlib.h>
+#include "minibwa.h"
+#include "kseq.h"
+KSEQ_INIT(gzFile, gzread)
+
+/* Probe for Task 4 (mb_reconcile_alt grouping).  mb_map() runs the full SE
+ * pipeline, including the gated mb_reconcile_alt pass, so the per-hit `parent`
+ * field already reflects the FINAL liftover grouping: a hit is its own group
+ * representative iff parent == id, otherwise it is a subordinate whose parent is
+ * the id of its group's representative.  For every hit of every read it prints:
+ *
+ *   <qname> ctg=<ctg> qs=<qs> qe=<qe> rev=<rev> is_alt=<is_alt> \
+ *           id=<id> parent=<parent> subsc=<subsc> dp_max2=<dp_max2|.>
+ *
+ * so the test can assert directly on the grouping (parent == id => independent
+ * representative; parent == <other id> => merged subordinate) without depending
+ * on SAM-flag heuristics that are insensitive to the merge decision. */
+int main(int argc, char *argv[])
+{
+	mb_opt_t opt;
+	mb_opt_init(&opt);
+
+	if (argc < 3) {
+		fprintf(stderr, "Usage: ex-group-check <idxPrefix> <query.fa>\n");
+		return 1;
+	}
+
+	gzFile f = gzopen(argv[2], "r");
+	assert(f);
+	kseq_t *ks = kseq_init(f);
+
+	mb_idx_t *idx = mb_idx_load(argv[1], 0);
+	assert(idx);
+	/* The loaders do not resolve .alt; the caller does.  This example needs it. */
+	mb_idx_set_alt_auto(idx, argv[1]);
+
+	while (kseq_read(ks) >= 0) {
+		mb_hit_t *hit;
+		int32_t j, n_hit;
+		hit = mb_map(&opt, idx, ks->seq.l, ks->seq.s, 0, &n_hit, 0, ks->name.s);
+		for (j = 0; j < n_hit; ++j) {
+			mb_hit_t *h = &hit[j];
+			printf("%s\tctg=%s\tqs=%d\tqe=%d\trev=%d\tis_alt=%d\tid=%d\tparent=%d\tsubsc=%d\tdp_max2=",
+				ks->name.s, mb_idx_ctg_name(idx, h->tid),
+				h->qs, h->qe, h->rev, h->is_alt, h->id, h->parent, h->subsc);
+			if (h->p) printf("%d\n", h->p->dp_max2);
+			else      printf(".\n");
+			free(h->p);
+		}
+		free(hit);
+	}
+	mb_idx_destroy(idx);
+	kseq_destroy(ks);
+	gzclose(f);
+	return 0;
+}

api-test/ex-lift-check.c

@@ -0,0 +1,49 @@
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include "minibwa.h"
+#include "l2bit.h"
+/* 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> [<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);
+    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;}
+    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);
+    printf("lift(100)=%d pri=%llu  lift(400)=%d pri=%llu\n",ok1,(unsigned long long)pp1,ok2,(unsigned long long)pp2);
+    return 0;
+}

api-test/ex-place-check.c

@@ -0,0 +1,80 @@
+#include <stdlib.h>
+#include <assert.h>
+#include <stdio.h>
+#include <string.h>
+#include <zlib.h>
+#include "minibwa.h"
+#include "mbpriv.h"   /* mb_idx_s -> l2b for mb_hit_place */
+#include "kseq.h"
+KSEQ_INIT(gzFile, gzread)
+
+/* Probe for Task 2 (mb_hit_place).  For every hit of every read it prints:
+ *   <qname> <ctg> ts=<ts> te=<te> rev=<rev> is_alt=<is_alt> cg=<cigar|.> ||
+ *     place: pri=<pri_ctg> lst=<lifted_st> len=<lifted_en> prev=<rev> liftable=<liftable>
+ * so the test script can assert on the placement (pri contig name, lifted_st,
+ * folded rev, liftable) without depending on internal tid numbering.
+ *
+ * Optional 3rd arg <hitspec> = "all" (default) or "coarse": when "coarse" the
+ * probe NULLs h->p before calling mb_hit_place to exercise the pre-DP [ts,te]
+ * path against the very same hit (so coarse vs exact can be compared directly). */
+int main(int argc, char *argv[])
+{
+	mb_opt_t opt;
+	int coarse = 0;
+	mb_opt_init(&opt);
+
+	if (argc < 3) {
+		fprintf(stderr, "Usage: ex-place-check <idxPrefix> <query.fa> [all|coarse]\n");
+		return 1;
+	}
+	if (argc >= 4 && strcmp(argv[3], "coarse") == 0) coarse = 1;
+
+	gzFile f = gzopen(argv[2], "r");
+	assert(f);
+	kseq_t *ks = kseq_init(f);
+
+	mb_idx_t *idx = mb_idx_load(argv[1], 0);
+	assert(idx);
+	/* The loaders do not resolve .alt; the caller does.  This example needs it. */
+	mb_idx_set_alt_auto(idx, argv[1]);
+
+	while (kseq_read(ks) >= 0) {
+		mb_hit_t *hit;
+		int32_t i, j, n_hit;
+		hit = mb_map(&opt, idx, ks->seq.l, ks->seq.s, 0, &n_hit, 0, ks->name.s);
+		for (j = 0; j < n_hit; ++j) {
+			mb_hit_t *h = &hit[j];
+			mb_extra_t *saved = h->p;
+			mb_place_t pl;
+			const char *pri_name;
+
+			printf("%s\t%s\tts=%ld\tte=%ld\trev=%d\tis_alt=%d\tcg=",
+				ks->name.s, mb_idx_ctg_name(idx, h->tid),
+				(long)h->ts, (long)h->te, h->rev, h->is_alt);
+			if (h->p && h->p->n_cigar > 0) {
+				for (i = 0; i < h->p->n_cigar; ++i)
+					printf("%d%c", h->p->cigar[i]>>4, MB_CIGAR_STR[h->p->cigar[i]&0xf]);
+			} else putchar('.');
+
+			if (coarse) h->p = 0;           /* force the pre-DP [ts,te] path */
+			pl = mb_hit_place(idx->l2b, h);
+			/* Restore before printing so ctg name lookup via idx still works if
+			 * needed.  Free `saved` directly -- it is the one allocation we own;
+			 * h->p after the restore is the same pointer, but freeing via `saved`
+			 * makes ownership unambiguous and avoids any confusion with the
+			 * subsequent free(hit) that releases only the flat hit array. */
+			h->p = saved;
+			free(saved);
+			saved = h->p = NULL;            /* prevent any accidental re-use */
+
+			pri_name = (pl.pri_tid >= 0) ? mb_idx_ctg_name(idx, pl.pri_tid) : "?";
+			printf("\t||\tplace:\tpri=%s\tlst=%ld\tlen=%ld\tprev=%d\tliftable=%d\n",
+				pri_name, (long)pl.lifted_st, (long)pl.lifted_en, pl.rev, pl.liftable);
+		}
+		free(hit);
+	}
+	mb_idx_destroy(idx);
+	kseq_destroy(ks);
+	gzclose(f);
+	return 0;
+}

format.c

@@ -274,7 +274,8 @@ void mb_fmt_sam(void *km, kstring_t *s, const l2b_t *l2b, const mb_bseq1_t *t, i
 		if (t->qual) sam_write_sq(s, t->qual, t->l_seq, 0, 0);
 		else kom_sprintf_lite(s, "*");
 	} else {
-		if ((flag & 0x900) == 0 || (opt->flag & MB_F_SUPP_SOFT)) {
+		if ((flag & 0x900) == 0 || (opt->flag & MB_F_SUPP_SOFT)
+		    || ((opt->flag & MB_F_ALT_RECORDS) && r->is_alt)) {
 			sam_write_sq(s, t->seq, t->l_seq, r->rev, r->rev);
 			kom_sprintf_lite(s, "\t");
 			if (t->qual) sam_write_sq(s, t->qual, t->l_seq, r->rev, 0);