Mutect2 add mitochondria mode and update tests (#967)

* new mitochondria mode added, tests updated to allow for temp fix for test data

* add cram test

* bam/bam_idx renamed to input and input_index

Co-authored-by: GCJMackenzie <gavin.mackenzie@nibsc.org>

modules/gatk4/mutect2/main.nf

@@ -19,9 +19,11 @@ process GATK4_MUTECT2 {
     }
 
     input:
-    tuple val(meta) , path(bam) , path(bai) , val(which_norm)
-    val run_single
-    val run_pon
+    tuple val(meta) , path(input) , path(input_index) , val(which_norm)
+    val  run_single
+    val  run_pon
+    val  run_mito
+    val  interval_label
     path fasta
     path fastaidx
     path dict
@@ -39,35 +41,34 @@ process GATK4_MUTECT2 {
 
     script:
     def prefix   = options.suffix ? "${meta.id}${options.suffix}" : "${meta.id}"
-    def inputsList = []
-    def normalsList = []
-    def inputsCommand = ''
-    def panelsCommand = ''
-    def normalsCommand = ''
+    def panels_command = ''
+    def normals_command = ''
 
-    bam.each() {a -> inputsList.add(" -I " + a ) }
-    inputsCommand = inputsList.join( ' ')
+    def inputs_command = '-I ' + input.join( ' -I ')
 
     if(run_pon) {
-        panelsCommand = ''
-        normalsCommand = ''
+        panels_command = ''
+        normals_command = ''
 
     } else if(run_single) {
-        panelsCommand = " --germline-resource $germline_resource --panel-of-normals $panel_of_normals"
-        normalsCommand = ''
+        panels_command = " --germline-resource $germline_resource --panel-of-normals $panel_of_normals"
+        normals_command = ''
+
+    } else if(run_mito){
+        panels_command = "-L ${interval_label} --mitochondria-mode"
+        normals_command = ''
 
     } else {
-        panelsCommand = " --germline-resource $germline_resource --panel-of-normals $panel_of_normals --f1r2-tar-gz ${prefix}.f1r2.tar.gz"
-        which_norm.each() {a -> normalsList.add(" -normal " + a ) }
-        normalsCommand = normalsList.join( ' ')
+        panels_command = " --germline-resource $germline_resource --panel-of-normals $panel_of_normals --f1r2-tar-gz ${prefix}.f1r2.tar.gz"
+        normals_command = '-normal ' + which_norm.join( ' -normal ')
     }
 
     """
     gatk Mutect2 \\
         -R ${fasta} \\
-        ${inputsCommand} \\
-        ${normalsCommand} \\
-        ${panelsCommand} \\
+        ${inputs_command} \\
+        ${normals_command} \\
+        ${panels_command} \\
         -O ${prefix}.vcf.gz \\
         $options.args
 

modules/gatk4/mutect2/meta.yml

@@ -22,23 +22,34 @@ input:
       description: |
         Groovy Map containing sample information
         e.g. [ id:'test']
-  - bam:
+  - input:
       type: list
-      description: list of BAM files
-      pattern: "*.bam"
-  - bai:
+      description: list of BAM files, also able to take CRAM as an input
+      pattern: "*.{bam/cram}"
+  - input_index:
       type: list
-      description: list of BAM file indexes
-      pattern: "*.bam.bai"
+      description: list of BAM file indexes, also able to take CRAM indexes as an input
+      pattern: "*.{bam.bai/cram.crai}"
   - which_norm:
       type: list
       description: optional list of sample headers contained in the normal sample bam files (these are required for tumor_normal_pair mode)
+      pattern: "testN"
   - run_single:
       type: boolean
       description: Specify whether or not to run in tumor_single mode instead of tumor_normal_pair mode (will be ignored if run_pon is also true)
+      pattern: "true/false"
   - run_pon:
       type: boolean
       description: Specify whether or not to run in panel_of_normal mode instead of tumor_normal_pair mode
+      pattern: "true/false"
+  - run_mito:
+      type: boolean
+      description: Specify whether or not to run in mitochondria-mode instead of tumor_normal_pair mode
+      pattern: "true/false"
+  - interval_label:
+      type: string
+      description: Specify the label used for mitochondrial chromosome when mutect2 is run in mitochondria mode.
+      pattern: "chrM"
   - fasta:
       type: file
       description: The reference fasta file

tests/modules/gatk4/mutect2/main.nf

@@ -3,6 +3,8 @@
 nextflow.enable.dsl = 2
 
 include { GATK4_MUTECT2 } from '../../../../modules/gatk4/mutect2/main.nf' addParams( options: [:] )
+// used to run with the mitochondria mode setting as this increases sensitivity, allowing for some tumor_normal variants to be detected while the old test data is still in use, will be removed when new test data for sarek is available.
+include { GATK4_MUTECT2 as GATK4_TEMPFIX_MUTECT2 } from '../../../../modules/gatk4/mutect2/main.nf' addParams( options: [args: '--mitochondria-mode'] )
 
 workflow test_gatk4_mutect2_tumor_normal_pair {
     input = [ [ id:'test'], // meta map
@@ -12,6 +14,8 @@ workflow test_gatk4_mutect2_tumor_normal_pair {
             ]
     run_single = false
     run_pon = false
+    run_mito = false
+    interval_label = []
     fasta = file(params.test_data['homo_sapiens']['genome']['genome_fasta'], checkIfExists: true)
     fastaidx = file(params.test_data['homo_sapiens']['genome']['genome_fasta_fai'], checkIfExists: true)
     dict = file(params.test_data['homo_sapiens']['genome']['genome_dict'], checkIfExists: true)
@@ -20,7 +24,7 @@ workflow test_gatk4_mutect2_tumor_normal_pair {
     panel_of_normals = file(params.test_data['homo_sapiens']['genome']['mills_and_1000g_indels_vcf_gz'], checkIfExists: true)
     panel_of_normals_idx = file(params.test_data['homo_sapiens']['genome']['mills_and_1000g_indels_vcf_gz_tbi'], checkIfExists: true)
 
-    GATK4_MUTECT2 ( input , run_single , run_pon , fasta , fastaidx , dict , germline_resource, germline_resource_idx , panel_of_normals , panel_of_normals_idx )
+    GATK4_TEMPFIX_MUTECT2 ( input , run_single , run_pon , run_mito , interval_label , fasta , fastaidx , dict , germline_resource, germline_resource_idx , panel_of_normals , panel_of_normals_idx )
 }
 
 workflow test_gatk4_mutect2_tumor_single {
@@ -31,6 +35,8 @@ workflow test_gatk4_mutect2_tumor_single {
             ]
     run_single = true
     run_pon = false
+    run_mito = false
+    interval_label = []
     fasta = file(params.test_data['homo_sapiens']['genome']['genome_fasta'], checkIfExists: true)
     fastaidx = file(params.test_data['homo_sapiens']['genome']['genome_fasta_fai'], checkIfExists: true)
     dict = file(params.test_data['homo_sapiens']['genome']['genome_dict'], checkIfExists: true)
@@ -39,7 +45,28 @@ workflow test_gatk4_mutect2_tumor_single {
     panel_of_normals = file(params.test_data['homo_sapiens']['genome']['mills_and_1000g_indels_vcf_gz'], checkIfExists: true)
     panel_of_normals_idx = file(params.test_data['homo_sapiens']['genome']['mills_and_1000g_indels_vcf_gz_tbi'], checkIfExists: true)
 
-    GATK4_MUTECT2 ( input , run_single , run_pon , fasta , fastaidx , dict , germline_resource, germline_resource_idx , panel_of_normals , panel_of_normals_idx )
+    GATK4_MUTECT2 ( input , run_single , run_pon , run_mito , interval_label , fasta , fastaidx , dict , germline_resource, germline_resource_idx , panel_of_normals , panel_of_normals_idx )
+}
+
+workflow test_gatk4_mutect2_cram_input {
+    input = [ [ id:'test'], // meta map
+              [ file(params.test_data['homo_sapiens']['illumina']['test2_paired_end_recalibrated_sorted_cram'], checkIfExists: true)],
+              [ file(params.test_data['homo_sapiens']['illumina']['test2_paired_end_recalibrated_sorted_cram_crai'], checkIfExists: true)],
+              []
+            ]
+    run_single = true
+    run_pon = false
+    run_mito = false
+    interval_label = []
+    fasta = file(params.test_data['homo_sapiens']['genome']['genome_fasta'], checkIfExists: true)
+    fastaidx = file(params.test_data['homo_sapiens']['genome']['genome_fasta_fai'], checkIfExists: true)
+    dict = file(params.test_data['homo_sapiens']['genome']['genome_dict'], checkIfExists: true)
+    germline_resource = file(params.test_data['homo_sapiens']['genome']['gnomad_r2_1_1_vcf_gz'], checkIfExists: true)
+    germline_resource_idx = file(params.test_data['homo_sapiens']['genome']['gnomad_r2_1_1_vcf_gz_tbi'], checkIfExists: true)
+    panel_of_normals = file(params.test_data['homo_sapiens']['genome']['mills_and_1000g_indels_vcf_gz'], checkIfExists: true)
+    panel_of_normals_idx = file(params.test_data['homo_sapiens']['genome']['mills_and_1000g_indels_vcf_gz_tbi'], checkIfExists: true)
+
+    GATK4_MUTECT2 ( input , run_single , run_pon , run_mito , interval_label , fasta , fastaidx , dict , germline_resource, germline_resource_idx , panel_of_normals , panel_of_normals_idx )
 }
 
 workflow test_gatk4_mutect2_generate_pon {
@@ -50,6 +77,8 @@ workflow test_gatk4_mutect2_generate_pon {
             ]
     run_single = false
     run_pon = true
+    run_mito = false
+    interval_label = []
     fasta = file(params.test_data['homo_sapiens']['genome']['genome_fasta'], checkIfExists: true)
     fastaidx = file(params.test_data['homo_sapiens']['genome']['genome_fasta_fai'], checkIfExists: true)
     dict = file(params.test_data['homo_sapiens']['genome']['genome_dict'], checkIfExists: true)
@@ -58,5 +87,27 @@ workflow test_gatk4_mutect2_generate_pon {
     panel_of_normals = []
     panel_of_normals_idx = []
 
-    GATK4_MUTECT2 ( input , run_single , run_pon , fasta , fastaidx , dict , germline_resource, germline_resource_idx , panel_of_normals , panel_of_normals_idx )
+    GATK4_MUTECT2 ( input , run_single , run_pon, run_mito , interval_label , fasta , fastaidx , dict , germline_resource, germline_resource_idx , panel_of_normals , panel_of_normals_idx )
+}
+
+// mitochondria mode would ideally have some mitochondria test data, but since the mitochondria settings only increase detection sensitivity, we can use the chr22 data as a stand in as it is already a small dataset, the extra variants detected compared to generate_pon shows the mode is working.
+workflow test_gatk4_mutect2_mitochondria {
+    input = [ [ id:'test'], // meta map
+              [ file(params.test_data['homo_sapiens']['illumina']['test_paired_end_recalibrated_sorted_bam'], checkIfExists: true)],
+              [ file(params.test_data['homo_sapiens']['illumina']['test_paired_end_recalibrated_sorted_bam_bai'], checkIfExists: true)],
+              []
+            ]
+    run_single = false
+    run_pon = false
+    run_mito = true
+    interval_label = 'chr22'
+    fasta = file(params.test_data['homo_sapiens']['genome']['genome_fasta'], checkIfExists: true)
+    fastaidx = file(params.test_data['homo_sapiens']['genome']['genome_fasta_fai'], checkIfExists: true)
+    dict = file(params.test_data['homo_sapiens']['genome']['genome_dict'], checkIfExists: true)
+    germline_resource = []
+    germline_resource_idx = []
+    panel_of_normals = []
+    panel_of_normals_idx = []
+
+    GATK4_MUTECT2 ( input , run_single , run_pon, run_mito , interval_label , fasta , fastaidx , dict , germline_resource, germline_resource_idx , panel_of_normals , panel_of_normals_idx )
 }

tests/modules/gatk4/mutect2/test.yml

@@ -7,7 +7,7 @@
     - path: output/gatk4/test.f1r2.tar.gz
     - path: output/gatk4/test.vcf.gz
     - path: output/gatk4/test.vcf.gz.stats
-      md5sum: 6ecb874e6a95aa48233587b876c2a7a9
+      md5sum: 887d54e393510f1d0aa2c33bc6155161
     - path: output/gatk4/test.vcf.gz.tbi
 
 - name: gatk4 mutect2 test_gatk4_mutect2_tumor_single
@@ -18,7 +18,18 @@
   files:
     - path: output/gatk4/test.vcf.gz
     - path: output/gatk4/test.vcf.gz.stats
-      md5sum: e7ef613f7d158b8a0adf44abe5db2029
+      md5sum: 106c5828b02b906c97922618b6072169
+    - path: output/gatk4/test.vcf.gz.tbi
+
+- name: gatk4 mutect2 test_gatk4_mutect2_cram_input
+  command: nextflow run tests/modules/gatk4/mutect2 -entry test_gatk4_mutect2_cram_input -c tests/config/nextflow.config
+  tags:
+    - gatk4
+    - gatk4/mutect2
+  files:
+    - path: output/gatk4/test.vcf.gz
+    - path: output/gatk4/test.vcf.gz.stats
+      md5sum: 106c5828b02b906c97922618b6072169
     - path: output/gatk4/test.vcf.gz.tbi
 
 - name: gatk4 mutect2 test_gatk4_mutect2_generate_pon
@@ -31,3 +42,14 @@
     - path: output/gatk4/test.vcf.gz.stats
       md5sum: 4f77301a125913170b8e9e7828b4ca3f
     - path: output/gatk4/test.vcf.gz.tbi
+
+- name: gatk4 mutect2 test_gatk4_mutect2_mitochondria
+  command: nextflow run tests/modules/gatk4/mutect2 -entry test_gatk4_mutect2_mitochondria -c tests/config/nextflow.config
+  tags:
+    - gatk4
+    - gatk4/mutect2
+  files:
+    - path: output/gatk4/test.vcf.gz
+    - path: output/gatk4/test.vcf.gz.stats
+      md5sum: fc6ea14ca2da346babe78161beea28c9
+    - path: output/gatk4/test.vcf.gz.tbi