* samtools-0.1.12-4 (r884)

 * fixed a long-existing flaw in the INDEL calling model

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index 239fb8d731040b4c2eea951040c2eaf8334c5c7e..3a8782b9e473b2fdc76c779468cfa98547d1caa8 100644 (file)
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -113,7 +113,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
        extern void ks_introsort_uint32_t(int, uint32_t*);
        int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score1, *score2, max_ref2;
        int N, K, l_run, ref_type, n_alt;
-       char *inscns = 0, *ref2, *query;
+       char *inscns = 0, *ref2, *query, **ref_sample;
        khash_t(rg) *hash = (khash_t(rg)*)rghash;
        if (ref == 0 || bca == 0) return -1;
        // mark filtered reads
@@ -188,6 +188,50 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                        if (ref[i] == 0) break;
                right = i;
        }
+       /* The following block fixes a long-existing flaw in the INDEL
+        * calling model: the interference of nearby SNPs. However, it also
+        * reduces the power because sometimes, substitutions caused by
+        * indels are not distinguishable from true mutations. Multiple
+        * sequence realignment helps to increase the power.
+        */
+       { // construct per-sample consensus
+               int L = right - left + 1;
+               uint32_t *cns;
+               char *ref0, *r;
+               ref_sample = calloc(n, sizeof(void*));
+               cns = calloc(L, 4);
+               ref0 = calloc(L, 1);
+               for (i = 0; i < right - left; ++i)
+                       ref0[i] = bam_nt16_table[(int)ref[i+left]];
+               for (s = 0; s < n; ++s) {
+                       r = ref_sample[s] = calloc(L, 1);
+                       memset(cns, 0, sizeof(int) * L);
+                       // collect ref and non-ref counts
+                       for (i = 0; i < n_plp[s]; ++i) {
+                               bam_pileup1_t *p = plp[s] + i;
+                               bam1_t *b = p->b;
+                               uint32_t *cigar = bam1_cigar(b);
+                               uint8_t *seq = bam1_seq(b);
+                               int x = b->core.pos, y = 0;
+                               for (k = 0; k < b->core.n_cigar; ++k) {
+                                       int op = cigar[k]&0xf;
+                                       int j, l = cigar[k]>>4;
+                                       if (op == BAM_CMATCH) {
+                                               for (j = 0; j < l; ++j)
+                                                       if (x + j >= left && x + j < right)
+                                                               cns[x+j-left] += (bam1_seqi(seq, y+j) == ref0[x+j-left])? 1 : 0x10000;
+                                               x += l; y += l;
+                                       } else if (op == BAM_CDEL || op == BAM_CREF_SKIP) x += l;
+                                       else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
+                               }
+                       }
+                       // determine the consensus
+                       for (i = 0; i < right - left; ++i)
+                               r[i] = (cns[i] && (double)(cns[i]&0xffff) / ((cns[i]&0xffff)+(cns[i]>>16&0xffff)) < 0.7)? 15 : ref0[i];
+//                     for (i = 0; i < right - left; ++i) fputc("=ACMGRSVTWYHKDBN"[(int)r[i]], stderr); fputc('\n', stderr);
+               }
+               free(ref0); free(cns);
+       }
        { // the length of the homopolymer run around the current position
                int c = bam_nt16_table[(int)ref[pos + 1]];
                if (c == 15) l_run = 1;
@@ -251,23 +295,19 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                else ir = est_indelreg(pos, ref, -types[t], 0);
                if (ir > bca->indelreg) bca->indelreg = ir;
 //             fprintf(stderr, "%d, %d, %d\n", pos, types[t], ir);
-               // write ref2
-               for (k = 0, j = left; j <= pos; ++j)
-                       ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
-               if (types[t] <= 0) j += -types[t];
-               else for (l = 0; l < types[t]; ++l)
-                                ref2[k++] = inscns[t*max_ins + l];
-               if (types[0] < 0) { // mask deleted sequences to avoid a particular error in the model.
-                       int jj, tmp = types[t] >= 0? -types[0] : -types[0] + types[t];
-                       for (jj = 0; jj < tmp && j < right && ref[j]; ++jj, ++j)
-                               ref2[k++] = 4;
-               }
-               for (; j < right && ref[j]; ++j)
-                       ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
-               for (; k < max_ref2; ++k) ref2[k] = 4;
-               if (j < right) right = j;
-               // align each read to ref2
+               // realignment
                for (s = K = 0; s < n; ++s) {
+                       // write ref2
+                       for (k = 0, j = left; j <= pos; ++j)
+                               ref2[k++] = bam_nt16_nt4_table[(int)ref_sample[s][j-left]];
+                       if (types[t] <= 0) j += -types[t];
+                       else for (l = 0; l < types[t]; ++l)
+                                        ref2[k++] = inscns[t*max_ins + l];
+                       for (; j < right && ref[j]; ++j)
+                               ref2[k++] = bam_nt16_nt4_table[(int)ref_sample[s][j-left]];
+                       for (; k < max_ref2; ++k) ref2[k] = 4;
+                       if (j < right) right = j;
+                       // align each read to ref2
                        for (i = 0; i < n_plp[s]; ++i, ++K) {
                                bam_pileup1_t *p = plp[s] + i;
                                int qbeg, qend, tbeg, tend, sc;
@@ -406,6 +446,8 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
        }
        free(score1); free(score2);
        // free
+       for (i = 0; i < n; ++i) free(ref_sample[i]);
+       free(ref_sample);
        free(types); free(inscns);
        return n_alt > 0? 0 : -1;
 }

diff --git a/bamtk.c b/bamtk.c
index da54ffd9d65c5a4afea3d54426cd36830eda3897..613e76cf0d901e407eae38689e833838d9c43bf7 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,7 @@
 #endif
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.12-3 (r878)"
+#define PACKAGE_VERSION "0.1.12-4 (r884)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);

diff --git a/bcftools/vcfutils.pl b/bcftools/vcfutils.pl
index 5df8b0c05e1275e7edc4b8063efa5c585ee7a80a..e38f1ca5ef660eca6aa223d4ab4cfbd359408a2f 100755 (executable)
--- a/bcftools/vcfutils.pl
+++ b/bcftools/vcfutils.pl
@@ -215,7 +215,7 @@ Note: This command discards indels. Output: QUAL #non-indel #SNPs #transitions #
 }
 
 sub varFilter {
-  my %opts = (d=>2, D=>10000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000);
+  my %opts = (d=>2, D=>10000000, a=>2, W=>10, Q=>10, w=>10, p=>undef, 1=>1e-4, 2=>1e-100, 3=>0, 4=>1e-4, G=>0, S=>1000);
   getopts('pd:D:W:Q:w:a:1:2:3:4:G:S:', \%opts);
   die(qq/
 Usage:   vcfutils.pl varFilter [options] <in.vcf>
@@ -289,7 +289,7 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
        $flt = 1 if ($flt == 0 && $mq >= 0 && $mq < $opts{Q});
        $flt = 7 if ($flt == 0 && /PV4=([^,]+),([^,]+),([^,]+),([^,;\t]+)/
                                 && ($1<$opts{1} || $2<$opts{2} || $3<$opts{3} || $4<$opts{4}));
-       $flt = 8 if ($flt == 0 && ((/MXGQ=(\d+)/ && $1 >= $opts{G}) || (/MXSP=(\d+)/ && $1 < $opts{S})));
+       $flt = 8 if ($flt == 0 && ((/MXGQ=(\d+)/ && $1 < $opts{G}) || (/MXSP=(\d+)/ && $1 >= $opts{S})));
 
        my $score = $t[5] * 100 + $dp_alt;
        my $rlen = length($t[3]) - 1; # $indel_score<0 for SNPs
@@ -312,7 +312,10 @@ Note: Some of the filters rely on annotations generated by SAMtools\/BCFtools.
          } else { # SNP or MNP
                for my $x (@staging) {
                  next if (($x->[0]&3) != 3 || $x->[4] + $x->[2] + $ow < $t[1]);
-                 $flt = 5;
+                 if ($x->[4] + length($x->[7]) - 1 == $t[1] && substr($x->[7], -1, 1) eq substr($t[4], 0, 1)
+                         && length($x->[7]) - length($x->[6]) == 1) {
+                       $x->[1] = 5;
+                 } else { $flt = 5; }
                  last;
                }
                # check MNP