Introduced LIBPATH variable so this could be overridden to allow samtools to build...

[samtools.git] / bam_maqcns.c

diff --git a/bam_maqcns.c b/bam_maqcns.c
index c6c9af7dee454ed91cd75363470473405d7c97ea..464288ae7da0ac74f1443bbb55c03684927cb391 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;
 }
 
@@ -160,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) {
@@ -168,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[]
@@ -210,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;