Add ApplyBQSR, SamtoolStats, BamQC

main.nf

@@ -277,10 +277,13 @@ include { BUILD_INDICES } from './modules/local/subworkflow/build_indices'
 ================================================================================
 */
 
-include { GATK_BASERECALIBRATOR  as BASERECALIBRATOR }  from './modules/nf-core/software/gatk_baserecalibrator'
-include { GATK_GATHERBQSRREPORTS as GATHERBQSRREPORTS } from './modules/nf-core/software/gatk_gatherbqsrreports'
-include { GATK_MARKDUPLICATES    as MARKDUPLICATES }    from './modules/nf-core/software/gatk_markduplicates'
-include { MULTIQC }                                     from './modules/nf-core/software/multiqc'
+include { GATK_BASERECALIBRATOR  as BASERECALIBRATOR }      from './modules/nf-core/software/gatk_baserecalibrator'
+include { GATK_GATHERBQSRREPORTS as GATHERBQSRREPORTS }     from './modules/nf-core/software/gatk_gatherbqsrreports'
+include { GATK_MARKDUPLICATES    as MARKDUPLICATES }        from './modules/nf-core/software/gatk_markduplicates'
+include { GATK_APPLYBQSR         as APPLYBQSR }             from './modules/nf-core/software/gatk_applybqsr'
+include { SAMTOOLS_INDEX         as SAMTOOLS_INDEX_MAPPED } from './modules/nf-core/software/samtools_index'
+include { SAMTOOLS_INDEX         as SAMTOOLS_INDEX_RECAL }  from './modules/nf-core/software/samtools_index'
+include { MULTIQC }                                         from './modules/nf-core/software/multiqc'
 
 /*
 ================================================================================
@@ -369,20 +372,20 @@ workflow {
 
     BWAMEM2_MEM(QC_TRIM.out.reads, bwa, fasta, fai, params.modules['bwamem2_mem'])
 
-    BWAMEM2_MEM.out.map{ meta, bam, bai ->
+    BWAMEM2_MEM.out.map{ meta, bam -> //, bai ->
         patient = meta.patient
         sample  = meta.sample
         gender  = meta.gender
         status  = meta.status
-        [patient, sample, gender, status, bam, bai]
+        [patient, sample, gender, status, bam] //, bai]
     }.groupTuple(by: [0,1])
         .branch{
             single:   it[4].size() == 1
             multiple: it[4].size() > 1
         }.set{ bam }
 
     bam_single = bam.single.map {
-        patient, sample, gender, status, bam, bai ->
+        patient, sample, gender, status, bam -> //, bai ->
 
         def meta = [:]
         meta.patient = patient
@@ -391,11 +394,11 @@ workflow {
         meta.status = status[0]
         meta.id = sample
 
-        [meta, bam[0], bai[0]]
+        [meta, bam[0]] // , bai[0]]
     }
 
     bam_multiple = bam.multiple.map {
-        patient, sample, gender, status, bam, bai ->
+        patient, sample, gender, status, bam -> //, bai ->
 
         def meta = [:]
         meta.patient = patient
@@ -404,12 +407,12 @@ workflow {
         meta.status = status[0]
         meta.id = sample
 
-        [meta, bam, bai]
+        [meta, bam]
     }
 
     // multipleBam = multipleBam.mix(multipleBamSentieon)
 
-    bam_mapped = bam_single.mix(MERGE_BAM_MAPPED(bam_multiple))
+    bam_mapped = bam_single.mix(SAMTOOLS_INDEX_MAPPED(MERGE_BAM_MAPPED(bam_multiple)))
 
     report_markduplicates = Channel.empty()
     bam_markduplicates    = bam_mapped
@@ -447,8 +450,70 @@ workflow {
         }
 
         GATHERBQSRREPORTS(recaltable)
+        // if ('baserecalibrator' in skip_qc) baseRecalibratorReport.close()
     }
 
+    // (recalTableTSV, recalTableSampleTSV) = recalTableTSV.mix(recalTableTSVnoInt).into(2)
+// // Create TSV files to restart from this step
+// if (params.skip_markduplicates) {
+//     recalTableTSV.map { idPatient, idSample ->
+//         status = status_map[idPatient, idSample]
+//         gender = gender_map[idPatient]
+//         bam = "${params.outdir}/Preprocessing/${idSample}/Mapped/${idSample}.bam"
+//         bai = "${params.outdir}/Preprocessing/${idSample}/Mapped/${idSample}.bam.bai"
+//         recalTable = "${params.outdir}/Preprocessing/${idSample}/Mapped/${idSample}.recal.table"
+//         "${idPatient}\t${gender}\t${status}\t${idSample}\t${bam}\t${bai}\t${recalTable}\n"
+//     }.collectFile(
+//         name: 'mapped_no_duplicates_marked.tsv', sort: true, storeDir: "${params.outdir}/Preprocessing/TSV"
+//     )
+
+//     recalTableSampleTSV
+//         .collectFile(storeDir: "${params.outdir}/Preprocessing/TSV/") {
+//             idPatient, idSample ->
+//             status = status_map[idPatient, idSample]
+//             gender = gender_map[idPatient]
+//             bam = "${params.outdir}/Preprocessing/${idSample}/Mapped/${idSample}.bam"
+//             bai = "${params.outdir}/Preprocessing/${idSample}/Mapped/${idSample}.bam.bai"
+//             recalTable = "${params.outdir}/Preprocessing/${idSample}/Mapped/${idSample}.recal.table"
+//             ["mapped_no_duplicates_marked_${idSample}.tsv", "${idPatient}\t${gender}\t${status}\t${idSample}\t${bam}\t${bai}\t${recalTable}\n"]
+//     }
+// } else {
+//     recalTableTSV.map { idPatient, idSample ->
+//     status = status_map[idPatient, idSample]
+//     gender = gender_map[idPatient]
+//     bam = "${params.outdir}/Preprocessing/${idSample}/DuplicatesMarked/${idSample}.md.bam"
+//     bai = "${params.outdir}/Preprocessing/${idSample}/DuplicatesMarked/${idSample}.md.bam.bai"
+//     recalTable = "${params.outdir}/Preprocessing/${idSample}/DuplicatesMarked/${idSample}.recal.table"
+
+//         "${idPatient}\t${gender}\t${status}\t${idSample}\t${bam}\t${bai}\t${recalTable}\n"
+//     }.collectFile(
+//         name: 'duplicates_marked.tsv', sort: true, storeDir: "${params.outdir}/Preprocessing/TSV"
+//     )
+
+//     recalTableSampleTSV
+//         .collectFile(storeDir: "${params.outdir}/Preprocessing/TSV/") {
+//             idPatient, idSample ->
+//             status = status_map[idPatient, idSample]
+//             gender = gender_map[idPatient]
+//             bam = "${params.outdir}/Preprocessing/${idSample}/DuplicatesMarked/${idSample}.md.bam"
+//             bai = "${params.outdir}/Preprocessing/${idSample}/DuplicatesMarked/${idSample}.md.bam.bai"
+//             recalTable = "${params.outdir}/Preprocessing/${idSample}/DuplicatesMarked/${idSample}.recal.table"
+//             ["duplicates_marked_${idSample}.tsv", "${idPatient}\t${gender}\t${status}\t${idSample}\t${bam}\t${bai}\t${recalTable}\n"]
+//     }
+
+//}
+    bam_applybqsr = MARKDUPLICATES.out.bam.join(GATHERBQSRREPORTS.out.table) //by:[0]
+    bam_applybqsr = bam_applybqsr.combine(BUILD_INDICES.out.intervals)
+
+// if (step == 'recalibrate') bamApplyBQSR = input_sample
+
+    APPLYBQSR(bam_applybqsr, dict, fasta, fai)
+
+    APPLYBQSR.out.dump()
+    // (bam_recalibrated_to_merge, bam_recalibrated_to_index) = bam_recalibrated_to_merge.groupTuple(by:[0, 1]).into(2)
+    MERGE_BAM_RECAL(APPLYBQSR.out)
+    SAMTOOLS_INDEX_RECAL(MERGE_BAM_RECAL.out)
+
     OUTPUT_DOCUMENTATION(
         output_docs,
         output_docs_images)

modules/local/process/bwamem2_mem.nf

@@ -16,7 +16,7 @@ process BWAMEM2_MEM {
         val options
 
     output:
-        tuple val(meta), path("*.bam"), path("*.bai")
+        tuple val(meta), path("*.bam")//, path("*.bai")
 
     script:
     CN = params.sequencing_center ? "CN:${params.sequencing_center}\\t" : ""
@@ -31,7 +31,7 @@ process BWAMEM2_MEM {
         ${fasta} ${reads} | \
     samtools sort --threads ${task.cpus} -m 2G - > ${meta.id}.bam
 
-    samtools index ${meta.id}.bam
+    # samtools index ${meta.id}.bam
 
     echo \$(bwa-mem2 version 2>&1) > bwa-mem2.version.txt
     """

modules/local/process/merge_bam.nf

@@ -4,14 +4,16 @@ process MERGE_BAM {
     tag "${meta.id}"
 
     input:
-        tuple val(meta), path(bam), path(bai)
+        tuple val(meta), path(bam)//, path(bai) optional: true
 
     output:
-        tuple val(meta), path("${meta.sample}.bam"), path("${meta.sample}.bam.bai")
+        tuple val(meta), path("${meta.sample}.bam")//, path("${meta.sample}.bam.bai") optional: true
 
+    //     when: !(params.no_intervals)
+//    samtools merge --threads ${task.cpus} ${idSample}.bam ${bam}
+//      samtools merge --threads ${task.cpus} ${idSample}.recal.bam ${bam}
     script:
     """
     samtools merge --threads ${task.cpus} ${meta.sample}.bam ${bam}
-    samtools index ${meta.sample}.bam
     """
 }

modules/nf-core/software/gatk_applybqsr.nf

@@ -0,0 +1,28 @@
+process GATK_APPLYBQSR {
+    label 'memory_singleCPU_2_task'
+    label 'cpus_2'
+
+    tag "${meta.id}-${interval.baseName}"
+
+    input:
+        tuple val(meta), path(bam), path(bai), path(recalibrationReport), file(interval)
+        path dict
+        path fasta
+        path fai
+
+    output:
+        tuple val(meta), path("${prefix}${meta.sample}.recal.bam") 
+
+    script:
+    prefix = params.no_intervals ? "" : "${interval.baseName}_"
+    options_intervals = params.no_intervals ? "" : "-L ${interval}"
+    """
+    gatk --java-options -Xmx${task.memory.toGiga()}g \
+        ApplyBQSR \
+        -R ${fasta} \
+        --input ${bam} \
+        --output ${prefix}${meta.sample}.recal.bam \
+        ${options_intervals} \
+        --bqsr-recal-file ${recalibrationReport}
+    """
+}

modules/nf-core/software/samtools_index.nf

@@ -0,0 +1,19 @@
+process SAMTOOLS_INDEX {
+   label 'cpus_8'
+
+    tag "${meta.id}"
+
+    input:
+        tuple val(meta), path(bam)
+
+    output:
+        tuple val(meta), path(bam), path("*.bai")
+//    samtools index ${idSample}.bam
+
+    //    samtools index ${idSample}.recal.bam
+
+    script:
+    """
+    samtools index $bam
+    """
+}