Merge dev change

pyproject.toml

@@ -1,7 +1,7 @@
 
 [project]
 name = "sentieon_cli"
-version = "1.6.2"
+version = "1.6.3"
 description = "Pipeline implementations for the Sentieon software"
 authors = [
     {name = "Don Freed", email = "don.freed@sentieon.com"},

sentieon_cli/sentieon_pangenome.py

@@ -29,6 +29,7 @@
     DNAscope,
     Driver,
     GCBias,
+    GVCFtyper,
     InsertSizeMetricAlgo,
     LocusCollector,
     MeanQualityByCycle,
@@ -128,7 +129,13 @@ class SentieonPangenome(BasePangenome):
             },
             "pangenome_ref_name": {
                 "default": "GRCh38",
-                "help": "Reference name in the pangenome (GRCh38)",
+                "help": (
+                    "Reference name in the pangenome (GRCh38). The "
+                    "'extract.<pangenome_ref_name>.model' member of the model "
+                    "bundle is preferred; if it is absent and the reference "
+                    "name is 'GRCh38', the pipeline falls back to "
+                    "'extract.model'."
+                ),
             },
             "pangenome_contig_prefix": {
                 "default": "GRCh38#0#",
@@ -240,6 +247,7 @@ def __init__(self) -> None:
         self.call_svs = False
         self.gvcf = False
         self.pangenome_ref_name = "GRCh38"
+        self.extract_model_name = "extract.model"
         self.pangenome_contig_prefix = "GRCh38#0#"
         self.skip_metrics = False
         self.skip_multiqc = False
@@ -414,6 +422,13 @@ def validate_expansion(self) -> None:
         if self.expansion_catalog is None:
             return
 
+        if self.tech.upper() == "ULTIMA":
+            self.logger.error(
+                "`--expansion_catalog` is not supported with single-end "
+                "(Ultima) input."
+            )
+            sys.exit(2)
+
         if len(self.sample_input) > 1:
             self.logger.error(
                 "`--expansion_catalog` accepts only a single `--sample_input` "
@@ -503,9 +518,21 @@ def validate_bundle(self) -> None:
             )
 
         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 so that
+        # existing bundles continue to work.
+        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
+
         if (
             "dnascope.model" not in bundle_members
-            or "extract.model" not in bundle_members
+            or self.extract_model_name not in bundle_members
             or "minimap2.model" not in bundle_members
         ):
             self.logger.error(
@@ -613,6 +640,9 @@ def build_first_dag(self) -> DAG:
         self.ploidy_json = pathlib.Path(
             str(self.output_vcf).replace(".vcf.gz", "_ploidy.json")
         )
+        out_gvcf = pathlib.Path(
+            str(self.output_vcf).replace(".vcf.gz", ".g.vcf.gz")
+        )
 
         # Intermediate file paths
         bwa_bam = self.tmp_dir.joinpath("sample-bwa.bam")
@@ -675,7 +705,7 @@ def build_first_dag(self) -> DAG:
                     ext_fastq,
                     self.sample_input,
                     self.reference,
-                    self.model_bundle.joinpath("extract.model"),
+                    self.model_bundle.joinpath(self.extract_model_name),
                     self.tmp_dir,
                     rw_bam,
                     threads=self.cores,
@@ -824,10 +854,18 @@ def build_first_dag(self) -> DAG:
         dag.add_job(dnascope_job, dnascope_dependencies)
         apply_dependencies.add(dnascope_job)
 
+        # When --gvcf is set, the model-apply / transfer outputs are
+        # gVCFs; GVCFtyper produces the final VCF at self.output_vcf.
+        small_variants_out = out_gvcf if self.gvcf else self.output_vcf
+
         # transfer annotations from the pop_vcf
         if self.pop_vcf:
             transfer_jobs, concat_job = build_transfer_jobs(
-                transfer_vcf if not self.skip_model_apply else self.output_vcf,
+                (
+                    transfer_vcf
+                    if not self.skip_model_apply
+                    else (small_variants_out)
+                ),
                 self.pop_vcf,
                 raw_vcf,
                 self.tmp_dir,
@@ -842,10 +880,21 @@ def build_first_dag(self) -> DAG:
             dag.add_job(concat_job, set(transfer_jobs))
             apply_dependencies.add(concat_job)
 
+        small_variants_last_job: Optional[Job] = None
         if not self.skip_model_apply:
             # DNAModelApply
-            apply_job = self.build_dnamodelapply_job(transfer_vcf)
+            apply_job = self.build_dnamodelapply_job(
+                transfer_vcf, small_variants_out
+            )
             dag.add_job(apply_job, apply_dependencies)
+            small_variants_last_job = apply_job
+        elif self.pop_vcf:
+            small_variants_last_job = concat_job
+
+        # Genotype the gVCF to also produce a regular VCF at output_vcf
+        if self.gvcf and small_variants_last_job is not None:
+            gvcftyper_job = self.build_gvcftyper_job(self.output_vcf, out_gvcf)
+            dag.add_job(gvcftyper_job, {small_variants_last_job})
 
         if self.call_svs and sv_vcf and self.has_cnv_model:
             self._add_cnv_jobs(
@@ -889,7 +938,7 @@ def build_alignment_job(
                 self.r1_fastq,
                 self.r2_fastq,
                 "@RG\\t" + "\\t".join([f"{x[0]}:{x[1]}" for x in rg.items()]),
-                self.model_bundle.joinpath("extract.model"),
+                self.model_bundle.joinpath(self.extract_model_name),
                 self.model_bundle.joinpath("bwa.model"),
                 self.cores,
                 unzip=unzip,
@@ -1169,10 +1218,11 @@ def build_dnascope_job(
             self.cores,
         )
 
-    def build_dnamodelapply_job(self, in_vcf) -> Job:
-        if not self.output_vcf:
-            self.logger.error("output_vcf is required")
-            sys.exit(2)
+    def build_dnamodelapply_job(
+        self,
+        in_vcf: pathlib.Path,
+        out_vcf: pathlib.Path,
+    ) -> Job:
         if not self.model_bundle:
             self.logger.error("model_bundle is required")
             sys.exit(2)
@@ -1185,7 +1235,7 @@ def build_dnamodelapply_job(self, in_vcf) -> Job:
             DNAModelApply(
                 model=self.model_bundle.joinpath("dnascope.model"),
                 vcf=in_vcf,
-                output=self.output_vcf,
+                output=out_vcf,
             )
         )
         return Job(
@@ -1194,6 +1244,33 @@ def build_dnamodelapply_job(self, in_vcf) -> Job:
             self.cores,
         )
 
+    def build_gvcftyper_job(
+        self,
+        out_vcf: pathlib.Path,
+        in_gvcf: pathlib.Path,
+    ) -> Job:
+        """Genotype a gVCF to produce a VCF"""
+        if not self.reference:
+            self.logger.error("reference is required")
+            sys.exit(2)
+
+        driver = Driver(
+            reference=self.reference,
+            thread_count=self.cores,
+            interval=self.bed,
+        )
+        driver.add_algo(
+            GVCFtyper(
+                output=out_vcf,
+                vcf=in_gvcf,
+            )
+        )
+        return Job(
+            Pipeline(Command(*driver.build_cmd())),
+            "gvcftyper",
+            self.cores,
+        )
+
     def build_segdup_job(
         self,
         out_segdup: pathlib.Path,

tests/unit/test_sentieon_pangenome.py

@@ -130,6 +130,47 @@ def test_skip_model_apply(self):
         args_str = " ".join([str(arg) for arg in concat_job.shell.nodes[0].args])
         assert str(pipeline.output_vcf) in args_str
 
+    def test_gvcf_adds_gvcftyper_and_routes_outputs(self):
+        """--gvcf produces a .g.vcf.gz from model-apply and runs
+        GVCFtyper to produce the final VCF at output_vcf."""
+        pipeline = self.create_pipeline()
+        pipeline.gvcf = True
+        dag = pipeline.build_first_dag()
+
+        all_jobs = list(dag.waiting_jobs.keys()) + list(dag.ready_jobs.keys())
+        job_names = [job.name for job in all_jobs]
+        assert "model-apply" in job_names
+        assert "gvcftyper" in job_names
+
+        expected_gvcf = str(pipeline.output_vcf).replace(
+            ".vcf.gz", ".g.vcf.gz"
+        )
+
+        # DNAscope emits gVCF and model-apply writes to .g.vcf.gz
+        for name in ("dnascope-raw", "model-apply"):
+            job = next(j for j in all_jobs if j.name == name)
+            cmd_str = str(job.shell)
+            if name == "dnascope-raw":
+                assert "--emit_mode gvcf" in cmd_str
+            else:
+                assert expected_gvcf in cmd_str
+
+        # GVCFtyper reads the gVCF and writes the final output VCF
+        gvcftyper_job = next(j for j in all_jobs if j.name == "gvcftyper")
+        cmd_str = str(gvcftyper_job.shell)
+        assert "--algo GVCFtyper" in cmd_str
+        assert expected_gvcf in cmd_str
+        assert str(pipeline.output_vcf) in cmd_str
+
+    def test_no_gvcftyper_without_gvcf(self):
+        """Without --gvcf, no GVCFtyper job is added"""
+        pipeline = self.create_pipeline()
+        dag = pipeline.build_first_dag()
+
+        all_jobs = list(dag.waiting_jobs.keys()) + list(dag.ready_jobs.keys())
+        job_names = [job.name for job in all_jobs]
+        assert "gvcftyper" not in job_names
+
     def test_call_svs(self):
         """Test that PangenomeSV is added when --call_svs is enabled"""
         pipeline = self.create_pipeline()