Changes:

- Add nf-test skeleton code - Change Metator input to BAM files - Update conda-checking code in bin3c modules - Add paramters to choose minimum contig size and minimum map % identity when binning
sanger-tol · Dec 20, 2024 · f6db65a · f6db65a
1 parent 9693dd7
commit f6db65a
Show file tree

Hide file tree

Showing 22 changed files with 281 additions and 91 deletions.
diff --git a/.gitignore b/.gitignore
@@ -8,4 +8,5 @@ testing*
 *.pyc
 null/
 co2footprint*
-
+.nf-test/
+.nf-test.log
diff --git a/assets/test_full_input.yaml b/assets/test_full_input.yaml
diff --git a/assets/test_input.yaml b/assets/test_input.yaml
diff --git a/bin/bin_summary.R b/bin/bin_summary.R
@@ -156,11 +156,3 @@ summary <- map(input_types, \(x) split_and_read(input, x)) |>
     reduce(\(x, y) left_join(x, y, by = "bin"))
 
 write_tsv(summary, glue::glue("{input$prefix}.bin_summary.tsv"))
-
-writeLines(
-    c("BIN_SUMMARY:",
-        paste0("    R: ", paste0(R.Version()[c("major","minor")], collapse = ".")),
-        paste0("    tidyverse: ", packageVersion("tidyverse"))
-    ),
-    "versions.yml"
-)
diff --git a/conf/modules.config b/conf/modules.config
@@ -20,8 +20,9 @@ process {
     ]
 
     withName: 'BIN3C_MKMAP' {
-        ext.prefix = { "${meta.id}_${meta.assembler}_bin3c" }
-        tag        = { "${meta.id}_${meta.assembler}" }
+        ext.args   = { "--min-reflen ${params.minimum_contig_size}" }
+        ext.prefix = { "${meta.id}_${meta.assembler}_bin3c"         }
+        tag        = { "${meta.id}_${meta.assembler}"               }
     }
 
     withName: 'BIN3C_CLUSTER' {
@@ -76,7 +77,7 @@ process {
     }
 
     withName: 'COVERM_CONTIG' {
-        ext.args   = { "-m metabat" }
+        ext.args   = { "-m metabat --min-read-percent-identity ${params.min_hifi_perc_identity}" }
         ext.prefix = { "${meta.id}_${meta.assembler}_depth" }
         tag        = { "${meta.id}_${meta.assembler}" }
         publishDir = [
@@ -178,7 +179,7 @@ process {
     }
 
     withName: 'MAXBIN2' {
-        ext.args   = { "" }
+        ext.args   = { "-min_contig_length ${params.minimum_contig_size}" }
         ext.prefix = { "${meta.id}_${meta.assembler}_maxbin2" }
         tag        = { "${meta.id}_${meta.assembler}" }
         publishDir = [
@@ -189,7 +190,7 @@ process {
     }
 
     withName: 'METABAT2_METABAT2' {
-        ext.args   = { "" }
+        ext.args   = { "--minContig ${params.minimum_contig_size}" }
         ext.prefix = { "${meta.id}_${meta.assembler}_metabat2" }
         tag        = { "${meta.id}_${meta.assembler}" }
         publishDir = [
@@ -199,8 +200,13 @@ process {
         ]
     }
 
+    withName: 'METATOR_PROCESS_INPUT_BAM' {
+        ext.prefix = { "${meta.id}_${meta.assembler}" } // metator already appends "_metator" to files
+        tag        = { "${meta.id}_${meta.assembler}" }
+    }
+
     withName: 'METATOR_PIPELINE' {
-        ext.args   = { "--start fastq --aligner bwa" }
+        ext.args   = { "--start bam" }
         ext.prefix = { "${meta.id}_${meta.assembler}" } // metator already appends "_metator" to files
         tag        = { "${meta.id}_${meta.assembler}" }
         publishDir = [
@@ -211,7 +217,7 @@ process {
     }
 
     withName: 'MINIMAP2_ALIGN' {
-        ext.args   = { "" }
+        ext.args   = { "-x map-hifi" }
         ext.args2  = { "" }
         ext.prefix = { "${meta2.id}_${meta2.assembler}" }
         tag        = { "${meta2.id}_${meta2.assembler}" }
@@ -255,7 +261,7 @@ process {
         ]
     }
 
-    withName: 'SORT_HIC_BAM' {
+    withName: 'SAMTOOLS_SORT_HIC_BAM' {
         ext.args   = { "-n" }
         ext.prefix = { "${meta.id}_${meta.assembler}" }
         tag        = { "${meta.id}_${meta.assembler}" }

diff --git a/conf/test.config b/conf/test.config
@@ -23,6 +23,6 @@ params {
     config_profile_description = 'Minimal test dataset to check pipeline function'
 
     // Input data
-    input         = "${baseDir}/assets/test_input.yaml"
+    input         = params.pipelines_testdata_base_path + "test_input.yaml"
     enable_gtdbtk = false
 }
diff --git a/modules/local/bin3c/cluster/main.nf b/modules/local/bin3c/cluster/main.nf
@@ -13,6 +13,9 @@ process BIN3C_CLUSTER {
     path("versions.yml")                  , emit: versions
 
     script:
+    if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) {
+        error "ERROR: Bin3C is only avaliable as a Docker or Singularity container. If you need to run with conda, run with --enable_bin3c false"
+    }
     def args   = task.ext.args   ?: ''
     def prefix = task.ext.prefix ?: "${meta.id}"
     def contigs_basename = contigs.getBaseName()

diff --git a/modules/local/bin3c/mkmap/main.nf b/modules/local/bin3c/mkmap/main.nf
@@ -13,6 +13,9 @@ process BIN3C_MKMAP {
     path("versions.yml")           , emit: versions
 
     script:
+    if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) {
+        error "ERROR: Bin3C is only avaliable as a Docker or Singularity container. If you need to run with conda, run with --enable_bin3c false"
+    }
     def args   = task.ext.args   ?: ''
     def prefix = task.ext.prefix ?: "${meta.id}"
     // if(!hic_enzymes) error("Error: no enzymes entry found in Hi-C meta object!")

diff --git a/modules/local/bin_summary/main.nf b/modules/local/bin_summary/main.nf
@@ -31,5 +31,23 @@ process BIN_SUMMARY {
         ${checkm_input} \\
         ${tax_input} \\
         ${args}
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        R: \$(Rscript -e 'cat(paste0(R.Version()[c("major","minor")], collapse = "."))')
+        tidyverse: \$(Rscript -e 'cat(packageVersion("tidyverse"))')
+    END_VERSIONS
+    """
+
+    stub:
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    touch ${prefix}.bin_summary.tsv
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        R: \$(Rscript -e 'cat(paste0(R.Version()[c("major","minor")], collapse = "."))')
+        tidyverse: \$(Rscript -e 'cat(packageVersion("tidyverse"))')
+    END_VERSIONS
     """
 }
diff --git a/modules/local/metator/process_input_bam/environment.yml b/modules/local/metator/process_input_bam/environment.yml
@@ -0,0 +1,7 @@
+channels:
+  - conda-forge
+  - bioconda
+dependencies:
+  - bioconda::htslib=1.21
+  - bioconda::samtools=1.21
+  - bioconda::bioawk=1.0
diff --git a/modules/local/metator/process_input_bam/main.nf b/modules/local/metator/process_input_bam/main.nf
@@ -0,0 +1,64 @@
+process METATOR_PROCESS_INPUT_BAM {
+    tag "${meta.id}"
+    label 'process_low'
+
+    conda "${moduleDir}/environment.yml"
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'oras://community.wave.seqera.io/library/htslib_samtools_bioawk:3ff2c81f84424e4c' :
+        'community.wave.seqera.io/library/htslib_samtools_bioawk:420f5543dfc64992' }"
+
+    input:
+    tuple val(meta), path(bam), val(direction)
+
+    output:
+    tuple val(meta), path("*.bam"), emit: filtered_bam
+    path "versions.yml"           , emit: versions
+
+    script:
+    def args   = task.ext.args   ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    switch(direction) {
+        case "fwd": flag = "0x40"; break
+        case "rev": flag = "0x80" ; break
+        default:
+            error("ERROR: METATOR_PROCESS_INPUT_BAM direction was not 'fwd' or 'rev'!")
+            break
+    }
+    """
+    samtools view --threads ${task.cpus-1} -f ${flag} -o temp.bam ${bam}
+    samtools view -H --threads ${task.cpus-1} temp.bam > temp_header
+
+    samtools view --threads ${task.cpus-1} temp.bam |\\
+        bioawk -c sam '{ \$flag = and(\$flag , 3860 ) ; print \$0 }' |\\
+        cat temp_header - |\\
+        samtools sort --threads ${task.cpus-1} -n -o ${prefix}.${direction}.bam
+
+    rm temp.bam temp_header
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//')
+        bioawk: 1.0
+    END_VERSIONS
+    """
+
+    stub:
+    def args   = task.ext.args   ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    switch(direction) {
+        case "fwd": flag = "0x40"; break
+        case "rev": flag = "0x80" ; break
+        default:
+            error("ERROR: METATOR_PROCESS_INPUT_BAM direction was not 'fwd' or 'rev'!")
+            break
+    }
+    """
+    touch ${prefix}.${direction}.bam
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//')
+        bioawk: 1.0
+    END_VERSIONS
+    """
+}
diff --git a/modules/nf-core/prokka/main.nf b/modules/nf-core/prokka/main.nf
diff --git a/nextflow.config b/nextflow.config
@@ -22,6 +22,8 @@ params {
     enable_maxbin2             = true
     enable_bin3c               = true
     enable_metator             = true
+    minimum_contig_size        = 3000
+    min_hifi_perc_identity     = 97
 
     // Bin refinement options
     enable_bin_refinement      = true
@@ -74,7 +76,7 @@ params {
     help_full                    = false
     show_hidden                  = false
     version                      = false
-    pipelines_testdata_base_path = 'https://raw.githubusercontent.com/nf-core/test-datasets/'
+    pipelines_testdata_base_path = 'https://raw.githubusercontent.com/prototaxites/lrm_testdata/refs/heads/main/'
 
     // Config options
     config_profile_name        = null

diff --git a/nextflow_schema.json b/nextflow_schema.json
@@ -83,6 +83,16 @@
                     "type": "boolean",
                     "default": true,
                     "description": "Enable binning with MetaTor."
+                },
+                "minimum_contig_size": {
+                    "type": "integer",
+                    "default": 3000,
+                    "description": "The minimum size of contig to be considered when binning. This option applies to MetaBat2, MaxBin2, and Bin3C, but not Metator."
+                },
+                "min_hifi_perc_identity": {
+                    "type": "integer",
+                    "default": 97,
+                    "description": "The minimum required percent identity of a read mapping to the assembly to be counted during coverage estimation for MetaBat2 and MaxBin2."
                 }
             }
         },
@@ -134,8 +144,8 @@
                     "description": "Enable QC using CheckM2."
                 },
                 "checkm2_db_version": {
-                    "type": "number",
-                    "default": 5571251.0,
+                    "type": "string",
+                    "default": "5571251",
                     "description": "If no local CheckM2 database provided, the Zenodo record ID of a CheckM2 database to download."
                 },
                 "checkm2_local_db": {
@@ -330,7 +340,7 @@
                     "type": "string",
                     "fa_icon": "far fa-check-circle",
                     "description": "Base URL or local path to location of pipeline test dataset files",
-                    "default": "https://raw.githubusercontent.com/nf-core/test-datasets/",
+                    "default": "https://raw.githubusercontent.com/prototaxites/lrm_testdata/refs/heads/main/",
                     "hidden": true
                 }
             }

diff --git a/nf-test.config b/nf-test.config
@@ -0,0 +1,10 @@
+config {
+    testsDir "."
+    workDir ".nf-test"
+    configFile "tests/nextflow.config"
+    profile "test"
+
+    plugins {
+        load "[email protected]"
+    }
+}
diff --git a/subworkflows/local/binning.nf b/subworkflows/local/binning.nf
@@ -5,6 +5,7 @@ include { MAXBIN2                     } from '../../modules/nf-core/maxbin2/main
 include { GAWK as GAWK_MAXBIN2_DEPTHS } from '../../modules/nf-core/gawk/main'
 include { METABAT2_METABAT2           } from '../../modules/nf-core/metabat2/metabat2/main'
 include { METATOR_PIPELINE            } from '../../modules/local/metator/pipeline/main'
+include { METATOR_PROCESS_INPUT_BAM   } from '../../modules/local/metator/process_input_bam/main'
 
 workflow BINNING {
     take:
@@ -49,11 +50,7 @@ workflow BINNING {
         ch_bins = ch_bins.mix(ch_maxbin2_bins)
     }
 
-    //
-    // LOGIC: Bin3C is not available in conda
-    //        only run if we are not running with the conda profile
-    //
-    if(params.enable_bin3c && !(workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1)) {
+    if(params.enable_bin3c) {
         ch_bin3c_mkmap_input = assemblies
             | combine(hic_bam, by: 0)
 
@@ -71,18 +68,39 @@ workflow BINNING {
         ch_bins = ch_bins.mix(ch_bin3c_bins)
     }
 
-    // NOTE: currently hi-c input to metator is supplied as fastq rather than BAM
-    // metator aligns fwd and rev reads independently
-    // TODO: process mapped bam from READ_MAPPING subworkflow into separate fwd/rev
-    // bam files
     if(params.enable_metator) {
+        // Old code (keep for now) - run Metator in mapping mode,
+        // providing hi-c reads as input rather than bam files
+        /*
         ch_assemblies_combine = assemblies
             | map {meta, contigs -> [ meta.subMap('id'), meta, contigs ]}
 
         ch_metator_input = ch_assemblies_combine
             | combine(hic_reads, by: 0)
             | map { meta_join, meta_assembly, contigs, hic -> [meta_assembly, contigs, hic, []]}
 
+        METATOR_PIPELINE(ch_metator_input, hic_enzymes)
+        ch_versions = ch_versions.mix(METATOR_PIPELINE.out.versions)
+        */
+
+        // Metator expects us to have aligned forward and reverse reads
+        // independently of one another - munge the bam file
+        // to filter out forward and reverse reads and remove mate information
+        // from SAM flags: bitwise and(flag, 3860)
+        ch_directions = Channel.of("fwd", "rev")
+        ch_hic_bam_to_process = hic_bam
+            | combine(ch_directions)
+
+        METATOR_PROCESS_INPUT_BAM(ch_hic_bam_to_process)
+        ch_versions = ch_versions.mix(METATOR_PROCESS_INPUT_BAM.out.versions)
+
+        ch_metator_input = METATOR_PROCESS_INPUT_BAM.out.filtered_bam
+            | groupTuple(by: 0, size: 2)
+            | combine(assemblies, by: 0)
+            | map { meta, bams, contigs ->
+                [ meta, contigs, bams.sort(), [] ]
+            }
+
         METATOR_PIPELINE(ch_metator_input, hic_enzymes)
         ch_versions = ch_versions.mix(METATOR_PIPELINE.out.versions)