Concatenating germline vcfs

.github/workflows/pytest-workflow.yml

@@ -51,6 +51,8 @@ jobs:
             tags: snpeff
           - profile: "conda"
             tags: vep
+          - profile: "conda"
+            tags: concatenate_vcfs
           - profile: "singularity"
             tags: merge
     env:

CHANGELOG.md

@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 - [#864](https://github.com/nf-core/sarek/pull/864) - Added possibilities to export assembled haplotypes and locally realigned reads
+- [#792](https://github.com/nf-core/sarek/pull/792) - Added the option `--concatenate_vcfs` for concatenating the germline vcf-files. Per default, the resulting vcf-files will be placed under `<outDir>/variant_calling/concat`.
 
 ### Changed
 

conf/modules/post_variant_calling.config

@@ -0,0 +1,44 @@
+/*
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Config file for defining DSL2 per module options and publishing paths
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Available keys to override module options:
+        ext.args   = Additional arguments appended to command in module.
+        ext.args2  = Second set of arguments appended to command in module (multi-tool modules).
+        ext.args3  = Third set of arguments appended to command in module (multi-tool modules).
+        ext.prefix = File name prefix for output files.
+        ext.when   = When to run the module.
+----------------------------------------------------------------------------------------
+*/
+
+// POSTPROCESSING VCFS
+// Like, for instance, concatenating the unannotated, germline vcf-files
+
+process {
+    withName: 'GERMLINE_VCFS_CONCAT'{
+        publishDir       = [
+            //specify to avoid publishing, overwritten otherwise
+            enabled: false
+        ]
+    }
+
+    withName: 'GERMLINE_VCFS_CONCAT_SORT'{
+        ext.prefix       = { "${meta.id}.germline" }
+        publishDir       = [
+            mode: params.publish_dir_mode,
+                path: { "${params.outdir}/variant_calling/concat/${meta.id}/" }
+        ]
+    }
+
+    withName: 'TABIX_EXT_VCF_.*' {
+        ext.prefix       = { "${input.baseName}" }
+    }
+
+    withName: 'TABIX_GERMLINE_VCFS_CONCAT_SORT'{
+        ext.prefix       = { "${meta.id}.germline" }
+        publishDir       = [
+            mode: params.publish_dir_mode,
+                path: { "${params.outdir}/variant_calling/concat/${meta.id}/" }
+        ]
+    }
+}

modules.json

@@ -15,6 +15,10 @@
                         "git_sha": "6301e29d77e7ec7ce98b55b8a361b316a9a91bfe",
                         "installed_by": ["modules"]
                     },
+                    "bcftools/concat": {
+                        "branch": "master",
+                        "git_sha": "5e34754d42cd2d5d248ca8673c0a53cdf5624905"
+                    },
                     "bcftools/sort": {
                         "branch": "master",
                         "git_sha": "78cf39939fbe160a1410c44a6c5946f9a4c56e7e",

modules/local/add_info_to_vcf/main.nf

@@ -0,0 +1,40 @@
+process ADD_INFO_TO_VCF {
+    tag "$meta.id"
+
+    conda (params.enable_conda ? "anaconda::gawk=5.1.0" : null)
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/gawk:5.1.0' :
+        'quay.io/biocontainers/gawk:5.1.0' }"
+
+    input:
+    tuple val(meta), path(vcf_gz)
+
+    output:
+    tuple val(meta), path("*.added_info.vcf"), emit: vcf
+    path "versions.yml"                      , emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    """
+    input="input.vcf"
+    output="${vcf_gz.baseName.minus(".vcf")}.added_info.vcf"
+    zcat $vcf_gz > \$input
+    ## Add info header lines
+    grep -E "^##" \$input > \$output
+    ## Add description of new INFO value
+    echo '##INFO=<ID=SOURCE,Number=1,Type=String,Description="Name of vcf-file from whence the variant came">' >> \$output
+    ## Add column header
+    grep -E "^#CHROM" \$input >> \$output
+    ## Add SOURCE value to INFO column of variant calls
+    if grep -Ev "^#" \$input; then
+        grep -Ev "^#" \$input | awk 'BEGIN{FS=OFS="\t"} { \$8=="." ? \$8="SOURCE=$vcf_gz" : \$8=\$8";SOURCE=$vcf_gz"; print }' >> \$output
+    fi
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        gawk: \$(awk -Wversion | sed '1!d; s/.*Awk //; s/,.*//')
+    END_VERSIONS
+    """
+}

nextflow.config

@@ -63,6 +63,7 @@ params {
     cf_minqual                  = 0     // ControlFreec default values
     cf_window                   = null  // by default we are not using this in Control-FREEC
     cnvkit_reference            = null  // by default the reference is build from the fasta file
+    concatenate_vcfs            = false // by default we don't concatenate the germline-vcf-files
     ignore_soft_clipped_bases   = false // no --dont-use-soft-clipped-bases for GATK Mutect2
     wes                         = false // Set to true, if data is exome/targeted sequencing data. Used to use correct models in various variant callers
     joint_germline              = false // g.vcf & joint germline calling are not run by default if HaplotypeCaller is selected
@@ -316,6 +317,8 @@ includeConfig 'conf/modules/mutect2.config'
 includeConfig 'conf/modules/strelka.config'
 includeConfig 'conf/modules/tiddit.config'
 
+includeConfig 'conf/modules/post_variant_calling.config'
+
 //annotate
 includeConfig 'conf/modules/annotate.config'
 

nextflow_schema.json

@@ -226,6 +226,12 @@
             "default": "",
             "fa_icon": "fas fa-toolbox",
             "properties": {
+                "concatenate_vcfs": {
+                    "type": "boolean",
+                    "fa_icon": "fas fa-merge",
+                    "description": "Option for concatenating germline vcf-files.",
+                    "help_text": "Concatenating the germline vcf-files from each applied variant-caller into one vcf-file using bfctools concat."
+                },
                 "only_paired_variant_calling": {
                     "type": "boolean",
                     "fa_icon": "fas fa-forward",

subworkflows/local/bam_variant_calling_germline_all/main.nf

@@ -37,7 +37,6 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
     //TODO: Temporary until the if's can be removed and printing to terminal is prevented with "when" in the modules.config
     deepvariant_vcf     = Channel.empty()
     freebayes_vcf       = Channel.empty()
-    genotype_gvcf       = Channel.empty()
     haplotypecaller_vcf = Channel.empty()
     manta_vcf           = Channel.empty()
     mpileup_vcf         = Channel.empty()
@@ -95,7 +94,6 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
             dict
         )
 
-        mpileup_germline = BAM_VARIANT_CALLING_MPILEUP.out.mpileup
         mpileup_vcf = BAM_VARIANT_CALLING_MPILEUP.out.vcf
         ch_versions = ch_versions.mix(BAM_VARIANT_CALLING_MPILEUP.out.versions)
     }
@@ -116,7 +114,7 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
             []
         )
 
-        ch_versions     = ch_versions.mix(BAM_VARIANT_CALLING_CNVKIT.out.versions)
+        ch_versions = ch_versions.mix(BAM_VARIANT_CALLING_CNVKIT.out.versions)
     }
 
     // DEEPVARIANT
