Fixed the VCF header to pass validation

[samtools.git] / bcftools / call1.c

diff --git a/bcftools/call1.c b/bcftools/call1.c
index b913c9a4990e1f59bbdde8667d6496f842c02d21..2308607ab6d9de176570282f32a4797a275b6f7c 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -26,6 +26,7 @@ KSTREAM_INIT(gzFile, gzread, 16384)
 #define VC_CALL_GT 2048
 #define VC_ADJLD   4096
 #define VC_NO_INDEL 8192
+#define VC_ANNO_MAX 16384
 
 typedef struct {
        int flag, prior_type, n1;
@@ -200,6 +201,42 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p
        return is_var;
 }
 
+static void write_header(bcf_hdr_t *h)
+{
+       kstring_t str;
+       str.l = h->l_txt? h->l_txt - 1 : 0;
+       str.m = str.l + 1; str.s = h->txt;
+       if (!strstr(str.s, "##INFO=<ID=DP,"))
+               kputs("##INFO=<ID=DP,Number=1,Type=Integer,Description=\"Raw read depth\">\n", &str);
+       if (!strstr(str.s, "##INFO=<ID=DP4,"))
+               kputs("##INFO=<ID=DP4,Number=4,Type=Integer,Description=\"# high-quality ref-forward bases, ref-reverse, alt-forward and alt-reverse bases\">\n", &str);
+       if (!strstr(str.s, "##INFO=<ID=MQ,"))
+               kputs("##INFO=<ID=MQ,Number=1,Type=Integer,Description=\"Root-mean-square mapping quality of covering reads\">\n", &str);
+       if (!strstr(str.s, "##INFO=<ID=FQ,"))
+               kputs("##INFO=<ID=FQ,Number=1,Type=Float,Description=\"Phred probability that sample chromosomes are not all the same\">\n", &str);
+       if (!strstr(str.s, "##INFO=<ID=AF1,"))
+               kputs("##INFO=<ID=AF1,Number=1,Type=Float,Description=\"Max-likelihood estimate of the site allele frequency of the first ALT allele\">\n", &str);
+       if (!strstr(str.s, "##INFO=<ID=CI95,"))
+               kputs("##INFO=<ID=CI95,Number=2,Type=Float,Description=\"Equal-tail Bayesian credible interval of the site allele frequency at the 95% level\">\n", &str);
+       if (!strstr(str.s, "##INFO=<ID=PV4,"))
+               kputs("##INFO=<ID=PV4,Number=4,Type=Float,Description=\"P-values for strand bias, baseQ bias, mapQ bias and tail distance bias\">\n", &str);
+    if (!strstr(str.s, "##INFO=<ID=INDEL,"))
+        kputs("##INFO=<ID=INDEL,Number=0,Type=Flag,Descriptin=\"Indicates that the variant is an INDEL.\">\n", &str);
+    if (!strstr(str.s, "##INFO=<ID=GT,"))
+        kputs("##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n", &str);
+    if (!strstr(str.s, "##FORMAT=<ID=GQ,"))
+        kputs("##FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\"Genotype Quality\">\n", &str);
+    if (!strstr(str.s, "##FORMAT=<ID=GL,"))
+        kputs("##FORMAT=<ID=GL,Number=3,Type=Float,Description=\"Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)\">\n", &str);
+       if (!strstr(str.s, "##FORMAT=<ID=DP,"))
+               kputs("##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"# high-quality bases\">\n", &str);
+       if (!strstr(str.s, "##FORMAT=<ID=SP,"))
+               kputs("##FORMAT=<ID=SP,Number=1,Type=Integer,Description=\"Phred-scaled strand bias P-value\">\n", &str);
+       if (!strstr(str.s, "##FORMAT=<ID=PL,"))
+               kputs("##FORMAT=<ID=PL,Number=-1,Type=Integer,Description=\"List of Phred-scaled genotype likelihoods, number of values is (#ALT+1)*(#ALT+2)/2\">\n", &str);
+       h->l_txt = str.l + 1; h->txt = str.s;
+}
+
 double bcf_ld_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4]);
 
 int bcfview(int argc, char *argv[])
@@ -207,6 +244,7 @@ int bcfview(int argc, char *argv[])
        extern int bcf_2qcall(bcf_hdr_t *h, bcf1_t *b);
        extern void bcf_p1_indel_prior(bcf_p1aux_t *ma, double x);
        extern int bcf_fix_gt(bcf1_t *b);
+       extern int bcf_anno_max(bcf1_t *b);
        bcf_t *bp, *bout = 0;
        bcf1_t *b, *blast;
        int c;
@@ -221,7 +259,7 @@ int bcfview(int argc, char *argv[])
        tid = begin = end = -1;
        memset(&vc, 0, sizeof(viewconf_t));
        vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; vc.indel_frac = -1.;
-       while ((c = getopt(argc, argv, "N1:l:cHAGvbSuP:t:p:QgLi:I")) >= 0) {
+       while ((c = getopt(argc, argv, "N1:l:cHAGvbSuP:t:p:QgLi:IM")) >= 0) {
                switch (c) {
                case '1': vc.n1 = atoi(optarg); break;
                case 'l': vc.fn_list = strdup(optarg); break;
@@ -236,6 +274,7 @@ int bcfview(int argc, char *argv[])
                case 'H': vc.flag |= VC_HWE; break;
                case 'g': vc.flag |= VC_CALL_GT | VC_CALL; break;
                case 'I': vc.flag |= VC_NO_INDEL; break;
+               case 'M': vc.flag |= VC_ANNO_MAX; break;
                case 't': vc.theta = atof(optarg); break;
                case 'p': vc.pref = atof(optarg); break;
                case 'i': vc.indel_frac = atof(optarg); break;
@@ -285,7 +324,10 @@ int bcfview(int argc, char *argv[])
        bp = vcf_open(argv[optind], moder);
        h = vcf_hdr_read(bp);
        bout = vcf_open("-", modew);
-       if (!(vc.flag & VC_QCALL)) vcf_hdr_write(bout, h);
+       if (!(vc.flag & VC_QCALL)) {
+               if (vc.flag & VC_CALL) write_header(h);
+               vcf_hdr_write(bout, h);
+       }
        if (vc.flag & VC_CALL) {
                p1 = bcf_p1_init(h->n_smpl);
                if (vc.prior_file) {
@@ -369,11 +411,12 @@ int bcfview(int argc, char *argv[])
                        }
                        bcf_cpy(blast, b);
                }
+               if (vc.flag & VC_ANNO_MAX) bcf_anno_max(b);
                if (vc.flag & VC_NO_GENO) { // do not output GENO fields
                        b->n_gi = 0;
                        b->fmt[0] = '\0';
-               }
-               bcf_fix_gt(b);
+                       b->l_str = b->fmt - b->str + 1;
+               } else bcf_fix_gt(b);
                vcf_write(bout, h, b);
        }
        if (vc.prior_file) free(vc.prior_file);