Merge pull request #32 from Juke34/bwamem2

add bwamem2

Author: Juke34

README.md

@@ -50,6 +50,7 @@ You can choose to run one or several aligner in parallel.
 | bowtie2 | ✅ | ✅ | ⚠️ | ⚠️ |
 | bwaaln | ✅ | ✅ R1 and R2 independently aligned then merged with bwa sampe | ✅ | ✅ |
 | bwamem | ✅ | ✅ | ⚠️ | ⚠️ |
+| bwamem2 | ✅ | ✅ | ⚠️ | ⚠️ |
 | bwasw | ✅ | ✅ | ⚠️ | ⚠️ |
 | graphmap2 | ⚠️ | ⚠️ R1 and R2 independently aligned then merged with cat | ✅ | ✅ |
 | hisat2 | ✅ | ✅ | ⚠️ | ⚠️ |
@@ -83,6 +84,7 @@ It is then translated to the correct option in the following aligners:
 | bowtie2 | --fr / --rf / --ff |  ISF ISR IU / OSF OSR OU / MSF MSR MU| read orientation |
 | bwaaln | 🚫 | 🚫 | 🚫 |
 | bwamem | 🚫 | 🚫 | 🚫 |
+| bwamem2 | 🚫 | 🚫 | 🚫 |
 | bwasw | 🚫 | 🚫 | 🚫 |
 | graphmap2 | 🚫 | 🚫 | 🚫 |
 | hisat2 | --rna-strandness [ F / R / FR / RF ] | SF / SR / ISF OSF MSF / ISR OSR MSR | strand information |
@@ -117,6 +119,7 @@ If you provide an annotation file the pipeline will pass automatically the file
 | bowtie2 | 🚫 |
 | bwaaln | 🚫 |
 | bwamem | 🚫 |
+| bwamem2 | 🚫 |
 | bwasw | 🚫 |
 | graphmap2 | GTF (--gtf)  |
 | hisat2 | 🚫 |
@@ -316,7 +319,7 @@ On success you should get a message looking like this:
         --reads                     path to the reads file or folder
         --reads_extension           extension of the reads files (default: .fastq.gz)
         --genome                    path to the genome file
-        --aligner                   aligner(s) to use among this list (comma or space separated) [bbmap, bowtie, bowtie2, bwaaln, bwamem, bwasw, graphmap2, hisat2, kallisto, minimap2, novoalign, nucmer, ngmlr, star, subread, sublong]
+        --aligner                   aligner(s) to use among this list (comma or space separated) [bbmap, bowtie, bowtie2, bwaaln, bwamem, bwamem2, bwasw, graphmap2, hisat2, kallisto, minimap2, novoalign, nucmer, ngmlr, star, subread, sublong]
         --outdir                    path to the output directory (default: alignment_results)
         --annotation                [Optional][used by graphmap2, STAR, subread] Absolute path to the annotation file (gtf or gff3)
 
@@ -340,6 +343,7 @@ On success you should get a message looking like this:
         --bowtie2_options           additional options for bowtie2
         --bwaaln_options            additional options for bwaaln
         --bwamem_options            additional options for bwamem
+        --bwamem2_options            additional options for bwamem2
         --bwasw_options             additional options for bwasw
         --graphmap2_options         additional options for graphmap2
         --hisat2_options            additional options for hisat2

aline.nf

@@ -38,13 +38,14 @@ params.annotation = ""
 params.trimming_fastp = false
 
 // Aligner params
-align_tools = [ 'bbmap', 'bowtie', 'bowtie2', 'bwaaln', 'bwamem', 'bwasw', 'graphmap2', 'hisat2', 'kallisto', 'minimap2', 'novoalign', 'nucmer', 'ngmlr', 'star', 'subread', 'sublong' ]
+align_tools = [ 'bbmap', 'bowtie', 'bowtie2', 'bwaaln', 'bwamem', 'bwamem2', 'bwasw', 'graphmap2', 'hisat2', 'kallisto', 'minimap2', 'novoalign', 'nucmer', 'ngmlr', 'star', 'subread', 'sublong' ]
 params.aligner = ''
 params.bbmap_options = ''
 params.bowtie_options = ''
 params.bowtie2_options = ''
 params.bwaaln_options = ''
 params.bwamem_options = ''
+params.bwamem2_options = ''
 params.bwasw_options = ''
 params.graphmap2_options = '' // owler option is possible
 params.hisat2_options = ''
