From 10cb6b05077ad193dcf00adf814b780dbf2ed564 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 30 Dec 2013 16:18:45 -0500 Subject: [PATCH] r428: allow to change the default chain_drop_ratio --- bwamem.c | 40 +++++++++++++++++++++++----------------- fastmap.c | 4 +++- main.c | 2 +- 3 files changed, 27 insertions(+), 19 deletions(-) diff --git a/bwamem.c b/bwamem.c index a609e0d..9b05528 100644 --- a/bwamem.c +++ b/bwamem.c @@ -232,12 +232,32 @@ static void mem_insert_seed(const mem_opt_t *opt, int64_t l_pac, kbtree_t(chn) * } } +int mem_chain_weight(const mem_chain_t *c) +{ + int64_t end; + int j, w = 0, tmp; + for (j = 0, end = 0; j < c->n; ++j) { + const mem_seed_t *s = &c->seeds[j]; + if (s->qbeg >= end) w += s->len; + else if (s->qbeg + s->len > end) w += s->qbeg + s->len - end; + end = end > s->qbeg + s->len? end : s->qbeg + s->len; + } + tmp = w; + for (j = 0, end = 0; j < c->n; ++j) { + const mem_seed_t *s = &c->seeds[j]; + if (s->rbeg >= end) w += s->len; + else if (s->rbeg + s->len > end) w += s->rbeg + s->len - end; + end = end > s->qbeg + s->len? end : s->qbeg + s->len; + } + return w < tmp? w : tmp; +} + void mem_print_chain(const bntseq_t *bns, mem_chain_v *chn) { int i, j; for (i = 0; i < chn->n; ++i) { mem_chain_t *p = &chn->a[i]; - err_printf("CHAIN(%d) n=%d", i, p->n); + err_printf("CHAIN(%d) n=%d w=%d", i, p->n, mem_chain_weight(p)); for (j = 0; j < p->n; ++j) { bwtint_t pos; int is_rev, ref_id; @@ -294,22 +314,8 @@ int mem_chain_flt(const mem_opt_t *opt, int n_chn, mem_chain_t *chains) a = malloc(sizeof(flt_aux_t) * n_chn); for (i = 0; i < n_chn; ++i) { mem_chain_t *c = &chains[i]; - int64_t end; - int w = 0, tmp; - for (j = 0, end = 0; j < c->n; ++j) { - const mem_seed_t *s = &c->seeds[j]; - if (s->qbeg >= end) w += s->len; - else if (s->qbeg + s->len > end) w += s->qbeg + s->len - end; - end = end > s->qbeg + s->len? end : s->qbeg + s->len; - } - tmp = w; - for (j = 0, end = 0; j < c->n; ++j) { - const mem_seed_t *s = &c->seeds[j]; - if (s->rbeg >= end) w += s->len; - else if (s->rbeg + s->len > end) w += s->rbeg + s->len - end; - end = end > s->qbeg + s->len? end : s->qbeg + s->len; - } - w = w < tmp? w : tmp; + int w; + w = mem_chain_weight(c); a[i].beg = c->seeds[0].qbeg; a[i].end = c->seeds[c->n-1].qbeg + c->seeds[c->n-1].len; a[i].w = w; a[i].p = c; a[i].p2 = 0; diff --git a/fastmap.c b/fastmap.c index 1884ebc..40cea8c 100644 --- a/fastmap.c +++ b/fastmap.c @@ -30,7 +30,7 @@ int main_mem(int argc, char *argv[]) int64_t n_processed = 0; opt = mem_opt_init(); - while ((c = getopt(argc, argv, "paMCSPHk:c:v:s:r:t:R:A:B:O:E:U:w:L:d:T:Q:")) >= 0) { + while ((c = getopt(argc, argv, "paMCSPHk:c:v:s:r:t:R:A:B:O:E:U:w:L:d:T:Q:D:")) >= 0) { if (c == 'k') opt->min_seed_len = atoi(optarg); else if (c == 'w') opt->w = atoi(optarg); else if (c == 'A') opt->a = atoi(optarg); @@ -49,6 +49,7 @@ int main_mem(int argc, char *argv[]) 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 == 'D') opt->chain_drop_ratio = atof(optarg); else if (c == 'C') copy_comment = 1; else if (c == 'Q') { opt->mapQ_coef_len = atoi(optarg); @@ -75,6 +76,7 @@ int main_mem(int argc, char *argv[]) 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); + fprintf(stderr, " -D FLOAT drop chains shorter than FLOAT fraction of the longest overlapping chain [%.2f]\n", opt->chain_drop_ratio); fprintf(stderr, " -S skip mate rescue\n"); fprintf(stderr, " -P skip pairing; mate rescue performed unless -S also in use\n"); fprintf(stderr, " -A INT score for a sequence match [%d]\n", opt->a); diff --git a/main.c b/main.c index 3595cbc..bb11578 100644 --- a/main.c +++ b/main.c @@ -4,7 +4,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.7.5a-r427" +#define PACKAGE_VERSION "0.7.5a-r428" #endif int bwa_fa2pac(int argc, char *argv[]);