Fixing mutect2 PR

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/Bgzip.scala

@@ -28,7 +28,7 @@ class Bgzip(val parent: Configurable) extends BiopetCommandLineFunction {
   var input: List[File] = Nil
 
   @Output(doc = "Compressed output file", required = false)
-  var output: File = null
+  var output: File = _
 
   var f: Boolean = config("f", default = false)
   executable = config("exe", default = "bgzip", freeVar = false)
@@ -39,7 +39,7 @@ class Bgzip(val parent: Configurable) extends BiopetCommandLineFunction {
     if (output == null && !outputAsStdout) Logging.addError("Output is missing for Bgzip")
   }
 
-  def cmdLine =
+  def cmdLine: String =
     required(executable) +
       conditional(f, "-f") +
       " -c " + repeat(input) +

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/bcftools/BcftoolsReheader.scala

+/**
+  * Biopet is built on top of GATK Queue for building bioinformatic
+  * pipelines. It is mainly intended to support LUMC SHARK cluster which is running
+  * SGE. But other types of HPC that are supported by GATK Queue (such as PBS)
+  * should also be able to execute Biopet tools and pipelines.
+  *
+  * Copyright 2014 Sequencing Analysis Support Core - Leiden University Medical Center
+  *
+  * Contact us at: sasc@lumc.nl
+  *
+  * A dual licensing mode is applied. The source code within this project is freely available for non-commercial use under an AGPL
+  * license; For commercial users or users who do not want to follow the AGPL
+  * license, please contact us to obtain a separate license.
+  */
+package nl.lumc.sasc.biopet.extensions.bcftools
+
+import java.io.File
+
+import nl.lumc.sasc.biopet.utils.config.Configurable
+import org.broadinstitute.gatk.utils.commandline.{Output, Input}
+
+class BcftoolsReheader(val parent: Configurable) extends Bcftools {
+
+  @Input(doc = "Input vcf file", required = false)
+  var input: File = _
+
+  @Input(doc = "File specifying how sample names should be renamed", required = true)
+  var renameSamples: File = _
+
+  @Output(doc = "Output vcf file", required = false)
+  var output: File = _
+
+  def cmdLine: String =
+    required(executable) +
+      required("reheader") +
+      required("--samples", renameSamples) +
+      optional("--output", output) +
+      optional(input)
+}
+
+object BcftoolsReheader {
+  def apply(parent: Configurable, renameSamples: File): BcftoolsReheader = {
+    val reheader = new BcftoolsReheader(parent)
+    reheader.renameSamples = renameSamples
+    reheader
+  }
+}

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/gatk/MuTect2.scala

@@ -14,20 +14,6 @@ class MuTect2(val parent: Configurable) extends CommandLineGATK with ScatterGath
 
   def analysis_type: String = "MuTect2"
 
-  /** Getter and setter for tumor sample bam file. */
-  def tumorSampleBam = input_file.find(file => isFileWithTag(file, "tumor")).getOrElse(null)
-  def tumorSampleBam_= (value:File):Unit = {
-    input_file = input_file.filterNot(file => isFileWithTag(file, "tumor"))
-    input_file :+= TaggedFile(value, "tumor")
-  }
-
-  /** Getter and setter for normal sample bam file. */
-  def normalSampleBam = input_file.find(file => isFileWithTag(file, "normal")).getOrElse(null)
-  def normalSampleBam_= (value:File):Unit = {
-    input_file = input_file.filterNot(file => isFileWithTag(file, "normal"))
-    input_file :+= TaggedFile(value, "normal")
-  }
-
   /** vcf file with info from cosmic db TODO desc  */
   @Input(fullName = "cosmic", shortName = "cosmic", required = false)
   var cosmic: Option[File] = config("cosmic")
@@ -46,7 +32,7 @@ class MuTect2(val parent: Configurable) extends CommandLineGATK with ScatterGath
   var contaminationFile: Option[File] = config("contamination_file")
 
   /** Output file of the program. */
-  @Output(fullName = "out", shortName = "o", required = true)
+  @Output(fullName = "out", shortName = "o", required = false)
   @Gather(classOf[CatVariantsGatherer])
   var outputVcf: File = _
 
@@ -197,7 +183,6 @@ class MuTect2(val parent: Configurable) extends CommandLineGATK with ScatterGath
 
 
   override def cmdLine = super.cmdLine +
-    required("--out", outputVcf) +
     optional("--cosmic", cosmic) +
     optional("--dbsnp", dbsnp) +
     optional("--normal_panel", ponFile) +
@@ -225,7 +210,8 @@ class MuTect2(val parent: Configurable) extends CommandLineGATK with ScatterGath
     conditional(annotateNDA, "--annotateNDA") +
     conditional(enableClusteredReadPositionFilter, "--enable_clustered_read_position_filter") +
     conditional(enableStrandArtifactFilter, "--enable_strand_artifact_filter") +
-    conditional(useNewAFCalculator, "--useNewAFCalculator")
+    conditional(useNewAFCalculator, "--useNewAFCalculator") +
+    (if (outputAsStdout) "" else required("--out", outputVcf))
 }
 
 object MuTect2 {

biopet-utils/src/main/scala/nl/lumc/sasc/biopet/utils/IoUtils.scala

@@ -93,4 +93,10 @@ object IoUtils {
       case e: IOException => false
     }
   }
