Reading VCF files?

[ap--]

Hello,

I wanted to try oxbow and was testing on the following VCF file, but it doesn't seem to work. I'm not sure if this should be supported already, or if I am missing something.

I built the oxbow wheel from current main:

$ git rev-parse HEAD
e3d2a1751901430a16438134b87bc16f21d90269

files

$ wget https://ftp.ncbi.nlm.nih.gov/snp/latest_release/VCF/GCF_000001405.40.gz 
$ wget https://ftp.ncbi.nlm.nih.gov/snp/latest_release/VCF/GCF_000001405.40.gz.tbi
$ ls -l GCF_000001405.40.* 
-rw-r--r--  1 andreaspoehlmann  staff  26611209012 Oct 16 10:51 GCF_000001405.40.gz
-rw-r--r--  1 andreaspoehlmann  staff      3118040 Oct 16 10:58 GCF_000001405.40.gz.tbi
$ md5sum GCF_000001405.40.*
a1082ca70e15eb63301dfc33b19d0ae7  GCF_000001405.40.gz
76959b1691e8e62cd650664b00b7ea02  GCF_000001405.40.gz.tbi

code

# read_vcf.py                                                                       
import importlib.metadata
import oxbow as ox
import polars as pl

print("oxbow.__version__", importlib.metadata.version("oxbow"))

ipc = ox.read_vcf("GCF_000001405.40.gz", index="GCF_000001405.40.gz.tbi")
df = pl.read_ipc(ipc)
print(df)

error

$ python read_vcf.py
oxbow.__version__ 0.2.0
thread '<unnamed>' panicked at 'called `Result::unwrap()` on an `Err` value: ExternalError(Custom { kind: InvalidData, error: InvalidInfo(InvalidField(InvalidValue(Other(Other("RS")), InvalidInteger(ParseIntError { kind: PosOverflow })))) })', src/lib.rs:117:49
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
Traceback (most recent call last):
  File "/Users/andreaspoehlmann/development/oxbow-test/read_vcf.py", line 8, in <module>
    ipc = ox.read_vcf("GCF_000001405.40.gz", index="GCF_000001405.40.gz.tbi")
pyo3_runtime.PanicException: called `Result::unwrap()` on an `Err` value: ExternalError(Custom { kind: InvalidData, error: InvalidInfo(InvalidField(InvalidValue(Other(Other("RS")), InvalidInteger(ParseIntError { kind: PosOverflow })))) })

system info

$ python -VV
Python 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:41:52) [Clang 15.0.7 ]
$ uname -a
Darwin F2WR4P9QNH 23.0.0 Darwin Kernel Version 23.0.0: Fri Sep 15 14:43:05 PDT 2023; root:xnu-10002.1.13~1/RELEASE_ARM64_T6020 arm64
$ system_profiler SPHardwareDataType | grep -e "Model\|Memory\|Cores"
      Model Name: MacBook Pro
      Model Identifier: Mac14,5
      Total Number of Cores: 12 (8 performance and 4 efficiency)
      Memory: 64 GB

Cheers, Andreas 😃

[GarrettNg]

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.

[ap--]

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?

[GarrettNg]

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.

[ap--]

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