From 8a1d52bcbed9b33702993499cb44469f41c07fc7 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 6 Aug 2021 21:40:43 -0400 Subject: [PATCH] r1094: for --split-prefix update max_dp at the end --- align.c | 4 ++-- main.c | 2 +- map.c | 8 ++++++++ mmpriv.h | 1 + 4 files changed, 12 insertions(+), 3 deletions(-) diff --git a/align.c b/align.c index e587bbc..ddbb0bd 100644 --- a/align.c +++ b/align.c @@ -931,7 +931,7 @@ static int32_t mm_recal_max_dp(const mm_reg1_t *r, double b2, int32_t match_sc) return (int32_t)(match_sc * (r->mlen - b2 * n_mis - gap_cost) + .499); } -static void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b) +void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b) { int32_t max = -1, max2 = -1, i, max_i = -1; double div, b2; @@ -1011,7 +1011,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m kfree(km, qseq0[0]); kfree(km, ez.cigar); mm_filter_regs(opt, qlen, n_regs_, regs); - if (!(opt->flag&MM_F_SR) && qlen >= opt->rank_min_len) { + if (!(opt->flag&MM_F_SR) && !opt->split_prefix && qlen >= opt->rank_min_len) { mm_update_dp_max(qlen, *n_regs_, regs, opt->rank_frac, opt->a, opt->b); mm_filter_regs(opt, qlen, n_regs_, regs); } diff --git a/main.c b/main.c index 81a6980..9869015 100644 --- a/main.c +++ b/main.c @@ -7,7 +7,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.21-dev-r1090-dirty" +#define MM_VERSION "2.21-dev-r1094-dirty" #ifdef __linux__ #include diff --git a/map.c b/map.c index f03b7ea..9c4ad64 100644 --- a/map.c +++ b/map.c @@ -494,6 +494,14 @@ static void merge_hits(step_t *s) } } } + if (!(opt->flag&MM_F_SR) && s->seq[k].l_seq >= opt->rank_min_len) + mm_update_dp_max(s->seq[k].l_seq, s->n_reg[k], s->reg[k], opt->rank_frac, opt->a, opt->b); + for (j = 0; j < s->n_reg[k]; ++j) { + mm_reg1_t *r = &s->reg[k][j]; + if (r->p) r->p->dp_max2 = 0; // reset ->dp_max2 as mm_set_parent() doesn't clear it; necessary with mm_update_dp_max() + r->subsc = 0; // this may not be necessary + r->n_sub = 0; // n_sub will be an underestimate as we don't see all the chains now, but it can't be accurate anyway + } mm_hit_sort(km, &s->n_reg[k], s->reg[k], opt->alt_drop); mm_set_parent(km, opt->mask_level, opt->mask_len, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop); if (!(opt->flag & MM_F_ALL_CHAINS)) { diff --git a/mmpriv.h b/mmpriv.h index db9143d..933e41f 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -95,6 +95,7 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs); void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r, float alt_diff_frac); 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_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b); 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);