@@ -158,13 +159,13 @@ if ("bbmap" in aligner_list && !params.relax){
 }
 
 // --- bwa aln tool ---
-//if ("bwaaln" in aligner_list && !params.relax){
-//    if (params.read_type == "pacbio" || params.read_type == "ont"){
-//        log.warn("""Error: bwaaln is not suitable for long reads.
-//However, if you know what you are doing you can activate the AliNe --relax parameter to use it anyway.\n""")
-//        stop_pipeline = true
-//    }
-//}
+if ("bwaaln" in aligner_list && !params.relax){
+    if (params.read_type == "pacbio" || params.read_type == "ont"){
+        log.warn("""Error: bwaaln is not suitable for long reads.
+However, if you know what you are doing you can activate the AliNe --relax parameter to use it anyway.\n""")
+        stop_pipeline = true
+    }
+}
 
 // --- bwa mem tool ---
 if ("bwamem" in aligner_list && !params.relax){
@@ -180,20 +181,29 @@ if ("bwamem" in aligner_list && !params.relax){
     }
 }
 
-// --- bwa sw tool ---
-if ("bwasw" in aligner_list && !params.relax){
+// --- bwa mem2 tool ---
+if ("bwamem2" in aligner_list && !params.relax){
     if (params.read_type == "pacbio"){
-        if ( !params.bwasw_options.contains(" pacbio") ){
-            params.replace("bwamem_options", "${params.bwasw_options} -x pacbio")
+        if ( !params.bwamem2_options.contains(" pacbio") ){
+            params.replace("bwamem2_options", "${params.bwamem2_options} -x pacbio")
         }
     }
     if (params.read_type == "ont"){
-        if ( !params.bwasw_options.contains(" ont2d") ){
-            params.replace("bwamem_options", "${params.bwasw_options} -x ont2d")
+        if ( !params.bwamem2_options.contains(" ont2d") ){
+            params.replace("bwamem2_options", "${params.bwamem2_options} -x ont2d")
         }
     }
 }
 
+// --- bwa sw tool ---
+if ("bwasw" in aligner_list && !params.relax){
+    if (params.read_type == "pacbio" || params.read_type == "ont"){
+        log.warn("""Error: bwasw is not suitable for long reads.
+However, if you know what you are doing you can activate the AliNe --relax parameter to use it anyway.\n""")
+        stop_pipeline = true
+    }
+}
+
 // --- graphmap2 tool ---
 if ("graphmap2" in aligner_list ){
     if (annotation_file && !params.graphmap2_options.contains("--gtf ") ){
@@ -353,11 +363,12 @@ include {bbmap_index; bbmap} from "$baseDir/modules/bbmap.nf"
 include {bowtie_index; bowtie} from "$baseDir/modules/bowtie.nf"
 include {bowtie2_index; bowtie2} from "$baseDir/modules/bowtie2.nf"
 include {bwa_index; bwaaln; bwamem; bwasw} from "$baseDir/modules/bwa.nf"
+include {bwamem2_index; bwamem2} from "$baseDir/modules/bwamem2.nf"
 include {seqkit_convert} from "$baseDir/modules/seqkit.nf"
 include {graphmap2_index; graphmap2} from "$baseDir/modules/graphmap2.nf"
 include {fastp} from "$baseDir/modules/fastp.nf"
 include {fastqc as fastqc_raw; fastqc as fastqc_fastp; fastqc as fastqc_ali_bbmap; fastqc as fastqc_ali_bowtie ; fastqc as fastqc_ali_bowtie2 ; 
-         fastqc as fastqc_ali_bwaaln; fastqc as fastqc_ali_bwamem; fastqc as fastqc_ali_bwasw; fastqc as fastqc_ali_graphmap2 ; 
+         fastqc as fastqc_ali_bwaaln; fastqc as fastqc_ali_bwamem; fastqc as fastqc_ali_bwamem2; fastqc as fastqc_ali_bwasw; fastqc as fastqc_ali_graphmap2 ; 
          fastqc as fastqc_ali_hisat2; fastqc as fastqc_ali_kallisto; fastqc as fastqc_ali_minimap2; fastqc as fastqc_ali_ngmlr; 
          fastqc as fastqc_ali_novoalign ; fastqc as fastqc_ali_nucmer; fastqc as fastqc_ali_star; fastqc as fastqc_ali_subread ; 
          fastqc as fastqc_ali_sublong } from "$baseDir/modules/fastqc.nf"
@@ -370,16 +381,16 @@ include {nucmer} from "$baseDir/modules/mummer4.nf"
 include {novoalign_index; novoalign} from "$baseDir/modules/novoalign.nf"
 include {salmon_index; salmon_guess_lib; set_tuple_withUserLib} from "$baseDir/modules/salmon.nf" 
 include {samtools_sam2bam_nucmer; samtools_sam2bam as samtools_sam2bam_bowtie; samtools_sam2bam as samtools_sam2bam_bowtie2; samtools_sam2bam as samtools_sam2bam_bwaaln; 
-         samtools_sam2bam as samtools_sam2bam_bwamem; samtools_sam2bam as samtools_sam2bam_bwasw; samtools_sam2bam as samtools_sam2bam_graphmap2; 
+         samtools_sam2bam as samtools_sam2bam_bwamem; samtools_sam2bam as samtools_sam2bam_bwamem2; samtools_sam2bam as samtools_sam2bam_bwasw; samtools_sam2bam as samtools_sam2bam_graphmap2; 
          samtools_sam2bam as samtools_sam2bam_hisat2; samtools_sam2bam as samtools_sam2bam_minimap2; 
          samtools_sam2bam as samtools_sam2bam_ngmlr; samtools_sam2bam as samtools_sam2bam_novoalign } from "$baseDir/modules/samtools.nf"
 include {samtools_sort as samtools_sort_bbmap; samtools_sort as samtools_sort_bowtie; samtools_sort as samtools_sort_bowtie2; samtools_sort as samtools_sort_bwaaln; 
-         samtools_sort as samtools_sort_bwamem; samtools_sort as samtools_sort_bwasw; samtools_sort as samtools_sort_graphmap2; 
+         samtools_sort as samtools_sort_bwamem; samtools_sort as samtools_sort_bwamem2; samtools_sort as samtools_sort_bwasw; samtools_sort as samtools_sort_graphmap2; 
          samtools_sort as samtools_sort_hisat2; samtools_sort as samtools_sort_minimap2; samtools_sort as samtools_sort_ngmlr; 
          samtools_sort as samtools_sort_novoalign;  samtools_sort as samtools_sort_nucmer;
          samtools_sort as samtools_sort_sublong } from "$baseDir/modules/samtools.nf"
 include {samtools_stats as samtools_stats_ali_bbmap; samtools_stats as samtools_stats_ali_bowtie ; samtools_stats as samtools_stats_ali_bowtie2 ; 
-         samtools_stats as samtools_stats_ali_bwaaln; samtools_stats as samtools_stats_ali_bwamem; samtools_stats as samtools_stats_ali_bwasw; samtools_stats as samtools_stats_ali_graphmap2 ; 
+         samtools_stats as samtools_stats_ali_bwaaln; samtools_stats as samtools_stats_ali_bwamem; samtools_stats as samtools_stats_ali_bwamem2; samtools_stats as samtools_stats_ali_bwasw; samtools_stats as samtools_stats_ali_graphmap2 ; 
          samtools_stats as samtools_stats_ali_hisat2; samtools_stats as samtools_stats_ali_kallisto; samtools_stats as samtools_stats_ali_minimap2; samtools_stats as samtools_stats_ali_ngmlr; 
          samtools_stats as samtools_stats_ali_novoalign ; samtools_stats as samtools_stats_ali_nucmer; samtools_stats as samtools_stats_ali_star; samtools_stats as samtools_stats_ali_subread ; 
          samtools_stats as samtools_stats_ali_sublong } from "$baseDir/modules/samtools.nf"
@@ -670,7 +681,7 @@ workflow align {
             }
         }
 
-        // ------------------- BWA -----------------
+        // ------------------- BWA ALN/MEM/SW -----------------
         if ("bwaaln" in aligner_list || "bwamem" in aligner_list || "bwasw" in aligner_list){
             // index
             bwa_index(genome.collect(), "alignment/bwa/indicies")
@@ -739,6 +750,31 @@ workflow align {
             }
         }
 
+        // ------------------- BWA MEM2 -----------------
+        if ("bwamem2" in aligner_list){
+            // index
+            bwamem2_index(genome.collect(), "alignment/bwamem2/indicies")
+            // align
+            bwamem2(reads, genome.collect(), bwamem2_index.out.collect(), "alignment/bwamem2/") 
+            logs.concat(bwamem2.out.bwamem2_summary).set{logs} // save log
+            // convert sam to bam
+            samtools_sam2bam_bwamem2(bwamem2.out.tuple_sample_sam)
+            // sort
+            samtools_sort_bwamem2(samtools_sam2bam_bwamem2.out.tuple_sample_bam, "alignment/bwamem2/")
+            samtools_sort_bwamem2.out.tuple_sample_sortedbam.set{bwamem2_ali} // set name
+            // save aligned reads
+            sorted_bam.concat(bwamem2_ali).set{sorted_bam} 
+            // stat on aligned reads
+            if(params.fastqc){
+                fastqc_ali_bwamem2(bwamem2_ali, "fastqc/bwamem2", "bwamem2")
+                logs.concat(fastqc_ali_bwamem2.out).set{logs} // save log
+            }
+            if(params.samtools_stats){
+                samtools_stats_ali_bwamem2(bwamem2_ali, genome.collect(), "samtools_stats/bwamem2", "bwamem2")
+                logs.concat(samtools_stats_ali_bwamem2.out).set{logs} // save log
+            }
+        }
+
         // ------------------- GRAPHMAP2 -----------------
         if ("graphmap2" in aligner_list ){
             // index
@@ -1058,6 +1094,7 @@ def helpMSG() {
         --bowtie2_options           additional options for bowtie2
         --bwaaln_options            additional options for bwaaln
         --bwamem_options            additional options for bwamem
+        --bwamem2_options            additional options for bwamem2
         --bwasw_options             additional options for bwasw
         --graphmap2_options         additional options for graphmap2
         --hisat2_options            additional options for hisat2
@@ -1106,6 +1143,11 @@ def printAlignerOptions(aligner_list, annotation_file, star_index_options) {
     bwamem parameters
         bwamem_options             : ${params.bwamem_options}
     """} 
+    if ("bwamem2" in aligner_list){
+        sentence += """
+    bwamem2 parameters
+        bwamem2_options             : ${params.bwamem2_options}
+    """} 
     if ("bwasw" in aligner_list){
         sentence += """
     bwasw parameters

config/softwares.config

@@ -14,6 +14,9 @@ process {
     withLabel: 'bwa' {
         container = 'quay.io/biocontainers/bwa:0.7.17--he4a0461_11'
     }
+    withLabel: 'bwamem2' {
+        container = 'quay.io/biocontainers/bwa-mem2:2.2.1--he70b90d_8'
+    }
     withLabel: 'fastp' {
         container = 'quay.io/biocontainers/fastp:0.23.4--h125f33a_5'
     }

modules/bwamem2.nf

@@ -0,0 +1,63 @@
+/* Module related to bwamem2
+https://github.com/bwa-mem2/bwa-mem2?tab=readme-ov-file
+
+info:
+We are happy to announce that the index size on disk is down by 8 times and in memory by 4 times due to moving to only one type of FM-index (2bit.64 instead of 2bit.64 and 8bit.32) and 8x compression of suffix array. 
+For example, for human genome, index size on disk is down to ~10GB from ~80GB and memory footprint is down to ~10GB from ~40GB. 
+There is a substantial reduction in index IO time due to the reduction and hardly any performance impact on read mapping. 
+Due to this change in index structure (in commit #4b59796, 10th October 2020), you will need to rebuild the index.
+Added MC flag in the output sam file in commit a591e22. Output should match original bwa-mem version 0.7.17.
+*/ 
+
+/*
+* To index with BWA MEM2
+*/
+process bwamem2_index {
+    label 'bwamem2'
+    tag "$genome_fasta"
+    publishDir "${params.outdir}/${outpath}", mode: 'copy'
+
+    input:
+        path(genome_fasta)
+        val outpath
+
+    output:
+        path("*")
+
+    script:
+        """
+        bwa-mem2 index $genome_fasta -p ${genome_fasta.baseName}
+        """
+}
+
+/*
+* To align with BWA MEM2
+*/
+process bwamem2 {
+    label 'bwamem2'
+    tag "$sample"
+    publishDir "${params.outdir}/${outpath}", pattern: "*bwamem2.log", mode: 'copy'
+
+    input:
+        tuple val(sample), path(reads), val(readtype), val(read_length)
+        path genome
+        path bwa_index_files
+        val outpath
+
+    output:
+        tuple val(sample), path ("*bwamem2.sam"), emit: tuple_sample_sam
+        path "*bwamem2.log",  emit: bwamem2_summary
+
+    script:
+        fileName = reads[0].baseName.replace('.fastq','')
+      
+        if (params.read_type == "short_paired"){
+            """
+            bwa-mem2 mem ${params.bwamem2_options} -t ${task.cpus} ${genome.baseName} ${reads[0]} ${reads[1]} > ${fileName}_bwamem2.sam 2> ${fileName}_bwamem2.log 
+            """
+        } else {
+            """
+            bwa-mem2 mem ${params.bwamem2_options} -t ${task.cpus} ${genome.baseName} ${reads} > ${fileName}_bwamem2.sam 2> ${fileName}_bwamem2.log 
+            """
+        }
+}

profiles/test_illumina_paired.config

@@ -7,7 +7,7 @@
 params {
     reads = "$baseDir/test/illumina/"
     genome = "$baseDir/test/yeast.fa"
-    aligner = 'bbmap,bowtie,bowtie2,bwaaln,bwamem,bwasw,graphmap2,hisat2,minimap2,nucmer,star,subread'
+    aligner = 'bbmap,bowtie,bowtie2,bwaaln,bwamem,bwamem2,bwasw,graphmap2,hisat2,minimap2,nucmer,star,subread'
     star_options = "--genomeSAindexNbases 9" // the default 14 is too large for the genome size=1351857
     multiqc_config = "$baseDir/config/multiqc_conf.yml"
 }

profiles/test_illumina_single.config

@@ -8,7 +8,7 @@ params {
     reads = "$baseDir/test/illumina/"
     genome = "$baseDir/test/yeast.fa"
     params.read_type = "short_single"
-    aligner = 'bbmap,bowtie,bowtie2,bwaaln,bwamem,bwasw,graphmap2,hisat2,kallisto,minimap2,ngmlr,nucmer,star,subread,sublong'
+    aligner = 'bbmap,bowtie,bowtie2,bwaaln,bwamem,bwamem2,bwasw,graphmap2,hisat2,kallisto,minimap2,ngmlr,nucmer,star,subread,sublong'
     trimming_fastp = true
     fastqc = true
     samtools_stats = true

profiles/test_ont.config

@@ -8,7 +8,7 @@ params {
     reads = "$baseDir/test/nanopore"
     genome = "$baseDir/test/yeast.fa"
     read_type = "ont"
-    aligner = 'bbmap,bowtie2,bwaaln,bwamem,bwasw,graphmap2,hisat2,kallisto,minimap2,ngmlr,nucmer,star,subread,sublong'
+    aligner = 'bowtie2,bwamem,graphmap2,hisat2,kallisto,minimap2,ngmlr,nucmer,star,subread,sublong'
     library_type = 'U'
     star_options = '--outFilterMismatchNmax 100 --seedSearchLmax 30   --seedSearchStartLmax 30 --seedPerReadNmax 100000   --seedPerWindowNmax 100 --alignTranscriptsPerReadNmax 100000 --alignTranscriptsPerWindowNmax 10000'
     multiqc_config = "$baseDir/config/multiqc_conf.yml"

profiles/test_pacbio.config

@@ -8,7 +8,7 @@ params {
     reads = "$baseDir/test/pacbio/"
     genome = "$baseDir/test/yeast.fa"
     read_type = "pacbio"
-    aligner = 'bbmap,bowtie,bowtie2,bwaaln,bwamem,bwasw,graphmap2,hisat2,kallisto,minimap2,ngmlr,nucmer,star,subread,sublong'
+    aligner = 'bbmap,bowtie,bowtie2,,bwamem,bwamem2,graphmap2,hisat2,kallisto,minimap2,ngmlr,nucmer,star,subread,sublong'
     star_options = '--outFilterMismatchNmax 100 --seedSearchLmax 30   --seedSearchStartLmax 30 --seedPerReadNmax 100000 --seedPerWindowNmax 100 --alignTranscriptsPerReadNmax 100000 --alignTranscriptsPerWindowNmax 10000'
     star_index_options = '--genomeSAindexNbases 9'
     multiqc_config = "$baseDir/config/multiqc_conf.yml"