Add coverm/contig to nf-core/modules

CHANGELOG.md

@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### `Added`
 
+- [#997](https://github.com/nf-core/mag/pull/997) - Add CoverM as the default contig depth/coverage calculator, replacing `jgi_summarize_bam_contig_depths` from MetaBAT2. New parameters added: `--depth_calculator` (default: `coverm`), `--coverm_mapper` (default: `bowtie2`), and `--coverm_methods` (default: `metabat`) (by @vinisalazar)
+
 ### `Changed`
 
 - [#1011](https://github.com/nf-core/mag/pull/1011) - Reverted CheckM2 database download workaround from #966 (by @dialvarezs)
@@ -21,6 +23,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### `Dependencies`
 
+| Tool   | Previous version | New version |
+| ------ | ---------------- | ----------- |
+| CoverM |                  | 0.7.0       |
+
 ### `Deprecated`
 
 - [#1018](https://github.com/nf-core/mag/pull/1018) - Remove `mag_depths_plot` local module (by @dialvarezs)

CITATIONS.md

@@ -48,6 +48,10 @@
 
   > De Coster W, D'Hert S, Schultz DT, Cruts M, Van Broeckhoven C. NanoPack: visualizing and processing long-read sequencing data. Bioinformatics. 2018 Aug 1;34(15):2666-2669. doi: 10.1093/bioinformatics/bty149
 
+- [CoverM](https://doi.org/10.1093/bioinformatics/btaf147)
+
+  > Aroney, S. T. N., Newell, R. J. P., Nissen, J. N., Camargo, A. P., Tyson, G. W., & Woodcroft, B. J. (2025). CoverM: read alignment statistics for metagenomics. Bioinformatics, 41(4), btaf147. doi: 10.1093/bioinformatics/btaf147
+
 - [COMEBin](https://doi.org/10.1038/s41467-023-44290-z)
 
   > Wang Z, You R, Han H, Liu W, Sun F, Zhu S. Effective binning of metagenomic contigs using contrastive multi-view representation learning. Nat Commun. 2024 Jan 17;15(1):585. doi: 10.1038/s41467-023-44290-z. PMID: 38233391; PMCID: PMC10794208.

README.md

@@ -45,7 +45,7 @@ The pipeline then:
 - performs assembly using [MEGAHIT](https://github.com/voutcn/megahit) and [SPAdes](http://cab.spbu.ru/software/spades/), and checks their quality using [Quast](http://quast.sourceforge.net/quast) and [ALE](https://github.com/sc932/ALE) (if short read data is used)
 - (optionally) performs ancient DNA assembly validation using [PyDamage](https://github.com/maxibor/pydamage) and contig consensus sequence recalling with [Freebayes](https://github.com/freebayes/freebayes) and [BCFtools](http://samtools.github.io/bcftools/bcftools.html)
 - predicts protein-coding genes for the assemblies using [Prodigal](https://github.com/hyattpd/Prodigal), and bins with [Prokka](https://github.com/tseemann/prokka) and optionally [MetaEuk](https://www.google.com/search?channel=fs&client=ubuntu-sn&q=MetaEuk)
-- performs metagenome binning using [MetaBAT2](https://bitbucket.org/berkeleylab/metabat/src/master/), [MaxBin2](https://sourceforge.net/projects/maxbin2/), [CONCOCT](https://github.com/BinPro/CONCOCT), [COMEBin](https://github.com/ziyewang/COMEBin), [MetaBinner](https://github.com/ziyewang/MetaBinner), and/or [SemiBin2](https://github.com/BigDataBiology/SemiBin)
+- estimates contig coverage using [CoverM](https://github.com/wwood/CoverM) and performs metagenome binning using [MetaBAT2](https://bitbucket.org/berkeleylab/metabat/src/master/), [MaxBin2](https://sourceforge.net/projects/maxbin2/), [CONCOCT](https://github.com/BinPro/CONCOCT), [COMEBin](https://github.com/ziyewang/COMEBin), [MetaBinner](https://github.com/ziyewang/MetaBinner), and/or [SemiBin2](https://github.com/BigDataBiology/SemiBin)
 - checks the quality of the genome bins using [Busco](https://busco.ezlab.org/), [CheckM](https://ecogenomics.github.io/CheckM/), or [CheckM2](https://github.com/chklovski/CheckM2) and optionally [GUNC](https://grp-bork.embl-community.io/gunc/)
 - Performs ancient DNA validation and repair with [pyDamage](https://github.com/maxibor/pydamage) and [freebayes](https://github.com/freebayes/freebayes)
 - optionally refines bins with [DAS Tool](https://github.com/cmks/DAS_Tool)

bin/get_mag_depths.py

@@ -48,7 +48,8 @@ def main(args=None):
     # load contig depths for all samples into dict (could use pandas as well)
     sample_names = []
     dict_contig_depths = {}
-    with gzip.open(args.depths, "rt") as infile:
+    opener = gzip.open if args.depths.endswith('.gz') else open
+    with opener(args.depths, "rt") as infile:
         reader = csv.DictReader(infile, delimiter="\t")
         # process header to extract sample names from column names
         depth_columns = []

conf/modules.config

@@ -340,25 +340,6 @@ process {
         ]
     }
 
-    withName: MINIMAP2_ASSEMBLY_ALIGN {
-        ext.args   = {
-            meta.lr_platform == "OXFORD_NANOPORE"
-                ? "-x map-ont"
-                : meta.lr_platform == "OXFORD_NANOPORE_HQ"
-                    ? "-x lr:hq"
-                    : meta.lr_platform == "PACBIO_CLR"
-                        ? "-x map-pb"
-                        : meta.lr_platform == "PACBIO_HIFI" ? "-x map-hifi" : ""
-        }
-        ext.prefix = { "${meta.assembler}-${meta.id}-${meta2.id}" }
-        publishDir = [
-            path: { "${params.outdir}/Assembly/${meta2.assembler}/QC/${meta2.id}" },
-            mode: params.publish_dir_mode,
-            pattern: "*.{bam,bai}",
-            enabled: params.save_assembly_mapped_reads,
-        ]
-    }
-
     withName: SAMTOOLS_HOSTREMOVED_UNMAPPED {
         ext.args   = '-f 4'
         ext.prefix = { "${meta.id}_${meta.run}_hostremoved" }
@@ -808,23 +789,72 @@ process {
         ]
     }
 
-    withName: METABAT2_JGISUMMARIZEBAMCONTIGDEPTHS_SHORTREAD {
+    withName: 'COVERM_CONTIG_SHORTREAD' {
         ext.prefix = { "${meta.assembler}-${meta.id}-depth" }
-        ext.args   = params.shortread_percentidentity ? "--percentIdentity ${params.shortread_percentidentity}" : ''
+        ext.args   = [
+            '--methods metabat',
+            params.shortread_coverage_mapper != 'bowtie2' ? "--mapper ${params.shortread_coverage_mapper}" : '',
+            params.shortread_percentidentity ? "--min-read-percent-identity ${params.shortread_percentidentity}" : '',
+        ].join(' ').trim()
         publishDir = [
-            path: { "${params.outdir}/GenomeBinning/depths/contigs" },
-            mode: params.publish_dir_mode,
-            pattern: '*-depth.txt.gz',
+            [
+                path: { "${params.outdir}/GenomeBinning/depths/contigs" },
+                mode: params.publish_dir_mode,
+                pattern: '*.depth.tsv',
+            ],
+            [
+                path: { "${params.outdir}/Assembly/${meta.assembler}/QC/${meta.id}" },
+                mode: params.publish_dir_mode,
+                pattern: "*.bam",
+                enabled: params.save_assembly_mapped_reads,
+            ],
         ]
     }
 
-    withName: METABAT2_JGISUMMARIZEBAMCONTIGDEPTHS_LONGREAD {
+    withName: 'COVERM_CONTIG_LONGREAD' {
         ext.prefix = { "${meta.assembler}-${meta.id}-depth" }
-        ext.args   = params.longread_percentidentity ? "--percentIdentity ${params.longread_percentidentity}" : ''
+        ext.args   = {
+            def mapper = meta.lr_platform in ['OXFORD_NANOPORE', 'OXFORD_NANOPORE_HQ']
+                ? 'minimap2-ont'
+                : meta.lr_platform == 'PACBIO_CLR'
+                    ? 'minimap2-pb'
+                    : meta.lr_platform == 'PACBIO_HIFI' ? 'minimap2-hifi' : ''
+            [
+                '--methods metabat',
+                mapper ? "--mapper ${mapper}" : '',
+                params.longread_percentidentity ? "--min-read-percent-identity ${params.longread_percentidentity}" : '',
+            ].join(' ').trim()
+        }
+        publishDir = [
+            [
+                path: { "${params.outdir}/GenomeBinning/depths/contigs" },
+                mode: params.publish_dir_mode,
+                pattern: '*.depth.tsv',
+            ],
+            [
+                path: { "${params.outdir}/Assembly/${meta.assembler}/QC/${meta.id}" },
+                mode: params.publish_dir_mode,
+                pattern: "*.bam",
+                enabled: params.save_assembly_mapped_reads,
+            ],
+        ]
+    }
+
+    withName: 'SAMTOOLS_COVERM_SHORTREAD_INDEX' {
+        publishDir = [
+            path: { "${params.outdir}/Assembly/${meta.assembler}/QC/${meta.id}" },
+            mode: params.publish_dir_mode,
+            pattern: "*.bai",
+            enabled: params.save_assembly_mapped_reads,
+        ]
+    }
+
+    withName: 'SAMTOOLS_COVERM_LONGREAD_INDEX' {
         publishDir = [
-            path: { "${params.outdir}/GenomeBinning/depths/contigs" },
+            path: { "${params.outdir}/Assembly/${meta.assembler}/QC/${meta.id}" },
             mode: params.publish_dir_mode,
-            pattern: '*-depth.txt.gz',
+            pattern: "*.bai",
+            enabled: params.save_assembly_mapped_reads,
         ]
     }
 

conf/test_coverm_shortread.config

@@ -0,0 +1,26 @@
+/*
+========================================================================================
+    Nextflow config file for testing CoverM native short-read coverage mapping
+========================================================================================
+*/
+
+includeConfig 'test.config'
+
+params {
+    config_profile_name        = 'Test CoverM short-read profile'
+    config_profile_description = 'Minimal test dataset using CoverM native short-read mapping'
+
+    shortread_coverage_mapper  = 'minimap2-sr'
+    save_assembly_mapped_reads = true
+    skip_clipping              = true
+    keep_phix                  = true
+    skip_semibin               = true
+    skip_prokka                = true
+    skip_quast                 = true
+    skip_binqc                 = true
+    run_busco                  = false
+    run_checkm                 = false
+    run_checkm2                = false
+    skip_gtdbtk                = true
+    cat_db                     = null
+}

conf/test_full.config

@@ -98,11 +98,11 @@ process {
         memory        = { 2.GB * task.attempt }
         errorStrategy = { task.exitStatus in [1, 255] ? 'ignore' : 'retry' }
     }
-    withName: 'NFCORE_MAG:MAG:BINNING:METABAT2_JGISUMMARIZEBAMCONTIGDEPTHS_LONGREAD' {
+    withName: 'NFCORE_MAG:MAG:BINNING_PREPARATION:LONGREAD_BINNING_PREPARATION:COVERM_CONTIG_LONGREAD' {
         cpus   = { 3 * task.attempt }
         memory = { 1.GB * task.attempt }
     }
-    withName: 'NFCORE_MAG:MAG:BINNING:METABAT2_JGISUMMARIZEBAMCONTIGDEPTHS_SHORTREAD' {
+    withName: 'NFCORE_MAG:MAG:BINNING_PREPARATION:SHORTREAD_BINNING_PREPARATION:COVERM_CONTIG_SHORTREAD' {
         cpus   = { 3 * task.attempt }
         memory = { 1.GB * task.attempt }
     }
@@ -118,11 +118,7 @@ process {
         cpus   = { 18 * task.attempt }
         memory = { 1.GB * task.attempt }
     }
-    withName: 'NFCORE_MAG:MAG:BINNING_PREPARATION:LONGREAD_BINNING_PREPARATION:MINIMAP2_ASSEMBLY_ALIGN' {
-        cpus   = { 13 * task.attempt }
-        memory = { 17.GB * task.attempt }
-    }
-    withName: 'NFCORE_MAG:MAG:BINNING_PREPARATION:LONGREAD_BINNING_PREPARATION:MINIMAP2_ASSEMBLY_INDEX' {
+    withName: 'NFCORE_MAG:MAG:BINNING_PREPARATION:LONGREAD_BINNING_PREPARATION:SAMTOOLS_COVERM_LONGREAD_INDEX' {
         cpus   = { 1 * task.attempt }
         memory = { 2.GB * task.attempt }
     }

docs/output.md

@@ -361,15 +361,17 @@ Protein-coding genes are predicted for each assembly.
 
 ### Contig sequencing depth
 
-Sequencing depth per contig and sample is generated by MetaBAT2's `jgi_summarize_bam_contig_depths --outputDepth`. The values correspond to `(sum of exactly aligned bases) / ((contig length)-2*75)`. For example, for two reads aligned exactly with `10` and `9` bases on a 1000 bp long contig the depth is calculated by `(10+9)/(1000-2*75)` (1000bp length of contig minus 75bp from each end, which is excluded).
+Sequencing depth per contig and sample is generated by [CoverM](https://github.com/wwood/CoverM) using the `metabat` method (equivalent to MetaBAT2's `jgi_summarize_bam_contig_depths --outputDepth`). The values correspond to `(sum of exactly aligned bases) / ((contig length)-2*75)`. For example, for two reads aligned exactly with `10` and `9` bases on a 1000 bp long contig the depth is calculated by `(10+9)/(1000-2*75)` (1000bp length of contig minus 75bp from each end, which is excluded).
+
+The legacy `jgi_summarize_bam_contig_depths` tool from MetaBAT2 can still be used by setting `--depth_calculator metabat2`.
 
 These depth files are used for downstream binning steps.
 
 <details markdown="1">
 <summary>Output files</summary>
 
 - `GenomeBinning/depths/contigs/`
-  - `[assembler]-[sample/group]-depth.txt.gz`: Sequencing depth for each contig and sample or group, only for short reads.
+  - `[assembler]-[sample/group]-depth.tsv`: Sequencing depth for each contig and sample or group.
 
 </details>
 

docs/usage.md

@@ -437,12 +437,24 @@ This may result in a different set or none of contigs being evaluated in pyDamag
 
 In order to run the binning tools included in the pipeline, MAG must first align reads back to the assemblies, and estimate the coverage of each contig.
 
-During the coverage estimation step, these alignments are by default filtered to retain alignments that have a percentage identity of 97% (i.e., of the base pairs that match between the read and the contig, 97% are identical). This value is a good default for short read Illumina data, however for certain long read technologies, the error rates in the reads can be much higher.
-For example, older Oxford Nanopore chemistries can have error rates approaching
-15% - 20%.
+The pipeline uses [CoverM](https://github.com/wwood/CoverM) for coverage estimation. CoverM outputs a MetaBAT2-compatible depth file (using `--methods metabat`) which is compatible with all downstream binning tools.
+
+### Mapper selection
+
+The `--shortread_coverage_mapper` parameter controls which aligner is used for mapping short reads back to assembled contigs:
+
+- `bowtie2` (default): Runs Bowtie2 manually to produce BAM files, which are then passed to CoverM.
+- CoverM-native short-read mappers (`bwa-mem`, `bwa-mem2`, `minimap2-sr`, `strobealign`): CoverM handles both mapping and coverage calculation directly from reads and exports cached BAM files for downstream binners.
+
+Long-read coverage mapping is handled by CoverM automatically from the samplesheet `long_reads_platform` metadata, using `minimap2-ont` for Oxford Nanopore reads, `minimap2-pb` for PacBio CLR reads, and `minimap2-hifi` for PacBio HiFi reads.
+
+### Alignment identity filtering
+
+During the coverage estimation step, alignments are by default filtered to retain alignments that have a percentage identity of 97% (i.e., of the base pairs that match between the read and the contig, 97% are identical). This value is a good default for short read Illumina data, however for certain long read technologies, the error rates in the reads can be much higher.
+For example, older Oxford Nanopore chemistries can have error rates approaching 15% - 20%.
 
 If you are having trouble with the coverage estimation steps (for example, the output depths for each bin are all at or near zero), it may be worth manually adjusting this parameter, if it is appropriate for your data.
-You can do this by adjusting the `--longread_percentidentity` and `--shortread_percentidentity` parameters for long reads and short reads, respectively.
+You can do this by adjusting the `--longread_percentidentity` and `--shortread_percentidentity` parameters for long reads and short reads, respectively. These set CoverM's `--min-read-percent-identity` flag.
 For older ONT data, you may wish to look at values of around 85% to improve coverage estimation.
 
 ## A note on bin refinement

modules.json

@@ -131,6 +131,11 @@
                         "git_sha": "e753770db613ce014b3c4bc94f6cba443427b726",
                         "installed_by": ["fasta_binning_concoct"]
                     },
+                    "coverm/contig": {
+                        "branch": "master",
+                        "git_sha": "38c52b07d406618f6ba4bb21e80172e0230c5623",
+                        "installed_by": ["modules"]
+                    },
                     "dastool/dastool": {
                         "branch": "master",
                         "git_sha": "e753770db613ce014b3c4bc94f6cba443427b726",

modules/local/mag_depths_convert/main.nf

@@ -15,21 +15,25 @@ process CONVERT_DEPTHS {
     path "versions.yml"                                   , emit: versions
 
     script:
+    def depth_unzipped = depth.toString() - '.gz'
     """
-    gunzip -f ${depth}
+    # Decompress only if gzipped
+    if [[ "${depth}" == *.gz ]]; then
+        gunzip -f "${depth}"
+    fi
 
     # Determine the number of abundance columns
-    n_abund=\$(awk 'NR==1 {print int((NF-3)/2)}' ${depth.toString() - '.gz'})
+    n_abund=\$(awk 'NR==1 {print int((NF-3)/2)}' ${depth_unzipped})
 
     # Get column names
-    read -r header<${depth.toString() - '.gz'}
+    read -r header<${depth_unzipped}
     header=(\$header)
 
     # Generate abundance files for each read set
     for i in \$(seq 1 \$n_abund); do
         col=\$((i*2+2))
-        name=\$( echo \${header[\$col-1]} | sed s/\\.bam\$// )
-        bioawk -t '{if (NR > 1) {print \$1, \$'"\$col"'}}' ${depth.toString() - '.gz'} > \${name}.abund
+        name=\$( basename "\${header[\$col-1]}" | sed s/\\.bam\$// )
+        bioawk -t '{if (NR > 1) {print \$1, \$'"\$col"'}}' ${depth_unzipped} > \${name}.abund
     done
 
     cat <<-END_VERSIONS > versions.yml