* samtools-0.1.3-9 (r256)

[samtools.git] / bam_maqcns.c

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 2270c21ac2ce01f1cdf8d15cbfdfa998543f2977..482f33ca5bbfec2af6162518d21936e86105c058 100644 (file)
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -14,7 +14,7 @@ typedef struct __bmc_aux_t {
 typedef struct {
        float esum[4], fsum[4];
        uint32_t c[4];
-       uint32_t mapQ_max;
+       uint32_t rms_mapQ;
 } glf_call_aux_t;
 
 char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
@@ -129,6 +129,7 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
        int i, j, k, w[8], c, n;
        glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
        float p[16], min_p = 1e30;
+       uint64_t rms;
 
        g->ref_base = ref_base;
        if (_n == 0) return g;
@@ -141,13 +142,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;
        }
@@ -155,7 +157,7 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
        // generate esum and fsum
        b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
        for (k = 0; k != 8; ++k) w[k] = 0;
-       b->mapQ_max = 0;
+       rms = 0;
        for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
                uint32_t info = bm->aux->info[j];
                if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
@@ -166,8 +168,9 @@ 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];
                }
-               if (b->mapQ_max < (info&0x7f)) b->mapQ_max = info&0x7f;
+               rms += (int)(info&0x7f) * (info&0x7f);
        }
+       b->rms_mapQ = (uint8_t)(sqrt((double)rms / n) + .499);
        // rescale ->c[]
        for (j = c = 0; j != 4; ++j) c += b->c[j];
        if (c > 255) {
@@ -208,7 +211,7 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
        }
 
        // convert necessary information to glf1_t
-       g->ref_base = ref_base; g->max_mapQ = b->mapQ_max;
+       g->ref_base = ref_base; g->max_mapQ = b->rms_mapQ;
        g->depth = n > 16777215? 16777215 : n;
        for (j = 0; j != 4; ++j)
                for (k = j; k < 4; ++k)
@@ -281,19 +284,23 @@ void bam_maqindel_ret_destroy(bam_maqindel_ret_t *mir)
 
 #define MINUS_CONST 0x10000000
 
-bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref)
+bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
+                                                                int _n_types, int *_types)
 {
        int i, j, n_types, *types, left, right;
        bam_maqindel_ret_t *ret = 0;
-       for (i = 0; i < n; ++i) {
-               const bam_pileup1_t *p = pl + i;
-               if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
+       // if there is no proposed indel, check if there is an indel from the alignment
+       if (_n_types == 0) {
+               for (i = 0; i < n; ++i) {
+                       const bam_pileup1_t *p = pl + i;
+                       if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
+               }
+               if (i == n) return 0; // no indel
        }
-       if (i == n) return 0; // no indel
        { // calculate how many types of indels are available (set n_types and types)
                int m;
                uint32_t *aux;
-               aux = (uint32_t*)calloc(n+1, 4);
+               aux = (uint32_t*)calloc(n + _n_types + 1, 4);
                m = 0;
                aux[m++] = MINUS_CONST; // zero indel is always a type
                for (i = 0; i < n; ++i) {
@@ -301,9 +308,12 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                        if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
                                aux[m++] = MINUS_CONST + p->indel;
                }
+               if (_n_types) // then also add this to aux[]
+                       for (i = 0; i < _n_types; ++i)
+                               if (_types[i]) aux[m++] = MINUS_CONST + _types[i];
                ks_introsort(uint32_t, m, aux);
-               n_types = 1;
-               for (i = 1; i < m; ++i)
+               // squeeze out identical types
+               for (i = 1, n_types = 1; i < m; ++i)
                        if (aux[i] != aux[i-1]) ++n_types;
                types = (int*)calloc(n_types, sizeof(int));
                j = 0;
@@ -451,14 +461,11 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                                        ret->s[1][k+1] = ref[pos + k + 1];
                        } else ret->s[1][0] = '*';
                        // write count
-                       for (j = 0; j < n; ++j) {
-                               if (score[max1_i*n+j] < 0 && score[max2_i*n+j] < 0) ++ret->cnt_anti;
-                               else {
-                                       int diff = score[max1_i*n+j] - score[max2_i*n+j];
-                                       if (diff > mi->ambi_thres) ++ret->cnt1;
-                                       else if (diff < -mi->ambi_thres) ++ret->cnt2;
-                                       else ++ret->cnt_ambi;
-                               }
+                       for (i = 0; i < n; ++i) {
+                               const bam_pileup1_t *p = pl + i;
+                               if (p->indel == ret->indel1) ++ret->cnt1;
+                               else if (p->indel == ret->indel2) ++ret->cnt2;
+                               else ++ret->cnt_anti;
                        }
                        // write gl[]
                        ret->gl[0] = ret->gl[1] = 0;