SNP calling from the GL field

diff --git a/bcftools/bcf.h b/bcftools/bcf.h
index e6ba2252282f52af80bcdcc89aab92c716563a8a..d0f9b02bb5b554df977c2202902605c0642cf223 100644 (file)
--- a/bcftools/bcf.h
+++ b/bcftools/bcf.h
@@ -66,6 +66,8 @@ extern "C" {
        int vcf_read(bcf_t *bp, bcf_hdr_t *h, bcf1_t *b);
 
        int bcf_shrink_alt(int n_smpl, bcf1_t *b, int n);
        int vcf_read(bcf_t *bp, bcf_hdr_t *h, bcf1_t *b);
 
        int bcf_shrink_alt(int n_smpl, bcf1_t *b, int n);

+       int bcf_gl2pl(int n_smpl, bcf1_t *b);
+

        void *bcf_build_refhash(bcf_hdr_t *h);
        void bcf_str2id_destroy(void *_hash);
        int bcf_str2id_add(void *_hash, const char *str);
        void *bcf_build_refhash(bcf_hdr_t *h);
        void bcf_str2id_destroy(void *_hash);
        int bcf_str2id_add(void *_hash, const char *str);

diff --git a/bcftools/bcfutils.c b/bcftools/bcfutils.c
index 6908bc478bf7081409f9ef9f57b332ef9031ce39..dd1e89f16cf1c48ac8f90ceae35bae312d76ba54 100644 (file)
--- a/bcftools/bcfutils.c
+++ b/bcftools/bcfutils.c
@@ -82,3 +82,29 @@ int bcf_shrink_alt(int n_smpl, bcf1_t *b, int n)
        bcf_sync(n_smpl, b);
        return 0;
 }
        bcf_sync(n_smpl, b);
        return 0;
 }

+
+int bcf_gl2pl(int n_smpl, bcf1_t *b)
+{
+       char *p;
+       int i;
+       bcf_ginfo_t *g;
+       float *d0;
+       uint8_t *d1;
+       if (strstr(b->fmt, "PL")) return -1;
+       if ((p = strstr(b->fmt, "GL")) == 0) return -1;
+       *p = 'P';
+       for (i = 0; i < b->n_gi; ++i)
+               if (b->gi[i].fmt == bcf_str2int("GL", 2))
+                       break;
+       g = b->gi + i;
+       g->fmt = bcf_str2int("PL", 2);
+       g->len /= 4; // 4 == sizeof(float)
+       d0 = (float*)g->data; d1 = (uint8_t*)g->data;
+       for (i = 0; i < n_smpl * g->len; ++i) {
+               int x = (int)(-10. * d0[i] + .499);
+               if (x > 255) x = 255;
+               if (x < 0) x = 0;
+               d1[i] = x;
+       }
+       return 0;
+}

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 1060513cfb104072c271aee4646905351b06f99b..995d500eb69ad028a1b0ed2a75da5bde88d92bee 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -287,6 +287,7 @@ int bcfview(int argc, char *argv[])
                }
                if (vc.flag & VC_CALL) {
                        bcf_p1rst_t pr;
                }
                if (vc.flag & VC_CALL) {
                        bcf_p1rst_t pr;

+                       bcf_gl2pl(h->n_smpl, b);

                        bcf_p1_cal(b, p1, &pr); // pr.g[3] is not calculated here
                        if (vc.flag&VC_HWE) bcf_p1_cal_g3(p1, pr.g);
                        if ((n_processed + 1) % 50000 == 0) bcf_p1_dump_afs(p1);
                        bcf_p1_cal(b, p1, &pr); // pr.g[3] is not calculated here
                        if (vc.flag&VC_HWE) bcf_p1_cal_g3(p1, pr.g);
                        if ((n_processed + 1) % 50000 == 0) bcf_p1_dump_afs(p1);