From 43b498a37e89acde2efbaabcff436a3060d7d6d9 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@me.com>
Date: Fri, 9 May 2014 14:56:59 -0400
Subject: [PATCH] r759: bugfix - frac_rep not working

Also added commented code for a 3rd round seeding. Not used.
---
 bwamem.c      | 22 +++++++++++++++++++++-
 bwamem.h      |  2 +-
 bwamem_pair.c |  1 +
 main.c        |  2 +-
 4 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/bwamem.c b/bwamem.c
index 5f51030..28d8c88 100644
--- a/bwamem.c
+++ b/bwamem.c
@@ -131,7 +131,27 @@ static void mem_collect_intv(const mem_opt_t *opt, const bwt_t *bwt, int len, co
 		if (end - start < split_len || p->x[2] > opt->split_width) continue;
 		bwt_smem1(bwt, len, seq, (start + end)>>1, p->x[2]+1, &a->mem1, a->tmpv);
 		for (i = 0; i < a->mem1.n; ++i)
-			kv_push(bwtintv_t, a->mem, a->mem1.a[i]);
+			if ((a->mem1.a[i].info>>32) - (uint32_t)a->mem1.a[i].info >= opt->min_seed_len)
+				kv_push(bwtintv_t, a->mem, a->mem1.a[i]);
+/*
+		{ // third pass:
+			int max_len = 0, max_i = -1;
+			for (i = 0; i < a->mem1.n; ++i) {
+				bwtintv_t *q = &a->mem1.a[i];
+				int start = q->info>>32, end = (int32_t)q->info;
+				if (end - start > max_len)
+					max_len = end - start, max_i = i;
+			}
+			if (max_i >= 0 && max_len > .9 * (end - start) && max_len >= split_len) {
+				bwtintv_t t = a->mem1.a[max_i];
+				int start = t.info>>32, end = (int32_t)t.info;
+				bwt_smem1(bwt, len, seq, (start + end)>>1, t.x[2]+1, &a->mem1, a->tmpv);
+				for (i = 0; i < a->mem1.n; ++i)
+					if ((a->mem1.a[i].info>>32) - (uint32_t)a->mem1.a[i].info >= opt->min_seed_len)
+						kv_push(bwtintv_t, a->mem, a->mem1.a[i]);
+			}
+		}
+*/
 	}
 	// sort
 	ks_introsort(mem_intv, a->mem.n, a->mem.a);
diff --git a/bwamem.h b/bwamem.h
index d4e2d1d..98e135d 100644
--- a/bwamem.h
+++ b/bwamem.h
@@ -65,7 +65,7 @@ typedef struct {
 	int secondary;  // index of the parent hit shadowing the current hit; <0 if primary
 	int seedlen0;   // length of the starting seed
 	int n_comp;     // number of sub-alignments chained together
-	int frac_rep;
+	float frac_rep;
 	uint64_t hash;
 } mem_alnreg_t;
 
diff --git a/bwamem_pair.c b/bwamem_pair.c
index 202e050..2901dfe 100644
--- a/bwamem_pair.c
+++ b/bwamem_pair.c
@@ -289,6 +289,7 @@ int mem_sam_pe(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, co
 		if (n_sub > 0) q_pe -= (int)(4.343 * log(n_sub+1) + .499);
 		if (q_pe < 0) q_pe = 0;
 		if (q_pe > 60) q_pe = 60;
+		q_pe = (int)(q_pe * (1. - .5 * (a[0].a[0].frac_rep + a[1].a[0].frac_rep)) + .499);
 		// the following assumes no split hits
 		if (o > score_un) { // paired alignment is preferred
 			mem_alnreg_t *c[2];
diff --git a/main.c b/main.c
index 9727770..ad4e750 100644
--- a/main.c
+++ b/main.c
@@ -4,7 +4,7 @@
 #include "utils.h"
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.7.8-r758-dirty"
+#define PACKAGE_VERSION "0.7.8-r759-dirty"
 #endif
 
 int bwa_fa2pac(int argc, char *argv[]);