From 091c8d2918d888c664dd8baaf9f2a84b2d2b6cd3 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Thu, 12 Feb 2009 11:15:32 +0000
Subject: [PATCH]  * samtools-0.1.2-6  * improve indel caller by setting
 maximum window size

---
 bam_maqcns.c         | 11 ++++++++---
 bamtk.c              |  2 +-
 misc/indel_filter.pl |  2 +-
 3 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 3c3d818..2270c21 100644
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -4,6 +4,8 @@
 #include "ksort.h"
 KSORT_INIT_GENERIC(uint32_t)
 
+#define MAX_WINDOW 33
+
 typedef struct __bmc_aux_t {
 	int max;
 	uint32_t *info;
@@ -323,6 +325,8 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
 				if (seg.tend > right) right = seg.tend;
 			}
 		}
+		if (pos - left > MAX_WINDOW) left = pos - MAX_WINDOW;
+		if (right - pos> MAX_WINDOW) right = pos + MAX_WINDOW;
 	}
 	{ // the core part
 		char *ref2, *inscns = 0;
@@ -383,7 +387,7 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
 				bam_segreg(pos, c, cigar, &seg);
 				for (ps = s = 0, l = seg.qbeg; c->pos + l < right && l < seg.qend; ++l) {
 					int cq = bam1_seqi(bam1_seq(p->b), l), ct;
-					ct = c->pos + l >= left? ref2[c->pos + l - left] : 15; // "<" should not happen if there is no bug
+					ct = c->pos + l >= left? ref2[c->pos + l - left] : 15; // "<" will happen if reads are too long
 					if (cq < 15 && ct < 15) {
 						s += cq == ct? 1 : -mi->mm_penalty;
 						if (cq != ct) ps += bam1_qual(p->b)[l];
@@ -460,8 +464,9 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
 			ret->gl[0] = ret->gl[1] = 0;
 			for (j = 0; j < n; ++j) {
 				int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
-				ret->gl[0] += s1 < s2? 0 : s1 - s2 < mi->q_indel? s1 - s2 : mi->q_indel;
-				ret->gl[1] += s2 < s1? 0 : s2 - s1 < mi->q_indel? s2 - s1 : mi->q_indel;
+				//printf("%d, %d\n", s1, s2);
+				if (s1 > s2) ret->gl[0] += s1 - s2 < mi->q_indel? s1 - s2 : mi->q_indel;
+				else ret->gl[1] += s2 - s1 < mi->q_indel? s2 - s1 : mi->q_indel;
 			}
 		}
 		free(score); free(pscore); free(ref2); free(inscns);
diff --git a/bamtk.c b/bamtk.c
index f9269a3..b0eec6e 100644
--- a/bamtk.c
+++ b/bamtk.c
@@ -3,7 +3,7 @@
 #include "bam.h"
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.2-5"
+#define PACKAGE_VERSION "0.1.2-6"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);
diff --git a/misc/indel_filter.pl b/misc/indel_filter.pl
index 7427f06..d357fcc 100755
--- a/misc/indel_filter.pl
+++ b/misc/indel_filter.pl
@@ -25,7 +25,7 @@ while (<>) {
 	next if ($t[3] eq '*/*');
 	next if ($t[5] == 0);
   }
-  next if ($t[8] + $t[9] + $t[10] + $t[11] > $opts{D});
+  next if ($t[7] > $opts{D});
   @$curr = ($t[0], $t[1], $t[5], $_);
   my $do_swap = 1;
   if (defined $last->[0]) {
-- 
2.39.5