diff --git a/hit.c b/hit.c index 33f9f3f..9dfa7a4 100644 --- a/hit.c +++ b/hit.c @@ -408,7 +408,33 @@ void mm_seg_free(void *km, int n_segs, mm_seg_t *segs) kfree(km, segs); } -void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr) +static void mm_set_inv_mapq(void *km, int n_regs, mm_reg1_t *regs) +{ + int i, n_aux; + uint64_t *aux; + if (n_regs < 3) return; + for (i = 0; i < n_regs; ++i) + if (regs[i].inv) break; + if (i == n_regs) return; // no inversion hits + + aux = (uint64_t*)kmalloc(km, n_regs * 8); + for (i = n_aux = 0; i < n_regs; ++i) + if (regs[i].parent == i || regs[i].parent < 0) + aux[n_aux++] = (uint64_t)regs[i].as << 32 | i; + radix_sort_64(aux, aux + n_aux); + + for (i = 1; i < n_aux - 1; ++i) { + mm_reg1_t *inv = ®s[(int32_t)aux[i]]; + if (inv->inv) { + mm_reg1_t *l = ®s[(int32_t)aux[i-1]]; + mm_reg1_t *r = ®s[(int32_t)aux[i+1]]; + inv->mapq = l->mapq < r->mapq? l->mapq : r->mapq; + } + } + kfree(km, aux); +} + +void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr) { static const float q_coef = 40.0f; int64_t sum_sc = 0; @@ -451,4 +477,5 @@ void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, in if (r->p && r->p->dp_max > r->p->dp_max2 && r->mapq == 0) r->mapq = 1; } else r->mapq = 0; } + mm_set_inv_mapq(km, n_regs, regs); } diff --git a/main.c b/main.c index 6a0f333..4d1e45a 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "getopt.h" -#define MM_VERSION "2.8-r710-dirty" +#define MM_VERSION "2.8-r711-dirty" #ifdef __linux__ #include diff --git a/map.c b/map.c index ad55cf8..4dcff71 100644 --- a/map.c +++ b/map.c @@ -340,7 +340,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** if (n_segs == 1) { // uni-segment regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a); - mm_set_mapq(n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr); + mm_set_mapq(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr); n_regs[0] = n_regs0, regs[0] = regs0; } else { // multi-segment mm_seg_t *seg; @@ -349,7 +349,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** for (i = 0; i < n_segs; ++i) { mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b); // update mm_reg1_t::parent regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a); - mm_set_mapq(n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr); + mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr); } mm_seg_free(b->km, n_segs, seg); if (n_segs == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR)) diff --git a/mmpriv.h b/mmpriv.h index 36c2033..3279695 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -78,7 +78,7 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int void mm_filter_regs(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs); void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs, mm128_t *a); void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r); -void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr); +void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr); void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);