r428: allow to change the default chain_drop_ratio

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;
-		int w = 0, tmp;
+		w = mem_chain_weight(c);
-		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;
 		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;

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);

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[]);