X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fbcf.c;h=24728dbe1515e7cc1171a8cc08f842bafa167d5e;hb=0ccd9d36ebf35ce620a8248ecf4336c84065e6c0;hp=19969eb5d4a4f15ad0e4e9f2486123794e8b7410;hpb=6418571502ba9c6bd4eb930c7ae0309f7f251b4e;p=samtools.git diff --git a/bcftools/bcf.c b/bcftools/bcf.c index 19969eb..24728db 100644 --- a/bcftools/bcf.c +++ b/bcftools/bcf.c @@ -4,12 +4,6 @@ #include "kstring.h" #include "bcf.h" -void bcf_hdr_clear(bcf_hdr_t *b) -{ - free(b->name); free(b->sname); free(b->txt); free(b->ns); free(b->sns); - memset(b, 0, sizeof(bcf_hdr_t)); -} - bcf_t *bcf_open(const char *fn, const char *mode) { bcf_t *b; @@ -19,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; } @@ -100,17 +93,23 @@ int bcf_hdr_sync(bcf_hdr_t *b) return 0; } -int bcf_sync(int n_smpl, bcf1_t *b) +int bcf_sync(bcf1_t *b) { char *p, *tmp[5]; - int i, n; + int i, n, n_smpl = b->n_smpl; 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, "[bcf_sync] incorrect number of fields (%d != 5)\n", 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]; @@ -138,12 +137,14 @@ int bcf_sync(int n_smpl, 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("GL", 2)) { + } 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; } b->gi[i].data = realloc(b->gi[i].data, n_smpl * b->gi[i].len); } @@ -172,6 +173,7 @@ int bcf_read(bcf_t *bp, const bcf_hdr_t *h, bcf1_t *b) int i, l = 0; if (b == 0) return -1; if (bgzf_read(bp->fp, &b->tid, 4) == 0) return -1; + b->n_smpl = h->n_smpl; bgzf_read(bp->fp, &b->pos, 4); bgzf_read(bp->fp, &b->qual, 4); bgzf_read(bp->fp, &b->l_str, 4); @@ -182,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(h->n_smpl, 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; @@ -195,7 +197,7 @@ int bcf_destroy(bcf1_t *b) int i; if (b == 0) return -1; free(b->str); - for (i = 0; i < b->n_gi; ++i) + for (i = 0; i < b->m_gi; ++i) free(b->gi[i].data); free(b->gi); free(b); @@ -219,7 +221,7 @@ void bcf_fmt_core(const bcf_hdr_t *h, bcf1_t *b, kstring_t *s) fmt_str(b->str, s); kputc('\t', s); fmt_str(b->ref, s); kputc('\t', s); fmt_str(b->alt, s); kputc('\t', s); - ksprintf(s, "%.3g", b->qual); /*kputw(b->qual, s);*/ kputc('\t', s); + ksprintf(s, "%.3g", b->qual); kputc('\t', s); fmt_str(b->flt, s); kputc('\t', s); fmt_str(b->info, s); if (b->fmt[0]) { @@ -228,27 +230,68 @@ 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) { + int ploidy = b->ploidy ? b->ploidy[j] : 2; kputc('\t', s); for (i = 0; i < b->n_gi; ++i) { if (i) kputc(':', s); if (b->gi[i].fmt == bcf_str2int("PL", 2)) { uint8_t *d = (uint8_t*)b->gi[i].data + j * x; int k; - for (k = 0; k < x; ++k) { - if (k > 0) kputc(',', s); - kputw(d[k], s); - } - } else if (b->gi[i].fmt == bcf_str2int("DP", 2)) { + if ( ploidy==1 ) + for (k=0; kn_alleles; k++) + { + if (k>0) kputc(',', s); + kputw(d[(k+1)*(k+2)/2-1], s); + } + else + for (k = 0; k < x; ++k) { + if (k > 0) kputc(',', s); + kputw(d[k], s); + } + } 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]; - 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 ( ploidy==1 ) + { + if ( y>>7&1 ) + kputc('.', s); + else + kputc('0' + (y>>3&7), s); + } + else + { + if ( y>>7&1 ) + kputsn("./.", 3, s); + else { + kputc('0' + (y>>3&7), s); + kputc("/|"[y>>6&1], s); + kputc('0' + (y&7), s); + } + } + } else if (b->gi[i].fmt == bcf_str2int("GL", 2)) { + float *d = (float*)b->gi[i].data + j * x; + int k; + //printf("- %lx\n", d); + for (k = 0; k < x; ++k) { + if (k > 0) kputc(',', s); + ksprintf(s, "%.2f", d[k]); + } + } else kputc('.', s); // custom fields } } } @@ -260,3 +303,94 @@ 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 remove_tag(char *str, const char *tag, char delim) +{ + char *tmp = str, *p; + int len_diff = 0, ori_len = strlen(str); + while ( *tmp && (p = strstr(tmp,tag)) ) + { + if ( p>str ) + { + if ( *(p-1)!=delim ) { tmp=p+1; continue; } // shared substring + p--; + } + char *q=p+1; + while ( *q && *q!=delim ) q++; + if ( p==str && *q ) q++; // the tag is first, don't move the delim char + len_diff += q-p; + if ( ! *q ) { *p = 0; break; } // the tag was last, no delim follows + else + memmove(p,q,ori_len-(int)(p-str)-(int)(q-p)); // *q==delim + } + if ( len_diff==ori_len ) + str[0]='.', str[1]=0, len_diff--; + + return len_diff; +} + + +void rm_info(kstring_t *s, const char *key) +{ + char *p = s->s; + int n = 0; + while ( n<4 ) + { + if ( !*p ) n++; + p++; + } + char *q = p+1; + while ( *q && q-s->sl ) q++; + + int nrm = remove_tag(p, key, ';'); + if ( nrm ) + memmove(q-nrm, q, s->s+s->l-q+1); + s->l -= nrm; +} + +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; +}