Permalink
Browse files

r323: added Z-dropoff, a variant of blast's X-drop

  • Loading branch information...
1 parent d6096c3 commit e0991d6a459daa27031df17405e0d62f0395a6f1 @lh3 committed Mar 5, 2013
Showing with 13 additions and 7 deletions.
  1. +3 −2 bwamem.c
  2. +1 −0 bwamem.h
  3. +5 −1 fastmap.c
  4. +2 −2 ksw.c
  5. +1 −1 ksw.h
  6. +1 −1 main.c
View
@@ -44,6 +44,7 @@ mem_opt_t *mem_opt_init()
o = calloc(1, sizeof(mem_opt_t));
o->flag = 0;
o->a = 1; o->b = 4; o->q = 6; o->r = 1; o->w = 60; o->max_w = 500;
+ o->zdrop = 100;
o->pen_unpaired = 9;
o->pen_clip = 5;
o->min_seed_len = 19;
@@ -560,7 +561,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int
rs = malloc(tmp);
for (i = 0; i < tmp; ++i) rs[i] = rseq[tmp - 1 - i];
for (aw[0] = opt->w;; aw[0] <<= 1) {
- a->score = ksw_extend(s->qbeg, qs, tmp, rs, 5, opt->mat, opt->q, opt->r, aw[0], 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->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", tmps, a->score, aw[0], max_off[0]); fflush(stdout);
if (a->score == tmps || aw[0]<<1 > opt->max_w || max_off[0] < (aw[0]>>1) + (aw[0]>>2)) break;
tmps = a->score;
@@ -577,7 +578,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int
re = s->rbeg + s->len - rmax[0];
assert(re >= 0);
for (aw[1] = opt->w;; aw[1] <<= 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], 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->zdrop, sc0, &qle, &tle, &gtle, &gscore, &max_off[1]);
if (bwa_verbose >= 4) printf("R\t%d < %d; w=%d; max_off=%d\n", tmps, a->score, aw[1], max_off[1]); fflush(stdout);
if (a->score == tmps || aw[1]<<1 > opt->max_w || max_off[1] < (aw[1]>>1) + (aw[1]>>2)) break;
tmps = a->score;
View
@@ -23,6 +23,7 @@ typedef struct {
int pen_clip; // clipping penalty. This score is not deducted from the DP score.
int w; // band width
int max_w; // max band width
+ int zdrop; // Z-dropoff
int flag; // see MEM_F_* macros
int min_seed_len; // minimum seed length
View
@@ -26,9 +26,10 @@ int main_mem(int argc, char *argv[])
void *ko = 0, *ko2 = 0;
opt = mem_opt_init();
- while ((c = getopt(argc, argv, "paMCPHk:c:v:s:r:t:R:A:B:O:E:U:w:L:")) >= 0) {
+ while ((c = getopt(argc, argv, "paMCPHk:c:v:s:r:t:R:A:B:O:E:U:w:L:d:W:")) >= 0) {
if (c == 'k') opt->min_seed_len = atoi(optarg);
else if (c == 'w') opt->w = atoi(optarg);
+ else if (c == 'W') opt->max_w = atoi(optarg);
else if (c == 'A') opt->a = atoi(optarg);
else if (c == 'B') opt->b = atoi(optarg);
else if (c == 'O') opt->q = atoi(optarg);
@@ -42,6 +43,7 @@ int main_mem(int argc, char *argv[])
else if (c == 'p') opt->flag |= MEM_F_PE;
else if (c == 'M') opt->flag |= MEM_F_NO_MULTI;
else if (c == 'c') opt->max_occ = atoi(optarg);
+ else if (c == 'd') opt->zdrop = atoi(optarg);
else if (c == 'v') bwa_verbose = atoi(optarg);
else if (c == 'r') opt->split_factor = atof(optarg);
else if (c == 'C') copy_comment = 1;
@@ -57,6 +59,8 @@ int main_mem(int argc, char *argv[])
fprintf(stderr, " -t INT number of threads [%d]\n", opt->n_threads);
fprintf(stderr, " -k INT minimum seed length [%d]\n", opt->min_seed_len);
fprintf(stderr, " -w INT band width for banded alignment [%d]\n", opt->w);
+ fprintf(stderr, " -W INT max band width [%d]\n", opt->max_w);
+ fprintf(stderr, " -d INT off-diagnal X-dropoff [%d]\n", opt->zdrop);
fprintf(stderr, " -r FLOAT look for internal seeds inside a seed longer than {-k} * FLOAT [%g]\n", opt->split_factor);
// fprintf(stderr, " -s INT look for internal seeds inside a seed with less than INT occ [%d]\n", opt->split_width);
fprintf(stderr, " -c INT skip seeds with more than INT occurrences [%d]\n", opt->max_occ);
View
4 ksw.c
@@ -359,7 +359,7 @@ typedef struct {
int32_t h, e;
} 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 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 zdrop, int h0, int *_qle, int *_tle, int *_gtle, int *_gscore, int *_max_off)
{
eh_t *eh; // score array
int8_t *qp; // query profile
@@ -426,7 +426,7 @@ int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target,
max_ie = gscore > h1? max_ie : i;
gscore = gscore > h1? gscore : h1;
}
- if (m == 0) break;
+ if (m == 0 || max - m - abs((i - max_i) - (j - max_j)) * gape > zdrop) break; // drop to zero, or below Z-dropoff
if (m > max) {
max = m, max_i = i, max_j = mj;
max_off = max_off > abs(mj - i)? max_off : abs(mj - i);
View
2 ksw.h
@@ -102,7 +102,7 @@ extern "C" {
*
* @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 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 zdrop, int h0, int *qle, int *tle, int *gtle, int *gscore, int *max_off);
#ifdef __cplusplus
}
View
2 main.c
@@ -4,7 +4,7 @@
#include "utils.h"
#ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.7.0-r320"
+#define PACKAGE_VERSION "0.7.0-r323-beta"
#endif
static int usage()

0 comments on commit e0991d6

Please sign in to comment.