From 5d39e3a77ff00e7b75877546fdde593bb18b5cfd Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sat, 12 Sep 2026 12:46:02 -0700 Subject: [PATCH 1/2] Detect the pangenome reference from the GBZ metadata A wrong `--pangenome_ref_name` was silent: `vg haplotypes --set-reference` only warns, the sampled graph loses its reference paths, `vg convert -Q` writes no rGFA tags and `pgutil lift` exits 0 with every read unmapped. Read the reference name and contig prefix from the GBWT metadata in the GBZ header instead (new `pangenome_meta` module, well under a second per graph). The backbone reference is the unfragmented reference sample. Make `--pangenome_ref_name` and `--pangenome_contig_prefix` hidden overrides that must agree with the graph, check the backbone contigs and build against the reference FASTA, replace the `grch38.gbz` filename-suffix check, apply the GRCh38 contig-length check only to GRCh38 backbones, and add `check-sample-gbz` / `check-sample-gfa` jobs to both DAGs so a sampled graph or GFA without the reference fails the run. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01Lo4R3WSNwCYMcKzvGKCqSs --- sentieon_cli/base_pangenome.py | 342 +++++++- sentieon_cli/command_strings.py | 61 +- sentieon_cli/hybrid_pangenome.py | 98 +-- sentieon_cli/pangenome_meta.py | 782 ++++++++++++++++++ sentieon_cli/sentieon_pangenome.py | 104 +-- tests/data/pangenome/generate.sh | 40 + tests/data/pangenome/rgfa_no_header.gfa | 3 + tests/data/pangenome/rgfa_no_tags.gfa | 4 + tests/data/pangenome/rgfa_ok.gfa | 6 + tests/data/pangenome/rgfa_wrong_reference.gfa | 6 + tests/data/pangenome/tiny.gbz | Bin 0 -> 2184 bytes tests/data/pangenome/tiny.gfa | 28 + tests/data/pangenome/tiny.hapl | Bin 0 -> 760 bytes tests/data/pangenome/tiny_chm13.gbz | Bin 0 -> 2184 bytes tests/data/pangenome/tiny_chm13.gfa | 28 + tests/data/pangenome/tiny_no_backbone.gbz | Bin 0 -> 2240 bytes tests/data/pangenome/tiny_no_backbone.gfa | 30 + tests/data/pangenome/tiny_no_ref_tag.gbz | Bin 0 -> 2152 bytes tests/data/pangenome/tiny_no_ref_tag.gfa | 28 + tests/data/pangenome/tiny_sampled.gbz | Bin 0 -> 2184 bytes tests/unit/test_hybrid_pangenome.py | 70 ++ tests/unit/test_pangenome_meta.py | 394 +++++++++ tests/unit/test_pipeline_validation.py | 331 +++++++- tests/unit/test_sentieon_pangenome.py | 79 ++ ...test_sentieon_pangenome_kmc_integration.py | 7 + 25 files changed, 2279 insertions(+), 162 deletions(-) create mode 100644 sentieon_cli/pangenome_meta.py create mode 100755 tests/data/pangenome/generate.sh create mode 100644 tests/data/pangenome/rgfa_no_header.gfa create mode 100644 tests/data/pangenome/rgfa_no_tags.gfa create mode 100644 tests/data/pangenome/rgfa_ok.gfa create mode 100644 tests/data/pangenome/rgfa_wrong_reference.gfa create mode 100644 tests/data/pangenome/tiny.gbz create mode 100644 tests/data/pangenome/tiny.gfa create mode 100644 tests/data/pangenome/tiny.hapl create mode 100644 tests/data/pangenome/tiny_chm13.gbz create mode 100644 tests/data/pangenome/tiny_chm13.gfa create mode 100644 tests/data/pangenome/tiny_no_backbone.gbz create mode 100644 tests/data/pangenome/tiny_no_backbone.gfa create mode 100644 tests/data/pangenome/tiny_no_ref_tag.gbz create mode 100644 tests/data/pangenome/tiny_no_ref_tag.gfa create mode 100644 tests/data/pangenome/tiny_sampled.gbz create mode 100644 tests/unit/test_pangenome_meta.py diff --git a/sentieon_cli/base_pangenome.py b/sentieon_cli/base_pangenome.py index bf7f90f..9b6c82d 100644 --- a/sentieon_cli/base_pangenome.py +++ b/sentieon_cli/base_pangenome.py @@ -2,10 +2,11 @@ A base class for pangenome pipelines """ +import argparse import copy import pathlib import sys -from typing import Iterable, List, Optional, Sequence +from typing import Dict, Iterable, List, Optional, Sequence, Set from . import command_strings as cmds from .dag import DAG @@ -21,7 +22,14 @@ WgsMetricsAlgo, ) from .job import Job +from .pangenome_meta import ( + PangenomeMetadataError, + PangenomeReference, + detect_pangenome_reference, + read_hapl_header, +) from .pipeline import BasePipeline +from .shard import GRCH38_CONTIGS, detect_reference_build from .stages.base import StageContext from .stages.cnv import PangenomeCNVResult, PangenomeCNVStage from .stages.expansion import ExpansionHunterResult, ExpansionHunterStage @@ -36,6 +44,19 @@ ) from .util import path_arg, require_versions, sample_sex_arg +# The pangenome reference of the graphs the pipelines shipped with, used +# when a graph's own metadata cannot be read +DEFAULT_PANGENOME_REF_NAME = "GRCh38" +DEFAULT_PANGENOME_CONTIG_PREFIX = "GRCh38#0#" + +# The linear reference build that each pangenome backbone reference +# implies. Contig names alone cannot tell the builds apart, as both name +# their chromosomes chr1..chrM. +PANGENOME_REFERENCE_BUILDS: Dict[str, str] = { + "GRCh38": "hg38", + "CHM13": "chm13", +} + class BasePangenome(BasePipeline): """A pipeline base class for short reads""" @@ -173,6 +194,23 @@ class BasePangenome(BasePipeline): f"passed to T1K (default: {DEFAULT_T1K_KIR_LOCUS})." ), }, + # Hidden arguments. The pangenome reference name and contig + # prefix are read from the `--gbz` file; these override the + # detected values for a graph the detection cannot handle. + "pangenome_ref_name": { + "help": argparse.SUPPRESS, + }, + "pangenome_contig_prefix": { + "help": argparse.SUPPRESS, + }, + "skip_contig_checks": { + "help": argparse.SUPPRESS, + "action": "store_true", + }, + "skip_pangenome_name_checks": { + "help": argparse.SUPPRESS, + "action": "store_true", + }, } ) @@ -199,6 +237,17 @@ def __init__(self) -> None: self.t1k_kir_coord: Optional[pathlib.Path] = None self.t1k_kir_locus: str = DEFAULT_T1K_KIR_LOCUS self.has_cnv_model = False + # The pangenome reference, resolved from the graph's metadata by + # `resolve_pangenome_reference`. The two arguments are hidden + # overrides, so they stay `None` unless the user supplies them. + self.pangenome_ref_name: Optional[str] = None + self.pangenome_contig_prefix: Optional[str] = None + self.pangenome_reference: Optional[PangenomeReference] = None + self.skip_contig_checks = False + self.skip_pangenome_name_checks = False + # Parsed by `validate`; declared here for the shared checks + self.fai_data: Dict[str, Dict[str, int]] = {} + self.pop_vcf_contigs: Dict[str, Optional[int]] = {} # Stashed by `build_dag` for the second, sex-aware DAG self.ploidy_json: Optional[pathlib.Path] = None self.sr_alignments: List[pathlib.Path] = [] @@ -230,6 +279,297 @@ def validate_t1k(self) -> None: require_versions(T1K_MIN_VERSIONS, skip=self.skip_version_check) + def resolve_pangenome_reference(self) -> None: + """Resolve the pangenome reference name and contig prefix. + + Both values are read from the `--gbz` file's own metadata. The + hidden `--pangenome_ref_name` and `--pangenome_contig_prefix` + arguments override the detected values, and disagreeing with the + graph is an error unless `--skip_pangenome_name_checks` is set. + + The check matters because a wrong reference name is silent: `vg + haplotypes --set-reference` only warns, writes a sampled graph + with no reference paths, `vg convert -Q` then writes a GFA with + no `SN:Z:` rGFA tags, and `pgutil lift` exits 0 with every read + unmapped. + + Called from `validate` once the reference index is parsed and + before `validate_bundle`, which selects the model bundle's + `extract..model` member by the resolved name. + """ + gbz = self.required(self.gbz, "gbz") + supplied_name = self.pangenome_ref_name + supplied_prefix = self.pangenome_contig_prefix + + detected: Optional[PangenomeReference] = None + try: + detected = detect_pangenome_reference(gbz) + except (PangenomeMetadataError, OSError) as err: + self._undetected_pangenome_reference(gbz, err) + + if detected is None: + self.pangenome_ref_name = ( + supplied_name or DEFAULT_PANGENOME_REF_NAME + ) + self.pangenome_contig_prefix = ( + supplied_prefix or DEFAULT_PANGENOME_CONTIG_PREFIX + ) + else: + self._check_pangenome_overrides(detected) + self.pangenome_reference = detected + self.pangenome_ref_name = supplied_name or detected.ref_name + self.pangenome_contig_prefix = ( + supplied_prefix or detected.contig_prefix + ) + + self.logger.info( + "Using pangenome reference '%s' (contig prefix '%s')", + self.pangenome_ref_name, + self.pangenome_contig_prefix, + ) + + if detected is not None: + self.check_pangenome_reference_fasta(detected) + self.check_pangenome_hapl(detected) + + def _undetected_pangenome_reference( + self, gbz: pathlib.Path, err: Exception + ) -> None: + """Handle a `--gbz` file whose reference could not be detected. + + Dry runs (and the unit tests behind them) parse placeholder files, + so they fall back to the historical defaults. A real run needs + both overrides to continue without the detected values. + """ + if self.dry_run: + self.logger.debug( + "Could not read the pangenome metadata of '%s': %s", gbz, err + ) + return + if self.pangenome_ref_name and self.pangenome_contig_prefix: + self.logger.warning( + "Could not read the pangenome metadata of '%s': %s. Using " + "the supplied `--pangenome_ref_name` and " + "`--pangenome_contig_prefix`.", + gbz, + err, + ) + return + self.logger.error( + "Could not read the pangenome metadata of '%s': %s. Supply " + "both `--pangenome_ref_name` and `--pangenome_contig_prefix` " + "to run without the detected values.", + gbz, + err, + ) + sys.exit(2) + + def _check_pangenome_overrides(self, detected: PangenomeReference) -> None: + """Compare the supplied overrides with the detected values""" + overrides = ( + ( + "--pangenome_ref_name", + self.pangenome_ref_name, + detected.ref_name, + ), + ( + "--pangenome_contig_prefix", + self.pangenome_contig_prefix, + detected.contig_prefix, + ), + ) + mismatched = False + for flag, supplied, found in overrides: + if supplied is None or supplied == found: + continue + mismatched = True + if self.skip_pangenome_name_checks: + self.logger.warning( + "The supplied `%s` value '%s' does not match the value " + "detected in the pangenome graph, '%s'. Using the " + "supplied value.", + flag, + supplied, + found, + ) + else: + self.logger.error( + "The supplied `%s` value '%s' does not match the value " + "detected in the pangenome graph, '%s'. Drop the " + "argument to use the detected value.", + flag, + supplied, + found, + ) + if mismatched and not self.skip_pangenome_name_checks: + sys.exit(2) + + def _pangenome_check_failed(self, msg: str, *args: object) -> None: + """Fail a pangenome consistency check. + + `--skip_pangenome_name_checks` downgrades every one of these + checks to a warning. + """ + if self.skip_pangenome_name_checks: + self.logger.warning(msg, *args) + return + self.logger.error(msg, *args) + sys.exit(2) + + def check_pangenome_reference_fasta( + self, detected: PangenomeReference + ) -> None: + """Check the reference FASTA against the pangenome backbone. + + Every contig of the backbone reference must be in the reference + index, and the build the backbone implies must match the build + detected from the index. Both checks are needed: contig names + alone cannot tell GRCh38 from CHM13. + """ + missing = [ctg for ctg in detected.contigs if ctg not in self.fai_data] + if missing: + self._pangenome_check_failed( + "%d contig(s) of the pangenome reference '%s' are missing " + "from the reference FASTA index: %s", + len(missing), + detected.ref_name, + ", ".join(missing[:10]), + ) + + expected_build = PANGENOME_REFERENCE_BUILDS.get(detected.ref_name) + fai_build = detect_reference_build(self.fai_data) + if expected_build is None or fai_build is None: + self.logger.info( + "The pangenome reference is '%s' and the reference FASTA " + "build is '%s'; the two were not compared", + detected.ref_name, + fai_build, + ) + return + if expected_build != fai_build: + self._pangenome_check_failed( + "The pangenome reference '%s' expects the '%s' reference " + "build, but the `--reference` FASTA is '%s'", + detected.ref_name, + expected_build, + fai_build, + ) + + def check_pangenome_hapl( + self, detected: Optional[PangenomeReference] + ) -> None: + """Check the `--hapl` file against the `--gbz` graph. + + A mismatched top-level chain count is only a warning: the count + equals the backbone's path count for every graph seen so far, but + it is a property of the snarl decomposition rather than a + guarantee. + """ + hapl = self.required(self.hapl, "hapl") + try: + header = read_hapl_header(hapl) + except (PangenomeMetadataError, OSError) as err: + if self.dry_run: + self.logger.debug( + "Could not read the haplotype file '%s': %s", hapl, err + ) + return + self._pangenome_check_failed("%s", err) + return + + if detected is None: + return + if header.top_level_chains != detected.n_paths: + self.logger.warning( + "The `--hapl` file '%s' has %d top-level chains, but the " + "pangenome reference '%s' has %d paths. The `--gbz` and " + "`--hapl` files may not be a matching pair.", + hapl, + header.top_level_chains, + detected.ref_name, + detected.n_paths, + ) + + def validate_grch38_contigs(self) -> None: + """Check the reference and pop VCF contig lengths against GRCh38. + + The lengths only describe GRCh38; other pangenome references are + covered by the reference FASTA checks in + `resolve_pangenome_reference` instead. + """ + if self.skip_contig_checks: + return + if self.pangenome_ref_name != DEFAULT_PANGENOME_REF_NAME: + self.logger.info( + "The pangenome reference is '%s', so the GRCh38 " + "contig-length check is skipped", + self.pangenome_ref_name, + ) + return + + # Check the fai file contigs + mismatch_contigs: Set[str] = set() + for ctg, length in GRCH38_CONTIGS.items(): + fai_length = self.fai_data.get(ctg, {}).get("length", -1) + if length != fai_length: + mismatch_contigs.add(ctg) + if mismatch_contigs: + self.logger.error( + "Reference contigs with unexpected lengths: %s", + ", ".join(mismatch_contigs), + ) + sys.exit(2) + + # Check the pop VCF file contigs + if self.dry_run: + return + mismatch_contigs = set() + for ctg, length in GRCH38_CONTIGS.items(): + if length != self.pop_vcf_contigs.get(ctg, -1): + mismatch_contigs.add(ctg) + if mismatch_contigs: + self.logger.error( + "Pop VCF contigs with unexpected lengths: %s", + ", ".join(mismatch_contigs), + ) + sys.exit(2) + + def ref_name(self) -> str: + """The resolved pangenome reference name. + + `resolve_pangenome_reference` always sets it; this narrows the + optional attribute for the consumers of the graph. + """ + return self.required(self.pangenome_ref_name, "pangenome_ref_name") + + def contig_prefix(self) -> str: + """The resolved pangenome contig prefix""" + return self.required( + self.pangenome_contig_prefix, "pangenome_contig_prefix" + ) + + def build_check_gbz_job(self, sample_gbz: pathlib.Path) -> Job: + """Check the sampled pangenome before anything consumes it""" + return Job( + cmds.cmd_check_pangenome_gbz( + sample_gbz, self.ref_name(), self.contig_prefix() + ), + "check-sample-gbz", + 1, + task_name="pangenome", + ) + + def build_check_gfa_job(self, sample_gfa: pathlib.Path) -> Job: + """Check the converted GFA before anything consumes it""" + return Job( + cmds.cmd_check_pangenome_gfa( + sample_gfa, self.ref_name(), self.contig_prefix() + ), + "check-sample-gfa", + 1, + task_name="pangenome", + ) + def output_path(self, suffix: str) -> pathlib.Path: """A path next to the output VCF, with `.vcf.gz` replaced""" output_vcf = self.required(self.output_vcf, "output_vcf") diff --git a/sentieon_cli/command_strings.py b/sentieon_cli/command_strings.py index 9b23910..c208463 100644 --- a/sentieon_cli/command_strings.py +++ b/sentieon_cli/command_strings.py @@ -319,6 +319,63 @@ def cmd_pyexec_sad_lad_update( return Pipeline(Command(*cmd)) +PANGENOME_META_MODULE = "sentieon_cli.pangenome_meta" + +# `sentieon_cli/__init__.py` imports every pipeline, so `pangenome_meta` is +# already in `sys.modules` when `-m` re-executes it as `__main__`. The +# duplicate module is harmless -- the script holds no state -- but runpy +# warns about it on the stderr of every check job, so runpy's own warnings +# are silenced rather than left to clutter the job logs. +PANGENOME_META_WARNING_FILTER = "ignore::RuntimeWarning:runpy" + + +def cmd_pangenome_meta( + subcommand: str, + input_flag: str, + input_file: pathlib.Path, + reference_name: str, + contig_prefix: str, +) -> Pipeline: + """Run a `sentieon_cli.pangenome_meta` check on a pangenome file""" + cmd = [ + sys.executable, + "-W", + PANGENOME_META_WARNING_FILTER, + "-m", + PANGENOME_META_MODULE, + subcommand, + input_flag, + str(input_file), + "--reference_name", + reference_name, + "--contig_prefix", + contig_prefix, + ] + return Pipeline(Command(*cmd)) + + +def cmd_check_pangenome_gbz( + gbz: pathlib.Path, + reference_name: str, + contig_prefix: str, +) -> Pipeline: + """Check that a sampled pangenome kept its reference paths""" + return cmd_pangenome_meta( + "check-gbz", "--gbz", gbz, reference_name, contig_prefix + ) + + +def cmd_check_pangenome_gfa( + gfa: pathlib.Path, + reference_name: str, + contig_prefix: str, +) -> Pipeline: + """Check that a pangenome GFA carries the reference's rGFA tags""" + return cmd_pangenome_meta( + "check-gfa", "--gfa", gfa, reference_name, contig_prefix + ) + + def cmd_pyexec_indel2cnv( out_vcf: pathlib.Path, in_vcf: pathlib.Path, @@ -1275,8 +1332,8 @@ def cmd_bwa_extract( def cmd_vg_convert_gfa( output_gfa: pathlib.Path, input_gbz: pathlib.Path, - reference_name="GRCh38", - threads=1, + reference_name: str = "GRCh38", + threads: int = 1, ) -> Pipeline: """Convert GBZ to GFA""" cmd = [ diff --git a/sentieon_cli/hybrid_pangenome.py b/sentieon_cli/hybrid_pangenome.py index 78a1ab3..ea83b2b 100644 --- a/sentieon_cli/hybrid_pangenome.py +++ b/sentieon_cli/hybrid_pangenome.py @@ -27,7 +27,6 @@ from .job import Job from .logging import get_logger from .shard import ( - GRCH38_CONTIGS, detect_reference_build, determine_shards_from_fai, parse_fai, @@ -151,16 +150,6 @@ class HybridPangenome(BasePangenome): ), "type": path_arg(exists=True, is_file=True), }, - "pangenome_contig_prefix": { - "default": "GRCh38#0#", - "help": ( - "Prefix to strip from pangenome contig names (GRCh38#0#)" - ), - }, - "pangenome_ref_name": { - "default": "GRCh38", - "help": "Reference name in the pangenome (GRCh38).", - }, "rgsm": { "help": ( "Overwrite the SM tag of the input readgroups for " @@ -185,18 +174,10 @@ class HybridPangenome(BasePangenome): "type": path_arg(exists=True, is_file=True), }, # Hidden arguments - "skip_contig_checks": { - "help": argparse.SUPPRESS, - "action": "store_true", - }, "skip_model_apply": { "help": argparse.SUPPRESS, "action": "store_true", }, - "skip_pangenome_name_checks": { - "help": argparse.SUPPRESS, - "action": "store_true", - }, "skip_pop_vcf_id_check": { "help": argparse.SUPPRESS, "action": "store_true", @@ -223,8 +204,6 @@ def __init__(self) -> None: self.call_cnvs = False self.lr_align_input = False self.lr_input_ref: Optional[pathlib.Path] = None - self.pangenome_contig_prefix = "GRCh38#0#" - self.pangenome_ref_name = "GRCh38" self.rgsm: Optional[str] = None self.extract_model_name = "extract.model" self.gfa2fa_with_vg = False @@ -235,9 +214,7 @@ def __init__(self) -> None: self.sr_readgroups: List[List[Dict[str, str]]] = [] self.lr_readgroups: List[List[Dict[str, str]]] = [] self.sample_sm = "" - self.skip_contig_checks = False self.skip_model_apply = False - self.skip_pangenome_name_checks = False self.skip_pop_vcf_id_check = False self.skip_small_variants = False self.skip_svs = False @@ -264,7 +241,10 @@ def validate(self) -> None: self.shards = determine_shards_from_fai( self.fai_data, 10 * 1000 * 1000 ) - self.pop_vcf_contigs: Dict[str, Optional[int]] = {} + # Before `validate_bundle`, which picks the bundle's + # `extract..model` member by the resolved name + self.resolve_pangenome_reference() + self.pop_vcf_contigs = {} if self.pop_vcf: self.pop_vcf_contigs = vcf_contigs(self.pop_vcf, self.dry_run) self.logger.debug("VCF contigs are: %s", self.pop_vcf_contigs) @@ -322,51 +302,7 @@ def validate(self) -> None: "across decoy and unplaced contigs." ) - if not self.skip_pangenome_name_checks: - if not str(self.gbz).endswith("grch38.gbz"): - self.logger.error( - "The `--gbz` file does not have the expected suffix. " - "Check that you are using a GRCh38 pangenome." - ) - sys.exit(2) - - if not str(self.hapl).endswith("grch38.hapl"): - self.logger.error( - "The `--hapl` file does not have the expected suffix. " - "Check that you are using a GRCh38 pangenome." - ) - sys.exit(2) - - if not self.skip_contig_checks: - # Check the fai file contigs - mismatch_contigs: Set[str] = set() - for ctg, length in GRCH38_CONTIGS.items(): - d = self.fai_data.get(ctg, {}) - fai_length = d.get("length", -1) - if length != fai_length: - mismatch_contigs.add(ctg) - if mismatch_contigs: - mismatch_contigs_s = ", ".join(mismatch_contigs) - self.logger.error( - "Reference contigs with unexpected lengths: %s", - mismatch_contigs_s, - ) - sys.exit(2) - - # Check the pop VCF file contigs - if not self.dry_run: - mismatch_contigs = set() - for ctg, length in GRCH38_CONTIGS.items(): - vcf_length = self.pop_vcf_contigs.get(ctg, -1) - if length != vcf_length: - mismatch_contigs.add(ctg) - if mismatch_contigs: - mismatch_contigs_s = ", ".join(mismatch_contigs) - self.logger.error( - "Pop VCF contigs with unexpected lengths: %s", - mismatch_contigs_s, - ) - sys.exit(2) + self.validate_grch38_contigs() # After the version gate and the contig checks, so the tool choice # is only logged for runs that proceed. The `reference_build` it @@ -860,9 +796,19 @@ def build_dag(self) -> DAG: haplotypes_job = self.build_haplotypes_job(sample_pangenome, kmer_file) dag.add_job(haplotypes_job, haplotype_dependencies) + # Confirm the sampled pangenome kept its reference paths before + # anything consumes it + check_gbz_job = self.build_check_gbz_job(sample_pangenome) + dag.add_job(check_gbz_job, {haplotypes_job}) + # convert the sample pangenome to GFA gfa_job = self.build_gfa_job(hap_raw_gfa, sample_pangenome) - dag.add_job(gfa_job, {haplotypes_job}) + dag.add_job(gfa_job, {check_gbz_job}) + + # Confirm the GFA carries the reference's rGFA tags, without which + # `pgutil lift` leaves every read unmapped + check_gfa_job = self.build_check_gfa_job(hap_raw_gfa) + dag.add_job(check_gfa_job, {gfa_job}) # Graph update without the SV BED update_raw_job = self.build_graph_update_job( @@ -871,7 +817,7 @@ def build_dag(self) -> DAG: calling_lr, name="graph-update-raw", ) - dag.add_job(update_raw_job, {gfa_job} | realign_jobs) + dag.add_job(update_raw_job, {check_gfa_job} | realign_jobs) # Call SVs from the long reads and collect graph update regions longreadsv_result = LongReadSVStage( @@ -1224,7 +1170,7 @@ def build_haplotypes_job( "--include-reference", "--diploid-sampling", "--set-reference", - self.pangenome_ref_name, + self.ref_name(), ], ), "vg-haplotypes", @@ -1242,7 +1188,7 @@ def build_gfa_job( output_gfa, input_gbz, threads=self.cores, - reference_name=self.pangenome_ref_name, + reference_name=self.ref_name(), ), "vg-convert-gfa", 0, @@ -1269,7 +1215,7 @@ def build_graph_update_job( out_gfa, gfa_file=in_gfa, target_bed=bed, - prefix=self.pangenome_contig_prefix, + prefix=self.contig_prefix(), ) ) return Job( @@ -1334,7 +1280,7 @@ def build_minimap2_lift_job( self.model_bundle.joinpath("minimap2.model"), threads=self.cores, mm2_xargs=["--secondary=yes"], - lift_prefix=self.pangenome_contig_prefix, + lift_prefix=self.contig_prefix(), ), "mm2-lift", self.cores, @@ -1395,7 +1341,7 @@ def build_pangenomesv_job( out_vcf, gfa_file=pangenome_gfa, min_af=PANGENOME_SV_MIN_AF, - prefix=self.pangenome_contig_prefix, + prefix=self.contig_prefix(), ) ) return Job( diff --git a/sentieon_cli/pangenome_meta.py b/sentieon_cli/pangenome_meta.py new file mode 100644 index 0000000..000eb42 --- /dev/null +++ b/sentieon_cli/pangenome_meta.py @@ -0,0 +1,782 @@ +""" +Read pangenome graph metadata without loading the graph. + +A GBZ file is a `simple-sds` serialization of:: + + GBZ header (16 B) | GBZ tags | GBWT | GBWTGraph sequences ... + +and the GBWT inside it is:: + + GBWT header (48 B) | GBWT tags | BWT | DA samples | Option + +Every component carries its own length, so the multi-GB BWT body is +skipped with a seek and only the few MB holding the path names are read. +Parsing a whole-genome HPRC graph this way takes well under a second. + +The metadata answers the two questions the pangenome pipelines ask before +they run: + +* `pangenome_ref_name`, the sample name that `vg haplotypes + --set-reference` and `vg convert -Q` expect. A wrong name is not an + error for `vg`: it warns, writes a sampled graph with an empty + `reference_samples` tag and no reference paths, and then converts a GFA + with no `SN:Z:` rGFA tags. `pgutil lift` exits 0 with every read + unmapped, so the run completes with silently ruined accuracy. +* `pangenome_contig_prefix`, the `##` prefix of the + reference path names, stripped by `pgutil lift` and by the `PangenomeSV` + and `PGHapUpdate` algos. + +The backbone reference is the sample in the `reference_samples` tag whose +paths are unfragmented: one path per contig, all with fragment 0. Any +other listed sample is chopped into fragments by the graph construction. + +The module also runs as a script, for the run-time checks the pipelines +wire into their DAGs:: + + python -m sentieon_cli.pangenome_meta detect --gbz graph.gbz + python -m sentieon_cli.pangenome_meta check-gbz --gbz sample.gbz + --reference_name GRCh38 --contig_prefix 'GRCh38#0#' + python -m sentieon_cli.pangenome_meta check-gfa --gfa sample.gfa + --reference_name GRCh38 --contig_prefix 'GRCh38#0#' +""" + +from __future__ import annotations + +import argparse +import pathlib +import struct +import sys +from collections import defaultdict +from dataclasses import dataclass, field +from typing import BinaryIO, Dict, List, Optional, Sequence, Tuple, Union + +from .logging import get_logger + +# Spelled out rather than taken from `__name__`: run as a script with +# `python -m sentieon_cli.pangenome_meta`, `__name__` is "__main__" and the +# logger would no longer be a child of the package logger, so its messages +# would miss the console handler and the run log. +logger = get_logger("sentieon_cli.pangenome_meta") + +# Bytes in a `simple-sds` word +WORD = 8 + +# Magic numbers of the serialized structures +GBZ_TAG = 0x205A4247 +GBWT_TAG = 0x6B376B37 +METADATA_TAG = 0x6B375E7A +HAPL_MAGIC = 0x4C504148 + +# How much of a GFA `check_gfa` reads before giving up on the `SN:Z:` tags +DEFAULT_MAX_GFA_LINES = 100000 + +# A GBWT path name record: (sample, contig, phase, fragment) +PathName = Tuple[int, int, int, int] + +PathLike = Union[str, pathlib.Path] + + +class PangenomeMetadataError(ValueError): + """A pangenome file could not be parsed, or does not describe one + unambiguous reference sample.""" + + +def unpack_bits(words: bytes, count: int, width: int) -> List[int]: + """Unpack `count` little-endian `width`-bit integers from 64-bit words. + + `simple-sds` packs an integer vector into a stream of 64-bit words with + no padding, so an element may straddle a word boundary. The running + buffer keeps the decode linear in `count`. + """ + mask = (1 << width) - 1 + out: List[int] = [] + buf = 0 + nbuf = 0 + offset = 0 + while len(out) < count: + while nbuf < width: + end = offset + WORD + buf |= int.from_bytes(words[offset:end], "little") << nbuf + nbuf += 64 + offset = end + out.append(buf & mask) + buf >>= width + nbuf -= width + return out + + +class _Reader: + """A cursor over the `simple-sds` serialization of a GBZ or GBWT. + + Every method reads (or skips) exactly one serialized structure and + leaves the cursor on the next one. + """ + + def __init__(self, fh: BinaryIO, path: PathLike) -> None: + self.fh = fh + self.path = str(path) + self.size = fh.seek(0, 2) + fh.seek(0) + + def _fail(self, msg: str) -> PangenomeMetadataError: + return PangenomeMetadataError(f"{self.path}: {msg}") + + def _read(self, nbytes: int) -> bytes: + data = self.fh.read(nbytes) + if len(data) != nbytes: + raise self._fail("the file ends before its contents do") + return data + + def _checked(self, nbytes: int) -> int: + """Reject a length that cannot fit in the file""" + if nbytes < 0 or nbytes > self.size: + raise self._fail( + f"implausible length {nbytes} in the file structure; the " + "file may be truncated or may not be a pangenome file" + ) + return nbytes + + # --- simple-sds primitives ------------------------------------------- + def u64(self) -> int: + return int(struct.unpack(" Tuple[int, int]: + low, high = struct.unpack(" None: + self.fh.seek(self._checked(nbytes), 1) + + def vec_u8(self) -> bytes: + count = self._checked(self.u64()) + data = self._read(count) + self._read((-count) % WORD) # padding to a whole word + return data + + def skip_vec_u8(self) -> None: + count = self.u64() + self.skip(count + (-count) % WORD) + + def vec_u64_words(self) -> bytes: + return self._read(self._checked(self.u64() * WORD)) + + def skip_vec_u64(self) -> None: + self.skip(self.u64() * WORD) + + def skip_option(self) -> None: + """Skip an `Option<...>` stored as a word count plus its words""" + self.skip(self.u64() * WORD) + + def raw_vector(self) -> Tuple[int, bytes]: + nbits = self.u64() + return nbits, self.vec_u64_words() + + def skip_raw_vector(self) -> None: + self.u64() + self.skip_vec_u64() + + def int_vector(self) -> List[int]: + count = self._checked(self.u64()) + width = self.u64() + if not 0 < width <= 64: + raise self._fail(f"unsupported integer width {width}") + nbits, words = self.raw_vector() + if nbits != count * width: + raise self._fail("inconsistent integer vector") + return unpack_bits(words, count, width) + + def skip_int_vector(self) -> None: + self.u64() + self.u64() + self.skip_raw_vector() + + def bit_vector(self) -> Tuple[int, bytes]: + self.u64() # the number of set bits + nbits, words = self.raw_vector() + self.skip_option() # rank support + self.skip_option() # select support + self.skip_option() # select-zero support + return nbits, words + + def skip_bit_vector(self) -> None: + self.u64() + self.skip_raw_vector() + self.skip_option() + self.skip_option() + self.skip_option() + + def sparse_vector_ones(self) -> List[int]: + """The positions of the set bits of an Elias-Fano sparse vector""" + self.u64() # universe size + _nbits, words = self.bit_vector() + count = self._checked(self.u64()) + width = self.u64() + lbits, lwords = self.raw_vector() + if lbits != count * width: + raise self._fail("inconsistent sparse vector") + low = unpack_bits(lwords, count, width) + + # The high part is unary-coded: the i-th set bit at position `p` + # contributes the bucket `p - i`. + result: List[int] = [] + index = 0 + for offset in range(0, len(words), WORD): + end = offset + WORD + word = int.from_bytes(words[offset:end], "little") + base = offset * 8 + while word: + lsb = word & -word + position = base + lsb.bit_length() - 1 + word ^= lsb + bucket = position - index + result.append((bucket << width) | low[index]) + index += 1 + if index == count: + return result + return result + + def skip_sparse_vector(self) -> None: + self.u64() + self.skip_bit_vector() + self.skip_int_vector() + + # --- gbwt-rs structures ---------------------------------------------- + def string_array(self) -> List[str]: + """A `StringArray`: offsets, a byte alphabet and packed indices""" + starts = self.sparse_vector_ones() + alphabet = self.vec_u8() + packed = self.int_vector() + try: + data = bytes(alphabet[index] for index in packed) + except IndexError as err: + raise self._fail( + "a string array indexes its alphabet out of range" + ) from err + starts.append(len(data)) + return [ + data[begin:end].decode("utf-8", "replace") + for begin, end in zip(starts, starts[1:]) + ] + + def tags(self) -> Dict[str, str]: + """A `Tags`: a string array of alternating keys and values""" + array = self.string_array() + return {array[2 * i]: array[2 * i + 1] for i in range(len(array) // 2)} + + def dictionary(self) -> List[str]: + """A `Dictionary`: a string array plus a sorted-id index""" + strings = self.string_array() + self.skip_int_vector() # sorted_ids + return strings + + +@dataclass +class GBZMetadata: + """The headers, tags and path names of a GBZ file""" + + gbz_version: int + gbz_tags: Dict[str, str] + gbwt_version: int + gbwt_tags: Dict[str, str] + sample_names: List[str] + contig_names: List[str] + path_names: List[PathName] + sample_count: int = 0 + haplotype_count: int = 0 + contig_count: int = 0 + sequence_count: int = 0 + bytes_read: int = 0 + + @property + def reference_samples(self) -> List[str]: + """The sample names listed by the `reference_samples` GBWT tag""" + return self.reference_samples_tag.split() + + @property + def reference_samples_tag(self) -> str: + """The raw `reference_samples` GBWT tag, `''` when it is absent""" + return self.gbwt_tags.get("reference_samples", "") + + def sample_of(self, path: PathName) -> str: + """The sample name of a path record""" + return self.sample_names[path[0]] + + def contig_of(self, path: PathName) -> str: + """The contig name of a path record""" + return self.contig_names[path[1]] + + def full_path_name(self, path: PathName) -> str: + """The `##` name of a path record""" + return f"{self.sample_of(path)}#{path[2]}#{self.contig_of(path)}" + + +@dataclass +class HaplHeader: + """The header of a `vg haplotypes` (.hapl) file""" + + version: int + top_level_chains: int + construction_jobs: int + total_subchains: int + total_kmers: int + k: int + tags: Dict[str, str] = field(default_factory=dict) + + +@dataclass +class PangenomeReference: + """The backbone reference sample of a pangenome graph""" + + ref_name: str + contig_prefix: str + contigs: List[str] + reference_samples: List[str] + n_paths: int + + +@dataclass +class _SampleSummary: + """What the path names say about one sample""" + + n_paths: int + phases: List[int] + contigs: List[str] + fragmented: bool + + +def read_gbz_metadata(path: PathLike) -> GBZMetadata: + """Read the tags and path names of a GBZ file. + + Only the file's headers, tags and metadata block are read; the BWT and + the graph sequences are skipped with seeks. + + Raises: + PangenomeMetadataError: the file is not a GBZ, is truncated, or + carries no GBWT metadata. + """ + try: + with open(path, "rb") as fh: + return _read_gbz_metadata(fh, path) + except PangenomeMetadataError: + raise + except (struct.error, IndexError, ValueError, EOFError) as err: + raise PangenomeMetadataError( + f"{path}: the file could not be parsed as a GBZ pangenome " + f"({err})" + ) from err + + +def _read_gbz_metadata(fh: BinaryIO, path: PathLike) -> GBZMetadata: + reader = _Reader(fh, path) + + tag, gbz_version = reader.u32x2() + if tag != GBZ_TAG: + raise PangenomeMetadataError(f"{path}: not a GBZ pangenome file") + reader.u64() # GBZ flags + gbz_tags = reader.tags() + + tag, gbwt_version = reader.u32x2() + if tag != GBWT_TAG: + raise PangenomeMetadataError( + f"{path}: the GBWT header was not found where the GBZ format " + "puts it" + ) + sequence_count = reader.u64() + reader.u64() # size + reader.u64() # offset + reader.u64() # alphabet size + reader.u64() # GBWT flags + gbwt_tags = reader.tags() + + # The BWT body: a sparse-vector index and a (possibly compressed) + # Vec. Both are skipped without reading their contents. + reader.skip_sparse_vector() + reader.skip_vec_u8() + reader.skip_vec_u64() # the document array samples + + # Option + if reader.u64() == 0: + raise PangenomeMetadataError( + f"{path}: the GBWT carries no metadata, so its path names are " + "unknown" + ) + tag, _meta_version = reader.u32x2() + if tag != METADATA_TAG: + raise PangenomeMetadataError( + f"{path}: the GBWT metadata header was not found" + ) + sample_count = reader.u64() + haplotype_count = reader.u64() + contig_count = reader.u64() + reader.u64() # metadata flags + n_paths = reader._checked(reader.u64()) + raw = reader._read(n_paths * 16) + path_names: List[PathName] = [ + (rec[0], rec[1], rec[2], rec[3]) + for rec in struct.iter_unpack(" HaplHeader: + """Read the header of a `vg haplotypes` (.hapl) file. + + Raises: + PangenomeMetadataError: the file is not a `vg haplotypes` file or + is too short to hold a header. + """ + try: + with open(path, "rb") as fh: + reader = _Reader(fh, path) + magic, version = reader.u32x2() + if magic != HAPL_MAGIC: + raise PangenomeMetadataError( + f"{path}: not a vg haplotypes file" + ) + return HaplHeader( + version=version, + top_level_chains=reader.u64(), + construction_jobs=reader.u64(), + total_subchains=reader.u64(), + total_kmers=reader.u64(), + k=reader.u64(), + # Tags were added to the header in version 6 + tags=reader.tags() if version >= 6 else {}, + ) + except PangenomeMetadataError: + raise + except (struct.error, IndexError, ValueError, EOFError) as err: + raise PangenomeMetadataError( + f"{path}: the file could not be parsed as a vg haplotypes " + f"file ({err})" + ) from err + + +def _summarize_samples(meta: GBZMetadata) -> Dict[str, _SampleSummary]: + """Group the path names of a GBZ by sample. + + A sample is "fragmented" when any of its paths carries a non-zero + fragment, or when one (contig, phase) pair has more than one path. + """ + by_sample: Dict[int, List[PathName]] = defaultdict(list) + for path in meta.path_names: + by_sample[path[0]].append(path) + + summaries: Dict[str, _SampleSummary] = {} + for sample_id, paths in by_sample.items(): + if sample_id >= len(meta.sample_names): + continue + per_contig: Dict[Tuple[int, int], int] = defaultdict(int) + for path in paths: + per_contig[(path[1], path[2])] += 1 + summaries[meta.sample_names[sample_id]] = _SampleSummary( + n_paths=len(paths), + phases=sorted({path[2] for path in paths}), + contigs=sorted({meta.contig_of(path) for path in paths}), + fragmented=( + any(path[3] != 0 for path in paths) + or any(count > 1 for count in per_contig.values()) + ), + ) + return summaries + + +def detect_pangenome_reference(gbz_path: PathLike) -> PangenomeReference: + """Detect the backbone reference sample of a pangenome graph. + + The GBWT `reference_samples` tag lists the graph's reference samples, + for example `"GRCh38 CHM13"`. Exactly one of them is the backbone: its + paths are unfragmented, one per contig. The others were chopped into + fragments when the graph was built. + + Raises: + PangenomeMetadataError: the file is not a GBZ, has no metadata, + names no reference sample, or names no single unfragmented + one. + """ + meta = read_gbz_metadata(gbz_path) + ref_samples = meta.reference_samples + if not ref_samples: + raise PangenomeMetadataError( + f"{gbz_path}: the GBWT 'reference_samples' tag is missing or " + "empty, so the graph names no reference sample" + ) + + summaries = _summarize_samples(meta) + present = [name for name in ref_samples if name in summaries] + if not present: + raise PangenomeMetadataError( + f"{gbz_path}: no sample named by the 'reference_samples' tag " + f"({' '.join(ref_samples)}) has any path in the graph" + ) + + backbone = [name for name in present if not summaries[name].fragmented] + if not backbone: + raise PangenomeMetadataError( + f"{gbz_path}: every reference sample ({', '.join(present)}) is " + "split into fragments, so the backbone reference cannot be " + "identified" + ) + if len(backbone) > 1: + raise PangenomeMetadataError( + f"{gbz_path}: more than one reference sample is unfragmented " + f"({', '.join(backbone)}), so the backbone reference is " + "ambiguous" + ) + + name = backbone[0] + summary = summaries[name] + if len(summary.phases) != 1: + phases = ", ".join(str(phase) for phase in summary.phases) + raise PangenomeMetadataError( + f"{gbz_path}: the reference sample '{name}' has paths on more " + f"than one haplotype ({phases}), so its contig prefix is " + "ambiguous" + ) + + reference = PangenomeReference( + ref_name=name, + contig_prefix=f"{name}#{summary.phases[0]}#", + contigs=summary.contigs, + reference_samples=ref_samples, + n_paths=summary.n_paths, + ) + logger.debug( + "Detected pangenome reference '%s' (contig prefix '%s', %d paths) " + "in '%s'", + reference.ref_name, + reference.contig_prefix, + reference.n_paths, + gbz_path, + ) + return reference + + +def check_sample_gbz( + gbz_path: PathLike, + ref_name: str, + contig_prefix: str, +) -> None: + """Check that a sampled pangenome kept its reference paths. + + `vg haplotypes --set-reference NAME` only warns when `NAME` is not a + reference sample of the input graph. The sampled graph it writes then + has an empty `reference_samples` tag and no reference paths at all, + which every downstream step turns into empty output. + + Raises: + PangenomeMetadataError: the file could not be read, its + `reference_samples` tag is not exactly `ref_name`, or the + reference sample has no path under `contig_prefix`. + """ + meta = read_gbz_metadata(gbz_path) + observed = meta.reference_samples + if observed != [ref_name]: + raise PangenomeMetadataError( + f"{gbz_path}: the sampled pangenome has " + f"reference_samples='{meta.reference_samples_tag}', expected " + f"'{ref_name}'. The pangenome reference name may be wrong: " + "`vg haplotypes --set-reference` drops every reference path " + "when the name it is given is not a reference sample of the " + "input graph." + ) + + reference_paths = [ + path for path in meta.path_names if meta.sample_of(path) == ref_name + ] + if not any( + meta.full_path_name(path).startswith(contig_prefix) + for path in reference_paths + ): + raise PangenomeMetadataError( + f"{gbz_path}: none of the {len(reference_paths)} paths of the " + f"reference sample '{ref_name}' starts with the contig prefix " + f"'{contig_prefix}'. The pangenome reference name or contig " + "prefix may be wrong." + ) + logger.info( + "The sampled pangenome '%s' has %d reference paths for '%s'", + gbz_path, + len(reference_paths), + ref_name, + ) + + +def _gfa_header_reference_samples(header: str) -> List[str]: + """The samples named by the `RS:Z:` tag of a GFA header line""" + for gfa_field in header.rstrip("\n").split("\t")[1:]: + if gfa_field.startswith("RS:Z:"): + return gfa_field.removeprefix("RS:Z:").split() + return [] + + +def check_gfa( + gfa_path: PathLike, + ref_name: str, + contig_prefix: str, + max_lines: int = DEFAULT_MAX_GFA_LINES, +) -> None: + """Check that a GFA carries the rGFA tags of the reference. + + `vg convert -f -Q NAME` writes an rGFA whose `S` lines carry + `SN:Z:##`, and lists the graph's remaining + reference samples in the header's `RS:Z:` tag. When `NAME` does not + match a reference sample, `vg` writes no `SN:Z:` tags at all and + leaves `NAME` in `RS:Z:`; `pgutil lift` then leaves every read + unmapped. + + Only the first `max_lines` lines are read: a whole-genome GFA runs to + many GB, and its first segment lines already carry the tags. + + Raises: + PangenomeMetadataError: the header lists `ref_name` among its + reference samples, or no segment line in the head of the file + carries an `SN:Z:` tag. + """ + sn_tag = f"SN:Z:{contig_prefix}" + header_checked = False + lines_read = 0 + with open(gfa_path, "r") as fh: + for line in fh: + if lines_read >= max_lines: + break + lines_read += 1 + if line.startswith("H\t") and not header_checked: + header_checked = True + rs_samples = _gfa_header_reference_samples(line) + if ref_name in rs_samples: + raise PangenomeMetadataError( + f"{gfa_path}: the GFA header still lists " + f"'{ref_name}' among its reference samples " + f"(RS:Z:{' '.join(rs_samples)}), so the rGFA tags " + "were written for a different reference. The " + "pangenome reference name may be wrong." + ) + elif line.startswith("S\t") and sn_tag in line: + logger.info( + "The GFA '%s' carries '%s' rGFA tags", gfa_path, sn_tag + ) + return + + raise PangenomeMetadataError( + f"{gfa_path}: no segment line in the first {lines_read} lines " + f"carries an '{sn_tag}' tag. `vg convert -Q` writes no rGFA tags " + "when its reference name does not match a reference sample of the " + "graph, and `pgutil lift` then leaves every read unmapped." + ) + + +def _build_parser() -> argparse.ArgumentParser: + """The argument parser of the `pangenome_meta` script""" + parser = argparse.ArgumentParser( + prog="python -m sentieon_cli.pangenome_meta", + description="Inspect and check pangenome graph metadata", + ) + subparsers = parser.add_subparsers(dest="command", required=True) + + detect = subparsers.add_parser( + "detect", + help="Print the reference sample detected in a pangenome GBZ", + ) + detect.add_argument( + "--gbz", required=True, type=pathlib.Path, help="The GBZ file" + ) + + check_gbz = subparsers.add_parser( + "check-gbz", + help="Check that a sampled GBZ kept its reference paths", + ) + check_gbz.add_argument( + "--gbz", + required=True, + type=pathlib.Path, + help="The sampled GBZ file", + ) + check_gbz.add_argument( + "--reference_name", + required=True, + help="The expected pangenome reference sample name", + ) + check_gbz.add_argument( + "--contig_prefix", + required=True, + help="The expected pangenome contig prefix", + ) + + check_gfa_parser = subparsers.add_parser( + "check-gfa", + help="Check that a GFA carries the reference's rGFA tags", + ) + check_gfa_parser.add_argument( + "--gfa", required=True, type=pathlib.Path, help="The GFA file" + ) + check_gfa_parser.add_argument( + "--reference_name", + required=True, + help="The expected pangenome reference sample name", + ) + check_gfa_parser.add_argument( + "--contig_prefix", + required=True, + help="The expected pangenome contig prefix", + ) + check_gfa_parser.add_argument( + "--max_lines", + type=int, + default=DEFAULT_MAX_GFA_LINES, + help=( + "Stop looking for the rGFA tags after this many lines " + f"(default: {DEFAULT_MAX_GFA_LINES})" + ), + ) + return parser + + +def main(argv: Optional[Sequence[str]] = None) -> int: + """Run a pangenome metadata check; return 0 on success and 1 on + failure""" + args = _build_parser().parse_args(argv) + try: + if args.command == "detect": + reference = detect_pangenome_reference(args.gbz) + print(f"pangenome_ref_name\t{reference.ref_name}") + print(f"pangenome_contig_prefix\t{reference.contig_prefix}") + print( + "reference_samples\t" + " ".join(reference.reference_samples) + ) + print(f"n_paths\t{reference.n_paths}") + print(f"n_contigs\t{len(reference.contigs)}") + print("contigs\t" + ",".join(reference.contigs)) + elif args.command == "check-gbz": + check_sample_gbz(args.gbz, args.reference_name, args.contig_prefix) + else: + check_gfa( + args.gfa, + args.reference_name, + args.contig_prefix, + max_lines=args.max_lines, + ) + except (PangenomeMetadataError, OSError) as err: + logger.error("%s", err) + return 1 + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/sentieon_cli/sentieon_pangenome.py b/sentieon_cli/sentieon_pangenome.py index 65c2375..27e258d 100644 --- a/sentieon_cli/sentieon_pangenome.py +++ b/sentieon_cli/sentieon_pangenome.py @@ -54,7 +54,6 @@ vcf_id, ) from .shard import ( - GRCH38_CONTIGS, determine_shards_from_fai, parse_fai, vcf_contigs, @@ -113,22 +112,6 @@ class SentieonPangenome(BasePangenome): "help": "Generate a gVCF output file.", "action": "store_true", }, - "pangenome_ref_name": { - "default": "GRCh38", - "help": ( - "Reference name in the pangenome (GRCh38). The " - "'extract..model' member of the model " - "bundle is preferred; if it is absent and the reference " - "name is 'GRCh38', the pipeline falls back to " - "'extract.model'." - ), - }, - "pangenome_contig_prefix": { - "default": "GRCh38#0#", - "help": ( - "Prefix to strip from pangenome contig names (GRCh38#0#)" - ), - }, "skip_metrics": { "help": "Skip metrics collection and multiQC", "action": "store_true", @@ -138,14 +121,6 @@ class SentieonPangenome(BasePangenome): "action": "store_true", }, # Hidden arguments - "skip_contig_checks": { - "help": argparse.SUPPRESS, - "action": "store_true", - }, - "skip_pangenome_name_checks": { - "help": argparse.SUPPRESS, - "action": "store_true", - }, "skip_pop_vcf_id_check": { "help": argparse.SUPPRESS, "action": "store_true", @@ -171,13 +146,9 @@ def __init__(self) -> None: self.bed: Optional[pathlib.Path] = None self.call_svs = False self.gvcf = False - self.pangenome_ref_name = "GRCh38" self.extract_model_name = "extract.model" - self.pangenome_contig_prefix = "GRCh38#0#" self.skip_metrics = False self.skip_multiqc = False - self.skip_contig_checks: bool = False - self.skip_pangenome_name_checks: bool = False self.skip_pop_vcf_id_check: bool = False self.skip_model_apply = False self.skip_small_variants = False @@ -201,7 +172,10 @@ def validate(self) -> None: self.shards = determine_shards_from_fai( self.fai_data, 10 * 1000 * 1000 ) - self.pop_vcf_contigs: Dict[str, Optional[int]] = {} + # Before `validate_bundle`, which picks the bundle's + # `extract..model` member by the resolved name + self.resolve_pangenome_reference() + self.pop_vcf_contigs = {} if self.pop_vcf: self.pop_vcf_contigs = vcf_contigs(self.pop_vcf, self.dry_run) self.logger.debug("VCF contigs are: %s", self.pop_vcf_contigs) @@ -253,51 +227,7 @@ def validate(self) -> None: "across decoy and unplaced contigs." ) - if not self.skip_pangenome_name_checks: - if not str(self.gbz).endswith("grch38.gbz"): - self.logger.error( - "The `--gbz` file does not have the expected suffix. " - "Check that you are using a GRCh38 pangenome." - ) - sys.exit(2) - - if not str(self.hapl).endswith("grch38.hapl"): - self.logger.error( - "The `--hapl` file does not have the expected suffix. " - "Check that you re using a GRCh38 pangenome." - ) - sys.exit(2) - - if not self.skip_contig_checks: - # Check the fai file contigs - mismatch_contigs: Set[str] = set() - for ctg, length in GRCH38_CONTIGS.items(): - d = self.fai_data.get(ctg, {}) - fai_length = d.get("length", -1) - if length != fai_length: - mismatch_contigs.add(ctg) - if mismatch_contigs: - mismatch_contigs_s = ", ".join(mismatch_contigs) - self.logger.error( - "Reference contigs with unexpected lengths: %s", - mismatch_contigs_s, - ) - sys.exit(2) - - # Check the pop VCF file contigs - if not self.dry_run: - mismatch_contigs = set() - for ctg, length in GRCH38_CONTIGS.items(): - vcf_length = self.pop_vcf_contigs.get(ctg, -1) - if length != vcf_length: - mismatch_contigs.add(ctg) - if mismatch_contigs: - mismatch_contigs_s = ", ".join(mismatch_contigs) - self.logger.error( - "Pop VCF contigs with unexpected lengths: %s", - mismatch_contigs_s, - ) - sys.exit(2) + self.validate_grch38_contigs() def validate_segdup(self) -> None: if self.segdup_caller is None: @@ -589,11 +519,21 @@ def build_dag(self) -> DAG: haplotypes_job = self.build_haplotypes_job(sample_pangenome, kmer_file) dag.add_job(haplotypes_job, haplotype_dependencies) + # Confirm the sampled pangenome kept its reference paths before + # anything consumes it + check_gbz_job = self.build_check_gbz_job(sample_pangenome) + dag.add_job(check_gbz_job, {haplotypes_job}) + # convert the sample pangenome gfa_job = self.build_gfa_job(sample_gfa, sample_pangenome) fasta_job = self.build_fasta_job(sample_fasta, sample_pangenome) - dag.add_job(gfa_job, {haplotypes_job}) - dag.add_job(fasta_job, {haplotypes_job}) + dag.add_job(gfa_job, {check_gbz_job}) + dag.add_job(fasta_job, {check_gbz_job}) + + # Confirm the GFA carries the reference's rGFA tags, without which + # `pgutil lift` leaves every read unmapped + check_gfa_job = self.build_check_gfa_job(sample_gfa) + dag.add_job(check_gfa_job, {gfa_job}) # minimap2 alignment of the extracted fastq dnascope_dependencies = set() @@ -603,7 +543,7 @@ def build_dag(self) -> DAG: sample_fasta, sample_gfa, ) - dag.add_job(mm2_job, mm2_dependencies | {gfa_job, fasta_job}) + dag.add_job(mm2_job, mm2_dependencies | {check_gfa_job, fasta_job}) dnascope_dependencies.add(mm2_job) # With fastq input, perform dedup and metrics @@ -711,7 +651,7 @@ def build_dag(self) -> DAG: PangenomeSV( sv_vcf, gfa_file=gfa_file, - prefix=self.pangenome_contig_prefix, + prefix=self.contig_prefix(), ) ) call = DNAscopeStage( @@ -827,7 +767,7 @@ def build_haplotypes_job( "--include-reference", "--diploid-sampling", "--set-reference", - self.pangenome_ref_name, + self.ref_name(), ], ), "vg-haplotypes", @@ -845,7 +785,7 @@ def build_gfa_job( output_gfa, input_gbz, threads=self.cores, - reference_name=self.pangenome_ref_name, + reference_name=self.ref_name(), ), "vg-convert-gfa", 0, @@ -899,7 +839,7 @@ def build_minimap2_lift_job( "@RG\\t" + "\\t".join([f"{x[0]}:{x[1]}" for x in rg2.items()]), mm2_model, threads=self.cores, - lift_prefix=self.pangenome_contig_prefix, + lift_prefix=self.contig_prefix(), ), "mm2-lift", self.cores, diff --git a/tests/data/pangenome/generate.sh b/tests/data/pangenome/generate.sh new file mode 100755 index 0000000..1d51556 --- /dev/null +++ b/tests/data/pangenome/generate.sh @@ -0,0 +1,40 @@ +#!/bin/bash +# Regenerate the binary pangenome fixtures in this directory. +# +# The `.gfa` files are the checked-in sources; `vg` (tested with v1.72.0) +# turns them into the GBZ graphs and the `vg haplotypes` index that +# `tests/unit/test_pangenome_meta.py` parses. Everything here is a toy +# graph of 8 segments and two 30 bp contigs, so the outputs stay a few KB. +# +# The graphs differ only in their walks and in the `RS:Z:` header tag, +# which `vg` turns into the GBWT `reference_samples` tag: +# +# tiny.gfa GRCh38 is the unfragmented backbone, CHM13 is +# split into two fragments on chr1 +# tiny_chm13.gfa the same graph with the two roles swapped +# tiny_no_backbone.gfa both reference samples are fragmented +# tiny_no_ref_tag.gfa no `RS:Z:` tag, so no reference sample at all +# +# `tiny_sampled.gbz` stands in for the graph `vg haplotypes +# --set-reference GRCh38` writes: one reference sample, as +# `check_sample_gbz` requires. +# +# The `rgfa_*.gfa` files are hand-written heads of a `vg convert -f -Q` +# rGFA and need no regeneration. +set -euo pipefail +cd "$(dirname "$0")" + +for name in tiny tiny_chm13 tiny_no_backbone tiny_no_ref_tag; do + vg gbwt -G "${name}.gfa" --gbz-format -g "${name}.gbz" +done + +vg gbwt -Z tiny.gbz --set-reference GRCh38 -g tiny_sampled.gbz + +vg index -j tiny.dist tiny.gbz +vg gbwt -Z tiny.gbz -r tiny.ri +vg haplotypes -d tiny.dist -r tiny.ri -H tiny.hapl tiny.gbz +rm -f tiny.dist tiny.ri + +# What the fixtures should contain +vg gbwt -Z tiny.gbz --tags +vg paths -x tiny.gbz -L diff --git a/tests/data/pangenome/rgfa_no_header.gfa b/tests/data/pangenome/rgfa_no_header.gfa new file mode 100644 index 0000000..afa1dce --- /dev/null +++ b/tests/data/pangenome/rgfa_no_header.gfa @@ -0,0 +1,3 @@ +S 1 ACGTACGTAC SN:Z:GRCh38#0#chr1 SO:i:0 SR:i:0 +S 2 GGGGGGGGGG SN:Z:GRCh38#0#chr1 SO:i:10 SR:i:0 +L 1 + 2 + 0M diff --git a/tests/data/pangenome/rgfa_no_tags.gfa b/tests/data/pangenome/rgfa_no_tags.gfa new file mode 100644 index 0000000..c4dcc80 --- /dev/null +++ b/tests/data/pangenome/rgfa_no_tags.gfa @@ -0,0 +1,4 @@ +H VN:Z:1.1 +S 1 ACGTACGTAC +S 2 GGGGGGGGGG +L 1 + 2 + 0M diff --git a/tests/data/pangenome/rgfa_ok.gfa b/tests/data/pangenome/rgfa_ok.gfa new file mode 100644 index 0000000..c6a8256 --- /dev/null +++ b/tests/data/pangenome/rgfa_ok.gfa @@ -0,0 +1,6 @@ +H VN:Z:1.1 RS:Z:CHM13 +S 1 ACGTACGTAC SN:Z:GRCh38#0#chr1 SO:i:0 SR:i:0 +S 2 GGGGGGGGGG SN:Z:GRCh38#0#chr1 SO:i:10 SR:i:0 +S 3 TTTTTTTTTT SR:i:1 +L 1 + 2 + 0M +L 1 + 3 + 0M diff --git a/tests/data/pangenome/rgfa_wrong_reference.gfa b/tests/data/pangenome/rgfa_wrong_reference.gfa new file mode 100644 index 0000000..46a1708 --- /dev/null +++ b/tests/data/pangenome/rgfa_wrong_reference.gfa @@ -0,0 +1,6 @@ +H VN:Z:1.1 RS:Z:CHM13 GRCh38 +S 1 ACGTACGTAC +S 2 GGGGGGGGGG +S 3 TTTTTTTTTT +L 1 + 2 + 0M +L 1 + 3 + 0M diff --git a/tests/data/pangenome/tiny.gbz b/tests/data/pangenome/tiny.gbz new file mode 100644 index 0000000000000000000000000000000000000000..d24bff4ddc989cf03dfaccdca0397c3aef83487b GIT binary patch literal 2184 zcmaJ?QAkr!7(VCTJGY$HGK*nR(}H}dIkh%}z}u=VP>HdD1!d`M)3v3pVH83yLQjGe zicyiE4+UWbSrR=3y;i=&^dWj!k@e6+_*B%ncmKEKxOFe?yZ?W_@B9CA&$-9Ik<%`P z{kuBoQ6tA9 zzML6mUBt5xE>j~mid){i{GK>A9@|)6G4XscFp<6U3MRUJBHDL=lWcG~w@ zY3sc+8P2}=wi?e}TN&1)_X4k;nk{-0-h|=3jr>A{m~R!uR>D1mn4kQ99?*CxfKDYs z1tIDMsnAZ!yWE~8pFa?64|R$6v+nh<2IArU{Gb?hvWF+{COf`HJej4LDl5;qJU?>_ z3n!jWWV7tlN5)2oznRcOh<uscvLi{!eHke1>YgkcJb#BM0=jB zRDvP;*h~GD5u$(0cV?6ZO2BqWUr90a?O2Esdejv?o0q))&9tFV2m2EH4t1877U%z0 z3OexRxc^LkQIUsNiuyu-wT}ih9-Lb29=WNnh$c?b_|>!Cgu%f+ljnu+#ZG(~&=>M+ z<7$yaedG}X@%>=j@R6U_#)i&_kwd)J&lp3Da|`p6*J&IBLR0L*FzU%M{P2E^V=(`n z$4$E}WC~py#prKh)}L{v%Lc{wWM1H59OonIebqlk1e;%p)2@Gbui7pk$Qs5 literal 0 HcmV?d00001 diff --git a/tests/data/pangenome/tiny.gfa b/tests/data/pangenome/tiny.gfa new file mode 100644 index 0000000..519c6e5 --- /dev/null +++ b/tests/data/pangenome/tiny.gfa @@ -0,0 +1,28 @@ +H VN:Z:1.1 RS:Z:GRCh38 CHM13 +S s1 ACGTACGTAC +S s2 GGGGGGGGGG +S s3 TTTTTTTTTT +S s4 CCCCCCCCCC +S s5 AAAAAAAAAA +S s6 ACACACACAC +S s7 TGTGTGTGTG +S s8 GATTACAGAT +L s1 + s2 + 0M +L s1 + s3 + 0M +L s2 + s4 + 0M +L s3 + s4 + 0M +L s5 + s6 + 0M +L s5 + s7 + 0M +L s6 + s8 + 0M +L s7 + s8 + 0M +W GRCh38 0 chr1 0 30 >s1>s2>s4 +W GRCh38 0 chr2 0 30 >s5>s6>s8 +W CHM13 0 chr1 0 20 >s1>s3 +W CHM13 0 chr1 20 30 >s4 +W CHM13 0 chr2 0 30 >s5>s7>s8 +W HG01 1 chr1 0 30 >s1>s3>s4 +W HG01 2 chr1 0 30 >s1>s2>s4 +W HG01 1 chr2 0 30 >s5>s6>s8 +W HG01 2 chr2 0 30 >s5>s7>s8 +W HG02 1 chr1 0 30 >s1>s2>s4 +W HG02 2 chr2 0 30 >s5>s6>s8 diff --git a/tests/data/pangenome/tiny.hapl b/tests/data/pangenome/tiny.hapl new file mode 100644 index 0000000000000000000000000000000000000000..7bf79bb9bda60fdc654433fedcf1845354a8811d GIT binary patch literal 760 zcmeZp4DexNU|?VZVh}*5A)*WnvQP#Gl!mDT=>-8s2*n@>rEyU(eK7U7_##jP3=EBo zO-#+qEi4n0l2cOC(sT0)KvM2dc~&T$3^fqu9$}}Ii;`7%!_P3j7P%fFy{;jc$>Q<~)Qu2D zVpJs5g+W+>l|@%Yx0N@syyzy2tcxzftD?TM^FC{QZk>aB&ilU4^S z9rURZt|OEPUlNuA=ztDEKlt8oY$w0l<85ec+OyZ^Z$8-*i}#*NB>VbL52Viw4v&mx z##lG;OoU6+h=uIt*U!JF4qi@Fl-s7iT;K6VjUQ-Lt1iCnOS>w=EvE8A)ipEz&vIMu zjj2fX>9^J7^7)l9Eq*Kb;!&p{QMvrb+Y@WZVq&Oje0ZBs_V?2*;m@n<%NaA zPp)ONtm`9Vx30oo!^~da%8bH{xtor8TQjtl3-N_-K}7yN~G zC3m7e{IG#IKZxsj#P32|*XcHVunY6)G1$1bFg|geE+N1*#U6xFPqZsYg!$30LF1e! zbCFFHhpvrm^p|nQx>`ofL;O7D9`F#y{fPRXNRM>}4>owX&xpPTpaWmDAwXVlLo7MW zo#+Sa<|HEWRLvMSkC?BNJs0<-^f8Nhu|MI*{f#=&d+4w@LRYF!hy7hnh&sC|P5r$4 zH|yLbvD;c87>NY5KqMH{A`va1aVTRKPE2|pn64@l+QY;oFS4F=sE7TQ?=xDwKX{=7 vU$nu)dg8zA2nA3N`=*ur`FzfLf2cEg7JB3t@4EqPU=$P(x_n-t4u<&`Z@z+2 literal 0 HcmV?d00001 diff --git a/tests/data/pangenome/tiny_chm13.gfa b/tests/data/pangenome/tiny_chm13.gfa new file mode 100644 index 0000000..40451f7 --- /dev/null +++ b/tests/data/pangenome/tiny_chm13.gfa @@ -0,0 +1,28 @@ +H VN:Z:1.1 RS:Z:CHM13 GRCh38 +S s1 ACGTACGTAC +S s2 GGGGGGGGGG +S s3 TTTTTTTTTT +S s4 CCCCCCCCCC +S s5 AAAAAAAAAA +S s6 ACACACACAC +S s7 TGTGTGTGTG +S s8 GATTACAGAT +L s1 + s2 + 0M +L s1 + s3 + 0M +L s2 + s4 + 0M +L s3 + s4 + 0M +L s5 + s6 + 0M +L s5 + s7 + 0M +L s6 + s8 + 0M +L s7 + s8 + 0M +W CHM13 0 chr1 0 30 >s1>s2>s4 +W CHM13 0 chr2 0 30 >s5>s6>s8 +W GRCh38 0 chr1 0 20 >s1>s3 +W GRCh38 0 chr1 20 30 >s4 +W GRCh38 0 chr2 0 30 >s5>s7>s8 +W HG01 1 chr1 0 30 >s1>s3>s4 +W HG01 2 chr1 0 30 >s1>s2>s4 +W HG01 1 chr2 0 30 >s5>s6>s8 +W HG01 2 chr2 0 30 >s5>s7>s8 +W HG02 1 chr1 0 30 >s1>s2>s4 +W HG02 2 chr2 0 30 >s5>s6>s8 diff --git a/tests/data/pangenome/tiny_no_backbone.gbz b/tests/data/pangenome/tiny_no_backbone.gbz new file mode 100644 index 0000000000000000000000000000000000000000..453f4602592b561b7770b6396737d8bfbacd850f GIT binary patch literal 2240 zcmaJ?O-NKx6u$SpH$VPNjXw*Ks0pJiFnwyy7zDni90QY(H(@O_bIh@E%rvo|O+-PP z1f!IqBB2&p&<_nI(JE*Y^s|b|MYOOeBU(hepnKn(&!%hBdko(_=X~Ee=iYbTo5x$)=OUb=M!aM^xchO51$X1 zKjvp=cRZR*r`ewOjExh&h_HeX{aC4J5^DpGVaQQE~du4{%?DEXP(j>6A{ zqJU}KFQ^5JbRQQqO=uXaq<--nu^$m~t~=+jB8LKN$%ej!Y^YK;R*F#3iH4jD^9I%F zlFSNqFsD_-2XF3S)!YCc0>!ccI`CC{!^Ih=1^6*frzi$VdMFL8%C$U+}VR8Ao3DVFS@O;xbS6I@?CQ?6}M$ZTMhk z=a(_q_}*dv>ia!`fGqrv3t`k#ZR9H;t^9muz7KK7gZYE5m3ZiH@@9CnY{Y%oO%Cu7 z$9EO=Kad{t03K}cF#kkf0?>i4+7KYGS=$pI7Ebho=U7HW_|^9?IzxmU09~#={k}Jj5_L9Gp2}49FV<1T*>x)%vF)%Gwj&YCwym%wz{5`N z8?U+Nx}cx2?)Qz0BFjmKdYCtyM>LwFK}82X^1y}~&lCU7BNRYA%uN&dq4$s1>s2 +W GRCh38 0 chr1 20 30 >s4 +W GRCh38 0 chr2 0 20 >s5>s6 +W GRCh38 0 chr2 20 30 >s8 +W CHM13 0 chr1 0 20 >s1>s3 +W CHM13 0 chr1 20 30 >s4 +W CHM13 0 chr2 0 30 >s5>s7>s8 +W HG01 1 chr1 0 30 >s1>s3>s4 +W HG01 2 chr1 0 30 >s1>s2>s4 +W HG01 1 chr2 0 30 >s5>s6>s8 +W HG01 2 chr2 0 30 >s5>s7>s8 +W HG02 1 chr1 0 30 >s1>s2>s4 +W HG02 2 chr2 0 30 >s5>s6>s8 diff --git a/tests/data/pangenome/tiny_no_ref_tag.gbz b/tests/data/pangenome/tiny_no_ref_tag.gbz new file mode 100644 index 0000000000000000000000000000000000000000..a19ffcec2e298de409e0c7cef7226d418ca17358 GIT binary patch literal 2152 zcmah~OK1~O6rFDxn?$2l301JIg57AESWW6GU&k(t8!X8-KhJ3HShxx)#i_4>Pka$8I4pxS(Nd;47ci`1LH zcia5<;@evB^Oe;JEqC8|_0(=LR*ID@?X9GjB8&NMBHYATteJ=7ZzSwykOzGpBX0w} zgfrQEfvre+*jbOwb1m*-w8}o9DPLyqrRA`cBHqEvw=Y#^gFox z9}8>0WadL1k0hvqXg_SpWI>c<_RAhc#Um*KQin1(PZlt$#1b8MaB;@w~ze=8PVsEr2rVLkKNQC_67Z8z7MB};9|CM`2oVn@54cK$wyt* zvpr}V*aV!>O&=8M;N0z^IC#sSmY)Bs6!O65@!&~kQ{jhJnd+dw_W4moOwKN4CU0n~ ztcl|kzkDiIvN$+r{Jiiz_=wLN^aX!CZedQ;haVW$I#bp_ywQ7o%*wMce8Ba3Rty;T z66VM6(+mPu;RXX?)Z;Pqus+r`7uStsU#e)U{qH-zRn^Bs70)R9p6mgA^pE=w^;f7z z>@#@4;Nd(2Yd1#WChYKiIbr5uwvM=iEHJo*Q1s#Y)zkd0m`O_;Ej@ zPWll#Y;2L|)@LyA-cCfFy(H85Ec{dG-y;eSG~F;wUDHj&&`eX)HD;QGt7oRW9;-Ly zDeXyNniawQje0n5wLa60`-2yG;PV(f>?i(s1>s2>s4 +W GRCh38 0 chr2 0 30 >s5>s6>s8 +W CHM13 0 chr1 0 20 >s1>s3 +W CHM13 0 chr1 20 30 >s4 +W CHM13 0 chr2 0 30 >s5>s7>s8 +W HG01 1 chr1 0 30 >s1>s3>s4 +W HG01 2 chr1 0 30 >s1>s2>s4 +W HG01 1 chr2 0 30 >s5>s6>s8 +W HG01 2 chr2 0 30 >s5>s7>s8 +W HG02 1 chr1 0 30 >s1>s2>s4 +W HG02 2 chr2 0 30 >s5>s6>s8 diff --git a/tests/data/pangenome/tiny_sampled.gbz b/tests/data/pangenome/tiny_sampled.gbz new file mode 100644 index 0000000000000000000000000000000000000000..c0e19be0c0eea49a89ebf40486ab527b1fac1709 GIT binary patch literal 2184 zcmaJ?O-NKx6u$SqH{;|_O=%%hm|Bqq=5Ni|%GaxA78)_U0SU^|F-IrOUyMSCKxidM zW*89(x(PxHlAu}zE!vc`$jViV3aqsd6m{Qw=d zQ}n76ZX)D_uL!dNbU?qN?^|{qyUAZwT~k|EfA~l+)NrA#J>GG#Gtt$3si*gH-@xF| z@Ca)qo`-OS8u5|c@cQ|WVcc@x_V2CqDMWsCwG<2guaTV zj(c~kkI%lZCB9r+9Wmnf!!MsWEe4apq@8;Q<%<$xzT3#oAuJ)p{PGq43B2Z@LENs)~@vC7!(rZ*o74_H%^OE(q@muE{}ypQeU z=1;eSQQ~hUtRX}{JUUk%+fL3503GI2uIPWwul>=a4|TkPc^G#dxuhE7LX*Gp(LH#C z&SW-C$QE?&lMUT3&XJ@GAr=C7I)WNyjOMKK>`IK7xUn%IoSMBbz z?xwIaB9_?}ybd>5=AW!pT-h zIiYJJ8~shqnZw0{-hTN#IRiYz)qCYL>9Nn?!A76B?}@$vpaWmEAwb@m+V;eNbfO>Z zTQL!lCotpQJZipd_Ix~;wCAkm#rcFE_c!Xq9-_m{2wk>5{RK;o5Ooewnu@1YzgbHq zXAcZBY+0sZTH&x^S%zsy$YIydk5?~vCbco+ap$-!vW0Z0hx4}HXDoAn@InW^YJ-RU r#DCdA3ZNd&O(Xf&^A%TwA}xt4&?CQk-}T@Cqo9J&t>;zhV48mcb%uUo literal 0 HcmV?d00001 diff --git a/tests/unit/test_hybrid_pangenome.py b/tests/unit/test_hybrid_pangenome.py index 4d35072..fbda6eb 100644 --- a/tests/unit/test_hybrid_pangenome.py +++ b/tests/unit/test_hybrid_pangenome.py @@ -118,6 +118,8 @@ def create_pipeline(self): pipeline.tmp_dir = self.mock_dir # State normally set by validate() + pipeline.pangenome_ref_name = "GRCh38" + pipeline.pangenome_contig_prefix = "GRCh38#0#" pipeline.fai_data = {"chr1": {"length": 1000}} pipeline.shards = [MagicMock()] pipeline.shards[0].contig = "chr1" @@ -515,6 +517,74 @@ def test_pangenomesv_command(self): # SV calling is not restricted to the small-variant BED assert f"--interval {self.mock_bed}" not in cmd_str + def test_pangenome_check_jobs(self): + """The graph checks sit between vg and everything downstream. + + A wrong `--pangenome_ref_name` is silent in `vg`: the sampled + graph loses its reference paths and the GFA loses its rGFA tags, + and `pgutil lift` then leaves every read unmapped. + """ + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + job_names, all_jobs = self._get_all_job_names(dag) + + assert "check-sample-gbz" in job_names + assert "check-sample-gfa" in job_names + + assert self._get_dep_names(dag, all_jobs, "check-sample-gbz") == { + "vg-haplotypes" + } + assert self._get_dep_names(dag, all_jobs, "vg-convert-gfa") == { + "check-sample-gbz" + } + assert self._get_dep_names(dag, all_jobs, "check-sample-gfa") == { + "vg-convert-gfa" + } + # The raw graph update is the only consumer of the converted GFA + assert self._get_dep_names(dag, all_jobs, "graph-update-raw") == { + "check-sample-gfa" + } + + # The checks share the log group of the jobs they guard + haplotypes = self._get_job(all_jobs, "vg-haplotypes") + check_gbz = self._get_job(all_jobs, "check-sample-gbz") + check_gfa = self._get_job(all_jobs, "check-sample-gfa") + assert check_gbz.task_name == haplotypes.task_name + assert check_gfa.task_name == haplotypes.task_name + + def test_pangenome_check_commands(self): + """The checks run the resolved name and prefix""" + pipeline = self.create_pipeline() + pipeline.pangenome_ref_name = "CHM13" + pipeline.pangenome_contig_prefix = "CHM13#0#" + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + gbz_cmd = str(self._get_job(all_jobs, "check-sample-gbz").shell) + assert "sentieon_cli.pangenome_meta check-gbz" in gbz_cmd + assert f"--gbz {self.mock_dir}/sample_pangenome.gbz" in gbz_cmd + assert "--reference_name CHM13" in gbz_cmd + # The `#` characters are shell-quoted in the rendered command + assert "--contig_prefix 'CHM13#0#'" in gbz_cmd + + gfa_cmd = str(self._get_job(all_jobs, "check-sample-gfa").shell) + assert "sentieon_cli.pangenome_meta check-gfa" in gfa_cmd + assert f"--gfa {self.mock_dir}/sample-hap.raw.gfa" in gfa_cmd + assert "--reference_name CHM13" in gfa_cmd + assert "--contig_prefix 'CHM13#0#'" in gfa_cmd + + def test_the_resolved_reference_name_reaches_vg(self): + """`vg haplotypes` and `vg convert` use the resolved name""" + pipeline = self.create_pipeline() + pipeline.pangenome_ref_name = "CHM13" + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + haplotypes = str(self._get_job(all_jobs, "vg-haplotypes").shell) + assert "--set-reference CHM13" in haplotypes + convert = str(self._get_job(all_jobs, "vg-convert-gfa").shell) + assert "-Q CHM13" in convert + def test_pangenome_contig_prefix(self): """`--pangenome_contig_prefix` reaches every graph consumer""" pipeline = self.create_pipeline() diff --git a/tests/unit/test_pangenome_meta.py b/tests/unit/test_pangenome_meta.py new file mode 100644 index 0000000..aa37e42 --- /dev/null +++ b/tests/unit/test_pangenome_meta.py @@ -0,0 +1,394 @@ +""" +Unit tests for the pangenome metadata reader and its run-time checks. + +The binary fixtures in `tests/data/pangenome` are built by the +`generate.sh` script next to them; see that script for the `vg` commands +and for what each graph contains. +""" + +import pathlib +import struct +import subprocess +import sys + +import pytest + +from sentieon_cli.pangenome_meta import ( + HAPL_MAGIC, + PangenomeMetadataError, + check_gfa, + check_sample_gbz, + detect_pangenome_reference, + main, + read_gbz_metadata, + read_hapl_header, + unpack_bits, +) + +DATA = pathlib.Path(__file__).resolve().parents[1] / "data" / "pangenome" + +# The tiny graphs hold two reference samples over two 30 bp contigs +TINY_CONTIGS = ["chr1", "chr2"] + + +class TestUnpackBits: + """The packed-integer decoder behind every `simple-sds` vector""" + + def test_byte_wide_values(self): + words = struct.pack("` writes""" + with pytest.raises(PangenomeMetadataError) as excinfo: + check_sample_gbz( + DATA / "tiny_no_ref_tag.gbz", "GRCh38", "GRCh38#0#" + ) + message = str(excinfo.value) + assert "reference_samples=''" in message + assert "reference name may be wrong" in message + + def test_a_reference_sample_without_the_prefix_fails(self): + with pytest.raises( + PangenomeMetadataError, match="contig prefix 'GRCh38#1#'" + ): + check_sample_gbz( + DATA / "tiny_sampled.gbz", "GRCh38", "GRCh38#1#" + ) + + +class TestCheckGfa: + """The run-time check on the GFA `vg convert -Q` writes""" + + def test_an_rgfa_with_the_reference_tags_passes(self): + check_gfa(DATA / "rgfa_ok.gfa", "GRCh38", "GRCh38#0#") + + def test_a_header_that_still_lists_the_reference_fails(self): + with pytest.raises(PangenomeMetadataError) as excinfo: + check_gfa( + DATA / "rgfa_wrong_reference.gfa", "GRCh38", "GRCh38#0#" + ) + message = str(excinfo.value) + assert "RS:Z:CHM13 GRCh38" in message + assert "reference name may be wrong" in message + + def test_a_file_without_sn_tags_fails(self): + with pytest.raises(PangenomeMetadataError, match="SN:Z:GRCh38#0#"): + check_gfa(DATA / "rgfa_no_tags.gfa", "GRCh38", "GRCh38#0#") + + def test_a_file_without_a_header_line_is_tolerated(self): + check_gfa(DATA / "rgfa_no_header.gfa", "GRCh38", "GRCh38#0#") + + def test_the_wrong_contig_prefix_fails(self): + """The header passes, but no segment carries the prefix""" + with pytest.raises(PangenomeMetadataError, match="SN:Z:GRCh38#1#"): + check_gfa(DATA / "rgfa_ok.gfa", "GRCh38", "GRCh38#1#") + + def test_a_reference_name_in_the_header_fails_first(self): + """`RS:Z:` naming the reference is checked before the tags""" + with pytest.raises(PangenomeMetadataError, match="RS:Z:CHM13"): + check_gfa(DATA / "rgfa_ok.gfa", "CHM13", "CHM13#0#") + + def test_only_the_head_of_the_file_is_read(self): + """The tags are past `max_lines`, so the check gives up""" + with pytest.raises( + PangenomeMetadataError, match="first 1 lines" + ): + check_gfa( + DATA / "rgfa_ok.gfa", "GRCh38", "GRCh38#0#", max_lines=1 + ) + + +class TestCli: + """The `python -m sentieon_cli.pangenome_meta` entry point""" + + def test_detect_prints_the_reference(self, capsys): + assert main(["detect", "--gbz", str(DATA / "tiny.gbz")]) == 0 + out = capsys.readouterr().out + assert "pangenome_ref_name\tGRCh38" in out + assert "pangenome_contig_prefix\tGRCh38#0#" in out + assert "n_contigs\t2" in out + + def test_detect_fails_on_a_graph_without_a_backbone(self): + argv = ["detect", "--gbz", str(DATA / "tiny_no_backbone.gbz")] + assert main(argv) == 1 + + def test_check_gbz_passes(self): + argv = [ + "check-gbz", + "--gbz", + str(DATA / "tiny_sampled.gbz"), + "--reference_name", + "GRCh38", + "--contig_prefix", + "GRCh38#0#", + ] + assert main(argv) == 0 + + def test_check_gbz_fails(self): + argv = [ + "check-gbz", + "--gbz", + str(DATA / "tiny.gbz"), + "--reference_name", + "GRCh38", + "--contig_prefix", + "GRCh38#0#", + ] + assert main(argv) == 1 + + def test_check_gbz_fails_on_a_missing_file(self): + argv = [ + "check-gbz", + "--gbz", + str(DATA / "does-not-exist.gbz"), + "--reference_name", + "GRCh38", + "--contig_prefix", + "GRCh38#0#", + ] + assert main(argv) == 1 + + def test_check_gfa_passes(self): + argv = [ + "check-gfa", + "--gfa", + str(DATA / "rgfa_ok.gfa"), + "--reference_name", + "GRCh38", + "--contig_prefix", + "GRCh38#0#", + ] + assert main(argv) == 0 + + def test_check_gfa_fails(self): + argv = [ + "check-gfa", + "--gfa", + str(DATA / "rgfa_wrong_reference.gfa"), + "--reference_name", + "GRCh38", + "--contig_prefix", + "GRCh38#0#", + ] + assert main(argv) == 1 + + def test_check_gfa_honors_max_lines(self): + argv = [ + "check-gfa", + "--gfa", + str(DATA / "rgfa_ok.gfa"), + "--reference_name", + "GRCh38", + "--contig_prefix", + "GRCh38#0#", + "--max_lines", + "1", + ] + assert main(argv) == 1 + + def test_the_module_runs_as_a_script(self): + """The DAG jobs run the module with `python -m`""" + failing = subprocess.run( + [ + sys.executable, + "-m", + "sentieon_cli.pangenome_meta", + "check-gbz", + "--gbz", + str(DATA / "tiny.gbz"), + "--reference_name", + "GRCh38", + "--contig_prefix", + "GRCh38#0#", + ], + capture_output=True, + text=True, + ) + assert failing.returncode == 1 + assert "reference name may be wrong" in failing.stderr + + passing = subprocess.run( + [ + sys.executable, + "-m", + "sentieon_cli.pangenome_meta", + "check-gbz", + "--gbz", + str(DATA / "tiny_sampled.gbz"), + "--reference_name", + "GRCh38", + "--contig_prefix", + "GRCh38#0#", + ], + capture_output=True, + text=True, + ) + assert passing.returncode == 0 diff --git a/tests/unit/test_pipeline_validation.py b/tests/unit/test_pipeline_validation.py index f628345..c1bae30 100644 --- a/tests/unit/test_pipeline_validation.py +++ b/tests/unit/test_pipeline_validation.py @@ -8,7 +8,8 @@ import sys import os import json -from unittest.mock import patch +from types import SimpleNamespace +from unittest.mock import MagicMock, patch # Add the parent directory to the path so we can import sentieon_cli sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))) @@ -16,9 +17,21 @@ from sentieon_cli.dnascope import DNAscopePipeline from sentieon_cli.dnascope_longread import DNAscopeLRPipeline from sentieon_cli.dnascope_hybrid import DNAscopeHybridPipeline +from sentieon_cli.hybrid_pangenome import HybridPangenome +from sentieon_cli.pangenome_meta import ( + PangenomeMetadataError, + PangenomeReference, +) +from sentieon_cli.sentieon_pangenome import SentieonPangenome +from sentieon_cli.shard import GRCH38_CONTIGS from sentieon_cli.util import set_bwt_max_mem from tests.utils.test_helpers import create_mock_args +# A reference index that satisfies the GRCh38 contig-length check +GRCH38_CONTIGS_AS_FAI = { + ctg: {"length": length} for ctg, length in GRCH38_CONTIGS.items() +} + class TestDNAscopePipelineValidation: """Test DNAscope pipeline validation logic""" @@ -514,5 +527,321 @@ def test_partial_index_files_missing_raises(self): assert excinfo.value.code == 2 +PANGENOME_PIPELINES = [SentieonPangenome, HybridPangenome] + +GRCH38_REFERENCE = PangenomeReference( + ref_name="GRCh38", + contig_prefix="GRCh38#0#", + contigs=["chr1", "chr2"], + reference_samples=["GRCh38", "CHM13"], + n_paths=2, +) +CHM13_REFERENCE = PangenomeReference( + ref_name="CHM13", + contig_prefix="CHM13#0#", + contigs=["chr1", "chr2"], + reference_samples=["CHM13", "GRCh38"], + n_paths=2, +) + + +@pytest.mark.parametrize( + "cls", PANGENOME_PIPELINES, ids=lambda c: c.__name__ +) +class TestResolvePangenomeReference: + """`BasePangenome.resolve_pangenome_reference` in both pipelines. + + The reference name and contig prefix come from the graph's own + metadata, so a graph built against the wrong reference cannot reach + `vg haplotypes --set-reference`, which would silently drop every + reference path. + """ + + def _pipeline(self, cls, tmp_path, **attributes): + pipeline = cls() + pipeline.logger = MagicMock() + pipeline.gbz = tmp_path / "graph.gbz" + pipeline.hapl = tmp_path / "graph.hapl" + pipeline.gbz.touch() + pipeline.hapl.touch() + pipeline.fai_data = {"chr1": {"length": 10}, "chr2": {"length": 10}} + for key, value in attributes.items(): + setattr(pipeline, key, value) + return pipeline + + def _detect(self, reference): + return patch( + "sentieon_cli.base_pangenome.detect_pangenome_reference", + return_value=reference, + ) + + def _build(self, build): + return patch( + "sentieon_cli.base_pangenome.detect_reference_build", + return_value=build, + ) + + def _hapl(self, chains=2): + return patch( + "sentieon_cli.base_pangenome.read_hapl_header", + return_value=SimpleNamespace(top_level_chains=chains), + ) + + def test_the_detected_reference_populates_both_attributes( + self, cls, tmp_path + ): + pipeline = self._pipeline(cls, tmp_path) + with self._detect(GRCH38_REFERENCE), self._build("hg38"), ( + self._hapl() + ): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "GRCh38" + assert pipeline.pangenome_contig_prefix == "GRCh38#0#" + assert pipeline.pangenome_reference is GRCH38_REFERENCE + + def test_a_chm13_graph_resolves_to_the_chm13_prefix( + self, cls, tmp_path + ): + pipeline = self._pipeline(cls, tmp_path) + with self._detect(CHM13_REFERENCE), self._build("chm13"), ( + self._hapl() + ): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "CHM13" + assert pipeline.pangenome_contig_prefix == "CHM13#0#" + + def test_a_matching_override_is_accepted(self, cls, tmp_path): + pipeline = self._pipeline( + cls, + tmp_path, + pangenome_ref_name="GRCh38", + pangenome_contig_prefix="GRCh38#0#", + ) + with self._detect(GRCH38_REFERENCE), self._build("hg38"), ( + self._hapl() + ): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "GRCh38" + + def test_a_mismatched_name_override_exits(self, cls, tmp_path): + pipeline = self._pipeline(cls, tmp_path, pangenome_ref_name="hg38") + with self._detect(GRCH38_REFERENCE), self._build("hg38"), ( + self._hapl() + ): + with pytest.raises(SystemExit) as excinfo: + pipeline.resolve_pangenome_reference() + assert excinfo.value.code == 2 + pipeline.logger.error.assert_called() + + def test_a_mismatched_prefix_override_exits(self, cls, tmp_path): + pipeline = self._pipeline( + cls, tmp_path, pangenome_contig_prefix="GRCh38#1#" + ) + with self._detect(GRCH38_REFERENCE), self._build("hg38"), ( + self._hapl() + ): + with pytest.raises(SystemExit) as excinfo: + pipeline.resolve_pangenome_reference() + assert excinfo.value.code == 2 + + def test_the_skip_flag_honors_a_mismatched_override( + self, cls, tmp_path + ): + pipeline = self._pipeline( + cls, + tmp_path, + pangenome_ref_name="hg38", + skip_pangenome_name_checks=True, + ) + with self._detect(GRCH38_REFERENCE), self._build("hg38"), ( + self._hapl() + ): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "hg38" + assert pipeline.pangenome_contig_prefix == "GRCh38#0#" + pipeline.logger.warning.assert_called() + pipeline.logger.error.assert_not_called() + + def test_a_contig_missing_from_the_fai_exits(self, cls, tmp_path): + pipeline = self._pipeline(cls, tmp_path) + pipeline.fai_data = {"chr1": {"length": 10}} + with self._detect(GRCH38_REFERENCE), self._build("hg38"), ( + self._hapl() + ): + with pytest.raises(SystemExit) as excinfo: + pipeline.resolve_pangenome_reference() + assert excinfo.value.code == 2 + message = pipeline.logger.error.call_args[0] + assert "missing from the reference FASTA index" in message[0] + assert "chr2" in message + + def test_the_skip_flag_downgrades_a_missing_contig(self, cls, tmp_path): + pipeline = self._pipeline( + cls, tmp_path, skip_pangenome_name_checks=True + ) + pipeline.fai_data = {"chr1": {"length": 10}} + with self._detect(GRCH38_REFERENCE), self._build("hg38"), ( + self._hapl() + ): + pipeline.resolve_pangenome_reference() + pipeline.logger.warning.assert_called() + + def test_a_reference_build_mismatch_exits(self, cls, tmp_path): + """A CHM13 graph with an hg38 reference FASTA. + + Contig names alone cannot catch this: both builds name their + chromosomes chr1..chrM. + """ + pipeline = self._pipeline(cls, tmp_path) + with self._detect(CHM13_REFERENCE), self._build("hg38"), ( + self._hapl() + ): + with pytest.raises(SystemExit) as excinfo: + pipeline.resolve_pangenome_reference() + assert excinfo.value.code == 2 + assert "reference build" in pipeline.logger.error.call_args[0][0] + + def test_an_unknown_reference_build_is_only_logged(self, cls, tmp_path): + pipeline = self._pipeline(cls, tmp_path) + with self._detect(GRCH38_REFERENCE), self._build(None), ( + self._hapl() + ): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "GRCh38" + pipeline.logger.error.assert_not_called() + + def test_a_hapl_chain_count_mismatch_only_warns(self, cls, tmp_path): + pipeline = self._pipeline(cls, tmp_path) + with self._detect(GRCH38_REFERENCE), self._build("hg38"), ( + self._hapl(chains=195) + ): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "GRCh38" + pipeline.logger.warning.assert_called() + + def test_an_unreadable_hapl_exits(self, cls, tmp_path): + pipeline = self._pipeline(cls, tmp_path) + hapl_error = PangenomeMetadataError("not a vg haplotypes file") + with self._detect(GRCH38_REFERENCE), self._build("hg38"), patch( + "sentieon_cli.base_pangenome.read_hapl_header", + side_effect=hapl_error, + ): + with pytest.raises(SystemExit) as excinfo: + pipeline.resolve_pangenome_reference() + assert excinfo.value.code == 2 + + def test_a_dry_run_falls_back_to_the_defaults(self, cls, tmp_path): + """Mock graph files cannot be parsed, but dry runs still work""" + pipeline = self._pipeline(cls, tmp_path, dry_run=True) + error = PangenomeMetadataError("not a GBZ pangenome file") + with patch( + "sentieon_cli.base_pangenome.detect_pangenome_reference", + side_effect=error, + ): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "GRCh38" + assert pipeline.pangenome_contig_prefix == "GRCh38#0#" + assert pipeline.pangenome_reference is None + + def test_a_dry_run_keeps_the_overrides(self, cls, tmp_path): + pipeline = self._pipeline( + cls, + tmp_path, + dry_run=True, + pangenome_ref_name="CHM13", + pangenome_contig_prefix="CHM13#0#", + ) + error = PangenomeMetadataError("not a GBZ pangenome file") + with patch( + "sentieon_cli.base_pangenome.detect_pangenome_reference", + side_effect=error, + ): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "CHM13" + assert pipeline.pangenome_contig_prefix == "CHM13#0#" + + def test_an_unreadable_graph_exits_outside_a_dry_run( + self, cls, tmp_path + ): + pipeline = self._pipeline(cls, tmp_path) + error = PangenomeMetadataError("not a GBZ pangenome file") + with patch( + "sentieon_cli.base_pangenome.detect_pangenome_reference", + side_effect=error, + ): + with pytest.raises(SystemExit) as excinfo: + pipeline.resolve_pangenome_reference() + assert excinfo.value.code == 2 + + def test_both_overrides_survive_an_unreadable_graph(self, cls, tmp_path): + pipeline = self._pipeline( + cls, + tmp_path, + pangenome_ref_name="CHM13", + pangenome_contig_prefix="CHM13#0#", + ) + error = PangenomeMetadataError("not a GBZ pangenome file") + with patch( + "sentieon_cli.base_pangenome.detect_pangenome_reference", + side_effect=error, + ), self._hapl(): + pipeline.resolve_pangenome_reference() + assert pipeline.pangenome_ref_name == "CHM13" + assert pipeline.pangenome_contig_prefix == "CHM13#0#" + pipeline.logger.warning.assert_called() + + +@pytest.mark.parametrize( + "cls", PANGENOME_PIPELINES, ids=lambda c: c.__name__ +) +class TestGrch38ContigChecks: + """The GRCh38 contig-length check only runs for GRCh38 pangenomes""" + + def _pipeline(self, cls, **attributes): + pipeline = cls() + pipeline.logger = MagicMock() + pipeline.pangenome_ref_name = "GRCh38" + pipeline.fai_data = dict(GRCH38_CONTIGS_AS_FAI) + pipeline.pop_vcf_contigs = dict(GRCH38_CONTIGS) + for key, value in attributes.items(): + setattr(pipeline, key, value) + return pipeline + + def test_matching_contigs_pass(self, cls): + self._pipeline(cls).validate_grch38_contigs() + + def test_a_wrong_contig_length_exits(self, cls): + pipeline = self._pipeline(cls) + pipeline.fai_data["chr1"] = {"length": 1} + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_grch38_contigs() + assert excinfo.value.code == 2 + + def test_skip_contig_checks_bypasses_the_check(self, cls): + pipeline = self._pipeline(cls, skip_contig_checks=True) + pipeline.fai_data["chr1"] = {"length": 1} + pipeline.validate_grch38_contigs() + + def test_a_non_grch38_pangenome_skips_the_check(self, cls): + """A CHM13 reference has different contig lengths""" + pipeline = self._pipeline(cls, pangenome_ref_name="CHM13") + pipeline.fai_data = {"chr1": {"length": 1}} + pipeline.pop_vcf_contigs = {} + pipeline.validate_grch38_contigs() + assert "GRCh38" in pipeline.logger.info.call_args[0][0] + + def test_a_wrong_pop_vcf_contig_exits(self, cls): + pipeline = self._pipeline(cls) + pipeline.pop_vcf_contigs = {} + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_grch38_contigs() + assert excinfo.value.code == 2 + + def test_a_dry_run_skips_the_pop_vcf_check(self, cls): + pipeline = self._pipeline(cls, dry_run=True) + pipeline.pop_vcf_contigs = {} + pipeline.validate_grch38_contigs() + + if __name__ == "__main__": pytest.main([__file__]) diff --git a/tests/unit/test_sentieon_pangenome.py b/tests/unit/test_sentieon_pangenome.py index ca6e9fa..9d96fe2 100644 --- a/tests/unit/test_sentieon_pangenome.py +++ b/tests/unit/test_sentieon_pangenome.py @@ -75,6 +75,10 @@ def create_pipeline(self): pipeline.skip_pop_vcf_id_check = True pipeline.tmp_dir = self.mock_dir + # The pangenome reference, normally resolved by validate() + pipeline.pangenome_ref_name = "GRCh38" + pipeline.pangenome_contig_prefix = "GRCh38#0#" + # Mock parsing fai pipeline.fai_data = {"chr1": {"length": 1000}} pipeline.shards = [MagicMock()] @@ -607,6 +611,81 @@ def test_skip_small_variants_with_call_svs(self): assert "--algo PangenomeSV" in cmd_str assert "--algo DNAscope" not in cmd_str + def test_pangenome_check_jobs(self): + """The graph checks sit between vg and everything downstream. + + A wrong `--pangenome_ref_name` is silent in `vg`: the sampled + graph loses its reference paths and the GFA loses its rGFA tags, + and `pgutil lift` then leaves every read unmapped. + """ + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + job_names, all_jobs = self._get_all_job_names(dag) + + assert "check-sample-gbz" in job_names + assert "check-sample-gfa" in job_names + + check_gbz = next(j for j in all_jobs if j.name == "check-sample-gbz") + check_gfa = next(j for j in all_jobs if j.name == "check-sample-gfa") + haplotypes = next(j for j in all_jobs if j.name == "vg-haplotypes") + gfa = next(j for j in all_jobs if j.name == "vg-convert-gfa") + fasta = next(j for j in all_jobs if j.name == "vg-paths-fasta") + mm2 = next(j for j in all_jobs if j.name == "mm2-lift") + + assert self._get_deps(dag, check_gbz) == {haplotypes} + assert self._get_deps(dag, gfa) == {check_gbz} + assert self._get_deps(dag, fasta) == {check_gbz} + assert self._get_deps(dag, check_gfa) == {gfa} + assert check_gfa in self._get_deps(dag, mm2) + assert "vg-convert-gfa" not in self._get_dep_names(dag, mm2) + + # The checks share the log group of the jobs they guard + assert check_gbz.task_name == haplotypes.task_name + assert check_gfa.task_name == gfa.task_name + + def test_pangenome_check_commands(self): + """The checks run the resolved name and prefix""" + pipeline = self.create_pipeline() + pipeline.pangenome_ref_name = "CHM13" + pipeline.pangenome_contig_prefix = "CHM13#0#" + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + gbz_cmd = str( + next(j for j in all_jobs if j.name == "check-sample-gbz").shell + ) + assert "sentieon_cli.pangenome_meta check-gbz" in gbz_cmd + # runpy's duplicate-module warning is kept out of the job logs + assert "-W ignore::RuntimeWarning:runpy" in gbz_cmd + assert f"--gbz {self.mock_dir}/sample_pangenome.gbz" in gbz_cmd + assert "--reference_name CHM13" in gbz_cmd + # The `#` characters are shell-quoted in the rendered command + assert "--contig_prefix 'CHM13#0#'" in gbz_cmd + + gfa_cmd = str( + next(j for j in all_jobs if j.name == "check-sample-gfa").shell + ) + assert "sentieon_cli.pangenome_meta check-gfa" in gfa_cmd + assert f"--gfa {self.mock_dir}/sample-hap.gfa" in gfa_cmd + assert "--reference_name CHM13" in gfa_cmd + assert "--contig_prefix 'CHM13#0#'" in gfa_cmd + + def test_the_resolved_reference_name_reaches_vg(self): + """`vg haplotypes` and `vg convert` use the resolved name""" + pipeline = self.create_pipeline() + pipeline.pangenome_ref_name = "CHM13" + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + haplotypes = str( + next(j for j in all_jobs if j.name == "vg-haplotypes").shell + ) + assert "--set-reference CHM13" in haplotypes + convert = str( + next(j for j in all_jobs if j.name == "vg-convert-gfa").shell + ) + assert "-Q CHM13" in convert + def _get_all_job_names(self, dag): """Helper to get all job names from a DAG""" all_jobs = list(dag.waiting_jobs.keys()) + list(dag.ready_jobs.keys()) diff --git a/tests/unit/test_sentieon_pangenome_kmc_integration.py b/tests/unit/test_sentieon_pangenome_kmc_integration.py index eaad07e..9b6d8b8 100644 --- a/tests/unit/test_sentieon_pangenome_kmc_integration.py +++ b/tests/unit/test_sentieon_pangenome_kmc_integration.py @@ -21,6 +21,10 @@ def test_validate_kmc_patch_success(self): pipeline.reference = "ref" pipeline.pop_vcf = "pop" pipeline.gbz = "ref.grch38.gbz" + pipeline.hapl = "ref.grch38.hapl" + # The mock graph files cannot be parsed; a dry run falls back to + # the default pangenome reference + pipeline.dry_run = True # Mock dependencies pipeline.logger = MagicMock() @@ -56,6 +60,9 @@ def test_validate_kmc_patch_failure(self): pipeline.validate_ref = MagicMock() pipeline.collect_readgroups = MagicMock() pipeline.bed = "bed" + pipeline.gbz = "ref.grch38.gbz" + pipeline.hapl = "ref.grch38.hapl" + pipeline.dry_run = True pipeline.logger = MagicMock() From fb193d34de276d1aee34e0451500e430a3357c05 Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sat, 12 Sep 2026 14:10:21 -0700 Subject: [PATCH 2/2] Check pop VCF contig lengths against the reference for any backbone `validate_grch38_contigs` compared the reference index and the `--pop_vcf` with a hard-coded table of the 25 GRCh38 chromosome lengths, so it only ran for a GRCh38 backbone and never looked at the 170 unplaced contigs the HPRC graph also carries. Drop the table and compare the two files instead: every contig of the backbone reference detected in the GBZ must carry the same length in the `--pop_vcf` header as in the reference FASTA index. A CHM13 graph is now checked the same way a GRCh38 one is, a contig absent from the pop VCF or without a length is an error, and the message names the first ten offenders with both lengths. Only the hidden `--skip_contig_checks` flag bypasses the check; `--skip_pangenome_name_checks` does not. The pop VCF header is now parsed in pure Python. `vcf_contigs` and `vcf_id` shelled out to `bcftools view -h`, and `vcf_contigs` returned a fake three-contig dictionary in dry runs, which left the dry run checking nothing and the real run checking everything. `read_vcf_header` opens the file with `gzip` when its first two bytes say so -- the extension is not trusted -- and stops at the first line that is not a header line, so the body is never read: the header of the 10 GB population VCF parses in about 1 ms. Dry runs now parse it too, and fall back only when the file cannot be read at all, which is what the unit tests and the placeholder files behind them need. The `SentieonVcfID` comparison follows the same rule, so a dry run catches a pop VCF that does not belong with the model bundle. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01RWPD3suVQRpMZbJHijXhao --- sentieon_cli/base_pangenome.py | 173 ++++++++++--- sentieon_cli/dnascope_hybrid.py | 2 +- sentieon_cli/dnascope_longread.py | 2 +- sentieon_cli/hybrid_pangenome.py | 23 +- sentieon_cli/sentieon_pangenome.py | 20 +- sentieon_cli/shard.py | 74 ++---- sentieon_cli/util.py | 93 +++++-- tests/unit/test_pipeline_validation.py | 233 +++++++++++++++--- ...test_sentieon_pangenome_kmc_integration.py | 7 +- tests/unit/test_vcf_header.py | 175 +++++++++++++ 10 files changed, 628 insertions(+), 174 deletions(-) create mode 100644 tests/unit/test_vcf_header.py diff --git a/sentieon_cli/base_pangenome.py b/sentieon_cli/base_pangenome.py index 9b6c82d..5b2f328 100644 --- a/sentieon_cli/base_pangenome.py +++ b/sentieon_cli/base_pangenome.py @@ -6,7 +6,7 @@ import copy import pathlib import sys -from typing import Dict, Iterable, List, Optional, Sequence, Set +from typing import Dict, Iterable, List, Optional, Sequence from . import command_strings as cmds from .dag import DAG @@ -29,7 +29,7 @@ read_hapl_header, ) from .pipeline import BasePipeline -from .shard import GRCH38_CONTIGS, detect_reference_build +from .shard import detect_reference_build, parse_vcf_contigs from .stages.base import StageContext from .stages.cnv import PangenomeCNVResult, PangenomeCNVStage from .stages.expansion import ExpansionHunterResult, ExpansionHunterStage @@ -42,7 +42,14 @@ T1KResult, T1KStage, ) -from .util import path_arg, require_versions, sample_sex_arg +from .util import ( + VcfHeaderError, + path_arg, + read_vcf_header, + require_versions, + sample_sex_arg, + vcf_id_from_header, +) # The pangenome reference of the graphs the pipelines shipped with, used # when a graph's own metadata cannot be read @@ -245,8 +252,11 @@ def __init__(self) -> None: self.pangenome_reference: Optional[PangenomeReference] = None self.skip_contig_checks = False self.skip_pangenome_name_checks = False - # Parsed by `validate`; declared here for the shared checks + # Parsed by `validate`; declared here for the shared checks. + # `pop_vcf_header` stays `None` when the `--pop_vcf` header + # could not be read, which only a dry run survives. self.fai_data: Dict[str, Dict[str, int]] = {} + self.pop_vcf_header: Optional[List[str]] = None self.pop_vcf_contigs: Dict[str, Optional[int]] = {} # Stashed by `build_dag` for the second, sex-aware DAG self.ploidy_json: Optional[pathlib.Path] = None @@ -490,50 +500,145 @@ def check_pangenome_hapl( detected.n_paths, ) - def validate_grch38_contigs(self) -> None: - """Check the reference and pop VCF contig lengths against GRCh38. + def load_pop_vcf_header(self, pop_vcf: Optional[pathlib.Path]) -> None: + """Parse the `--pop_vcf` header. - The lengths only describe GRCh38; other pangenome references are - covered by the reference FASTA checks in - `resolve_pangenome_reference` instead. + Sets `pop_vcf_header` and `pop_vcf_contigs`. The header is read + in dry runs too, so a dry run checks the pop VCF the way the real + run will. Dry runs (and the unit tests behind them) are given + placeholder files, so an unreadable header is only logged there + and leaves both attributes empty; a real run cannot continue + without the annotations the pop VCF carries. """ - if self.skip_contig_checks: + self.pop_vcf_header = None + self.pop_vcf_contigs = {} + if pop_vcf is None: return - if self.pangenome_ref_name != DEFAULT_PANGENOME_REF_NAME: - self.logger.info( - "The pangenome reference is '%s', so the GRCh38 " - "contig-length check is skipped", - self.pangenome_ref_name, + try: + header = read_vcf_header(pop_vcf) + except (VcfHeaderError, OSError) as err: + if self.dry_run: + self.logger.debug( + "Could not read the header of the `--pop_vcf` '%s': %s", + pop_vcf, + err, + ) + return + self.logger.error( + "Could not read the header of the `--pop_vcf` '%s': %s", + pop_vcf, + err, ) - return + sys.exit(2) + self.pop_vcf_header = header + self.pop_vcf_contigs = parse_vcf_contigs(header) + self.logger.debug("VCF contigs are: %s", self.pop_vcf_contigs) - # Check the fai file contigs - mismatch_contigs: Set[str] = set() - for ctg, length in GRCH38_CONTIGS.items(): - fai_length = self.fai_data.get(ctg, {}).get("length", -1) - if length != fai_length: - mismatch_contigs.add(ctg) - if mismatch_contigs: + def check_pop_vcf_id(self, bundle_vcf_id: str) -> None: + """Compare the `--pop_vcf` ID with the one the bundle expects. + + Called from `validate_bundle` in both pipelines, dry runs + included, so a dry run catches a pop VCF that does not belong + with the model bundle. Only a dry run reaches this without a + header: `load_pop_vcf_header` has already failed the real run. + """ + if self.pop_vcf_header is None: + self.logger.debug( + "The `--pop_vcf` header could not be read, so its " + "SentieonVcfID was not compared with the model bundle's" + ) + return + pop_vcf_id = vcf_id_from_header(self.pop_vcf_header) + if bundle_vcf_id != pop_vcf_id: self.logger.error( - "Reference contigs with unexpected lengths: %s", - ", ".join(mismatch_contigs), + "The ID of the `--pop_vcf` does not match the model " + "bundle. The bundle expects '%s' and the `--pop_vcf` " + "carries '%s'.", + bundle_vcf_id, + pop_vcf_id, ) sys.exit(2) - # Check the pop VCF file contigs - if self.dry_run: + def validate_pangenome_contig_lengths(self) -> None: + """Check the `--pop_vcf` contig lengths against the reference. + + `DNAModelApply` transfers the pop VCF's annotations by position, + so a pop VCF built against another reference silently annotates + the wrong sites. Every contig of the pangenome backbone + reference must therefore carry the same length in the + `--pop_vcf` header as in the reference FASTA index, whatever the + build: the lengths come from the two files rather than from a + table of GRCh38 lengths. + + The backbone contigs come from the graph's own metadata. When + that could not be read -- a dry run with a placeholder graph, or + a real run driven by the hidden `--pangenome_ref_name` override + -- the check falls back to the contigs the index and the pop VCF + share. That fallback cannot report a contig missing from the pop + VCF, but it still catches a pop VCF from another build. + + `--skip_pangenome_name_checks` does not downgrade this check; + only the hidden `--skip_contig_checks` flag bypasses it. + """ + if self.skip_contig_checks: return - mismatch_contigs = set() - for ctg, length in GRCH38_CONTIGS.items(): - if length != self.pop_vcf_contigs.get(ctg, -1): - mismatch_contigs.add(ctg) - if mismatch_contigs: + + if not self.pop_vcf_contigs: + message = ( + "No `--pop_vcf` contigs are known, so their lengths were " + "not checked against the reference FASTA" + ) + if self.dry_run: + self.logger.debug(message) + else: + self.logger.warning(message) + return + + if self.pangenome_reference is not None: + contigs: List[str] = list(self.pangenome_reference.contigs) + else: + contigs = sorted(set(self.fai_data) & set(self.pop_vcf_contigs)) + + checked = 0 + mismatched: List[str] = [] + for ctg in contigs: + fai_length = self.fai_data.get(ctg, {}).get("length") + if fai_length is None: + # `check_pangenome_reference_fasta` has already reported + # a backbone contig that the reference index is missing + continue + checked += 1 + if ctg not in self.pop_vcf_contigs: + found = "absent" + else: + pop_length = self.pop_vcf_contigs[ctg] + if pop_length == fai_length: + continue + found = "no length" if pop_length is None else str(pop_length) + mismatched.append( + f"{ctg} (reference {fai_length}, pop VCF {found})" + ) + + if mismatched: self.logger.error( - "Pop VCF contigs with unexpected lengths: %s", - ", ".join(mismatch_contigs), + "%d of the %d contig(s) of the pangenome reference '%s' do " + "not have the reference FASTA's length in the `--pop_vcf`: " + "%s. The `--pop_vcf` must be built against the same " + "reference as the `--gbz` pangenome.", + len(mismatched), + checked, + self.pangenome_ref_name, + ", ".join(mismatched[:10]), ) sys.exit(2) + self.logger.info( + "The `--pop_vcf` and the reference FASTA agree on the length of " + "all %d contig(s) of the pangenome reference '%s'", + checked, + self.pangenome_ref_name, + ) + def ref_name(self) -> str: """The resolved pangenome reference name. diff --git a/sentieon_cli/dnascope_hybrid.py b/sentieon_cli/dnascope_hybrid.py index d8e1f38..8ad702b 100644 --- a/sentieon_cli/dnascope_hybrid.py +++ b/sentieon_cli/dnascope_hybrid.py @@ -377,7 +377,7 @@ def validate(self) -> None: self.fai_data = parse_fai(pathlib.Path(str(self.reference) + ".fai")) self.pop_vcf_contigs: Dict[str, Optional[int]] = {} if self.pop_vcf: - self.pop_vcf_contigs = vcf_contigs(self.pop_vcf, self.dry_run) + self.pop_vcf_contigs = vcf_contigs(self.pop_vcf) self.logger.debug("VCF contigs are: %s", self.pop_vcf_contigs) self.shards = determine_shards_from_fai( self.fai_data, 10 * 1000 * 1000 diff --git a/sentieon_cli/dnascope_longread.py b/sentieon_cli/dnascope_longread.py index 0bd9eaf..032f653 100644 --- a/sentieon_cli/dnascope_longread.py +++ b/sentieon_cli/dnascope_longread.py @@ -303,7 +303,7 @@ def validate(self) -> None: self.fai_data = parse_fai(pathlib.Path(str(self.reference) + ".fai")) self.pop_vcf_contigs: Dict[str, Optional[int]] = {} if self.pop_vcf: - self.pop_vcf_contigs = vcf_contigs(self.pop_vcf, self.dry_run) + self.pop_vcf_contigs = vcf_contigs(self.pop_vcf) self.logger.debug("VCF contigs are: %s", self.pop_vcf_contigs) self.shards = determine_shards_from_fai( self.fai_data, 10 * 1000 * 1000 diff --git a/sentieon_cli/hybrid_pangenome.py b/sentieon_cli/hybrid_pangenome.py index ea83b2b..d3ae092 100644 --- a/sentieon_cli/hybrid_pangenome.py +++ b/sentieon_cli/hybrid_pangenome.py @@ -30,7 +30,6 @@ detect_reference_build, determine_shards_from_fai, parse_fai, - vcf_contigs, ) from .shell_pipeline import Command, Pipeline from .stages.alignment import ( @@ -61,7 +60,6 @@ path_arg, require_versions, total_memory, - vcf_id, ) HYBRID_PANGENOME_MIN_VERSIONS = { @@ -244,10 +242,7 @@ def validate(self) -> None: # Before `validate_bundle`, which picks the bundle's # `extract..model` member by the resolved name self.resolve_pangenome_reference() - self.pop_vcf_contigs = {} - if self.pop_vcf: - self.pop_vcf_contigs = vcf_contigs(self.pop_vcf, self.dry_run) - self.logger.debug("VCF contigs are: %s", self.pop_vcf_contigs) + self.load_pop_vcf_header(self.pop_vcf) self.validate_bundle() self.validate_output_vcf() @@ -302,7 +297,7 @@ def validate(self) -> None: "across decoy and unplaced contigs." ) - self.validate_grch38_contigs() + self.validate_pangenome_contig_lengths() # After the version gate and the contig checks, so the tool choice # is only logged for runs that proceed. The `reference_build` it @@ -486,18 +481,8 @@ def validate_bundle(self) -> None: self.has_cnv_model = "cnv.model" in bundle_members bundle_vcf_id = bundle_info.get("SentieonVcfID") - if ( - bundle_vcf_id - and not self.skip_pop_vcf_id_check - and not self.dry_run - ): - assert self.pop_vcf is not None - pop_vcf_id = vcf_id(self.pop_vcf) - if bundle_vcf_id != pop_vcf_id: - self.logger.error( - "The ID of the `--pop_vcf` does not match the model bundle" - ) - sys.exit(2) + if bundle_vcf_id and not self.skip_pop_vcf_id_check: + self.check_pop_vcf_id(bundle_vcf_id) def collect_readgroups(self) -> None: """Collect readgroup tags from the inputs""" diff --git a/sentieon_cli/sentieon_pangenome.py b/sentieon_cli/sentieon_pangenome.py index 27e258d..ed445e8 100644 --- a/sentieon_cli/sentieon_pangenome.py +++ b/sentieon_cli/sentieon_pangenome.py @@ -51,12 +51,10 @@ path_arg, require_versions, total_memory, - vcf_id, ) from .shard import ( determine_shards_from_fai, parse_fai, - vcf_contigs, ) SENT_PANGENOME_MIN_VERSIONS = { @@ -175,10 +173,7 @@ def validate(self) -> None: # Before `validate_bundle`, which picks the bundle's # `extract..model` member by the resolved name self.resolve_pangenome_reference() - self.pop_vcf_contigs = {} - if self.pop_vcf: - self.pop_vcf_contigs = vcf_contigs(self.pop_vcf, self.dry_run) - self.logger.debug("VCF contigs are: %s", self.pop_vcf_contigs) + self.load_pop_vcf_header(self.pop_vcf) self.validate_bundle() self.validate_fastq_rg() @@ -227,7 +222,7 @@ def validate(self) -> None: "across decoy and unplaced contigs." ) - self.validate_grch38_contigs() + self.validate_pangenome_contig_lengths() def validate_segdup(self) -> None: if self.segdup_caller is None: @@ -284,7 +279,7 @@ def validate_cnv(self) -> None: require_versions(CNV_MIN_VERSIONS, skip=self.skip_version_check) def validate_bundle(self) -> None: - pop_vcf = self.required(self.pop_vcf, "pop_vcf") + self.required(self.pop_vcf, "pop_vcf") gbz = self.required(self.gbz, "gbz") bundle_info_bytes = ar_load( str(self.model_bundle) + "/bundle_info.json" @@ -357,13 +352,8 @@ def validate_bundle(self) -> None: "PangenomeSV will still run." ) - if not self.skip_pop_vcf_id_check and not self.dry_run: - pop_vcf_id = vcf_id(pop_vcf) - if bundle_vcf_id != pop_vcf_id: - self.logger.error( - "The ID of the `--pop_vcf` does not match the model bundle" - ) - sys.exit(2) + if not self.skip_pop_vcf_id_check: + self.check_pop_vcf_id(bundle_vcf_id) def validate_fastq_rg(self) -> None: if len(self.r1_fastq) != len(self.r2_fastq): diff --git a/sentieon_cli/shard.py b/sentieon_cli/shard.py index 9e509e5..8d04858 100644 --- a/sentieon_cli/shard.py +++ b/sentieon_cli/shard.py @@ -4,45 +4,16 @@ import pathlib import re -import subprocess as sp -from typing import Dict, List, NamedTuple, Optional +from typing import Dict, Iterable, List, NamedTuple, Optional from importlib.resources import files from .logging import get_logger +from .util import VcfHeaderError, read_vcf_header logger = get_logger(__name__) -GRCH38_CONTIGS: Dict[str, int] = { - "chr1": 248956422, - "chr2": 242193529, - "chr3": 198295559, - "chr4": 190214555, - "chr5": 181538259, - "chr6": 170805979, - "chr7": 159345973, - "chr8": 145138636, - "chr9": 138394717, - "chr10": 133797422, - "chr11": 135086622, - "chr12": 133275309, - "chr13": 114364328, - "chr14": 107043718, - "chr15": 101991189, - "chr16": 90338345, - "chr17": 83257441, - "chr18": 80373285, - "chr19": 58617616, - "chr20": 64444167, - "chr21": 46709983, - "chr22": 50818468, - "chrX": 156040895, - "chrY": 57227415, - "chrM": 16569, -} - - # Signature contigs used to recognize the common human reference builds. # A build is only recognized when every signature contig is present with # the exact expected length. @@ -192,29 +163,16 @@ def determine_shards_from_fai( return shards -def vcf_contigs( - in_vcf: pathlib.Path, dry_run=False +def parse_vcf_contigs( + header: Iterable[str], ) -> Dict[str, Optional[int]]: - """Report the contigs in the input VCF""" - if dry_run: - return { - "chr1": 100, - "chr2": 200, - "chr3": 300, - } + """The ID and length of every `##contig` line of a VCF header. + + A `##contig` line without a `length` maps to `None`. + """ kvpat = re.compile(r'(.*?)=(".*?"|.*?)(?:,|$)') - cmd = ["bcftools", "view", "-h", str(in_vcf)] - p = sp.run(cmd, capture_output=True, text=True) - if p.returncode != 0: - logger.error( - "`%s` failed with return code %d: %s", - " ".join(cmd), - p.returncode, - p.stderr.strip(), - ) - return {} contigs: Dict[str, Optional[int]] = {} - for line in p.stdout.split("\n"): + for line in header: if not line.startswith("##contig"): continue s = line.index("<") @@ -224,3 +182,17 @@ def vcf_contigs( length: Optional[str] = d.get("length", None) contigs[ctg] = int(length) if length else None return contigs + + +def vcf_contigs(in_vcf: pathlib.Path) -> Dict[str, Optional[int]]: + """Report the contigs in the input VCF. + + Returns an empty dictionary, and logs why, when the header could not + be read. + """ + try: + header = read_vcf_header(in_vcf) + except (VcfHeaderError, OSError) as err: + logger.error("%s", err) + return {} + return parse_vcf_contigs(header) diff --git a/sentieon_cli/util.py b/sentieon_cli/util.py index 9dd6794..7d79905 100644 --- a/sentieon_cli/util.py +++ b/sentieon_cli/util.py @@ -4,6 +4,7 @@ import argparse from enum import Enum +import gzip from importlib.metadata import PackageNotFoundError, version import multiprocessing as mp import os @@ -13,7 +14,16 @@ import subprocess as sp import sys import tempfile -from typing import Callable, Dict, List, Mapping, Optional, Tuple +from typing import ( + IO, + Callable, + Dict, + Iterable, + List, + Mapping, + Optional, + Tuple, +) import packaging.version @@ -33,6 +43,9 @@ _UNSAFE = re.compile(r"[^A-Za-z0-9._-]") +# The first two bytes of a gzip (and so of a BGZF) file +GZIP_MAGIC = b"\x1f\x8b" + class SampleSex(Enum): FEMALE = 1 @@ -360,25 +373,77 @@ def parse_rg_line(rg_line: str) -> Dict[str, str]: return parsed -def vcf_id(in_vcf: pathlib.Path) -> Optional[str]: - """Collect the SentieonVcfID header""" - cmd = ["bcftools", "view", "-h", str(in_vcf)] - p = sp.run(cmd, capture_output=True, text=True) - if p.returncode != 0: - logger.error( - "`%s` failed with return code %d: %s", - " ".join(cmd), - p.returncode, - p.stderr.strip(), - ) - return None - for line in p.stdout.split("\n"): +class VcfHeaderError(ValueError): + """A VCF header could not be read""" + + +def _open_vcf_text(in_vcf: pathlib.Path) -> IO[str]: + """Open a VCF for reading, transparently decompressing it. + + The first two bytes decide: a BGZF file is a gzip file, and a VCF is + routinely named without regard to how it is compressed. + """ + with open(in_vcf, "rb") as raw: + compressed = raw.read(len(GZIP_MAGIC)) == GZIP_MAGIC + if compressed: + return gzip.open(in_vcf, "rt") + return open(in_vcf, "rt") + + +def read_vcf_header(in_vcf: pathlib.Path) -> List[str]: + """Read the header lines of a VCF file. + + Lines are read until the first one that is not a header line, so the + body of the file is never touched: the header of a whole-genome + population VCF parses in milliseconds, with no `bcftools` on the + PATH and no subprocess. + + Raises: + VcfHeaderError: the file holds no VCF header, or its header + could not be decompressed or decoded. + OSError: the file does not exist or could not be read. + """ + header: List[str] = [] + try: + with _open_vcf_text(in_vcf) as fh: + for line in fh: + if not line.startswith("#"): + break + header.append(line.rstrip("\n")) + if line.startswith("#CHROM"): + break + except (gzip.BadGzipFile, EOFError, UnicodeDecodeError) as err: + raise VcfHeaderError( + f"{in_vcf}: the VCF header could not be read ({err})" + ) from err + if not header: + raise VcfHeaderError(f"{in_vcf}: the file has no VCF header") + return header + + +def vcf_id_from_header(header: Iterable[str]) -> Optional[str]: + """The `##SentieonVcfID` of VCF header lines, `None` when absent""" + for line in header: if line.startswith("##SentieonVcfID="): i = line.index("=") return line[i + 1 :] # noqa: E203 return None +def vcf_id(in_vcf: pathlib.Path) -> Optional[str]: + """Collect the SentieonVcfID header. + + Returns `None` when the header has no ID and when it could not be + read at all; an unreadable header is logged. + """ + try: + header = read_vcf_header(in_vcf) + except (VcfHeaderError, OSError) as err: + logger.error("%s", err) + return None + return vcf_id_from_header(header) + + def check_kmc_patch(kmc_cmd: str = "kmc") -> bool: """Check if the KMC version supports the required patch""" # Create a temporary directory for the test diff --git a/tests/unit/test_pipeline_validation.py b/tests/unit/test_pipeline_validation.py index c1bae30..ba5a546 100644 --- a/tests/unit/test_pipeline_validation.py +++ b/tests/unit/test_pipeline_validation.py @@ -23,15 +23,9 @@ PangenomeReference, ) from sentieon_cli.sentieon_pangenome import SentieonPangenome -from sentieon_cli.shard import GRCH38_CONTIGS -from sentieon_cli.util import set_bwt_max_mem +from sentieon_cli.util import VcfHeaderError, set_bwt_max_mem from tests.utils.test_helpers import create_mock_args -# A reference index that satisfies the GRCh38 contig-length check -GRCH38_CONTIGS_AS_FAI = { - ctg: {"length": length} for ctg, length in GRCH38_CONTIGS.items() -} - class TestDNAscopePipelineValidation: """Test DNAscope pipeline validation logic""" @@ -794,53 +788,222 @@ def test_both_overrides_survive_an_unreadable_graph(self, cls, tmp_path): @pytest.mark.parametrize( "cls", PANGENOME_PIPELINES, ids=lambda c: c.__name__ ) -class TestGrch38ContigChecks: - """The GRCh38 contig-length check only runs for GRCh38 pangenomes""" +class TestPangenomeContigLengthChecks: + """`BasePangenome.validate_pangenome_contig_lengths` in both pipelines. - def _pipeline(self, cls, **attributes): + The pop VCF is compared with the reference FASTA rather than with a + table of GRCh38 lengths, so the check covers every backbone. + """ + + def _pipeline(self, cls, reference=GRCH38_REFERENCE, **attributes): pipeline = cls() pipeline.logger = MagicMock() - pipeline.pangenome_ref_name = "GRCh38" - pipeline.fai_data = dict(GRCH38_CONTIGS_AS_FAI) - pipeline.pop_vcf_contigs = dict(GRCH38_CONTIGS) + pipeline.pangenome_reference = reference + pipeline.pangenome_ref_name = ( + None if reference is None else reference.ref_name + ) + pipeline.fai_data = { + "chr1": {"length": 1000}, + "chr2": {"length": 2000}, + } + pipeline.pop_vcf_contigs = {"chr1": 1000, "chr2": 2000} for key, value in attributes.items(): setattr(pipeline, key, value) return pipeline - def test_matching_contigs_pass(self, cls): - self._pipeline(cls).validate_grch38_contigs() + def test_matching_lengths_pass(self, cls): + pipeline = self._pipeline(cls) + pipeline.validate_pangenome_contig_lengths() + pipeline.logger.error.assert_not_called() + # The INFO line reports how many backbone contigs were compared + assert pipeline.logger.info.call_args[0][1] == 2 - def test_a_wrong_contig_length_exits(self, cls): + def test_a_mismatching_length_exits(self, cls): pipeline = self._pipeline(cls) - pipeline.fai_data["chr1"] = {"length": 1} + pipeline.pop_vcf_contigs["chr2"] = 2001 with pytest.raises(SystemExit) as excinfo: - pipeline.validate_grch38_contigs() + pipeline.validate_pangenome_contig_lengths() assert excinfo.value.code == 2 + reported = pipeline.logger.error.call_args[0][-1] + assert "chr2" in reported + assert "2000" in reported and "2001" in reported + + def test_a_contig_missing_from_the_pop_vcf_exits(self, cls): + pipeline = self._pipeline(cls) + del pipeline.pop_vcf_contigs["chr2"] + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_pangenome_contig_lengths() + assert excinfo.value.code == 2 + assert "absent" in pipeline.logger.error.call_args[0][-1] + + def test_a_contig_without_a_length_exits(self, cls): + pipeline = self._pipeline(cls) + pipeline.pop_vcf_contigs["chr2"] = None + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_pangenome_contig_lengths() + assert excinfo.value.code == 2 + assert "no length" in pipeline.logger.error.call_args[0][-1] def test_skip_contig_checks_bypasses_the_check(self, cls): pipeline = self._pipeline(cls, skip_contig_checks=True) - pipeline.fai_data["chr1"] = {"length": 1} - pipeline.validate_grch38_contigs() - - def test_a_non_grch38_pangenome_skips_the_check(self, cls): - """A CHM13 reference has different contig lengths""" - pipeline = self._pipeline(cls, pangenome_ref_name="CHM13") - pipeline.fai_data = {"chr1": {"length": 1}} - pipeline.pop_vcf_contigs = {} - pipeline.validate_grch38_contigs() - assert "GRCh38" in pipeline.logger.info.call_args[0][0] - - def test_a_wrong_pop_vcf_contig_exits(self, cls): + pipeline.pop_vcf_contigs["chr1"] = 1 + pipeline.validate_pangenome_contig_lengths() + pipeline.logger.error.assert_not_called() + + def test_the_name_check_flag_does_not_bypass_the_check(self, cls): + """Only `--skip_contig_checks` turns this check off""" + pipeline = self._pipeline(cls, skip_pangenome_name_checks=True) + pipeline.pop_vcf_contigs["chr1"] = 1 + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_pangenome_contig_lengths() + assert excinfo.value.code == 2 + + def test_a_chm13_backbone_is_checked(self, cls): + """The old GRCh38-only check skipped every other backbone""" + pipeline = self._pipeline(cls, reference=CHM13_REFERENCE) + pipeline.pop_vcf_contigs["chr1"] = 248956422 + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_pangenome_contig_lengths() + assert excinfo.value.code == 2 + assert "CHM13" in pipeline.logger.error.call_args[0] + + def test_a_matching_chm13_backbone_passes(self, cls): + pipeline = self._pipeline(cls, reference=CHM13_REFERENCE) + pipeline.validate_pangenome_contig_lengths() + pipeline.logger.error.assert_not_called() + + def test_empty_pop_vcf_contigs_warn_and_skip(self, cls): + pipeline = self._pipeline(cls, pop_vcf_contigs={}) + pipeline.validate_pangenome_contig_lengths() + pipeline.logger.warning.assert_called() + pipeline.logger.error.assert_not_called() + + def test_empty_pop_vcf_contigs_only_debug_in_a_dry_run(self, cls): + pipeline = self._pipeline(cls, pop_vcf_contigs={}, dry_run=True) + pipeline.validate_pangenome_contig_lengths() + pipeline.logger.warning.assert_not_called() + pipeline.logger.debug.assert_called() + + def test_an_undetected_backbone_checks_the_shared_contigs(self, cls): + """Without graph metadata, the fai/pop VCF intersection is used""" + pipeline = self._pipeline(cls, reference=None) + pipeline.pangenome_ref_name = "GRCh38" + pipeline.pop_vcf_contigs = {"chr1": 1001, "chr3": 30} + with pytest.raises(SystemExit) as excinfo: + pipeline.validate_pangenome_contig_lengths() + assert excinfo.value.code == 2 + reported = pipeline.logger.error.call_args[0][-1] + assert "chr1" in reported and "chr2" not in reported + + def test_a_contig_missing_from_the_fai_is_not_reported_twice(self, cls): + """`check_pangenome_reference_fasta` already reported it""" + pipeline = self._pipeline(cls) + del pipeline.fai_data["chr2"] + del pipeline.pop_vcf_contigs["chr2"] + pipeline.validate_pangenome_contig_lengths() + pipeline.logger.error.assert_not_called() + assert pipeline.logger.info.call_args[0][1] == 1 + + +@pytest.mark.parametrize( + "cls", PANGENOME_PIPELINES, ids=lambda c: c.__name__ +) +class TestPopVcfHeader: + """`load_pop_vcf_header` and `check_pop_vcf_id` in both pipelines""" + + HEADER = ( + "##fileformat=VCFv4.2\n" + "##SentieonVcfID=population-test-20260101\n" + "##contig=\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + ) + + def _pipeline(self, cls, **attributes): + pipeline = cls() + pipeline.logger = MagicMock() + for key, value in attributes.items(): + setattr(pipeline, key, value) + return pipeline + + def _pop_vcf(self, tmp_path, text=None): + pop_vcf = tmp_path / "pop.vcf.gz" + pop_vcf.write_text(self.HEADER if text is None else text) + return pop_vcf + + def test_the_header_and_the_contigs_are_parsed(self, cls, tmp_path): pipeline = self._pipeline(cls) - pipeline.pop_vcf_contigs = {} + pipeline.load_pop_vcf_header(self._pop_vcf(tmp_path)) + assert pipeline.pop_vcf_contigs == {"chr1": 1000} + assert pipeline.pop_vcf_header is not None + + def test_no_pop_vcf_leaves_both_empty(self, cls): + pipeline = self._pipeline(cls) + pipeline.load_pop_vcf_header(None) + assert pipeline.pop_vcf_contigs == {} + assert pipeline.pop_vcf_header is None + + def test_a_dry_run_falls_back_on_an_unparsable_pop_vcf( + self, cls, tmp_path + ): + """Unit tests and dry runs supply placeholder files""" + pipeline = self._pipeline(cls, dry_run=True) + placeholder = tmp_path / "placeholder.vcf.gz" + placeholder.touch() + pipeline.load_pop_vcf_header(placeholder) + assert pipeline.pop_vcf_header is None + assert pipeline.pop_vcf_contigs == {} + pipeline.logger.debug.assert_called() + pipeline.logger.error.assert_not_called() + + def test_an_unparsable_pop_vcf_exits_outside_a_dry_run( + self, cls, tmp_path + ): + pipeline = self._pipeline(cls) + placeholder = tmp_path / "placeholder.vcf.gz" + placeholder.touch() + with pytest.raises(SystemExit) as excinfo: + pipeline.load_pop_vcf_header(placeholder) + assert excinfo.value.code == 2 + + def test_a_missing_pop_vcf_exits_outside_a_dry_run(self, cls, tmp_path): + pipeline = self._pipeline(cls) + with pytest.raises(SystemExit) as excinfo: + pipeline.load_pop_vcf_header(tmp_path / "absent.vcf.gz") + assert excinfo.value.code == 2 + + def test_a_matching_vcf_id_passes(self, cls, tmp_path): + pipeline = self._pipeline(cls) + pipeline.load_pop_vcf_header(self._pop_vcf(tmp_path)) + pipeline.check_pop_vcf_id("population-test-20260101") + pipeline.logger.error.assert_not_called() + + def test_a_mismatching_vcf_id_exits_in_a_dry_run(self, cls, tmp_path): + """The ID check no longer waits for a real run""" + pipeline = self._pipeline(cls, dry_run=True) + pipeline.load_pop_vcf_header(self._pop_vcf(tmp_path)) with pytest.raises(SystemExit) as excinfo: - pipeline.validate_grch38_contigs() + pipeline.check_pop_vcf_id("population-other-20240101") assert excinfo.value.code == 2 - def test_a_dry_run_skips_the_pop_vcf_check(self, cls): + def test_an_unread_header_skips_the_vcf_id_check(self, cls): + pipeline = self._pipeline(cls, dry_run=True) + pipeline.check_pop_vcf_id("population-test-20260101") + pipeline.logger.error.assert_not_called() + pipeline.logger.debug.assert_called() + + def test_an_unparsable_header_does_not_raise_in_a_dry_run( + self, cls, tmp_path + ): + """The whole fall-back path, as a `--dry_run` walks it""" pipeline = self._pipeline(cls, dry_run=True) - pipeline.pop_vcf_contigs = {} - pipeline.validate_grch38_contigs() + with patch( + "sentieon_cli.base_pangenome.read_vcf_header", + side_effect=VcfHeaderError("no VCF header"), + ): + pipeline.load_pop_vcf_header(tmp_path / "pop.vcf.gz") + pipeline.check_pop_vcf_id("population-test-20260101") + pipeline.validate_pangenome_contig_lengths() + pipeline.logger.error.assert_not_called() if __name__ == "__main__": diff --git a/tests/unit/test_sentieon_pangenome_kmc_integration.py b/tests/unit/test_sentieon_pangenome_kmc_integration.py index 9b6d8b8..b4f5fe7 100644 --- a/tests/unit/test_sentieon_pangenome_kmc_integration.py +++ b/tests/unit/test_sentieon_pangenome_kmc_integration.py @@ -22,8 +22,9 @@ def test_validate_kmc_patch_success(self): pipeline.pop_vcf = "pop" pipeline.gbz = "ref.grch38.gbz" pipeline.hapl = "ref.grch38.hapl" - # The mock graph files cannot be parsed; a dry run falls back to - # the default pangenome reference + # The mock graph and pop VCF files cannot be parsed; a dry run + # falls back to the default pangenome reference and to an + # empty set of pop VCF contigs pipeline.dry_run = True # Mock dependencies @@ -39,7 +40,6 @@ def test_validate_kmc_patch_success(self): with patch("sentieon_cli.util.check_version", return_value=True), \ patch("sentieon_cli.sentieon_pangenome.parse_fai", return_value={}), \ patch("sentieon_cli.sentieon_pangenome.determine_shards_from_fai", return_value=[]), \ - patch("sentieon_cli.sentieon_pangenome.vcf_contigs", return_value={}), \ patch("sentieon_cli.sentieon_pangenome.check_kmc_patch", return_value=True) as mock_check, \ patch("sys.exit") as mock_exit: @@ -69,7 +69,6 @@ def test_validate_kmc_patch_failure(self): with patch("sentieon_cli.util.check_version", return_value=True), \ patch("sentieon_cli.sentieon_pangenome.parse_fai", return_value={}), \ patch("sentieon_cli.sentieon_pangenome.determine_shards_from_fai", return_value=[]), \ - patch("sentieon_cli.sentieon_pangenome.vcf_contigs", return_value={}), \ patch("sentieon_cli.sentieon_pangenome.check_kmc_patch", return_value=False) as mock_check, \ patch("sys.exit") as mock_exit: diff --git a/tests/unit/test_vcf_header.py b/tests/unit/test_vcf_header.py new file mode 100644 index 0000000..859455a --- /dev/null +++ b/tests/unit/test_vcf_header.py @@ -0,0 +1,175 @@ +""" +Unit tests for the pure-Python VCF header reader +""" + +import gzip +import pathlib +import struct +import zlib + +import pytest + +from sentieon_cli.shard import parse_vcf_contigs, vcf_contigs +from sentieon_cli.util import ( + VcfHeaderError, + read_vcf_header, + vcf_id, + vcf_id_from_header, +) + +HEADER = """##fileformat=VCFv4.2 +##SentieonVcfID=population-test-20260101 +##contig= +##contig= +##contig= +#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO +""" + +BODY = "chr1\t1\t.\tA\tC\t50\tPASS\t.\n#not-a-header-line\n" + + +def bgzf(data: bytes) -> bytes: + """A single BGZF block, the way a `bgzip`-compressed VCF starts. + + BGZF is gzip with an extra `BC` subfield, so the magic bytes the + reader sniffs are the same and `gzip` decompresses it unchanged. + """ + deflater = zlib.compressobj(6, zlib.DEFLATED, -zlib.MAX_WBITS) + deflated = deflater.compress(data) + deflater.flush() + bsize = 18 + len(deflated) + 8 - 1 + header = struct.pack( + " pathlib.Path: + """A plain-text VCF, named `.vcf.gz` to defeat the extension""" + path = tmp_path / "plain.vcf.gz" + path.write_text(HEADER + BODY) + return path + + +@pytest.fixture +def gzip_vcf(tmp_path: pathlib.Path) -> pathlib.Path: + """A gzip-compressed VCF, named `.vcf` to defeat the extension""" + path = tmp_path / "compressed.vcf" + with gzip.open(path, "wt") as fh: + fh.write(HEADER + BODY) + return path + + +@pytest.fixture +def bgzf_vcf(tmp_path: pathlib.Path) -> pathlib.Path: + """A BGZF-compressed VCF, as `bgzip` would write it""" + path = tmp_path / "bgzf.vcf.gz" + path.write_bytes(bgzf((HEADER + BODY).encode())) + return path + + +class TestReadVcfHeader: + """`read_vcf_header` decompresses by magic bytes, not by name""" + + def test_a_plain_text_vcf_is_read(self, plain_vcf): + assert read_vcf_header(plain_vcf) == HEADER.splitlines() + + def test_a_gzip_vcf_is_read(self, gzip_vcf): + assert read_vcf_header(gzip_vcf) == HEADER.splitlines() + + def test_a_bgzf_vcf_is_read(self, bgzf_vcf): + assert read_vcf_header(bgzf_vcf) == HEADER.splitlines() + + def test_the_body_is_never_read(self, plain_vcf): + """Reading stops at `#CHROM`, so a `#` line below it is not a + header line""" + header = read_vcf_header(plain_vcf) + assert header[-1].startswith("#CHROM") + assert "#not-a-header-line" not in header + + def test_a_header_without_chrom_stops_at_the_body(self, tmp_path): + path = tmp_path / "no_chrom.vcf" + path.write_text("##fileformat=VCFv4.2\nchr1\t1\t.\tA\tC\n") + assert read_vcf_header(path) == ["##fileformat=VCFv4.2"] + + def test_a_file_with_no_header_raises(self, tmp_path): + path = tmp_path / "headerless.vcf" + path.write_text("chr1\t1\t.\tA\tC\t50\tPASS\t.\n") + with pytest.raises(VcfHeaderError): + read_vcf_header(path) + + def test_an_empty_file_raises(self, tmp_path): + path = tmp_path / "empty.vcf.gz" + path.touch() + with pytest.raises(VcfHeaderError): + read_vcf_header(path) + + def test_a_missing_file_raises_oserror(self, tmp_path): + with pytest.raises(OSError): + read_vcf_header(tmp_path / "absent.vcf.gz") + + def test_truncated_gzip_raises(self, tmp_path, gzip_vcf): + path = tmp_path / "truncated.vcf.gz" + path.write_bytes(gzip_vcf.read_bytes()[:20]) + with pytest.raises((VcfHeaderError, OSError)): + read_vcf_header(path) + + +class TestVcfContigs: + """`vcf_contigs` keeps its return type and its `None` lengths""" + + def test_contigs_and_lengths_are_parsed(self, bgzf_vcf): + assert vcf_contigs(bgzf_vcf) == { + "chr1": 248956422, + "chr2": 242193529, + "chrUn_KI270302v1": None, + } + + def test_a_contig_without_a_length_is_none(self): + contigs = parse_vcf_contigs(["##contig="]) + assert contigs == {"chr1": None} + + def test_a_quoted_assembly_does_not_confuse_the_parser(self): + contigs = parse_vcf_contigs( + ['##contig='] + ) + assert contigs == {"chr1": 10} + + def test_a_missing_file_returns_an_empty_dict(self, tmp_path): + assert vcf_contigs(tmp_path / "absent.vcf.gz") == {} + + def test_a_headerless_file_returns_an_empty_dict(self, tmp_path): + path = tmp_path / "headerless.vcf" + path.write_text("chr1\t1\t.\tA\tC\n") + assert vcf_contigs(path) == {} + + +class TestVcfId: + """`vcf_id` reports `None` the same way `vcf_contigs` reports `{}`""" + + def test_the_id_is_read(self, gzip_vcf): + assert vcf_id(gzip_vcf) == "population-test-20260101" + + def test_a_header_without_an_id(self, tmp_path): + path = tmp_path / "no_id.vcf" + path.write_text("##fileformat=VCFv4.2\n") + assert vcf_id(path) is None + + def test_a_missing_file_returns_none(self, tmp_path): + assert vcf_id(tmp_path / "absent.vcf.gz") is None + + def test_an_empty_file_returns_none(self, tmp_path): + path = tmp_path / "empty.vcf.gz" + path.touch() + assert vcf_id(path) is None + + def test_vcf_id_from_header(self): + assert vcf_id_from_header(HEADER.splitlines()) == ( + "population-test-20260101" + ) + assert vcf_id_from_header([]) is None + + +if __name__ == "__main__": + pytest.main([__file__])