diff --git a/Dockerfile b/Dockerfile index 800b254..c17c07d 100644 --- a/Dockerfile +++ b/Dockerfile @@ -26,9 +26,9 @@ ARG EXPANSIONHUNTER_SHA256=ebf3ec0ace6e6e3bbce12c26463da5d9f8e16374eff1ad10f0f1a ARG T1K_VERSION=1.0.9 ARG T1K_GIT_TAG=v1.0.9 ARG T1K_COMMIT=9376b555c1d8d2f8ca357c2656f49f450462dbc3 -ARG SEGDUP_CALLER_VERSION=0.5.1 -ARG SEGDUP_CALLER_GIT_TAG=v0.5.1 -ARG SEGDUP_CALLER_COMMIT=0406ea78b7ff53d7a169b9aa945d7fcc7257ff12 +ARG SEGDUP_CALLER_VERSION=0.9.0 +ARG SEGDUP_CALLER_GIT_TAG=v0.9.0 +ARG SEGDUP_CALLER_COMMIT=37c299fe16970215712b7a4ebdf1932f923e903a # Pinned Poetry toolchain; should match the version used to generate poetry.lock. ARG POETRY_VERSION=2.3.4 ARG POETRY_PLUGIN_EXPORT_VERSION=1.9.0 @@ -291,8 +291,8 @@ RUN sentieon driver --help && \ kmc --help && \ ExpansionHunter --help && \ perl -c "$(command -v run-t1k)" && \ - segdup-caller --version && \ - python -c "import sys; from packaging.version import Version; \ + SENTIEON_LICENSE=placeholder segdup-caller --version && \ + SENTIEON_LICENSE=placeholder python -c "import sys; from packaging.version import Version; \ from sentieon_cli.sentieon_pangenome import SEGDUP_MIN_VERSION; \ from sentieon_cli.util import check_version; \ sys.exit(0 if all(check_version(c, v) for c, v in SEGDUP_MIN_VERSION.items()) else 1)" && \ diff --git a/docs/dnascope-hybrid.md b/docs/dnascope-hybrid.md deleted file mode 100644 index 86442bb..0000000 --- a/docs/dnascope-hybrid.md +++ /dev/null @@ -1,135 +0,0 @@ -# DNAscope Hybrid - -Sentieon DNAscope Hybrid is a pipeline for germline variant calling using combined short-read and long-read data from a single sample. The DNAscope Hybrid pipeline is able to utilize the strengths of both short and long-read technologies to generate variant callsets that are more accurate than either short-read or long-read data alone. - -The pipeline supports input data in the following formats; both short-read and long-read input are required: -* Unaligned short-read data in gzipped FASTQ format. -* Aligned short-reads in BAM or CRAM format. -* Unaligned long-read data in the uBAM or uCRAM format. -* Aligned long-read data in BAM or CRAM format. - -By default, the pipeline will generate the following output files: -* Small variants (SNVs and indels) in the VCF format. -* Structural variants in the VCF format. -* Copy-number variants in the VCF format. - -If unaligned reads are used as input, the pipeline will also output aligned reads in BAM or CRAM format. - -The DNAscope Hybrid pipeline is implemented using 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.01 or higher. -- [Python] version 3.8 or higher. -- [bcftools] version 1.10 or higher. -- [bedtools] -- [MultiQC] version 1.18 or higher for metrics report generation. -- [samtools] version 1.16 or higher. -- [mosdepth] version 0.2.6 or higher for coverage metrics collection from long-read data. - -The `sentieon`, `python`, `bcftools`, `bedtools`, `samtools`, `multiqc`, and `mosdepth` executables will be accessed through the user's `PATH` environment variable. - -## Input data requirements - -### The Reference genome - -DNAscope LongRead will call variants present in the sample relative to a high quality reference genome in FASTA format. Besides the reference genome file, a samtools fasta index file (.fai) needs to be present. Short-read alignment also requires bwa index files. - -We recommend aligning to a reference genome without alternate contigs. If alternate contigs are present in the genome and the pipeline is performing short-read alignment, please also supply a ".alt" file to activate [alt-aware alignment] in bwa. - -## Usage - -### Germline variant calling from aligned short and long-read data - -A single command is run to call SNVs, indels, SVs and CNVs from aligned short and long reads: -```sh -sentieon-cli dnascope-hybrid \ - -r REFERENCE \ - --sr_aln SR_ALN [SR_ALN ...] \ - --lr_aln LR_ALN [LR_ALN ...] \ - -m MODEL_BUNDLE \ - [-b DIPLOID_BED] \ - [-d DBSNP] \ - [--dry_run] \ - [--gvcf] \ - [--sr_duplicate_marking MARKDUP] \ - [-t NUMBER_THREADS] \ - sample.vcf.gz -``` - -The DNAscope Hybrid pipeline requires the following arguments: -- `-r REFERENCE`: the location of the reference FASTA file. A reference fasta index, ".fai" file, is also required. -- `--sr_aln`: the input short-read data in BAM or CRAM format. One or more files can be supplied by passing multiple files after the `--sr_aln` argument. -- `--lr_aln`: the input long-read data in BAM or CRAM format. One or more files can be supplied by passing multiple files after the `--lr_aln` argument. -- `-m MODEL_BUNDLE`: the location of the model bundle. Model bundle files can be found in the [sentieon-models] repository. -- `sample.vcf.gz`: the location of the output VCF file for SNVs and indels. The pipeline requires the output file end with the suffix, ".vcf.gz". - -The DNAscope Hybrid pipeline accepts the following optional arguments: -- `-b DIPLOID_BED`: interval in the reference to restrict diploid variant calling, in BED file format. Supplying this file will limit diploid variant calling to the intervals inside the BED file. -- `-d DBSNP`: the location of the Single Nucleotide Polymorphism database (dbSNP) used to label known variants in VCF (`.vcf`) or bgzip compressed VCF (`.vcf.gz`) format. Only one file is supported. Supplying this file will annotate variants with their dbSNP refSNP ID numbers. A VCF index file is required. -- `--dry_run`: print the pipeline commands, but do not actually execute them. -- `--gvcf`: generate a gVCF output file along with the VCF. -- `--sr_duplicate_marking`: setting for duplicate marking. `markdup` will mark duplicate reads. `rmdup` will remove duplicate reads. `none` will skip duplicate marking. The default setting is `markdup`. -- `-t NUMBER_THREADS`: number of computing threads that will be used by the software to run parallel processes. The argument is optional; if omitted, the pipeline will use as many threads as the server has. -- `-h`: print the command-line help and exit. - -### Germline variant calling from unaligned short and long-read data - -A single command is run to call SNVs, indels, SVs and CNVs from unaligned short and long reads: -```sh -sentieon-cli dnascope-hybrid \ - -r REFERENCE \ - --sr_r1_fastq SR_R1_FQ [SR_R1_FQ ...] \ - --sr_r2_fastq SR_R2_FQ [SR_R2_FQ ...] \ - --sr_readgroups SR_READGROUP [SR_READGROUP ...] \ - --lr_aln LR_ALN [LR_ALN ...] \ - --lr_align_input \ - -m MODEL_BUNDLE \ - [-b DIPLOID_BED] \ - [--bam_format] \ - [-d DBSNP] \ - [--dry_run] \ - [--gvcf] \ - [--sr_duplicate_marking MARKDUP] \ - [-t NUMBER_THREADS] \ - sample.vcf.gz -``` - -The DNAscope Hybrid pipeline requires the following arguments: -- `--sr_r1_fastq`: the input R1 short-read data in gzipped FASTQ format. One or more files can be supplied by passing multiple files after the `--sr_r1_fastq` argument. -- `--sr_r2_fastq`: the input R2 short-read data in gzipped FASTQ format. One or more files can be supplied by passing multiple files after the `--sr_r2_fastq` argument. -- `--sr_readgroups`: readgroup information for each FASTQ. The pipeline will expect the same number of arguments to `--sr_r1_fastq` and `--sr_readgroups`. An example argument is, `--sr_readgroups "@RG\tID:HG002-1\tSM:HG002\tLB:HG002-LB-1\tPL:ILLUMINA"` -- `--lr_aln`: the input long-read data in uBAM or uCRAM format. One or more files can be supplied by passing multiple files after the `--lr_aln` argument. -- `--lr_align_input`: directs the pipeline to align the input long-reads. - -The DNAscope Hybrid pipeline accepts the following optional arguments: -- `--bam_format`: use BAM format instead of CRAM for output aligned files. -- `--lr_input_ref`: a reference fasta used for decoding the input long-read file(s). Required with long-read uCRAM or CRAM input. Can be different from the fasta used with the `-r` argument. - -## Pipeline output - -### List of output files - -The following files are output by the DNAscope Hybrid pipeline: -- `sample.vcf.gz`: SNV and indel variant calls across the regions of the genome as defined in the `-b DIPLOID_BED` file. -- `sample.sv.vcf.gz`: structural variant calls from the Sentieon LongReadSV tool. -- `sample.cnv.vcf.gz`: copy-number variant calls from the Sentieon CNVscope tool. -- `sample_deduped.cram`: aligned, coordinate-sorted and duplicate-marked short-read data from the input FASTQ files. -- `sample_mm2_sorted_*.cram`: aligned and coordinate-sorted long-reads from the input uBAM, uCRAM, BAM, or CRAM files. -- `sample_metrics`: a directory containing QC metrics for the analyzed sample. - -## Troubleshooting - -### The pipeline complains, "Input ... has a different RG-SM tag" - -This error will occur if the pipeline detects that the input files have (or will have) different readgroup `SM` tags. To fix this error, please use the `--rgsm` argument to adjust the `SM` tags of the input files during variant calling. Note that with this argument, all reads in the input files will be used during variant calling. - - -[Python]: https://www.python.org/ -[bcftools]: http://samtools.github.io/bcftools/bcftools.html -[bedtools]: https://bedtools.readthedocs.io/en/latest/ -[MultiQC]: https://multiqc.info/ -[mosdepth]: https://github.com/brentp/mosdepth -[samtools]: https://www.htslib.org/ -[alt-aware alignment]: https://github.com/lh3/bwa/blob/master/README-alt.md -[sentieon-models]: https://github.com/Sentieon/sentieon-models diff --git a/docs/dnascope-longread.md b/docs/dnascope-longread.md deleted file mode 100644 index 9ee2f4b..0000000 --- a/docs/dnascope-longread.md +++ /dev/null @@ -1,153 +0,0 @@ -# DNAscope LongRead - -Sentieon DNAscope LongRead is a pipeline for alignment and germline variant calling (SNVs, SVs, CNVs, and indels) from long-read sequence data. The DNAscope LongRead pipeline is able to take advantage of longer read lengths to perform quick and accurate variant calling using specially calibrated machine learning models. - -The pipeline will accept as input aligned reads in BAM or CRAM format, or un-aligned reads in FASTQ, uBAM, or uCRAM format. The pipeline will output variants in the VCF (or gVCF) formats and aligned reads in BAM or CRAM formats. - -DNAscope LongRead is implemented using 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 202308.01 or higher. -- [Python] version 3.8 or higher. -- [bcftools] version 1.10 or higher for the SNV and indel calling pipeline. -- [bedtools] for the SNV and indel calling pipeline. -- [samtools] version 1.16 or higher for alignment of read data in uBAM or uCRAM format or re-alignment of previously aligned reads. -- [mosdepth] version 0.2.6 or higher for coverage metrics of long-read data. -- [hificnv] version 1.0.0 or higher for CNV calling. - -The `sentieon`, `python`, `bcftools`, `bedtools`, `samtools`, `hificnv`, and `mosdepth` executables will be accessed through the user's `PATH` environment variable. - -## Input data requirements - -### The Reference genome - -DNAscope LongRead 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. - -We recommend aligning to a reference genome without alternate contigs. - - -## Usage - -### Alignment and variant calling from FASTQ - -A single command is run to align and call SNVs, indels, and structural variants from PacBio HiFi or ONT reads in the FASTQ format: -```sh -sentieon-cli dnascope-longread [-h] \ - -r REFERENCE \ - --fastq INPUT_FASTQ ... \ - --readgroups READGROUP ... \ - -m MODEL_BUNDLE \ - [-d DBSNP] \ - [-b DIPLOID_BED] \ - [-t NUMBER_THREADS] \ - [-g] \ - --tech HiFi|ONT \ - [--haploid_bed HAPLOID_BED] \ - [--cnv_excluded_regions CNV_EXCLUDE_BED] \ - sample.vcf.gz -``` - -With FASTQ input, the DNAscope LongRead pipeline requires the following arguments: -- `-r REFERENCE`: the location of the reference FASTA file. A reference fasta index, ".fai" file, is also required. -- `--fastq INPUT_FASTQ`: the input sample file in FASTQ format. One or more files can be supplied by passing multiple files after the `--fastq` argument. -- `--readgroups READGROUP`: readgroup information for the read data. The `--readgroups` argument is required if the input data is in the FASTQ format. This argument expects complete readgroup strings and these strings will be passed to `minimap2` through the `-R` argument. An example argument is `--readgroups '@RG\tID:foo\tSM:bar'`. -- `-m MODEL_BUNDLE`: the location of the model bundle. Model bundle files can be found in the [sentieon-models] repository. -- `--tech HiFi|ONT`: Sequencing technology used to generate the reads. Supported arguments are `ONT` or `HiFi`. -- `sample.vcf.gz`: the location of the output VCF file for SNVs and indels. 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 Sentieon LongRead pipeline accepts the following optional arguments: -- `-d DBSNP`: the location of the Single Nucleotide Polymorphism database (dbSNP) used to label known variants in VCF (`.vcf`) or bgzip compressed VCF (`.vcf.gz`) format. Only one file is supported. Supplying this file will annotate variants with their dbSNP refSNP ID numbers. A VCF index file is required. -- `-b DIPLOID_BED`: interval in the reference to restrict diploid variant calling, in BED file format. Supplying this file will limit diploid variant calling to the intervals inside the BED file. -- `--haploid_bed HAPLOID_BED`: interval in the reference to restrict haploid variant calling, in BED file format. Supplying this file will perform haploid variant calling across the intervals inside the BED file. -- `--cnv_excluded_regions`: a BED file of excluded CNV regions passed to hificnv. See the hificnv documentation for more details, https://github.com/PacificBiosciences/HiFiCNV. -- `-t NUMBER_THREADS`: number of computing threads that will be used by the software to run parallel processes. The argument is optional; if omitted, the pipeline will use as many threads as the server has. -- `-g`: output variants in the gVCF format, in addition to the VCF output file. The tool will output a bgzip compressed gVCF file with a corresponding index file. -- `-h`: print the command-line help and exit. -- `--dry_run`: print the pipeline commands, but do not actually execute them. - -### Alignment and variant calling from uBAM, uCRAM, BAM, or CRAM - -A single command is run to align and call SNVs, indels, and structural variants from PacBio HiFi or ONT reads in the uBAM, uCRAM, BAM, or CRAM formats: -```sh -sentieon-cli dnascope-longread [-h] \ - -r REFERENCE \ - -i SAMPLE_INPUT ... \ - --align \ - -m MODEL_BUNDLE \ - [-d DBSNP] \ - [-b DIPLOID_BED] \ - [-t NUMBER_THREADS] \ - [-g] \ - --tech HiFi|ONT \ - [--haploid_bed HAPLOID_BED] \ - [--cnv_excluded_regions CNV_EXCLUDE_BED] \ - [--input_ref INPUT_REF] \ - sample.vcf.gz -``` - -With uBAM, uCRAM, BAM, or CRAM input, the DNAscope LongRead pipeline requires the following new arguments: -- `-i SAMPLE_INPUT`: the input sample file in uBAM or uCRAM format. One or more files can be supplied by passing multiple files after the `-i` argument. -- `--align`: re-align the input read data to the reference genome using Sentieon minimap2. - -The DNAscope LongRead pipeline accepts the following new optional arguments: -- `--input_ref INPUT_REF`: a reference fasta used for decoding the input file(s). Required with uCRAM or CRAM input. Can be different from the fasta used with the `-r` argument. - -### Variant calling from BAM or CRAM - -A single command is run to call SNVs, indels, and structural variants from PacBio HiFi or ONT reads in the BAM, or CRAM formats: -```sh -sentieon-cli dnascope-longread [-h] \ - -r REFERENCE \ - -i SAMPLE_INPUT ... \ - -m MODEL_BUNDLE \ - [-d DBSNP] \ - [-b DIPLOID_BED] \ - [-t NUMBER_THREADS] \ - [-g] \ - --tech HiFi|ONT \ - [--haploid_bed HAPLOID_BED] \ - [--cnv_excluded_regions CNV_EXCLUDE_BED] \ - sample.vcf.gz -``` - -Not supplying the `--align` argument will direct the pipeline to call variants directly from the input reads. - -## Pipeline output - -### List of output files - -The following files are output when processing FASTQ data or uBAM, uCRAM, BAM, or CRAM files with the `--align` argument: -- `sample.vcf.gz`: SNV and indel variant calls across the regions of the genome as defined in the `-b DIPLOID_BED` file. -- `sample.sv.vcf.gz`: structural variant calls from the Sentieon LongReadSV tool. -- `sample_mm2_sorted_fq_*.cram`: aligned and coordinate-sorted reads from the input FASTQ files. -- `sample_mm2_sorted_*.cram`: aligned and coordinate-sorted reads from the input uBAM, uCRAM, BAM, or CRAM files. -- `sample.hificnv`: the base name of HiFiCNV output files. - -## Other considerations - -### Diploid and haploid variant calling - -The default pipeline is recommended for use with samples from diploid organisms. For samples with both diploid and haploid chromosomes, the `-b DIPLOID_BED` option can be used to limit diploid variant calling to diploid chromosomes and the `--haploid_bed HAPLOID_BED` argument can be used to perform haploid variant calling across haploid chromosomes. Diploid and haploid variants will be output to separate VCF files. - -Diploid and haploid BED files for the human hg38 reference genome (with male samples) can be found in the [/data](/data) folder in this repository. - -### Modification - -Scripts in this repository are made available under the [BSD 2-Clause license](/LICENSE). - -The Python scripts in the `sentieon_cli/scripts` folder perform low-level manipulation of intermediate gVCF and VCF files generated by the pipeline. Due to the low-level data handling performed by these scripts, modification of these files by users is discouraged. - -## References -**[Sentieon DNAscope LongRead – A highly Accurate, Fast, and Efficient Pipeline for Germline Variant Calling from PacBio HiFi reads]** - A preprint describing the DNAscope LongRead pipeline for calling variants from PacBio HiFi data. - - -[Python]: https://www.python.org/ -[bcftools]: http://samtools.github.io/bcftools/bcftools.html -[bedtools]: https://bedtools.readthedocs.io/en/latest/ -[mosdepth]: https://github.com/brentp/mosdepth -[samtools]: https://www.htslib.org/ -[sentieon-models]: https://github.com/Sentieon/sentieon-models -[hificnv]: https://github.com/PacificBiosciences/HiFiCNV - -[Sentieon DNAscope LongRead – A highly Accurate, Fast, and Efficient Pipeline for Germline Variant Calling from PacBio HiFi reads]: https://www.biorxiv.org/content/10.1101/2022.06.01.494452v1 diff --git a/docs/dnascope.md b/docs/dnascope.md deleted file mode 100644 index 8e70936..0000000 --- a/docs/dnascope.md +++ /dev/null @@ -1,176 +0,0 @@ -# DNAscope - -Sentieon DNAscope is a pipeline for alignment and germline variant calling (SNVs, SVs and indels) from short-read DNA sequence data. The DNAscope pipeline uses a combination of traditional statistical approaches and machine learning to achieve high variant calling accuracy. The DNAscope pipeline supports samples sequenced using whole-genome or targeted (hybrid-capture) enrichment library preps. - -The pipeline accepts as input aligned reads in BAM or CRAM format, or un-aligned reads in FASTQ, uBAM, or uCRAM format. The pipeline will output variants in the VCF (or gVCF) formats and aligned reads in BAM or CRAM formats. - -DNAscope is implemented using 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 202308 or higher. -- [samtools] version 1.16 or higher for alignment of reads in uBAM or uCRAM format or re-alignment of previously aligned reads. -- [MultiQC] version 1.18 or higher for metrics report generation. - -The `sentieon`, `samtools`, and `multiqc` executables will be accessed through the user's `PATH` environment variable. - -## Input data requirements - -### The Reference genome - -DNAscope 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. Read alignment also requires bwa index files. - -We recommend aligning to a reference genome without alternate contigs. If alternate contigs are present in the genome, please also supply a ".alt" file to activate [alt-aware alignment] in bwa. - -## Usage - -### Alignment and variant calling from FASTQ - -A single command is run to align, preprocess, and call SNVs, indels, and structural variants from FASTQ: -```sh -sentieon-cli dnascope [-h] \ - -r REFERENCE \ - --r1_fastq R1_FASTQ ... \ - --r2_fastq R2_FASTQ ... \ - --readgroups READGROUPS ... \ - -m MODEL_BUNDLE \ - [-d DBSNP] \ - [-b BED] \ - [--interval_padding INTERVAL_PADDING] \ - [-t NUMBER_THREADS] \ - [--pcr_free] \ - [-g] \ - [--duplicate_marking DUP_MARKING] \ - [--assay ASSAY] \ - [--consensus] \ - [--dry_run] \ - [--bam_format] \ - sample.vcf.gz -``` - -With FASTQ input, the DNAscope pipeline requires the following arguments: -- `-r REFERENCE`: the location of the reference FASTA file. A reference fasta index, ".fai" file, and bwa index files, are also required. -- `--r1_fastq R1_FASTQ`: the R1 input FASTQ. Can be used multiple times. `--r1_fastq` files without a corresponding `--r2_fastq` are assumed to be single-ended. Be aware that the pipeline performs single-sample processing, and all fastq are expected to be from the same sample. -- `--r2_fastq R2_FASTQ`: the R2 input FASTQ. Can be used multiple times. -- `--readgroups READGROUPS`: readgroup information for each FASTQ. The pipeline will expect the same number of arguments to `--r1_fastq` and `--readgroups`. An example argument is, `--readgroups "@RG\tID:HG002-1\tSM:HG002\tLB:HG002-LB-1\tPL:ILLUMINA"` -- `-m MODEL_BUNDLE`: the location of the model bundle. Model bundle files can be found in the [sentieon-models] repository. -- `sample.vcf.gz`: the location of the output VCF file for SNVs and indels. 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 DNAscope pipeline accepts the following optional arguments: -- `-d DBSNP`: the location of the Single Nucleotide Polymorphism database (dbSNP) used to label known variants in VCF (`.vcf`) or bgzip compressed VCF (`.vcf.gz`) format. Only one file is supported. Supplying this file will annotate variants with their dbSNP refSNP ID numbers. A VCF index file is required. -- `-b BED`: interval in the reference to restrict variant calling, in BED file format. Supplying this file will limit variant calling to the intervals inside the BED file. If a BED file is not supplied, the software will process the whole genome. -- `--interval_padding INTERVAL_PADDING`: adds INTERVAL_PADDING bases padding to the edges of the input intervals. The default value is 0. -- `-t NUMBER_THREADS`: number of computing threads that will be used by the software to run parallel processes. The argument is optional; if omitted, the pipeline will use as many threads as the server has. -- `--pcr_free`: Call variants using `--pcr_indel_model NONE`, which is appropriate for libraries prepared with a PCR-free library prep. Deduplication is still performed to identify optical duplicates. -- `-g`: output variants in the gVCF format, in addition to the VCF output file. The tool will output a bgzip compressed gVCF file with a corresponding index file. -- `--duplicate_marking DUP_MARKING`: setting for duplicate marking. `markdup` will mark duplicate reads. `rmdup` will remove duplicate reads. `none` will skip duplicate marking. The default setting is `markdup`. -- `--assay ASSAY`: assay setting for metrics collection `WGS` or `WES`. The default setting is `WGS`. -- `--consensus`: generate consensus reads during duplicate marking. -- `-h`: print the command-line help and exit. -- `--dry_run`: print the pipeline commands, but do not actually execute them. -- `--bam_format`: use BAM format instead of CRAM for output aligned files. - -### Alignment and variant calling from uBAM or uCRAM -A single command is run to align, preprocess, and call SNVs, indels, and structural variants from uBAM or uCRAM files: -```sh -sentieon-cli dnascope [-h] \ - -r REFERENCE \ - -i SAMPLE_INPUT ... \ - --align \ - [--input_ref INPUT_REF] \ - -m MODEL_BUNDLE \ - [-d DBSNP] \ - [-b BED] \ - [--interval_padding INTERVAL_PADDING] \ - [-t NUMBER_THREADS] \ - [--pcr_free] \ - [-g] \ - [--duplicate_marking DUP_MARKING] \ - [--assay ASSAY] \ - [--consensus] \ - [--dry_run] \ - [--bam_format] \ - sample.vcf.gz -``` - -With uBAM or uCRAM input, the DNAscope pipeline requires the following new arguments: -- `-i SAMPLE_INPUT`: the input sample file in uBAM or uCRAM format. One or more files can be supplied by passing multiple files after the `-i` argument. -- `--align`: directs the pipeline to align the input reads. - -The DNAscope pipeline accepts the following new optional arguments: -- `--input_ref INPUT_REF`: a reference fasta used for decoding the input file(s). Required with uCRAM input. Can be different from the fasta used with the `-r` argument. - -### Alignment and variant calling from sorted BAM or CRAM -A single command is run to align, preprocess, and call SNVs, indels, and structural variants from BAM or CRAM files: -```sh -sentieon-cli dnascope [-h] \ - -r REFERENCE \ - -i SAMPLE_INPUT ... \ - --collate_align \ - [--input_ref INPUT_REF] \ - -m MODEL_BUNDLE \ - [-d DBSNP] \ - [-b BED] \ - [--interval_padding INTERVAL_PADDING] \ - [-t NUMBER_THREADS] \ - [--pcr_free] \ - [-g] \ - [--duplicate_marking DUP_MARKING] \ - [--assay ASSAY] \ - [--consensus] \ - [--dry_run] \ - [--bam_format] \ - sample.vcf.gz -``` - -With BAM or CRAM input, the DNAscope pipeline requires the following new arguments: -- `--collate_align`: directs the pipeline to collate and then align the input reads. - -### Variant calling from sorted BAM or CRAM -A single command is run to preprocess, and call SNVs, indels, and structural variants from BAM or CRAM files: -```sh -sentieon-cli dnascope [-h] \ - -r REFERENCE \ - -i SAMPLE_INPUT ... \ - -m MODEL_BUNDLE \ - [-d DBSNP] \ - [-b BED] \ - [--interval_padding INTERVAL_PADDING] \ - [-t NUMBER_THREADS] \ - [--pcr_free] \ - [-g] \ - [--duplicate_marking DUP_MARKING] \ - [--assay ASSAY] \ - [--consensus] \ - [--dry_run] \ - [--bam_format] \ - sample.vcf.gz -``` - -Not supplying the `--align` and `--collate_align` arguments will direct the pipeline to call variants directly from the input reads. - -## Pipeline output - -### List of output files - -The following files are output when processing WGS FASTQ with default arguments: -- `sample.vcf.gz`: SNV and indel variant calls across the regions of the genome as defined in the `-b BED` file. -- `sample_deduped.cram` or `sample_deduped.bam`: aligned, coordinate-sorted and duplicate-marked read data from the input FASTQ files. -- `sample_svs.vcf.gz`: structural variant calls from DNAscope and SVSolver. -- `sample_metrics`: a directory containing QC metrics for the analyzed sample. - - `sample_metrics/coverage*`: coverage metrics for the processed sample. Only available for WGS samples. - - `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. Only available for WGS samples. - - `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. Only available for WGS samples. - - `sample_metrics/{sample}.txt.hybrid-selection.txt`: Metrics from the HsMetricAlgo algo. - - `sample_metrics/multiqc_report.html`: collected QC metrics aggregated by MultiQC. - -[samtools]: https://www.htslib.org/ -[MultiQC]: https://multiqc.info/ -[alt-aware alignment]: https://github.com/lh3/bwa/blob/master/README-alt.md -[sentieon-models]: https://github.com/Sentieon/sentieon-models diff --git a/pyproject.toml b/pyproject.toml index 734c3b4..8f866cb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,7 @@ [project] name = "sentieon_cli" -version = "1.6.3" +version = "1.7.0" description = "Pipeline implementations for the Sentieon software" authors = [ {name = "Don Freed", email = "don.freed@sentieon.com"}, diff --git a/sentieon_cli/command_strings.py b/sentieon_cli/command_strings.py index 86d00bf..4bdebde 100644 --- a/sentieon_cli/command_strings.py +++ b/sentieon_cli/command_strings.py @@ -368,9 +368,11 @@ def hybrid_stage1( def hybrid_stage3( - out_bam: pathlib.Path, + out_aln: pathlib.Path, + reference: pathlib.Path, driver: BaseDriver, cores: int, + util_sort_args: str = "", ) -> Pipeline: sort_cmd = Command( "sentieon", @@ -380,8 +382,11 @@ def hybrid_stage3( "-", "-t", str(cores), + "--reference", + str(reference), "-o", - str(out_bam), + str(out_aln), + *util_sort_args.split(), ) return Pipeline(Command(*driver.build_cmd()), sort_cmd) @@ -1235,6 +1240,7 @@ def cmd_segdup_caller( input_vcf: Optional[pathlib.Path] = None, sex: Optional[str] = None, genes: Optional[str] = None, + overrides: Optional[List[str]] = None, ) -> Pipeline: cmd = [ "segdup-caller", @@ -1251,6 +1257,8 @@ def cmd_segdup_caller( cmd.extend(["--sex", sex]) if genes: cmd.extend(["--genes", genes]) + for override in overrides or []: + cmd.extend(["--set", override]) cmd.extend(["--outdir", str(out_segdup)]) return Pipeline(Command(*cmd)) diff --git a/sentieon_cli/dnascope_hybrid.py b/sentieon_cli/dnascope_hybrid.py index 958095b..e1a9e28 100644 --- a/sentieon_cli/dnascope_hybrid.py +++ b/sentieon_cli/dnascope_hybrid.py @@ -916,7 +916,10 @@ def call_variants( rm_cmd = ["rm", str(stage1_bam), str(stage1_hap_bam)] rm_job3 = Job(Pipeline(Command(*rm_cmd, fail_ok=True)), "rm-tmp3", 0) - stage3_bam = self.tmp_dir.joinpath("hybrid_stage3.bam") + suffix = "bam" if self.bam_format else "cram" + stage3_aln = pathlib.Path( + str(self.output_vcf).replace(".vcf.gz", f"_sr_realigned.{suffix}") + ) driver = Driver( reference=self.reference, thread_count=self.cores, @@ -933,7 +936,8 @@ def call_variants( ) third_stage_job = Job( cmds.hybrid_stage3( - stage3_bam, + stage3_aln, + reference=self.reference, driver=driver, cores=self.cores, ), @@ -949,7 +953,7 @@ def call_variants( reference=self.reference, thread_count=self.cores, replace_rg=rg_info.replace_rg_args[0], - input=lr_aln + [stage3_bam], + input=lr_aln + [stage3_aln], interval=stage2_bed, read_filter=rg_info.ultima_read_filter + rg_info.lr_rg_read_filter, ) diff --git a/sentieon_cli/sentieon_pangenome.py b/sentieon_cli/sentieon_pangenome.py index 53393c3..55616b8 100644 --- a/sentieon_cli/sentieon_pangenome.py +++ b/sentieon_cli/sentieon_pangenome.py @@ -68,7 +68,7 @@ } SEGDUP_MIN_VERSION = { - "segdup-caller": packaging.version.Version("0.5.1"), + "segdup-caller": packaging.version.Version("0.7.0"), } EXPANSION_MIN_VERSION = { @@ -743,8 +743,23 @@ def build_first_dag(self) -> DAG: if self.r1_fastq: dnascope_bams.append(out_bwa_aln) dnascope_bams.append(out_mm2_aln) + + # 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" + 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 @@ -1058,6 +1073,7 @@ def build_dedup_job( 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") @@ -1090,7 +1106,7 @@ def build_dedup_job( input=input_bam, read_filter=read_filters, ) - driver2.add_algo(Dedup(output_bam, score_file)) + driver2.add_algo(Dedup(output_bam, score_file, metrics=metrics)) dedup_job = Job( Pipeline(Command(*driver2.build_cmd())), @@ -1287,6 +1303,13 @@ def build_segdup_job( sys.exit(2) sex = "male" if self.sample_sex == SampleSex.MALE else "female" + + # segdup-caller's default `main.min_map_qual` of 45 is too strict + # for Ultima alignments. + overrides = [] + if self.tech.upper() == "ULTIMA": + overrides.append("main.min_map_qual=30") + return Job( cmds.cmd_segdup_caller( out_segdup, @@ -1296,6 +1319,7 @@ def build_segdup_job( input_vcf=input_vcf, sex=sex, genes=genes, + overrides=overrides, ), "segdup-caller", self.cores, diff --git a/sentieon_cli/util.py b/sentieon_cli/util.py index c8c3f8f..68021ea 100644 --- a/sentieon_cli/util.py +++ b/sentieon_cli/util.py @@ -54,9 +54,15 @@ def check_version( return True cmd_list.append("--version") - cmd_version_str = ( - sp.check_output(cmd_list).decode("utf-8", "ignore").strip() - ) + try: + cmd_version_str = ( + sp.check_output(cmd_list).decode("utf-8", "ignore").strip() + ) + except (sp.CalledProcessError, OSError) as e: + logger.error( + "Error: could not determine the version of '%s': %s", cmd, e + ) + return False if cmd_list[0] == "sentieon": cmd_version_str = cmd_version_str.split("-")[-1] elif cmd_list[0] == "pbsv": diff --git a/tests/test_util.py b/tests/test_util.py index 0395b8c..c835f55 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -1,5 +1,9 @@ import os +import subprocess as sp import sys +from unittest.mock import patch + +import packaging.version sys.path.insert( 0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) @@ -28,3 +32,63 @@ def test_split_numa_nodes(): "48-55", "56-63", ] + + +def test_check_version_parses_version(): + """A well-behaved `--version` is compared against the minimum""" + min_version = packaging.version.Version("0.9.0") + with ( + patch.object( + sentieon_cli.util.shutil, + "which", + return_value="/bin/segdup-caller", + ), + patch.object( + sentieon_cli.util.sp, + "check_output", + return_value=b"segdup-caller 0.9.0\n", + ), + ): + assert sentieon_cli.util.check_version("segdup-caller", min_version) + + with ( + patch.object( + sentieon_cli.util.shutil, + "which", + return_value="/bin/segdup-caller", + ), + patch.object( + sentieon_cli.util.sp, + "check_output", + return_value=b"segdup-caller 0.8.0\n", + ), + ): + assert not sentieon_cli.util.check_version( + "segdup-caller", min_version + ) + + +def test_check_version_handles_failed_command(): + """A non-zero `--version` reports False rather than raising. + + segdup-caller >=0.6.0 exits non-zero on `--version` when + SENTIEON_LICENSE is unset. + """ + min_version = packaging.version.Version("0.9.0") + for err in ( + sp.CalledProcessError(1, ["segdup-caller", "--version"]), + OSError("Exec format error"), + ): + with ( + patch.object( + sentieon_cli.util.shutil, + "which", + return_value="/bin/segdup-caller", + ), + patch.object( + sentieon_cli.util.sp, "check_output", side_effect=err + ), + ): + assert not sentieon_cli.util.check_version( + "segdup-caller", min_version + ) diff --git a/tests/unit/test_sentieon_pangenome.py b/tests/unit/test_sentieon_pangenome.py index c85a799..e6c1914 100644 --- a/tests/unit/test_sentieon_pangenome.py +++ b/tests/unit/test_sentieon_pangenome.py @@ -14,6 +14,7 @@ sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))) from sentieon_cli.sentieon_pangenome import SentieonPangenome +from sentieon_cli.base_pangenome import SampleSex from sentieon_cli.dag import DAG @@ -91,6 +92,57 @@ def create_pipeline(self): return pipeline + def create_fastq_pipeline(self): + """Create a fastq-input SentieonPangenome pipeline for testing. + + The default fixture uses BAM/CRAM input, which does not exercise the + dedup/metrics branch of ``build_first_dag``. This variant provides + FASTQ input so the bwa/mm2 dedup jobs are created. + """ + mock_r1 = self.mock_dir / "sample_R1.fastq.gz" + mock_r2 = self.mock_dir / "sample_R2.fastq.gz" + for fq in (mock_r1, mock_r2): + fq.touch() + + pipeline = self.create_pipeline() + pipeline.sample_input = [] + pipeline.r1_fastq = [mock_r1] + pipeline.r2_fastq = [mock_r2] + pipeline.fastq_readgroup = {"ID": "rg1", "SM": "sample1"} + pipeline.skip_metrics = False + pipeline.skip_multiqc = True + return pipeline + + def test_dedup_metrics_output(self): + """Dedup on the primary (bwa) alignment emits a --metrics file that + lands in the metrics directory scanned by MultiQC.""" + pipeline = self.create_fastq_pipeline() + dag = pipeline.build_first_dag() + + job_names, all_jobs = self._get_all_job_names(dag) + assert "dedup-bwa" in job_names + assert "dedup-mm2" in job_names + + bwa_dedup = next(j for j in all_jobs if j.name == "dedup-bwa") + bwa_cmd = str(bwa_dedup.shell) + assert "--algo Dedup" in bwa_cmd + assert "--metrics" in bwa_cmd + assert "output_metrics/output.txt.dedup_metrics.txt" in bwa_cmd + + # The mm2 dedup does not emit a duplicate-metrics file + mm2_dedup = next(j for j in all_jobs if j.name == "dedup-mm2") + assert "--metrics" not in str(mm2_dedup.shell) + + def test_dedup_metrics_skipped(self): + """No Dedup --metrics output when metrics collection is skipped.""" + pipeline = self.create_fastq_pipeline() + pipeline.skip_metrics = True + dag = pipeline.build_first_dag() + + _, all_jobs = self._get_all_job_names(dag) + bwa_dedup = next(j for j in all_jobs if j.name == "dedup-bwa") + assert "--metrics" not in str(bwa_dedup.shell) + def test_model_apply_default(self): """Test that model apply job is created by default""" pipeline = self.create_pipeline() @@ -171,6 +223,27 @@ def test_no_gvcftyper_without_gvcf(self): job_names = [job.name for job in all_jobs] assert "gvcftyper" not in job_names + def build_segdup_cmd(self, tech): + """Build the segdup-caller command string for a given platform""" + pipeline = self.create_pipeline() + pipeline.tech = tech + pipeline.sample_sex = SampleSex.FEMALE + job = pipeline.build_segdup_job( + self.mock_dir / "output_segdups", + self.mock_bam, + self.mock_vcf, + None, + ) + return str(job.shell) + + def test_segdup_ultima_lowers_min_map_qual(self): + """Ultima input overrides segdup-caller's default min_map_qual""" + assert "--set main.min_map_qual=30" in self.build_segdup_cmd("Ultima") + + def test_segdup_no_override_for_short_reads(self): + """Non-Ultima input leaves the segdup-caller defaults alone""" + assert "--set" not in self.build_segdup_cmd("Illumina") + def test_call_svs(self): """Test that PangenomeSV is added when --call_svs is enabled""" pipeline = self.create_pipeline()