Fix mapAD multiref mapping

jch-13 · jch-13 · commit 8c34dae7c872 · 2024-06-07T10:45:26.000+02:00
diff --git a/conf/modules.config b/conf/modules.config
@@ -586,6 +586,30 @@ process {
         ]
     }
 
+    withName: ".*FASTQ_ALIGN_MAPAD:BAM_SORT_STATS_SAMTOOLS:SAMTOOLS_SORT" {
+        tag = { "${meta.reference}:${meta.genomic_region}|${meta.sample_id}_${meta.library_id}" }
+        ext.prefix = { "${meta.sample_id}_${meta.library_id}_${meta.reference}_${meta.genomic_region}_mapAD" }
+        publishDir = [
+            enabled: false
+        ]
+    }
+
+    withName: ".*FASTQ_ALIGN_MAPAD:BAM_SORT_STATS_SAMTOOLS:SAMTOOLS_INDEX" {
+        tag = { "${meta.reference}:${meta.genomic_region}|${meta.sample_id}_${meta.library_id}" }
+        ext.prefix = { "${meta.sample_id}_${meta.library_id}_${meta.reference}_${meta.genomic_region}_mapAD" }
+        publishDir = [
+            enabled: false
+        ]
+    }
+
+    withName: ".*FASTQ_ALIGN_MAPAD:BAM_SORT_STATS_SAMTOOLS:SAMTOOLS_STATS" {
+        tag = { "${meta.reference}:${meta.genomic_region}|${meta.sample_id}_${meta.library_id}" }
+        ext.prefix = { "${meta.sample_id}_${meta.library_id}_${meta.reference}_${meta.genomic_region}_mapAD" }
+        publishDir = [
+            enabled: false
+        ]
+    }
+
     withName: SAMTOOLS_SORT_DEDUPPED {
         tag = { "${meta.reference}|${meta.sample_id}_${meta.library_id}" }
         ext.prefix = { "${meta.sample_id}_${meta.library_id}_${meta.reference}_dedupped" }
diff --git a/nextflow.config b/nextflow.config
@@ -123,7 +123,6 @@ params {
     mapping_bowtie2_trim3       = 0
     mapping_bowtie2_maxins      = 500
     mapping_mapad_p             = 0.03
-    mapping_mapad_lib           = 'single_stranded'
     mapping_mapad_f             = 0.5
     mapping_mapad_t             = 0.5
     mapping_mapad_d             = 0.02
diff --git a/nextflow_schema.json b/nextflow_schema.json
@@ -611,14 +611,6 @@
                     "help_text": "Configures the `mapad map -p` parameter, defining how many mismatches are allowed in a read. Default is set following recommendations from [Heide et al. 2024](https://localhost) who tested when aligning to human reference genomes. \n\n> Modifies mapAD parameter: `-p`",
                     "fa_icon": "fas fa-sort-numeric-down"
                 },
-                "mapping_mapad_lib": {
-                    "type": "string",
-                    "default": "single_stranded",
-                    "description": "Specify mapAD's library-preparation parameter, which adapts its aDNA-specific mismatch scoring correspondingly.",
-                    "help_text": "Configures whether mapAD's aDNA-scoring model is set to double- or single-stranded library preparation mode. In single-stranded mode, we expect elevated C->T error rates at both ends of the reads, whereas in double-stranded mode we expect increased C->T error rates at the 5‘ end, but elevated G->A error rates at the 3’ end. Default is set to `single_stranded`.\n\n> Modifies mapAD parameter: `-l/--library`",
-                    "fa_icon": "fas fa-exchange-alt",
-                    "enum": ["single_stranded", "double_stranded"]
-                },
                 "mapping_mapad_f": {
                     "type": "number",
                     "default": 0.5,
diff --git a/subworkflows/local/map.nf b/subworkflows/local/map.nf
@@ -17,7 +17,7 @@ include { SAMTOOLS_FLAGSTAT as SAMTOOLS_FLAGSTAT_MAPPED } from '../../modules/nf
 workflow MAP {
     take:
     reads // [ [meta], [read1, reads2] ] or [ [meta], [read1] ]
-    index // [ [meta], [ index ] ]
+    index // [ [meta], [ index ], [ fasta ] ]
 
     main:
     ch_versions       = Channel.empty()
@@ -41,7 +41,7 @@ workflow MAP {
             .groupTuple()
 
         ch_input_for_mapping = sharded_reads
-            .combine(index)
+            .combine(index.map{ meta, index, fasta -> [ meta, index ] })
             .multiMap {
                 meta, reads, meta2, index ->
                     new_meta = meta.clone()
@@ -52,7 +52,7 @@ workflow MAP {
 
     } else {
         ch_input_for_mapping = reads
-            .combine(index)
+            .combine(index.map{ meta, index, fasta -> [ meta, index ] })
             .multiMap {
                 meta, reads, meta2, index ->
                     new_meta = meta.clone()
@@ -63,7 +63,7 @@ workflow MAP {
     }
 
     if ( params.mapping_tool == 'bwaaln' ) {
-        ch_index_for_mapping = index
+        ch_index_for_mapping = index.map{ meta, index, fasta -> [ meta, index ] }
         ch_reads_for_mapping = reads
 
         FASTQ_ALIGN_BWAALN ( ch_reads_for_mapping, ch_index_for_mapping )
@@ -80,7 +80,7 @@ workflow MAP {
 
     } else if ( params.mapping_tool == 'bwamem' ) {
         ch_input_for_mapping = reads
-                            .combine( index )
+                            .combine( index.map{ meta, index, fasta -> [ meta, index ] } )
                             .multiMap {
                                 meta, reads, meta2, index ->
                                     new_meta = meta + [ reference: meta2.id ]
@@ -98,7 +98,7 @@ workflow MAP {
 
     } else if ( params.mapping_tool == 'bowtie2' ) {
         ch_input_for_mapping = reads
-                            .combine( index )
+                            .combine( index.map{ meta, index, fasta -> [ meta, index ] } )
                             .multiMap {
                                 meta, reads, meta2, index ->
                                     new_meta = meta + [ reference: meta2.id ]
@@ -115,12 +115,23 @@ workflow MAP {
         ch_mapped_lane_bai = params.fasta_largeref ? SAMTOOLS_INDEX_BT2.out.csi : SAMTOOLS_INDEX_BT2.out.bai
 
     } else if ( params.mapping_tool == 'mapad' ) {
+        ch_input_for_mapping = reads
+                            .combine( index )
+                            .multiMap {
+                                meta, reads, meta2, index, fasta ->
+                                    new_meta = meta + [ reference: meta2.id ]
+                                    reads: [ new_meta, reads ]
+                                    index: [ meta2, index ]
+                                    fasta: [ meta2, fasta ]
+                                    strandedness: meta.strandedness!="single"
+                            }
+
         FASTQ_ALIGN_MAPAD (
-                reads,
-                index,
-                params.fasta,
+                ch_input_for_mapping.reads,
+                ch_input_for_mapping.index,
+                ch_input_for_mapping.fasta,
                 params.mapping_mapad_p,
-                params.mapping_mapad_lib=="double_stranded" ? true : false,
+                ch_input_for_mapping.strandedness,
                 params.mapping_mapad_f,
                 params.mapping_mapad_t,
                 params.mapping_mapad_d,
diff --git a/workflows/eager.nf b/workflows/eager.nf
@@ -182,7 +182,7 @@ workflow EAGER {
     ch_reference_for_mapping = REFERENCE_INDEXING.out.reference
             .map{
                 meta, fasta, fai, dict, index, circular_target ->
-                [ meta, index ]
+                [ meta, index, fasta ]
             }
 
     MAP ( ch_reads_for_mapping, ch_reference_for_mapping )

Original file line number	Diff line number	Diff line change
`@@ -182,7 +182,7 @@ workflow EAGER {`
`182`	`182`	`ch_reference_for_mapping = REFERENCE_INDEXING.out.reference`
`183`	`183`	`.map{`
`184`	`184`	`meta, fasta, fai, dict, index, circular_target ->`
`185`		`- [ meta, index ]`
	`185`	`+ [ meta, index, fasta ]`
`186`	`186`	`}`
`187`	`187`
`188`	`188`	`MAP ( ch_reads_for_mapping, ch_reference_for_mapping )`