diff --git a/align.c b/align.c index 7f0fab1..05af553 100644 --- a/align.c +++ b/align.c @@ -1,6 +1,7 @@ #include #include #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; } diff --git a/map.c b/map.c index 43e360a..b01d4e7 100644 --- a/map.c +++ b/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; diff --git a/minimap.h b/minimap.h index dbb0017..e9dafb1 100644 --- a/minimap.h +++ b/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; diff --git a/mmpriv.h b/mmpriv.h index 66d8187..1314b99 100644 --- a/mmpriv.h +++ b/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