Reading VCF files? about oxbow HOT 4 CLOSED

The noodles developer has previously mentioned a strong desire to maintain spec compliance, so in this case, your error was a feature and not a bug. We've found that many VCF files in the wild aren't spec compliant, so the situation isn't ideal regarding those files unfortunately.

from oxbow.

Hi Andreas, thank you for the detailed report. I was able to replicate your issue.

The error indicates that RS values in the INFO column are causing integer overflow.

To check this, I dug into the VCF file itself with pysam, another tool.

# pysam_check.py
import pysam
with pysam.VariantFile("./GCF_000001405.40.gz") as f:
    records = []
    for record in f.fetch():
        records.append(record)
    print(records)

This produced the following output and warning, which seems consistent with the error message you got.

NC_000001.11    6259533 rs2148352434    C       T       .       .       RS=.;dbSNPBuildID=156;SSR=0;GENEINFO=GPR153:387509;VC=SNV;INT;R5;GNO;FREQ=1000Genomes:0.9998,0.0001562

[W::vcf_parse_info] Extreme INFO/RS value encountered and set to missing at NC_000001.11:6259533

The ID value rs2148352424 seems to represent the actual RS value, which pysam is converting to RS=. in the INFO column since it's too large.

Oxbow is using the noodles library to parse some file formats, including VCF files. Noodles aims to be spec compliant, and appears to be following the VCF 4.4 spec, which uses 32-bit signed integers.

Any attempts to parse numbers greater than 2^31-1 would cause errors.

Next, I looked through the contents of the VCF records to confirm if there are indeed such large numbers present.

# pysam_print_vcf.py
import pysam
with pysam.VariantFile("./GCF_000001405.40.gz") as f:
    for record in f.fetch():
        print(record)

$ python pysam_print_vcf.py > vcf_out  # manually interrupted since it's long running
$ grep "RS=\." vcf_out

Output snippet:

...
NC_000001.11    23815613        rs2148423122    A       G       .       .       RS=.;dbSNPBuildID=156;SSR=0;GENEINFO=HMGCL:3155;VC=SNV;INT
NC_000001.11    23815620        rs2148423131    T       C       .       .       RS=.;dbSNPBuildID=156;SSR=0;GENEINFO=HMGCL:3155;VC=SNV;INT
NC_000001.11    23815640        rs2148423145    T       C       .       .       RS=.;dbSNPBuildID=156;SSR=0;GENEINFO=HMGCL:3155;VC=SNV;INT
...

There are many RS values that exceed the 32-bit signed integer limit.

Unfortunately, this means that you won't be able to parse this particular VCF file with oxbow due to upstream spec conformance.

If you want to view the data by setting the too-large values as missing, you can use another tool like pysam or our VCF dataframe flattener, VCFormer, which wraps pysam.

from oxbow.

Thank you for the quick reply @GarrettNg!
I'll check out vcformer. Do you think it would make sense to report this upstream with noodles?

from oxbow.

Makes sense. I'll see if I can report the incorrect format to ncbi.

from oxbow.