* samtools-0.1.9-6 (r803)

[samtools.git] / bam2bcf.c

diff --git a/bam2bcf.c b/bam2bcf.c
index be8cf1fe529a735800947b80b2a52c7fbaadb2d4..1ba288ccd949e3b506b7c786f8f609480af05cb2 100644 (file)
--- a/bam2bcf.c
+++ b/bam2bcf.c
@@ -20,7 +20,7 @@ bcf_callaux_t *bcf_call_init(double theta, int min_baseQ)
        if (theta <= 0.) theta = CALL_DEFTHETA;
        bca = calloc(1, sizeof(bcf_callaux_t));
        bca->capQ = 60;
-       bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 80;
+       bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 100;
        bca->min_baseQ = min_baseQ;
        bca->e = errmod_init(1. - theta);
        return bca;
@@ -53,11 +53,13 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
        memset(r, 0, sizeof(bcf_callret1_t));
        for (i = n = 0; i < _n; ++i) {
                const bam_pileup1_t *p = pl + i;
-               int q, b, mapQ, baseQ, is_diff, min_dist;
+               int q, b, mapQ, baseQ, is_diff, min_dist, seqQ;
                // set base
                if (p->is_del || p->is_refskip || (p->b->core.flag&BAM_FUNMAP)) continue;
                baseQ = q = is_indel? p->aux&0xff : (int)bam1_qual(p->b)[p->qpos]; // base/indel quality
+               seqQ = is_indel? (p->aux>>8&0xff) : 99;
                if (q < bca->min_baseQ) continue;
+               if (q > seqQ) q = seqQ;
                mapQ = p->b->core.qual < bca->capQ? p->b->core.qual : bca->capQ;
                if (q > mapQ) q = mapQ;
                if (q > 63) q = 63;
@@ -67,7 +69,7 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
                        b = bam_nt16_nt4_table[b? b : ref_base]; // b is the 2-bit base
                        is_diff = (ref4 < 4 && b == ref4)? 0 : 1;
                } else {
-                       b = p->aux>>8&0x3f;
+                       b = p->aux>>16&0x3f;
                        is_diff = (b != 0);
                }
                bca->bases[n++] = q<<5 | (int)bam1_strand(p->b)<<4 | b;
@@ -181,9 +183,11 @@ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bc
                kputc(ref[pos], &s);
                for (i = 1; i < 4; ++i) {
                        if (bc->a[i] < 0) break;
-                       if (i > 1) kputc(',', &s);
+                       if (i > 1) {
+                               kputc(',', &s); kputc(ref[pos], &s);
+                       }
                        if (bca->indel_types[bc->a[i]] < 0) { // deletion
-                               for (j = -bca->indel_types[bc->a[i]]; j < bca->indelreg; ++i)
+                               for (j = -bca->indel_types[bc->a[i]]; j < bca->indelreg; ++j)
                                        kputc(ref[pos+1+j], &s);
                        } else { // insertion; cannot be a reference unless a bug
                                char *inscns = &bca->inscns[bc->a[i] * bca->maxins];