Restores mapping to nf-core; Keeps gatk *Spark

conf/base.config

@@ -72,4 +72,4 @@ process {
     container = {(params.annotation_cache && params.vep_cache) ? 'nfcore/sarek:2.5.2' : "nfcore/sarekvep:2.5.2.${params.genome}"}
     errorStrategy = {task.exitStatus == 143 ? 'retry' : 'ignore'}
   }
-}
+}

conf/igenomes.config

@@ -201,5 +201,21 @@ params {
       bwaIndex               = "${params.igenomes_base}/Zea_mays/Ensembl/AGPv3/Sequence/BWAIndex/genome.fa.{amb,ann,bwt,pac,sa}"
       fasta                  = "${params.igenomes_base}/Zea_mays/Ensembl/AGPv3/Sequence/WholeGenomeFasta/genome.fa"
     }
+    'GRCh38_PGP_UK' {
+      acLoci                 = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/ASCAT/1000G_phase3_GRCh38_maf0.3.loci"
+      acLociGC               = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/ASCAT/1000G_phase3_GRCh38_maf0.3.loci.gc"
+      chrDir                 = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Sequence/Chromosomes"
+      chrLength              = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Sequence/Length/Homo_sapiens_assembly38.len"
+      dbsnp                  = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GATKBundle/dbsnp_146.hg38.vcf.gz"
+      dbsnpIndex             = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GATKBundle/dbsnp_146.hg38.vcf.gz.tbi"
+      fasta                  = "s3://lifebit-featured-datasets/pipelines/nf-core-sarek/resources/GRCh38_full_analysis_set_plus_decoy_hla_noChr.fa"
+      germlineResource       = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GermlineResource/gnomAD.r2.1.1.GRCh38.PASS.AC.AF.only.vcf.gz"
+      germlineResourceIndex  = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GermlineResource/gnomAD.r2.1.1.GRCh38.PASS.AC.AF.only.vcf.gz.tbi"
+      intervals              = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/intervals/wgs_calling_regions.hg38.bed"
+      knownIndels            = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GATKBundle/{Mills_and_1000G_gold_standard.indels.hg38,beta/Homo_sapiens_assembly38.known_indels}.vcf.gz"
+      knownIndelsIndex       = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh38/Annotation/GATKBundle/{Mills_and_1000G_gold_standard.indels.hg38,beta/Homo_sapiens_assembly38.known_indels}.vcf.gz.tbi"
+      snpeffDb               = "GRCh38.86"
+      vepCacheVersion        = "95"
+    }
   }
 }

main.nf

