r438: fixed a rare case that leads to missing hits

It is a bug in chaining.
This commit is contained in:
Heng Li 2017-09-25 14:59:34 -04:00
parent 8f25cfa36e
commit 9538c985aa
2 changed files with 5 additions and 4 deletions

View File

@ -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;
}

4
main.c
View File

@ -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 <sys/resource.h>
@ -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