Add subworkflow for gridss somatic structural variant calling

subworkflows/nf-core/bam_somatic_sv_calling_gridss/main.nf

@@ -0,0 +1,96 @@
+include { SAMTOOLS_INDEX as SAMTOOLS_INDEX_NORMAL                                         } from '../../../modules/nf-core/samtools/index/main'
+include { SAMTOOLS_INDEX as SAMTOOLS_INDEX_TUMOR                                          } from '../../../modules/nf-core/samtools/index/main'
+include { GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS as GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS_NORMAL } from '../../../modules/nf-core/gridss/extractoverlappingfragments/main'
+include { GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS as GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS_TUMOR  } from '../../../modules/nf-core/gridss/extractoverlappingfragments/main'
+include { GRIDSS_PREPROCESS as GRIDSS_PREPROCESS_NORMAL                                   } from '../../../modules/nf-core/gridss/preprocess/main'
+include { GRIDSS_PREPROCESS as GRIDSS_PREPROCESS_TUMOR                                    } from '../../../modules/nf-core/gridss/preprocess/main'
+include { GRIDSS_ASSEMBLE                                                                 } from '../../../modules/nf-core/gridss/assemble/main'
+include { GRIDSS_CALL                                                                     } from '../../../modules/nf-core/gridss/call/main'
+include { GRIDSS_SOMATICFILTER                                                            } from '../../../modules/nf-core/gridss/somaticfilter/main'
+
+
+workflow BAM_SOMATIC_SV_CALLING_GRIDSS {
+
+    take:
+    ch_input_bam            // channel: mandatory [ val(meta), [ tumor_bam, normal_bam ], [tumor_bai, normal_bai] ]
+    ch_fasta_fai_bwaindex   // channel: mandatory [ val(meta2), fasta, fai, bwa_index ]
+    ch_pondir               // channel: mandatory [ val(meta3), [ pondir ] ]
+    ch_gridss_config        // channel: optional  [ val(meta4), [ gridss_config ] ]
+    ch_target_bed           // channel: optional  [ val(meta5), [ target_bed ] ]
+    is_targeted             // boolean: mandatory true for panel/exome data, false to run untargeted
+
+    main:
+
+    def ch_input_normal = ch_input_bam
+        .map { meta, bams, bais ->
+            tuple(meta, bams[1], bais[1])
+        }
+    def ch_input_tumor  = ch_input_bam
+        .map { meta, bams, bais ->
+            tuple(meta, bams[0], bais[0])
+        }
+
+    // Targeted (panel/exome) data is subset to the fragments overlapping the target
+    // regions before preprocessing, untargeted data goes straight to preprocessing
+    if (is_targeted) {
+
+        GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS_NORMAL(
+            ch_input_normal,
+            ch_target_bed
+        )
+        GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS_TUMOR(
+            ch_input_tumor,
+            ch_target_bed
+        )
+
+        SAMTOOLS_INDEX_NORMAL ( GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS_NORMAL.out.bam )
+        SAMTOOLS_INDEX_TUMOR  ( GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS_TUMOR.out.bam )
+
+        ch_input_normal = GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS_NORMAL.out.bam
+            .join(SAMTOOLS_INDEX_NORMAL.out.index)
+        ch_input_tumor  = GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS_TUMOR.out.bam
+            .join(SAMTOOLS_INDEX_TUMOR.out.index)
+    }
+
+    GRIDSS_PREPROCESS_NORMAL(
+        ch_input_normal,
+        ch_fasta_fai_bwaindex
+    )
+
+    GRIDSS_PREPROCESS_TUMOR(
+        ch_input_tumor,
+        ch_fasta_fai_bwaindex
+    )
+
+    def ch_input_assemble = ch_input_tumor
+        .join(ch_input_normal)
+        .join(GRIDSS_PREPROCESS_TUMOR.out.preprocess_dir)
+        .join(GRIDSS_PREPROCESS_NORMAL.out.preprocess_dir)
+        .map { meta, tbam, tbai, nbam, nbai, preprocess_tumor, preprocess_normal ->
+            tuple(meta, [tbam, nbam], [tbai, nbai], [preprocess_tumor, preprocess_normal])
+        }
+
+    GRIDSS_ASSEMBLE(
+        ch_input_assemble,
+        ch_fasta_fai_bwaindex,
+        ch_gridss_config
+    )
+
+    def ch_input_call = ch_input_assemble
+        .join(GRIDSS_ASSEMBLE.out.assemble_dir)
+
+    GRIDSS_CALL (
+        ch_input_call,
+        ch_fasta_fai_bwaindex,
+        ch_gridss_config
+    )
+
+    GRIDSS_SOMATICFILTER ( GRIDSS_CALL.out.vcf, ch_pondir )
+
+
+    emit:
+    gridss_vcf                  = GRIDSS_CALL.out.vcf                      // channel: [ val(meta), join_call_vcf ]
+    all_somatic_vcf             = GRIDSS_SOMATICFILTER.out.all_sv          // channel: [ val(meta), somatic_vcf ]
+    high_confidence_somatic_vcf = GRIDSS_SOMATICFILTER.out.high_conf_sv    // channel: [ val(meta), high_conf_sv_vcf ]
+
+}

