From 9538c985aa6342c82d8781e3a912373c6938c4fd Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 25 Sep 2017 14:59:34 -0400 Subject: [PATCH] r438: fixed a rare case that leads to missing hits It is a bug in chaining. --- chain.c | 5 +++-- main.c | 4 ++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/chain.c b/chain.c index 9c3ffeb..6c86bcb 100644 --- a/chain.c +++ b/chain.c @@ -73,7 +73,8 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m } if (p[j] >= 0) t[p[j]] = i; } - f[i] = max_f, p[i] = max_j, v[i] = max_f_past; // v[] keeps the max score in the previous chain + f[i] = max_f, p[i] = max_j; + v[i] = max_f_past > max_f? max_f_past : max_f; // v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak } // find the ending positions of chains @@ -91,7 +92,7 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m for (i = n_u = 0; i < n; ++i) { if (t[i] == 0 && v[i] >= min_sc) { j = i; - while (j >= 0 && f[j] < v[j]) j = p[j]; // find the point that maximizes f[] + while (j >= 0 && f[j] < v[j]) j = p[j]; // find the peak that maximizes f[] if (j < 0) j = i; // TODO: this should really be assert(j>=0) u[n_u++] = (uint64_t)f[j] << 32 | j; } diff --git a/main.c b/main.c index 16af19b..b89917c 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "getopt.h" -#define MM_VERSION "2.2-r437-dirty" +#define MM_VERSION "2.2-r438-dirty" #ifdef __linux__ #include @@ -109,7 +109,7 @@ int main(int argc, char *argv[]) else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc else if (c == 0 && long_idx == 4) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname else if (c == 0 && long_idx == 5) opt.flag |= MM_F_NO_SELF; // --no-self - else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED; // --print-seed + else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ; // --print-aln-seq