Skip to content

Gbz name validation #38

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
DonFreed merged 2 commits into from
Sep 13, 2026
Merged

Gbz name validation #38

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5d39e3a
Detect the pangenome reference from the GBZ metadata
DonFreed Sep 12, 2026
fb193d3
Check pop VCF contig lengths against the reference for any backbone
DonFreed Sep 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
449 changes: 447 additions & 2 deletions sentieon_cli/base_pangenome.py
View file
Open in desktop

Large diffs are not rendered by default.

61 changes: 59 additions & 2 deletions sentieon_cli/command_strings.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 = [
Expand Down
2 changes: 1 addition & 1 deletion sentieon_cli/dnascope_hybrid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion sentieon_cli/dnascope_longread.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
117 changes: 24 additions & 93 deletions sentieon_cli/hybrid_pangenome.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 (
Expand Down Expand Up @@ -62,7 +60,6 @@
path_arg,
require_versions,
total_memory,
vcf_id,
)

HYBRID_PANGENOME_MIN_VERSIONS = {
Expand Down Expand Up @@ -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 "
Expand All @@ -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",
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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.<pangenome_ref_name>.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()
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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"""
Expand Down Expand Up @@ -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(
Expand All @@ -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(
Expand Down Expand Up @@ -1224,7 +1155,7 @@ def build_haplotypes_job(
"--include-reference",
"--diploid-sampling",
"--set-reference",
self.pangenome_ref_name,
self.ref_name(),
],
),
"vg-haplotypes",
Expand All @@ -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,
Expand All @@ -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(
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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(
Expand Down
Loading
Loading