Add blobtk depth for read coverage

conf/base.config

@@ -133,4 +133,9 @@ process {
         // 100k reads per hour for PacBio, 50m for HiC
         time   = { 1.h  * Math.ceil( meta.read_count / (meta.datatype == "pacbio" ? 100000 : 50000000)) * task.attempt }
     }
+
+    withName: BLOBTK_DEPTH {
+        cpus   = { log_increase_cpus(4, 2*task.attempt, 1, 2) }
+        memory = { 1.GB  * Math.ceil( meta.read_count / 1000000 ) * task.attempt }
+    }
 }

conf/modules.config

@@ -163,7 +163,7 @@ process {
         ext.args   = { (meta.datatype == "pacbio" ? "-y pbccs " : "") + "-O bam" }
     }
 
-    withName: '.*:CONVERT_STATS:SAMTOOLS_.*' {
+    withName: '.*:CONVERT_STATS:SAMTOOLS_.*|BLOBTK_DEPTH' {
         beforeScript = { "export REF_PATH=spoof"}
         publishDir = [
             path: { "${params.outdir}/read_mapping/${meta.datatype}" },
@@ -172,4 +172,8 @@ process {
         ]
     }
 
+    withName: BLOBTK_DEPTH {
+        ext.prefix = { "${cram.baseName}" }
+    }
+
 }

docs/output.md

@@ -89,6 +89,19 @@ The filtered PacBio reads are aligned with `MINIMAP2_ALIGN`. The sorted and merg
 
 If provided using the `--header` option, all output alignments (`*.cram` or `*.bam`) will include any additional metadata supplied as a SAM header template, replacing the existing _@HD_ and _@SD_ entries (note that this behaviour can be altered by modifying the `ext.args` for `SAMTOOLS_REHEADER` in `modules.config`).
 
+### Read coverage
+
+Read coverage of the output alignment file is calculated with [blotk_depth](https://github.com/genomehubs/blobtk/wiki/blobtk-depth).
+
+<details markdown="1">
+<summary>Output files</summary>
+
+- `read_mapping`
+  - `<sequence-type>`
+    - `<gca_accession>.unmasked.<sequence-type>.<sample_id>.[cr|b]am.regions.bed.gz`: Read coverage in BED format
+
+</details>
+
 ### Statistics
 
 The output alignments, along with the index, are used to calculate mapping statistics. Output files are generated using `SAMTOOLS_STATS`, `SAMTOOLS_FLAGSTAT` and `SAMTOOLS_IDXSTATS`.

modules/local/blobtk_depth.nf

@@ -0,0 +1,34 @@
+process BLOBTK_DEPTH {
+    tag "${meta.id}"
+    label 'process_single'
+
+    if (workflow.profile.tokenize(',').intersect(['conda', 'mamba']).size() >= 1) {
+        exit 1, "BLOBTOOLKIT_DEPTH module does not support Conda. Please use Docker / Singularity / Podman instead."
+    }
+    container "docker.io/genomehubs/blobtk:0.6.5"
+
+    input:
+    tuple val(meta), path(cram), path(bai)
+
+    output:
+    tuple val(meta), path('*.regions.bed.gz') , emit: bed
+    path "versions.yml"                       , emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def args = task.ext.args ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    blobtk depth \\
+        -c ${cram} \\
+        $args \\
+        -O ${prefix}.regions.bed.gz \\
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        blobtk: \$(blobtk --version | cut -d' ' -f2)
+    END_VERSIONS
+    """
+}

subworkflows/local/convert_stats.nf

@@ -2,6 +2,7 @@
 // Convert BAM to CRAM, create index and calculate statistics
 //
 
+include { BLOBTK_DEPTH                      } from '../../modules/local/blobtk_depth'
 include { CRUMBLE                           } from '../../modules/nf-core/crumble/main'
 include { SAMTOOLS_VIEW as SAMTOOLS_CRAM    } from '../../modules/nf-core/samtools/view/main'
 include { SAMTOOLS_VIEW as SAMTOOLS_REINDEX } from '../../modules/nf-core/samtools/view/main'
@@ -81,15 +82,20 @@ workflow CONVERT_STATS {
 
         // Set the BAM and BAI channels for emission
         ch_bam = SAMTOOLS_REINDEX.out.bam
-        ch_bai = SAMTOOLS_REINDEX.out.bai
+        ch_csi = SAMTOOLS_REINDEX.out.csi
 
         // If using BAM for stats, use the reindexed BAM
         if ( !("cram" in outfmt_options) ) {
-            ch_data_for_stats = ch_bam.join ( ch_bai )
+            ch_data_for_stats = ch_bam.join ( ch_csi )
         }
     }
 
 
+    // Calculate read depth
+    BLOBTK_DEPTH ( ch_data_for_stats )
+    ch_versions = ch_versions.mix( BLOBTK_DEPTH.out.versions.first() )
+
+
     // Calculate statistics
     SAMTOOLS_STATS ( ch_data_for_stats, fasta )
     ch_versions = ch_versions.mix( SAMTOOLS_STATS.out.versions.first() )