subworkflows/nf-core/bam_somatic_sv_calling_gridss/meta.yml

@@ -0,0 +1,94 @@
+# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/subworkflows/yaml-schema.json
+name: "bam_somatic_sv_calling_gridss"
+description: |
+  Call somatic structural variants from a tumor/normal BAM pair with GRIDSS.
+  Runs the full GRIDSS workflow (optional targeted fragment extraction, preprocess,
+  joint assembly, joint calling) followed by gridss_somatic_filter to separate
+  somatic events from germline and panel-of-normals artefacts.
+keywords:
+  - gridss
+  - structural variants
+  - somatic
+  - tumor
+  - normal
+  - bam
+  - vcf
+components:
+  - gridss/extractoverlappingfragments
+  - gridss/preprocess
+  - gridss/assemble
+  - gridss/call
+  - gridss/somaticfilter
+  - samtools/index
+input:
+  - ch_input_bam:
+      type: file
+      description: |
+        Channel containing the tumor/normal BAM pair and their indices.
+        The tumor entry must be first and the normal entry second in each list;
+        this is the sample order used for the joint GRIDSS call, so
+        gridss_somatic_filter must be given `--tumourordinal 1 --normalordinal 2`
+        via `ext.args` on GRIDSS_SOMATICFILTER.
+        Structure: [ val(meta), [ path(tumor_bam), path(normal_bam) ], [ path(tumor_bai), path(normal_bai) ] ]
+      pattern: "*.{bam,bam.bai}"
+  - ch_fasta_fai_bwaindex:
+      type: file
+      description: |
+        Channel containing the reference genome, its samtools index and the BWA index
+        directory. All BWA index files must share the basename of the FASTA.
+        Structure: [ val(meta2), path(fasta), path(fasta_fai), path(bwa_index) ]
+      pattern: "*.{fa,fna,fasta,fai}"
+  - ch_pondir:
+      type: directory
+      description: |
+        Directory containing the panel-of-normals BED/BEDPE files used by
+        gridss_somatic_filter to remove false positive somatic events.
+        Generate it with gridss.GeneratePonBedpe. Pass [ [], [] ] to filter without a PON.
+        Structure: [ val(meta3), path(pondir) ]
+  - ch_gridss_config:
+      type: file
+      description: |
+        Optional GRIDSS configuration file (Java .properties format) applied to the
+        assemble and call steps. Pass [ [], [] ] to use the GRIDSS defaults.
+        Structure: [ val(meta4), path(gridss_config) ]
+      pattern: "*.properties"
+  - ch_target_bed:
+      type: file
+      description: |
+        BED file of target regions for panel/exome data. Only used when `is_targeted`
+        is true, in which case reads overlapping the targets are extracted with
+        gridss_extract_overlapping_fragments and re-indexed before preprocessing.
+        Pass [ [], [] ] when running untargeted.
+        Structure: [ val(meta5), path(target_bed) ]
+      pattern: "*.bed"
+  - is_targeted:
+      type: boolean
+      description: |
+        Set to true for panel/exome data to subset the BAM files to the fragments
+        overlapping `ch_target_bed` before preprocessing, false to run GRIDSS on the
+        full BAM files. When true, `ext.prefix` must disambiguate the tumor and normal
+        outputs of GRIDSS_EXTRACTOVERLAPPINGFRAGMENTS, which otherwise both default to
+        `${meta.id}.subset`.
+output:
+  - gridss_vcf:
+      type: file
+      description: |
+        Unfiltered joint (tumor + normal) structural variant calls from GRIDSS.
+        Structure: [ val(meta), path(vcf) ]
+      pattern: "*.sv.gridss.vcf.gz"
+  - all_somatic_vcf:
+      type: file
+      description: |
+        All somatic structural variants retained by gridss_somatic_filter.
+        Structure: [ val(meta), path(vcf) ]
+      pattern: "*.all_somatic.vcf.bgz"
+  - high_confidence_somatic_vcf:
+      type: file
+      description: |
+        High confidence subset of the somatic structural variant calls.
+        Structure: [ val(meta), path(vcf) ]
+      pattern: "*.high_confidence_somatic.vcf.bgz"
+authors:
+  - "@imsarath"
+maintainers:
+  - "@imsarath"

subworkflows/nf-core/bam_somatic_sv_calling_gridss/tests/main.nf.test

