Bad floating point parsing in AF field

[jackgoldsmith4]

I created a genomicsDB array using the VCF importer and this VCF. I then used a GenomicsDBFeatureReader to iterate over the VariantContext objects generated from that array, but the numbers I got for the AF INFO field (type = "float") were sometimes slightly different than the original values from the VCF. For example, in the second VariantContext, the original value for the AF field is 0.007252, but the number in the VariantContext object is 0.00725199.

[kdatta]

This is an artifact of the print function used to write the VCF file which rounded the values to 7th decimal place. The VariantContext object in memory contains the actual value i.e. 0.00725199

[kgururaj]

This is always going to happen - FP representation in ASCII is fuzzy. An "equality" comparator for FP numbers should be a "nearness" comparator with a tolerance threshold.

[danking]

Floating point numbers can be represented precisely in decimal notation. IEEE 754 defines the size of a 32-bit floating point mantissa as 24-bits. log_10(2^24) is approximately 7.22, ergo, you need at least 8 decimal digits to distinguish between 0.00725199 and 0.00725200.

A deterministic, lossless format should preserve equality of floating-point numbers.

[kgururaj]

Yep, the binary IEEE representation of FP numbers preserves the exact value (for the most part). It's the conversion to ASCII base 10 representation that leads to these weird issues. Each software system does this slightly differently.

If you use the BCF2Codec, such discrepancies might reduce.

[danking]

To be very clear, the two character strings in question 0.007252 and 0.00725199 are not parsed by sscanf as the same number, they're parsed as two distinct floating point numbers. These two numbers are separated from one another by 20 other floating point numbers.

0.007252f    : 0 01110111 11011011010001000101111
0.00725199   : 0 01110111 11011011010001000011010

Note the last six bits. Here's a few of the intermediate numbers:

0.0072519904 : 0 01110111 11011011010001000011011
0.007251991  : 0 01110111 11011011010001000011100
0.0072519914 : 0 01110111 11011011010001000011101

All of these numbers have the same decimal representation when rounded to eight digits after the decimal point, namely: 0.00725199.

However, all is not lost. We can reliably interchange these numbers by using eight significant digits. This is quite easy in C! Anywhere genomicsdb prints a floating point number, use the format specifier %0.8e and Hail will be able to load exactly the value that genomicsdb is using internally.

This leaves open the question of why the character string 7.252e-03 (from the VCF file) is eventually represented in genomicsdb as a floating-point number with a printed representation of 0.00725199.

[kgururaj]

GenomicsDB uses htslib to print out the values in the VCF format. The Java API uses htsjdk VCFCodec to parse and print the values. The htslib library uses %g to print out FP values - 6 significant decimals.

It would probably be better to print 9 significant decimals for single precision FP. However, once this value is consumed by downstream software (Java libs etc), we have no control.

FYI I could not reproduce this issue on our systems.

[danking]

Hey @kgururaj, I created a simple example that exhibits the behavior. If you clone that to /private/tmp/ and rename the cloned folder to foo (so that the gdb workspace is located at /private/tmp/foo/bar), then change directories to /private/tmp/foo/gdb-fp and execute sbt run, you should get the same output I've pasted below.

The floating-point number 0.007252 should not get printed as 0.00725199 by any system. I'm not sure at which step the value of the AF field was corrupted. AFAIK, these are the only possible points at which it could have been corrupted:

1. vcf2tiledb loader-config.json
2. the C/C++ code that connects GenomicsDB to HTSJDK
3. the Java API of GenomicsDB
4. the HTSJDK VCFCodec

>> comparing output from gdb to vcf <<
Class of the AF field in the attributes map: java.lang.String
af.toDouble, then printed as a String by Scala: 0.00725199

