Adding bwa index

public/biopet-framework/src/main/scala/nl/lumc/sasc/biopet/extensions/bwa/BwaIndex.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 that are
+ * not part of GATK Queue 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.bwa
+
+import java.io.File
+
+import nl.lumc.sasc.biopet.core.Reference
+import nl.lumc.sasc.biopet.core.config.Configurable
+import org.broadinstitute.gatk.utils.commandline.{Input, Output}
+
+/**
+ * Extension for bwa aln
+ *
+ * Based on version 0.7.12-r1039
+ *
+ * Created by pjvan_thof on 1/16/15.
+ */
+class BwaIndex(val root: Configurable) extends Bwa {
+  @Input(doc = "Fastq file", required = true)
+  var reference: File = _
+
+  @Output(doc = "Index files for bwa", required = false)
+  private var output: List[File] = Nil
+
+  var a: Option[String] = config("a", freeVar = false)
+  var p: Option[String] = config("p", freeVar = false)
+  var b: Option[Int] = config("e", freeVar = false)
+  var _6: Boolean = config("6", default = false, freeVar = false)
+
+  override def defaultCoreMemory = 4.0
+
+  override def beforeGraph() {
+    super.beforeGraph()
+    List(".sa", ".pac")
+      .foreach(ext => output ::= new File(reference.getAbsolutePath + ext))
+    output = output.distinct
+  }
+
+  /** Returns command to execute */
+  def cmdLine = required(executable) +
+    required("index") +
+    optional("-a", a) +
+    optional("-p", p) +
+    optional("-b", b) +
+    conditional(_6, "-6") +
+    required(reference)
+}

public/biopet-framework/src/main/scala/nl/lumc/sasc/biopet/extensions/samtools/SamtoolsFaidx.scala

@@ -31,6 +31,7 @@ class SamtoolsFaidx(val root: Configurable) extends Samtools {
   def output = _output
 
   override def beforeGraph: Unit = {
+    super.beforeGraph
     _output = new File(input.getParentFile, input.getName + ".fai")
   }
 

public/biopet-framework/src/main/scala/nl/lumc/sasc/biopet/pipelines/GenerateIndexes.scala

@@ -19,9 +19,10 @@ import java.io.File
 
 import nl.lumc.sasc.biopet.core.{PipelineCommand, BiopetQScript}
 import nl.lumc.sasc.biopet.core.config.Configurable
+import nl.lumc.sasc.biopet.extensions.bwa.BwaIndex
 import nl.lumc.sasc.biopet.extensions.picard.CreateSequenceDictionary
 import nl.lumc.sasc.biopet.extensions.samtools.SamtoolsFaidx
-import nl.lumc.sasc.biopet.extensions.{Md5sum, Zcat, Curl}
+import nl.lumc.sasc.biopet.extensions.{Ln, Md5sum, Zcat, Curl}
 import nl.lumc.sasc.biopet.utils.ConfigUtils
 import org.broadinstitute.gatk.queue.QScript
 import org.broadinstitute.gatk.utils.commandline
@@ -47,15 +48,15 @@ class GenerateIndexes(val root: Configurable) extends QScript with BiopetQScript
       val speciesDir = new File(outputDir, speciesName)
       for ((genomeName, c) <- speciesConfig) {
         val genomeConfig = ConfigUtils.any2map(c)
-        val fastaUrl = genomeConfig.getOrElse("fasta_url",
-          throw new IllegalArgumentException(s"No fasta_url found for $speciesName - $genomeName")).toString
+        val fastaUri = genomeConfig.getOrElse("fasta_uri",
+          throw new IllegalArgumentException(s"No fasta_uri found for $speciesName - $genomeName")).toString
 
         val genomeDir = new File(speciesDir, genomeName)
         val fastaFile = new File(genomeDir, "reference.fa")
 
         val curl = new Curl(this)
-        curl.url = fastaUrl
-        if (fastaUrl.endsWith(".gz")) {
+        curl.url = fastaUri
+        if (fastaUri.endsWith(".gz")) {
           curl.output = new File(genomeDir, "reference.fa.gz")
           curl.isIntermediate = true
           add(Zcat(this, curl.output, fastaFile))
@@ -72,9 +73,28 @@ class GenerateIndexes(val root: Configurable) extends QScript with BiopetQScript
         createDict.output = new File(genomeDir, fastaFile.getName.stripSuffix(".fa") + ".dict")
         createDict.species = Some(speciesName)
         createDict.genomeAssembly = Some(genomeName)
-        createDict.uri = Some(fastaUrl)
+        createDict.uri = Some(fastaUri)
         add(createDict)
 
+        def createLinks(dir: File): File = {
+          val newFastaFile = new File(dir, fastaFile.getName)
+          val newFai = new File(dir, faidx.output.getName)
+          val newDict = new File(dir, createDict.output.getName)
+
+          add(Ln(this, faidx.output, newFai))
+          add(Ln(this, createDict.output, newDict))
+          val lnFasta = Ln(this, fastaFile, newFastaFile)
+          lnFasta.deps ++= List(newFai, newDict)
+          add(lnFasta)
+          newFastaFile
+        }
+
+        // Bwa index
+
+        val bwaIndex = new BwaIndex(this)
+        bwaIndex.reference = createLinks(new File(genomeDir, "bwa"))
+        add(bwaIndex)
+
         //TODO: other indexes
       }
     }