@@ -793,6 +793,12 @@ process MapReads {
     extra = status == 1 ? "-B 3" : ""
     convertToFastq = hasExtension(inputFile1, "bam") ? "gatk --java-options -Xmx${task.memory.toGiga()}g SamToFastq --INPUT=${inputFile1} --FASTQ=/dev/stdout --INTERLEAVE=true --NON_PF=true | \\" : ""
     input = hasExtension(inputFile1, "bam") ? "-p /dev/stdin - 2> >(tee ${inputFile1}.bwa.stderr.log >&2)" : "${inputFile1} ${inputFile2}"
+    // Pseudo-code: Add soft-coded memory allocation to the two tools:
+    bwa_memory  = task.memory.toGiga() * 0.60
+    sort_memory = task.memory.toGiga() * 0.40
+    // Pseudo-code: Add soft-coded memory allocation to the two tools:
+    bwa_cpus  = (task.cpus * 0.60).toInteger()
+    sort_cpus = (task.cpus * 0.40).toInteger()
     """
         ${convertToFastq}
         bwa mem -K 100000000 -R \"${readGroup}\" ${extra} -t ${task.cpus} -M ${fasta} \
@@ -940,10 +946,12 @@ process IndexBamFile {
 
 // STEP 2: MARKING DUPLICATES
 
-process MarkDuplicates {
-    label 'cpus_16'
+process MarkDuplicatesSpark {
+    label 'cpus_max'
+    label 'memory_max'
 
     tag {idPatient + "-" + idSample}
+    echo true
 
     publishDir params.outdir, mode: params.publishDirMode,
         saveAs: {
@@ -956,23 +964,21 @@ process MarkDuplicates {
         set idPatient, idSample, file("${idSample}.bam") from mergedBam
 
     output:
-        set idPatient, idSample, file("${idSample}.md.bam"), file("${idSample}.md.bai") into duplicateMarkedBams
+        set idPatient, idSample, file("${idSample}.md.bam"), file("${idSample}.md.bam.bai") into duplicateMarkedBams
         file ("${idSample}.bam.metrics") into markDuplicatesReport
 
     when: params.knownIndels
 
     script:
     markdup_java_options = task.memory.toGiga() > 8 ? params.markdup_java_options : "\"-Xms" +  (task.memory.toGiga() / 2).trunc() + "g -Xmx" + (task.memory.toGiga() - 1) + "g\""
     """
-    gatk --java-options ${markdup_java_options} \
-        MarkDuplicates \
-        --MAX_RECORDS_IN_RAM 50000 \
-        --INPUT ${idSample}.bam \
-        --METRICS_FILE ${idSample}.bam.metrics \
-        --TMP_DIR . \
-        --ASSUME_SORT_ORDER coordinate \
-        --CREATE_INDEX true \
-        --OUTPUT ${idSample}.md.bam
+    gatk MarkDuplicatesSpark \
+        --input ${idSample}.bam \
+        --output ${idSample}.md.bam \
+        --tmp-dir . \
+        --verbosity DEBUG \
+        --create-output-bam-index true \
+        --spark-runner LOCAL --spark-master local[${task.cpus}]
     """
 }
 
@@ -1044,10 +1050,11 @@ process SentieonDedup {
 
 // STEP 3: CREATING RECALIBRATION TABLES
 
-process BaseRecalibrator {
+process BaseRecalibratorSpark {
     label 'cpus_1'
 
-    tag {idPatient + "-" + idSample + "-" + intervalBed.baseName}
+    tag {idPatient + "-" + idSample + "-" + intervalBed.simpleName}
+    echo true
 
     input:
         set idPatient, idSample, file(bam), file(bai), file(intervalBed) from bamBaseRecalibrator
@@ -1073,15 +1080,16 @@ process BaseRecalibrator {
     // TODO: --use-original-qualities ???
     """
     gatk --java-options -Xmx${task.memory.toGiga()}g \
-        BaseRecalibrator \
-        -I ${bam} \
-        -O ${prefix}${idSample}.recal.table \
-        --tmp-dir /tmp \
-        -R ${fasta} \
+    BaseRecalibratorSpark \
+        --input ${bam} \
+        --output ${prefix}${idSample}.recal.table \
+        --tmp-dir . \
+        --reference ${fasta} \
         ${intervalsOptions} \
         ${dbsnpOptions} \
         ${knownOptions} \
-        --verbosity INFO
+        --verbosity INFO \
+        --spark-runner LOCAL --spark-master local[${task.cpus}]
     """
 }
 
@@ -1164,7 +1172,7 @@ bamApplyBQSR = bamApplyBQSR.dump(tag:'BAM + BAI + RECAL TABLE + INT')
 
 // STEP 4: RECALIBRATING
 
-process ApplyBQSR {
+process ApplyBQSRSpark {
     label 'memory_singleCPU_2_task'
     label 'cpus_2'
 
@@ -1184,12 +1192,14 @@ process ApplyBQSR {
     intervalsOptions = params.no_intervals ? "" : "-L ${intervalBed}"
     """
     gatk --java-options -Xmx${task.memory.toGiga()}g \
-        ApplyBQSR \
-        -R ${fasta} \
+        ApplyBQSRSpark \
+        --reference ${fasta} \
         --input ${bam} \
         --output ${prefix}${idSample}.recal.bam \
         ${intervalsOptions} \
-        --bqsr-recal-file ${recalibrationReport}
+        --bqsr-recal-file ${recalibrationReport} \
+        --verbosity DEBUG \
+        --spark-runner LOCAL --spark-master local[${task.cpus}] &> applyBQSRspark.log.txt
     """
 }
 
@@ -1284,7 +1294,7 @@ bamRecalSentieonSampleTSV
         ["recalibrated_sentieon_${idSample}.tsv", "${idPatient}\t${gender}\t${status}\t${idSample}\t${bam}\t${bai}\n"]
 }
 
-// STEP 4.5: MERGING THE RECALIBRATED BAM FILES
+// STEP 4.5.1: MERGING THE RECALIBRATED BAM FILES
 
 process MergeBamRecal {
     label 'cpus_8'
@@ -1300,6 +1310,7 @@ process MergeBamRecal {
         set idPatient, idSample, file("${idSample}.recal.bam"), file("${idSample}.recal.bam.bai") into bamRecal
         set idPatient, idSample, file("${idSample}.recal.bam") into bamRecalQC
         set idPatient, idSample into bamRecalTSV
+        file("${idSample}.recal.bam") into (bamGenomeChronicler, bamGenomeChroniclerToPrint)
 
     when: !(params.no_intervals)
 
@@ -1309,6 +1320,40 @@ process MergeBamRecal {
     samtools index ${idSample}.recal.bam
     """
 }
+bamGenomeChroniclerToPrint.view()
+
+// TODO: Bind this with HaplotypeCaller output and migrate process chuck after HaplotypeCasller + VEP
+Channel.fromPath(params.vepFile)
+       .ifEmpty { exit 1, "--vepFile not specified or no file found at that destination with the suffix .html. Please make sure to provide the file path correctly}" }
+       .set { vepGenomeChronicler }
+
+
+// STEP 4.5.2: RUNNING GenomeChronicler FOR THE RECALIBRATED BAM FILES
+// TODO: Update this when there is a different VEP html report for each bam
+process RunGenomeChronicler {
+  tag "$bam"
+  publishDir "$params.outdir/GenomeChronicler", mode: 'copy'
+  echo true
+
+  input:
+  file(bam) from bamGenomeChronicler
+  each file(vep) from vepGenomeChronicler
+
+  output:
+  file("results_${bam.simpleName}") into chronicler_results
+
+  script:
+
+  optional_argument = vep.endsWith("no_vepFile.txt") ? '' : "--vepFile ${vep}"
+
+  """
+  genomechronicler \
+  --resultsDir '/GenomeChronicler' \
+  --bamFile $bam $optional_argument &> STDERR.txt
+  cp -r /GenomeChronicler/results/results_${bam.simpleName} .
+  mv STDERR.txt results_${bam.simpleName}/
+  """
+}
 
 // STEP 4.5': INDEXING THE RECALIBRATED BAM FILES
 
