* fixed a silly bug in printing MNP
* restrict to at most 1 alternative allele
if (bc->a[i] < 0) break;
if (i > 1) kputc(',', &s);
for (j = 0; j < bca->indelreg; ++j)
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
}
kputc('\0', &s);
} else { // a SNP
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;
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;
bca->indelreg = 0;
for (s = 0; s < n; ++s) {
for (i = 0; i < n_plp[s]; ++i) {
bca->indelreg = 0;
for (s = 0; s < n; ++s) {
for (i = 0; i < n_plp[s]; ++i) {
else stop = 1;
x |= c << 2*k;
}
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;
}
}
if (k == B2B_MAX_MNP && j >= 2) aux[m++] = x;
}
}
free(aux);
// collect types
ks_introsort(uint32_t, n_types, cnt);
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;
// calculate MNP length
for (i = B2B_MAX_MNP - 1, k = 0; i >= 0; --i) {
int c = types[0] >> 2*i & 3;
x = (1<<2*bca->indelreg) - 1;
for (t = 0; t < n_types; ++t) types[t] &= x;
ref_seq &= x;
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;
}
{ // calculate left and right boundary
left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
#endif
#ifndef PACKAGE_VERSION
#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[]);
#endif
int bam_taf2baf(int argc, char *argv[]);