@@ -128,8 +126,8 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
             fasta_fai
         )
 
-        deepvariant_vcf = Channel.empty().mix(BAM_VARIANT_CALLING_DEEPVARIANT.out.deepvariant_vcf)
-        ch_versions     = ch_versions.mix(BAM_VARIANT_CALLING_DEEPVARIANT.out.versions)
+        deepvariant_vcf      = Channel.empty().mix(BAM_VARIANT_CALLING_DEEPVARIANT.out.deepvariant_vcf)
+        ch_versions          = ch_versions.mix(BAM_VARIANT_CALLING_DEEPVARIANT.out.versions)
     }
 
     // FREEBAYES
@@ -147,8 +145,8 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
             fasta_fai
         )
 
-        freebayes_vcf   = BAM_VARIANT_CALLING_FREEBAYES.out.freebayes_vcf
-        ch_versions     = ch_versions.mix(BAM_VARIANT_CALLING_FREEBAYES.out.versions)
+        freebayes_vcf     = BAM_VARIANT_CALLING_FREEBAYES.out.freebayes_vcf
+        ch_versions       = ch_versions.mix(BAM_VARIANT_CALLING_FREEBAYES.out.versions)
     }
 
     // HAPLOTYPECALLER
@@ -184,8 +182,9 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
                         known_sites_snps_tbi,
                         intervals_bed_combined_haplotypec)
 
