* samtools-0.1.2-6

 * improve indel caller by setting maximum window size

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 3c3d81892ed95d06e298ee57539b36e3c5a26bf4..2270c21ac2ce01f1cdf8d15cbfdfa998543f2977 100644 (file)
--- 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 f9269a36205a6267fa3d2dc6c7d58fdecfcfa424..b0eec6ee881e871ddc34cf43aa69ddcc955f40c4 100644 (file)
--- 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 7427f068941cc37d75a1a4eefadc580ca352ee47..d357fccb8f117b830e293707b19fb9b172e92b8b 100755 (executable)
--- 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]) {