r374: fix - clipping penalty not always working

This only happens to gaps where mem underestimates the bandwidth without
considering the clipping penalty.
This commit is contained in:
Heng Li 2013-04-10 01:09:37 -04:00
parent 53bb846407
commit 3d8a8c1e37
4 changed files with 6 additions and 6 deletions

View File

@ -568,7 +568,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int
for (i = 0; i < MAX_BAND_TRY; ++i) { for (i = 0; i < MAX_BAND_TRY; ++i) {
int prev = a->score; int prev = a->score;
aw[0] = opt->w << i; aw[0] = opt->w << i;
a->score = ksw_extend(s->qbeg, qs, tmp, rs, 5, opt->mat, opt->q, opt->r, aw[0], opt->zdrop, s->len * opt->a, &qle, &tle, &gtle, &gscore, &max_off[0]); a->score = ksw_extend(s->qbeg, qs, tmp, rs, 5, opt->mat, opt->q, opt->r, aw[0], opt->pen_clip, opt->zdrop, s->len * opt->a, &qle, &tle, &gtle, &gscore, &max_off[0]);
if (bwa_verbose >= 4) { printf("L\t%d < %d; w=%d; max_off=%d\n", prev, a->score, aw[0], max_off[0]); fflush(stdout); } if (bwa_verbose >= 4) { printf("L\t%d < %d; w=%d; max_off=%d\n", prev, a->score, aw[0], max_off[0]); fflush(stdout); }
if (a->score == prev || max_off[0] < (aw[0]>>1) + (aw[0]>>2)) break; if (a->score == prev || max_off[0] < (aw[0]>>1) + (aw[0]>>2)) break;
} }
@ -591,7 +591,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int
for (i = 0; i < MAX_BAND_TRY; ++i) { for (i = 0; i < MAX_BAND_TRY; ++i) {
int prev = a->score; int prev = a->score;
aw[1] = opt->w << i; aw[1] = opt->w << i;
a->score = ksw_extend(l_query - qe, query + qe, rmax[1] - rmax[0] - re, rseq + re, 5, opt->mat, opt->q, opt->r, aw[1], opt->zdrop, sc0, &qle, &tle, &gtle, &gscore, &max_off[1]); a->score = ksw_extend(l_query - qe, query + qe, rmax[1] - rmax[0] - re, rseq + re, 5, opt->mat, opt->q, opt->r, aw[1], opt->pen_clip, opt->zdrop, sc0, &qle, &tle, &gtle, &gscore, &max_off[1]);
if (bwa_verbose >= 4) { printf("R\t%d < %d; w=%d; max_off=%d\n", prev, a->score, aw[1], max_off[1]); fflush(stdout); } if (bwa_verbose >= 4) { printf("R\t%d < %d; w=%d; max_off=%d\n", prev, a->score, aw[1], max_off[1]); fflush(stdout); }
if (a->score == prev || max_off[1] < (aw[1]>>1) + (aw[1]>>2)) break; if (a->score == prev || max_off[1] < (aw[1]>>1) + (aw[1]>>2)) break;
} }

4
ksw.c
View File

@ -359,7 +359,7 @@ typedef struct {
int32_t h, e; int32_t h, e;
} eh_t; } eh_t;
int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int zdrop, int h0, int *_qle, int *_tle, int *_gtle, int *_gscore, int *_max_off) int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int end_bonus, int zdrop, int h0, int *_qle, int *_tle, int *_gtle, int *_gscore, int *_max_off)
{ {
eh_t *eh; // score array eh_t *eh; // score array
int8_t *qp; // query profile int8_t *qp; // query profile
@ -381,7 +381,7 @@ int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target,
k = m * m; k = m * m;
for (i = 0, max = 0; i < k; ++i) // get the max score for (i = 0, max = 0; i < k; ++i) // get the max score
max = max > mat[i]? max : mat[i]; max = max > mat[i]? max : mat[i];
max_gap = (int)((double)(qlen * max - gapo) / gape + 1.); max_gap = (int)((double)(qlen * max + end_bonus - gapo) / gape + 1.);
max_gap = max_gap > 1? max_gap : 1; max_gap = max_gap > 1? max_gap : 1;
w = w < max_gap? w : max_gap; w = w < max_gap? w : max_gap;
// DP loop // DP loop

2
ksw.h
View File

@ -102,7 +102,7 @@ extern "C" {
* *
* @return best semi-local alignment score * @return best semi-local alignment score
*/ */
int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int zdrop, int h0, int *qle, int *tle, int *gtle, int *gscore, int *max_off); int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int end_bonus, int zdrop, int h0, int *qle, int *tle, int *gtle, int *gscore, int *max_off);
#ifdef __cplusplus #ifdef __cplusplus
} }

2
main.c
View File

@ -3,7 +3,7 @@
#include "utils.h" #include "utils.h"
#ifndef PACKAGE_VERSION #ifndef PACKAGE_VERSION
#define PACKAGE_VERSION "0.7.3-r373-beta" #define PACKAGE_VERSION "0.7.3-r374-beta"
#endif #endif
int bwa_fa2pac(int argc, char *argv[]); int bwa_fa2pac(int argc, char *argv[]);