Merge Don't encode negative ints as reserved values (PR #766)

Add NEWS item describing the bug.

Author: daviesrob

NEWS

@@ -1,6 +1,16 @@
 Noteworthy changes in release a.b
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+* Fixed bug where some 8 or 16-bit negative integers were stored using values
+  reserved by the BCF specification.  These numbers are now promoted to the
+  next size up, so -120 to -128 are stored using at least 16 bits, and -32760
+  to -32768 are stored using 32 bits.
+
+  Note that while BCF files affected by this bug are technically incorrect,
+  it is still possible to read them.  When converting to VCF format,
+  HTSlib (and therefore bcftools) will interpret the values as intended
+  and write out the correct negative numbers. (#766, thanks to John Marshall)
+
 Noteworthy changes in release 1.9 (18th July 2018)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

htslib/vcf.h

@@ -953,8 +953,8 @@ static inline void bcf_enc_size(kstring_t *s, int size, int type)
 
 static inline int bcf_enc_inttype(long x)
 {
-    if (x <= INT8_MAX && x > bcf_int8_missing) return BCF_BT_INT8;
-    if (x <= INT16_MAX && x > bcf_int16_missing) return BCF_BT_INT16;
+    if (x <= INT8_MAX && x >= INT8_MIN + 8) return BCF_BT_INT8;
+    if (x <= INT16_MAX && x >= INT16_MIN + 8) return BCF_BT_INT16;
     return BCF_BT_INT32;
 }
 
@@ -966,10 +966,10 @@ static inline void bcf_enc_int1(kstring_t *s, int32_t x)
     } else if (x == bcf_int32_missing) {
         bcf_enc_size(s, 1, BCF_BT_INT8);
         kputc(bcf_int8_missing, s);
-    } else if (x <= INT8_MAX && x > bcf_int8_missing) {
+    } else if (x <= INT8_MAX && x >= INT8_MIN + 8) {
         bcf_enc_size(s, 1, BCF_BT_INT8);
         kputc(x, s);
-    } else if (x <= INT16_MAX && x > bcf_int16_missing) {
+    } else if (x <= INT16_MAX && x >= INT16_MIN + 8) {
         int16_t z = x;
         bcf_enc_size(s, 1, BCF_BT_INT16);
         kputsn((char*)&z, 2, s);

test/test-vcf-api.c

@@ -62,6 +62,7 @@ void write_bcf(char *fname)
     bcf_hdr_append(hdr, "##phasing=partial");
     bcf_hdr_append(hdr, "##INFO=<ID=NS,Number=1,Type=Integer,Description=\"Number of Samples With Data\">");
     bcf_hdr_append(hdr, "##INFO=<ID=DP,Number=1,Type=Integer,Description=\"Total Depth\">");
+    bcf_hdr_append(hdr, "##INFO=<ID=NEG,Number=.,Type=Integer,Description=\"Test Negative Numbers\">");
     bcf_hdr_append(hdr, "##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency\">");
     bcf_hdr_append(hdr, "##INFO=<ID=AA,Number=1,Type=String,Description=\"Ancestral Allele\">");
     bcf_hdr_append(hdr, "##INFO=<ID=DB,Number=0,Type=Flag,Description=\"dbSNP membership, build 129\">");
@@ -82,7 +83,7 @@ void write_bcf(char *fname)
 
 
     // Add a record
-    // 20     14370   rs6054257 G      A       29   PASS   NS=3;DP=14;AF=0.5;DB;H2           GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
+    // 20     14370   rs6054257 G      A       29   PASS   NS=3;DP=14;NEG=-127;AF=0.5;DB;H2           GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
     // .. CHROM
     rec->rid = bcf_hdr_name2id(hdr, "20");
     // .. POS
@@ -101,6 +102,8 @@ void write_bcf(char *fname)
     bcf_update_info_int32(hdr, rec, "NS", &tmpi, 1);
     tmpi = 14;
     bcf_update_info_int32(hdr, rec, "DP", &tmpi, 1);
+    tmpi = -127;
+    bcf_update_info_int32(hdr, rec, "NEG", &tmpi, 1);
     float tmpf = 0.5;
     bcf_update_info_float(hdr, rec, "AF", &tmpf, 1);
     bcf_update_info_flag(hdr, rec, "DB", NULL, 1);
@@ -133,7 +136,7 @@ void write_bcf(char *fname)
     bcf_update_format_string(hdr, rec, "TS", (const char**)tmp_str, 3);
     if ( bcf_write1(fp, hdr, rec)!=0 ) error("Failed to write to %s\n", fname);
 
-    // 20     1110696 . A      G,T     67   .   NS=2;DP=10;AF=0.333,.;AA=T;DB GT 2 1   ./.
+    // 20     1110696 . A      G,T     67   .   NS=2;DP=10;NEG=-128;AF=0.333,.;AA=T;DB GT 2 1   ./.
     bcf_clear1(rec);
     rec->rid = bcf_hdr_name2id(hdr, "20");
     rec->pos = 1110695;
@@ -143,6 +146,8 @@ void write_bcf(char *fname)
     bcf_update_info_int32(hdr, rec, "NS", &tmpi, 1);
     tmpi = 10;
     bcf_update_info_int32(hdr, rec, "DP", &tmpi, 1);
+    tmpi = -128;
+    bcf_update_info_int32(hdr, rec, "NEG", &tmpi, 1);
     float *tmpfa = (float*)malloc(2*sizeof(float));
     tmpfa[0] = 0.333;
     bcf_float_set_missing(tmpfa[1]);
@@ -291,6 +296,7 @@ void test_get_info_values(const char *fname)
     while (bcf_read(fp, hdr, line) == 0)
     {
         float *afs = 0;
+        int32_t *negs = NULL;
         int count = 0;
         int ret = bcf_get_info_float(hdr, line, "AF", &afs, &count);
 
@@ -312,6 +318,21 @@ void test_get_info_values(const char *fname)
         }
 
         free(afs);
+
+        int32_t expected = (line->pos == 14369)? -127 : -128;
+        count = 0;
+        ret = bcf_get_info_int32(hdr, line, "NEG", &negs, &count);
+        if (ret != 1 || negs[0] != expected)
+        {
+            if (ret < 0)
+                fprintf(stderr, "NEG should be %d, got error ret=%d\n", expected, ret);
+            else if (ret == 0)
+                fprintf(stderr, "NEG should be %d, got no entries\n", expected);
+            else
+                fprintf(stderr, "NEG should be %d, got %d entries (first is %d)\n", expected, ret, negs[0]);
+            exit(1);
+        }
+        free(negs);
     }
 
     bcf_destroy(line);

