Add applyvqsr (#1101)

* initial commit to setup branch

* workflow finished

* Update nextflow.config

* tumour to tumor, getpileup passed as nomral and tumor

* paired_somatic renamed to tumor_normal_somatic

* Apply suggestions from code review

Co-authored-by: Maxime U. Garcia <maxime.garcia@scilifelab.se>

* Update subworkflows/nf-core/gatk_tumor_normal_somatic_variant_calling/main.nf

Co-authored-by: Maxime U. Garcia <maxime.garcia@scilifelab.se>

* updated index names in meta.yml

* changed index file names in main script and test

* Apply suggestions from code review

Co-authored-by: Maxime U. Garcia <maxime.garcia@scilifelab.se>

* Apply suggestions from code review

* fixed bug from changes

* Apply suggestions from code review

* modified yml to allow new subworkflow testing

* Update test.yml

* Update test.yml

* add applyvqsr

* added memory options, new test data used

* Update main.nf

* Update main.nf

Co-authored-by: GCJMackenzie <gavin.mackenzie@nibsc.org>
Co-authored-by: Maxime U. Garcia <maxime.garcia@scilifelab.se>

modules/gatk4/applyvqsr/main.nf

@@ -0,0 +1,55 @@
+process GATK4_APPLYVQSR {
+    tag "$meta.id"
+    label 'process_low'
+
+    conda (params.enable_conda ? "bioconda::gatk4=4.2.3.0" : null)
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/gatk4:4.2.3.0--hdfd78af_0' :
+        'quay.io/biocontainers/gatk4:4.2.3.0--hdfd78af_0' }"
+
+    input:
+    tuple val(meta), path(vcf), path(tbi), path(recal), path(recalidx), path(tranches)
+    path fasta
+    path fai
+    path dict
+    val allelespecific
+    val truthsensitivity
+    val mode
+
+    output:
+    tuple val(meta), path("*.vcf.gz")     , emit: vcf
+    tuple val(meta), path("*.tbi")        , emit: tbi
+    path "versions.yml"                   , emit: versions
+
+    script:
+    def args = task.ext.args ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    refCommand = fasta ? "-R ${fasta} " : ''
+    alleleSpecificCommand = allelespecific ? '-AS' : ''
+    truthSensitivityCommand = truthsensitivity ? "--truth-sensitivity-filter-level ${truthsensitivity}" : ''
+    modeCommand = mode ? "--mode ${mode} " : 'SNP'
+
+    def avail_mem = 3
+    if (!task.memory) {
+        log.info '[GATK ApplyVQSR] Available memory not known - defaulting to 3GB. Specify process memory requirements to change this.'
+    } else {
+        avail_mem = task.memory.giga
+    }
+    """
+    gatk --java-options "-Xmx${avail_mem}g" ApplyVQSR \\
+        ${refCommand} \\
+        -V ${vcf} \\
+        -O ${prefix}.vcf.gz \\
+        ${alleleSpecificCommand} \\
+        ${truthSensitivityCommand} \\
+        --tranches-file $tranches \\
+        --recal-file $recal \\
+        ${modeCommand} \\
+        $args
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        gatk4: \$(echo \$(gatk --version 2>&1) | sed 's/^.*(GATK) v//; s/ .*\$//')
+    END_VERSIONS
+    """
+}

modules/gatk4/applyvqsr/meta.yml

