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 44eed4e468bcbad77923719ac486966e4eaef457..f36b0ee2ab443affe0635866a8d593c5cb54fdf7 100644 (file)
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -214,25 +214,22 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
        }
 
        { // fix p[k<<2|k]
-               float max1, max2, min;
+               float max1, max2, min1, min2;
                int max_k, min_k;
                max_k = min_k = -1;
-               max1 = max2 = -1.0; min = 1e30;
+               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];
-                       }
+                       } else if (b->esum[k] > max2) max2 = b->esum[k];
                }
                for (k = 0; k < 4; ++k) {
-                       if (min > p[k<<2|k]) {
-                               min = p[k<<2|k];
-                               min_k = 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)
-                       p[max_k<<2|max_k] = min > 1.0? min - 1.0 : 0.0;
+               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
@@ -442,7 +439,7 @@ 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;
+                               //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]);
                        }
@@ -502,6 +499,15 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                                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);
        }