-        haplotypecaller_vcf  = BAM_VARIANT_CALLING_HAPLOTYPECALLER.out.filtered_vcf
-        ch_versions          = ch_versions.mix(BAM_VARIANT_CALLING_HAPLOTYPECALLER.out.versions)
+        haplotypecaller_vcf      = BAM_VARIANT_CALLING_HAPLOTYPECALLER.out.filtered_vcf
+        ch_versions              = ch_versions.mix(BAM_VARIANT_CALLING_HAPLOTYPECALLER.out.versions)
+
     }
 
     // MANTA
@@ -197,8 +196,9 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
             fasta_fai
         )
 
-        manta_vcf   = BAM_VARIANT_CALLING_GERMLINE_MANTA.out.manta_vcf
-        ch_versions = ch_versions.mix(BAM_VARIANT_CALLING_GERMLINE_MANTA.out.versions)
+
+        manta_vcf     = BAM_VARIANT_CALLING_GERMLINE_MANTA.out.manta_vcf
+        ch_versions   = ch_versions.mix(BAM_VARIANT_CALLING_GERMLINE_MANTA.out.versions)
     }
 
     // STRELKA
@@ -210,8 +210,8 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
             fasta_fai
         )
 
-        strelka_vcf = BAM_VARIANT_CALLING_SINGLE_STRELKA.out.strelka_vcf
-        ch_versions = ch_versions.mix(BAM_VARIANT_CALLING_SINGLE_STRELKA.out.versions)
+        strelka_vcf     = BAM_VARIANT_CALLING_SINGLE_STRELKA.out.strelka_vcf
+        ch_versions     = ch_versions.mix(BAM_VARIANT_CALLING_SINGLE_STRELKA.out.versions)
     }
 
     //TIDDIT
@@ -222,14 +222,13 @@ workflow BAM_VARIANT_CALLING_GERMLINE_ALL {
             bwa
         )
 
-        tiddit_vcf = BAM_VARIANT_CALLING_SINGLE_TIDDIT.out.tiddit_vcf
-        ch_versions = ch_versions.mix(BAM_VARIANT_CALLING_SINGLE_TIDDIT.out.versions)
+        tiddit_vcf     = BAM_VARIANT_CALLING_SINGLE_TIDDIT.out.tiddit_vcf
+        ch_versions    = ch_versions.mix(BAM_VARIANT_CALLING_SINGLE_TIDDIT.out.versions)
     }
 
     emit:
     deepvariant_vcf
     freebayes_vcf
-    genotype_gvcf
     haplotypecaller_vcf
     manta_vcf
     mpileup_vcf

subworkflows/local/bam_variant_calling_haplotypecaller/main.nf

@@ -134,7 +134,20 @@ workflow BAM_VARIANT_CALLING_HAPLOTYPECALLER {
                     known_sites_indels.concat(known_sites_snps).flatten().unique().collect(),
                     known_sites_indels_tbi.concat(known_sites_snps_tbi).flatten().unique().collect())
 