@@ -0,0 +1,88 @@
+name: gatk4_applyvqsr
+description: |
+  Apply a score cutoff to filter variants based on a recalibration table.
+  AplyVQSR performs the second pass in a two-stage process called Variant Quality Score Recalibration (VQSR).
+  Specifically, it applies filtering to the input variants based on the recalibration table produced
+  in the first step by VariantRecalibrator and a target sensitivity value.
+keywords:
+  - gatk4
+  - applyvqsr
+  - VQSR
+tools:
+  - gatk4:
+      description: |
+        Developed in the Data Sciences Platform at the Broad Institute, the toolkit offers a wide variety of tools
+        with a primary focus on variant discovery and genotyping. Its powerful processing engine
+        and high-performance computing features make it capable of taking on projects of any size.
+      homepage: https://gatk.broadinstitute.org/hc/en-us
+      documentation: https://gatk.broadinstitute.org/hc/en-us/categories/360002369672s
+      doi: 10.1158/1538-7445.AM2017-3590
+      licence: ['Apache-2.0']
+
+input:
+  - meta:
+      type: map
+      description: |
+        Groovy Map containing sample information
+        e.g. [ id:'test']
+  - vcf:
+      type: file
+      description: VCF file to be recalibrated, this should be the same file as used for the first stage VariantRecalibrator.
+      pattern: "*.vcf"
+  - tbi:
+      type: file
+      description: Tbi index for the input vcf file.
+      pattern: "*.vcf.tbi"
+  - recal:
+      type: file
+      description: Recalibration file produced when the input vcf was run through VariantRecalibrator in stage 1.
+      pattern: "*.recal"
+  - recalidx:
+      type: file
+      description: Index file for the recalibration file.
+      pattern: ".recal.idx"
+  - tranches:
+      type: boolean
+      description: Tranches file produced when the input vcf was run through VariantRecalibrator in stage 1.
+      pattern: ".tranches"
+  - fasta:
+      type: file
+      description: The reference fasta file
+      pattern: "*.fasta"
+  - fai:
+      type: file
+      description: Index of reference fasta file
+      pattern: "*.fasta.fai"
+  - dict:
+      type: file
+      description: GATK sequence dictionary
+      pattern: "*.dict"
+  - allelespecific:
+      type: boolean
+      description: Whether or not to run ApplyVQSR in allele specific mode, this should be kept the same as the stage 1 VariantRecalibrator run.
+      pattern: "{true,false}"
+  - truthsensitivity:
+      type: double
+      description: Value to be used as the truth sensitivity cutoff score.
+      pattern: "99.0"
+  - mode:
+      type: String
+      description: Specifies which recalibration mode to employ, should be the same as the stage 1 VariantRecalibrator run. (SNP is default, BOTH is intended for testing only)
+      pattern: "{SNP,INDEL,BOTH}"
+
+output:
+  - vcf:
+      type: file
+      description: compressed vcf file containing the recalibrated variants.
+      pattern: "*.vcf.gz"
+  - tbi:
+      type: file
+      description: Index of recalibrated vcf file.
+      pattern: "*vcf.gz.tbi"
+  - versions:
+      type: file
+      description: File containing software versions.
+      pattern: "versions.yml"
+
+authors:
+  - "@GCJMackenzie"

tests/config/pytest_modules.yml

