* samtools-0.1.9-8 (r809)

 * fixed a bug in the indel caller

diff --git a/bam2bcf.c b/bam2bcf.c
index 1ba288ccd949e3b506b7c786f8f609480af05cb2..8537e232191189ee489ba7ce681008c12a39b829 100644 (file)
--- a/bam2bcf.c
+++ b/bam2bcf.c
@@ -124,7 +124,10 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/,
                if (((ref4 < 4 && j < 4) || (ref4 == 4 && j < 5)) && i >= 0)
                        call->unseen = j, call->a[j++] = qsum[i]&3;
                call->n_alleles = j;
-       } else call->n_alleles = j;
+       } else {
+               call->n_alleles = j;
+               if (call->n_alleles == 1) return -1; // no reliable supporting read. stop doing anything
+       }
        // set the PL array
        if (call->n < n) {
                call->n = n;

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index 2d0b5c6f8a66a4c5c006dbec560718ab1f7a0377..c28b6859dd0101dc5f9ded34a42ded1add58f781 100644 (file)
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -110,7 +110,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                                          const void *rghash)
 {
        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, *score, N, K, l_run, ref_type;
+       int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type, n_alt;
        char *inscns = 0, *ref2, *query;
        khash_t(rg) *hash = (khash_t(rg)*)rghash;
        if (ref == 0 || bca == 0) return -1;
@@ -342,13 +342,14 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                        memcpy(&bca->inscns[t * bca->maxins], &inscns[(sumq[t]&0x3f) * max_ins], bca->maxins);
                }
                // update p->aux
-               for (s = K = 0; s < n; ++s) {
-                       for (i = 0; i < n_plp[s]; ++i, ++K) {
+               for (s = n_alt = 0; s < n; ++s) {
+                       for (i = 0; i < n_plp[s]; ++i) {
                                bam_pileup1_t *p = plp[s] + i;
                                int x = types[p->aux>>16&0x3f];
                                for (j = 0; j < 4; ++j)
                                        if (x == bca->indel_types[j]) break;
                                p->aux = j<<16 | (j == 4? 0 : (p->aux&0xffff));
+                               if ((p->aux>>16&0x3f) > 0) ++n_alt;
 //                             fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d\n", pos, s, i, bam1_qname(p->b), p->aux>>16&63, p->aux&0xff);
                        }
                }               
@@ -356,5 +357,5 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
        free(score);
        // free
        free(types); free(inscns);
-       return 0;
+       return n_alt > 0? 0 : -1;
 }

diff --git a/bam_plcmd.c b/bam_plcmd.c
index 0fbc9303c25744ebde8d1b0e269c539d5c4f5325..1dc819a6d8870782dd94eb5d446821c4622d4b94 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -739,12 +739,13 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                        if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref, rghash) >= 0) {
                                for (i = 0; i < gplp.n; ++i)
                                        bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], -1, bca, bcr + i);
-                               bcf_call_combine(gplp.n, bcr, -1, &bc);
-                               b = calloc(1, sizeof(bcf1_t));
-                               bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0,
-                                                        (conf->flag&MPLP_FMT_SP), bca, ref);
-                               bcf_write(bp, bh, b);
-                               bcf_destroy(b);
+                               if (bcf_call_combine(gplp.n, bcr, -1, &bc) >= 0) {
+                                       b = calloc(1, sizeof(bcf1_t));
+                                       bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0,
+                                                                (conf->flag&MPLP_FMT_SP), bca, ref);
+                                       bcf_write(bp, bh, b);
+                                       bcf_destroy(b);
+                               }
                        }
                } else {
                        printf("%s\t%d\t%c", h->target_name[tid], pos + 1, (ref && pos < ref_len)? ref[pos] : 'N');

diff --git a/bamtk.c b/bamtk.c
index d2817b9902a5bccf985bb8360775029bef946267..04a7899bceade9d1ef37579c40007b9a2348ad51 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,7 @@
 #endif
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.9-8 (r807)"
+#define PACKAGE_VERSION "0.1.9-8 (r809)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);