X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fbcf.c;h=0524408da6cf29f4deb6cee09c69d3dce61c4841;hb=ccb7838cb53bf0ca6917a77f7991d940057c12db;hp=35c3d943e04a62fa65387394d141946d7e16b47b;hpb=82e06ed4289c179c7a0da9b8ffa31b15836ae763;p=samtools.git diff --git a/bcftools/bcf.c b/bcftools/bcf.c index 35c3d94..0524408 100644 --- a/bcftools/bcf.c +++ b/bcftools/bcf.c @@ -13,7 +13,6 @@ bcf_t *bcf_open(const char *fn, const char *mode) } else { b->fp = strcmp(fn, "-")? bgzf_open(fn, mode) : bgzf_fdopen(fileno(stdin), mode); } - b->fp->owned_file = 1; return b; } @@ -101,10 +100,16 @@ int bcf_sync(bcf1_t *b) ks_tokaux_t aux; // set ref, alt, flt, info, fmt b->ref = b->alt = b->flt = b->info = b->fmt = 0; - for (p = b->str, n = 0; p < b->str + b->l_str; ++p) - if (*p == 0 && p+1 != b->str + b->l_str) tmp[n++] = p + 1; + for (p = b->str, n = 0; p < b->str + b->l_str; ++p) { + if (*p == 0 && p+1 != b->str + b->l_str) { + if (n == 5) { + ++n; + break; + } else tmp[n++] = p + 1; + } + } if (n != 5) { - fprintf(stderr, "[%s] incorrect number of fields (%d != 5). Corrupted file?\n", __func__, n); + fprintf(stderr, "[%s] incorrect number of fields (%d != 5) at %d:%d\n", __func__, n, b->tid, b->pos); return -1; } b->ref = tmp[0]; b->alt = tmp[1]; b->flt = tmp[2]; b->info = tmp[3]; b->fmt = tmp[4]; @@ -132,10 +137,12 @@ int bcf_sync(bcf1_t *b) for (i = 0; i < b->n_gi; ++i) { if (b->gi[i].fmt == bcf_str2int("PL", 2)) { b->gi[i].len = b->n_alleles * (b->n_alleles + 1) / 2; - } else if (b->gi[i].fmt == bcf_str2int("DP", 2) || b->gi[i].fmt == bcf_str2int("HQ", 2)) { + } else if (b->gi[i].fmt == bcf_str2int("DP", 2) || b->gi[i].fmt == bcf_str2int("HQ", 2) || b->gi[i].fmt == bcf_str2int("DV", 2)) { b->gi[i].len = 2; } else if (b->gi[i].fmt == bcf_str2int("GQ", 2) || b->gi[i].fmt == bcf_str2int("GT", 2)) { b->gi[i].len = 1; + } else if (b->gi[i].fmt == bcf_str2int("SP", 2)) { + b->gi[i].len = 4; } else if (b->gi[i].fmt == bcf_str2int("GL", 2)) { b->gi[i].len = b->n_alleles * (b->n_alleles + 1) / 2 * 4; } @@ -177,7 +184,7 @@ int bcf_read(bcf_t *bp, const bcf_hdr_t *h, bcf1_t *b) } bgzf_read(bp->fp, b->str, b->l_str); l = 12 + b->l_str; - bcf_sync(b); + if (bcf_sync(b) < 0) return -2; for (i = 0; i < b->n_gi; ++i) { bgzf_read(bp->fp, b->gi[i].data, b->gi[i].len * h->n_smpl); l += b->gi[i].len * h->n_smpl; @@ -223,7 +230,31 @@ void bcf_fmt_core(const bcf_hdr_t *h, bcf1_t *b, kstring_t *s) } x = b->n_alleles * (b->n_alleles + 1) / 2; if (b->n_gi == 0) return; + int iPL = -1; + if ( b->n_alleles > 2 ) { + for (i=0; in_gi; i++) { + if ( b->gi[i].fmt == bcf_str2int("PL", 2) ) { + iPL = i; + break; + } + } + } for (j = 0; j < h->n_smpl; ++j) { + + // Determine GT with maximum PL (multiple ALT sites only) + int imax=-1; + if ( iPL!=-1 ) { + uint8_t *d = (uint8_t*)b->gi[iPL].data + j * x; + int k,identical=1; + imax=0; + for (k=1; kn_gi; ++i) { if (i) kputc(':', s); @@ -234,19 +265,37 @@ void bcf_fmt_core(const bcf_hdr_t *h, bcf1_t *b, kstring_t *s) if (k > 0) kputc(',', s); kputw(d[k], s); } - } else if (b->gi[i].fmt == bcf_str2int("DP", 2)) { + } else if (b->gi[i].fmt == bcf_str2int("DP", 2) || b->gi[i].fmt == bcf_str2int("DV", 2)) { kputw(((uint16_t*)b->gi[i].data)[j], s); } else if (b->gi[i].fmt == bcf_str2int("GQ", 2)) { kputw(((uint8_t*)b->gi[i].data)[j], s); + } else if (b->gi[i].fmt == bcf_str2int("SP", 2)) { + kputw(((int32_t*)b->gi[i].data)[j], s); } else if (b->gi[i].fmt == bcf_str2int("GT", 2)) { - int y = ((uint8_t*)b->gi[i].data)[j]; - if (y>>7&1) { - kputsn("./.", 3, s); - } else { - kputc('0' + (y>>3&7), s); - kputc("/|"[y>>6&1], s); - kputc('0' + (y&7), s); - } + int y = ((uint8_t*)b->gi[i].data)[j]; + if ( y>>7&1 ) + kputsn("./.", 3, s); + else if ( imax==-1 ) + { + kputc('0' + (y>>3&7), s); + kputc("/|"[y>>6&1], s); + kputc('0' + (y&7), s); + } + else + { + // Arguably, the while loop will be faster than two sqrts + int n = 0; + int row = 1; + while ( n>6&1], s); + kputw(row, s); + } } else if (b->gi[i].fmt == bcf_str2int("GL", 2)) { float *d = (float*)b->gi[i].data + j * x; int k; @@ -255,7 +304,7 @@ void bcf_fmt_core(const bcf_hdr_t *h, bcf1_t *b, kstring_t *s) if (k > 0) kputc(',', s); ksprintf(s, "%.2f", d[k]); } - } + } else kputc('.', s); // custom fields } } } @@ -267,3 +316,50 @@ char *bcf_fmt(const bcf_hdr_t *h, bcf1_t *b) bcf_fmt_core(h, b, &s); return s.s; } + +int bcf_append_info(bcf1_t *b, const char *info, int l) +{ + int shift = b->fmt - b->str; + int l_fmt = b->l_str - shift; + char *ori = b->str; + if (b->l_str + l > b->m_str) { // enlarge if necessary + b->m_str = b->l_str + l; + kroundup32(b->m_str); + b->str = realloc(b->str, b->m_str); + } + memmove(b->str + shift + l, b->str + shift, l_fmt); // move the FORMAT field + memcpy(b->str + shift - 1, info, l); // append to the INFO field + b->str[shift + l - 1] = '\0'; + b->fmt = b->str + shift + l; + b->l_str += l; + if (ori != b->str) bcf_sync(b); // synchronize when realloc changes the pointer + return 0; +} + +int bcf_cpy(bcf1_t *r, const bcf1_t *b) +{ + char *t1 = r->str; + bcf_ginfo_t *t2 = r->gi; + int i, t3 = r->m_str, t4 = r->m_gi; + *r = *b; + r->str = t1; r->gi = t2; r->m_str = t3; r->m_gi = t4; + if (r->m_str < b->m_str) { + r->m_str = b->m_str; + r->str = realloc(r->str, r->m_str); + } + memcpy(r->str, b->str, r->m_str); + bcf_sync(r); // calling bcf_sync() is simple but inefficient + for (i = 0; i < r->n_gi; ++i) + memcpy(r->gi[i].data, b->gi[i].data, r->n_smpl * r->gi[i].len); + return 0; +} + +int bcf_is_indel(const bcf1_t *b) +{ + char *p; + if (strlen(b->ref) > 1) return 1; + for (p = b->alt; *p; ++p) + if (*p != ',' && p[1] != ',' && p[1] != '\0') + return 1; + return 0; +}