From 58c9c4d6a6150c0ccf16d73b49232bea3a8b049c Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sat, 15 Aug 2026 15:27:30 -0700 Subject: [PATCH 01/10] Add the hybrid-pangenome pipeline Combine the human pangenome with sample short-read fastq and aligned long-read data for variant calling: k-mer counting across both read sets, personalized diploid graph construction with vg, graph updates with PGHapUpdateAlgo using the long-read alignments and LongReadSV calls, re-alignment of extracted short reads to the personalized pangenome with liftover, short-read dedup and metrics, and PangenomeSV + DNAscope calling with LR readgroup attributes. The dedup and metrics job builders move from SentieonPangenome to BasePangenome so both pangenome pipelines share them; the emitted dnascope-pangenome commands are unchanged. Co-Authored-By: Claude Fable 5 --- sentieon_cli/__init__.py | 7 + sentieon_cli/base_pangenome.py | 145 ++++- sentieon_cli/command_strings.py | 116 ++++ sentieon_cli/driver.py | 20 + sentieon_cli/hybrid_pangenome.py | 911 ++++++++++++++++++++++++++++ sentieon_cli/sentieon_pangenome.py | 144 +---- tests/unit/test_hybrid_pangenome.py | 474 +++++++++++++++ 7 files changed, 1676 insertions(+), 141 deletions(-) create mode 100644 sentieon_cli/hybrid_pangenome.py create mode 100644 tests/unit/test_hybrid_pangenome.py diff --git a/sentieon_cli/__init__.py b/sentieon_cli/__init__.py index 59c9971..bffbb96 100644 --- a/sentieon_cli/__init__.py +++ b/sentieon_cli/__init__.py @@ -4,6 +4,7 @@ from .dnascope import DNAscopePipeline from .dnascope_hybrid import DNAscopeHybridPipeline from .dnascope_longread import DNAscopeLRPipeline +from .hybrid_pangenome import HybridPangenome from .job import Job from .sentieon_pangenome import SentieonPangenome from .util import __version__ @@ -76,6 +77,12 @@ def main(): pipeline.add_arguments(dnascope_pangenome_subparser) dnascope_pangenome_subparser.set_defaults(pipeline=pipeline.main) + # Hybrid pangenome + pipeline = HybridPangenome() + hybrid_pangenome_subparser = subparsers.add_parser("hybrid-pangenome") + pipeline.add_arguments(hybrid_pangenome_subparser) + hybrid_pangenome_subparser.set_defaults(pipeline=pipeline.main) + args = parser.parse_args() args.loglevel = resolve_loglevel(args) # Job ids must be unique for the whole run, which may execute more than diff --git a/sentieon_cli/base_pangenome.py b/sentieon_cli/base_pangenome.py index 7284de6..c126bbd 100644 --- a/sentieon_cli/base_pangenome.py +++ b/sentieon_cli/base_pangenome.py @@ -6,13 +6,28 @@ from enum import Enum import json import pathlib -from typing import List, Optional +import sys +from typing import List, Optional, Tuple from importlib.resources import files from . import command_strings as cmds +from .driver import ( + AlignmentStat, + BaseDistributionByCycle, + CoverageMetrics, + Dedup, + Driver, + GCBias, + InsertSizeMetricAlgo, + LocusCollector, + MeanQualityByCycle, + QualDistribution, + WgsMetricsAlgo, +) from .job import Job from .pipeline import BasePipeline +from .shell_pipeline import Command, Pipeline from .util import path_arg @@ -155,6 +170,134 @@ def build_ploidy_job( ) return ploidy_job + def build_dedup_job( + self, + output_bam, + input_bam: List[pathlib.Path], + tag: str, + left_align_rgid: Optional[str] = None, + metrics: Optional[pathlib.Path] = None, + ) -> Tuple[Job, Job]: + """Build deduplication job""" + score_file = self.tmp_dir.joinpath(f"sample-{tag}-score.txt.gz") + + read_filters = [] + if left_align_rgid: + read_filters.append( + f"IndelLeftAlignReadTransform,rgid={left_align_rgid}" + ) + + # LocusCollector + Dedup + driver = Driver( + reference=self.reference, + thread_count=self.cores, + input=input_bam, + read_filter=read_filters, + ) + driver.add_algo(LocusCollector(score_file)) + + lc_job = Job( + Pipeline(Command(*driver.build_cmd())), + f"locuscollector-{tag}", + self.cores, + task_name="dedup", + ) + + driver2 = Driver( + reference=self.reference, + thread_count=self.cores, + input=input_bam, + read_filter=read_filters, + ) + driver2.add_algo(Dedup(output_bam, score_file, metrics=metrics)) + + dedup_job = Job( + Pipeline(Command(*driver2.build_cmd())), + f"dedup-{tag}", + self.cores, + task_name="dedup", + ) + + return lc_job, dedup_job + + def build_metrics_job( + self, + sample_input: List[pathlib.Path], + ) -> Tuple[Job, Job]: + """Build a metrics job""" + if not self.output_vcf: + self.logger.error("output_vcf is required") + sys.exit(2) + + # Create the metrics directory + sample_name = self.output_vcf.name.replace(".vcf.gz", "") + metric_base = sample_name + ".txt" + metrics_dir = pathlib.Path( + str(self.output_vcf).replace(".vcf.gz", "_metrics") + ) + if not self.dry_run: + metrics_dir.mkdir(exist_ok=True) + + is_metrics = metrics_dir.joinpath(metric_base + ".insert_size.txt") + mqbc_metrics = metrics_dir.joinpath( + metric_base + ".mean_qual_by_cycle.txt" + ) + bdbc_metrics = metrics_dir.joinpath( + metric_base + ".base_distribution_by_cycle.txt" + ) + qualdist_metrics = metrics_dir.joinpath( + metric_base + ".qual_distribution.txt" + ) + as_metrics = metrics_dir.joinpath(metric_base + ".alignment_stat.txt") + coverage_metrics = metrics_dir.joinpath("coverage") + + # WGS metrics + wgs_metrics = metrics_dir.joinpath(metric_base + ".wgs.txt") + gc_metrics = metrics_dir.joinpath(metric_base + ".gc_bias.txt") + gc_summary = metrics_dir.joinpath(metric_base + ".gc_bias_summary.txt") + + driver = Driver( + reference=self.reference, + thread_count=self.cores, + input=sample_input, + ) + + driver.add_algo(InsertSizeMetricAlgo(is_metrics)) + driver.add_algo(MeanQualityByCycle(mqbc_metrics)) + driver.add_algo(BaseDistributionByCycle(bdbc_metrics)) + driver.add_algo(QualDistribution(qualdist_metrics)) + driver.add_algo(AlignmentStat(as_metrics)) + driver.add_algo(GCBias(gc_metrics, summary=gc_summary)) + driver.add_algo(WgsMetricsAlgo(wgs_metrics, include_unpaired="true")) + driver.add_algo(CoverageMetrics(coverage_metrics)) + + metrics_job = Job( + Pipeline(Command(*driver.build_cmd())), + "metrics", + 0, + task_name="metrics", + ) + + rehead_script = pathlib.Path( + str( + files("sentieon_cli.scripts").joinpath("rehead_wgs_metrics.py") + ) + ) + rehead_job = Job( + Pipeline( + Command( + sys.executable, + str(rehead_script), + "--metrics_file", + str(wgs_metrics), + ) + ), + "Rehead metrics", + 0, + task_name="metrics", + ) + return (metrics_job, rehead_job) + def get_sex(self, ploidy_json: pathlib.Path) -> None: """Retrieve the sample sex""" if self.dry_run: diff --git a/sentieon_cli/command_strings.py b/sentieon_cli/command_strings.py index 2595fe1..e09877a 100644 --- a/sentieon_cli/command_strings.py +++ b/sentieon_cli/command_strings.py @@ -1415,6 +1415,122 @@ def cmd_minimap2_lift( return Pipeline(mm2_cmd, lift_cmd, sort_cmd) +def cmd_hybrid_kmc( + output_prefix: pathlib.Path, + fastq: List[pathlib.Path], + aln: List[pathlib.Path], + reference: pathlib.Path, + tmp_dir: pathlib.Path, + k: int = 29, + memory: int = 30, + threads: int = 1, + unzip: str = "gzip", +) -> Pipeline: + """Count k-mers across fastq and aligned reads. + + KMC accepts a single input format per run, so the fastq files are + converted to FASTA and the aligned reads are extracted with + `samtools fasta`, and one FASTA stream is fed to the patched KMC + through stdin. + """ + cat_args: List[Union[str, InputProcSub]] = [] + if fastq: + fq_fasta = Pipeline( + Command( + unzip, + "-dc", + *[str(x) for x in fastq], + ), + Command("awk", 'NR%4==1{print ">"substr($0,2)} NR%4==2{print}'), + ) + cat_args.append(InputProcSub(fq_fasta)) + for aln_file in aln: + cat_args.append( + InputProcSub( + Pipeline( + Command( + "samtools", + "fasta", + "--reference", + str(reference), + "-@", + str(threads), + str(aln_file), + ) + ) + ) + ) + cat_cmd = Command("cat", *cat_args) + kmc_cmd = Command( + "kmc", + f"-k{k}", + f"-m{memory}", + "-okff", + f"-t{threads}", + "-fa", + "/dev/stdin", + str(output_prefix), + str(tmp_dir), + ) + return Pipeline(cat_cmd, kmc_cmd) + + +# Extract graph update regions from the phase sets of LongReadSV calls. +LONGREAD_SV_BED_AWK = """!/^#/ { + n=split($10,a,":"); + ps=a[n]; + if (ps != "." && ps ~ /^chr[^_]+_[0-9]+_[0-9]+$/) { + split(ps,b,"_"); + print b[1], b[2], b[3]; + } +}""" + + +def cmd_longread_sv_bed( + out_bed: pathlib.Path, + sv_vcf: pathlib.Path, +) -> Pipeline: + """Generate a BED file of graph update regions from LongReadSV calls""" + zcat_cmd = Command("zcat", str(sv_vcf)) + awk_cmd = Command("awk", "-F\t", LONGREAD_SV_BED_AWK, "OFS=\t") + sort_cmd = Command("sort", "-k1,1", "-k2,2n") + merge_cmd = Command("bedtools", "merge") + return Pipeline( + zcat_cmd, + awk_cmd, + sort_cmd, + merge_cmd, + file_output=out_bed, + ) + + +def cmd_pgutil_gfa2fa( + out_fasta: pathlib.Path, + ref_fai: pathlib.Path, + gfa_file: pathlib.Path, +) -> Pipeline: + """Generate FASTA sequences from a pangenome graph""" + cmd = Command( + "sentieon", + "pgutil", + "gfa2fa", + "-F", + str(ref_fai), + "-g", + str(gfa_file), + "-o", + str(out_fasta), + ) + return Pipeline(cmd) + + +def cmd_samtools_faidx( + fasta: pathlib.Path, +) -> Pipeline: + """Index a FASTA file""" + return Pipeline(Command("samtools", "faidx", str(fasta))) + + def cmd_bcftools_merge_trim( shard_vcf: pathlib.Path, raw_vcf: pathlib.Path, diff --git a/sentieon_cli/driver.py b/sentieon_cli/driver.py index c2bbb7f..e39fe80 100644 --- a/sentieon_cli/driver.py +++ b/sentieon_cli/driver.py @@ -542,6 +542,26 @@ def __init__( self.prefix = prefix +class PGHapUpdateAlgo(BaseAlgo): + """algo PGHapUpdateAlgo""" + + name = "PGHapUpdateAlgo" + + def __init__( + self, + output: pathlib.Path, + gfa_file: pathlib.Path, + target_bed: Optional[pathlib.Path] = None, + min_map_qual: Optional[int] = None, + prefix: Optional[str] = None, + ): + self.output = output + self.gfa_file = gfa_file + self.target_bed = target_bed + self.min_map_qual = min_map_qual + self.prefix = prefix + + class BaseDriver: """A base class for the Sentieon driver""" diff --git a/sentieon_cli/hybrid_pangenome.py b/sentieon_cli/hybrid_pangenome.py new file mode 100644 index 0000000..24ef897 --- /dev/null +++ b/sentieon_cli/hybrid_pangenome.py @@ -0,0 +1,911 @@ +""" +The Sentieon hybrid-pangenome pipeline + +Combines the human pangenome with sample short-read and long-read data +for highly accurate variant calling. +""" + +import argparse +import copy +import json +import pathlib +import shutil +import sys +from typing import Dict, List, Optional, Set + +import packaging.version + +from . import command_strings as cmds +from .archive import ar_load +from .base_pangenome import BasePangenome +from .dag import DAG +from .driver import ( + DNAModelApply, + DNAscope, + Driver, + LongReadSV, + PangenomeSV, + PGHapUpdateAlgo, +) +from .job import Job +from .logging import get_logger +from .shard import ( + GRCH38_CONTIGS, + determine_shards_from_fai, + parse_fai, + vcf_contigs, +) +from .shell_pipeline import Command, Pipeline +from .transfer import build_transfer_jobs +from .util import ( + __version__, + check_kmc_patch, + check_version, + parse_rg_line, + path_arg, + total_memory, + vcf_id, +) + +HYBRID_PANGENOME_MIN_VERSIONS = { + "kmc": None, + "sentieon driver": packaging.version.Version("202503.04"), + "vg": None, + "bcftools": packaging.version.Version("1.22"), + "samtools": packaging.version.Version("1.16"), + "bedtools": None, +} + +# LongReadSV settings for finding graph update regions +LONGREADSV_MIN_SV_SIZE = 20 + +# PangenomeSV settings +PANGENOME_SV_MIN_AF = 0.1 + +logger = get_logger(__name__) + + +class HybridPangenome(BasePangenome): + """The Sentieon hybrid-pangenome pipeline""" + + params = copy.deepcopy(BasePangenome.params) + params.update( + { + # Required arguments + "lr_aln": { + "nargs": "*", + "help": ( + "Long-read BAM or CRAM files aligned to the linear " + "reference genome." + ), + "required": True, + "type": path_arg(exists=True, is_file=True), + }, + "pop_vcf": { + "flags": ["--pop_vcf"], + "help": ( + "A VCF containing annotations for use with DNAModelApply." + ), + "type": path_arg(exists=True, is_file=True), + "required": True, + }, + "readgroup": { + "help": "Readgroup information for the fastq files.", + }, + # Additional arguments + "bed": { + "flags": ["-b", "--bed"], + "help": ( + "Region BED file. Supplying this file will limit " + "small-variant calling to the intervals inside the BED " + "file." + ), + "type": path_arg(exists=True, is_file=True), + }, + "pangenome_ref_name": { + "default": "GRCh38", + "help": "Reference name in the pangenome (GRCh38).", + }, + "rgsm": { + "help": ( + "Overwrite the SM tag of the input readgroups for " + "compatibility" + ), + }, + "skip_metrics": { + "help": "Skip metrics collection and multiQC", + "action": "store_true", + }, + "skip_multiqc": { + "help": "Skip multiQC report generation", + "action": "store_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", + }, + "skip_small_variants": { + "help": argparse.SUPPRESS, + "action": "store_true", + }, + "skip_svs": { + "help": argparse.SUPPRESS, + "action": "store_true", + }, + } + ) + + positionals = BasePangenome.positionals + + def __init__(self) -> None: + super().__init__() + self.lr_aln: List[pathlib.Path] = [] + self.pop_vcf: Optional[pathlib.Path] = None + self.readgroup: Optional[str] = None + self.bed: Optional[pathlib.Path] = None + self.pangenome_ref_name = "GRCh38" + self.rgsm: Optional[str] = None + self.extract_model_name = "extract.model" + self.skip_metrics = False + self.skip_multiqc = False + 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 + + def validate(self) -> None: + """Validate pipeline inputs""" + self.validate_ref() + self.fai_data = parse_fai(pathlib.Path(str(self.reference) + ".fai")) + 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) + + self.validate_bundle() + self.validate_output_vcf() + + if not self.r1_fastq or not self.readgroup: + self.logger.error( + "Please supply the short reads with the `--r1_fastq`, " + "`--r2_fastq`, and `--readgroup` arguments" + ) + sys.exit(2) + if len(self.r1_fastq) != len(self.r2_fastq): + self.logger.error( + "The number of input `--r1_fastq` files does not equal the " + "number of `--r2_fastq` files" + ) + sys.exit(2) + if not self.lr_aln: + self.logger.error( + "Please supply the long-read alignments with the `--lr_aln` " + "argument" + ) + sys.exit(2) + + self.validate_bwa_index() + self.collect_readgroups() + self.validate_readgroups() + + if not self.skip_version_check: + for cmd, min_version in HYBRID_PANGENOME_MIN_VERSIONS.items(): + if not check_version(cmd, min_version): + sys.exit(2) + + if not check_kmc_patch("kmc"): + self.logger.error( + "Error: The 'kmc' executable in the PATH does not " + "support reading from stdin. Please ensure " + "you are using the patched version of KMC from " + "https://github.com/Sentieon/KMC/releases." + ) + sys.exit(2) + + if self.bed is None: + self.logger.info( + "A BED file is recommended to avoid small-variant calling " + "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) + + def validate_bundle(self) -> None: + """Validate the model bundle""" + bundle_info_bytes = ar_load( + str(self.model_bundle) + "/bundle_info.json" + ) + if isinstance(bundle_info_bytes, list): + bundle_info_bytes = b"{}" + bundle_info = json.loads(bundle_info_bytes.decode()) + + req_version_s = bundle_info.get("minScriptVersion") + if req_version_s: + req_version = packaging.version.Version(req_version_s) + if req_version > packaging.version.Version(__version__): + self.logger.error( + "The model bundle requires version %s or later of the " + "sentieon-cli.", + req_version, + ) + sys.exit(2) + + bundle_pipeline = bundle_info.get("pipeline") + if bundle_pipeline and bundle_pipeline != "Hybrid pangenome": + self.logger.error("The model bundle is for a different pipeline.") + sys.exit(2) + + bundle_members = set(ar_load(str(self.model_bundle))) + + # Prefer a reference-specific extract model. Fall back to the generic + # 'extract.model' only for the default 'GRCh38' reference. + extract_candidate = f"extract.{self.pangenome_ref_name}.model" + if extract_candidate in bundle_members: + self.extract_model_name = extract_candidate + elif self.pangenome_ref_name == "GRCh38": + self.extract_model_name = "extract.model" + else: + self.extract_model_name = extract_candidate + + required_members = { + "bwa.model", + "dnascope.model", + "longreadsv.model", + "minimap2.model", + self.extract_model_name, + } + missing_members = required_members - bundle_members + if missing_members: + self.logger.error( + "Expected model files not found in the model bundle file: %s", + ", ".join(sorted(missing_members)), + ) + sys.exit(2) + + 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) + + def collect_readgroups(self) -> None: + """Collect readgroup tags from the inputs""" + assert self.readgroup is not None + try: + parsed_rg = parse_rg_line(self.readgroup.replace(r"\t", "\t")) + except ValueError as e: + self.logger.error( + "Invalid --readgroup value '%s': %s", self.readgroup, e + ) + sys.exit(2) + if not parsed_rg.get("ID"): + self.logger.error( + "Readgroup '%s' does not have a RGID tag", + self.readgroup, + ) + sys.exit(2) + if not parsed_rg.get("SM"): + self.logger.error( + "Readgroup '%s' does not have a RGSM tag", + self.readgroup, + ) + sys.exit(2) + self.fastq_readgroup: Dict[str, str] = parsed_rg + + self.lr_readgroups: List[List[Dict[str, str]]] = [] + for aln in self.lr_aln: + self.lr_readgroups.append([]) + for rg_line in cmds.get_rg_lines(aln, self.dry_run): + self.lr_readgroups[-1].append(parse_rg_line(rg_line)) + + def validate_readgroups(self) -> None: + """Confirm that all readgroups have a consistent SM tag""" + rg_sm = self.fastq_readgroup.get("SM") + for aln, aln_rgs in zip(self.lr_aln, self.lr_readgroups): + for aln_rg in aln_rgs: + if not aln_rg.get("ID"): + self.logger.error( + "Found a readgroup without an ID tag in '%s': %s", + aln, + str(aln_rg), + ) + sys.exit(2) + sm = aln_rg.get("SM") + if not sm: + self.logger.error( + "Found a readgroup without a SM tag in '%s': %s", + aln, + str(aln_rg), + ) + sys.exit(2) + if self.dry_run or self.rgsm: + continue + if sm != rg_sm: + self.logger.error( + "Input readgroup '%s' has a different RG-SM tag " + "from the `--readgroup` argument.\n" + "found='%s' expected='%s'. Please set the `--rgsm` " + "argument to override the SM tag in the input files", + str(aln_rg), + sm, + rg_sm, + ) + sys.exit(2) + self.sample_sm: str = self.rgsm if self.rgsm else str(rg_sm) + + def configure(self) -> None: + """Configure pipeline parameters""" + pass + + def find_unzip(self) -> str: + """The decompression tool for fastq input""" + unzip = "igzip" + if not shutil.which(unzip): + self.logger.info( + "igzip is recommended for decompression, but is not " + "available. Falling back to gzip." + ) + unzip = "gzip" + return unzip + + def build_dag(self) -> DAG: + """Build the DAG for the hybrid-pangenome pipeline""" + if not self.reference: + self.logger.error("reference is required") + sys.exit(2) + if not self.model_bundle: + self.logger.error("model_bundle is required") + sys.exit(2) + if not self.output_vcf: + self.logger.error("output_vcf is required") + sys.exit(2) + if not self.pop_vcf: + self.logger.error("pop_vcf is required") + sys.exit(2) + + self.logger.info("Building the hybrid-pangenome DAG") + dag = DAG() + + ref_fai = pathlib.Path(str(self.reference) + ".fai") + + # Output files + suffix = "bam" if self.bam_format else "cram" + out_bwa_aln = pathlib.Path( + str(self.output_vcf).replace(".vcf.gz", f"_bwa_deduped.{suffix}") + ) + out_lift_aln = pathlib.Path( + str(self.output_vcf).replace(".vcf.gz", f"_lift_deduped.{suffix}") + ) + sv_vcf = pathlib.Path( + str(self.output_vcf).replace(".vcf.gz", "_sv.vcf.gz") + ) + + # Intermediate file paths + bwa_bam = self.tmp_dir.joinpath("sample-bwa.bam") + lift_bam = self.tmp_dir.joinpath("sample-lift.bam") + ext_fastq = self.tmp_dir.joinpath("sample-extract.fq.gz") + kmer_prefix = self.tmp_dir.joinpath("sample") + kmer_file = pathlib.Path(str(kmer_prefix) + ".kff") + sample_pangenome = self.tmp_dir.joinpath("sample_pangenome.gbz") + hap_raw_gfa = self.tmp_dir.joinpath("sample-hap.raw.gfa") + pangenome_raw_gfa = self.tmp_dir.joinpath("sample-pangenome-raw.gfa") + pangenome_gfa = self.tmp_dir.joinpath("sample-pangenome.gfa") + pangenome_fasta = self.tmp_dir.joinpath("sample-pangenome.fa") + longread_sv_vcf = self.tmp_dir.joinpath("sample-longread-sv.vcf.gz") + sv_bed = self.tmp_dir.joinpath("sample-sv.bed") + raw_vcf = self.tmp_dir.joinpath("sample-dnascope.vcf.gz") + transfer_vcf = self.tmp_dir.joinpath("sample-dnascope_transfer.vcf.gz") + + total_mem_gb = total_memory() / (1024.0**3) + + # KMC k-mer counting across the short and long reads + kmc_job = self.build_hybrid_kmc_job(kmer_prefix) + dag.add_job(kmc_job) + haplotype_dependencies: Set[Job] = {kmc_job} + + # BWA alignment and extraction + bwa_job = self.build_alignment_job(bwa_bam, ext_fastq) + dag.add_job(bwa_job) + # Do not run vg-haplotypes with bwa in low-mem environments + if total_mem_gb < 70: + haplotype_dependencies.add(bwa_job) + + # vg haplotypes - create a sample-specific pangenome + haplotypes_job = self.build_haplotypes_job(sample_pangenome, kmer_file) + dag.add_job(haplotypes_job, haplotype_dependencies) + + # convert the sample pangenome to GFA + gfa_job = self.build_gfa_job(hap_raw_gfa, sample_pangenome) + dag.add_job(gfa_job, {haplotypes_job}) + + # Graph update without the SV BED + update_raw_job = self.build_graph_update_job( + pangenome_raw_gfa, + hap_raw_gfa, + name="graph-update-raw", + ) + dag.add_job(update_raw_job, {gfa_job}) + + # Call SVs from the long reads and collect graph update regions + longreadsv_job = self.build_longreadsv_job(longread_sv_vcf) + dag.add_job(longreadsv_job) + sv_bed_job = Job( + cmds.cmd_longread_sv_bed(sv_bed, longread_sv_vcf), + "longread-sv-bed", + 0, + task_name="pangenome-update", + ) + dag.add_job(sv_bed_job, {longreadsv_job}) + + # Graph update with the SV BED + update_job = self.build_graph_update_job( + pangenome_gfa, + pangenome_raw_gfa, + bed=sv_bed, + name="graph-update", + ) + dag.add_job(update_job, {update_raw_job, sv_bed_job}) + + # FASTA generation from the updated graph + gfa2fa_job = Job( + cmds.cmd_pgutil_gfa2fa(pangenome_fasta, ref_fai, pangenome_gfa), + "gfa2fa", + 0, + task_name="pangenome", + ) + dag.add_job(gfa2fa_job, {update_job}) + faidx_job = Job( + cmds.cmd_samtools_faidx(pangenome_fasta), + "faidx", + 0, + task_name="pangenome", + ) + dag.add_job(faidx_job, {gfa2fa_job}) + + # minimap2 alignment of the extracted reads with liftover + mm2_job = self.build_minimap2_lift_job( + lift_bam, + ext_fastq, + pangenome_fasta, + pangenome_gfa, + ) + dag.add_job(mm2_job, {bwa_job, faidx_job, update_job}) + + # Deduplicate the short-read alignments. The `--lr_aln` input is + # assumed to be deduplicated already. + # Emit Dedup metrics for the primary (bwa) short-read alignment so + # they land in the metrics directory scanned by MultiQC. + dedup_metrics: Optional[pathlib.Path] = None + if not self.skip_metrics: + metrics_dir = pathlib.Path( + str(self.output_vcf).replace(".vcf.gz", "_metrics") + ) + if not self.dry_run: + metrics_dir.mkdir(exist_ok=True) + sample_name = self.output_vcf.name.replace(".vcf.gz", "") + dedup_metrics = metrics_dir.joinpath( + sample_name + ".txt.dedup_metrics.txt" + ) + + bwa_lc_job, bwa_dedup_job = self.build_dedup_job( + out_bwa_aln, [bwa_bam], "bwa", metrics=dedup_metrics + ) + lift_lc_job, lift_dedup_job = self.build_dedup_job( + out_lift_aln, + [lift_bam], + "lift", + ) + dag.add_job(bwa_lc_job, {bwa_job}) + dag.add_job(bwa_dedup_job, {bwa_lc_job}) + dag.add_job(lift_lc_job, {mm2_job}) + dag.add_job(lift_dedup_job, {lift_lc_job}) + + # Alignment metrics from the deduplicated short reads + if not self.skip_metrics: + metrics_job, rehead_job = self.build_metrics_job( + [out_bwa_aln, out_lift_aln], + ) + dag.add_job(metrics_job, {bwa_dedup_job, lift_dedup_job}) + dag.add_job(rehead_job, {metrics_job}) + if not self.skip_multiqc: + multiqc_job = self.multiqc() + if multiqc_job: + dag.add_job(multiqc_job, {rehead_job}) + + # Variant calling with the original, lifted, and long reads + calling_bams = [out_bwa_aln, out_lift_aln] + list(self.lr_aln) + replace_rg = self.build_replace_rg() + calling_dependencies = {bwa_dedup_job, lift_dedup_job} + + if not self.skip_svs: + pangenomesv_job = self.build_pangenomesv_job( + sv_vcf, + calling_bams, + pangenome_gfa, + replace_rg, + ) + dag.add_job(pangenomesv_job, calling_dependencies | {update_job}) + + if self.skip_small_variants: + return dag + + dnascope_job = self.build_dnascope_job( + raw_vcf, + calling_bams, + replace_rg, + ) + dag.add_job(dnascope_job, calling_dependencies) + + # transfer annotations from the pop_vcf + transfer_target = ( + transfer_vcf if not self.skip_model_apply else self.output_vcf + ) + transfer_jobs, concat_job = build_transfer_jobs( + transfer_target, + self.pop_vcf, + raw_vcf, + self.tmp_dir, + self.shards, + self.pop_vcf_contigs, + self.fai_data, + self.dry_run, + self.cores, + ) + for job in transfer_jobs: + dag.add_job(job, {dnascope_job}) + dag.add_job(concat_job, set(transfer_jobs)) + + if not self.skip_model_apply: + apply_job = self.build_dnamodelapply_job( + transfer_vcf, self.output_vcf + ) + dag.add_job(apply_job, {concat_job}) + + return dag + + def build_hybrid_kmc_job(self, kmer_prefix: pathlib.Path) -> Job: + """Count k-mers across the short-read fastq and long-read + alignments""" + assert self.reference is not None + kmc_job = Job( + cmds.cmd_hybrid_kmc( + kmer_prefix, + self.r1_fastq + self.r2_fastq, + self.lr_aln, + self.reference, + self.tmp_dir, + memory=self.kmer_memory, + threads=self.cores, + unzip=self.find_unzip(), + ), + "kmc", + 0, # run in the background + task_name="kmer-counting", + ) + return kmc_job + + def build_alignment_job( + self, + sample_bam: pathlib.Path, + sample_fastq: pathlib.Path, + ) -> Job: + """Build the bwa alignment and read extraction job""" + assert self.reference is not None + assert self.model_bundle is not None + + rg = copy.deepcopy(self.fastq_readgroup) + rg["SM"] = self.sample_sm + rg["LR"] = "0" + bwa_job = Job( + cmds.cmd_bwa_extract( + sample_bam, + sample_fastq, + self.reference, + self.r1_fastq, + self.r2_fastq, + "@RG\\t" + "\\t".join([f"{x[0]}:{x[1]}" for x in rg.items()]), + self.model_bundle.joinpath(self.extract_model_name), + self.model_bundle.joinpath("bwa.model"), + self.cores, + unzip=self.find_unzip(), + ), + "bwa-extract", + self.cores, + task_name="alignment", + ) + return bwa_job + + def build_haplotypes_job( + self, output_gbz: pathlib.Path, kmer_file: pathlib.Path + ) -> Job: + """Build vg haplotypes job""" + assert self.hapl is not None + assert self.gbz is not None + + haplotypes_job = Job( + cmds.cmd_vg_haplotypes( + output_gbz, + kmer_file, + self.hapl, + self.gbz, + threads=self.cores, + xargs=[ + "--include-reference", + "--diploid-sampling", + "--set-reference", + self.pangenome_ref_name, + ], + ), + "vg-haplotypes", + self.cores, + task_name="pangenome", + ) + return haplotypes_job + + def build_gfa_job( + self, output_gfa: pathlib.Path, input_gbz: pathlib.Path + ) -> Job: + """Build vg convert to GFA job""" + gfa_job = Job( + cmds.cmd_vg_convert_gfa( + output_gfa, + input_gbz, + threads=self.cores, + reference_name=self.pangenome_ref_name, + ), + "vg-convert-gfa", + 0, + task_name="pangenome", + ) + return gfa_job + + def build_graph_update_job( + self, + out_gfa: pathlib.Path, + in_gfa: pathlib.Path, + bed: Optional[pathlib.Path] = None, + name: str = "graph-update", + ) -> Job: + """Update the personalized graph using the long-read alignments""" + driver = Driver( + reference=self.reference, + thread_count=self.cores, + input=list(self.lr_aln), + ) + driver.add_algo( + PGHapUpdateAlgo(out_gfa, gfa_file=in_gfa, target_bed=bed) + ) + return Job( + Pipeline(Command(*driver.build_cmd())), + name, + self.cores, + task_name="pangenome-update", + ) + + def build_longreadsv_job(self, out_vcf: pathlib.Path) -> Job: + """Call SVs from the long reads for the graph update""" + assert self.model_bundle is not None + driver = Driver( + reference=self.reference, + thread_count=self.cores, + input=list(self.lr_aln), + ) + driver.add_algo( + LongReadSV( + out_vcf, + model=self.model_bundle.joinpath("longreadsv.model"), + min_sv_size=LONGREADSV_MIN_SV_SIZE, + ) + ) + return Job( + Pipeline(Command(*driver.build_cmd())), + "longreadsv", + self.cores, + task_name="pangenome-update", + ) + + def build_minimap2_lift_job( + self, + out_bam: pathlib.Path, + ext_fastq: pathlib.Path, + pangenome_fasta: pathlib.Path, + pangenome_gfa: pathlib.Path, + ) -> Job: + """Align the extracted reads to the personalized pangenome and lift + the alignments back to the linear reference""" + assert self.reference is not None + assert self.model_bundle is not None + + rg = copy.deepcopy(self.fastq_readgroup) + rg["ID"] = rg["ID"] + "-pg" + rg["SM"] = self.sample_sm + rg["LR"] = "2" + mm2_job = Job( + cmds.cmd_minimap2_lift( + out_bam, + pangenome_fasta, + ext_fastq, + pangenome_gfa, + self.reference, + "@RG\\t" + "\\t".join([f"{x[0]}:{x[1]}" for x in rg.items()]), + self.model_bundle.joinpath("minimap2.model"), + threads=self.cores, + mm2_xargs=["--secondary=yes"], + ), + "mm2-lift", + self.cores, + task_name="pangenome-alignment", + ) + return mm2_job + + def build_replace_rg(self) -> List[List[str]]: + """`--replace_rg` arguments setting the LR readgroup attribute for + the variant calling stages. + + The bwa (LR:0) and lifted (LR:2) alignments are generated by the + pipeline with the LR attribute already in place; only the input + long-read readgroups are rewritten (LR:1), as they cannot be + assumed to carry the attribute. + """ + replace_rg: List[List[str]] = [[], []] # bwa and lifted alignments + for aln_rgs in self.lr_readgroups: + replace_rg.append([]) + for aln_rg in aln_rgs: + rg_id = aln_rg.get("ID") + sm = self.rgsm if self.rgsm else aln_rg.get("SM") + replace_rg[-1].append(f"{rg_id}=ID:{rg_id}\\tSM:{sm}\\tLR:1") + return replace_rg + + def build_pangenomesv_job( + self, + out_vcf: pathlib.Path, + input_bams: List[pathlib.Path], + pangenome_gfa: pathlib.Path, + replace_rg: List[List[str]], + ) -> Job: + """Call SVs with the original, lifted, and long reads""" + driver = Driver( + reference=self.reference, + thread_count=self.cores, + input=input_bams, + replace_rg=replace_rg, + ) + driver.add_algo( + PangenomeSV( + out_vcf, + gfa_file=pangenome_gfa, + min_af=PANGENOME_SV_MIN_AF, + ) + ) + return Job( + Pipeline(Command(*driver.build_cmd())), + "pangenome-sv", + self.cores, + task_name="sv-calling", + ) + + def build_dnascope_job( + self, + out_vcf: pathlib.Path, + input_bams: List[pathlib.Path], + replace_rg: List[List[str]], + ) -> Job: + """Call small variants with the original, lifted, and long reads""" + assert self.model_bundle is not None + pcr_indel_model = "NONE" if self.pcr_free else "CONSERVATIVE" + driver = Driver( + reference=self.reference, + thread_count=self.cores, + input=input_bams, + interval=self.bed, + replace_rg=replace_rg, + ) + driver.add_algo( + DNAscope( + out_vcf, + model=self.model_bundle.joinpath("dnascope.model"), + pcr_indel_model=pcr_indel_model, + dbsnp=self.dbsnp, + ) + ) + return Job( + Pipeline(Command(*driver.build_cmd())), + "dnascope-raw", + self.cores, + task_name="variant-calling", + ) + + def build_dnamodelapply_job( + self, + in_vcf: pathlib.Path, + out_vcf: pathlib.Path, + ) -> Job: + """Apply the DNAscope model""" + assert self.model_bundle is not None + driver = Driver( + reference=self.reference, + thread_count=self.cores, + ) + driver.add_algo( + DNAModelApply( + model=self.model_bundle.joinpath("dnascope.model"), + vcf=in_vcf, + output=out_vcf, + ) + ) + return Job( + Pipeline(Command(*driver.build_cmd())), + "model-apply", + self.cores, + task_name="model-apply", + ) diff --git a/sentieon_cli/sentieon_pangenome.py b/sentieon_cli/sentieon_pangenome.py index bb985db..c6816ae 100644 --- a/sentieon_cli/sentieon_pangenome.py +++ b/sentieon_cli/sentieon_pangenome.py @@ -20,24 +20,14 @@ from .base_pangenome import BasePangenome, SampleSex from .dag import DAG from .driver import ( - AlignmentStat, - BaseDistributionByCycle, CNVModelApply, CNVscope, - CoverageMetrics, - Dedup, DNAModelApply, DNAscope, Driver, - GCBias, GVCFtyper, - InsertSizeMetricAlgo, - LocusCollector, - MeanQualityByCycle, PangenomeSV, - QualDistribution, ReadWriter, - WgsMetricsAlgo, ) from .job import Job from .logging import get_logger @@ -774,7 +764,10 @@ def build_first_dag(self) -> DAG: out_bwa_aln, [bwa_bam], "bwa", metrics=dedup_metrics ) mm2_lc_job, mm2_dedup_job = self.build_dedup_job( - out_mm2_aln, [mm2_bam], "mm2", left_align=True + out_mm2_aln, + [mm2_bam], + "mm2", + left_align_rgid=f"{self.fastq_readgroup['ID']}-mm2", ) dag.add_job(bwa_lc_job, bwa_lc_dependencies) dag.add_job(bwa_dedup_job, {bwa_lc_job}) @@ -1084,135 +1077,6 @@ def build_minimap2_lift_job( ) return mm2_job - def build_dedup_job( - self, - output_bam, - input_bam: List[pathlib.Path], - tag: str, - left_align=False, - metrics: Optional[pathlib.Path] = None, - ) -> Tuple[Job, Job]: - """Build deduplication job""" - score_file = self.tmp_dir.joinpath(f"sample-{tag}-score.txt.gz") - - read_filters = [] - if left_align: - read_filters.append( - "IndelLeftAlignReadTransform," - f"rgid={self.fastq_readgroup['ID']}-mm2" - ) - - # LocusCollector + Dedup - driver = Driver( - reference=self.reference, - thread_count=self.cores, - input=input_bam, - read_filter=read_filters, - ) - driver.add_algo(LocusCollector(score_file)) - - lc_job = Job( - Pipeline(Command(*driver.build_cmd())), - f"locuscollector-{tag}", - self.cores, - task_name="dedup", - ) - - driver2 = Driver( - reference=self.reference, - thread_count=self.cores, - input=input_bam, - read_filter=read_filters, - ) - driver2.add_algo(Dedup(output_bam, score_file, metrics=metrics)) - - dedup_job = Job( - Pipeline(Command(*driver2.build_cmd())), - f"dedup-{tag}", - self.cores, - task_name="dedup", - ) - - return lc_job, dedup_job - - def build_metrics_job( - self, - sample_input: List[pathlib.Path], - ) -> Tuple[Job, Job]: - """Build a metrics job""" - if not self.output_vcf: - self.logger.error("output_vcf is required") - sys.exit(2) - - # Create the metrics directory - sample_name = self.output_vcf.name.replace(".vcf.gz", "") - metric_base = sample_name + ".txt" - metrics_dir = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_metrics") - ) - if not self.dry_run: - metrics_dir.mkdir(exist_ok=True) - - is_metrics = metrics_dir.joinpath(metric_base + ".insert_size.txt") - mqbc_metrics = metrics_dir.joinpath( - metric_base + ".mean_qual_by_cycle.txt" - ) - bdbc_metrics = metrics_dir.joinpath( - metric_base + ".base_distribution_by_cycle.txt" - ) - qualdist_metrics = metrics_dir.joinpath( - metric_base + ".qual_distribution.txt" - ) - as_metrics = metrics_dir.joinpath(metric_base + ".alignment_stat.txt") - coverage_metrics = metrics_dir.joinpath("coverage") - - # WGS metrics - wgs_metrics = metrics_dir.joinpath(metric_base + ".wgs.txt") - gc_metrics = metrics_dir.joinpath(metric_base + ".gc_bias.txt") - gc_summary = metrics_dir.joinpath(metric_base + ".gc_bias_summary.txt") - - driver = Driver( - reference=self.reference, - thread_count=self.cores, - input=sample_input, - ) - - driver.add_algo(InsertSizeMetricAlgo(is_metrics)) - driver.add_algo(MeanQualityByCycle(mqbc_metrics)) - driver.add_algo(BaseDistributionByCycle(bdbc_metrics)) - driver.add_algo(QualDistribution(qualdist_metrics)) - driver.add_algo(AlignmentStat(as_metrics)) - driver.add_algo(GCBias(gc_metrics, summary=gc_summary)) - driver.add_algo(WgsMetricsAlgo(wgs_metrics, include_unpaired="true")) - driver.add_algo(CoverageMetrics(coverage_metrics)) - - metrics_job = Job( - Pipeline(Command(*driver.build_cmd())), - "metrics", - 0, - task_name="metrics", - ) - - rehead_script = pathlib.Path( - str( - files("sentieon_cli.scripts").joinpath("rehead_wgs_metrics.py") - ) - ) - rehead_job = Job( - Pipeline( - Command( - sys.executable, - str(rehead_script), - "--metrics_file", - str(wgs_metrics), - ) - ), - "Rehead metrics", - 0, - task_name="metrics", - ) - return (metrics_job, rehead_job) - def build_dnascope_job( self, out_vcf: pathlib.Path, diff --git a/tests/unit/test_hybrid_pangenome.py b/tests/unit/test_hybrid_pangenome.py new file mode 100644 index 0000000..d2b99c7 --- /dev/null +++ b/tests/unit/test_hybrid_pangenome.py @@ -0,0 +1,474 @@ +""" +Unit tests for the HybridPangenome pipeline logic +""" + +import os +import pathlib +import sys +import tempfile +from unittest.mock import MagicMock + +# Add the parent directory to the path to import sentieon_cli +sys.path.insert( + 0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) +) + +from sentieon_cli.hybrid_pangenome import HybridPangenome +from sentieon_cli.command_strings import LONGREAD_SV_BED_AWK +from sentieon_cli.dag import DAG + + +class TestHybridPangenome: + """Test the hybrid-pangenome pipeline logic""" + + def setup_method(self): + """Setup test fixtures""" + self.temp_dir = tempfile.mkdtemp() + self.mock_dir = pathlib.Path(self.temp_dir) + + self.mock_vcf = self.mock_dir / "output.vcf.gz" + self.mock_ref = self.mock_dir / "reference.fa" + self.mock_lr_bam = self.mock_dir / "longreads.bam" + self.mock_bundle = self.mock_dir / "model.bundle" + self.mock_gbz = self.mock_dir / "pangenome.grch38.gbz" + self.mock_hapl = self.mock_dir / "pangenome.grch38.hapl" + self.mock_pop_vcf = self.mock_dir / "population.vcf.gz" + self.mock_dbsnp = self.mock_dir / "dbsnp.vcf.gz" + self.mock_bed = self.mock_dir / "autosomes.bed" + self.mock_r1 = self.mock_dir / "sample_R1.fastq.gz" + self.mock_r2 = self.mock_dir / "sample_R2.fastq.gz" + + for file_path in [ + self.mock_ref, + self.mock_lr_bam, + self.mock_bundle, + self.mock_gbz, + self.mock_hapl, + self.mock_pop_vcf, + self.mock_dbsnp, + self.mock_bed, + self.mock_r1, + self.mock_r2, + ]: + file_path.touch() + + with open(str(self.mock_ref) + ".fai", "w") as f: + f.write("chr1\t1000\t0\t80\t81\n") + + def create_pipeline(self): + """Create a HybridPangenome pipeline for testing""" + pipeline = HybridPangenome() + + pipeline.logger = MagicMock() + + # Configure arguments + pipeline.output_vcf = self.mock_vcf + pipeline.reference = self.mock_ref + pipeline.model_bundle = self.mock_bundle + pipeline.r1_fastq = [self.mock_r1] + pipeline.r2_fastq = [self.mock_r2] + pipeline.lr_aln = [self.mock_lr_bam] + pipeline.gbz = self.mock_gbz + pipeline.hapl = self.mock_hapl + pipeline.pop_vcf = self.mock_pop_vcf + pipeline.dbsnp = self.mock_dbsnp + pipeline.bed = self.mock_bed + pipeline.cores = 2 + pipeline.dry_run = True + pipeline.skip_version_check = True + pipeline.skip_multiqc = True + pipeline.tmp_dir = self.mock_dir + + # State normally set by validate() + pipeline.fai_data = {"chr1": {"length": 1000}} + pipeline.shards = [MagicMock()] + pipeline.shards[0].contig = "chr1" + pipeline.shards[0].start = 1 + pipeline.shards[0].stop = 1000 + pipeline.pop_vcf_contigs = {"chr1": 1000} + pipeline.fastq_readgroup = {"ID": "rg1", "SM": "sample1"} + pipeline.lr_readgroups = [[{"ID": "lr-rg1", "SM": "sample1"}]] + pipeline.sample_sm = "sample1" + + return pipeline + + 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()) + return [job.name for job in all_jobs], all_jobs + + def _get_job(self, all_jobs, name): + return next(j for j in all_jobs if j.name == name) + + def test_dag_jobs_present(self): + """All expected jobs are in the DAG""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + assert isinstance(dag, DAG) + + job_names, _ = self._get_all_job_names(dag) + for name in ( + "kmc", + "bwa-extract", + "vg-haplotypes", + "vg-convert-gfa", + "graph-update-raw", + "longreadsv", + "longread-sv-bed", + "graph-update", + "gfa2fa", + "faidx", + "mm2-lift", + "locuscollector-bwa", + "dedup-bwa", + "locuscollector-lift", + "dedup-lift", + "metrics", + "pangenome-sv", + "dnascope-raw", + "model-apply", + ): + assert name in job_names, f"missing job: {name}" + + def test_kmc_command(self): + """KMC counts the short-read fastq and long-read fasta together""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "kmc").shell) + assert "-k29" in cmd_str + assert "-m30" in cmd_str + assert "-okff" in cmd_str + assert "-fa /dev/stdin" in cmd_str + assert "samtools fasta" in cmd_str + assert str(self.mock_lr_bam) in cmd_str + assert str(self.mock_r1) in cmd_str + assert str(self.mock_r2) in cmd_str + + def test_bwa_readgroup_lr0(self): + """The bwa alignment carries the LR:0 readgroup attribute""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "bwa-extract").shell) + # The user's RGID is used unchanged + assert r"ID:rg1\t" in cmd_str + assert "ID:rg1-bwa" not in cmd_str + assert "LR:0" in cmd_str + assert "pgutil extract" in cmd_str + assert "bwa.model" in cmd_str + # The alignment is written to the temporary directory for dedup + bwa_bam = self.mock_dir / "sample-bwa.bam" + assert f"-b {bwa_bam}" in cmd_str + + def test_mm2_lift_readgroup_lr2(self): + """The lifted alignment carries the LR:2 readgroup attribute""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "mm2-lift").shell) + assert "ID:rg1-pg" in cmd_str + assert "LR:2" in cmd_str + assert "--secondary=yes" in cmd_str + assert "pgutil lift" in cmd_str + assert "--prefix" in cmd_str + assert "GRCh38#0#" in cmd_str + assert "minimap2.model" in cmd_str + # `util sort` writes the sorted alignment to the temporary directory + assert "util sort" in cmd_str + lift_bam = self.mock_dir / "sample-lift.bam" + assert f"-o {lift_bam}" in cmd_str + + def test_dedup_commands(self): + """The short-read alignments are deduplicated, with Dedup metrics + from the bwa alignment only""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + bwa_bam = self.mock_dir / "sample-bwa.bam" + lift_bam = self.mock_dir / "sample-lift.bam" + out_bwa = str(self.mock_vcf).replace(".vcf.gz", "_bwa_deduped.cram") + out_lift = str(self.mock_vcf).replace(".vcf.gz", "_lift_deduped.cram") + dedup_metrics = str( + self.mock_dir / "output_metrics" / "output.txt.dedup_metrics.txt" + ) + + lc_cmd = str(self._get_job(all_jobs, "locuscollector-bwa").shell) + assert "--algo LocusCollector" in lc_cmd + assert str(bwa_bam) in lc_cmd + + bwa_cmd = str(self._get_job(all_jobs, "dedup-bwa").shell) + assert "--algo Dedup" in bwa_cmd + assert str(bwa_bam) in bwa_cmd + assert out_bwa in bwa_cmd + assert f"--metrics {dedup_metrics}" in bwa_cmd + assert "IndelLeftAlignReadTransform" not in bwa_cmd + + lift_cmd = str(self._get_job(all_jobs, "dedup-lift").shell) + assert "--algo Dedup" in lift_cmd + assert str(lift_bam) in lift_cmd + assert out_lift in lift_cmd + assert "IndelLeftAlignReadTransform" not in lift_cmd + # Dedup metrics are only collected from the bwa alignment + assert "--metrics" not in lift_cmd + + # The long-read input is assumed to be deduplicated already + for job_name in ("locuscollector-bwa", "dedup-bwa", "dedup-lift"): + assert str(self.mock_lr_bam) not in str( + self._get_job(all_jobs, job_name).shell + ) + + def test_metrics_job(self): + """Metrics are collected from the deduplicated alignments""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + job_names, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "metrics").shell) + assert ( + str(self.mock_vcf).replace(".vcf.gz", "_bwa_deduped.cram") + in cmd_str + ) + assert ( + str(self.mock_vcf).replace(".vcf.gz", "_lift_deduped.cram") + in cmd_str + ) + assert "--algo GCBias" in cmd_str + assert "--algo WgsMetricsAlgo" in cmd_str + assert "Rehead metrics" in job_names + + def test_skip_metrics(self): + """skip_metrics removes the metrics jobs and the Dedup metrics""" + pipeline = self.create_pipeline() + pipeline.skip_metrics = True + dag = pipeline.build_dag() + job_names, all_jobs = self._get_all_job_names(dag) + + assert "metrics" not in job_names + assert "Rehead metrics" not in job_names + assert "multiqc" not in job_names + # Deduplication still runs + assert "dedup-bwa" in job_names + assert "dedup-lift" in job_names + assert "--metrics" not in str( + self._get_job(all_jobs, "dedup-bwa").shell + ) + + def test_graph_update_commands(self): + """PGHapUpdateAlgo runs without and then with the SV BED""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + raw_cmd = str(self._get_job(all_jobs, "graph-update-raw").shell) + assert "--algo PGHapUpdateAlgo" in raw_cmd + assert "--gfa_file " in raw_cmd + assert "sample-hap.raw.gfa" in raw_cmd + assert "sample-pangenome-raw.gfa" in raw_cmd + assert "--target_bed" not in raw_cmd + assert str(self.mock_lr_bam) in raw_cmd + + update_cmd = str(self._get_job(all_jobs, "graph-update").shell) + assert "--algo PGHapUpdateAlgo" in update_cmd + assert "--gfa_file " in update_cmd + assert "sample-pangenome-raw.gfa" in update_cmd + assert "--target_bed " in update_cmd + assert "sample-sv.bed" in update_cmd + + def test_longreadsv_and_bed(self): + """LongReadSV runs on the long reads and the awk BED script is + retained verbatim""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + sv_cmd = str(self._get_job(all_jobs, "longreadsv").shell) + assert "--algo LongReadSV" in sv_cmd + assert "longreadsv.model" in sv_cmd + assert "--min_sv_size 20" in sv_cmd + assert str(self.mock_lr_bam) in sv_cmd + + bed_cmd = str(self._get_job(all_jobs, "longread-sv-bed").shell) + assert 'n=split($10,a,":")' in bed_cmd + assert "sort -k1,1 -k2,2n" in bed_cmd + assert "bedtools merge" in bed_cmd + assert "sample-sv.bed" in bed_cmd + # The awk script matches the validated implementation + assert "chr[^_]+_[0-9]+_[0-9]+" in LONGREAD_SV_BED_AWK + + def test_gfa2fa_and_faidx(self): + """The updated graph is converted to an indexed FASTA""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + gfa2fa_cmd = str(self._get_job(all_jobs, "gfa2fa").shell) + assert "pgutil gfa2fa" in gfa2fa_cmd + assert str(self.mock_ref) + ".fai" in gfa2fa_cmd + assert "sample-pangenome.gfa" in gfa2fa_cmd + assert "sample-pangenome.fa" in gfa2fa_cmd + + faidx_cmd = str(self._get_job(all_jobs, "faidx").shell) + assert "samtools faidx" in faidx_cmd + assert "sample-pangenome.fa" in faidx_cmd + + def test_calling_inputs_and_replace_rg(self): + """PangenomeSV and DNAscope use the bwa, lifted, and long reads, + rewriting only the long-read readgroups with LR:1""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + bwa_aln = str(self.mock_vcf).replace(".vcf.gz", "_bwa_deduped.cram") + lift_aln = str(self.mock_vcf).replace(".vcf.gz", "_lift_deduped.cram") + replace_arg = r"lr-rg1=ID:lr-rg1\tSM:sample1\tLR:1" + + for job_name in ("pangenome-sv", "dnascope-raw"): + cmd_str = str(self._get_job(all_jobs, job_name).shell) + assert bwa_aln in cmd_str + assert lift_aln in cmd_str + assert str(self.mock_lr_bam) in cmd_str + assert replace_arg in cmd_str + # The long-read input is preceded by its --replace_rg argument + assert cmd_str.index(replace_arg) < cmd_str.index( + str(self.mock_lr_bam) + ) + # The pipeline-generated alignments are not rewritten + assert cmd_str.count("--replace_rg") == 1 + + def test_pangenomesv_command(self): + """PangenomeSV uses the updated graph and min_af""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "pangenome-sv").shell) + assert "--algo PangenomeSV" in cmd_str + assert "--gfa_file" in cmd_str + assert "sample-pangenome.gfa" in cmd_str + assert "--min_af 0.1" in cmd_str + sv_vcf = str(self.mock_vcf).replace(".vcf.gz", "_sv.vcf.gz") + assert sv_vcf in cmd_str + # SV calling is not restricted to the small-variant BED + assert f"--interval {self.mock_bed}" not in cmd_str + + def test_dnascope_command(self): + """DNAscope runs with the model, interval, and pcr_indel_model""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "dnascope-raw").shell) + assert "--algo DNAscope" in cmd_str + assert "dnascope.model" in cmd_str + assert "--pcr_indel_model CONSERVATIVE" in cmd_str + assert f"--interval {self.mock_bed}" in cmd_str + assert f"--dbsnp {self.mock_dbsnp}" in cmd_str + + def test_pcr_free(self): + """--pcr_free calls DNAscope with `--pcr_indel_model NONE`""" + pipeline = self.create_pipeline() + pipeline.pcr_free = True + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "dnascope-raw").shell) + assert "--pcr_indel_model NONE" in cmd_str + + def test_bam_format(self): + """--bam_format switches the deduplicated outputs to BAM""" + pipeline = self.create_pipeline() + pipeline.bam_format = True + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + bwa_bam = str(self.mock_vcf).replace(".vcf.gz", "_bwa_deduped.bam") + lift_bam = str(self.mock_vcf).replace(".vcf.gz", "_lift_deduped.bam") + + assert bwa_bam in str(self._get_job(all_jobs, "dedup-bwa").shell) + assert lift_bam in str(self._get_job(all_jobs, "dedup-lift").shell) + + cmd_str = str(self._get_job(all_jobs, "dnascope-raw").shell) + assert bwa_bam in cmd_str + assert lift_bam in cmd_str + + def test_model_apply_writes_output_vcf(self): + """DNAModelApply produces the final output VCF""" + pipeline = self.create_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "model-apply").shell) + assert "--algo DNAModelApply" in cmd_str + assert str(self.mock_vcf) in cmd_str + + def test_skip_model_apply(self): + """With skip_model_apply the transfer writes the final VCF""" + pipeline = self.create_pipeline() + pipeline.skip_model_apply = True + dag = pipeline.build_dag() + job_names, all_jobs = self._get_all_job_names(dag) + + assert "model-apply" not in job_names + concat_job = self._get_job(all_jobs, "merge-trim-concat") + assert str(pipeline.output_vcf) in str(concat_job.shell) + + def test_skip_svs(self): + """skip_svs removes PangenomeSV but not the graph jobs""" + pipeline = self.create_pipeline() + pipeline.skip_svs = True + dag = pipeline.build_dag() + job_names, _ = self._get_all_job_names(dag) + + assert "pangenome-sv" not in job_names + assert "dnascope-raw" in job_names + assert "graph-update" in job_names + + def test_skip_small_variants(self): + """skip_small_variants removes DNAscope, transfer, and model-apply""" + pipeline = self.create_pipeline() + pipeline.skip_small_variants = True + dag = pipeline.build_dag() + job_names, _ = self._get_all_job_names(dag) + + assert "dnascope-raw" not in job_names + assert "model-apply" not in job_names + assert "merge-trim-concat" not in job_names + assert "pangenome-sv" in job_names + + def test_multiple_lr_inputs(self): + """Each long-read input gets its own --replace_rg arguments""" + lr_bam2 = self.mock_dir / "longreads2.bam" + lr_bam2.touch() + + pipeline = self.create_pipeline() + pipeline.lr_aln = [self.mock_lr_bam, lr_bam2] + pipeline.lr_readgroups = [ + [{"ID": "lr-rg1", "SM": "sample1"}], + [{"ID": "lr-rg2", "SM": "sample1"}], + ] + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "dnascope-raw").shell) + assert r"lr-rg1=ID:lr-rg1\tSM:sample1\tLR:1" in cmd_str + assert r"lr-rg2=ID:lr-rg2\tSM:sample1\tLR:1" in cmd_str + assert str(lr_bam2) in cmd_str + + def test_rgsm_overrides_sm(self): + """--rgsm overrides the SM tag in the rewritten readgroups""" + pipeline = self.create_pipeline() + pipeline.rgsm = "override_sm" + pipeline.sample_sm = "override_sm" + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "dnascope-raw").shell) + assert r"lr-rg1=ID:lr-rg1\tSM:override_sm\tLR:1" in cmd_str + + bwa_cmd = str(self._get_job(all_jobs, "bwa-extract").shell) + assert "SM:override_sm" in bwa_cmd From 74f25ba364b864b82116967f2fc414ac1a540e33 Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sat, 15 Aug 2026 19:33:16 -0700 Subject: [PATCH 02/10] Support aligned short-read and unaligned long-read input in hybrid-pangenome Aligned short-read BAM/CRAM input via `--sr_aln`, mutually exclusive with fastq input. Aligned input is assumed deduplicated: no dedup or metrics jobs run, read extraction and k-mer counting share a single pass over the input (driver ReadWriter | pgutil extract | kmc, with the long-read FASTA streams concatenated), and the lifted alignment is written directly to the output. The calling stages rewrite the input readgroups with `--replace_rg` to set LR:0. Unaligned long-read (uBAM/uCRAM) input via `--lr_align_input` and `--lr_input_ref`. Each input is realigned with minimap2 using the new `minimap2_lr.model` bundle member; k-mer counting reads the original input files in parallel with the realignment. Also fix `cmd_samtools_fastq_minimap2` and `cmd_samtools_fastq_bwa` to decode the input file with `input_ref` rather than the alignment target reference, matching the flag's documented purpose. This changes behavior for dnascope-longread and dnascope-hybrid runs where the two references differ. Co-Authored-By: Claude Fable 5 --- sentieon_cli/command_strings.py | 174 ++++++-- sentieon_cli/hybrid_pangenome.py | 491 ++++++++++++++++----- tests/unit/test_command_strings_realign.py | 138 ++++++ tests/unit/test_hybrid_pangenome.py | 320 ++++++++++++++ 4 files changed, 964 insertions(+), 159 deletions(-) create mode 100644 tests/unit/test_command_strings_realign.py diff --git a/sentieon_cli/command_strings.py b/sentieon_cli/command_strings.py index e09877a..6a091ec 100644 --- a/sentieon_cli/command_strings.py +++ b/sentieon_cli/command_strings.py @@ -10,7 +10,7 @@ import subprocess as sp import sys import typing -from typing import Any, Dict, List, Optional, Union +from typing import Any, Dict, List, Optional, Tuple, Union from .driver import BaseDriver, Driver, ReadWriter from .logging import get_logger @@ -512,12 +512,15 @@ def cmd_samtools_fastq_minimap2( fastq_taglist: str = "*", minimap2_args: str = "-YL", util_sort_args: str = "--cram_write_options version=3.0,compressor=rans", + minimap2_model: Optional[Union[pathlib.Path, str]] = None, ) -> Pipeline: """Re-align an input BAM/CRAM/uBAM/uCRAM file with minimap2""" + # `input_ref` decodes the input file, which may use a different + # reference from the alignment target ref_cmd: List[str] = [] if input_ref: - ref_cmd = ["--reference", str(reference)] + ref_cmd = ["--reference", str(input_ref)] cmd_list = [ Command( "samtools", @@ -530,6 +533,9 @@ def cmd_samtools_fastq_minimap2( str(input_aln), ) ] + mm2_model = ( + minimap2_model if minimap2_model else f"{model_bundle}/minimap2.model" + ) cmd_list.append( Command( "sentieon", @@ -540,7 +546,7 @@ def cmd_samtools_fastq_minimap2( "-a", minimap2_args, "-x", - f"{model_bundle}/minimap2.model", + str(mm2_model), str(reference), "/dev/stdin", ) @@ -597,9 +603,11 @@ def cmd_samtools_fastq_bwa( util_sort_args: str = "--cram_write_options version=3.0,compressor=rans", ) -> Pipeline: """Re-align an input BAM/CRAM/uBAM/uCRAM file with bwa""" + # `input_ref` decodes the input file, which may use a different + # reference from the alignment target ref_cmd: List[str] = [] if input_ref: - ref_cmd = ["--reference", str(reference)] + ref_cmd = ["--reference", str(input_ref)] if collate: collate_cmd = Command( @@ -879,19 +887,22 @@ def cmd_kmc( return Pipeline(Command(*cmd)) -def cmd_extract_kmc( - output_prefix: pathlib.Path, - out_fastq: pathlib.Path, +def _readwriter_extract_cmds( input_aln: List[pathlib.Path], reference: pathlib.Path, extract_model: pathlib.Path, - tmp_dir: pathlib.Path, rw_bam: pathlib.Path, + out_fastq: pathlib.Path, threads: int = 1, -) -> Pipeline: - # `rw_bam` must already be a symlink to /dev/stdout; - # ReadWriter writes to it, which resolves to the driver's stdout (= the - # pipe to pgutil extract). +) -> Tuple[Command, Command]: + """Read an aligned input once, writing the extracted reads to a fastq + file and the full read stream to stdout. + + `rw_bam` must already be a symlink to /dev/stdout; ReadWriter writes to + it, which resolves to the driver's stdout (= the pipe to pgutil + extract). Duplicate, secondary, and supplementary reads are excluded + with `--output_flag_filter 0xf00:0`. + """ driver = Driver( reference=reference, thread_count=threads, @@ -925,6 +936,27 @@ def cmd_extract_kmc( "-a", "-", ) + return driver_cmd, extract_cmd + + +def cmd_extract_kmc( + output_prefix: pathlib.Path, + out_fastq: pathlib.Path, + input_aln: List[pathlib.Path], + reference: pathlib.Path, + extract_model: pathlib.Path, + tmp_dir: pathlib.Path, + rw_bam: pathlib.Path, + threads: int = 1, +) -> Pipeline: + driver_cmd, extract_cmd = _readwriter_extract_cmds( + input_aln, + reference, + extract_model, + rw_bam, + out_fastq, + threads, + ) kmc_cmd = Command( "kmc", "-k29", @@ -1415,11 +1447,57 @@ def cmd_minimap2_lift( return Pipeline(mm2_cmd, lift_cmd, sort_cmd) +def _aln_fasta_procsubs( + aln: List[Tuple[pathlib.Path, pathlib.Path]], + threads: int = 1, +) -> List[InputProcSub]: + """FASTA process substitutions for `(alignment, decode reference)` + pairs""" + procsubs: List[InputProcSub] = [] + for aln_file, decode_ref in aln: + procsubs.append( + InputProcSub( + Pipeline( + Command( + "samtools", + "fasta", + "--reference", + str(decode_ref), + "-@", + str(threads), + str(aln_file), + ) + ) + ) + ) + return procsubs + + +def _kmc_stdin_cmd( + output_prefix: pathlib.Path, + tmp_dir: pathlib.Path, + k: int = 29, + memory: int = 30, + threads: int = 1, +) -> Command: + """The patched KMC reading a FASTA stream from stdin""" + return Command( + "kmc", + f"-k{k}", + f"-m{memory}", + "-okff", + f"-t{threads}", + "-fa", + "/dev/stdin", + str(output_prefix), + str(tmp_dir), + ) + + def cmd_hybrid_kmc( output_prefix: pathlib.Path, fastq: List[pathlib.Path], - aln: List[pathlib.Path], - reference: pathlib.Path, + aln: List[Tuple[pathlib.Path, pathlib.Path]], tmp_dir: pathlib.Path, k: int = 29, memory: int = 30, @@ -1431,7 +1509,8 @@ def cmd_hybrid_kmc( KMC accepts a single input format per run, so the fastq files are converted to FASTA and the aligned reads are extracted with `samtools fasta`, and one FASTA stream is fed to the patched KMC - through stdin. + through stdin. Each alignment is supplied as an + `(alignment, decode reference)` pair. """ cat_args: List[Union[str, InputProcSub]] = [] if fastq: @@ -1444,34 +1523,47 @@ def cmd_hybrid_kmc( Command("awk", 'NR%4==1{print ">"substr($0,2)} NR%4==2{print}'), ) cat_args.append(InputProcSub(fq_fasta)) - for aln_file in aln: - cat_args.append( - InputProcSub( - Pipeline( - Command( - "samtools", - "fasta", - "--reference", - str(reference), - "-@", - str(threads), - str(aln_file), - ) - ) - ) - ) + cat_args.extend(_aln_fasta_procsubs(aln, threads)) cat_cmd = Command("cat", *cat_args) - kmc_cmd = Command( - "kmc", - f"-k{k}", - f"-m{memory}", - "-okff", - f"-t{threads}", - "-fa", - "/dev/stdin", - str(output_prefix), - str(tmp_dir), + kmc_cmd = _kmc_stdin_cmd(output_prefix, tmp_dir, k, memory, threads) + return Pipeline(cat_cmd, kmc_cmd) + + +def cmd_hybrid_extract_kmc( + output_prefix: pathlib.Path, + out_fastq: pathlib.Path, + sr_aln: List[pathlib.Path], + lr_aln: List[Tuple[pathlib.Path, pathlib.Path]], + reference: pathlib.Path, + extract_model: pathlib.Path, + tmp_dir: pathlib.Path, + rw_bam: pathlib.Path, + k: int = 29, + memory: int = 30, + threads: int = 1, +) -> Pipeline: + """Extract reads from aligned short-read input and count k-mers across + the short and long reads in a single pass. + + The aligned short reads are read once: `pgutil extract` writes the + extracted reads to `out_fastq` and passes the full read stream on to + KMC. Each long-read alignment is supplied as an + `(alignment, decode reference)` pair. + """ + driver_cmd, extract_cmd = _readwriter_extract_cmds( + sr_aln, + reference, + extract_model, + rw_bam, + out_fastq, + threads, ) + cat_args: List[Union[str, InputProcSub]] = [ + InputProcSub(Pipeline(driver_cmd, extract_cmd)) + ] + cat_args.extend(_aln_fasta_procsubs(lr_aln, threads)) + cat_cmd = Command("cat", *cat_args) + kmc_cmd = _kmc_stdin_cmd(output_prefix, tmp_dir, k, memory, threads) return Pipeline(cat_cmd, kmc_cmd) diff --git a/sentieon_cli/hybrid_pangenome.py b/sentieon_cli/hybrid_pangenome.py index 24ef897..a5a0d62 100644 --- a/sentieon_cli/hybrid_pangenome.py +++ b/sentieon_cli/hybrid_pangenome.py @@ -11,7 +11,7 @@ import pathlib import shutil import sys -from typing import Dict, List, Optional, Set +from typing import Dict, List, Optional, Set, Tuple import packaging.version @@ -90,7 +90,10 @@ class HybridPangenome(BasePangenome): "required": True, }, "readgroup": { - "help": "Readgroup information for the fastq files.", + "help": ( + "Readgroup information for the fastq files. Required " + "with fastq input; cannot be used with `--sr_aln`." + ), }, # Additional arguments "bed": { @@ -102,6 +105,21 @@ class HybridPangenome(BasePangenome): ), "type": path_arg(exists=True, is_file=True), }, + "lr_align_input": { + "help": ( + "Align the `--lr_aln` input files to the linear " + "reference genome with minimap2. Use with unaligned " + "(uBAM or uCRAM) long-read input." + ), + "action": "store_true", + }, + "lr_input_ref": { + "help": ( + "A reference FASTA used to decode the `--lr_aln` input " + "files. Only used with `--lr_align_input`." + ), + "type": path_arg(exists=True, is_file=True), + }, "pangenome_ref_name": { "default": "GRCh38", "help": "Reference name in the pangenome (GRCh38).", @@ -120,6 +138,15 @@ class HybridPangenome(BasePangenome): "help": "Skip multiQC report generation", "action": "store_true", }, + "sr_aln": { + "nargs": "*", + "help": ( + "Aligned short-read BAM or CRAM files. Assumed " + "deduplicated or duplicate-marked. Cannot be used with " + "fastq input." + ), + "type": path_arg(exists=True, is_file=True), + }, # Hidden arguments "skip_contig_checks": { "help": argparse.SUPPRESS, @@ -156,11 +183,18 @@ def __init__(self) -> None: self.pop_vcf: Optional[pathlib.Path] = None self.readgroup: Optional[str] = None self.bed: Optional[pathlib.Path] = None + self.lr_align_input = False + self.lr_input_ref: Optional[pathlib.Path] = None self.pangenome_ref_name = "GRCh38" self.rgsm: Optional[str] = None self.extract_model_name = "extract.model" self.skip_metrics = False self.skip_multiqc = False + self.sr_aln: List[pathlib.Path] = [] + self.fastq_readgroup: Optional[Dict[str, str]] = 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 @@ -183,18 +217,8 @@ def validate(self) -> None: self.validate_bundle() self.validate_output_vcf() - if not self.r1_fastq or not self.readgroup: - self.logger.error( - "Please supply the short reads with the `--r1_fastq`, " - "`--r2_fastq`, and `--readgroup` arguments" - ) - sys.exit(2) - if len(self.r1_fastq) != len(self.r2_fastq): - self.logger.error( - "The number of input `--r1_fastq` files does not equal the " - "number of `--r2_fastq` files" - ) - sys.exit(2) + self.validate_sr_inputs() + if not self.lr_aln: self.logger.error( "Please supply the long-read alignments with the `--lr_aln` " @@ -202,7 +226,20 @@ def validate(self) -> None: ) sys.exit(2) - self.validate_bwa_index() + if self.lr_input_ref and not self.lr_align_input: + self.logger.warning( + "The `--lr_input_ref` argument is only used with " + "`--lr_align_input` and will be ignored" + ) + if self.lr_align_input and not self.lr_input_ref: + if any(str(aln).endswith(".cram") for aln in self.lr_aln): + self.logger.warning( + "CRAM input may fail to decode without the " + "`--lr_input_ref` argument" + ) + + if self.r1_fastq: + self.validate_bwa_index() self.collect_readgroups() self.validate_readgroups() @@ -272,6 +309,45 @@ def validate(self) -> None: ) sys.exit(2) + def validate_sr_inputs(self) -> None: + """Validate the short-read input arguments. + + Short reads are supplied either as fastq (with `--readgroup`) or as + aligned BAM/CRAM files (`--sr_aln`), but not both. + """ + if self.r1_fastq and self.sr_aln: + self.logger.error( + "Supplying both fastq (`--r1_fastq`) and aligned " + "(`--sr_aln`) short reads is not supported" + ) + sys.exit(2) + if not self.r1_fastq and not self.sr_aln: + self.logger.error( + "Please supply the short reads with the `--sr_aln` argument " + "or with the `--r1_fastq`, `--r2_fastq`, and `--readgroup` " + "arguments" + ) + sys.exit(2) + + if self.r1_fastq: + if not self.readgroup: + self.logger.error( + "The `--readgroup` argument is required with fastq input" + ) + sys.exit(2) + if len(self.r1_fastq) != len(self.r2_fastq): + self.logger.error( + "The number of input `--r1_fastq` files does not equal " + "the number of `--r2_fastq` files" + ) + sys.exit(2) + elif self.readgroup: + self.logger.error( + "The `--readgroup` argument cannot be used with aligned " + "short-read input (`--sr_aln`)" + ) + sys.exit(2) + def validate_bundle(self) -> None: """Validate the model bundle""" bundle_info_bytes = ar_load( @@ -310,12 +386,15 @@ def validate_bundle(self) -> None: self.extract_model_name = extract_candidate required_members = { - "bwa.model", "dnascope.model", "longreadsv.model", "minimap2.model", self.extract_model_name, } + if self.r1_fastq: + required_members.add("bwa.model") + if self.lr_align_input: + required_members.add("minimap2_lr.model") missing_members = required_members - bundle_members if missing_members: self.logger.error( @@ -340,29 +419,37 @@ def validate_bundle(self) -> None: def collect_readgroups(self) -> None: """Collect readgroup tags from the inputs""" - assert self.readgroup is not None - try: - parsed_rg = parse_rg_line(self.readgroup.replace(r"\t", "\t")) - except ValueError as e: - self.logger.error( - "Invalid --readgroup value '%s': %s", self.readgroup, e - ) - sys.exit(2) - if not parsed_rg.get("ID"): - self.logger.error( - "Readgroup '%s' does not have a RGID tag", - self.readgroup, - ) - sys.exit(2) - if not parsed_rg.get("SM"): - self.logger.error( - "Readgroup '%s' does not have a RGSM tag", - self.readgroup, - ) - sys.exit(2) - self.fastq_readgroup: Dict[str, str] = parsed_rg + if self.readgroup: + try: + parsed_rg = parse_rg_line(self.readgroup.replace(r"\t", "\t")) + except ValueError as e: + self.logger.error( + "Invalid --readgroup value '%s': %s", self.readgroup, e + ) + sys.exit(2) + if not parsed_rg.get("ID"): + self.logger.error( + "Readgroup '%s' does not have a RGID tag", + self.readgroup, + ) + sys.exit(2) + if not parsed_rg.get("SM"): + self.logger.error( + "Readgroup '%s' does not have a RGSM tag", + self.readgroup, + ) + sys.exit(2) + self.fastq_readgroup = parsed_rg - self.lr_readgroups: List[List[Dict[str, str]]] = [] + # Read the readgroups from the ORIGINAL long-read input; minimap2 + # realignment preserves the readgroup IDs with `addreplacerg` + self.sr_readgroups = [] + for aln in self.sr_aln: + self.sr_readgroups.append([]) + for rg_line in cmds.get_rg_lines(aln, self.dry_run): + self.sr_readgroups[-1].append(parse_rg_line(rg_line)) + + self.lr_readgroups = [] for aln in self.lr_aln: self.lr_readgroups.append([]) for rg_line in cmds.get_rg_lines(aln, self.dry_run): @@ -370,8 +457,15 @@ def collect_readgroups(self) -> None: def validate_readgroups(self) -> None: """Confirm that all readgroups have a consistent SM tag""" - rg_sm = self.fastq_readgroup.get("SM") - for aln, aln_rgs in zip(self.lr_aln, self.lr_readgroups): + rg_sm: Optional[str] = None + if self.fastq_readgroup: + rg_sm = self.fastq_readgroup.get("SM") + elif self.sr_readgroups and self.sr_readgroups[0]: + rg_sm = self.sr_readgroups[0][0].get("SM") + + aln_inputs = list(self.sr_aln) + list(self.lr_aln) + aln_readgroups = self.sr_readgroups + self.lr_readgroups + for aln, aln_rgs in zip(aln_inputs, aln_readgroups): for aln_rg in aln_rgs: if not aln_rg.get("ID"): self.logger.error( @@ -393,7 +487,7 @@ def validate_readgroups(self) -> None: if sm != rg_sm: self.logger.error( "Input readgroup '%s' has a different RG-SM tag " - "from the `--readgroup` argument.\n" + "from the other sample inputs.\n" "found='%s' expected='%s'. Please set the `--rgsm` " "argument to override the SM tag in the input files", str(aln_rg), @@ -401,7 +495,7 @@ def validate_readgroups(self) -> None: rg_sm, ) sys.exit(2) - self.sample_sm: str = self.rgsm if self.rgsm else str(rg_sm) + self.sample_sm = self.rgsm if self.rgsm else str(rg_sm) def configure(self) -> None: """Configure pipeline parameters""" @@ -466,19 +560,49 @@ def build_dag(self) -> DAG: raw_vcf = self.tmp_dir.joinpath("sample-dnascope.vcf.gz") transfer_vcf = self.tmp_dir.joinpath("sample-dnascope_transfer.vcf.gz") - total_mem_gb = total_memory() / (1024.0**3) - - # KMC k-mer counting across the short and long reads - kmc_job = self.build_hybrid_kmc_job(kmer_prefix) - dag.add_job(kmc_job) - haplotype_dependencies: Set[Job] = {kmc_job} + # Realign unaligned (uBAM or uCRAM) long-read input. K-mer counting + # reads the original input files, so it runs without waiting for + # the realignment. + calling_lr = list(self.lr_aln) + realign_jobs: Set[Job] = set() + if self.lr_align_input: + calling_lr, realign_jobs = self.lr_align_inputs() + for job in realign_jobs: + dag.add_job(job) + + haplotype_dependencies: Set[Job] = set() + mm2_dependencies: Set[Job] = set() + if self.r1_fastq: + # KMC k-mer counting across the short and long reads + kmc_job = self.build_hybrid_kmc_job(kmer_prefix) + dag.add_job(kmc_job) + haplotype_dependencies.add(kmc_job) + + # BWA alignment and extraction + bwa_job = self.build_alignment_job(bwa_bam, ext_fastq) + dag.add_job(bwa_job) + mm2_dependencies.add(bwa_job) + # Do not run vg-haplotypes with bwa in low-mem environments + total_mem_gb = total_memory() / (1024.0**3) + if total_mem_gb < 70: + haplotype_dependencies.add(bwa_job) + + # The lifted alignment is deduplicated before variant calling + lift_aln = lift_bam + else: + # Read extraction and k-mer counting in a single pass over the + # aligned short-read input + symlink_job, extract_kmc_job = self.build_extract_kmc_jobs( + kmer_prefix, ext_fastq + ) + dag.add_job(symlink_job) + dag.add_job(extract_kmc_job, {symlink_job}) + haplotype_dependencies.add(extract_kmc_job) + mm2_dependencies.add(extract_kmc_job) - # BWA alignment and extraction - bwa_job = self.build_alignment_job(bwa_bam, ext_fastq) - dag.add_job(bwa_job) - # Do not run vg-haplotypes with bwa in low-mem environments - if total_mem_gb < 70: - haplotype_dependencies.add(bwa_job) + # Aligned input is not deduplicated, so the lifted alignment is + # written directly to the output file + lift_aln = out_lift_aln # vg haplotypes - create a sample-specific pangenome haplotypes_job = self.build_haplotypes_job(sample_pangenome, kmer_file) @@ -492,13 +616,14 @@ def build_dag(self) -> DAG: update_raw_job = self.build_graph_update_job( pangenome_raw_gfa, hap_raw_gfa, + calling_lr, name="graph-update-raw", ) - dag.add_job(update_raw_job, {gfa_job}) + dag.add_job(update_raw_job, {gfa_job} | realign_jobs) # Call SVs from the long reads and collect graph update regions - longreadsv_job = self.build_longreadsv_job(longread_sv_vcf) - dag.add_job(longreadsv_job) + longreadsv_job = self.build_longreadsv_job(longread_sv_vcf, calling_lr) + dag.add_job(longreadsv_job, realign_jobs) sv_bed_job = Job( cmds.cmd_longread_sv_bed(sv_bed, longread_sv_vcf), "longread-sv-bed", @@ -511,6 +636,7 @@ def build_dag(self) -> DAG: update_job = self.build_graph_update_job( pangenome_gfa, pangenome_raw_gfa, + calling_lr, bed=sv_bed, name="graph-update", ) @@ -534,58 +660,67 @@ def build_dag(self) -> DAG: # minimap2 alignment of the extracted reads with liftover mm2_job = self.build_minimap2_lift_job( - lift_bam, + lift_aln, ext_fastq, pangenome_fasta, pangenome_gfa, ) - dag.add_job(mm2_job, {bwa_job, faidx_job, update_job}) + dag.add_job(mm2_job, mm2_dependencies | {faidx_job, update_job}) + + # Variant calling with the short reads, the lifted reads, and the + # long reads + if self.r1_fastq: + # Deduplicate the short-read alignments. The `--lr_aln` input + # is assumed to be deduplicated already. + # Emit Dedup metrics for the primary (bwa) short-read alignment + # so they land in the metrics directory scanned by MultiQC. + dedup_metrics: Optional[pathlib.Path] = None + if not self.skip_metrics: + metrics_dir = pathlib.Path( + str(self.output_vcf).replace(".vcf.gz", "_metrics") + ) + if not self.dry_run: + metrics_dir.mkdir(exist_ok=True) + sample_name = self.output_vcf.name.replace(".vcf.gz", "") + dedup_metrics = metrics_dir.joinpath( + sample_name + ".txt.dedup_metrics.txt" + ) - # Deduplicate the short-read alignments. The `--lr_aln` input is - # assumed to be deduplicated already. - # Emit Dedup metrics for the primary (bwa) short-read alignment so - # they land in the metrics directory scanned by MultiQC. - dedup_metrics: Optional[pathlib.Path] = None - if not self.skip_metrics: - metrics_dir = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_metrics") - ) - if not self.dry_run: - metrics_dir.mkdir(exist_ok=True) - sample_name = self.output_vcf.name.replace(".vcf.gz", "") - dedup_metrics = metrics_dir.joinpath( - sample_name + ".txt.dedup_metrics.txt" + bwa_lc_job, bwa_dedup_job = self.build_dedup_job( + out_bwa_aln, [bwa_bam], "bwa", metrics=dedup_metrics ) - - bwa_lc_job, bwa_dedup_job = self.build_dedup_job( - out_bwa_aln, [bwa_bam], "bwa", metrics=dedup_metrics - ) - lift_lc_job, lift_dedup_job = self.build_dedup_job( - out_lift_aln, - [lift_bam], - "lift", - ) - dag.add_job(bwa_lc_job, {bwa_job}) - dag.add_job(bwa_dedup_job, {bwa_lc_job}) - dag.add_job(lift_lc_job, {mm2_job}) - dag.add_job(lift_dedup_job, {lift_lc_job}) - - # Alignment metrics from the deduplicated short reads - if not self.skip_metrics: - metrics_job, rehead_job = self.build_metrics_job( - [out_bwa_aln, out_lift_aln], + lift_lc_job, lift_dedup_job = self.build_dedup_job( + out_lift_aln, + [lift_bam], + "lift", ) - dag.add_job(metrics_job, {bwa_dedup_job, lift_dedup_job}) - dag.add_job(rehead_job, {metrics_job}) - if not self.skip_multiqc: - multiqc_job = self.multiqc() - if multiqc_job: - dag.add_job(multiqc_job, {rehead_job}) - - # Variant calling with the original, lifted, and long reads - calling_bams = [out_bwa_aln, out_lift_aln] + list(self.lr_aln) + dag.add_job(bwa_lc_job, {bwa_job}) + dag.add_job(bwa_dedup_job, {bwa_lc_job}) + dag.add_job(lift_lc_job, {mm2_job}) + dag.add_job(lift_dedup_job, {lift_lc_job}) + + # Alignment metrics from the deduplicated short reads + if not self.skip_metrics: + metrics_job, rehead_job = self.build_metrics_job( + [out_bwa_aln, out_lift_aln], + ) + dag.add_job(metrics_job, {bwa_dedup_job, lift_dedup_job}) + dag.add_job(rehead_job, {metrics_job}) + if not self.skip_multiqc: + multiqc_job = self.multiqc() + if multiqc_job: + dag.add_job(multiqc_job, {rehead_job}) + + calling_bams = [out_bwa_aln, out_lift_aln] + calling_lr + calling_dependencies = {bwa_dedup_job, lift_dedup_job} + else: + # Aligned short reads are used as-is; duplicate, secondary, and + # supplementary reads are excluded from the lifted alignment + # during read extraction + calling_bams = list(self.sr_aln) + [out_lift_aln] + calling_lr + calling_dependencies = {mm2_job} + calling_dependencies |= realign_jobs replace_rg = self.build_replace_rg() - calling_dependencies = {bwa_dedup_job, lift_dedup_job} if not self.skip_svs: pangenomesv_job = self.build_pangenomesv_job( @@ -633,16 +768,27 @@ def build_dag(self) -> DAG: return dag + def lr_kmc_pairs(self) -> List[Tuple[pathlib.Path, pathlib.Path]]: + """`(alignment, decode reference)` pairs for the long-read k-mer + counting. + + K-mer counting reads the ORIGINAL `--lr_aln` input, which may use a + different reference from the alignment target. + """ + assert self.reference is not None + decode_ref = self.reference + if self.lr_align_input and self.lr_input_ref: + decode_ref = self.lr_input_ref + return [(aln, decode_ref) for aln in self.lr_aln] + def build_hybrid_kmc_job(self, kmer_prefix: pathlib.Path) -> Job: """Count k-mers across the short-read fastq and long-read alignments""" - assert self.reference is not None kmc_job = Job( cmds.cmd_hybrid_kmc( kmer_prefix, self.r1_fastq + self.r2_fastq, - self.lr_aln, - self.reference, + self.lr_kmc_pairs(), self.tmp_dir, memory=self.kmer_memory, threads=self.cores, @@ -654,6 +800,85 @@ def build_hybrid_kmc_job(self, kmer_prefix: pathlib.Path) -> Job: ) return kmc_job + def build_extract_kmc_jobs( + self, + kmer_prefix: pathlib.Path, + ext_fastq: pathlib.Path, + ) -> Tuple[Job, Job]: + """Extract reads from the aligned short-read input and count k-mers + across the short and long reads in a single pass""" + assert self.reference is not None + assert self.model_bundle is not None + + # ReadWriter cannot write to /dev/stdout directly; pre-create a + # symlink so the driver writes to a real path that resolves to its + # stdout (the pipe to pgutil extract). + rw_bam = self.tmp_dir.joinpath("extract-kmc-rw.bam") + symlink_job = Job( + Pipeline(Command("ln", "-sf", "/dev/stdout", str(rw_bam))), + "extract-kmc-symlink", + 1, + task_name="read-extraction", + ) + extract_kmc_job = Job( + cmds.cmd_hybrid_extract_kmc( + kmer_prefix, + ext_fastq, + self.sr_aln, + self.lr_kmc_pairs(), + self.reference, + self.model_bundle.joinpath(self.extract_model_name), + self.tmp_dir, + rw_bam, + memory=self.kmer_memory, + threads=self.cores, + ), + "extract-kmc", + self.cores, + task_name="read-extraction", + ) + return symlink_job, extract_kmc_job + + def lr_align_inputs(self) -> Tuple[List[pathlib.Path], Set[Job]]: + """Align the long-read input files to the linear reference genome + with minimap2""" + assert self.output_vcf is not None + assert self.reference is not None + assert self.model_bundle is not None + + res: List[pathlib.Path] = [] + realign_jobs: Set[Job] = set() + suffix = "bam" if self.bam_format else "cram" + for i, input_aln in enumerate(self.lr_aln): + out_aln = pathlib.Path( + str(self.output_vcf).replace( + ".vcf.gz", f"_mm2_sorted_{i}.{suffix}" + ) + ) + rg_lines = cmds.get_rg_lines(input_aln, self.dry_run) + realign_jobs.add( + Job( + cmds.cmd_samtools_fastq_minimap2( + out_aln, + input_aln, + self.reference, + self.model_bundle, + self.cores, + rg_lines, + self.sample_sm, + self.lr_input_ref, + minimap2_model=self.model_bundle.joinpath( + "minimap2_lr.model" + ), + ), + f"bam-realign-{i}", + self.cores, + task_name="alignment", + ) + ) + res.append(out_aln) + return (res, realign_jobs) + def build_alignment_job( self, sample_bam: pathlib.Path, @@ -663,6 +888,8 @@ def build_alignment_job( assert self.reference is not None assert self.model_bundle is not None + assert self.fastq_readgroup is not None + rg = copy.deepcopy(self.fastq_readgroup) rg["SM"] = self.sample_sm rg["LR"] = "0" @@ -733,6 +960,7 @@ def build_graph_update_job( self, out_gfa: pathlib.Path, in_gfa: pathlib.Path, + lr_aln: List[pathlib.Path], bed: Optional[pathlib.Path] = None, name: str = "graph-update", ) -> Job: @@ -740,7 +968,7 @@ def build_graph_update_job( driver = Driver( reference=self.reference, thread_count=self.cores, - input=list(self.lr_aln), + input=list(lr_aln), ) driver.add_algo( PGHapUpdateAlgo(out_gfa, gfa_file=in_gfa, target_bed=bed) @@ -752,13 +980,17 @@ def build_graph_update_job( task_name="pangenome-update", ) - def build_longreadsv_job(self, out_vcf: pathlib.Path) -> Job: + def build_longreadsv_job( + self, + out_vcf: pathlib.Path, + lr_aln: List[pathlib.Path], + ) -> Job: """Call SVs from the long reads for the graph update""" assert self.model_bundle is not None driver = Driver( reference=self.reference, thread_count=self.cores, - input=list(self.lr_aln), + input=list(lr_aln), ) driver.add_algo( LongReadSV( @@ -786,7 +1018,14 @@ def build_minimap2_lift_job( assert self.reference is not None assert self.model_bundle is not None - rg = copy.deepcopy(self.fastq_readgroup) + # With aligned input, the first readgroup of the first alignment + # file seeds the readgroup of the lifted alignment + rg_source = ( + self.fastq_readgroup + if self.fastq_readgroup + else self.sr_readgroups[0][0] + ) + rg = copy.deepcopy(rg_source) rg["ID"] = rg["ID"] + "-pg" rg["SM"] = self.sample_sm rg["LR"] = "2" @@ -808,22 +1047,38 @@ def build_minimap2_lift_job( ) return mm2_job + def _replace_rg_arg(self, aln_rg: Dict[str, str], lr: str) -> str: + """A `--replace_rg` argument setting the LR attribute of an input + readgroup""" + rg_id = aln_rg.get("ID") + sm = self.rgsm if self.rgsm else aln_rg.get("SM") + return f"{rg_id}=ID:{rg_id}\\tSM:{sm}\\tLR:{lr}" + def build_replace_rg(self) -> List[List[str]]: """`--replace_rg` arguments setting the LR readgroup attribute for the variant calling stages. - The bwa (LR:0) and lifted (LR:2) alignments are generated by the - pipeline with the LR attribute already in place; only the input - long-read readgroups are rewritten (LR:1), as they cannot be - assumed to carry the attribute. + The rows match the variant-calling input order: the short-read + alignments, the lifted alignment, then the long-read alignments. + Alignments generated by the pipeline (the bwa alignment with LR:0 + and the lifted alignment with LR:2) already carry the LR attribute + and take an empty row. User-supplied alignments are rewritten, as + they cannot be assumed to carry the attribute: `--sr_aln` input + with LR:0 and `--lr_aln` input with LR:1. """ - replace_rg: List[List[str]] = [[], []] # bwa and lifted alignments + replace_rg: List[List[str]] = [] + if self.r1_fastq: + replace_rg.append([]) # the bwa alignment + else: + for aln_rgs in self.sr_readgroups: + replace_rg.append( + [self._replace_rg_arg(rg, "0") for rg in aln_rgs] + ) + replace_rg.append([]) # the lifted alignment for aln_rgs in self.lr_readgroups: - replace_rg.append([]) - for aln_rg in aln_rgs: - rg_id = aln_rg.get("ID") - sm = self.rgsm if self.rgsm else aln_rg.get("SM") - replace_rg[-1].append(f"{rg_id}=ID:{rg_id}\\tSM:{sm}\\tLR:1") + replace_rg.append( + [self._replace_rg_arg(rg, "1") for rg in aln_rgs] + ) return replace_rg def build_pangenomesv_job( diff --git a/tests/unit/test_command_strings_realign.py b/tests/unit/test_command_strings_realign.py new file mode 100644 index 0000000..abfb56d --- /dev/null +++ b/tests/unit/test_command_strings_realign.py @@ -0,0 +1,138 @@ +""" +Unit tests for the BAM/CRAM re-alignment command builders +""" + +import os +import pathlib +import sys + +# Add the parent directory to the path to import sentieon_cli +sys.path.insert( + 0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) +) + +from sentieon_cli.command_strings import ( + cmd_samtools_fastq_bwa, + cmd_samtools_fastq_minimap2, +) + + +REFERENCE = pathlib.Path("/ref/hg38.fa") +INPUT_REF = pathlib.Path("/ref/hs37d5.fa") +MODEL_BUNDLE = pathlib.Path("/bundle/sample.bundle") +INPUT_ALN = pathlib.Path("/data/sample.cram") +OUT_ALN = pathlib.Path("/out/sample_mm2_sorted_0.cram") +RG_HEADER = pathlib.Path("/tmp/sample.hdr") + + +class TestSamtoolsFastqMinimap2: + """Test the minimap2 re-alignment command""" + + def build(self, **kwargs) -> str: + args = { + "out_aln": OUT_ALN, + "input_aln": INPUT_ALN, + "reference": REFERENCE, + "model_bundle": MODEL_BUNDLE, + "cores": 4, + "rg_lines": ["@RG\tID:rg1\tSM:sample1"], + "sample_name": "sample1", + } + args.update(kwargs) + return str(cmd_samtools_fastq_minimap2(**args)) + + def test_input_ref_decodes_the_input(self): + """`input_ref` decodes the input file, not the target reference""" + cmd_str = self.build(input_ref=INPUT_REF) + assert f"samtools fastq --reference {INPUT_REF}" in cmd_str + assert f"samtools fastq --reference {REFERENCE}" not in cmd_str + # The target reference is still used for alignment and sorting + assert f"--reference {REFERENCE}" in cmd_str + + def test_no_input_ref(self): + """No `--reference` is passed to samtools without `input_ref`""" + cmd_str = self.build() + assert "--reference" not in cmd_str.split("sentieon minimap2")[0] + assert "samtools fastq -@ 4" in cmd_str + + def test_default_minimap2_model(self): + """The minimap2 model defaults to the bundle's minimap2.model""" + cmd_str = self.build() + assert f"-x {MODEL_BUNDLE}/minimap2.model" in cmd_str + + def test_explicit_minimap2_model(self): + """An explicit minimap2 model overrides the default""" + model = MODEL_BUNDLE.joinpath("minimap2_lr.model") + cmd_str = self.build(minimap2_model=model) + assert f"-x {model}" in cmd_str + assert f"-x {MODEL_BUNDLE}/minimap2.model" not in cmd_str + + def test_sm_backfill(self): + """A missing SM tag is backfilled with the sample name""" + cmd_str = self.build(rg_lines=["@RG\tID:rg1"], sample_name="sample1") + assert "addreplacerg" in cmd_str + assert "SM:sample1" in cmd_str + + def test_sm_not_overwritten(self): + """An existing SM tag is left alone""" + cmd_str = self.build( + rg_lines=["@RG\tID:rg1\tSM:from_header"], + sample_name="sample1", + ) + assert "SM:from_header" in cmd_str + assert "SM:sample1" not in cmd_str + + def test_one_addreplacerg_per_readgroup(self): + """Each readgroup gets its own addreplacerg command""" + cmd_str = self.build( + rg_lines=["@RG\tID:rg1\tSM:s1", "@RG\tID:rg2\tSM:s1"], + ) + assert cmd_str.count("samtools addreplacerg") == 2 + + def test_output_and_sort(self): + """The realigned reads are sorted into the output file""" + cmd_str = self.build() + assert f"-o {OUT_ALN}" in cmd_str + assert "sentieon util sort" in cmd_str + assert "--sam2bam" in cmd_str + + +class TestSamtoolsFastqBwa: + """Test the bwa re-alignment command""" + + def build(self, **kwargs) -> str: + args = { + "out_aln": OUT_ALN, + "input_aln": INPUT_ALN, + "reference": REFERENCE, + "model_bundle": MODEL_BUNDLE, + "cores": 4, + "rg_header": RG_HEADER, + } + args.update(kwargs) + return str(cmd_samtools_fastq_bwa(**args)) + + def test_input_ref_without_collate(self): + """`input_ref` decodes the input file read by samtools fastq""" + cmd_str = self.build(input_ref=INPUT_REF) + assert f"samtools fastq --reference {INPUT_REF}" in cmd_str + assert "samtools collate" not in cmd_str + + def test_input_ref_with_collate(self): + """`input_ref` decodes the input file read by samtools collate""" + cmd_str = self.build(input_ref=INPUT_REF, collate=True) + assert f"samtools collate --reference {INPUT_REF}" in cmd_str + # The fastq stage reads the collated stream, not the input file + assert f"samtools fastq --reference {INPUT_REF}" not in cmd_str + + def test_no_input_ref(self): + """No `--reference` is passed to samtools without `input_ref`""" + cmd_str = self.build() + assert "--reference" not in cmd_str.split("sentieon bwa")[0] + assert "samtools fastq -@ 4" in cmd_str + + def test_bwa_model_and_header(self): + """bwa uses the bundle model and the readgroup header file""" + cmd_str = self.build() + assert f"-x {MODEL_BUNDLE}/bwa.model" in cmd_str + assert f"-H {RG_HEADER}" in cmd_str diff --git a/tests/unit/test_hybrid_pangenome.py b/tests/unit/test_hybrid_pangenome.py index d2b99c7..9953161 100644 --- a/tests/unit/test_hybrid_pangenome.py +++ b/tests/unit/test_hybrid_pangenome.py @@ -8,6 +8,8 @@ import tempfile from unittest.mock import MagicMock +import pytest + # Add the parent directory to the path to import sentieon_cli sys.path.insert( 0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) @@ -37,6 +39,8 @@ def setup_method(self): self.mock_bed = self.mock_dir / "autosomes.bed" self.mock_r1 = self.mock_dir / "sample_R1.fastq.gz" self.mock_r2 = self.mock_dir / "sample_R2.fastq.gz" + self.mock_sr_bam = self.mock_dir / "shortreads.bam" + self.mock_lr_ref = self.mock_dir / "lr_reference.fa" for file_path in [ self.mock_ref, @@ -49,6 +53,8 @@ def setup_method(self): self.mock_bed, self.mock_r1, self.mock_r2, + self.mock_sr_bam, + self.mock_lr_ref, ]: file_path.touch() @@ -92,6 +98,28 @@ def create_pipeline(self): return pipeline + def create_aligned_pipeline(self): + """Create a pipeline with aligned short-read input""" + pipeline = self.create_pipeline() + pipeline.r1_fastq = [] + pipeline.r2_fastq = [] + pipeline.readgroup = None + pipeline.fastq_readgroup = None + pipeline.sr_aln = [self.mock_sr_bam] + pipeline.sr_readgroups = [[{"ID": "sr-rg1", "SM": "sample1"}]] + return pipeline + + def create_lr_realign_pipeline(self, aligned_sr=False): + """Create a pipeline that realigns the long-read input""" + pipeline = ( + self.create_aligned_pipeline() + if aligned_sr + else self.create_pipeline() + ) + pipeline.lr_align_input = True + pipeline.lr_input_ref = self.mock_lr_ref + return pipeline + 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()) @@ -100,6 +128,11 @@ def _get_all_job_names(self, dag): def _get_job(self, all_jobs, name): return next(j for j in all_jobs if j.name == name) + def _get_dep_names(self, dag, all_jobs, name): + """Helper to get the dependency names of a job""" + job = self._get_job(all_jobs, name) + return {dep.name for dep in dag.waiting_jobs.get(job, set())} + def test_dag_jobs_present(self): """All expected jobs are in the DAG""" pipeline = self.create_pipeline() @@ -472,3 +505,290 @@ def test_rgsm_overrides_sm(self): bwa_cmd = str(self._get_job(all_jobs, "bwa-extract").shell) assert "SM:override_sm" in bwa_cmd + + # Aligned short-read input + + def test_aligned_dag_jobs(self): + """Aligned short-read input skips alignment, dedup, and metrics""" + pipeline = self.create_aligned_pipeline() + assert not pipeline.skip_metrics + dag = pipeline.build_dag() + job_names, _ = self._get_all_job_names(dag) + + for name in ( + "extract-kmc-symlink", + "extract-kmc", + "vg-haplotypes", + "graph-update", + "mm2-lift", + "pangenome-sv", + "dnascope-raw", + "model-apply", + ): + assert name in job_names, f"missing job: {name}" + + for name in ( + "kmc", + "bwa-extract", + "locuscollector-bwa", + "dedup-bwa", + "locuscollector-lift", + "dedup-lift", + "metrics", + "Rehead metrics", + "multiqc", + ): + assert name not in job_names, f"unexpected job: {name}" + + def test_aligned_extract_kmc_command(self): + """Read extraction and k-mer counting run in a single pass""" + pipeline = self.create_aligned_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + symlink_cmd = str(self._get_job(all_jobs, "extract-kmc-symlink").shell) + rw_bam = self.mock_dir / "extract-kmc-rw.bam" + assert f"ln -sf /dev/stdout {rw_bam}" in symlink_cmd + + cmd_str = str(self._get_job(all_jobs, "extract-kmc").shell) + # One pass over the aligned short reads, without duplicate, + # secondary, or supplementary reads + assert "--algo ReadWriter" in cmd_str + assert "--output_flag_filter 0xf00:0" in cmd_str + assert str(self.mock_sr_bam) in cmd_str + assert str(rw_bam) in cmd_str + # writing the extracted reads to the fastq and the reads to kmc + assert "pgutil extract" in cmd_str + ext_fastq = self.mock_dir / "sample-extract.fq.gz" + assert f"-o {ext_fastq}" in cmd_str + assert "-a -" in cmd_str + # concatenated with the long reads + assert "cat " in cmd_str + assert "samtools fasta" in cmd_str + assert str(self.mock_lr_bam) in cmd_str + assert "-fa /dev/stdin" in cmd_str + assert "-k29" in cmd_str + assert "-m30" in cmd_str + + def test_aligned_extract_kmc_memory(self): + """`--kmer_memory` is passed to the single-pass KMC""" + pipeline = self.create_aligned_pipeline() + pipeline.kmer_memory = 64 + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "extract-kmc").shell) + assert "-m64" in cmd_str + + def test_aligned_dag_dependencies(self): + """The pangenome is built from the extracted k-mers""" + pipeline = self.create_aligned_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + assert self._get_dep_names(dag, all_jobs, "extract-kmc") == { + "extract-kmc-symlink" + } + assert self._get_dep_names(dag, all_jobs, "vg-haplotypes") == { + "extract-kmc" + } + assert "extract-kmc" in self._get_dep_names(dag, all_jobs, "mm2-lift") + assert self._get_dep_names(dag, all_jobs, "dnascope-raw") == { + "mm2-lift" + } + + def test_aligned_lift_output(self): + """The lifted alignment is the final short-read output""" + pipeline = self.create_aligned_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + lift_cram = str(self.mock_vcf).replace(".vcf.gz", "_lift_deduped.cram") + cmd_str = str(self._get_job(all_jobs, "mm2-lift").shell) + assert f"-o {lift_cram}" in cmd_str + # The readgroup is seeded from the first input readgroup + assert "ID:sr-rg1-pg" in cmd_str + assert "LR:2" in cmd_str + + def test_aligned_replace_rg(self): + """The aligned short reads are rewritten with LR:0""" + pipeline = self.create_aligned_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + lift_cram = str(self.mock_vcf).replace(".vcf.gz", "_lift_deduped.cram") + sr_arg = r"sr-rg1=ID:sr-rg1\tSM:sample1\tLR:0" + lr_arg = r"lr-rg1=ID:lr-rg1\tSM:sample1\tLR:1" + + cmd_str = str(self._get_job(all_jobs, "dnascope-raw").shell) + assert cmd_str.count("--replace_rg") == 2 + # Each row precedes the input file it applies to + assert cmd_str.index(sr_arg) < cmd_str.index(str(self.mock_sr_bam)) + assert cmd_str.index(lr_arg) < cmd_str.index(str(self.mock_lr_bam)) + # The calling inputs are ordered short reads, lifted, long reads + assert ( + cmd_str.index(str(self.mock_sr_bam)) + < cmd_str.index(lift_cram) + < cmd_str.index(str(self.mock_lr_bam)) + ) + # The lifted alignment carries its LR tag already + assert "LR:2" not in cmd_str + + def test_aligned_bam_format(self): + """--bam_format switches the lifted output to BAM""" + pipeline = self.create_aligned_pipeline() + pipeline.bam_format = True + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + lift_bam = str(self.mock_vcf).replace(".vcf.gz", "_lift_deduped.bam") + assert lift_bam in str(self._get_job(all_jobs, "mm2-lift").shell) + assert lift_bam in str(self._get_job(all_jobs, "dnascope-raw").shell) + + # Short-read input validation + + def test_validate_sr_inputs_fastq(self): + """fastq input with a readgroup is accepted""" + pipeline = self.create_pipeline() + pipeline.readgroup = r"@RG\tID:rg1\tSM:sample1" + pipeline.validate_sr_inputs() + + def test_validate_sr_inputs_aligned(self): + """Aligned input without a readgroup is accepted""" + pipeline = self.create_aligned_pipeline() + pipeline.validate_sr_inputs() + + def test_validate_sr_inputs_both(self): + """fastq and aligned short reads cannot be combined""" + pipeline = self.create_pipeline() + pipeline.readgroup = r"@RG\tID:rg1\tSM:sample1" + pipeline.sr_aln = [self.mock_sr_bam] + with pytest.raises(SystemExit): + pipeline.validate_sr_inputs() + + def test_validate_sr_inputs_neither(self): + """Short reads are required""" + pipeline = self.create_pipeline() + pipeline.r1_fastq = [] + pipeline.r2_fastq = [] + with pytest.raises(SystemExit): + pipeline.validate_sr_inputs() + + def test_validate_sr_inputs_fastq_without_readgroup(self): + """fastq input requires a readgroup""" + pipeline = self.create_pipeline() + pipeline.readgroup = None + with pytest.raises(SystemExit): + pipeline.validate_sr_inputs() + + def test_validate_sr_inputs_fastq_length_mismatch(self): + """The r1 and r2 fastq lists must have the same length""" + pipeline = self.create_pipeline() + pipeline.readgroup = r"@RG\tID:rg1\tSM:sample1" + pipeline.r2_fastq = [] + with pytest.raises(SystemExit): + pipeline.validate_sr_inputs() + + def test_validate_sr_inputs_aligned_with_readgroup(self): + """`--readgroup` cannot be used with aligned input""" + pipeline = self.create_aligned_pipeline() + pipeline.readgroup = r"@RG\tID:rg1\tSM:sample1" + with pytest.raises(SystemExit): + pipeline.validate_sr_inputs() + + # Unaligned (uBAM/uCRAM) long-read input + + def test_lr_realign_job(self): + """The long-read input is realigned with minimap2""" + pipeline = self.create_lr_realign_pipeline() + dag = pipeline.build_dag() + job_names, all_jobs = self._get_all_job_names(dag) + + assert "bam-realign-0" in job_names + cmd_str = str(self._get_job(all_jobs, "bam-realign-0").shell) + # The input reference decodes the input file + assert f"samtools fastq --reference {self.mock_lr_ref}" in cmd_str + # The long-read minimap2 model aligns to the linear reference + assert "minimap2_lr.model" in cmd_str + assert str(self.mock_ref) in cmd_str + realigned = str(self.mock_vcf).replace(".vcf.gz", "_mm2_sorted_0.cram") + assert f"-o {realigned}" in cmd_str + assert self._get_dep_names(dag, all_jobs, "bam-realign-0") == set() + + def test_lr_realign_downstream(self): + """Downstream jobs consume the realigned long reads""" + pipeline = self.create_lr_realign_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + realigned = str(self.mock_vcf).replace(".vcf.gz", "_mm2_sorted_0.cram") + for job_name in ( + "longreadsv", + "graph-update-raw", + "graph-update", + "pangenome-sv", + "dnascope-raw", + ): + cmd_str = str(self._get_job(all_jobs, job_name).shell) + assert realigned in cmd_str, job_name + assert str(self.mock_lr_bam) not in cmd_str, job_name + + for job_name in ("longreadsv", "graph-update-raw", "dnascope-raw"): + assert "bam-realign-0" in self._get_dep_names( + dag, all_jobs, job_name + ), job_name + + # The readgroups of the realigned input are unchanged + cmd_str = str(self._get_job(all_jobs, "dnascope-raw").shell) + assert r"lr-rg1=ID:lr-rg1\tSM:sample1\tLR:1" in cmd_str + + def test_lr_realign_kmc_reads_original_input(self): + """K-mer counting reads the original long-read input""" + pipeline = self.create_lr_realign_pipeline() + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "kmc").shell) + assert f"samtools fasta --reference {self.mock_lr_ref}" in cmd_str + assert str(self.mock_lr_bam) in cmd_str + assert "_mm2_sorted_0" not in cmd_str + assert self._get_dep_names(dag, all_jobs, "kmc") == set() + + def test_lr_realign_without_input_ref(self): + """The target reference decodes the input without `lr_input_ref`""" + pipeline = self.create_lr_realign_pipeline() + pipeline.lr_input_ref = None + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + cmd_str = str(self._get_job(all_jobs, "kmc").shell) + assert f"samtools fasta --reference {self.mock_ref}" in cmd_str + + def test_aligned_sr_with_lr_realign(self): + """Aligned short reads combine with realigned long reads""" + pipeline = self.create_lr_realign_pipeline(aligned_sr=True) + dag = pipeline.build_dag() + job_names, all_jobs = self._get_all_job_names(dag) + + assert "bam-realign-0" in job_names + assert "extract-kmc" in job_names + assert "kmc" not in job_names + assert "dedup-bwa" not in job_names + + realigned = str(self.mock_vcf).replace(".vcf.gz", "_mm2_sorted_0.cram") + lift_cram = str(self.mock_vcf).replace(".vcf.gz", "_lift_deduped.cram") + cmd_str = str(self._get_job(all_jobs, "dnascope-raw").shell) + assert ( + cmd_str.index(str(self.mock_sr_bam)) + < cmd_str.index(lift_cram) + < cmd_str.index(realigned) + ) + assert self._get_dep_names(dag, all_jobs, "dnascope-raw") == { + "mm2-lift", + "bam-realign-0", + } + + # k-mer counting still reads the original long-read input + extract_cmd = str(self._get_job(all_jobs, "extract-kmc").shell) + assert str(self.mock_lr_bam) in extract_cmd + assert f"--reference {self.mock_lr_ref}" in extract_cmd From e6d526fce1091b6095257f4cf849ad73a43ed456 Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sun, 16 Aug 2026 00:23:54 -0700 Subject: [PATCH 03/10] Sort SV update regions in reference contig order Replace `sort -k1,1 -k2,2n` with `bedtools sort -faidx` in the LongReadSV BED generation so the regions sort correctly for references with an unusual contig order. Co-Authored-By: Claude Fable 5 --- sentieon_cli/command_strings.py | 3 ++- sentieon_cli/hybrid_pangenome.py | 2 +- tests/unit/test_hybrid_pangenome.py | 3 ++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/sentieon_cli/command_strings.py b/sentieon_cli/command_strings.py index 6a091ec..e55ab93 100644 --- a/sentieon_cli/command_strings.py +++ b/sentieon_cli/command_strings.py @@ -1581,11 +1581,12 @@ def cmd_hybrid_extract_kmc( def cmd_longread_sv_bed( out_bed: pathlib.Path, sv_vcf: pathlib.Path, + ref_fai: pathlib.Path, ) -> Pipeline: """Generate a BED file of graph update regions from LongReadSV calls""" zcat_cmd = Command("zcat", str(sv_vcf)) awk_cmd = Command("awk", "-F\t", LONGREAD_SV_BED_AWK, "OFS=\t") - sort_cmd = Command("sort", "-k1,1", "-k2,2n") + sort_cmd = Command("bedtools", "sort", "-faidx", str(ref_fai), "-i", "-") merge_cmd = Command("bedtools", "merge") return Pipeline( zcat_cmd, diff --git a/sentieon_cli/hybrid_pangenome.py b/sentieon_cli/hybrid_pangenome.py index a5a0d62..d8d60be 100644 --- a/sentieon_cli/hybrid_pangenome.py +++ b/sentieon_cli/hybrid_pangenome.py @@ -625,7 +625,7 @@ def build_dag(self) -> DAG: longreadsv_job = self.build_longreadsv_job(longread_sv_vcf, calling_lr) dag.add_job(longreadsv_job, realign_jobs) sv_bed_job = Job( - cmds.cmd_longread_sv_bed(sv_bed, longread_sv_vcf), + cmds.cmd_longread_sv_bed(sv_bed, longread_sv_vcf, ref_fai), "longread-sv-bed", 0, task_name="pangenome-update", diff --git a/tests/unit/test_hybrid_pangenome.py b/tests/unit/test_hybrid_pangenome.py index 9953161..fc7b6f3 100644 --- a/tests/unit/test_hybrid_pangenome.py +++ b/tests/unit/test_hybrid_pangenome.py @@ -327,7 +327,8 @@ def test_longreadsv_and_bed(self): bed_cmd = str(self._get_job(all_jobs, "longread-sv-bed").shell) assert 'n=split($10,a,":")' in bed_cmd - assert "sort -k1,1 -k2,2n" in bed_cmd + # Sort in reference contig order for bedtools merge + assert f"bedtools sort -faidx {self.mock_ref}.fai -i -" in bed_cmd assert "bedtools merge" in bed_cmd assert "sample-sv.bed" in bed_cmd # The awk script matches the validated implementation From 390c5cfe2d5444eaba6588954aaf18a6f03bd5ca Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sun, 16 Aug 2026 10:01:55 -0700 Subject: [PATCH 04/10] Thread the pangenome contig prefix through all graph consumers Add `--pangenome_contig_prefix` (default GRCh38#0#) and pass it to PGHapUpdateAlgo, pgutil lift, and PangenomeSV so a non-default pangenome reference name uses a consistent prefix throughout. Co-Authored-By: Claude Fable 5 --- sentieon_cli/hybrid_pangenome.py | 16 +++++++++++++++- tests/unit/test_hybrid_pangenome.py | 24 ++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/sentieon_cli/hybrid_pangenome.py b/sentieon_cli/hybrid_pangenome.py index d8d60be..8207b68 100644 --- a/sentieon_cli/hybrid_pangenome.py +++ b/sentieon_cli/hybrid_pangenome.py @@ -120,6 +120,12 @@ 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).", @@ -185,6 +191,7 @@ def __init__(self) -> None: self.bed: Optional[pathlib.Path] = None 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" @@ -971,7 +978,12 @@ def build_graph_update_job( input=list(lr_aln), ) driver.add_algo( - PGHapUpdateAlgo(out_gfa, gfa_file=in_gfa, target_bed=bed) + PGHapUpdateAlgo( + out_gfa, + gfa_file=in_gfa, + target_bed=bed, + prefix=self.pangenome_contig_prefix, + ) ) return Job( Pipeline(Command(*driver.build_cmd())), @@ -1040,6 +1052,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, ), "mm2-lift", self.cores, @@ -1100,6 +1113,7 @@ def build_pangenomesv_job( out_vcf, gfa_file=pangenome_gfa, min_af=PANGENOME_SV_MIN_AF, + prefix=self.pangenome_contig_prefix, ) ) return Job( diff --git a/tests/unit/test_hybrid_pangenome.py b/tests/unit/test_hybrid_pangenome.py index fc7b6f3..7307c4e 100644 --- a/tests/unit/test_hybrid_pangenome.py +++ b/tests/unit/test_hybrid_pangenome.py @@ -304,6 +304,8 @@ def test_graph_update_commands(self): assert "sample-pangenome-raw.gfa" in raw_cmd assert "--target_bed" not in raw_cmd assert str(self.mock_lr_bam) in raw_cmd + assert "--prefix" in raw_cmd + assert "GRCh38#0#" in raw_cmd update_cmd = str(self._get_job(all_jobs, "graph-update").shell) assert "--algo PGHapUpdateAlgo" in update_cmd @@ -311,6 +313,8 @@ def test_graph_update_commands(self): assert "sample-pangenome-raw.gfa" in update_cmd assert "--target_bed " in update_cmd assert "sample-sv.bed" in update_cmd + assert "--prefix" in update_cmd + assert "GRCh38#0#" in update_cmd def test_longreadsv_and_bed(self): """LongReadSV runs on the long reads and the awk BED script is @@ -385,11 +389,31 @@ def test_pangenomesv_command(self): assert "--gfa_file" in cmd_str assert "sample-pangenome.gfa" in cmd_str assert "--min_af 0.1" in cmd_str + assert "--prefix" in cmd_str + assert "GRCh38#0#" in cmd_str sv_vcf = str(self.mock_vcf).replace(".vcf.gz", "_sv.vcf.gz") assert sv_vcf in cmd_str # SV calling is not restricted to the small-variant BED assert f"--interval {self.mock_bed}" not in cmd_str + def test_pangenome_contig_prefix(self): + """`--pangenome_contig_prefix` reaches every graph consumer""" + pipeline = self.create_pipeline() + pipeline.pangenome_contig_prefix = "CHM13#0#" + dag = pipeline.build_dag() + _, all_jobs = self._get_all_job_names(dag) + + for job_name in ( + "mm2-lift", + "graph-update-raw", + "graph-update", + "pangenome-sv", + ): + cmd_str = str(self._get_job(all_jobs, job_name).shell) + assert "--prefix" in cmd_str, job_name + assert "CHM13#0#" in cmd_str, job_name + assert "GRCh38#0#" not in cmd_str, job_name + def test_dnascope_command(self): """DNAscope runs with the model, interval, and pcr_indel_model""" pipeline = self.create_pipeline() From f7b4a0f0e9762cf40e4738bf5904eb5c9f7ccdb3 Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sun, 16 Aug 2026 15:05:08 -0700 Subject: [PATCH 05/10] Harden input readgroup validation and the KMC patch check Reject BAM/CRAM inputs with no @RG header lines (previously the LR attribute was silently never applied, or build_dag crashed), require readgroup IDs to be unique across all inputs, and report malformed @RG header lines with the offending file instead of a traceback. check_kmc_patch now parses the k-mer count from the KMC output and fails on a zero count, catching KMC builds that read nothing from piped stdin while still exiting successfully. Co-Authored-By: Claude Fable 5 --- sentieon_cli/hybrid_pangenome.py | 107 ++++++++++++++++++++++--- sentieon_cli/util.py | 20 ++++- tests/unit/test_hybrid_pangenome.py | 120 ++++++++++++++++++++++++++++ tests/unit/test_kmc_patch.py | 62 +++++++++++++- 4 files changed, 290 insertions(+), 19 deletions(-) diff --git a/sentieon_cli/hybrid_pangenome.py b/sentieon_cli/hybrid_pangenome.py index 8207b68..0f33b8e 100644 --- a/sentieon_cli/hybrid_pangenome.py +++ b/sentieon_cli/hybrid_pangenome.py @@ -450,28 +450,51 @@ def collect_readgroups(self) -> None: # Read the readgroups from the ORIGINAL long-read input; minimap2 # realignment preserves the readgroup IDs with `addreplacerg` - self.sr_readgroups = [] - for aln in self.sr_aln: - self.sr_readgroups.append([]) - for rg_line in cmds.get_rg_lines(aln, self.dry_run): - self.sr_readgroups[-1].append(parse_rg_line(rg_line)) - - self.lr_readgroups = [] - for aln in self.lr_aln: - self.lr_readgroups.append([]) - for rg_line in cmds.get_rg_lines(aln, self.dry_run): - self.lr_readgroups[-1].append(parse_rg_line(rg_line)) + self.sr_readgroups = [ + self.parse_aln_readgroups(aln) for aln in self.sr_aln + ] + self.lr_readgroups = [ + self.parse_aln_readgroups(aln) for aln in self.lr_aln + ] + + def parse_aln_readgroups(self, aln: pathlib.Path) -> List[Dict[str, str]]: + """Parse the @RG lines from the header of an input alignment""" + parsed: List[Dict[str, str]] = [] + for rg_line in cmds.get_rg_lines(aln, self.dry_run): + try: + parsed.append(parse_rg_line(rg_line)) + except ValueError as e: + self.logger.error( + "Invalid readgroup line in '%s': '%s': %s", + aln, + rg_line, + e, + ) + sys.exit(2) + return parsed def validate_readgroups(self) -> None: """Confirm that all readgroups have a consistent SM tag""" + aln_inputs = list(self.sr_aln) + list(self.lr_aln) + aln_readgroups = self.sr_readgroups + self.lr_readgroups + + # Every input needs readgroups; the variant calling stages set the + # LR attribute of the input readgroups with `--replace_rg` + for aln, aln_rgs in zip(aln_inputs, aln_readgroups): + if not aln_rgs: + self.logger.error( + "No @RG lines found in the header of '%s'. Please add " + "readgroups to the input file", + aln, + ) + sys.exit(2) + rg_sm: Optional[str] = None if self.fastq_readgroup: rg_sm = self.fastq_readgroup.get("SM") elif self.sr_readgroups and self.sr_readgroups[0]: rg_sm = self.sr_readgroups[0][0].get("SM") - aln_inputs = list(self.sr_aln) + list(self.lr_aln) - aln_readgroups = self.sr_readgroups + self.lr_readgroups for aln, aln_rgs in zip(aln_inputs, aln_readgroups): for aln_rg in aln_rgs: if not aln_rg.get("ID"): @@ -502,8 +525,66 @@ def validate_readgroups(self) -> None: rg_sm, ) sys.exit(2) + + # The dry-run readgroups are synthetic and identical for every + # input, so they would always collide + if not self.dry_run: + self.validate_rg_ids(aln_inputs, aln_readgroups) + + if not self.rgsm and not rg_sm: + self.logger.error( + "Could not determine the sample name from the inputs. " + "Please set the `--rgsm` argument" + ) + sys.exit(2) self.sample_sm = self.rgsm if self.rgsm else str(rg_sm) + def validate_rg_ids( + self, + aln_inputs: List[pathlib.Path], + aln_readgroups: List[List[Dict[str, str]]], + ) -> None: + """Confirm that every input readgroup ID is unique. + + The inputs are merged by readgroup ID during variant calling, so a + collision would give one readgroup two conflicting `--replace_rg` + rewrites of the LR attribute. + """ + rg_sources: Dict[str, str] = {} + if self.fastq_readgroup: + rg_sources[self.fastq_readgroup["ID"]] = ( + "the `--readgroup` argument" + ) + for aln, aln_rgs in zip(aln_inputs, aln_readgroups): + for aln_rg in aln_rgs: + rg_id = aln_rg["ID"] + if rg_id in rg_sources: + self.logger.error( + "Duplicate readgroup ID '%s' found in %s and '%s'. " + "Readgroup IDs need to be unique across all inputs", + rg_id, + rg_sources[rg_id], + aln, + ) + sys.exit(2) + rg_sources[rg_id] = f"'{aln}'" + + # The lifted alignment reuses the first short-read readgroup ID + # with a '-pg' suffix + lift_rg = self.fastq_readgroup + if not lift_rg and self.sr_readgroups and self.sr_readgroups[0]: + lift_rg = self.sr_readgroups[0][0] + if lift_rg: + lift_id = lift_rg["ID"] + "-pg" + if lift_id in rg_sources: + self.logger.error( + "The readgroup ID '%s' in %s is reserved for the lifted " + "alignment. Please rename the input readgroup", + lift_id, + rg_sources[lift_id], + ) + sys.exit(2) + def configure(self) -> None: """Configure pipeline parameters""" pass diff --git a/sentieon_cli/util.py b/sentieon_cli/util.py index fe9cc73..f03efaa 100644 --- a/sentieon_cli/util.py +++ b/sentieon_cli/util.py @@ -288,7 +288,8 @@ def check_kmc_patch(kmc_cmd: str = "kmc") -> bool: temp_path = pathlib.Path(temp_dir) output_prefix = temp_path / "kmc_test" - # Test input sequence + # Test input sequence. Both reads have an N-free stretch longer + # than the k-mer size, so a working KMC counts k-mers from them. test_input = ( ">206B4ABXX100825:7:1:1360:6029/1\n" "TGATTTTNNNNNNNNNNNTGAAGAACGCACCCATGTTAAAGAGCATGACAAANNNANNACAAGGCTAAGNGGCGNG\n" # noqa: E501 @@ -309,14 +310,25 @@ def check_kmc_patch(kmc_cmd: str = "kmc") -> bool: ] try: - sp.run( + res = sp.run( cmd, input=test_input, text=True, check=True, - stdout=sp.DEVNULL, + stdout=sp.PIPE, stderr=sp.DEVNULL, ) - return True except (sp.CalledProcessError, FileNotFoundError): return False + + # KMC prints a stats block after a successful run + match = re.search( + r"Total no\. of k-mers\s*:\s*(\d+)", res.stdout or "" + ) + if not match: + logger.debug( + "Could not find the k-mer count in the `%s` output", + kmc_cmd, + ) + return False + return int(match.group(1)) > 0 diff --git a/tests/unit/test_hybrid_pangenome.py b/tests/unit/test_hybrid_pangenome.py index 7307c4e..42c836c 100644 --- a/tests/unit/test_hybrid_pangenome.py +++ b/tests/unit/test_hybrid_pangenome.py @@ -15,6 +15,7 @@ 0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) ) +from sentieon_cli import command_strings as cmds from sentieon_cli.hybrid_pangenome import HybridPangenome from sentieon_cli.command_strings import LONGREAD_SV_BED_AWK from sentieon_cli.dag import DAG @@ -817,3 +818,122 @@ def test_aligned_sr_with_lr_realign(self): extract_cmd = str(self._get_job(all_jobs, "extract-kmc").shell) assert str(self.mock_lr_bam) in extract_cmd assert f"--reference {self.mock_lr_ref}" in extract_cmd + + # Readgroup validation + # + # The readgroups are read from real input headers, so these tests set + # the parsed readgroups directly (or patch `get_rg_lines`) rather than + # relying on the synthetic readgroup of a dry run. + + def create_rg_pipeline(self): + """A pipeline with aligned inputs and readgroup checks enabled""" + pipeline = self.create_aligned_pipeline() + pipeline.dry_run = False + pipeline.sr_readgroups = [[{"ID": "sr-rg1", "SM": "sample1"}]] + pipeline.lr_readgroups = [[{"ID": "lr-rg1", "SM": "sample1"}]] + return pipeline + + def test_validate_readgroups_ok(self): + """Unique readgroups with a shared SM tag are accepted""" + pipeline = self.create_rg_pipeline() + pipeline.validate_readgroups() + assert pipeline.sample_sm == "sample1" + + def test_validate_readgroups_sr_without_rg_lines(self): + """An aligned short-read input without @RG lines is rejected""" + pipeline = self.create_rg_pipeline() + pipeline.sr_readgroups = [[]] + with pytest.raises(SystemExit): + pipeline.validate_readgroups() + + def test_validate_readgroups_lr_without_rg_lines(self): + """A long-read input without @RG lines is rejected""" + pipeline = self.create_rg_pipeline() + pipeline.lr_readgroups = [[]] + with pytest.raises(SystemExit): + pipeline.validate_readgroups() + + def test_validate_readgroups_no_rg_lines_with_rgsm(self): + """`--rgsm` does not bypass the missing readgroup check""" + pipeline = self.create_rg_pipeline() + pipeline.rgsm = "override_sm" + pipeline.sr_readgroups = [[]] + with pytest.raises(SystemExit): + pipeline.validate_readgroups() + + def test_validate_readgroups_duplicate_ids_across_inputs(self): + """A readgroup ID shared by the short and long reads is rejected""" + pipeline = self.create_rg_pipeline() + pipeline.sr_readgroups = [[{"ID": "rg1", "SM": "sample1"}]] + pipeline.lr_readgroups = [[{"ID": "rg1", "SM": "sample1"}]] + with pytest.raises(SystemExit): + pipeline.validate_readgroups() + + def test_validate_readgroups_duplicate_ids_within_input(self): + """A readgroup ID repeated inside one input is rejected""" + pipeline = self.create_rg_pipeline() + pipeline.sr_readgroups = [ + [{"ID": "rg1", "SM": "sample1"}, {"ID": "rg1", "SM": "sample1"}] + ] + with pytest.raises(SystemExit): + pipeline.validate_readgroups() + + def test_validate_readgroups_duplicate_with_fastq_readgroup(self): + """A long-read ID matching the `--readgroup` ID is rejected""" + pipeline = self.create_pipeline() + pipeline.dry_run = False + pipeline.fastq_readgroup = {"ID": "rg1", "SM": "sample1"} + pipeline.lr_readgroups = [[{"ID": "rg1", "SM": "sample1"}]] + with pytest.raises(SystemExit): + pipeline.validate_readgroups() + + def test_validate_readgroups_lift_id_collision(self): + """An input ID reserved for the lifted alignment is rejected""" + pipeline = self.create_pipeline() + pipeline.dry_run = False + pipeline.fastq_readgroup = {"ID": "rg1", "SM": "sample1"} + pipeline.lr_readgroups = [[{"ID": "rg1-pg", "SM": "sample1"}]] + with pytest.raises(SystemExit): + pipeline.validate_readgroups() + + def test_validate_readgroups_without_sample_name(self): + """A sample name is required for the output readgroups""" + pipeline = self.create_rg_pipeline() + pipeline.fastq_readgroup = None + pipeline.sr_aln = [] + pipeline.lr_aln = [] + pipeline.sr_readgroups = [] + pipeline.lr_readgroups = [] + with pytest.raises(SystemExit): + pipeline.validate_readgroups() + + def test_collect_readgroups_malformed_rg_line(self, monkeypatch): + """A malformed @RG line in an input header is rejected""" + pipeline = self.create_aligned_pipeline() + pipeline.dry_run = False + monkeypatch.setattr( + cmds, + "get_rg_lines", + lambda aln, dry_run: ["@RG\tID:sr-rg1\t"], + ) + with pytest.raises(SystemExit): + pipeline.collect_readgroups() + + def test_collect_readgroups_parses_input_headers(self, monkeypatch): + """The @RG lines of every input are parsed""" + pipeline = self.create_aligned_pipeline() + pipeline.dry_run = False + headers = { + str(self.mock_sr_bam): ["@RG\tID:sr-rg1\tSM:sample1"], + str(self.mock_lr_bam): ["@RG\tID:lr-rg1\tSM:sample1"], + } + monkeypatch.setattr( + cmds, + "get_rg_lines", + lambda aln, dry_run: headers[str(aln)], + ) + pipeline.collect_readgroups() + assert pipeline.sr_readgroups == [[{"ID": "sr-rg1", "SM": "sample1"}]] + assert pipeline.lr_readgroups == [[{"ID": "lr-rg1", "SM": "sample1"}]] + pipeline.validate_readgroups() + assert pipeline.sample_sm == "sample1" diff --git a/tests/unit/test_kmc_patch.py b/tests/unit/test_kmc_patch.py index 1589354..4c5f13a 100644 --- a/tests/unit/test_kmc_patch.py +++ b/tests/unit/test_kmc_patch.py @@ -3,20 +3,78 @@ """ import subprocess as sp -from unittest.mock import MagicMock, patch +from unittest.mock import patch from sentieon_cli.util import check_kmc_patch +def kmc_stdout(n_kmers: int) -> str: + """The stats block printed by KMC after a successful run""" + return ( + "1st stage: 0.068111s\n" + "2nd stage: 0.01257s\n" + "Total : 0.080681s\n" + "Tmp size : 0MB\n" + "\n" + "Stats:\n" + f" No. of k-mers below min. threshold : {n_kmers:>12}\n" + " No. of k-mers above max. threshold : 0\n" + f" No. of unique k-mers : {n_kmers:>12}\n" + " No. of unique counted k-mers : 0\n" + f" Total no. of k-mers : {n_kmers:>12}\n" + " Total no. of reads : 2\n" + " Total no. of super-k-mers : 3\n" + ) + + def test_check_kmc_patch_success(): - """Test check_kmc_patch returns True on success""" + """Test check_kmc_patch returns True when k-mers are counted""" with patch("subprocess.run") as mock_run: mock_run.return_value.returncode = 0 + mock_run.return_value.stdout = kmc_stdout(13) assert check_kmc_patch() is True mock_run.assert_called_once() args, kwargs = mock_run.call_args assert "kmc" == args[0][0] assert kwargs["input"] is not None + # The k-mer count is parsed from the captured stdout + assert kwargs["stdout"] == sp.PIPE + + +def test_check_kmc_patch_zero_kmers(): + """An unpatched KMC reads nothing from stdin but still exits 0""" + with patch("subprocess.run") as mock_run: + mock_run.return_value.returncode = 0 + mock_run.return_value.stdout = kmc_stdout(0) + assert check_kmc_patch() is False + + +def test_check_kmc_patch_missing_stats(): + """Test check_kmc_patch returns False without the k-mer count""" + with patch("subprocess.run") as mock_run: + mock_run.return_value.returncode = 0 + mock_run.return_value.stdout = "Stage 1: 100%\nStage 2: 100%\n" + assert check_kmc_patch() is False + + +def test_check_kmc_patch_no_stdout(): + """Test check_kmc_patch returns False without any output""" + with patch("subprocess.run") as mock_run: + mock_run.return_value.returncode = 0 + mock_run.return_value.stdout = None + assert check_kmc_patch() is False + + +def test_check_kmc_patch_custom_command(): + """Test check_kmc_patch runs the supplied kmc executable""" + with patch("subprocess.run") as mock_run: + mock_run.return_value.returncode = 0 + mock_run.return_value.stdout = kmc_stdout(13) + assert check_kmc_patch("/opt/kmc/kmc") is True + args, _kwargs = mock_run.call_args + assert args[0][0] == "/opt/kmc/kmc" + # KMC reads the test input from stdin + assert "/dev/stdin" in args[0] def test_check_kmc_patch_failure(): From feabf07f210ee1838107d2add3dbcbcb0084910b Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sun, 16 Aug 2026 22:14:19 -0700 Subject: [PATCH 06/10] Sort the HybridStage1 hap_bam output explicitly In Sentieon 202503.04 the HybridStage1 `--hap_bam` output is unsorted and `-` sends it to stdout. Send the hap BAM to stdout and pipe it to `sentieon util sort`, and move the algo's fastq output to a named fifo, `stage1_hap.fq`, read by bwa in the first-stage alignment. The haplotype job runs with zero scheduler threads so it always executes concurrently with the first-stage job reading the fifo. Requires sentieon driver 202503.04 or later. Co-Authored-By: Claude Fable 5 --- sentieon_cli/command_strings.py | 30 +++- sentieon_cli/dnascope_hybrid.py | 94 +++++++---- sentieon_cli/driver.py | 2 +- tests/unit/test_dnascope_hybrid_stage1.py | 187 ++++++++++++++++++++++ 4 files changed, 280 insertions(+), 33 deletions(-) create mode 100644 tests/unit/test_dnascope_hybrid_stage1.py diff --git a/sentieon_cli/command_strings.py b/sentieon_cli/command_strings.py index e55ab93..907c532 100644 --- a/sentieon_cli/command_strings.py +++ b/sentieon_cli/command_strings.py @@ -316,24 +316,46 @@ def cmd_pyexec_hybrid_anno( return Pipeline(Command(*cmd)) +def hybrid_stage1_hap( + out_hap_bam: pathlib.Path, + stage1_driver: BaseDriver, + cores: int, +) -> Pipeline: + """Sort the haplotype alignments written to stdout by HybridStage1""" + sort_cmd = Command( + "sentieon", + "util", + "sort", + "-i", + "-", + "-t", + str(cores), + "-o", + str(out_hap_bam), + # No `--sam2bam`: the algo writes unsorted BAM, not SAM, to stdout + ) + return Pipeline(Command(*stage1_driver.build_cmd()), sort_cmd) + + def hybrid_stage1( out_aln: pathlib.Path, reference: pathlib.Path, cores: int, readgroup: str, ins_driver: BaseDriver, - stage1_driver: BaseDriver, + hap_fastq_fifo: pathlib.Path, bwa_model: pathlib.Path, ) -> Pipeline: bwa_env = dict(os.environ) _ = bwa_env.pop("bwt_max_mem", None) - # Send the input of both fq commands to bwa with cat - fq1_cmd = Command(*stage1_driver.build_cmd()) + # Send both sets of reads to bwa with cat. The HybridStage1 driver of + # the `hybrid_stage1_hap` job writes its fastq output to the fifo, so + # the fifo is read first to drain it while that job runs. fq2_cmd = Command(*ins_driver.build_cmd()) cat_cmd = Command( "cat", - InputProcSub(Pipeline(fq1_cmd)), + str(hap_fastq_fifo), InputProcSub(Pipeline(fq2_cmd)), ) diff --git a/sentieon_cli/dnascope_hybrid.py b/sentieon_cli/dnascope_hybrid.py index 6be5285..012a5b0 100644 --- a/sentieon_cli/dnascope_hybrid.py +++ b/sentieon_cli/dnascope_hybrid.py @@ -51,7 +51,8 @@ CALLING_MIN_VERSIONS = { - "sentieon driver": packaging.version.Version("202503.01"), + # 202503.04 writes unsorted BAM to `--hap_bam`, supporting stdout + "sentieon driver": packaging.version.Version("202503.04"), "bedtools": None, "bcftools": packaging.version.Version("1.22"), "samtools": packaging.version.Version("1.16"), @@ -640,6 +641,8 @@ def build_dag(self) -> DAG: mapq0_slop_job, cat_merge_job, rm_job1, + stage1_fifo_job, + stage1_hap_job, stage1_job, rm_job2, second_stage_job, @@ -661,8 +664,11 @@ def build_dag(self) -> DAG: dag.add_job(mapq0_job, realign_jobs | sr_preprocessing_jobs) dag.add_job(mapq0_slop_job, {mapq0_job}) dag.add_job(cat_merge_job, {mapq0_slop_job, select_job}) - dag.add_job(stage1_job, {cat_merge_job}) - dag.add_job(second_stage_job, {stage1_job}) + # The haplotype job writes the fifo read by the first-stage job + dag.add_job(stage1_fifo_job) + dag.add_job(stage1_hap_job, {stage1_fifo_job, cat_merge_job}) + dag.add_job(stage1_job, {stage1_fifo_job, cat_merge_job}) + dag.add_job(second_stage_job, {stage1_job, stage1_hap_job}) dag.add_job(third_stage_job, {second_stage_job}) dag.add_job(call2_job, {third_stage_job}) dag.add_job(subset_job, {second_stage_job}) @@ -684,7 +690,7 @@ def build_dag(self) -> DAG: # Remove intermediate files during processing if not self.retain_tmpdir: dag.add_job(rm_job1, {cat_merge_job}) - dag.add_job(rm_job2, {stage1_job}) + dag.add_job(rm_job2, {stage1_job, stage1_hap_job}) dag.add_job(rm_job3, {second_stage_job}) dag.add_job(rm_job4, {third_stage_job}) dag.add_job(rm_job5, {concat_job}) @@ -698,27 +704,29 @@ def call_variants( rg_info: RgInfo, **_kwargs: Any, ) -> Tuple[ - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Job, - Optional[List[Job]], - Optional[Job], - Optional[Job], - Optional[Job], + Job, # call_job + Job, # select_job + Job, # mapq0_job + Job, # mapq0_slop_job + Job, # cat_merge_job + Job, # rm_job1 + Job, # stage1_fifo_job + Job, # stage1_hap_job + Job, # stage1_job + Job, # rm_job2 + Job, # second_stage_job + Job, # rm_job3 + Job, # third_stage_job + Job, # rm_job4 + Job, # call2_job + Job, # subset_job + Job, # concat_job + Job, # rm_job5 + Job, # anno_job + Optional[List[Job]], # transfer_jobs + Optional[Job], # transfer_concat_job + Optional[Job], # apply_job + Optional[Job], # norm_job ]: """ Call SNVs and indels using the DNAscope hybrid pipeline @@ -862,6 +870,17 @@ def call_variants( stage1_hap_bam = self.tmp_dir.joinpath("stage1_hap.bam") stage1_hap_bed = self.tmp_dir.joinpath("stage1_hap.bed") stage1_hap_vcf = self.tmp_dir.joinpath("stage1_hap.vcf") + + # The algo writes its fastq output to a fifo read by the bwa job + # and its unsorted haplotype alignments to stdout + stage1_fifo = self.tmp_dir.joinpath("stage1_hap.fq") + stage1_fifo_job = Job( + Pipeline(Command("mkfifo", str(stage1_fifo))), + "stage1-fifo", + 1, + task_name="hybrid-realignment", + ) + stage1_driver = Driver( reference=self.reference, thread_count=self.cores, @@ -872,13 +891,28 @@ def call_variants( ) stage1_driver.add_algo( HybridStage1( - "-", + stage1_fifo, model=self.model_bundle.joinpath("HybridStage1.model"), - hap_bam=stage1_hap_bam, + hap_bam="-", hap_bed=stage1_hap_bed, hap_vcf=stage1_hap_vcf, ) ) + stage1_hap_job = Job( + cmds.hybrid_stage1_hap( + stage1_hap_bam, + stage1_driver, + self.cores, + ), + "first-stage-hap", + # This job writes the fifo that the `first-stage` job reads, so + # the two need to run concurrently. Jobs requesting 0 threads + # start immediately instead of waiting for the thread budget; + # requesting `self.cores` here could let the scheduler serialize + # the two jobs and deadlock on the fifo. + 0, + task_name="hybrid-realignment", + ) stage1_bam = self.tmp_dir.joinpath("hybrid_stage1.bam") stage1_job = Job( @@ -888,7 +922,7 @@ def call_variants( cores=self.cores, readgroup=f"@RG\\tID:hybrid-18893\\tSM:{self.hybrid_rg_sm}", ins_driver=ins_driver, - stage1_driver=stage1_driver, + hap_fastq_fifo=stage1_fifo, bwa_model=self.model_bundle.joinpath("HybridStage1_bwa.model"), ), "first-stage", @@ -1086,6 +1120,8 @@ def call_variants( mapq0_slop_job, cat_merge_job, rm_job1, + stage1_fifo_job, + stage1_hap_job, stage1_job, rm_job2, second_stage_job, @@ -1143,6 +1179,8 @@ def call_variants( mapq0_slop_job, cat_merge_job, rm_job1, + stage1_fifo_job, + stage1_hap_job, stage1_job, rm_job2, second_stage_job, diff --git a/sentieon_cli/driver.py b/sentieon_cli/driver.py index e39fe80..c6872c3 100644 --- a/sentieon_cli/driver.py +++ b/sentieon_cli/driver.py @@ -410,7 +410,7 @@ def __init__( fa_file: Optional[pathlib.Path] = None, bed_file: Optional[pathlib.Path] = None, cut_indel: Optional[int] = None, - hap_bam: Optional[pathlib.Path] = None, + hap_bam: Optional[Union[pathlib.Path, str]] = None, hap_bed: Optional[pathlib.Path] = None, cut_len: Optional[int] = None, split_size: Optional[int] = None, diff --git a/tests/unit/test_dnascope_hybrid_stage1.py b/tests/unit/test_dnascope_hybrid_stage1.py new file mode 100644 index 0000000..6b72aeb --- /dev/null +++ b/tests/unit/test_dnascope_hybrid_stage1.py @@ -0,0 +1,187 @@ +""" +Unit tests for the DNAscope hybrid first-stage realignment jobs +""" + +import os +import pathlib +import sys +import tempfile +from unittest.mock import MagicMock, patch + +import packaging.version + +# Add the parent directory to the path to import sentieon_cli +sys.path.insert( + 0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) +) + +from sentieon_cli.dnascope_hybrid import ( + CALLING_MIN_VERSIONS, + DNAscopeHybridPipeline, +) + + +class TestDNAscopeHybridStage1: + """Test the first-stage haplotype alignment jobs""" + + def setup_method(self): + """Setup test fixtures""" + self.temp_dir = tempfile.mkdtemp() + self.mock_dir = pathlib.Path(self.temp_dir) + + self.mock_vcf = self.mock_dir / "output.vcf.gz" + self.mock_ref = self.mock_dir / "reference.fa" + self.mock_sr_aln = [self.mock_dir / "short.bam"] + self.mock_lr_aln = [self.mock_dir / "long.bam"] + self.mock_bundle = self.mock_dir / "model.bundle" + self.mock_bed = self.mock_dir / "interval.bed" + + for file_path in [ + self.mock_ref, + self.mock_sr_aln[0], + self.mock_lr_aln[0], + self.mock_bundle, + self.mock_bed, + ]: + file_path.touch() + + self.stage1_fifo = self.mock_dir / "stage1_hap.fq" + self.stage1_hap_bam = self.mock_dir / "stage1_hap.bam" + + def create_pipeline(self): + """Create a DNAscopeHybridPipeline for testing""" + with patch("sys.exit"): + pipeline = DNAscopeHybridPipeline() + + pipeline.logger = MagicMock() + + # Configure arguments + pipeline.output_vcf = self.mock_vcf + pipeline.reference = self.mock_ref + pipeline.model_bundle = self.mock_bundle + pipeline.bed = self.mock_bed + pipeline.cores = 2 + pipeline.dry_run = True + pipeline.skip_version_check = True + pipeline.tmp_dir = self.mock_dir + pipeline.pop_vcf = None + + # State normally set by validate() + pipeline.fai_data = {"chr1": {"length": 1000}} + pipeline.shards = [MagicMock()] + pipeline.shards[0].contig = "chr1" + pipeline.shards[0].start = 1 + pipeline.shards[0].stop = 1000 + pipeline.pop_vcf_contigs = {"chr1": 1000} + pipeline.lr_aln_readgroups = [[{"ID": "lr_rg1", "SM": "sample1"}]] + pipeline.sr_aln_readgroups = [[{"ID": "sr_rg1", "SM": "sample1"}]] + pipeline.hybrid_rg_sm = "sample1" + pipeline.hybrid_set_rg = False + pipeline.shortread_tech = "Illumina" + pipeline.longread_tech = "ONT" + pipeline.sr_aln = self.mock_sr_aln + pipeline.lr_aln = self.mock_lr_aln + + return pipeline + + def build_dag(self, pipeline): + """Build the DAG and index the jobs by name""" + with patch( + "sentieon_cli.dnascope_hybrid.check_version", return_value=True + ): + dag = pipeline.build_dag() + jobs = list(dag.waiting_jobs.keys()) + list(dag.ready_jobs.keys()) + return dag, {job.name: job for job in jobs} + + def dep_names(self, dag, job): + """The names of a job's dependencies""" + return {dep.name for dep in dag.waiting_jobs.get(job, set())} + + def test_fifo_job(self): + """The haplotype fastq fifo is created before the stage1 jobs""" + pipeline = self.create_pipeline() + dag, jobs = self.build_dag(pipeline) + + assert "stage1-fifo" in jobs + cmd_str = str(jobs["stage1-fifo"].shell) + assert cmd_str == f"mkfifo {self.stage1_fifo}" + assert self.dep_names(dag, jobs["stage1-fifo"]) == set() + + def test_stage1_hap_command(self): + """The haplotype alignments are sorted into stage1_hap.bam""" + pipeline = self.create_pipeline() + _dag, jobs = self.build_dag(pipeline) + + assert "first-stage-hap" in jobs + hap_job = jobs["first-stage-hap"] + cmd_str = str(hap_job.shell) + assert "--algo HybridStage1" in cmd_str + assert "HybridStage1.model" in cmd_str + # The unsorted haplotype BAM is written to stdout and sorted + assert "--hap_bam -" in cmd_str + assert "sentieon util sort" in cmd_str + assert f"-o {self.stage1_hap_bam}" in cmd_str + # The driver writes BAM, not SAM, so no conversion is needed + assert "--sam2bam" not in cmd_str + # The fastq output goes to the fifo + assert str(self.stage1_fifo) in cmd_str + + # The job must run concurrently with the fifo reader + assert hap_job.threads == 0 + + def test_stage1_reads_the_fifo(self): + """The bwa job reads the haplotype fastq from the fifo""" + pipeline = self.create_pipeline() + _dag, jobs = self.build_dag(pipeline) + + cmd_str = str(jobs["first-stage"].shell) + # The fifo is a plain `cat` argument, read before the proc sub + assert cmd_str.startswith(f"cat {self.stage1_fifo} ") + assert f"<({self.stage1_fifo}" not in cmd_str + # The haplotype driver moved to its own job; only the insertion + # driver remains in the bwa pipeline + assert "--hap_bam" not in cmd_str + assert "HybridStage1.model" not in cmd_str + assert "HybridStage1_ins.model" in cmd_str + assert "sentieon bwa mem" in cmd_str + assert "--sam2bam" in cmd_str + + def test_stage1_dependencies(self): + """The stage1 jobs wait for the fifo and the merged BED""" + pipeline = self.create_pipeline() + dag, jobs = self.build_dag(pipeline) + + assert self.dep_names(dag, jobs["first-stage-hap"]) == { + "stage1-fifo", + "concat-merge-bed", + } + assert self.dep_names(dag, jobs["first-stage"]) == { + "stage1-fifo", + "concat-merge-bed", + } + # The second stage reads both stage1 outputs + assert self.dep_names(dag, jobs["second-stage"]) == { + "first-stage", + "first-stage-hap", + } + # The haplotype VCF removed here is written by the haplotype job + assert self.dep_names(dag, jobs["rm-tmp2"]) == { + "first-stage", + "first-stage-hap", + } + + def test_second_stage_inputs(self): + """The second stage consumes the sorted haplotype BAM""" + pipeline = self.create_pipeline() + _dag, jobs = self.build_dag(pipeline) + + cmd_str = str(jobs["second-stage"].shell) + assert "--algo HybridStage2" in cmd_str + assert f"--input {self.stage1_hap_bam}" in cmd_str + assert f"--input {self.mock_dir / 'hybrid_stage1.bam'}" in cmd_str + + def test_driver_min_version(self): + """The unsorted `--hap_bam` output requires 202503.04""" + assert CALLING_MIN_VERSIONS[ + "sentieon driver" + ] >= packaging.version.Version("202503.04") From ecc766328f9320c2cb3e0e53e7c5e02000cb335d Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sun, 16 Aug 2026 23:03:17 -0700 Subject: [PATCH 07/10] Remove the temporary directory when a sentieon-pangenome run fails The two-DAG main() override performed its tmpdir cleanup inline after the last check_execution, so a failing run skipped it and leaked the directory. Wrap the body in try/finally, mirroring BasePipeline.main. Found by qualification (DEFECT_sentieon_pangenome_tmpdir_leak). Co-Authored-By: Claude Fable 5 --- sentieon_cli/sentieon_pangenome.py | 19 ++--- tests/unit/test_pipeline_lifecycle.py | 104 ++++++++++++++++++++++++++ 2 files changed, 114 insertions(+), 9 deletions(-) diff --git a/sentieon_cli/sentieon_pangenome.py b/sentieon_cli/sentieon_pangenome.py index c6816ae..938faf4 100644 --- a/sentieon_cli/sentieon_pangenome.py +++ b/sentieon_cli/sentieon_pangenome.py @@ -282,18 +282,19 @@ def main(self, args: argparse.Namespace) -> None: tmp_dir_str = tmp() self.tmp_dir = pathlib.Path(tmp_dir_str) - dag = self.build_first_dag() - executor = self.run(dag) - self.check_execution(dag, executor) - - if self.expansion_catalog or self.segdup_caller is not None: - self.get_sex(self.ploidy_json) - dag = self.build_second_dag() + try: + dag = self.build_first_dag() executor = self.run(dag) self.check_execution(dag, executor) - if not self.retain_tmpdir: - shutil.rmtree(tmp_dir_str) + if self.expansion_catalog or self.segdup_caller is not None: + self.get_sex(self.ploidy_json) + dag = self.build_second_dag() + executor = self.run(dag) + self.check_execution(dag, executor) + finally: + if not self.retain_tmpdir: + shutil.rmtree(tmp_dir_str) success = True finally: self.log_completion(success, start_time) diff --git a/tests/unit/test_pipeline_lifecycle.py b/tests/unit/test_pipeline_lifecycle.py index 593a03f..aa701e8 100644 --- a/tests/unit/test_pipeline_lifecycle.py +++ b/tests/unit/test_pipeline_lifecycle.py @@ -16,10 +16,12 @@ os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")), ) +from sentieon_cli import sentieon_pangenome # noqa: E402 from sentieon_cli.dag import DAG # noqa: E402 from sentieon_cli.exceptions import DagExecutionError # noqa: E402 from sentieon_cli.job import Job # noqa: E402 from sentieon_cli.pipeline import BasePipeline # noqa: E402 +from sentieon_cli.sentieon_pangenome import SentieonPangenome # noqa: E402 from sentieon_cli.shell_pipeline import Command, Pipeline # noqa: E402 @@ -138,6 +140,108 @@ def test_check_execution_names_the_failed_jobs(): pipeline.check_execution(DAG(), _StubExecutor([job])) +class _StubPangenome(SentieonPangenome): + """A SentieonPangenome with its run stages stubbed out. + + `fail_on` is the 1-based index of the `check_execution` call that + raises, mimicking a job failure in the first or the second DAG. + """ + + def __init__(self, fail_on=None): + super().__init__() + self.fail_on = fail_on + self.check_calls = 0 + self.dags_built = [] + + def validate_ref(self) -> None: + pass + + def validate(self) -> None: + pass + + def configure(self) -> None: + pass + + def build_first_dag(self) -> DAG: + self.dags_built.append("first") + self.ploidy_json = self.tmp_dir.joinpath("ploidy.json") + return DAG() + + def build_second_dag(self) -> DAG: + self.dags_built.append("second") + return DAG() + + def get_sex(self, ploidy_json) -> None: + pass + + def run(self, dag): + return _StubExecutor() + + def check_execution(self, dag, executor): + self.check_calls += 1 + if self.fail_on == self.check_calls: + raise DagExecutionError("a job failed") + + +def _stub_pangenome(monkeypatch, tmp_path, **kwargs): + """A stubbed pangenome pipeline writing temp dirs into `tmp_path`""" + monkeypatch.setenv("SENTIEON_TMPDIR", str(tmp_path)) + monkeypatch.setattr(sentieon_pangenome, "parse_fai", lambda fai: {}) + monkeypatch.setattr( + sentieon_pangenome, "determine_shards_from_fai", lambda fai, size: [] + ) + return _StubPangenome(**kwargs) + + +def test_pangenome_main_cleans_up_tmpdir_when_the_first_dag_fails( + tmp_path, monkeypatch +): + # SentieonPangenome overrides main() for its two-DAG flow; the temp + # directory has to be removed on the failure path too. + pipeline = _stub_pangenome(monkeypatch, tmp_path, fail_on=1) + + with pytest.raises(DagExecutionError): + pipeline.main(argparse.Namespace(loglevel="WARNING")) + + assert list(tmp_path.iterdir()) == [] # no leaked temp dir + + +def test_pangenome_main_cleans_up_tmpdir_when_the_second_dag_fails( + tmp_path, monkeypatch +): + pipeline = _stub_pangenome(monkeypatch, tmp_path, fail_on=2) + pipeline.segdup_caller = [] # request the second DAG + + with pytest.raises(DagExecutionError): + pipeline.main(argparse.Namespace(loglevel="WARNING")) + + assert pipeline.dags_built == ["first", "second"] + assert list(tmp_path.iterdir()) == [] + + +def test_pangenome_main_retains_tmpdir_on_failure_when_requested( + tmp_path, monkeypatch +): + pipeline = _stub_pangenome(monkeypatch, tmp_path, fail_on=1) + pipeline.retain_tmpdir = True + + with pytest.raises(DagExecutionError): + pipeline.main(argparse.Namespace(loglevel="WARNING")) + + assert len(list(tmp_path.iterdir())) == 1 # kept for inspection + + +def test_pangenome_main_cleans_up_tmpdir_on_success(tmp_path, monkeypatch): + pipeline = _stub_pangenome(monkeypatch, tmp_path) + pipeline.segdup_caller = [] # request the second DAG + + pipeline.main(argparse.Namespace(loglevel="WARNING")) + + assert pipeline.dags_built == ["first", "second"] + assert pipeline.check_calls == 2 + assert list(tmp_path.iterdir()) == [] + + def test_check_execution_flags_unexecuted_jobs(): pipeline = _DummyPipeline() pipeline.setup_logging(argparse.Namespace(loglevel="WARNING")) From 4f0750f2bdad6e518d3c3e042669c35a58069e49 Mon Sep 17 00:00:00 2001 From: Don Freed Date: Sun, 16 Aug 2026 23:17:44 -0700 Subject: [PATCH 08/10] Revert the k-mer count check in check_kmc_patch Restore the exit-status-only probe; the piped-KMC defect is being addressed in the KMC build itself. Co-Authored-By: Claude Fable 5 --- sentieon_cli/util.py | 20 +++--------- tests/unit/test_kmc_patch.py | 62 ++---------------------------------- 2 files changed, 6 insertions(+), 76 deletions(-) diff --git a/sentieon_cli/util.py b/sentieon_cli/util.py index f03efaa..fe9cc73 100644 --- a/sentieon_cli/util.py +++ b/sentieon_cli/util.py @@ -288,8 +288,7 @@ def check_kmc_patch(kmc_cmd: str = "kmc") -> bool: temp_path = pathlib.Path(temp_dir) output_prefix = temp_path / "kmc_test" - # Test input sequence. Both reads have an N-free stretch longer - # than the k-mer size, so a working KMC counts k-mers from them. + # Test input sequence test_input = ( ">206B4ABXX100825:7:1:1360:6029/1\n" "TGATTTTNNNNNNNNNNNTGAAGAACGCACCCATGTTAAAGAGCATGACAAANNNANNACAAGGCTAAGNGGCGNG\n" # noqa: E501 @@ -310,25 +309,14 @@ def check_kmc_patch(kmc_cmd: str = "kmc") -> bool: ] try: - res = sp.run( + sp.run( cmd, input=test_input, text=True, check=True, - stdout=sp.PIPE, + stdout=sp.DEVNULL, stderr=sp.DEVNULL, ) + return True except (sp.CalledProcessError, FileNotFoundError): return False - - # KMC prints a stats block after a successful run - match = re.search( - r"Total no\. of k-mers\s*:\s*(\d+)", res.stdout or "" - ) - if not match: - logger.debug( - "Could not find the k-mer count in the `%s` output", - kmc_cmd, - ) - return False - return int(match.group(1)) > 0 diff --git a/tests/unit/test_kmc_patch.py b/tests/unit/test_kmc_patch.py index 4c5f13a..1589354 100644 --- a/tests/unit/test_kmc_patch.py +++ b/tests/unit/test_kmc_patch.py @@ -3,78 +3,20 @@ """ import subprocess as sp -from unittest.mock import patch +from unittest.mock import MagicMock, patch from sentieon_cli.util import check_kmc_patch -def kmc_stdout(n_kmers: int) -> str: - """The stats block printed by KMC after a successful run""" - return ( - "1st stage: 0.068111s\n" - "2nd stage: 0.01257s\n" - "Total : 0.080681s\n" - "Tmp size : 0MB\n" - "\n" - "Stats:\n" - f" No. of k-mers below min. threshold : {n_kmers:>12}\n" - " No. of k-mers above max. threshold : 0\n" - f" No. of unique k-mers : {n_kmers:>12}\n" - " No. of unique counted k-mers : 0\n" - f" Total no. of k-mers : {n_kmers:>12}\n" - " Total no. of reads : 2\n" - " Total no. of super-k-mers : 3\n" - ) - - def test_check_kmc_patch_success(): - """Test check_kmc_patch returns True when k-mers are counted""" + """Test check_kmc_patch returns True on success""" with patch("subprocess.run") as mock_run: mock_run.return_value.returncode = 0 - mock_run.return_value.stdout = kmc_stdout(13) assert check_kmc_patch() is True mock_run.assert_called_once() args, kwargs = mock_run.call_args assert "kmc" == args[0][0] assert kwargs["input"] is not None - # The k-mer count is parsed from the captured stdout - assert kwargs["stdout"] == sp.PIPE - - -def test_check_kmc_patch_zero_kmers(): - """An unpatched KMC reads nothing from stdin but still exits 0""" - with patch("subprocess.run") as mock_run: - mock_run.return_value.returncode = 0 - mock_run.return_value.stdout = kmc_stdout(0) - assert check_kmc_patch() is False - - -def test_check_kmc_patch_missing_stats(): - """Test check_kmc_patch returns False without the k-mer count""" - with patch("subprocess.run") as mock_run: - mock_run.return_value.returncode = 0 - mock_run.return_value.stdout = "Stage 1: 100%\nStage 2: 100%\n" - assert check_kmc_patch() is False - - -def test_check_kmc_patch_no_stdout(): - """Test check_kmc_patch returns False without any output""" - with patch("subprocess.run") as mock_run: - mock_run.return_value.returncode = 0 - mock_run.return_value.stdout = None - assert check_kmc_patch() is False - - -def test_check_kmc_patch_custom_command(): - """Test check_kmc_patch runs the supplied kmc executable""" - with patch("subprocess.run") as mock_run: - mock_run.return_value.returncode = 0 - mock_run.return_value.stdout = kmc_stdout(13) - assert check_kmc_patch("/opt/kmc/kmc") is True - args, _kwargs = mock_run.call_args - assert args[0][0] == "/opt/kmc/kmc" - # KMC reads the test input from stdin - assert "/dev/stdin" in args[0] def test_check_kmc_patch_failure(): From 31ae7fd8aabf0e6b1a8402ae81175f45974159d9 Mon Sep 17 00:00:00 2001 From: Don Freed Date: Tue, 18 Aug 2026 07:54:48 -0700 Subject: [PATCH 09/10] Increment version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 8f866cb..d953ebd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,7 @@ [project] name = "sentieon_cli" -version = "1.7.0" +version = "2.0.0" description = "Pipeline implementations for the Sentieon software" authors = [ {name = "Don Freed", email = "don.freed@sentieon.com"}, From 5b60f8afb96cc4932afd19b01d12579f8157aca5 Mon Sep 17 00:00:00 2001 From: Don Freed Date: Tue, 18 Aug 2026 08:19:10 -0700 Subject: [PATCH 10/10] Use the SENTIEON_VERSION repo variable in CI-smoke Replace the hardcoded sentieon-version matrix value with ${{ vars.SENTIEON_VERSION }} in all three smoke jobs, matching the pattern already used by the Docker workflow. Each "Install sentieon" step now fails with an explicit error when the variable is unset, instead of a confusing tar failure on a 404 from S3. Co-Authored-By: Claude Opus 5 --- .github/workflows/main.yml | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 0e5e624..1cee559 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -15,7 +15,7 @@ jobs: max-parallel: 1 matrix: python-version: ["3.11", "3.14"] - sentieon-version: ["202503.02"] + sentieon-version: ["${{ vars.SENTIEON_VERSION }}"] os: [ubuntu-22.04] #, macos-latest, windows-latest] runs-on: ${{ matrix.os }} steps: @@ -41,6 +41,10 @@ jobs: pip install multiqc - name: Install sentieon run: | + if [ -z "$SENTIEON_VERSION" ]; then + echo "::error::Repository variable SENTIEON_VERSION is not set." >&2 + exit 1 + fi url="https://s3.amazonaws.com/sentieon-release/software/sentieon-genomics-$SENTIEON_VERSION.tar.gz" bin_dir=$(pwd)/sentieon-genomics-$SENTIEON_VERSION/bin machine_arch=$(uname -m) @@ -98,7 +102,7 @@ jobs: max-parallel: 1 matrix: python-version: ["3.11", "3.14"] - sentieon-version: ["202503.02"] + sentieon-version: ["${{ vars.SENTIEON_VERSION }}"] os: [ubuntu-22.04] #, macos-latest, windows-latest] runs-on: ${{ matrix.os }} steps: @@ -136,6 +140,10 @@ jobs: sudo chmod ugo+x /usr/local/bin/mosdepth - name: Install sentieon run: | + if [ -z "$SENTIEON_VERSION" ]; then + echo "::error::Repository variable SENTIEON_VERSION is not set." >&2 + exit 1 + fi url="https://s3.amazonaws.com/sentieon-release/software/sentieon-genomics-$SENTIEON_VERSION.tar.gz" bin_dir=$(pwd)/sentieon-genomics-$SENTIEON_VERSION/bin machine_arch=$(uname -m) @@ -220,7 +228,7 @@ jobs: max-parallel: 1 matrix: python-version: ["3.11", "3.14"] - sentieon-version: ["202503.02"] + sentieon-version: ["${{ vars.SENTIEON_VERSION }}"] os: [ubuntu-22.04] #, macos-latest, windows-latest] runs-on: ${{ matrix.os }} steps: @@ -261,6 +269,10 @@ jobs: pip install multiqc - name: Install sentieon run: | + if [ -z "$SENTIEON_VERSION" ]; then + echo "::error::Repository variable SENTIEON_VERSION is not set." >&2 + exit 1 + fi url="https://s3.amazonaws.com/sentieon-release/software/sentieon-genomics-$SENTIEON_VERSION.tar.gz" bin_dir=$(pwd)/sentieon-genomics-$SENTIEON_VERSION/bin machine_arch=$(uname -m)