* samtools-0.1.10-4 (r834)

 * fixed a silly bug in printing MNP
 * restrict to at most 1 alternative allele

diff --git a/bam2bcf.c b/bam2bcf.c
index de414d7176903dd323ced99f3b8de39396667517..45a0ccc36592edbcd8647e731a6476900c4b2c31 100644 (file)
--- a/bam2bcf.c
+++ b/bam2bcf.c
@@ -215,7 +215,7 @@ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bc
                        if (bc->a[i] < 0) break;
                        if (i > 1) kputc(',', &s);
                        for (j = 0; j < bca->indelreg; ++j)
-                               kputc("ACGT"[bca->indel_types[i]>>2*i&3], &s);
+                               kputc("ACGT"[bca->indel_types[i]>>2*j&3], &s);
                }
                kputc('\0', &s);
        } else { // a SNP

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index 790f59e5e886eed9f1d07b59c816b6188cff71f8..a930cba03b57f5aee2f59e869be04c7ccb530fad 100644 (file)
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -455,7 +455,7 @@ int bcf_call_mnp_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                aux = calloc(N + 1, 4);
                for (i = 0, x = 0; i < B2B_MAX_MNP; ++i)
                        x |= bam_nt16_nt4_table[bam_nt16_table[(int)ref[pos+i]]] << 2*i;
-               ref_seq = aux[m++] = x;
+               ref_seq = x;
                bca->indelreg = 0;
                for (s = 0; s < n; ++s) {
                        for (i = 0; i < n_plp[s]; ++i) {
@@ -472,7 +472,7 @@ int bcf_call_mnp_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                                        else stop = 1;
                                        x |= c << 2*k;
                                }
-//                             fprintf(stderr, "* %d, %d\n", x, j);
+//                             if (j >= 2 && k == B2B_MAX_MNP) fprintf(stderr, "* %c%c%c%c, %d\n", "ACGT"[x&3], "ACGT"[x>>2&3], "ACGT"[x>>4&3], "ACGT"[x>>6&3], j);
                                if (k == B2B_MAX_MNP && j >= 2) aux[m++] = x;
                        }
                }
@@ -495,12 +495,9 @@ int bcf_call_mnp_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                free(aux);
                // collect types
                ks_introsort(uint32_t, n_types, cnt);
-               for (t = 0; t < n_types; ++t)
-                       if ((cnt[t]&0xff) == ref_seq) break;
-               if (t == n_types - 1) --t; // we get at least two types
-               types = (int*)calloc(n_types, sizeof(int));
-               for (i = 0; t < n_types; ++i, ++t) types[i] = cnt[t]&0xff;
-               n_types = i;
+               types = (int*)calloc(2, sizeof(int));
+               types[0] = ref_seq; types[1] = cnt[n_types-1]&0xff;
+               ref_type = 0; n_types = 2;
                // calculate MNP length
                for (i = B2B_MAX_MNP - 1, k = 0; i >= 0; --i) {
                        int c = types[0] >> 2*i & 3;
@@ -513,10 +510,7 @@ int bcf_call_mnp_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                x = (1<<2*bca->indelreg) - 1;
                for (t = 0; t < n_types; ++t) types[t] &= x;
                ref_seq &= x;
-               for (t = 0; t < n_types; ++t)
-                       if (types[t] == ref_seq) break;
-               ref_type = t;
-//             fprintf(stderr, "%d, %d, [%x,%x] %d\n", pos, n_types, types[0], types[1], bca->indelreg);
+//             x = types[1]; fprintf(stderr, "%d, %d, %d, %c%c%c%c\n", pos, n_types, bca->indelreg, "ACGT"[x&3], "ACGT"[x>>2&3], "ACGT"[x>>4&3], "ACGT"[x>>6&3]);
        }
        { // calculate left and right boundary
                left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;

diff --git a/bamtk.c b/bamtk.c
index 89f7e80d8b12cf26fb4d450eb5f695c1ab2bb242..ba02dd44f674bf31c4806825e24253af91707b09 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,7 @@
 #endif
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.10-3 (r833)"
+#define PACKAGE_VERSION "0.1.10-4 (r834)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);