Hybrid pangenome

.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)

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 = "[email protected]"},

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

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: