Merge pull request #144 from biopet/refactorDelly

Refactor delly

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/delly/DellyCaller.scala → biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/delly/DellyCallerCall.scala

@@ -16,11 +16,18 @@ package nl.lumc.sasc.biopet.extensions.delly
 
 import java.io.File
 
-import nl.lumc.sasc.biopet.core.{Version, BiopetCommandLineFunction}
+import nl.lumc.sasc.biopet.core.{BiopetCommandLineFunction, Reference, Version}
 import nl.lumc.sasc.biopet.utils.config.Configurable
 import org.broadinstitute.gatk.utils.commandline.{Argument, Input, Output}
 
-class DellyCaller(val parent: Configurable) extends BiopetCommandLineFunction with Version {
+/**
+  *
+  * Created by imoustakas on 9-5-17.
+  */
+class DellyCallerCall(val parent: Configurable)
+    extends BiopetCommandLineFunction
+    with Version
+    with Reference {
   executable = config("exe", default = "delly")
 
   private lazy val versionexecutable: File = new File(executable)
@@ -29,33 +36,38 @@ class DellyCaller(val parent: Configurable) extends BiopetCommandLineFunction wi
   override def defaultCoreMemory = 4.0
 
   def versionCommand = versionexecutable.getAbsolutePath
-  def versionRegex = """DELLY \(Version: (.*)\)""".r
+  def versionRegex = """D(ELLY|elly) \(Version: (.*)\)""".r
   override def versionExitcode = List(0, 1)
-
   @Input(doc = "Input file (bam)")
   var input: File = _
 
-  @Output(doc = "Delly VCF output")
-  var outputvcf: File = _
+  @Input(doc = "Reference file", required = true)
+  var reference: File = _
 
-  @Argument(doc = "What kind of analysis to run: DEL,DUP,INV,TRA")
+  @Output(doc = "Delly BCF output")
+  var outputbcf: File = _
+
+  @Argument(doc = " SV type: DEL, DUP, INV, BND, INS")
   var analysistype: String = _
 
+  // Rest of available arguments for the call command of delly
+  var fileformat: Option[String] = config("fileformat")
+  var gsize: Option[Float] = config("gsize")
+  var keepdup: Boolean = config("keep-dup", default = false)
+  var buffersize: Option[Int] = config("buffer-size")
+  var outputdir: String = _
+
+  override def beforeGraph(): Unit = {
+    super.beforeGraph()
+    if (reference == null) reference = referenceFasta()
+  }
+
   def cmdLine =
     required(executable) +
-      required("-t", analysistype) +
-      required("-o", outputvcf) +
-      required(input) +
-      createEmptyOutputIfNeeded
-
-  // when no variants are found then the tool doesn't generate the output file either, in Biopet it's needed that the empty file would be there
-  private def createEmptyOutputIfNeeded =
-    s"""
-       |c=$$?
-       |if [ $$c -eq 0 ] && [ ! -f $outputvcf ]; then
-       |  echo '##fileformat=VCFv4.2' > $outputvcf
-       |  echo '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO' >> $outputvcf
-       |fi
-       |exit $$c""".stripMargin
+      required("call") +
+      required("-t", analysistype) + // SV type (DEL, DUP, INV, BND, INS)
+      required("-o", outputbcf) + // delly BCF output
+      required("-g", reference) + // reference file, on which the reads are aligned
+      required(input) // BAM file with the alignments
 
 }

shiva/src/main/scala/nl/lumc/sasc/biopet/pipelines/shiva/svcallers/Delly.scala

@@ -14,8 +14,8 @@
   */
 package nl.lumc.sasc.biopet.pipelines.shiva.svcallers
 
-import nl.lumc.sasc.biopet.extensions.delly.DellyCaller
-import nl.lumc.sasc.biopet.extensions.gatk.CatVariants
+import nl.lumc.sasc.biopet.extensions.bcftools.BcftoolsMerge
+import nl.lumc.sasc.biopet.extensions.delly.DellyCallerCall
 import nl.lumc.sasc.biopet.extensions.picard.SortVcf
 import nl.lumc.sasc.biopet.utils.Logging
 import nl.lumc.sasc.biopet.utils.config.Configurable
@@ -27,66 +27,80 @@ class Delly(val parent: Configurable) extends SvCaller {
   val del: Boolean = config("DEL", default = true)
   val dup: Boolean = config("DUP", default = true)
   val inv: Boolean = config("INV", default = true)
-  val tra: Boolean = config("TRA", default = true)
+  val bnd: Boolean = config("BND", default = true)
+  val ins: Boolean = config("INS", default = true)
 
   def biopetScript() {
+
     for ((sample, bamFile) <- inputBams) {
       val dellyDir = new File(outputDir, sample)
 
-      val catVariants = new CatVariants(this)
-      catVariants.outputFile = new File(dellyDir, sample + ".delly.vcf")
-      catVariants.isIntermediate = true
-      catVariants.writeHeaderToEmptyOutput = true
+      // Use bcftools merge to merge the bcf files. Output is an uncompressed vcf
+      val mergeVariants = new BcftoolsMerge(this)
+      mergeVariants.output = new File(dellyDir, sample + ".delly.vcf")
+      mergeVariants.m = Some("id")
+      mergeVariants.forcesamples = true
 
       if (del) {
-        val delly = new DellyCaller(this)
+        val delly = new DellyCallerCall(this)
         delly.input = bamFile
         delly.analysistype = "DEL"
         delly.isIntermediate = true
-        delly.outputvcf = new File(dellyDir, sample + ".delly.del.vcf")
+        delly.outputbcf = new File(dellyDir, sample + ".delly.del.bcf")
         add(delly)
 
-        catVariants.variant :+= delly.outputvcf
+        // bcf files must to be concatenated with bcftools merge
+        mergeVariants.input :+= delly.outputbcf
       }
       if (dup) {
-        val delly = new DellyCaller(this)
+        val delly = new DellyCallerCall(this)
         delly.input = bamFile
         delly.analysistype = "DUP"
         delly.isIntermediate = true
-        delly.outputvcf = new File(dellyDir, sample + ".delly.dup.vcf")
+        delly.outputbcf = new File(dellyDir, sample + ".delly.dup.bcf")
         add(delly)
 
-        catVariants.variant :+= delly.outputvcf
+        mergeVariants.input :+= delly.outputbcf
       }
       if (inv) {
-        val delly = new DellyCaller(this)
+        val delly = new DellyCallerCall(this)
         delly.input = bamFile
         delly.analysistype = "INV"
         delly.isIntermediate = true
-        delly.outputvcf = new File(dellyDir, sample + ".delly.inv.vcf")
+        delly.outputbcf = new File(dellyDir, sample + ".delly.inv.bcf")
+        add(delly)
+
+        mergeVariants.input :+= delly.outputbcf
+      }
+      if (ins) {
+        val delly = new DellyCallerCall(this)
+        delly.input = bamFile
+        delly.analysistype = "INS"
+        delly.isIntermediate = true
+        delly.outputbcf = new File(dellyDir, sample + ".delly.ins.bcf")
         add(delly)
 
-        catVariants.variant :+= delly.outputvcf
+        mergeVariants.input :+= delly.outputbcf
       }
-      if (tra) {
-        val delly = new DellyCaller(this)
+      if (bnd) {
+        val delly = new DellyCallerCall(this)
         delly.input = bamFile
-        delly.analysistype = "TRA"
+        delly.analysistype = "BND"
         delly.isIntermediate = true
-        delly.outputvcf = new File(dellyDir, sample + ".delly.tra.vcf")
+        delly.outputbcf = new File(dellyDir, sample + ".delly.tra.bcf")
         add(delly)
 
-        catVariants.variant :+= delly.outputvcf
+        mergeVariants.input :+= delly.outputbcf
       }
 
-      if (catVariants.variant.isEmpty)
+      if (mergeVariants.input.isEmpty)
         Logging.addError("At least 1 SV-type should be selected for Delly")
 
-      add(catVariants)
+      add(mergeVariants)
 
       val compressedVCF = new SortVcf(this)
-      compressedVCF.input = catVariants.outputFile
-      compressedVCF.output = catVariants.outputFile + ".gz"
+      compressedVCF.input = mergeVariants.output
+      compressedVCF.output = mergeVariants.output + ".gz"
       add(compressedVCF)
 
       addVCF(sample, compressedVCF.output)

shiva/src/test/scala/nl/lumc/sasc/biopet/pipelines/shiva/ShivaSvCallingTest.scala

@@ -24,7 +24,7 @@ import nl.lumc.sasc.biopet.extensions.breakdancer.{
   BreakdancerVCF
 }
 import nl.lumc.sasc.biopet.extensions.clever.CleverCaller
-import nl.lumc.sasc.biopet.extensions.delly.DellyCaller
+import nl.lumc.sasc.biopet.extensions.delly.DellyCallerCall
 import nl.lumc.sasc.biopet.extensions.pindel.{PindelCaller, PindelConfig, PindelVCF}
 import nl.lumc.sasc.biopet.utils.{ConfigUtils, Logging}
 import nl.lumc.sasc.biopet.utils.config.Config
@@ -117,7 +117,8 @@ class ShivaSvCallingTest extends TestNGSuite with Matchers {
       pipeline.functions.count(_.isInstanceOf[PindelVCF]) shouldBe (if (pindel) bams else 0)
 
       pipeline.functions.count(_.isInstanceOf[CleverCaller]) shouldBe (if (clever) bams else 0)
-      pipeline.functions.count(_.isInstanceOf[DellyCaller]) shouldBe (if (delly) bams * 4 else 0)
+      pipeline.functions.count(_.isInstanceOf[DellyCallerCall]) shouldBe (if (delly) bams * 5
+                                                                          else 0)
 
     }
   }
@@ -128,30 +129,37 @@ class ShivaSvCallingTest extends TestNGSuite with Matchers {
     for (del <- bool;
          dup <- bool;
          inv <- bool;
-         tra <- bool) yield Array(1, del, dup, inv, tra)
+         bnd <- bool;
+         ins <- bool) yield Array(1, del, dup, inv, bnd, ins)
   }
 
   @Test(dataProvider = "dellyOptions")
-  def testShivaDelly(bams: Int, del: Boolean, dup: Boolean, inv: Boolean, tra: Boolean): Unit = {
+  def testShivaDelly(bams: Int,
+                     del: Boolean,
+                     dup: Boolean,
+                     inv: Boolean,
+                     bnd: Boolean,
+                     ins: Boolean): Unit = {
     val outputDir = ShivaSvCallingTest.outputDir
     dirs :+= outputDir
 
     val map = Map("sv_callers" -> List("delly"),
                   "delly" ->
-                    Map("DEL" -> del, "DUP" -> dup, "INV" -> inv, "TRA" -> tra))
+                    Map("DEL" -> del, "DUP" -> dup, "INV" -> inv, "BND" -> bnd, "INS" -> ins))
     val pipeline = initPipeline(map, outputDir)
 
     pipeline.inputBams = Map("bam" -> ShivaSvCallingTest.inputTouch("bam" + ".bam"))
 
-    if (!del && !dup && !inv && !tra) intercept[IllegalStateException] {
+    if (!del && !dup && !inv && !bnd && !ins) intercept[IllegalStateException] {
       pipeline.init()
       pipeline.script()
     } else {
       pipeline.init()
       pipeline.script()
 
-      pipeline.functions.count(_.isInstanceOf[DellyCaller]) shouldBe
-        ((if (del) 1 else 0) + (if (dup) 1 else 0) + (if (inv) 1 else 0) + (if (tra) 1 else 0))
+      pipeline.functions.count(_.isInstanceOf[DellyCallerCall]) shouldBe
+        ((if (del) 1 else 0) + (if (dup) 1 else 0) + (if (inv) 1 else 0) +
+          (if (bnd) 1 else 0) + (if (ins) 1 else 0))
     }
   }
 
@@ -221,6 +229,7 @@ object ShivaSvCallingTest {
   copyFile("ref.dict")
   copyFile("ref.fa.fai")
 
+  // A dummy config file. Tools being used by this test must be included
   def config(outputDir: File) = Map(
     "skip_write_dependencies" -> true,
     "name_prefix" -> "test",
@@ -231,6 +240,7 @@ object ShivaSvCallingTest {
     "reference_fasta" -> (inputDir + File.separator + "ref.fa"),
     "gatk_jar" -> "test",
     "samtools" -> Map("exe" -> "test"),
+    "bcftools" -> Map("exe" -> "test"),
     "md5sum" -> Map("exe" -> "test"),
     "bgzip" -> Map("exe" -> "test"),
     "tabix" -> Map("exe" -> "test"),