* samtools-0.1.3-9 (r256)

[samtools.git] / bam_maqcns.c

diff --git a/bam_maqcns.c b/bam_maqcns.c
index ec11f3ede314d025713e9b268136d0172dc5e209..482f33ca5bbfec2af6162518d21936e86105c058 100644 (file)
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -142,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;
        }
@@ -283,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) {
@@ -303,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;
@@ -453,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;