@@ -504,6 +504,10 @@ gatk4/applybqsr:
   - modules/gatk4/applybqsr/**
   - tests/modules/gatk4/applybqsr/**
 
+gatk4/applyvqsr:
+  - modules/gatk4/applyvqsr/**
+  - tests/modules/gatk4/applyvqsr/**
+
 gatk4/baserecalibrator:
   - modules/gatk4/baserecalibrator/**
   - tests/modules/gatk4/baserecalibrator/**

tests/config/test_data.config

@@ -115,6 +115,9 @@ params {
                 transcriptome_fasta                            = "${test_data_dir}/genomics/homo_sapiens/genome/transcriptome.fasta"
                 genome2_fasta                                  = "${test_data_dir}/genomics/homo_sapiens/genome/genome2.fasta"
 		        genome_chain_gz                                = "${test_data_dir}/genomics/homo_sapiens/genome/genome.chain.gz"
+                genome_21_fasta                                = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/sequence/genome.fasta"
+                genome_21_fasta_fai                            = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/sequence/genome.fasta.fai"
+                genome_21_dict                                 = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/sequence/genome.dict"
 
                 dbsnp_146_hg38_vcf_gz                          = "${test_data_dir}/genomics/homo_sapiens/genome/vcf/dbsnp_146.hg38.vcf.gz"
                 dbsnp_146_hg38_vcf_gz_tbi                      = "${test_data_dir}/genomics/homo_sapiens/genome/vcf/dbsnp_146.hg38.vcf.gz.tbi"
@@ -122,6 +125,14 @@ params {
                 gnomad_r2_1_1_vcf_gz_tbi                       = "${test_data_dir}/genomics/homo_sapiens/genome/vcf/gnomAD.r2.1.1.vcf.gz.tbi"
                 mills_and_1000g_indels_vcf_gz                  = "${test_data_dir}/genomics/homo_sapiens/genome/vcf/mills_and_1000G.indels.vcf.gz"
                 mills_and_1000g_indels_vcf_gz_tbi              = "${test_data_dir}/genomics/homo_sapiens/genome/vcf/mills_and_1000G.indels.vcf.gz.tbi"
+                hapmap_3_3_hg38_21_vcf_gz                      = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/germlineresources/hapmap_3.3.hg38.vcf.gz"
+                hapmap_3_3_hg38_21_vcf_gz_tbi                  = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/germlineresources/hapmap_3.3.hg38.vcf.gz.tbi"
+                res_1000g_omni2_5_hg38_21_vcf_gz               = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/germlineresources/1000G_omni2.5.hg38.vcf.gz"
+                res_1000g_omni2_5_hg38_21_vcf_gz_tbi           = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/germlineresources/1000G_omni2.5.hg38.vcf.gz.tbi"
+                res_1000g_phase1_snps_hg38_21_vcf_gz           = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/germlineresources/1000G_phase1.snps.hg38.vcf.gz"
+                res_1000g_phase1_snps_hg38_21_vcf_gz_tbi       = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/germlineresources/1000G_phase1.snps.hg38.vcf.gz.tbi"
+                dbsnp_138_hg38_21_vcf_gz                       = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/germlineresources/dbsnp_138.hg38.vcf.gz"
+                dbsnp_138_hg38_21_vcf_gz_tbi                   = "${test_data_dir}/genomics/homo_sapiens/genome/chr21/germlineresources/dbsnp_138.hg38.vcf.gz.tbi"
 
                 syntheticvcf_short_vcf_gz                      = "${test_data_dir}/genomics/homo_sapiens/genome/vcf/syntheticvcf_short.vcf.gz"
                 syntheticvcf_short_vcf_gz_tbi                  = "${test_data_dir}/genomics/homo_sapiens/genome/vcf/syntheticvcf_short.vcf.gz.tbi"
@@ -197,6 +208,16 @@ params {
                 test_genomicsdb_tar_gz                         = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/test_genomicsdb.tar.gz"
                 test_pon_genomicsdb_tar_gz                     = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/test_pon_genomicsdb.tar.gz"
 
+                test2_haplotc_ann_vcf_gz                       = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/haplotypecaller_calls/test2_haplotc.ann.vcf.gz"
+                test2_haplotc_ann_vcf_gz_tbi                   = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/haplotypecaller_calls/test2_haplotc.ann.vcf.gz.tbi"
+
+                test2_recal                                    = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/variantrecalibrator/test2.recal"
+                test2_recal_idx                                = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/variantrecalibrator/test2.recal.idx"
+                test2_tranches                                 = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/variantrecalibrator/test2.tranches"
+                test2_allele_specific_recal                    = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/variantrecalibrator/test2_allele_specific.recal"
+                test2_allele_specific_recal_idx                = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/variantrecalibrator/test2_allele_specific.recal.idx"
+                test2_allele_specific_tranches                 = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/variantrecalibrator/test2_allele_specific.tranches"
+
                 test_test2_paired_mutect2_calls_vcf_gz         = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/paired_mutect2_calls/test_test2_paired_mutect2_calls.vcf.gz"
                 test_test2_paired_mutect2_calls_vcf_gz_tbi     = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/paired_mutect2_calls/test_test2_paired_mutect2_calls.vcf.gz.tbi"
                 test_test2_paired_mutect2_calls_vcf_gz_stats   = "${test_data_dir}/genomics/homo_sapiens/illumina/gatk/paired_mutect2_calls/test_test2_paired_mutect2_calls.vcf.gz.stats"

tests/modules/gatk4/applyvqsr/main.nf

@@ -0,0 +1,41 @@
+#!/usr/bin/env nextflow
+
+nextflow.enable.dsl = 2
+
+include { GATK4_APPLYVQSR } from '../../../../modules/gatk4/applyvqsr/main.nf'
+
+workflow test_gatk4_applyvqsr {
+    input = [ [ id:'test'], // meta map
+              file(params.test_data['homo_sapiens']['illumina']['test2_haplotc_ann_vcf_gz'], checkIfExists: true),
+              file(params.test_data['homo_sapiens']['illumina']['test2_haplotc_ann_vcf_gz_tbi'], checkIfExists: true),
+              file(params.test_data['homo_sapiens']['illumina']['test2_recal'], checkIfExists: true),
+              file(params.test_data['homo_sapiens']['illumina']['test2_recal_idx'], checkIfExists: true),
+              file(params.test_data['homo_sapiens']['illumina']['test2_tranches'], checkIfExists: true)
+            ]
+    fasta = file(params.test_data['homo_sapiens']['genome']['genome_21_fasta'], checkIfExists: true)
+    fai = file(params.test_data['homo_sapiens']['genome']['genome_21_fasta_fai'], checkIfExists: true)
+    dict = file(params.test_data['homo_sapiens']['genome']['genome_21_dict'], checkIfExists: true)
+    allelespecific = false
+    truthsensitivity = '99.0'
+    mode = 'SNP'
+
+    GATK4_APPLYVQSR ( input, fasta, fai, dict, allelespecific, truthsensitivity, mode )
+}
+
+workflow test_gatk4_applyvqsr_allele_specific {
+    input = [ [ id:'test'], // meta map
+              file(params.test_data['homo_sapiens']['illumina']['test2_haplotc_ann_vcf_gz'], checkIfExists: true),
+              file(params.test_data['homo_sapiens']['illumina']['test2_haplotc_ann_vcf_gz_tbi'], checkIfExists: true),
+              file(params.test_data['homo_sapiens']['illumina']['test2_allele_specific_recal'], checkIfExists: true),
+              file(params.test_data['homo_sapiens']['illumina']['test2_allele_specific_recal_idx'], checkIfExists: true),
+              file(params.test_data['homo_sapiens']['illumina']['test2_allele_specific_tranches'], checkIfExists: true)
+            ]
+    fasta = file(params.test_data['homo_sapiens']['genome']['genome_21_fasta'], checkIfExists: true)
+    fai = file(params.test_data['homo_sapiens']['genome']['genome_21_fasta_fai'], checkIfExists: true)
+    dict = file(params.test_data['homo_sapiens']['genome']['genome_21_dict'], checkIfExists: true)
+    allelespecific = true
+    truthsensitivity = '99.0'
+    mode = 'SNP'
+
+    GATK4_APPLYVQSR ( input, fasta, fai, dict, allelespecific, truthsensitivity, mode )
+}

tests/modules/gatk4/applyvqsr/nextflow.config

@@ -0,0 +1,5 @@
+process {
+
+    publishDir = { "${params.outdir}/${task.process.tokenize(':')[-1].tokenize('_')[0].toLowerCase()}" }
+
+}

tests/modules/gatk4/applyvqsr/test.yml

@@ -0,0 +1,17 @@
+- name: gatk4 applyvqsr test_gatk4_applyvqsr
+  command: nextflow run tests/modules/gatk4/applyvqsr -entry test_gatk4_applyvqsr -c tests/config/nextflow.config -c ./tests/modules/gatk4/applyvqsr/nextflow.config
+  tags:
+    - gatk4
+    - gatk4/applyvqsr
+  files:
+    - path: output/gatk4/test.vcf.gz
+    - path: output/gatk4/test.vcf.gz.tbi
+
+- name: gatk4 applyvqsr test_gatk4_applyvqsr_allele_specific
+  command: nextflow run tests/modules/gatk4/applyvqsr -entry test_gatk4_applyvqsr_allele_specific -c tests/config/nextflow.config -c ./tests/modules/gatk4/applyvqsr/nextflow.config
+  tags:
+    - gatk4
+    - gatk4/applyvqsr
+  files:
+    - path: output/gatk4/test.vcf.gz
+    - path: output/gatk4/test.vcf.gz.tbi