test/test-vcf-api.out

@@ -8,6 +8,7 @@
 ##phasing=partial
 ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
 ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
+##INFO=<ID=NEG,Number=.,Type=Integer,Description="Test Negative Numbers">
 ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
 ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
 ##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
@@ -20,9 +21,9 @@
 ##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
 ##FORMAT=<ID=TS,Number=1,Type=String,Description="Test String">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA00001	NA00002	NA00003
-20	14370	rs6054257	G	A	29	PASS	NS=3;DP=14;AF=0.5;DB;H2	GT:GQ:DP:HQ:TS	0|0:48:1:51,51:String1	1|0:48:8:51,51:SomeOtherString2	1/1:43:5:.,.:YetAnotherString3
-20	14370	.	G	A	29	PASS	NS=3;DP=14;AF=0.5;DB;H2	GT:DP:HQ:TS	0|0:1:51,51:String1	1|0:8:51,51:SomeOtherString2	1/1:5:.,.:YetAnotherString3
-20	14370	.	G	A	29	PASS	NS=3;DP=99;AF=0.5;DB;H2	GT:DP:HQ:TS	0|0:9:51,51:String1	1|0:9:51,51:SomeOtherString2	1/1:9:.,.:YetAnotherString3
-20	1110696	.	A	G,T	67	.	NS=2;DP=10;AF=0.333,.;AA=T;DB	GT	2	1	./.
-20	1110696	.	A	G,T	67	.	NS=2;DP=10;AF=0.333,.;AA=T;DB	GT	2	1	./.
-20	1110696	.	G	A	67	.	NS=2;DP=99;AF=0.333,.;AA=T;DB	GT:DP	2:9	1:9	./.:9
+20	14370	rs6054257	G	A	29	PASS	NS=3;DP=14;NEG=-127;AF=0.5;DB;H2	GT:GQ:DP:HQ:TS	0|0:48:1:51,51:String1	1|0:48:8:51,51:SomeOtherString2	1/1:43:5:.,.:YetAnotherString3
+20	14370	.	G	A	29	PASS	NS=3;DP=14;NEG=-127;AF=0.5;DB;H2	GT:DP:HQ:TS	0|0:1:51,51:String1	1|0:8:51,51:SomeOtherString2	1/1:5:.,.:YetAnotherString3
+20	14370	.	G	A	29	PASS	NS=3;DP=99;NEG=-127;AF=0.5;DB;H2	GT:DP:HQ:TS	0|0:9:51,51:String1	1|0:9:51,51:SomeOtherString2	1/1:9:.,.:YetAnotherString3
+20	1110696	.	A	G,T	67	.	NS=2;DP=10;NEG=-128;AF=0.333,.;AA=T;DB	GT	2	1	./.
+20	1110696	.	A	G,T	67	.	NS=2;DP=10;NEG=-128;AF=0.333,.;AA=T;DB	GT	2	1	./.
+20	1110696	.	G	A	67	.	NS=2;DP=99;NEG=-128;AF=0.333,.;AA=T;DB	GT:DP	2:9	1:9	./.:9

vcf.c

@@ -1742,13 +1742,13 @@ void bcf_enc_vint(kstring_t *s, int n, int32_t *a, int wsize)
             if (max < a[i]) max = a[i];
             if (min > a[i]) min = a[i];
         }
-        if (max <= INT8_MAX && min > bcf_int8_vector_end) {
+        if (max <= INT8_MAX && min >= INT8_MIN + 8) {
             bcf_enc_size(s, wsize, BCF_BT_INT8);
             for (i = 0; i < n; ++i)
                 if ( a[i]==bcf_int32_vector_end ) kputc(bcf_int8_vector_end, s);
                 else if ( a[i]==bcf_int32_missing ) kputc(bcf_int8_missing, s);
                 else kputc(a[i], s);
-        } else if (max <= INT16_MAX && min > bcf_int16_vector_end) {
+        } else if (max <= INT16_MAX && min >= INT16_MIN + 8) {
             uint8_t *p;
             bcf_enc_size(s, wsize, BCF_BT_INT16);
             ks_resize(s, s->l + n * sizeof(int16_t));