+
+  def writeLinesToFile(output: File, lines: List[String]): Unit = {
+    val writer = new PrintWriter(output)
+    lines.foreach(writer.println(_))
+    writer.close()
+  }
 }

biopet-utils/src/main/scala/nl/lumc/sasc/biopet/utils/Logging.scala

@@ -38,7 +38,9 @@ object Logging {
 
   private[biopet] val errors: ListBuffer[Exception] = ListBuffer()
 
-  def addError(error: String, debug: Option[String] = None, cause: Option[Exception] = None): Unit = {
+  def addError(error: String,
+               debug: Option[String] = None,
+               cause: Option[Exception] = None): Unit = {
     val msg = error + debug.map("; " + _).getOrElse("")
 
     cause match {

shiva/src/main/scala/nl/lumc/sasc/biopet/pipelines/shiva/ShivaVariantcalling.scala

@@ -161,7 +161,7 @@ class ShivaVariantcalling(val parent: Configurable)
       caller.outputDir = new File(outputDir, caller.name)
       caller.genders = genders
       caller match {
-        case c:SomaticVariantCaller => c.tnPairs = tumorSamples
+        case c: SomaticVariantCaller => c.tnPairs = tumorSamples
         case _ =>
       }
 

shiva/src/main/scala/nl/lumc/sasc/biopet/pipelines/shiva/variantcallers/somatic/MuTect2.scala

 package nl.lumc.sasc.biopet.pipelines.shiva.variantcallers.somatic
 
-import nl.lumc.sasc.biopet.extensions.gatk.BqsrGather
-import nl.lumc.sasc.biopet.extensions.{Awk, Tabix, gatk}
-import nl.lumc.sasc.biopet.utils.config.{ConfigValue, Configurable}
+import nl.lumc.sasc.biopet.extensions.bcftools.BcftoolsReheader
+import nl.lumc.sasc.biopet.extensions.gatk.{BqsrGather, CombineVariants}
+import nl.lumc.sasc.biopet.extensions._
+import nl.lumc.sasc.biopet.utils.IoUtils
+import nl.lumc.sasc.biopet.utils.config.Configurable
+import org.broadinstitute.gatk.queue.extensions.gatk.TaggedFile
 
 class MuTect2(val parent: Configurable) extends SomaticVariantCaller {
 
@@ -10,8 +13,6 @@ class MuTect2(val parent: Configurable) extends SomaticVariantCaller {
 
   override val mergeVcfResults: Boolean = false
 
-  // currently not relevant, at the moment only one somatic variant caller exists in Biopet
-  // and results from this won't be merged together with the results from other methods
   def defaultPrio: Int = -1
 
   lazy val runConEst: Boolean = config("run_contest", default = false)
@@ -20,19 +21,35 @@ class MuTect2(val parent: Configurable) extends SomaticVariantCaller {
 
     val outputFiles = for (pair <- tnPairs) yield {
 
-      val outputFile = new File(outputDir, s"${pair.tumorSample}-${pair.normalSample}.$name.vcf.gz")
+      val bqsrFile =
+        if (inputBqsrFiles.contains(pair.tumorSample) &&
+            inputBqsrFiles.contains(pair.normalSample)) {
+          val gather = new BqsrGather()
+          gather.inputBqsrFiles =
+            List(inputBqsrFiles(pair.tumorSample), inputBqsrFiles(pair.normalSample))
+          gather.outputBqsrFile =
+            new File(outputDir, s"${pair.tumorSample}-${pair.normalSample}.bqsr")
+          add(gather)
 
-      val muTect2 = {
-        gatk.MuTect2(this, inputBams(pair.tumorSample), inputBams(pair.normalSample), outputFile)
-      }
+          Some(gather.outputBqsrFile)
+        } else None
+
+      val outputFile =
+        new File(outputDir, s"${pair.tumorSample}-${pair.normalSample}.$name.vcf.gz")
+
+      val muTect2 = new gatk.MuTect2(this)
+      muTect2.input_file :+= TaggedFile(inputBams(pair.tumorSample), "tumor")
+      muTect2.input_file :+= TaggedFile(inputBams(pair.normalSample), "normal")
+      muTect2.BQSR = bqsrFile
 
       if (runConEst) {
         val namePrefix = outputFile.getAbsolutePath.stripSuffix(".vcf.gz")
         val contEstOutput: File = new File(s"$namePrefix.contamination.txt")
         val contEst = gatk.ContEst(this,
-          inputBams(pair.tumorSample),
-          inputBams(pair.normalSample),
-          contEstOutput)
+                                   inputBams(pair.tumorSample),
+                                   inputBams(pair.normalSample),
+                                   contEstOutput)
+        contEst.BQSR = bqsrFile
         add(contEst)
 
         val contaminationPerSample: File = new File(s"$namePrefix.contamination.short.txt")
@@ -43,7 +60,8 @@ class MuTect2(val parent: Configurable) extends SomaticVariantCaller {
         muTect2.contaminationFile = Some(contaminationPerSample)
       }
 
-      if (inputBqsrFiles.contains(pair.tumorSample) && inputBqsrFiles.contains(pair.normalSample)) {
+      if (inputBqsrFiles.contains(pair.tumorSample) && inputBqsrFiles
+            .contains(pair.normalSample)) {
         val gather = new BqsrGather()
         gather.inputBqsrFiles =
           List(inputBqsrFiles(pair.tumorSample), inputBqsrFiles(pair.normalSample))
@@ -53,12 +71,26 @@ class MuTect2(val parent: Configurable) extends SomaticVariantCaller {
         muTect2.BQSR = Some(gather.outputBqsrFile)
       }
 
-      // TODO: Add name change
+      val renameFile = new File(outputDir, s".rename.${pair.tumorSample}-${pair.normalSample}.txt")
+      IoUtils.writeLinesToFile(renameFile,
+                               List(
+                                 s"TUMOR ${pair.tumorSample}",
+                                 s"NORMAL ${pair.normalSample}"
+                               ))
 
-      add(muTect2)
+      val pipe = muTect2 | BcftoolsReheader(this, renameFile) | new Bgzip(this) > outputFile
+      pipe.threadsCorrection = -2
+      add(pipe)
       add(Tabix(this, outputFile))
 
       outputFile
     }
+
+    if (outputFiles.size > 1) {
+      add(CombineVariants(this, outputFiles, outputFile))
+    } else if (outputFiles.nonEmpty) {
+      add(Ln(this, outputFiles.head, outputFile))
+      add(Ln(this, outputFiles.head + ".tbi", outputFile + ".tbi"))
+    }
   }
 }