@@ -2830,7 +2875,7 @@ compressVCFsnpEffOut = compressVCFsnpEffOut.dump(tag:'VCF')
 process VEP {
     label 'VEP'
     label 'cpus_4'
-
+    echo true
     tag {"${idSample} - ${variantCaller} - ${vcf}"}
 
     publishDir params.outdir, mode: params.publishDirMode, saveAs: {
@@ -2861,6 +2906,8 @@ process VEP {
     cadd = (params.cadd_cache && params.cadd_WG_SNVs && params.cadd_InDels) ? "--plugin CADD,whole_genome_SNVs.tsv.gz,InDels.tsv.gz" : ""
     genesplicer = params.genesplicer ? "--plugin GeneSplicer,/opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/bin/genesplicer,/opt/conda/envs/nf-core-sarek-${workflow.manifest.version}/share/genesplicer-1.0-1/human,context=200,tmpdir=\$PWD/${reducedVCF}" : "--offline"
     """
+    hello_message="I am here at VEP"
+    echo $hello_message
     mkdir ${reducedVCF}
 
     vep \
@@ -2983,7 +3030,7 @@ compressVCFOutVEP = compressVCFOutVEP.dump(tag:'VCF')
 // STEP MULTIQC
 
 process MultiQC {
-    publishDir "${params.outdir}/Reports/MultiQC", mode: params.publishDirMode
+    publishDir "${params.outdir}/MultiQC", mode: params.publishDirMode
 
     input:
         file (multiqcConfig) from Channel.value(params.multiqc_config ? file(params.multiqc_config) : "")

nextflow.config

@@ -10,14 +10,15 @@ params {
 
   // Workflow flags
   annotateTools = null // Only with --step annotate
-  genome = 'GRCh38'
+  genome = 'GRCh38_PGP_UK'
   input = null // No default input
   noGVCF = null // g.vcf are produced by HaplotypeCaller
   noStrelkaBP = null // Strelka will use Manta candidateSmallIndels if available
   no_intervals = null // Intervals will be built from the fasta file
   skipQC = null // All QC tools are used
   step = 'mapping' // Starts with mapping
   tools = null // No default Variant Calling or Annotation tools
+  vepFile = 'https://raw.githubusercontent.com/cgpu/sarek/master/assets/no_vepFile.txt'
 
   // Workflow settings
   annotation_cache = null // Annotation cache disabled
@@ -89,6 +90,75 @@ params {
 // Developmental code should specify dev
 process.container = 'nfcore/sarek:2.5.2'
 
+
+process {
+
+    withName: MapReads {
+    container = "nfcore/sarek:2.5.1"
+    cpus = 30
+    memory = 210.GB
+    maxForks = 2
+  }
+
+  withName: BaseRecalibratorSpark {
+    container = "broadinstitute/gatk:4.1.4.0"
+    cpus = 8
+    memory = 32.GB
+    maxForks = 96
+  }
+
+    withName: MarkDuplicatesSpark {
+    container = "broadinstitute/gatk:4.1.3.0"
+    cpus = 28
+    memory = 180.GB
+    maxForks = 2
+  }
+
+  withName: RunGenomeChronicler {
+    container = "lifebitai/genomechronicler:pgp-uk-5513c6f"
+    cpus = 1
+    memory = 4.GB
+    maxForks = 96
+  }
+
+  withName: ApplyBQSRSpark {
+    container = "broadinstitute/gatk:4.1.4.0" 
+    cpus = 4
+    memory = 16.GB
+    maxForks = 64
+  }
+
+    withName: HaplotypeCaller {
+    container = "broadinstitute/gatk:4.1.4.0"
+    cpus = 2
+    memory = 8.GB
+    maxForks = 64
+  }
+
+  withName: Mutect2 {
+    container = "broadinstitute/gatk:4.1.4.0"
+    cpus = 2
+    memory = 8.GB
+    maxForks = 32
+  }
+
+  withName: PileupSummariesForMutect2 {
+    container = "broadinstitute/gatk:4.1.4.0"
+    cpus = 1
+    memory = 4.GB
+    maxForks = 96
+  }
+
+  withName: MultiQC {
+    container = "nfcore/sarek:2.5.1"
+    cpus = 8
+    memory = 32.GB
+    maxForks = 1
+    errorStrategy = 'retry'
+    maxRetries = 4
+  }
+}
+
 // Load base.config by default for all pipelines
 includeConfig 'conf/base.config'