backup

align.c

@@ -1,6 +1,7 @@
 #include <assert.h>
 #include <string.h>
 #include "minimap.h"
+#include "mmpriv.h"
 #include "ksw2.h"
 
 #ifndef kroundup32
@@ -29,6 +30,21 @@ static inline void mm_seq_rev(uint32_t len, uint8_t *seq)
 		t = seq[i], seq[i] = seq[len - 1 - i], seq[len - 1 - i] = t;
 }
 
+static void mm_reg_split(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
+{
+	if (n <= 0 || n >= r->cnt) return;
+	*r2 = *r;
+	r2->p = 0;
+	r2->cnt = r->cnt - n;
+	r2->score = (int32_t)(r->score * ((float)r2->cnt / r->cnt) + .499);
+	r2->as = r->as + n;
+	mm_reg_set_coor(r2, qlen, a);
+	r->cnt -= r2->cnt;
+	r->score -= r2->score;
+	mm_reg_set_coor(r, qlen, a);
+	r->split = r2->split = 1;
+}
+
 static void mm_update_extra(mm_extra_t *p, const uint8_t *qseq, const uint8_t *tseq, uint32_t n_cigar, uint32_t *cigar)
 {
 	uint32_t k, l, toff = 0, qoff = 0;
@@ -119,6 +135,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	ksw_gen_simple_mat(5, mat, opt->a, opt->b);
 	bw = (int)(opt->bw * 1.5 + 1.);
 
+	r2->cnt = 0;
 	mm_adjust_minier(mi, qseq0, &a[r->as], &rs, &qs);
 	mm_adjust_minier(mi, qseq0, &a[r->as + r->cnt - 1], &re, &qe);
 
@@ -171,8 +188,15 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			mm_update_extra(r->p, qseq, tseq, ez->n_cigar, ez->cigar);
 			if (ez->score == KSW_NEG_INF) { // truncated by Z-drop
+				int j;
+				for (j = i - 1; j >= 0; --j)
+					if ((int32_t)a[r->as + j].x < re + ez->max_t)
+						break;
 				r->p->score += ez->max;
-				abort();
+				re1 = re + (ez->max_t + 1);
+				qe1 = qe + (ez->max_q + 1);
+				mm_reg_split(r, r2, j + 1, qlen, a);
+				break;
 			} else {
 				r->p->score += ez->score;
 			}

map.c

@@ -168,18 +168,8 @@ mm_reg1_t *mm_gen_reg(int qlen, int n_u, uint64_t *u, mm128_t *a)
 		ri->parent = -1, ri->subsc = 0;
 		ri->score = u[i]>>32;
 		ri->cnt = (int32_t)u[i];
-		ri->rev = a[k].x>>63;
-		ri->rid = a[k].x<<1>>33;
-		ri->rs = (int32_t)a[k].x + 1 > (int32_t)(a[k].y>>32)? (int32_t)a[k].x + 1 - (int32_t)(a[k].y>>32) : 0;
-		ri->re = (int32_t)a[k + ri->cnt - 1].x + 1;
-		if (!ri->rev) {
-			ri->qs = (int32_t)a[k].y + 1 - (int32_t)(a[k].y>>32);
-			ri->qe = (int32_t)a[k + ri->cnt - 1].y + 1;
-		} else {
-			ri->qs = qlen - ((int32_t)a[k + ri->cnt - 1].y + 1);
-			ri->qe = qlen - ((int32_t)a[k].y + 1 - (int32_t)(a[k].y>>32));
-		}
 		ri->as = k;
+		mm_reg_set_coor(ri, qlen, a);
 		k += ri->cnt;
 	}
 	return r;

minimap.h

@@ -51,7 +51,7 @@ typedef struct {
 } mm_extra_t;
 
 typedef struct {
-	uint32_t cnt:31, rev:1;
+	uint32_t cnt:30, rev:1, split:1;
 	uint32_t rid:31, rep:1;
 	int32_t score;
 	int32_t qs, qe, rs, re;

mmpriv.h

@@ -22,4 +22,20 @@ void mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 }
 #endif
 
+static inline void mm_reg_set_coor(mm_reg1_t *r, int32_t qlen, const mm128_t *a)
+{
+	int32_t k = r->as;
+	r->rev = a[k].x>>63;
+	r->rid = a[k].x<<1>>33;
+	r->rs = (int32_t)a[k].x + 1 > (int32_t)(a[k].y>>32)? (int32_t)a[k].x + 1 - (int32_t)(a[k].y>>32) : 0;
+	r->re = (int32_t)a[k + r->cnt - 1].x + 1;
+	if (!r->rev) {
+		r->qs = (int32_t)a[k].y + 1 - (int32_t)(a[k].y>>32);
+		r->qe = (int32_t)a[k + r->cnt - 1].y + 1;
+	} else {
+		r->qs = qlen - ((int32_t)a[k + r->cnt - 1].y + 1);
+		r->qe = qlen - ((int32_t)a[k].y + 1 - (int32_t)(a[k].y>>32));
+	}
+}
+
 #endif