A bug fix, "samtools view" is now working again.

[samtools.git] / bam_maqcns.c

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 60a194525425d711a23a286d6f411358911173f1..f36b0ee2ab443affe0635866a8d593c5cb54fdf7 100644 (file)
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -108,6 +108,7 @@ bam_maqcns_t *bam_maqcns_init()
        bm->theta = 0.85;
        bm->n_hap = 2;
        bm->eta = 0.03;
+       bm->cap_mapQ = 60;
        return bm;
 }
 
@@ -142,13 +143,14 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
        }
        for (i = n = 0; i < _n; ++i) {
                const bam_pileup1_t *p = pl + i;
-               uint32_t q, x = 0;
+               uint32_t q, x = 0, qq;
                if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue;
                q = (uint32_t)bam1_qual(p->b)[p->qpos];
                x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
                if (p->b->core.qual < q) q = p->b->core.qual;
                x |= q << 24;
-               q = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos)];
+               qq = bam1_seqi(bam1_seq(p->b), p->qpos);
+               q = bam_nt16_nt4_table[qq? qq : ref_base];
                if (!p->is_del && q < 4) x |= 1 << 21 | q << 16;
                bm->aux->info[n++] = x;
        }
@@ -159,6 +161,7 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
        rms = 0;
        for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
                uint32_t info = bm->aux->info[j];
+               int tmp;
                if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
                k = info>>16&7;
                if (info>>24 > 0) {
@@ -167,7 +170,8 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
                        if (w[k] < 0xff) ++w[k];
                        ++b->c[k&3];
                }
-               rms += (int)(info&0x7f) * (info&0x7f);
+               tmp = (int)(info&0x7f) < bm->cap_mapQ? (int)(info&0x7f) : bm->cap_mapQ;
+               rms += tmp * tmp;
        }
        b->rms_mapQ = (uint8_t)(sqrt((double)rms / n) + .499);
        // rescale ->c[]
@@ -209,6 +213,25 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
                        if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
        }
 
+       { // fix p[k<<2|k]
+               float max1, max2, min1, min2;
+               int max_k, min_k;
+               max_k = min_k = -1;
+               max1 = max2 = -1.0; min1 = min2 = 1e30;
+               for (k = 0; k < 4; ++k) {
+                       if (b->esum[k] > max1) {
+                               max2 = max1; max1 = b->esum[k]; max_k = k;
+                       } else if (b->esum[k] > max2) max2 = b->esum[k];
+               }
+               for (k = 0; k < 4; ++k) {
+                       if (p[k<<2|k] < min1) {
+                               min2 = min1; min1 = p[k<<2|k]; min_k = k;
+                       } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
+               }
+               if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
+                       p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
+       }
+
        // convert necessary information to glf1_t
        g->ref_base = ref_base; g->max_mapQ = b->rms_mapQ;
        g->depth = n > 16777215? 16777215 : n;
@@ -396,7 +419,8 @@ 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; // "<" will happen if reads are too long
+                                       // in the following line, "<" will happen if reads are too long
+                                       ct = c->pos + l - seg.qbeg >= left? ref2[c->pos + l - seg.qbeg - left] : 15;
                                        if (cq < 15 && ct < 15) {
                                                s += cq == ct? 1 : -mi->mm_penalty;
                                                if (cq != ct) ps += bam1_qual(p->b)[l];
@@ -406,7 +430,7 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                                if (types[i] != 0) { // then try the other way to calculate the score
                                        for (ps = s = 0, l = seg.qbeg; c->pos + l + types[i] < right && l < seg.qend; ++l) {
                                                int cq = bam1_seqi(bam1_seq(p->b), l), ct;
-                                               ct = c->pos + l + types[i] >= left? ref2[c->pos + l + types[i] - left] : 15;
+                                               ct = c->pos + l - seg.qbeg + types[i] >= left? ref2[c->pos + l - seg.qbeg + types[i] - left] : 15;
                                                if (cq < 15 && ct < 15) {
                                                        s += cq == ct? 1 : -mi->mm_penalty;
                                                        if (cq != ct) ps += bam1_qual(p->b)[l];
@@ -415,8 +439,9 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                                }
                                if (score[i*n+j] < s) score[i*n+j] = s; // choose the higher of the two scores
                                if (pscore[i*n+j] > ps) pscore[i*n+j] = ps;
-                               if (types[i] != 0) score[i*n+j] -= mi->indel_err;
-                               //printf("%d, %d, %d, %d\n", i, types[i], j, score[i*n+j]);
+                               //if (types[i] != 0) score[i*n+j] -= mi->indel_err;
+                               //printf("%d, %d, %d, %d, %d, %d, %d\n", p->b->core.pos + 1, seg.qbeg, i, types[i], j,
+                               //         score[i*n+j], pscore[i*n+j]);
                        }
                }
                { // get final result
@@ -470,10 +495,19 @@ 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];
-                               //printf("%d, %d\n", s1, s2);
+                               //printf("%d, %d, %d, %d, %d\n", pl[j].b->core.pos+1, max1_i, max2_i, 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;
                        }
+                       // write cnt_ref and cnt_ambi
+                       if (max1_i != 0 && max2_i != 0) {
+                               for (j = 0; j < n; ++j) {
+                                       int diff1 = score[j] - score[max1_i * n + j];
+                                       int diff2 = score[j] - score[max2_i * n + j];
+                                       if (diff1 > 0 && diff2 > 0) ++ret->cnt_ref;
+                                       else if (diff1 == 0 || diff2 == 0) ++ret->cnt_ambi;
+                               }
+                       }
                }
                free(score); free(pscore); free(ref2); free(inscns);
        }