-        filtered_vcf = VCF_VARIANT_FILTERING_GATK.out.filtered_vcf.map{ meta, vcf-> [[patient:meta.patient, sample:meta.sample, status:meta.status, sex:meta.sex, id:meta.sample, num_intervals:meta.num_intervals, variantcaller:"haplotypecaller"], vcf]}
+        filtered_vcf = VCF_VARIANT_FILTERING_GATK.out.filtered_vcf.map{ meta, vcf-> [
+                [
+                    patient:meta.patient,
+                    sample:meta.sample,
+                    status:meta.status,
+                    sex:meta.sex,
+                    id:meta.sample,
+                    num_intervals:meta.num_intervals,
+                    variantcaller:"haplotypecaller"
+                ],
+                vcf
+            ]
+        }
+
         ch_versions = ch_versions.mix(GATK4_HAPLOTYPECALLER.out.versions)
         ch_versions = ch_versions.mix(MERGE_HAPLOTYPECALLER.out.versions)
         ch_versions = ch_versions.mix(VCF_VARIANT_FILTERING_GATK.out.versions)

subworkflows/local/bam_variant_calling_mpileup/main.nf

@@ -72,5 +72,5 @@ workflow BAM_VARIANT_CALLING_MPILEUP {
     emit:
     versions = ch_versions
     mpileup = Channel.empty().mix(CAT_MPILEUP.out.file_out, mpileup.no_intervals)
-    vcf = Channel.empty().mix(GATK4_MERGEVCFS.out.vcf,vcfs.no_intervals)
+    vcf = Channel.empty().mix(GATK4_MERGEVCFS.out.vcf, vcfs.no_intervals)
 }

subworkflows/local/bam_variant_calling_single_strelka/main.nf

@@ -76,7 +76,7 @@ workflow BAM_VARIANT_CALLING_SINGLE_STRELKA {
                         sex:            meta.sex,
                         status:         meta.status,
                         variantcaller:  "strelka"
-                    ],vcf]
+                    ], vcf]
                 }
 
     ch_versions = ch_versions.mix(MERGE_STRELKA.out.versions)

subworkflows/local/vcf_variant_filtering_gatk/main.nf

@@ -66,4 +66,3 @@ workflow VCF_VARIANT_FILTERING_GATK {
     versions = ch_versions
     filtered_vcf
 }
-

tests/config/pytest_tags.yml

@@ -263,3 +263,20 @@ vep:
   - modules/nf-core/ensemblvep/main.nf
   - modules/nf-core/tabix/bgziptabix/main.nf
   - subworkflows/nf-core/vcf_annotate_ensemblvep/main.nf
+
+## concatenate germline vcfs
+concatenate_vcfs:
+  - conf/modules/post_variant_calling.config
+  - modules/nf-core/deepvariant/main.nf # deepvariant
+  - modules/nf-core/tabix/tabix/main.nf
+  - modules/nf-core/freebayes/main.nf # freebayes
+  - modules/nf-core/gatk4/haplotypecaller/main.nf # haplotypecaller
+  - modules/nf-core/manta/germline/main.nf # manta
+  - modules/nf-core/bcftools/mpileup/main.nf # mpileup/bcftools
+  - modules/nf-core/bcftools/sort/main.nf
+  - modules/nf-core/tabix/bgziptabix/main.nf
+  - modules/nf-core/bcftools/concat/main.nf
+  - modules/nf-core/samtools/mpileup/main.nf
+  - modules/nf-core/gatk4/mergevcfs/main.nf # gatk4/mergevcfs
+  - modules/nf-core/strelka/germline/main.nf # strelka
+  - modules/nf-core/tiddit/sv/main.nf # tiddit

tests/csv/3.0/mapped_joint_bam.fixed.csv

@@ -0,0 +1,3 @@
+patient,status,sample,bam,bai
+testN,0,testN,https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam,https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam.bai
+testT,0,testT,https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/homo_sapiens/illumina/bam/test2.paired_end.sorted.bam,https://raw.githubusercontent.com/nf-core/test-datasets/modules/data/genomics/homo_sapiens/illumina/bam/test2.paired_end.sorted.bam.bai