* samtools-0.1.12-4 (r884)

[samtools.git] / bam2bcf_indel.c

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index 239fb8d731040b4c2eea951040c2eaf8334c5c7e..3a8782b9e473b2fdc76c779468cfa98547d1caa8 100644 (file)
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -113,7 +113,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
        extern void ks_introsort_uint32_t(int, uint32_t*);
        int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score1, *score2, max_ref2;
        int N, K, l_run, ref_type, n_alt;
-       char *inscns = 0, *ref2, *query;
+       char *inscns = 0, *ref2, *query, **ref_sample;
        khash_t(rg) *hash = (khash_t(rg)*)rghash;
        if (ref == 0 || bca == 0) return -1;
        // mark filtered reads
@@ -188,6 +188,50 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                        if (ref[i] == 0) break;
                right = i;
        }
+       /* The following block fixes a long-existing flaw in the INDEL
+        * calling model: the interference of nearby SNPs. However, it also
+        * reduces the power because sometimes, substitutions caused by
+        * indels are not distinguishable from true mutations. Multiple
+        * sequence realignment helps to increase the power.
+        */
+       { // construct per-sample consensus
+               int L = right - left + 1;
+               uint32_t *cns;
+               char *ref0, *r;
+               ref_sample = calloc(n, sizeof(void*));
+               cns = calloc(L, 4);
+               ref0 = calloc(L, 1);
+               for (i = 0; i < right - left; ++i)
+                       ref0[i] = bam_nt16_table[(int)ref[i+left]];
+               for (s = 0; s < n; ++s) {
+                       r = ref_sample[s] = calloc(L, 1);
+                       memset(cns, 0, sizeof(int) * L);
+                       // collect ref and non-ref counts
+                       for (i = 0; i < n_plp[s]; ++i) {
+                               bam_pileup1_t *p = plp[s] + i;
+                               bam1_t *b = p->b;
+                               uint32_t *cigar = bam1_cigar(b);
+                               uint8_t *seq = bam1_seq(b);
+                               int x = b->core.pos, y = 0;
+                               for (k = 0; k < b->core.n_cigar; ++k) {
+                                       int op = cigar[k]&0xf;
+                                       int j, l = cigar[k]>>4;
+                                       if (op == BAM_CMATCH) {
+                                               for (j = 0; j < l; ++j)
+                                                       if (x + j >= left && x + j < right)
+                                                               cns[x+j-left] += (bam1_seqi(seq, y+j) == ref0[x+j-left])? 1 : 0x10000;
+                                               x += l; y += l;
+                                       } else if (op == BAM_CDEL || op == BAM_CREF_SKIP) x += l;
+                                       else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
+                               }
+                       }
+                       // determine the consensus
+                       for (i = 0; i < right - left; ++i)
+                               r[i] = (cns[i] && (double)(cns[i]&0xffff) / ((cns[i]&0xffff)+(cns[i]>>16&0xffff)) < 0.7)? 15 : ref0[i];
+//                     for (i = 0; i < right - left; ++i) fputc("=ACMGRSVTWYHKDBN"[(int)r[i]], stderr); fputc('\n', stderr);
+               }
+               free(ref0); free(cns);
+       }
        { // the length of the homopolymer run around the current position
                int c = bam_nt16_table[(int)ref[pos + 1]];
                if (c == 15) l_run = 1;
@@ -251,23 +295,19 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                else ir = est_indelreg(pos, ref, -types[t], 0);
                if (ir > bca->indelreg) bca->indelreg = ir;
 //             fprintf(stderr, "%d, %d, %d\n", pos, types[t], ir);
-               // write ref2
-               for (k = 0, j = left; j <= pos; ++j)
-                       ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
-               if (types[t] <= 0) j += -types[t];
-               else for (l = 0; l < types[t]; ++l)
-                                ref2[k++] = inscns[t*max_ins + l];
-               if (types[0] < 0) { // mask deleted sequences to avoid a particular error in the model.
-                       int jj, tmp = types[t] >= 0? -types[0] : -types[0] + types[t];
-                       for (jj = 0; jj < tmp && j < right && ref[j]; ++jj, ++j)
-                               ref2[k++] = 4;
-               }
-               for (; j < right && ref[j]; ++j)
-                       ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
-               for (; k < max_ref2; ++k) ref2[k] = 4;
-               if (j < right) right = j;
-               // align each read to ref2
+               // realignment
                for (s = K = 0; s < n; ++s) {
+                       // write ref2
+                       for (k = 0, j = left; j <= pos; ++j)
+                               ref2[k++] = bam_nt16_nt4_table[(int)ref_sample[s][j-left]];
+                       if (types[t] <= 0) j += -types[t];
+                       else for (l = 0; l < types[t]; ++l)
+                                        ref2[k++] = inscns[t*max_ins + l];
+                       for (; j < right && ref[j]; ++j)
+                               ref2[k++] = bam_nt16_nt4_table[(int)ref_sample[s][j-left]];
+                       for (; k < max_ref2; ++k) ref2[k] = 4;
+                       if (j < right) right = j;
+                       // align each read to ref2
                        for (i = 0; i < n_plp[s]; ++i, ++K) {
                                bam_pileup1_t *p = plp[s] + i;
                                int qbeg, qend, tbeg, tend, sc;
@@ -406,6 +446,8 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
        }
        free(score1); free(score2);
        // free
+       for (i = 0; i < n; ++i) free(ref_sample[i]);
+       free(ref_sample);
        free(types); free(inscns);
        return n_alt > 0? 0 : -1;
 }