bugfix: wrong SP; missing DV in the VCF header

diff --git a/bam.h b/bam.h
index df345b6e47afcbcbdf545e9f954277b78825272b..11f8028cf6d96ece1390074c6fb155437ce5f766 100644 (file)
--- a/bam.h
+++ b/bam.h
@@ -40,7 +40,7 @@
   @copyright Genome Research Ltd.
  */
 
-#define BAM_VERSION "0.1.18-dev (r982:310)"
+#define BAM_VERSION "0.1.18-dev (r982:313)"
 
 #include <stdint.h>
 #include <stdlib.h>

diff --git a/bam2bcf.c b/bam2bcf.c
index 83f2e564b67b85c711b0a0b83d1e29fc52591eb5..6ac5dce182a2353ada163dd65bdb5ca530d59127 100644 (file)
--- a/bam2bcf.c
+++ b/bam2bcf.c
@@ -328,7 +328,7 @@ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bc
        if (bcr && fmt_flag) {
                uint16_t *dp = (fmt_flag & B2B_FMT_DP)? b->gi[1].data : 0;
                uint16_t *dv = (fmt_flag & B2B_FMT_DV)? b->gi[1 + ((fmt_flag & B2B_FMT_DP) != 0)].data : 0;
-               int32_t  *sp = (fmt_flag & B2B_FMT_DV)? b->gi[1 + ((fmt_flag & B2B_FMT_DP) != 0) + ((fmt_flag & B2B_FMT_DV) != 0)].data : 0;
+               int32_t  *sp = (fmt_flag & B2B_FMT_SP)? b->gi[1 + ((fmt_flag & B2B_FMT_DP) != 0) + ((fmt_flag & B2B_FMT_DV) != 0)].data : 0;
                for (i = 0; i < bc->n; ++i) {
                        bcf_callret1_t *p = bcr + i;
                        if (dp) dp[i] = p->depth  < 0xffff? p->depth  : 0xffff;

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 3cc464949eb69e630ee7905c96d4f8050b1b569c..f6717b37297216f4aba010611b7c72bda90892fa 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -33,6 +33,7 @@ KSTREAM_INIT(gzFile, gzread, 16384)
 #define VC_EM       0x10000
 #define VC_PAIRCALL 0x20000
 #define VC_QCNT     0x40000
+#define VC_INDEL_ONLY 0x80000
 
 typedef struct {
        int flag, prior_type, n1, n_sub, *sublist, n_perm;
@@ -283,6 +284,8 @@ static void write_header(bcf_hdr_t *h)
         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=DV,"))
+               kputs("##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"# high-quality non-reference 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,"))
@@ -318,7 +321,7 @@ int bcfview(int argc, char *argv[])
        memset(&vc, 0, sizeof(viewconf_t));
        vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; vc.indel_frac = -1.; vc.n_perm = 0; vc.min_perm_p = 0.01; vc.min_smpl_frac = 0; vc.min_lrt = 1;
        memset(qcnt, 0, 8 * 256);
-       while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:T:Y")) >= 0) {
+       while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:T:Yw")) >= 0) {
                switch (c) {
                case '1': vc.n1 = atoi(optarg); break;
                case 'l': vc.bed = bed_read(optarg); break;
@@ -335,6 +338,7 @@ int bcfview(int argc, char *argv[])
                case 'u': vc.flag |= VC_UNCOMP | VC_BCFOUT; break;
                case 'g': vc.flag |= VC_CALL_GT | VC_CALL; break;
                case 'I': vc.flag |= VC_NO_INDEL; break;
+               case 'w': vc.flag |= VC_INDEL_ONLY; break;
                case 'M': vc.flag |= VC_ANNO_MAX; break;
                case 'Y': vc.flag |= VC_QCNT; break;
                case 't': vc.theta = atof(optarg); break;
@@ -482,6 +486,7 @@ int bcfview(int argc, char *argv[])
                if (vc.flag & VC_FIX_PL) bcf_fix_pl(b);
                is_indel = bcf_is_indel(b);
                if ((vc.flag & VC_NO_INDEL) && is_indel) continue;
+               if ((vc.flag & VC_INDEL_ONLY) && !is_indel) continue;
                if ((vc.flag & VC_ACGT_ONLY) && !is_indel) {
                        int x;
                        if (b->ref[0] == 0 || b->ref[1] != 0) continue;