diff --git a/README.md b/README.md index 49642e9..b0db2eb 100644 --- a/README.md +++ b/README.md @@ -53,7 +53,7 @@ The `sentieon-cli` supports the following global arguments: - [**DNAscope**](https://support.sentieon.com/docs/sentieon_cli/#dnascope) - DNAscope pipeline implementation for germline SNV and indel calling from short read data. - [**DNAscope LongRead**](https://support.sentieon.com/docs/sentieon_cli/#dnascope-longread) - DNAscope LongRead pipeline implementations for germline SNV and indel calling from long read data. - [**DNAscope Hybrid**](https://support.sentieon.com/docs/sentieon_cli/#dnascope-hybrid) - DNAscope short-long-hybrid pipeline. -- [**Sentieon Pangenome**](https://support.sentieon.com/docs/sentieon_cli/#sentieon-pangenome) - Sentieon pangenome alignment and variant calling. Our recommended pipeline for short-read small variant calling. +- [**DNAscope Pangenome**](https://support.sentieon.com/docs/sentieon_cli/#sentieon-pangenome) - DNAscope pangenome alignment and variant calling. Our recommended pipeline for short-read small variant calling. ## License Unless otherwise indicated, files in this repository are licensed under a BSD 2-Clause License. diff --git a/docs/pangenome.md b/docs/pangenome.md deleted file mode 100644 index a71ca78..0000000 --- a/docs/pangenome.md +++ /dev/null @@ -1,244 +0,0 @@ -# Pangenome - -**Note: This is an experimental pangenome pipeline that is no longer recommended. For short-read pangenome alignment and variant calling, please use the [`sentieon-pangenome` pipeline](https://support.sentieon.com/docs/sentieon_cli/#sentieon-pangenome)** - -Pangenome is a pipeline for alignment and variant calling from short-read DNA sequence data using pangenome graphs. The Pangenome pipeline leverages graph-based reference representations to improve alignment and variant calling accuracy, particularly in complex genomic regions with high sequence diversity. - -The pipeline uses the open-source vg toolkit for pangenome alignment and structural variant calling. Sentieon tools are used for small variant calling, CNV detection, segdup variant calling, preprocessing and metrics collection. Specialized tools are used for HLA/KIR genotyping and repeat expansion analysis. - -The pipeline accepts as input unaligned reads in FASTQ format and outputs variants in VCF format, aligned reads in BAM format, and comprehensive quality control metrics. - -Pangenome is implemented using open-source tools and the Sentieon software package, which requires a valid license for use. Please contact info@sentieon.com for access to the Sentieon software and an evaluation license. - -## Prerequisites - -- Sentieon software package version 202503 or higher. -- [vg] toolkit for pangenome alignment operations. -- [KMC] version 3 or higher for k-mer counting. -- [samtools] version 1.16 or higher for alignment operations. -- [bcftools] version 1.10 or higher for VCF operations. -- [MultiQC] version 1.18 or higher for metrics report generation. -- [T1K] for HLA and KIR genotyping (optional). -- [ExpansionHunter] for repeat expansion calling (optional). -- [segdup-caller] for segmental duplication variant calling (optional). - -The executables will be accessed through the user's `PATH` environment variable. - -## Input data requirements - -### The Reference genome - -Pangenome will call variants present in the sample relative to a high quality reference genome sequence. Besides the reference genome file, a samtools fasta index file (.fai) needs to be present. - -### Pangenome graph files - -The pipeline requires several pangenome graph files: -- **GBZ file**: The pangenome graph in GBZ format -- **Haplotype file**: Haplotype information for the pangenome -- **Snarls file**: VG snarls file for the GBZ file (created using `vg snarls`) -- **XG file**: VG XG index file for the GBZ file (created using `vg convert`) - -### Model bundle - -A Sentieon model bundle containing machine learning models for variant calling is required. Model bundle files can be found in the [sentieon-models] repository. - -## Usage - -### Pangenome alignment and variant calling from FASTQ - -A single command is run to align reads to a pangenome graph, call small variants, structural variants, and CNVs: - -```sh -sentieon-cli pangenome [-h] \ - -r REFERENCE \ - --gbz GBZ \ - --hapl HAPL \ - --snarls SNARLS \ - --xg XG \ - -m MODEL_BUNDLE \ - --r1_fastq R1_FASTQ ... \ - --r2_fastq R2_FASTQ ... \ - --readgroups READGROUPS ... \ - [-d DBSNP] \ - [-t CORES] \ - [--kmer_memory KMER_MEMORY] \ - [--expansion_catalog EXPANSION_CATALOG] \ - [--t1k_hla_seq T1K_HLA_SEQ] \ - [--t1k_hla_coord T1K_HLA_COORD] \ - [--t1k_kir_seq T1K_KIR_SEQ] \ - [--t1k_kir_coord T1K_KIR_COORD] \ - [--segdup_caller_genes SEGDUP_CALLER_GENES] \ - [--dry_run] \ - sample.vcf.gz -``` - -The Pangenome pipeline requires the following arguments: -- `-r REFERENCE`: the location of the reference FASTA file. A reference fasta index, ".fai" file, is also required. -- `--gbz GBZ`: the pangenome graph file in GBZ format. -- `--hapl HAPL`: the haplotype file for the pangenome. -- `--snarls SNARLS`: the vg snarls file for the GBZ file. -- `--xg XG`: the XG file for the GBZ file. -- `-m MODEL_BUNDLE`: the location of the model bundle containing DNAscope and CNVscope models. -- `--r1_fastq R1_FASTQ`: the R1 input FASTQ. Can be used with multiple files. -- `--r2_fastq R2_FASTQ`: the R2 input FASTQ. Can be used with multiple files. -- `--readgroups READGROUPS`: readgroup information for each FASTQ. An example argument is, `--readgroups "@RG\tID:HG002-1\tSM:HG002\tLB:HG002-LB-1\tPL:ILLUMINA"`. -- `sample.vcf.gz`: the location of the output VCF file for small variants. The pipeline requires the output file end with the suffix, ".vcf.gz". The file path without the suffix will be used as the basename for other output files. - -The Pangenome pipeline accepts the following optional arguments: -- `-d DBSNP`: the location of the Single Nucleotide Polymorphism database (dbSNP) VCF used to label known variants. A VCF index file is required. -- `-t CORES`: number of computing threads/cores to use. Default is all available cores. -- `--kmer_memory KMER_MEMORY`: memory limit for KMC k-mer counting in GB. Default is 128 GB. -- `--expansion_catalog EXPANSION_CATALOG`: ExpansionHunter variant catalog for repeat expansion calling. -- `--t1k_hla_seq T1K_HLA_SEQ`: DNA HLA sequence FASTA file for T1K HLA calling. Requires `--t1k_hla_coord`. -- `--t1k_hla_coord T1K_HLA_COORD`: DNA HLA coordinate FASTA file for T1K HLA calling. Requires `--t1k_hla_seq`. -- `--t1k_kir_seq T1K_KIR_SEQ`: DNA KIR sequence FASTA file for T1K KIR calling. Requires `--t1k_kir_coord`. -- `--t1k_kir_coord T1K_KIR_COORD`: DNA KIR coordinate FASTA file for T1K KIR calling. Requires `--t1k_kir_seq`. -- `--segdup_caller_genes SEGDUP_CALLER_GENES`: Genes for segmental duplication calling. Example: 'CFH,CFHR3,CYP11B1,CYP2D6,GBA,NCF1,PMS2,SMN1,STRC'. -- `-h`: print the command-line help and exit. -- `--dry_run`: print the pipeline commands, but do not actually execute them. - -## Pipeline workflow - -The Pangenome pipeline consists of multiple stages executed across two DAGs: - -### First DAG -1. **K-mer counting**: `kmc` counts k-mers from input FASTQ files -2. **Sample-specific pangenome**: `vg haplotypes` creates a sample-specific pangenome using k-mer frequencies -3. **Pangenome alignment**: `vg giraffe` aligns reads to the sample-specific pangenome -4. **Read support computation**: `vg pack` computes read support for variants -5. **Structural variant calling**: `vg call` identifies structural variants from read support -6. **Linear alignment**: `vg surject` projects pangenome alignments back to linear reference -7. **Duplicate marking**: Sentieon `Dedup` marks duplicate reads and computes metrics -8. **Indel realignment**: Sentieon `Realigner` performs indel realignment -9. **Small variant calling**: Sentieon `DNAscope` calls small variants (SNVs/indels) -10. **Ploidy estimation**: Estimation of sample sex and ploidy for downstream processing -11. **HLA/KIR genotyping**: T1K genotypes HLA and KIR loci (optional) -12. **Segmental duplication calling**: segdup-caller identifies variants in segmental duplications (optional) - -### Second DAG -1. **CNV calling**: Sentieon CNVscope calls copy number variants with sex-specific processing -2. **Repeat expansion calling**: ExpansionHunter identifies repeat expansions (optional) - -## Pipeline output - -### List of output files - -The following files are output when processing WGS FASTQ: -- `sample.vcf.gz`: Small variant calls (SNVs and indels) from DNAscope. -- `sample_pangenome-aligned.bam`: Aligned, coordinate-sorted and duplicate-marked read data projected to the linear reference. -- `sample_svs.vcf.gz`: Structural variant calls from vg call. -- `sample_cnv.vcf.gz`: Copy number variant calls from CNVscope. -- `sample_ploidy.json`: Estimated sample sex and ploidy data. -- `sample_metrics`: A directory containing QC metrics for the analyzed sample. - - `sample_metrics/{sample}.txt.alignment_stat.txt`: Metrics from the AlignmentStat algo. - - `sample_metrics/{sample}.txt.base_distribution_by_cycle.txt`: Metrics from the BaseDistributionByCycle algo. - - `sample_metrics/{sample}.txt.dedup_metrics.txt`: Metrics from the Dedup algo. - - `sample_metrics/{sample}.txt.gc_bias*`: Metrics from the GCBias algo. - - `sample_metrics/{sample}.txt.insert_size.txt`: Metrics from the InsertSizeMetricAlgo algo. - - `sample_metrics/{sample}.txt.mean_qual_by_cycle.txt`: Metrics from the MeanQualityByCycle algo. - - `sample_metrics/{sample}.txt.qual_distribution.txt`: Metrics from the QualDistribution algo. - - `sample_metrics/{sample}.txt.wgs.txt`: Metrics from the WgsMetricsAlgo algo. - - `sample_metrics/multiqc_report.html`: Collected QC metrics aggregated by MultiQC. - -### Optional output files - -When optional analyses are enabled: -- `sample_hla/`: HLA genotyping results from T1K (when HLA files are provided). -- `sample_kir/`: KIR genotyping results from T1K (when KIR files are provided). -- `sample_segdups/`: Segmental duplication variant calls (when genes are specified). -- `sample_expansion*`: Repeat expansion calls from ExpansionHunter (when a catalog is provided). - -[vg]: https://github.com/vgteam/vg -[KMC]: https://github.com/refresh-bio/KMC -[samtools]: https://www.htslib.org/ -[bcftools]: http://samtools.github.io/bcftools/bcftools.html -[MultiQC]: https://multiqc.info/ -[T1K]: https://github.com/mourisl/T1K -[ExpansionHunter]: https://github.com/Illumina/ExpansionHunter -[segdup-caller]: https://github.com/Sentieon/segdup-caller -[sentieon-models]: https://github.com/Sentieon/sentieon-models - -## Example data processing - -### Obtaining the pangenome graph files - -Download vg the public human pangenome from the HPRC repository - -```sh -curl -LO "https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.gbz" -curl -LO "https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.hapl" -``` - -Generate a snarls file from the .gbz with `vg snarls`. -```sh -vg snarls hprc-v1.1-mc-grch38.gbz > hprc-v1.1-mc-grch38.snarls -``` - -Generate a xg file from the .gbz with `vg convert`. -```sh -vg convert -x --drop-haplotypes hprc-v1.1-mc-grch38.gbz > hprc-v1.1-mc-grch38.xg -``` - -### Download the hg38 reference genome - -Download the reference fasta file and samtools faidx file -```sh -curl -L 'https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz' | \ - gzip -dc > GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna -curl -LO 'https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.fai' -``` - -### Download an expansion catalog of clinically relevant expansions - -```sh -curl -L 'https://github.com/Illumina/ExpansionHunter/raw/refs/tags/v5.0.0/variant_catalog/hg38/variant_catalog.json' > hg38_variant_catalog.json -``` - -### Download t1k fasta files - -Download the two databases -```sh -perl t1k-build.pl -o hlaidx --download IPD-IMGT/HLA -perl t1k-build.pl -o kiridx --download IPD-KIR -``` - -Download the gencode gtf file for hg38 -```sh -curl -L 'https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.annotation.gtf.gz' | \ - gzip -dc > gencode.v38.annotation.gtf -``` - -Bulid the coordinate file -```sh -perl t1k-build.pl -o hlaidx -d hlaidx/hla.dat -g gencode.v38.annotation.gtf -perl t1k-build.pl -o kiridx -d kiridx/kir.dat -g gencode.v38.annotation.gtf -``` - -### Download a 30x WGS fastq - -```sh -curl -LO 'https://storage.googleapis.com/brain-genomics-public/research/sequencing/fastq/novaseq/wgs_pcr_free/30x/HG002.novaseq.pcr-free.30x.R1.fastq.gz' -curl -LO 'https://storage.googleapis.com/brain-genomics-public/research/sequencing/fastq/novaseq/wgs_pcr_free/30x/HG002.novaseq.pcr-free.30x.R2.fastq.gz' -``` - -### Run the pangenome pipeline - -```sh -sentieon-cli -v pangenome \ - --reference GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \ - --gbz hprc-v1.1-mc-grch38.gbz \ - --hapl hprc-v1.1-mc-grch38.hapl \ - --snarls hprc-v1.1-mc-grch38.snarls \ - --xg hprc-v1.1-mc-grch38.xg \ - --model_bundle pangenome.bundle \ - --r1_fastq HG002.novaseq.pcr-free.30x.R1.fastq.gz \ - --r2_fastq HG002.novaseq.pcr-free.30x.R2.fastq.gz \ - --readgroups '@RG\tID:HG002-1\tSM:HG002\tPL:ILLUMINA' \ - --expansion_catalog hg38_variant_catalog.json \ - --segdup_caller_genes 'CFH,CFHR3,CYP11B1,CYP2D6,GBA,NCF1,PMS2,SMN1,STRC' \ - --t1k_hla_seq hlaidx/data_dna_seq.fa \ - --t1k_hla_coord hlaidx/data_dna_coord.fa \ - --t1k_kir_seq kiridx/data_dna_seq.fa \ - --t1k_kir_coord kiridx/data_dna_coord.fa \ - HG002_pangenome_analysis.vcf.gz -``` diff --git a/sentieon_cli/__init__.py b/sentieon_cli/__init__.py index 54b80b5..3fcb34b 100644 --- a/sentieon_cli/__init__.py +++ b/sentieon_cli/__init__.py @@ -2,7 +2,6 @@ from .dnascope import DNAscopePipeline from .dnascope_hybrid import DNAscopeHybridPipeline from .dnascope_longread import DNAscopeLRPipeline -from .pangenome import PangenomePipeline from .sentieon_pangenome import SentieonPangenome from .util import __version__ @@ -51,17 +50,13 @@ def main(): pipeline.add_arguments(dnascope_hybrid_subparser) dnascope_hybrid_subparser.set_defaults(pipeline=pipeline.main) - # Pangenome parser - pipeline = PangenomePipeline() - pangenome_subparser = subparsers.add_parser("pangenome") - pipeline.add_arguments(pangenome_subparser) - pangenome_subparser.set_defaults(pipeline=pipeline.main) - - # Sentieon pangenome + # DNAscope pangenome pipeline = SentieonPangenome() - sentieon_pangenome_subparser = subparsers.add_parser("sentieon-pangenome") - pipeline.add_arguments(sentieon_pangenome_subparser) - sentieon_pangenome_subparser.set_defaults(pipeline=pipeline.main) + dnascope_pangenome_subparser = subparsers.add_parser( + "dnascope-pangenome", aliases=["sentieon-pangenome"] + ) + pipeline.add_arguments(dnascope_pangenome_subparser) + dnascope_pangenome_subparser.set_defaults(pipeline=pipeline.main) args = parser.parse_args() args.pipeline(args) diff --git a/sentieon_cli/pangenome.py b/sentieon_cli/pangenome.py deleted file mode 100644 index 5c02025..0000000 --- a/sentieon_cli/pangenome.py +++ /dev/null @@ -1,1233 +0,0 @@ -""" -Pangenome alignment and variant calling pipeline -""" - -import argparse -import copy -import itertools -import json -import os -import pathlib -import shutil -import sys -from typing import List, Optional, Tuple - -import packaging.version - -from importlib.resources import files - -from . import command_strings as cmds -from .archive import ar_load -from .base_pangenome import BasePangenome, SampleSex -from .dag import DAG -from .driver import ( - AlignmentStat, - BaseDistributionByCycle, - CNVscope, - CNVModelApply, - Dedup, - DNAModelApply, - DNAscope, - Driver, - GCBias, - InsertSizeMetricAlgo, - MeanQualityByCycle, - LocusCollector, - ReadWriter, - QualDistribution, - WgsMetricsAlgo, -) -from .job import Job -from .shell_pipeline import Command, Pipeline -from .util import __version__, check_version, parse_rg_line, path_arg, tmp - -PANGENOME_MIN_VERSIONS = { - "kmc": None, - "sentieon driver": packaging.version.Version("202503.01"), - "vg": None, - "bcftools": packaging.version.Version("1.22"), - "samtools": packaging.version.Version("1.16"), -} - -T1K_MIN_VERSION = { - "run-t1k": None, -} - -EXPANSION_MIN_VERSION = { - "ExpansionHunter": None, -} - -SEGDUP_MIN_VERSION = { - "segdup-caller": None, -} - - -SV_HDR_ATTR = { - "INFO": { - "AT": { - "Number": "R", - "Type": "String", - "Description": '"Allele Traversal as path in graph"', - }, - "DP": { - "Number": "1", - "Type": "Integer", - "Description": '"Total Depth"', - }, - }, - "FORMAT": { - "AD": { - "Number": ".", - "Type": "Integer", - "Description": ( - '"Allelic depths for the ref and alt alleles in the order ' - 'listed"' - ), - }, - "DP": { - "Number": "1", - "Type": "Integer", - "Description": '"Read Depth"', - }, - "GL": { - "Number": "G", - "Type": "Float", - "Description": ( - '"Genotype Likelihood, log10-scaled likelihoods of the data ' - "given the called genotype for each possible genotype " - "generated from the reference and alternate alleles given the " - 'sample ploidy"' - ), - }, - "GP": { - "Number": "1", - "Type": "Float", - "Description": ( - '"Genotype Probability, the log-scaled posterior probability ' - 'of the called genotype"' - ), - }, - "GQ": { - "Number": "1", - "Type": "Integer", - "Description": ( - '"Genotype Quality, the Phred-scaled probability estimate of ' - 'the called genotype"' - ), - }, - "GT": { - "Number": "1", - "Type": "String", - "Description": '"Genotype"', - }, - "MAD": { - "Number": "1", - "Type": "Integer", - "Description": '"Minimum site allele depth"', - }, - "XD": { - "Number": "1", - "Type": "Float", - "Description": ( - '"eXpected Depth, background coverage as used for the Poisson ' - 'model"' - ), - }, - }, - "FILTER": { - "lowad": { - "Description": ( - '"Variant does not meet minimum allele read support threshold ' - 'of 1"' - ), - }, - "lowdepth": {"Description": '"Variant has read depth less than 4"'}, - "PASS": {"Description": '"All filters passed"'}, - }, -} - - -GRCh38_CONTIGS_NOT_IN_PANGENOME = { - "chr14_KI270726v1_random", - "chr15_KI270727v1_random", - "chr1_KI270711v1_random", - "chr22_KI270737v1_random", - "chr22_KI270739v1_random", - "chrUn_GL000213v1", - "chrUn_KI270338v1", - "chrUn_KI270364v1", - "chrUn_KI270371v1", - "chrUn_KI270372v1", - "chrUn_KI270374v1", - "chrUn_KI270375v1", - "chrUn_KI270424v1", - "chrUn_KI270528v1", - "chrUn_KI270581v1", - "chrUn_KI270587v1", -} - - -class PangenomePipeline(BasePangenome): - """The Pangenome pipeline""" - - params = copy.deepcopy(BasePangenome.params) - params.update( - { - "readgroups": { - "nargs": "*", - "help": ( - "Readgroup information for the fastq files. Only the ID " - "and SM attributes are used." - ), - }, - # Required arguments - "snarls": { - "help": ( - "The vg snarls file for the .gbz file. Can be created " - "using vg." - ), - "required": True, - "type": path_arg(exists=True, is_file=True), - }, - "xg": { - "help": ( - "The xg file for the .gbz file. Can be created using vg" - ), - "required": True, - "type": path_arg(exists=True, is_file=True), - }, - # Additional arguments - "expansion_catalog": { - "help": ( - "An ExpansionHunter variant catalog. Required for " - "expansion calling." - ), - "type": path_arg(exists=True, is_file=True), - }, - "segdup_caller_genes": { - "type": str, - "help": ( - "Genes for SegDup calling. Ex: " - "'CFH,CFHR3,CYP11B1,CYP2D6,GBA,NCF1,PMS2,SMN1,STRC'. " - "Required for SegDup calling." - ), - }, - "skip_cnv": { - "help": "Skip CNV calling.", - "action": "store_true", - }, - "t1k_hla_seq": { - "help": ( - "The DNA HLA seq FASTA file for T1K. Required for HLA " - "calling." - ), - "type": path_arg(exists=True, is_file=True), - }, - "t1k_hla_coord": { - "help": ( - "The DNA HLA coord FASTA file for T1K. Required for HLA " - "calling." - ), - "type": path_arg(exists=True, is_file=True), - }, - "t1k_kir_seq": { - "help": ( - "The DNA KIR seq FASTA file for T1K. Required for KIA " - "calling." - ), - "type": path_arg(exists=True, is_file=True), - }, - "t1k_kir_coord": { - "help": ( - "The DNA KIR coord FASTA file for T1K. Required for KIA " - "calling." - ), - "type": path_arg(exists=True, is_file=True), - }, - } - ) - - positionals = BasePangenome.positionals - - def __init__(self) -> None: - super().__init__() - self.readgroups: List[str] = [] - self.snarls: Optional[pathlib.Path] = None - self.xg: Optional[pathlib.Path] = None - self.expansion_catalog: Optional[pathlib.Path] = None - self.segdup_caller_genes: Optional[str] = None - self.skip_cnv = False - self.t1k_hla_seq: Optional[pathlib.Path] = None - self.t1k_hla_coord: Optional[pathlib.Path] = None - self.t1k_kir_seq: Optional[pathlib.Path] = None - self.t1k_kir_coord: Optional[pathlib.Path] = None - - def main(self, args: argparse.Namespace) -> None: - """Run the pipeline""" - self.handle_arguments(args) - self.setup_logging(args) - self.validate() - - 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) - - self.get_sex(self.ploidy_json) - dag = self.build_second_dag(self.deduped_cram) - executor = self.run(dag) - self.check_execution(dag, executor) - - if not self.retain_tmpdir: - shutil.rmtree(tmp_dir_str) - - def build_dag(self) -> DAG: - return DAG() - - def validate_t1k(self) -> None: - if (self.t1k_hla_seq and not self.t1k_hla_coord) or ( - self.t1k_hla_coord and not self.t1k_hla_seq - ): - self.logger.error( - "For HLA calling, both the seq and coord fasta files need to " - "be supplied. Exiting" - ) - sys.exit(2) - - if (self.t1k_kir_seq and not self.t1k_kir_coord) or ( - self.t1k_kir_coord and not self.t1k_kir_seq - ): - self.logger.error( - "For KIR calling, both the seq and coord fasta files need to " - "be supplied. Exiting" - ) - sys.exit(2) - - if ( - self.t1k_hla_seq or self.t1k_kir_seq - ) and not self.skip_version_check: - for cmd, min_version in T1K_MIN_VERSION.items(): - if not check_version(cmd, min_version): - sys.exit(2) - - def validate_expansion(self) -> None: - if self.expansion_catalog and not self.skip_version_check: - for cmd, min_version in EXPANSION_MIN_VERSION.items(): - if not check_version(cmd, min_version): - sys.exit(2) - - def validate_segdup(self) -> None: - if self.segdup_caller_genes and not self.skip_version_check: - for cmd, min_version in SEGDUP_MIN_VERSION.items(): - if not check_version(cmd, min_version): - sys.exit(2) - - def validate(self) -> None: - """Validate pipeline inputs""" - if not self.output_vcf: - self.logger.error("output_vcf is required") - sys.exit(2) - if not self.reference: - self.logger.error("reference is required") - sys.exit(2) - - self.validate_bundle() - self.validate_fastq_rg(r1_required=True) - self.validate_output_vcf() - self.validate_ref() - self.validate_t1k() - self.validate_expansion() - self.validate_segdup() - - if not self.skip_version_check: - for cmd, min_version in PANGENOME_MIN_VERSIONS.items(): - if not check_version(cmd, min_version): - sys.exit(2) - - def validate_bundle(self) -> None: - 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()) - try: - req_version = packaging.version.Version( - bundle_info["minScriptVersion"] - ) - bundle_pipeline = bundle_info["pipeline"] - except KeyError: - self.logger.error( - "The model bundle does not have the expected attributes" - ) - sys.exit(2) - 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) - if bundle_pipeline != "DNAscope pangenome": - self.logger.error("The model bundle is for a different pipeline.") - sys.exit(2) - - bundle_members = set(ar_load(str(self.model_bundle))) - if ( - "bwa.model" not in bundle_members - or "cnv.model" not in bundle_members - or "dnascope.model" not in bundle_members - ): - self.logger.error( - "Expected model files not found in the model bundle file" - ) - sys.exit(2) - - def validate_fastq_rg(self, r1_required=False) -> None: - if not self.r1_fastq and r1_required: - self.logger.error("Please supply --r1_fastq arguments") - sys.exit(2) - - if len(self.r1_fastq) != len(self.readgroups): - self.logger.error( - "The number of readgroups does not equal the number of fastq " - "files" - ) - sys.exit(2) - - # Validate readgroups - rg_sample = None - for rg in self.readgroups: - rg_dict = parse_rg_line(rg.replace(r"\t", "\t")) - rg_sm = rg_dict.get("SM") - if not rg_sm: - self.logger.error( - "Found a readgroup without a SM tag: %s", - str(rg), - ) - sys.exit(2) - if rg_sample and rg_sample != rg_sm: - self.logger.error( - "Inconsistent readgroup sample information found in: %s", - str(rg), - ) - sys.exit(2) - rg_sample = rg_sm - if "ID" not in rg_dict: - self.logger.error( - "Found a readgroup without an ID tag: %s", - str(rg), - ) - sys.exit(2) - - def configure(self) -> None: - """Configure pipeline parameters""" - pass - - def build_first_dag(self) -> DAG: - """Build the first DAG for the pangenome pipeline""" - self.logger.info("Building the first pangenome DAG") - dag = DAG() - - # Find the sample name - sample_name = parse_rg_line(self.readgroups[0].replace(r"\t", "\t"))[ - "SM" - ] - - # Output files - suffix = "bam" if self.bam_format else "cram" - out_svs = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_svs.vcf.gz") - ) - self.deduped_cram = pathlib.Path( - str(self.output_vcf).replace( - ".vcf.gz", f"_pangenome-aligned.{suffix}" - ) - ) - self.ploidy_json = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_ploidy.json") - ) - out_hla = pathlib.Path(str(self.output_vcf).replace(".vcf.gz", "_hla")) - out_kir = pathlib.Path(str(self.output_vcf).replace(".vcf.gz", "_kir")) - out_segdup = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_segdups") - ) - - # Intermediate file paths - surject_paths_dict = self.tmp_dir.joinpath("GRCh38_pansn_paths.dict") - sv_header = self.tmp_dir.joinpath("sample_sv.hdr") - kmer_prefix = self.tmp_dir.joinpath("sample.fq") - kmer_file = pathlib.Path(str(kmer_prefix) + ".kff") - sample_pangenome = self.tmp_dir.joinpath("sample_pangenome.gbz") - sample_pack = self.tmp_dir.joinpath("sample_coverage.pack") - lc_score = self.tmp_dir.joinpath("sample_lc_score.txt.gz") - # - Ensure we have a bam file for t1k - rw_bam = self.tmp_dir.joinpath("sample_deduped.bam") - dnascope_tmp = self.tmp_dir.joinpath("sample_dnascope_tmp.vcf.gz") - - # Create files used to generate headers - self.create_paths_dict(surject_paths_dict) - self.create_sv_header(sv_header, sample_name) - - # KMC k-mer counting - kmc_job = self.build_kmc_job(kmer_prefix, self.cores) - dag.add_job(kmc_job) - - # vg haplotypes - create sample-specific pangenome - haplotypes_job = self.build_haplotypes_job(sample_pangenome, kmer_file) - dag.add_job(haplotypes_job, {kmc_job}) - - # vg giraffe - align reads to pangenome - sample_gams, giraffe_jobs = self.build_alignment_jobs(sample_pangenome) - alignment_dependencies = {haplotypes_job} - for job in giraffe_jobs: - dag.add_job(job, alignment_dependencies) - - # vg pack - compute read support - pack_job = self.build_pack_job(sample_pack, sample_gams) - dag.add_job(pack_job, set(giraffe_jobs)) - - # vg call - call SVs - call_job = self.build_call_job( - out_svs, sample_pack, sv_header, sample_name - ) - dag.add_job(call_job, {pack_job}) - - # surject and sort - move reads back to the linear reference - aligned_bams, surject_jobs = self.build_surject_jobs( - sample_gams, surject_paths_dict - ) - for job in surject_jobs: - dag.add_job(job, set(giraffe_jobs)) - - # duplicate marking - hdr_jobs, lc_job, dedup_job = self.build_dedup_job( - lc_score, self.deduped_cram, aligned_bams - ) - for job in hdr_jobs: - dag.add_job(job, set(surject_jobs)) - dag.add_job(lc_job, set(hdr_jobs)) - dag.add_job(dedup_job, {lc_job}) - - # metrics - collect some metrics after dedup - metrics_job, rehead_job = self.build_metrics_job( - rw_bam, self.deduped_cram - ) - dag.add_job(metrics_job, {dedup_job}) - dag.add_job(rehead_job, {metrics_job}) - - # multiqc - multiqc_job = self.multiqc() - if multiqc_job: - dag.add_job(multiqc_job, {rehead_job}) - - # small variant calling - dnascope_job, dnamodelapply_job = self.build_smallvar_job( - dnascope_tmp, self.deduped_cram - ) - dag.add_job(dnascope_job, {dedup_job}) - dag.add_job(dnamodelapply_job, {dnascope_job}) - - # estimate ploidy - ploidy_job = self.build_ploidy_job( - self.ploidy_json, - [rw_bam], - ) - dag.add_job(ploidy_job, {metrics_job}) - - # HLA and KIR calling with T1K - if self.t1k_hla_seq and self.t1k_hla_coord: - hla_job = self.build_t1k_job( - out_hla, - rw_bam, - self.t1k_hla_seq, - self.t1k_hla_coord, - "hla-wgs", - ) - dag.add_job(hla_job, {metrics_job}) - if self.t1k_kir_seq and self.t1k_kir_coord: - kir_job = self.build_t1k_job( - out_kir, - rw_bam, - self.t1k_kir_seq, - self.t1k_kir_coord, - "kir-wgs", - ) - dag.add_job(kir_job, {metrics_job}) - - # special-caller - if self.segdup_caller_genes: - segdup_job = self.build_segdup_job( - out_segdup, self.deduped_cram, genes=self.segdup_caller_genes - ) - dag.add_job(segdup_job, {dedup_job}) - - return dag - - def create_paths_dict( - self, - surject_paths_dict: pathlib.Path, - ctg_prefix="GRCh38#0#", - assembly="38", - species="Homo sapiens", - ): - """Create a .dict file for surjection""" - contigs: List[Tuple[str, str]] = [] - fai_file = str(self.reference) + ".fai" - if not self.dry_run: - with open(fai_file) as fh: - for line in fh: - line_split = line.rstrip().split("\t") - contigs.append((line_split[0], line_split[1])) - - with open(surject_paths_dict, "w") as ofh: - print("@HD\tVN:1.5", file=ofh) - for ctg_name, ctg_len in contigs: - if ctg_name in GRCh38_CONTIGS_NOT_IN_PANGENOME: - continue - flds = { - "SN": f"{ctg_prefix}{ctg_name}", - "LN": ctg_len, - "AS": assembly, - "SP": species, - } - print( - "@SQ\t" - + "\t".join([f"{x}:{y}" for x, y in flds.items()]), - file=ofh, - ) - - def create_sv_header(self, sv_header: pathlib.Path, sample_name: str): - """Create a header file for the SV VCF""" - contigs: List[Tuple[str, str]] = [] - fai_file = str(self.reference) + ".fai" - if not self.dry_run: - with open(fai_file) as fh: - for line in fh: - line_split = line.rstrip().split("\t") - contigs.append((line_split[0], line_split[1])) - - if not self.dry_run: - with open(sv_header, "w") as fh: - print(r"##fileformat=VCFv4.2", file=fh) - for contig, contig_l in contigs: - print( - f"##contig=", - file=fh, - ) - - for fld in ("FILTER", "INFO", "FORMAT"): - for hdr_id, hdr_d in SV_HDR_ATTR[fld].items(): - hdr_kv = [f"ID={hdr_id}"] - for k, v in hdr_d.items(): - hdr_kv.append(f"{k}={v}") - hdr_kv_str = ",".join(hdr_kv) - print(f"##{fld}=<{hdr_kv_str}>", file=fh) - - last_hdr_line = [ - "#CHROM", - "POS", - "ID", - "REF", - "ALT", - "QUAL", - "FILTER", - "INFO", - "FORMAT", - sample_name, - ] - print("\t".join(last_hdr_line), file=fh) - - def build_haplotypes_job( - self, output_gbz: pathlib.Path, kmer_file: pathlib.Path - ) -> Job: - """Build vg haplotypes job""" - if not self.hapl: - self.logger.error("hapl is required") - sys.exit(2) - if not self.gbz: - self.logger.error("gbz is required") - sys.exit(2) - - haplotypes_job = Job( - cmds.cmd_vg_haplotypes( - output_gbz, - kmer_file, - self.hapl, - self.gbz, - threads=self.cores, - ), - "vg-haplotypes", - self.cores, - ) - return haplotypes_job - - def build_alignment_jobs( - self, sample_pangenome: pathlib.Path - ) -> Tuple[List[pathlib.Path], List[Job]]: - """Build vg giraffe alignment jobs""" - # Process each pair of fastq files - sample_gams: List[pathlib.Path] = [] - giraffe_jobs: List[Job] = [] - for i, (r1, r2, rg) in enumerate( - itertools.zip_longest( - self.r1_fastq, self.r2_fastq, self.readgroups - ) - ): - sample_gam = self.tmp_dir.joinpath(f"sample_aligned_{i}.gam") - rg_dict = parse_rg_line(rg.replace(r"\t", "\t")) - giraffe_job = Job( - cmds.cmd_vg_giraffe( - sample_gam, - sample_pangenome, - fastq1=r1, - fastq2=r2, - readgroup=rg_dict["ID"], - sample=rg_dict["SM"], - threads=self.cores, - ), - "vg-giraffe", - self.cores, - ) - - sample_gams.append(sample_gam) - giraffe_jobs.append(giraffe_job) - return (sample_gams, giraffe_jobs) - - def build_pack_job( - self, - sample_pack: pathlib.Path, - sample_gams: List[pathlib.Path], - ) -> Job: - """Compute read support with vg pack""" - if not self.gbz: - self.logger.error("gbz is required") - sys.exit(2) - - pack_job = Job( - cmds.cmd_vg_pack( - sample_pack, - sample_gams, - self.gbz, - threads=self.cores, - ), - "vg-pack", - self.cores, - ) - return pack_job - - def build_call_job( - self, - out_svs: pathlib.Path, - sample_pack: pathlib.Path, - sv_header: pathlib.Path, - sample_name="", - ) -> Job: - """Call SVs with vg call""" - if not self.gbz: - self.logger.error("gbz is required") - sys.exit(2) - if not self.snarls: - self.logger.error("snarls is required") - sys.exit(2) - - call_job = Job( - cmds.cmd_sv_call( - out_svs, - sample_pack, - sv_header, - self.gbz, - self.snarls, - sample_name=sample_name, - ), - "vg-call", - self.cores, - ) - return call_job - - def build_surject_jobs( - self, - sample_gams: List[pathlib.Path], - surject_paths_dict: pathlib.Path, - ) -> Tuple[List[pathlib.Path], List[Job]]: - """Surject gam files on the linear reference""" - if not self.xg: - self.logger.error("xg is required") - sys.exit(2) - - aligned_bams: List[pathlib.Path] = [] - surject_jobs: List[Job] = [] - - for i, sample_gam in enumerate(sample_gams): - sample_bam = self.tmp_dir.joinpath(f"sample_aligned_{i}.bam") - surject_job = Job( - cmds.cmd_vg_surject( - sample_bam, - sample_gam, - self.xg, - surject_paths_dict, - threads=self.cores, - ), - f"vg_surject_{i}", - self.cores, - ) - - aligned_bams.append(sample_bam) - surject_jobs.append(surject_job) - return (aligned_bams, surject_jobs) - - def build_dedup_job( - self, - lc_score: pathlib.Path, - deduped_bam: pathlib.Path, - aligned_bams: List[pathlib.Path], - ) -> Tuple[List[Job], Job, Job]: - """Mark duplicates""" - if not self.output_vcf: - self.logger.error("output_vcf is required") - sys.exit(2) - - # Create .hdr files for each input file - hdr_jobs: List[Job] = [] - for i, bam in enumerate(aligned_bams): - hdr = pathlib.Path(str(bam) + ".hdr") - hdr_jobs.append( - Job( - cmds.strip_ctg_prefix( - hdr, - bam, - "GRCh38#0#", - ), - f"create_hdr_{i}", - 0, - ) - ) - - # Output metrics - metric_base = self.output_vcf.name.replace(".vcf.gz", "") + ".txt" - metrics_dir = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_metrics") - ) - try: - os.mkdir(metrics_dir) - except FileExistsError: - pass - 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") - gc_metrics = metrics_dir.joinpath(metric_base + ".gc_bias.txt") - gc_summary = metrics_dir.joinpath(metric_base + ".gc_bias_summary.txt") - dedup_metrics = metrics_dir.joinpath( - metric_base + ".dedup_metrics.txt" - ) - - driver = Driver( - reference=self.reference, - thread_count=self.cores, - input=aligned_bams, - ) - driver.add_algo(LocusCollector(lc_score)) - 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)) - - lc_job = Job( - Pipeline(Command(*driver.build_cmd())), - "locuscollector", - self.cores, - ) - - # Dedup - driver = Driver( - reference=self.reference, - thread_count=self.cores, - input=aligned_bams, - read_filter="IndelLeftAlignReadTransform", - ) - driver.add_algo( - Dedup( - deduped_bam, - lc_score, - cram_write_options="version=3.0,compressor=rans", - metrics=dedup_metrics, - ) - ) - dedup_job = Job( - Pipeline(Command(*driver.build_cmd())), "dedup", self.cores - ) - - return (hdr_jobs, lc_job, dedup_job) - - def build_metrics_job( - self, - deduped_bam: pathlib.Path, - deduped_cram: pathlib.Path, - ) -> Tuple[Job, Job]: - if not self.output_vcf: - self.logger.error("output_vcf is required") - sys.exit(2) - - # Output metrics - metric_base = self.output_vcf.name.replace(".vcf.gz", "") + ".txt" - metrics_dir = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_metrics") - ) - - # Output metrics - # - Run some metrics after duplicate marking - is_metrics = metrics_dir.joinpath(metric_base + ".insert_size.txt") - wgs_metrics = metrics_dir.joinpath(metric_base + ".wgs.txt") - - driver = Driver( - reference=self.reference, - thread_count=self.cores, - input=[deduped_cram], - ) - driver.add_algo(ReadWriter(deduped_bam)) - driver.add_algo(InsertSizeMetricAlgo(is_metrics)) - driver.add_algo(WgsMetricsAlgo(wgs_metrics, include_unpaired="true")) - metrics_job = Job( - Pipeline(Command(*driver.build_cmd())), - "metrics", - self.cores, - ) - - # Rehead WGS metrics so they are recognized by MultiQC - rehead_script = pathlib.Path( - str( - files("sentieon_cli.scripts").joinpath("rehead_wgs_metrics.py") - ) - ) - rehead_job = Job( - Pipeline( - Command( - "sentieon", - "pyexec", - str(rehead_script), - "--metrics_file", - str(wgs_metrics), - ) - ), - "Rehead metrics", - 0, - ) - return (metrics_job, rehead_job) - - def build_smallvar_job( - self, - dnascope_tmp: pathlib.Path, - realigned_cram: pathlib.Path, - ) -> Tuple[Job, Job]: - """Small variant calling""" - 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) - - model = self.model_bundle.joinpath("dnascope.model") - interval_ctgs = list(["chr" + str(x) for x in range(1, 23)]) - interval_ctgs.extend(["chrX", "chrY"]) - interval = ",".join(interval_ctgs) - - pcr_indel_model = "NONE" if self.pcr_free else "CONSERVATIVE" - driver = Driver( - reference=self.reference, - thread_count=self.cores, - input=[realigned_cram], - interval=interval, - ) - driver.add_algo( - DNAscope( - dnascope_tmp, - dbsnp=self.dbsnp, - pcr_indel_model=pcr_indel_model, - model=model, - ) - ) - dnascope_job = Job( - Pipeline(Command(*driver.build_cmd())), "dnascope", self.cores - ) - - driver = Driver( - reference=self.reference, - thread_count=self.cores, - ) - driver.add_algo( - DNAModelApply( - model, - dnascope_tmp, - self.output_vcf, - ) - ) - dnamodelapply_job = Job( - Pipeline(Command(*driver.build_cmd())), "model-apply", self.cores - ) - return (dnascope_job, dnamodelapply_job) - - def build_second_dag(self, input_aln: pathlib.Path) -> DAG: - """Build the second DAG for the pangenome pipeline""" - self.logger.info("Building the second pangenome DAG") - dag = DAG() - - # Output files - out_cnv = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_cnv.vcf.gz") - ) - out_expansions = pathlib.Path( - str(self.output_vcf).replace(".vcf.gz", "_expansion") - ) - - # Intermediate file paths - cnvscope_tmp = self.tmp_dir.joinpath("sample_cnvcope_tmp.vcf.gz") - cnvscope_autosomes = self.tmp_dir.joinpath( - "sample_cnvcope_autosomes.vcf.gz" - ) - cnvscope_haploid_tmp = self.tmp_dir.joinpath( - "sample_cnvcope_haploid_tmp.vcf.gz" - ) - cnvscope_haploid = self.tmp_dir.joinpath( - "sample_cnvcope_haploid.vcf.gz" - ) - - # CNV calling - if not self.skip_cnv: - ( - cnvscope_job, - cnvapply_job, - haploid_cnv_job, - haploid_apply_job, - concat_job, - ) = self.build_cnv_jobs( - out_cnv, - cnvscope_tmp, - cnvscope_autosomes, - cnvscope_haploid_tmp, - cnvscope_haploid, - input_aln, - ) - dag.add_job(cnvscope_job) - dag.add_job(cnvapply_job, {cnvscope_job}) - if haploid_cnv_job and haploid_apply_job and concat_job: - dag.add_job(haploid_cnv_job) - dag.add_job(haploid_apply_job, {haploid_cnv_job}) - dag.add_job(concat_job, {cnvapply_job, haploid_apply_job}) - - # Repeat expansions - if self.expansion_catalog: - expansion_job = self.build_expansion_job( - out_expansions, - input_aln, - self.expansion_catalog, - ) - dag.add_job(expansion_job) - - return dag - - def build_cnv_jobs( - self, - out_cnv: pathlib.Path, - cnvscope_tmp: pathlib.Path, - cnvscope_autosomes: pathlib.Path, - cnvscope_haploid_tmp: pathlib.Path, - cnvscope_haploid: pathlib.Path, - realigned_cram: pathlib.Path, - ) -> Tuple[Job, Job, Optional[Job], Optional[Job], Optional[Job]]: - """Call CNVs with CNVscope""" - if not self.model_bundle: - self.logger.error("model_bundle is required") - sys.exit(2) - - model = self.model_bundle.joinpath("cnv.model") - - # First pass will be autosomes for males - # or autosomes + chrX for females - first_contigs = list(["chr" + str(x) for x in range(1, 23)]) - if self.sample_sex != SampleSex.MALE: - first_contigs.append("chrX") - - driver = Driver( - reference=self.reference, - thread_count=self.cores, - input=[realigned_cram], - interval=",".join(first_contigs), - ) - driver.add_algo( - CNVscope( - cnvscope_tmp, - model=model, - ) - ) - cnvscope_job = Job( - Pipeline(Command(*driver.build_cmd())), - "cnvscope-autosomes", - self.cores, - ) - - # Autosome model apply - outfile = ( - cnvscope_autosomes - if self.sample_sex == SampleSex.MALE - else out_cnv - ) - driver = Driver( - reference=self.reference, - thread_count=self.cores, - ) - driver.add_algo( - CNVModelApply( - outfile, - model, - cnvscope_tmp, - ) - ) - cnvapply_job = Job( - Pipeline(Command(*driver.build_cmd())), - "cnvmodelapply-autosomes", - self.cores, - ) - - if self.sample_sex != SampleSex.MALE: - return (cnvscope_job, cnvapply_job, None, None, None) - - # Sex chrom CNVscope - driver = Driver( - reference=self.reference, - thread_count=self.cores, - input=[realigned_cram], - interval="chrX,chrY", - ) - driver.add_algo( - CNVscope( - cnvscope_haploid_tmp, - model=model, - ) - ) - haploid_cnv_job = Job( - Pipeline(Command(*driver.build_cmd())), - "cnvscope-haploid", - self.cores, - ) - - # Sex chrom CNVModelApply - driver = Driver( - reference=self.reference, - thread_count=self.cores, - ) - driver.add_algo( - CNVModelApply( - cnvscope_haploid, - model, - cnvscope_haploid_tmp, - ) - ) - haploid_apply_job = Job( - Pipeline(Command(*driver.build_cmd())), - "cnvmodelapply-haploid", - self.cores, - ) - - # Concate autosome and sex chromsome VCFs - concat_job = Job( - cmds.bcftools_concat( - out_cnv, [cnvscope_autosomes, cnvscope_haploid] - ), - "concat-calls", - 0, - ) - - return ( - cnvscope_job, - cnvapply_job, - haploid_cnv_job, - haploid_apply_job, - concat_job, - ) - - def build_expansion_job( - self, - out_expansions: pathlib.Path, - realigned_cram: pathlib.Path, - expansion_catalog: pathlib.Path, - ) -> Job: - """Identify repeat expansions""" - if not self.reference: - self.logger.error("reference is required") - sys.exit(2) - - expansion_job = Job( - cmds.cmd_expansion_hunter( - out_expansions, - realigned_cram, - reference=self.reference, - variant_catalog=expansion_catalog, - sex="male" if self.sample_sex == SampleSex.MALE else "female", - threads=self.cores, - ), - "expansion-hunter", - 0, - ) - return expansion_job - - def build_t1k_job( - self, - out_basename: pathlib.Path, - deduped_bam: pathlib.Path, - gene_seq: pathlib.Path, - gene_coord: pathlib.Path, - preset: str, - ) -> Job: - """Genotype diverse sequences with T1K""" - job = Job( - cmds.cmd_t1k( - out_basename, - deduped_bam, - gene_seq=gene_seq, - gene_coord=gene_coord, - preset=preset, - threads=1, - ), - "T1K-HLA-KIR", - 0, - ) - return job - - def build_segdup_job( - self, - out_segdup: pathlib.Path, - realigned_cram: pathlib.Path, - genes: str, - ) -> Job: - """Call variants in difficult SegDups""" - 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) - - segdup_job = Job( - cmds.cmd_segdup_caller( - out_segdup, - realigned_cram, - reference=self.reference, - sr_bundle=self.model_bundle, - genes=genes, - ), - "segdup-caller", - self.cores, - ) - return segdup_job