diff --git a/sentieon_cli/base_pangenome.py b/sentieon_cli/base_pangenome.py index bf7f90f..5b2f328 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 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 detect_reference_build, parse_vcf_contigs from .stages.base import StageContext from .stages.cnv import PangenomeCNVResult, PangenomeCNVStage from .stages.expansion import ExpansionHunterResult, ExpansionHunterStage @@ -34,7 +42,27 @@ 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 +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): @@ -173,6 +201,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 +244,20 @@ 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. + # `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 self.sr_alignments: List[pathlib.Path] = [] @@ -230,6 +289,392 @@ 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 load_pop_vcf_header(self, pop_vcf: Optional[pathlib.Path]) -> None: + """Parse the `--pop_vcf` header. + + 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. + """ + self.pop_vcf_header = None + self.pop_vcf_contigs = {} + if pop_vcf is None: + return + 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, + ) + 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) + + 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( + "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) + + 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 + + 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( + "%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. + + `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/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 78a1ab3..d3ae092 100644 --- a/sentieon_cli/hybrid_pangenome.py +++ b/sentieon_cli/hybrid_pangenome.py @@ -27,11 +27,9 @@ from .job import Job from .logging import get_logger from .shard import ( - GRCH38_CONTIGS, detect_reference_build, determine_shards_from_fai, parse_fai, - vcf_contigs, ) from .shell_pipeline import Command, Pipeline from .stages.alignment import ( @@ -62,7 +60,6 @@ path_arg, require_versions, total_memory, - vcf_id, ) HYBRID_PANGENOME_MIN_VERSIONS = { @@ -151,16 +148,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 +172,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 +202,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 +212,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,10 +239,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]] = {} - 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) + # Before `validate_bundle`, which picks the bundle's + # `extract..model` member by the resolved name + self.resolve_pangenome_reference() + self.load_pop_vcf_header(self.pop_vcf) self.validate_bundle() self.validate_output_vcf() @@ -322,51 +297,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_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 @@ -550,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""" @@ -860,9 +781,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 +802,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 +1155,7 @@ def build_haplotypes_job( "--include-reference", "--diploid-sampling", "--set-reference", - self.pangenome_ref_name, + self.ref_name(), ], ), "vg-haplotypes", @@ -1242,7 +1173,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 +1200,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 +1265,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 +1326,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..ed445e8 100644 --- a/sentieon_cli/sentieon_pangenome.py +++ b/sentieon_cli/sentieon_pangenome.py @@ -51,13 +51,10 @@ path_arg, require_versions, total_memory, - vcf_id, ) from .shard import ( - GRCH38_CONTIGS, determine_shards_from_fai, parse_fai, - vcf_contigs, ) SENT_PANGENOME_MIN_VERSIONS = { @@ -113,22 +110,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 +119,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 +144,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,10 +170,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]] = {} - 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) + # Before `validate_bundle`, which picks the bundle's + # `extract..model` member by the resolved name + self.resolve_pangenome_reference() + self.load_pop_vcf_header(self.pop_vcf) self.validate_bundle() self.validate_fastq_rg() @@ -253,51 +222,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_pangenome_contig_lengths() def validate_segdup(self) -> None: if self.segdup_caller is None: @@ -354,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" @@ -427,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): @@ -589,11 +509,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 +533,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 +641,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 +757,7 @@ def build_haplotypes_job( "--include-reference", "--diploid-sampling", "--set-reference", - self.pangenome_ref_name, + self.ref_name(), ], ), "vg-haplotypes", @@ -845,7 +775,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 +829,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/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/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 0000000..d24bff4 Binary files /dev/null and b/tests/data/pangenome/tiny.gbz differ 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 0000000..7bf79bb Binary files /dev/null and b/tests/data/pangenome/tiny.hapl differ diff --git a/tests/data/pangenome/tiny_chm13.gbz b/tests/data/pangenome/tiny_chm13.gbz new file mode 100644 index 0000000..f7a304c Binary files /dev/null and b/tests/data/pangenome/tiny_chm13.gbz differ 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 0000000..453f460 Binary files /dev/null and b/tests/data/pangenome/tiny_no_backbone.gbz differ diff --git a/tests/data/pangenome/tiny_no_backbone.gfa b/tests/data/pangenome/tiny_no_backbone.gfa new file mode 100644 index 0000000..753d81d --- /dev/null +++ b/tests/data/pangenome/tiny_no_backbone.gfa @@ -0,0 +1,30 @@ +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 20 >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 0000000..a19ffce Binary files /dev/null and b/tests/data/pangenome/tiny_no_ref_tag.gbz differ diff --git a/tests/data/pangenome/tiny_no_ref_tag.gfa b/tests/data/pangenome/tiny_no_ref_tag.gfa new file mode 100644 index 0000000..ceacf56 --- /dev/null +++ b/tests/data/pangenome/tiny_no_ref_tag.gfa @@ -0,0 +1,28 @@ +H VN:Z:1.1 +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_sampled.gbz b/tests/data/pangenome/tiny_sampled.gbz new file mode 100644 index 0000000..c0e19be Binary files /dev/null and b/tests/data/pangenome/tiny_sampled.gbz differ 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..ba5a546 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,7 +17,13 @@ from sentieon_cli.dnascope import DNAscopePipeline from sentieon_cli.dnascope_longread import DNAscopeLRPipeline from sentieon_cli.dnascope_hybrid import DNAscopeHybridPipeline -from sentieon_cli.util import set_bwt_max_mem +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.util import VcfHeaderError, set_bwt_max_mem from tests.utils.test_helpers import create_mock_args @@ -514,5 +521,490 @@ 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 TestPangenomeContigLengthChecks: + """`BasePangenome.validate_pangenome_contig_lengths` in both pipelines. + + 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_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_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_mismatching_length_exits(self, cls): + pipeline = self._pipeline(cls) + pipeline.pop_vcf_contigs["chr2"] = 2001 + 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 "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.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.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.check_pop_vcf_id("population-other-20240101") + assert excinfo.value.code == 2 + + 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) + 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__": 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..b4f5fe7 100644 --- a/tests/unit/test_sentieon_pangenome_kmc_integration.py +++ b/tests/unit/test_sentieon_pangenome_kmc_integration.py @@ -21,6 +21,11 @@ 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 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 pipeline.logger = MagicMock() @@ -35,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: @@ -56,13 +60,15 @@ 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() 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__])