* samtools-0.1.2-16

[samtools.git] / bam_maqcns.c

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 3c3d81892ed95d06e298ee57539b36e3c5a26bf4..6b2b5e03260da470801ce70b8ff18f788c8fd880 100644 (file)
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -4,6 +4,8 @@
 #include "ksort.h"
 KSORT_INIT_GENERIC(uint32_t)
 
+#define MAX_WINDOW 33
+
 typedef struct __bmc_aux_t {
        int max;
        uint32_t *info;
@@ -12,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 };
@@ -127,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;
@@ -153,7 +156,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);
@@ -164,8 +167,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) {
@@ -206,7 +210,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)
@@ -279,19 +283,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) {
@@ -299,9 +307,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;
@@ -323,6 +334,8 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                                if (seg.tend > right) right = seg.tend;
                        }
                }
+               if (pos - left > MAX_WINDOW) left = pos - MAX_WINDOW;
+               if (right - pos> MAX_WINDOW) right = pos + MAX_WINDOW;
        }
        { // the core part
                char *ref2, *inscns = 0;
@@ -383,7 +396,7 @@ 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; // "<" should not happen if there is no bug
+                                       ct = c->pos + l >= left? ref2[c->pos + l - left] : 15; // "<" will happen if reads are too long
                                        if (cq < 15 && ct < 15) {
                                                s += cq == ct? 1 : -mi->mm_penalty;
                                                if (cq != ct) ps += bam1_qual(p->b)[l];
@@ -460,8 +473,9 @@ 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];
-                               ret->gl[0] += s1 < s2? 0 : s1 - s2 < mi->q_indel? s1 - s2 : mi->q_indel;
-                               ret->gl[1] += s2 < s1? 0 : s2 - s1 < mi->q_indel? s2 - s1 : mi->q_indel;
+                               //printf("%d, %d\n", 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;
                        }
                }
                free(score); free(pscore); free(ref2); free(inscns);