feat: script fix_glnexus.py prepare GLNexus output for noodles (#356)

Author: holtgrewe

misc/fix_glnexus.py

@@ -0,0 +1,64 @@
+#!/usr/bin/env python3
+"""
+Fix GLNexus output for input to `mehari annotate seqvars`.
+
+First, note that multiallelic sites must be corrected by a call to
+`bcftools norm` already as an input to this script.
+
+What we will do is splitting the `FORMAT/RNC` field with a comma
+as noodles-vcf expects a list of characters for rather than a
+string of length 2 for `##FORMAT=<ID=RNC,Number=2,Type=Character>`.
+"""
+
+from pathlib import Path
+import sys
+from typing import Annotated
+
+import typer
+import vcfpy
+
+
+def main(
+    path_in: Annotated[Path, typer.Option(help="Path to input VCF")],
+    path_out: Annotated[Path, typer.Option(help="Path to output VCF")],
+    quiet: Annotated[bool, typer.Option(help="Disable verbose output")] = False,
+):
+    """
+    Fix GLNexus output to be compatible with `noodles-vcf`.
+
+    cf: https://github.com/varfish-org/mehari/issues/356
+    """
+    if not quiet:
+        print(f"Opening input file {path_in}", file=sys.stderr)
+    reader = vcfpy.Reader.from_path(path_in)
+    header = reader.header.copy()
+    for line in header.lines:
+        if line.key == "FORMAT" and line.mapping.get("ID") == "GQ":
+            line.mapping["Number"] = 1
+            line.mapping["Type"] = "Integer"
+    if not quiet:
+        print(f"Opening output file {path_out}", file=sys.stderr)
+    writer = vcfpy.Writer.from_path(path_out, header)
+    if not quiet:
+        print("Processing records...", file=sys.stderr)
+    with reader, writer:
+        for idx, record in enumerate(reader):
+            if idx % 10_000 == 0:
+                print(
+                    f"  at {idx} records {record.CHROM}:{record.POS}", file=sys.stderr
+                )
+            for call in record.calls:
+                if "RNC" in call.data:
+                    if (
+                        isinstance(call.data["RNC"], list)
+                        and len(call.data["RNC"]) == 1
+                    ):
+                        call.data["RNC"] = list(call.data["RNC"][0])
+            writer.write_record(record)
+    if not quiet:
+        print("... done", file=sys.stderr)
+        print("All done. Have a nice day!", file=sys.stderr)
+
+
+if __name__ == "__main__":
+    typer.run(main)

tests/data/annotate/seqvars/clair3-glnexus-min.multiallelic.vcf

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:f8adde9396d07ce7379922e57f7c9ca6e0b96ea7b4d1203075f3cb2d017c03f1
-size 5436
+oid sha256:acdcc6065fcab6156f630b1fe27afb9cedd09adac0f84215a879120ae84a3162
+size 5548

tests/data/annotate/seqvars/clair3-glnexus-min.original.vcf

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f8adde9396d07ce7379922e57f7c9ca6e0b96ea7b4d1203075f3cb2d017c03f1
+size 5436

tests/data/annotate/seqvars/clair3-glnexus-min.tsv

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d82cbc30da246c9ecbfcf41cf166713d0af30245d97d8b2d51a534354453745a
+size 3709

tests/data/annotate/seqvars/clair3-glnexus-min.vcf

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:6a9770c6d67d60c1aea6c5ac2e22c949bd7bced118e632151844535031547500
-size 5677
+oid sha256:b7ae1ba6ba0989506e92fbf1fef393fc4eb7e66cba2add3a8dbfb540a7571fbe
+size 5791