* vcfutils.pl: fixed a typo in help message

[samtools.git] / bcftools / bcf.c

diff --git a/bcftools/bcf.c b/bcftools/bcf.c
index e375a1fb8791cd6b3f34cfa070d53396b158133e..3531647e18775cd5049954385a1022d39a17ec6f 100644 (file)
--- 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;
@@ -100,21 +94,27 @@ 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;
-       if (n != 5) return -1;
+       if (n != 5) {
+               fprintf(stderr, "[%s] incorrect number of fields (%d != 5). Corrupted file?\n", __func__, n);
+               return -1;
+       }
        b->ref = tmp[0]; b->alt = tmp[1]; b->flt = tmp[2]; b->info = tmp[3]; b->fmt = tmp[4];
        // set n_alleles
-       for (p = b->alt, n = 1; *p; ++p)
-               if (*p == ',') ++n;
-       b->n_alleles = n + 1;
+       if (*b->alt == 0) b->n_alleles = 1;
+       else {
+               for (p = b->alt, n = 1; *p; ++p)
+                       if (*p == ',') ++n;
+               b->n_alleles = n + 1;
+       }
        // set n_gi and gi[i].fmt
        for (p = b->fmt, n = 1; *p; ++p)
                if (*p == ':') ++n;
@@ -137,7 +137,7 @@ int bcf_sync(int n_smpl, bcf1_t *b)
                } 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)) {
-                       b->gi[i].len = 4;
+                       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);
        }
@@ -166,6 +166,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);
@@ -176,7 +177,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);
+       bcf_sync(b);
        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;
@@ -189,7 +190,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);
@@ -213,7 +214,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);
-       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]) {
@@ -239,9 +240,21 @@ void bcf_fmt_core(const bcf_hdr_t *h, bcf1_t *b, kstring_t *s)
                                kputw(((uint8_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);
+                               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]);
+                               }
                        }
                }
        }
@@ -254,3 +267,41 @@ 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;
+       bcf_sync(b);
+//     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;
+}