chore: add clnvid cache rule (#1576)

#### Added

- rule and script for adding CLNVID to clinvar INFO field

Author: mathiasbio

.github/workflows/pytest_and_coveralls.yml

@@ -61,4 +61,4 @@ jobs:
           file: ./coverage.xml
           flags: unittests
           fail_ci_if_error: true
-          verbose: true
+          verbose: true

BALSAMIC/assets/scripts/add_clnvid_field.py

@@ -0,0 +1,82 @@
+import pysam
+import click
+
+
+def add_clnvid_header(output_handle) -> None:
+    """
+    Writes the INFO header line for the CLNVID field to the output VCF.
+    """
+    # Use String to be safe; ID column can contain non-numeric identifiers.
+    header_line = '##INFO=<ID=CLNVID,Number=1,Type=String,Description="ClinVar Variation ID taken from the VCF ID column">'
+    output_handle.write(f"{header_line}\n".encode("utf-8"))
+
+
+def process_vcf(input_path: str, output_path: str) -> None:
+    """
+    Processes a bgzipped VCF file using pysam, adds the CLNVID INFO field based on the ID column,
+    and writes to a new bgzipped VCF file in a tabix-compatible format.
+    """
+    saw_clnvid_header = False
+
+    with pysam.BGZFile(input_path, "r") as infile, pysam.BGZFile(
+        output_path, "w"
+    ) as outfile:
+        for raw_line in infile:
+            # Pass through meta headers; track if CLNVID header already exists
+            if raw_line.startswith(b"##"):
+                if b"##INFO=<ID=CLNVID" in raw_line:
+                    saw_clnvid_header = True
+                outfile.write(raw_line)
+                continue
+
+            # Column header line: add CLNVID header if missing, then write
+            if raw_line.startswith(b"#"):
+                if not saw_clnvid_header:
+                    add_clnvid_header(outfile)
+                    saw_clnvid_header = True
+                outfile.write(raw_line)
+                continue
+
+            # Variant line
+            line = raw_line.decode("utf-8").rstrip("\n")
+            fields = line.split("\t")
+
+            # Ensure we have at least up to INFO column
+            if len(fields) < 8:
+                # Malformed line, write back unchanged
+                outfile.write((line + "\n").encode("utf-8"))
+                continue
+
+            vcf_id = fields[2]
+            info = fields[7]
+
+            # Only add CLNVID when ID column is not '.'
+            if vcf_id != ".":
+                if info == "." or info == "":
+                    info = f"CLNVID={vcf_id}"
+                elif "CLNVID=" not in info:
+                    info = f"{info};CLNVID={vcf_id}"
+
+            # Write updated INFO back to fields
+            fields[7] = info
+            modified_line = "\t".join(fields) + "\n"
+            outfile.write(modified_line.encode("utf-8"))
+
+
+@click.command()
+@click.argument(
+    "input_path", type=click.Path(exists=True, readable=True, dir_okay=False)
+)
+@click.argument("output_path", type=click.Path(writable=True, dir_okay=False))
+def main(input_path: str, output_path: str) -> None:
+    """
+    Adds a CLNVID INFO field to each record in a bgzipped VCF file based on the ID column.
+
+    INPUT_PATH: Path to the input VCF file (.vcf.gz).
+    OUTPUT_PATH: Path to the output VCF file (.vcf.gz).
+    """
+    process_vcf(input_path, output_path)
+
+
+if __name__ == "__main__":
+    main()

BALSAMIC/commands/init/base.py

@@ -3,6 +3,7 @@
 import logging
 import subprocess
 import sys
+import shlex
 from datetime import datetime
 from pathlib import Path
 from typing import List, Optional, Union

BALSAMIC/commands/run/analysis.py

@@ -3,6 +3,7 @@
 import logging
 import os
 import subprocess
+import shlex
 import sys
 from pathlib import Path
 from typing import List

BALSAMIC/constants/cache.py

@@ -390,7 +390,7 @@ class DockerContainers(StrEnum):
             "url": "https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz",
             "file_type": FileType.VCF,
             "gzip": True,
-            "file_name": "clinvar.vcf",
+            "file_name": "clinvar_preprocess.vcf",
             "dir_name": "variants",
         },
         "somalier_sites": {

BALSAMIC/models/cache.py

@@ -1,4 +1,5 @@
 """Balsamic reference cache models."""
+import re
 import logging
 from pathlib import Path
 from typing import Dict, Optional, List, Union
@@ -148,6 +149,12 @@ def get_cadd_snv_file_paths(self) -> List[str]:
         """Return CADD SNV reference output files."""
         return [self.cadd_snv.file_path, f"{self.cadd_snv.file_path}.{FileType.TBI}"]
 
+    def get_processed_clinvar_file_path(self) -> str:
+        return re.sub(r"\.vcf$", "_processed.vcf.gz", self.clinvar.file_path)
+
+    def get_processed_clinvar_index_path(self) -> str:
+        return re.sub(r"\.vcf$", "_processed.vcf.gz.tbi", self.clinvar.file_path)
+
     def get_delly_file_paths(self) -> List[str]:
         """Return Delly associated output files."""
         return [
@@ -400,12 +407,13 @@ def get_reference_output_paths(self) -> List[str]:
             self.references.ascat_gc_correction.file_path,
             self.references.cadd_snv.file_path,
             self.references.simple_repeat.file_path,
-            f"{self.references.clinvar.file_path}.{FileType.GZ}",
             f"{self.references.cosmic.file_path}.{FileType.GZ}",
             f"{self.references.dbsnp.file_path}.{FileType.GZ}",
             self.references.rank_score.file_path,
             f"{self.references.somalier_sites.file_path}.{FileType.GZ}",
             self.references.wgs_calling_regions.file_path,
+            self.references.get_processed_clinvar_file_path(),
+            self.references.get_processed_clinvar_index_path(),
             *self.get_compressed_indexed_vcf_paths(),
             *self.references.get_1k_genome_file_paths(),
             *self.references.get_cadd_snv_file_paths(),
@@ -434,7 +442,7 @@ def get_analysis_references(
             ascat_gc_correction=self.references.ascat_gc_correction.file_path,
             cadd_snv=self.references.cadd_snv.file_path,
             simple_repeat=self.references.simple_repeat.file_path,
-            clinvar=f"{self.references.clinvar.file_path}.{FileType.GZ}",
+            clinvar=self.references.get_processed_clinvar_file_path(),
             cosmic=f"{self.references.cosmic.file_path}.{FileType.GZ}",
             dbsnp=f"{self.references.dbsnp.file_path}.{FileType.GZ}",
             delly_exclusion=self.references.delly_exclusion.file_path,

BALSAMIC/snakemake_rules/cache/reference_vcf.rule

@@ -1,5 +1,6 @@
 """Rules to process VCF reference files."""
 
+from BALSAMIC.utils.rule import get_script_path
 
 wildcard_constraints:
     vcf="|".join(cache_config.get_reference_file_paths_by_file_type(file_type=FileType.VCF)),
@@ -45,3 +46,42 @@ rule index_vcfs:
         """
         tabix -p vcf "{input.vcf_gz}" &> "{log}"
         """
+
+rule postprocess_clinvar:
+    """Add new info-field based on clinvar variant ID: CLNVID"""
+    input:
+        clinvar=f"{cache_config.references.clinvar.file_path}.{FileType.GZ}",
+    output:
+        clinvar=cache_config.references.get_processed_clinvar_file_path(),
+    params:
+        add_clnvid_script=get_script_path("add_clnvid_field.py"),
+    message:
+        "Adding INFO field CLNVID to ClinVar VCF: {input.clinvar}"
+    benchmark:
+        f"{cache_config.references.clinvar.file_path}.postprocess_clinvar.benchmark.{FileType.TSV}"
+    log:
+        f"{cache_config.references.clinvar.file_path}.postprocess_clinvar.{FileType.LOG}",
+    shell:
+        """
+        python {params.add_clnvid_script} {input.clinvar} {output.clinvar}
+        """
+
+rule index_postprocessed_clinvar_vcf:
+    """Index postprocessed clinvar VCF"""
+    input:
+        clinvar=cache_config.references.get_processed_clinvar_file_path(),
+        singularity_image=f"{config['containers_dir']}/{config['bioinfo_tools']['tabix']}.{FileType.SIF}",
+    output:
+        clinvar_tbi=cache_config.references.get_processed_clinvar_index_path(),
+    singularity:
+        f"{config['containers_dir']}/{config['bioinfo_tools']['tabix']}.{FileType.SIF}"
+    message:
+        "Index postprocessed ClinVar VCF: {input.clinvar}"
+    benchmark:
+        f"{cache_config.references.clinvar.file_path}.tabix_postprocess_clinvar.benchmark.{FileType.TSV}"
+    log:
+        f"{cache_config.references.clinvar.file_path}.tabix_postprocess_clinvar.{FileType.LOG}",
+    shell:
+        """
+        tabix -p vcf {input.clinvar}
+        """

BALSAMIC/workflows/reference.smk

@@ -12,7 +12,9 @@ from BALSAMIC.constants.rules import SNAKEMAKE_RULES
 from BALSAMIC.models.cache import CacheConfig, AnalysisReferences
 from BALSAMIC.utils.io import write_finish_file, write_json
 from BALSAMIC.utils.utils import get_relative_paths_dict
-
+from BALSAMIC.utils.rule import (
+    get_script_path,
+)
 LOG = logging.getLogger(__name__)
 
 # Balsamic cache configuration model

CHANGELOG.rst

@@ -22,6 +22,7 @@ Changed:
 * use median_target_coverage as threshold for wgs normals coverage instead of median coverage https://github.com/Clinical-Genomics/BALSAMIC/pull/1587
 
 
+
 Removed:
 ^^^^^^^^
 * removed immediate submit functionality https://github.com/Clinical-Genomics/BALSAMIC/pull/1604
@@ -44,6 +45,7 @@ Added:
 
 
 
+
 Changed:
 ^^^^^^^^
 

tests/conftest.py

@@ -2300,12 +2300,14 @@ def fixture_delly_exclusion_converted_file(session_tmp_path: Path) -> Path:
     return reference_file
 
 
-@pytest.fixture(scope="session", name="clinvar_file")
+@pytest.fixture(scope="session", name="clinvar_processed")
 def fixture_clinvar_file(session_tmp_path: Path) -> Path:
     """Return dummy clinvar file."""
-    clinvar_file: Path = Path(session_tmp_path, "variants", "clinvar.vcf.gz")
-    clinvar_file.touch()
-    return clinvar_file
+    clinvar_processed: Path = Path(
+        session_tmp_path, "variants", "clinvar_processed.vcf.gz"
+    )
+    clinvar_processed.touch()
+    return clinvar_processed
 
 
 @pytest.fixture(scope="session", name="cosmic_file")
@@ -2379,7 +2381,7 @@ def fixture_analysis_references_hg_data(
     cache_config: CacheConfig,
     analysis_references_data: Dict[str, Path],
     delly_exclusion_converted_file: Path,
-    clinvar_file: Path,
+    clinvar_processed: Path,
     cosmic_file: Path,
     dbsnp_file: Path,
     hc_vcf_1kg_file: Path,
@@ -2397,7 +2399,7 @@ def fixture_analysis_references_hg_data(
         ),
         "cadd_snv": Path(cache_config.references.cadd_snv.file_path),
         "simple_repeat": Path(cache_config.references.simple_repeat.file_path),
-        "clinvar": clinvar_file,
+        "clinvar": clinvar_processed,
         "cosmic": cosmic_file,
         "dbsnp": dbsnp_file,
         "delly_exclusion": Path(cache_config.references.delly_exclusion.file_path),