From 664c90281344b1e5711eee12d9f7b47ba969309a Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Wed, 10 Jun 2009 15:27:59 +0000
Subject: [PATCH]  * samtools-0.1.4-4 (r325)  * further avoid wrong consensus
 calls in repetitive regions.

---
 bam_maqcns.c | 19 ++++++++-----------
 bamtk.c      |  2 +-
 2 files changed, 9 insertions(+), 12 deletions(-)

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 44eed4e..464288a 100644
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -214,25 +214,22 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
 	}
 
 	{ // fix p[k<<2|k]
-		float max1, max2, min;
+		float max1, max2, min1, min2;
 		int max_k, min_k;
 		max_k = min_k = -1;
-		max1 = max2 = -1.0; min = 1e30;
+		max1 = max2 = -1.0; min1 = min2 = 1e30;
 		for (k = 0; k < 4; ++k) {
 			if (b->esum[k] > max1) {
 				max2 = max1; max1 = b->esum[k]; max_k = k;
-			} else if (b->esum[k] > max2) {
-				max2 = b->esum[k];
-			}
+			} else if (b->esum[k] > max2) max2 = b->esum[k];
 		}
 		for (k = 0; k < 4; ++k) {
-			if (min > p[k<<2|k]) {
-				min = p[k<<2|k];
-				min_k = k;
-			}
+			if (p[k<<2|k] < min1) {
+				min2 = min1; min1 = p[k<<2|k]; min_k = k;
+			} else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
 		}
-		if (max1 > max2 && min_k != max_k)
-			p[max_k<<2|max_k] = min > 1.0? min - 1.0 : 0.0;
+		if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
+			p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
 	}
 
 	// convert necessary information to glf1_t
diff --git a/bamtk.c b/bamtk.c
index 8441082..e901961 100644
--- a/bamtk.c
+++ b/bamtk.c
@@ -3,7 +3,7 @@
 #include "bam.h"
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.4-3 (r324)"
+#define PACKAGE_VERSION "0.1.4-4 (r325)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);
-- 
2.39.5