@@ -0,0 +1,168 @@
+nextflow_workflow {
+
+    name "Test Subworkflow BAM_SOMATIC_SV_CALLING_GRIDSS"
+    script "../main.nf"
+    workflow "BAM_SOMATIC_SV_CALLING_GRIDSS"
+    config "./nextflow.config"
+
+    tag "subworkflows"
+    tag "subworkflows_nfcore"
+    tag "subworkflows/bam_somatic_sv_calling_gridss"
+    tag "bwa/index"
+    tag "gridss"
+    tag "gridss/extractoverlappingfragments"
+    tag "gridss/preprocess"
+    tag "gridss/assemble"
+    tag "gridss/call"
+    tag "gridss/somaticfilter"
+    tag "samtools"
+    tag "samtools/index"
+
+    setup {
+        run("BWA_INDEX") {
+            script "../../../../modules/nf-core/bwa/index/main.nf"
+            process {
+                """
+                input[0] = [
+                    [ id:'genome' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true)
+                ]
+                """
+            }
+        }
+    }
+
+    test("human - bam - tumor_normal - untargeted") {
+
+        when {
+            workflow {
+                """
+                input[0] = channel.of([
+                    [ id:'test' ],
+                    [
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test2.paired_end.sorted.bam', checkIfExists: true),
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam', checkIfExists: true)
+                    ],
+                    [
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test2.paired_end.sorted.bam.bai', checkIfExists: true),
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam.bai', checkIfExists: true)
+                    ]
+                ])
+                input[1] = channel.of([
+                    [ id:'genome' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true),
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true)
+                ]).join(BWA_INDEX.out.index)
+                input[2] = [ [], [] ]
+                input[3] = [
+                    [ id:'gridss_config' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/gridss/gridss.properties', checkIfExists: true)
+                ]
+                input[4] = [ [], [] ]
+                input[5] = false
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert workflow.success },
+                { assert snapshot(
+                    path(workflow.out.gridss_vcf[0][1]).vcf.summary,
+                    workflow.out.all_somatic_vcf.collect { meta, vcf -> [ meta, file(vcf).name ] },
+                    workflow.out.high_confidence_somatic_vcf.collect { meta, vcf -> [ meta, file(vcf).name ] }
+                ).match() }
+            )
+        }
+    }
+
+    test("human - bam - tumor_normal - targeted") {
+
+        when {
+            workflow {
+                """
+                input[0] = channel.of([
+                    [ id:'test' ],
+                    [
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test2.paired_end.sorted.bam', checkIfExists: true),
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam', checkIfExists: true)
+                    ],
+                    [
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test2.paired_end.sorted.bam.bai', checkIfExists: true),
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam.bai', checkIfExists: true)
+                    ]
+                ])
+                input[1] = channel.of([
+                    [ id:'genome' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true),
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true)
+                ]).join(BWA_INDEX.out.index)
+                input[2] = [ [], [] ]
+                input[3] = [
+                    [ id:'gridss_config' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/gridss/gridss.properties', checkIfExists: true)
+                ]
+                input[4] = [
+                    [ id:'target_bed' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.bed', checkIfExists: true)
+                ]
+                input[5] = true
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert workflow.success },
+                { assert snapshot(
+                    path(workflow.out.gridss_vcf[0][1]).vcf.summary,
+                    workflow.out.all_somatic_vcf.collect { meta, vcf -> [ meta, file(vcf).name ] },
+                    workflow.out.high_confidence_somatic_vcf.collect { meta, vcf -> [ meta, file(vcf).name ] }
+                ).match() }
+            )
+        }
+    }
+
+    test("human - bam - tumor_normal - stub") {
+
+        options "-stub"
+
+        when {
+            workflow {
+                """
+                input[0] = channel.of([
+                    [ id:'test' ],
+                    [
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test2.paired_end.sorted.bam', checkIfExists: true),
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam', checkIfExists: true)
+                    ],
+                    [
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test2.paired_end.sorted.bam.bai', checkIfExists: true),
+                        file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/bam/test.paired_end.sorted.bam.bai', checkIfExists: true)
+                    ]
+                ])
+                input[1] = channel.of([
+                    [ id:'genome' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta', checkIfExists: true),
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/genome/genome.fasta.fai', checkIfExists: true)
+                ]).join(BWA_INDEX.out.index)
+                input[2] = [ [], [] ]
+                input[3] = [
+                    [ id:'gridss_config' ],
+                    file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/gridss/gridss.properties', checkIfExists: true)
+                ]
+                input[4] = [ [], [] ]
+                input[5] = false
+                """
+            }
+        }
+
+        then {
+            assertAll(
+                { assert workflow.success },
+                { assert snapshot(sanitizeOutput(workflow.out)).match() }
+            )
+        }
+    }
+
+}