diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a5767b4..9c1f46a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -13,7 +13,7 @@ jobs: strategy: matrix: NXF_VER: - - "23.10.0" + - "25.04.2" steps: - name: Check out pipeline code uses: actions/checkout@v4 diff --git a/conf/test.config b/conf/test.config index 2c2d694..c22bb46 100644 --- a/conf/test.config +++ b/conf/test.config @@ -29,4 +29,6 @@ params { email = "test@email.com" publish_dir_mode = "symlink" tracedir = "test_output/pipeline_info" + + deepvariant_model_type = 'WGS' } diff --git a/modules/local/deepvariant_checkpoint.nf b/modules/local/deepvariant_checkpoint.nf new file mode 100644 index 0000000..472ba77 --- /dev/null +++ b/modules/local/deepvariant_checkpoint.nf @@ -0,0 +1,66 @@ +process fetchDeepvariantCheckpoint { + label 'process_single' + tag "${model_type}" + conda 'conda-forge::wget=1.20.3' + storeDir "${projectDir}/ref_data/deepvariant_checkpoints" + // + // Pulls the EM-seq DeepVariant fine-tuned model for the requested --model_type (only WGS and WES) from + // Zenodo (doi.org/10.5281/zenodo.21416823) and stages it as a checkpoint directory + // WGS -> wgs.checkpoint.260717.{data-*,index} + // WES -> target_capture.checkpoint.260717.{data-*,index} + // + // The output dir contains the checkpoint (.index/.data-*) plus model.example_info.json. + // DV 1.10.0 inference resolves it by directory: get_model_example_info_json() looks in + // os.path.dirname(--customized_model) for 'model.example_info.json' + + input: + val(model_type) + + output: + tuple val(model_type), path("${label}"), val(prefix), emit: checkpoint + + script: + def spec = [ WGS: ['wgs', 'wgs.checkpoint.260717'], + WES: ['target_capture', 'target_capture.checkpoint.260717'] ][model_type] + if (!spec) + error "No EM-seq Zenodo checkpoint for model_type '${model_type}' (only WGS and WES)." + label = spec[0] + prefix = spec[1] + def base = 'https://zenodo.org/records/21416823/files' + """ + mkdir -p ${label} + + wget -q -O ${label}/${prefix}.data-00000-of-00001 "${base}/${prefix}.data-00000-of-00001?download=1" + wget -q -O ${label}/${prefix}.index "${base}/${prefix}.index?download=1" + wget -q -O ${label}/model.example_info.json "${base}/model.example_info.json?download=1" + """ +} + +workflow resolveDeepvariantCheckpoint { + main: + def mt = params.deepvariant_model_type + def valid = ['WGS', 'WES', 'PACBIO', 'ONT_R104', 'HYBRID_PACBIO_ILLUMINA'] + if (!mt) + error "Please set --deepvariant_model_type (one of: ${valid.join(', ')})." + if (!(mt in valid)) + error "deepvariant_model_type must be one of ${valid.join(', ')}, got '${mt}'." + if (params.deepvariant_custom_model && !(mt in ['WGS', 'WES'])) + error "deepvariant_custom_model is only available for WGS and WES — " + + "model_type '${mt}' has no custom model." + + if (params.run_deepvariant == false) { + ch_checkpoint = Channel.value([[], '', mt]) + } else if (params.deepvariant_checkpoint) { + ch_checkpoint = Channel.value([ file(params.deepvariant_checkpoint).parent, + file(params.deepvariant_checkpoint).name, + mt ]) + } else if (params.deepvariant_custom_model) { + ch_checkpoint = fetchDeepvariantCheckpoint(Channel.value(mt)) + .checkpoint.map { m, dir, prefix -> [dir, prefix, mt] } + } else { + ch_checkpoint = Channel.value([[], '', mt]) + } + + emit: + checkpoint = ch_checkpoint +} diff --git a/modules/local/deepvariant_custom.nf b/modules/local/deepvariant_custom.nf new file mode 100644 index 0000000..d994e5e --- /dev/null +++ b/modules/local/deepvariant_custom.nf @@ -0,0 +1,35 @@ +process deepvariantCustom { + label 'process_high' + tag "$library" + container "docker.io/google/deepvariant:1.10.0" + publishDir "${params.outdir}/deepvariantCustom", mode: params.publish_dir_mode + // For running a custom model the checkpoint directory must contain: + // - checkpoint files (.index, .data-*) + // - model.example_info.json (required by DV 1.10.0 ) + + input: + tuple val(library), path(bam), path(bai) + path(fasta) + path(fai) + tuple path(checkpoint_dir), val(checkpoint_name), val(model_type) + + output: + tuple val(library), path("${library}.deepvariant.vcf.gz"), + path("${library}.deepvariant.vcf.gz.tbi"), emit: vcf + + script: + def custom_model_opt = checkpoint_name + ? "--customized_model ${checkpoint_dir}/${checkpoint_name} --disable_small_model" + : '' + """ + export TMPDIR=/tmp + + run_deepvariant \\ + --model_type ${model_type} \\ + ${custom_model_opt} \\ + --ref ${fasta} \\ + --reads ${bam} \\ + --output_vcf ${library}.deepvariant.vcf.gz \\ + --num_shards ${task.cpus} +""" +} diff --git a/modules/local/fgbio_clip_bam.nf b/modules/local/fgbio_clip_bam.nf new file mode 100644 index 0000000..998a24c --- /dev/null +++ b/modules/local/fgbio_clip_bam.nf @@ -0,0 +1,29 @@ +process fgbioClipBam { + cpus 6 + memory '100 GB' + conda "bioconda::fgbio=3.1.2 bioconda::samtools=1.19" + publishDir "${params.outdir}/clipped_bams/", mode: 'symlink' + tag { sample_id } + + input: + tuple val(sample_id), path(bam), path(bai) + + output: + tuple val(sample_id), path("${sample_id}.clipped.bam"), path("${sample_id}.clipped.bam.bai") + + script: + """ + export _JAVA_OPTIONS="-Xmx90g" + samtools sort -n -u -@ ${task.cpus} ${bam} | \\ + fgbio ClipBam \\ + --input /dev/stdin \\ + --output /dev/stdout \\ + --ref ${params.fasta} \\ + --clip-overlapping-reads true \\ + --clipping-mode Hard \\ + --sort-order Coordinate | \\ + samtools view -b -@ ${task.cpus} -o ${sample_id}.clipped.bam + + samtools index ${sample_id}.clipped.bam + """ +} diff --git a/modules/local/strelka.nf b/modules/local/strelka.nf index c88fe1b..847ebb6 100644 --- a/modules/local/strelka.nf +++ b/modules/local/strelka.nf @@ -13,7 +13,7 @@ process strelka { tuple val(library), path("*.strelka_vcf.gz"), path("*.strelka_vcf.gz.tbi"), emit: vcf script: - def target_regions_opt = (target_bed && file(target_bed).exists()) ? '--callRegions target.bed.gz' : '' + def target_regions_opt = (target_bed && file(target_bed).exists()) ? '--callRegions target.bed.gz --exome' : '' """ if [ -n "${target_bed}" ] && [ -f "${target_bed}" ]; then ln -s "${target_bed}" target.bed diff --git a/nextflow.config b/nextflow.config index 3b8e092..d4c9a06 100644 --- a/nextflow.config +++ b/nextflow.config @@ -8,18 +8,23 @@ // Global default params, used in configs params { - tmpdir = '/tmp' - index_format = 'csi' - publish_dir_mode = 'copy' - outdir = 'results' - tracedir = "${params.outdir}/pipeline_info" - run_strelka = true - run_freebayes = true - freebayes_qual_filter = 15 // default quality filter for freebayes output. - run_happy = true - target_bed = null - happy_truth_vcf = null - happy_bed = null + tmpdir = '/tmp' + index_format = 'csi' + publish_dir_mode = 'copy' + outdir = 'results' + tracedir = "${params.outdir}/pipeline_info" + run_clip_bam = false + run_strelka = true + run_freebayes = true + run_deepvariant = true + deepvariant_model_type = null // WGS | WES | PACBIO | ONT_R104 | HYBRID_PACBIO_ILLUMINA + deepvariant_custom_model = true // Download an EM-seq fine-tuned model, only available for WGS and WES (target capture). + deepvariant_checkpoint = null // Optional local custom model prefix + freebayes_qual_filter = 15 + run_happy = true + target_bed = null + happy_truth_vcf = null + happy_bed = null } includeConfig 'conf/base.config' @@ -35,7 +40,8 @@ profiles { conda.cacheDir = "${projectDir}/.conda/envs" conda.enabled = true docker.enabled = false - singularity.enabled = false + singularity.enabled = true + singularity.autoMounts = true podman.enabled = false shifter.enabled = false charliecloud.enabled = false @@ -53,10 +59,12 @@ profiles { apptainer.enabled = false } - test { + test { includeConfig 'conf/test.config' conda.enabled = true conda.cacheDir = "${projectDir}/.conda/envs" + docker.enabled = true + singularity.enabled = false } } diff --git a/tests/main.nf.test b/tests/main.nf.test index 2f3c844..8222681 100644 --- a/tests/main.nf.test +++ b/tests/main.nf.test @@ -27,6 +27,8 @@ nextflow_pipeline { def happy_vcf_freebayes = path("${launchDir}/test_output/happyConcordanceFreebayes/200ng-V2-1.md_vs_truth.vcf.happy.vcf.gz").vcf.summary def happy_summary_freebayes = path("${launchDir}/test_output/happyConcordanceFreebayes/200ng-V2-1.md_vs_truth.vcf.happy.summary.csv").text.tokenize('\n') def snp_breakdown = path("${launchDir}/test_output/parseHappyVcf/200ng-V2-1.md.happy.snp_mutations.txt").text.tokenize('\n') + def deepvariant_vcf = path("${launchDir}/test_output/deepvariantCustom/200ng-V2-1.md.deepvariant.vcf.gz").vcf.summary + def happy_vcf_deepvariant = path("${launchDir}/test_output/happyConcordanceDeepvariant/200ng-V2-1.md_vs_truth.vcf.happy.vcf.gz").vcf.summary assertAll( { assert workflow.success }, @@ -40,7 +42,9 @@ nextflow_pipeline { ["happy_summary_strelka", happy_summary_strelka], ["happy_vcf_freebayes", happy_vcf_freebayes], ["happy_summary_freebayes", happy_summary_freebayes], - ["snp_breakdown", snp_breakdown] + ["snp_breakdown", snp_breakdown], + ["deepvariant_vcf", deepvariant_vcf], + ["happy_vcf_deepvariant", happy_vcf_deepvariant] ).match() } ) @@ -72,6 +76,8 @@ nextflow_pipeline { def happy_vcf_freebayes = path("${launchDir}/test_output/happyConcordanceFreebayes/200ng-V2-1.md_vs_truth.vcf.happy.vcf.gz").vcf.summary def happy_summary_freebayes = path("${launchDir}/test_output/happyConcordanceFreebayes/200ng-V2-1.md_vs_truth.vcf.happy.summary.csv").text.tokenize('\n') def snp_breakdown = path("${launchDir}/test_output/parseHappyVcf/200ng-V2-1.md.happy.snp_mutations.txt").text.tokenize('\n') + def deepvariant_vcf = path("${launchDir}/test_output/deepvariantCustom/200ng-V2-1.md.deepvariant.vcf.gz").vcf.summary + def happy_vcf_deepvariant = path("${launchDir}/test_output/happyConcordanceDeepvariant/200ng-V2-1.md_vs_truth.vcf.happy.vcf.gz").vcf.summary assertAll( { assert workflow.success }, @@ -85,7 +91,9 @@ nextflow_pipeline { ["happy_summary_strelka", happy_summary_strelka], ["happy_vcf_freebayes", happy_vcf_freebayes], ["happy_summary_freebayes", happy_summary_freebayes], - ["snp_breakdown", snp_breakdown] + ["snp_breakdown", snp_breakdown], + ["deepvariant_vcf", deepvariant_vcf], + ["happy_vcf_deepvariant", happy_vcf_deepvariant] ).match() } ) diff --git a/tests/main.nf.test.snap b/tests/main.nf.test.snap index e294762..78f1080 100644 --- a/tests/main.nf.test.snap +++ b/tests/main.nf.test.snap @@ -3,8 +3,8 @@ "content": [ { "tasksFailed": 0, - "tasksCount": 15, - "tasksSucceeded": 15 + "tasksCount": 20, + "tasksSucceeded": 20 }, [ "calcmd_bam", @@ -54,22 +54,30 @@ "snp_breakdown", [ "Count Ref Alt Type", - "1 G C Correct" + "1 G C NotPassed" ] + ], + [ + "deepvariant_vcf", + "VcfFile [chromosomes=[chr21], sampleCount=1, variantCount=60, phased=false, phasedAutodetect=false]" + ], + [ + "happy_vcf_deepvariant", + "VcfFile [chromosomes=[chr21], sampleCount=2, variantCount=4, phased=false, phasedAutodetect=false]" ] ], "meta": { - "nf-test": "0.8.4", - "nextflow": "23.04.1" + "nf-test": "0.9.0", + "nextflow": "24.04.2" }, - "timestamp": "2025-05-20T13:36:17.023260625" + "timestamp": "2026-07-17T17:08:12.351084775" }, "Pipeline test without optional inputs": { "content": [ { "tasksFailed": 0, - "tasksCount": 15, - "tasksSucceeded": 15 + "tasksCount": 19, + "tasksSucceeded": 19 }, [ "calcmd_bam", @@ -121,12 +129,20 @@ "Count Ref Alt Type", "1 G C Correct" ] + ], + [ + "deepvariant_vcf", + "VcfFile [chromosomes=[chr21], sampleCount=1, variantCount=60, phased=false, phasedAutodetect=false]" + ], + [ + "happy_vcf_deepvariant", + "VcfFile [chromosomes=[chr21], sampleCount=2, variantCount=4, phased=false, phasedAutodetect=false]" ] ], "meta": { - "nf-test": "0.8.4", - "nextflow": "23.04.1" + "nf-test": "0.9.0", + "nextflow": "24.04.2" }, - "timestamp": "2025-05-20T13:39:05.639276115" + "timestamp": "2026-07-17T17:11:27.137540031" } } \ No newline at end of file diff --git a/workflows/emseq_variant_calling.nf b/workflows/emseq_variant_calling.nf index 163881e..d430aea 100644 --- a/workflows/emseq_variant_calling.nf +++ b/workflows/emseq_variant_calling.nf @@ -25,6 +25,8 @@ tmpdir = Channel.value(params.tmpdir) // Optional input for regions to call, e.g. for target enrichment target_bed = Channel.value(params.target_bed ?: '') +// DeepVariant checkpoint — resolved to (dir, name): local override, else Zenodo fetch by mode. + // Hap.py concordance inputs, only if run_happy is true happy_truth_vcf = params.run_happy ? Channel.fromPath(params.happy_truth_vcf, checkIfExists: true) : Channel.empty() happy_truth_tbi = params.run_happy ? Channel.fromPath(params.happy_truth_vcf + '.tbi', checkIfExists: true) : Channel.empty() @@ -33,12 +35,16 @@ happy_bed = Channel.value(params.happy_bed ?: '') // optional regions to a // Local Modules: include { downloadRevelio } from '../modules/local/download_revelio.nf' include { calcMD } from '../modules/local/calc_md.nf' +include { fgbioClipBam } from '../modules/local/fgbio_clip_bam.nf' include { revelio } from '../modules/local/revelio.nf' include { strelka } from '../modules/local/strelka.nf' include { freebayes } from '../modules/local/freebayes.nf' include { happyConcordance as happyConcordanceStrelka; - happyConcordance as happyConcordanceFreebayes } from '../modules/local/happy_concordance.nf' + happyConcordance as happyConcordanceFreebayes; + happyConcordance as happyConcordanceDeepvariant; } from '../modules/local/happy_concordance.nf' include { parseHappyVcf } from '../modules/local/parse_happy_vcf.nf' +include { deepvariantCustom } from '../modules/local/deepvariant_custom.nf' +include { resolveDeepvariantCheckpoint } from '../modules/local/deepvariant_checkpoint.nf' workflow emseq_variant_calling { @@ -56,11 +62,18 @@ workflow emseq_variant_calling { index_format ) + // + // Module: Optionally clip overlapping reads with fgbio ClipBam before Revelio + // + if (params.run_clip_bam) { + fgbioClipBam(calcMD.out.calcmd_bam) + } + // // Module: Run Revelio to mask possibly converted bases by setting BQ to 0 // revelio ( - calcMD.out.calcmd_bam, + params.run_clip_bam ? fgbioClipBam.out : calcMD.out.calcmd_bam, downloadRevelio.out, tmpdir ) @@ -111,7 +124,7 @@ workflow emseq_variant_calling { fai, target_bed // call regions (optional) ) - + // // Module: Run hap.py to compare variants to a "truth" vcf // @@ -127,6 +140,31 @@ workflow emseq_variant_calling { ) } } + + // + // Module: Run DeepVariant with custom model + // + if (params.run_deepvariant) { + ck = resolveDeepvariantCheckpoint() + deepvariantCustom( + revelio.out.masked, + fasta, + fai, + ck.checkpoint + ) + + if (params.run_happy) { + + happyConcordanceDeepvariant( + deepvariantCustom.out.vcf, + happy_bed, // comparison regions (optional) + fasta, + fai, + happy_truth_vcf, + happy_truth_tbi + ) + } + } } workflow {