the line in question:
22	16050115	.	G	A	134.13	.	AC=20;AF=7.25200e-03;AN=2758;BaseQRankSum=5.04300e+00;ClippingRankSum=9.60000e-02;DP=113;FS=0.00000e+00;InbreedingCoeff=-7.70000e-02;MLEAC=14;MLEAF=5.07600e-03;MQ=2.99300e+01;MQ0=0;MQRankSum=3.18000e-01;QD=4.63000e+00;ReadPosRankSum=2.99000e-01;	GT:AD:DP:GQ:PL	./.	0/0:1,0:1:3:0,3,22	./.	0/0:1,0:1:3:0,3,22	0/0:1,0:1:3:0,3,22	0/0:3,0:3:9:0,9,84	./.	0/0:2,0:2:6:0,6,52	0/0:3,0:3:9:0,9,80	0/0:6,0:6:18:0,18,154	./.	./.	./.	./.	0/0:2,0:2:6:0,6,50	./.	0/0:3,0:3:9:0,9,81	0/0:2,0:2:6:0,6,60	0/0:2,0:2:6:0,6,49	0/0:3,0:3:9:0,9,69	0/0:1,0:1:3:0,3,35	0/0:1,0:1:3:0,3,22	./.	0/0:1,0:1:3:0,3,29	0/0:1,0:1:3:0,3,29	0/0:4,0:4:12:0,12,96	0/0:8,0:8:24:0,24,216	./.	0/0:3,0:3:9:0,9,97	./.	0/0:1,0:1:3:0,3,23	./.	./.	./.	0/0:4,0:4:12:0,12,114	0/0:1,0:1:3:0,3,22	./.	./.	./.	./.	0/0:3,0:3:9:0,9,68	0/0:3,0:3:9:0,9,81	0/0:2,0:2:6:0,6,58	./.	0/0:3,0:3:9:0,9,82	./.	0/0:6,0:6:18:0,18,160	0/0:2,0:2:6:0,6,71	0/0:2,0:2:6:0,6,55	./.	./.	./.	./.	./.	./.	0/0:3,0:3:9:0,9,76	0/0:1,0:1:3:0,3,28	./.	./.	0/0:3,0:3:9:0,9,79	./.	0/0:1,0:1:3:0,3,19	./.	./.	0/0:1,0:1:3:0,3,31	./.	./.	0/0:1,0:1:3:0,3,22	0/0:3,0:3:9:0,9,81	0/0:2,0:2:6:0,6,51	0/0:1,0:1:3:0,3,30	./.	0/0:1,0:1:3:0,3,29	0/0:1,0:1:3:0,3,22	./.	0/0:1,0:1:3:0,3,20	./.	./.	0/0:2,0:2:6:0,6,52	./.	./.	0/0:5,0:5:15:0,15,146	0/0:1,0:1:3:0,3,27	0/0:1,0:1:3:0,3,26	./.	0/0:3,0:3:9:0,9,86	0/0:1,0:1:3:0,3,25	./.	0/0:2,0:2:6:0,6,45	./.	0/0:1,0:1:3:0,3,22	0/0:1,0:1:3:0,3,22	./.	0/0:1,0:1:3:0,3,22	./.	0/0:1,0:1:3:0,3,22	./.	./.	./.	./.
original numeric text: 7.25200e-03
original numeric text parsed by Scala, then printed as a String by Scala: 0.007252

fromGDB == fromVCF?  false
fromGDB binary: 11111101111101101101000100001100111110000110000011111010010000
fromVCF binary: 11111101111101101101000100010111101101010010100001101011010110


>> now again with floats <<
af.toFloat, then printed as a String by Scala: 0.00725199
original numeric text parsed by Scala, then printed as a String by Scala: 0.007252
fromGDB == fromVCF?  false
fromGDB binary: 111011111011011010001000011010
fromVCF binary: 111011111011011010001000101111

[danking]

I'm not sure exactly whether it should be /private/tmp or /tmp. On OS X /tmp is a symlink to /private/tmp. My $CWD was /tmp when I created and populated the instance.