Make peakcalling from bam files possible

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/macs2/Macs2CallPeak.scala

@@ -49,7 +49,7 @@ class Macs2CallPeak(val parent: Configurable) extends Macs2 {
   var gsize: Option[Float] = config("gsize")
   var keepdup: Boolean = config("keep-dup", default = false)
   var buffersize: Option[Int] = config("buffer-size")
-  var outputdir: String = _
+  var outputdir: File = _
   var name: Option[String] = config("name")
   var bdg: Boolean = config("bdg", default = false)
   var verbose: Boolean = config("verbose", default = false)
@@ -76,12 +76,12 @@ class Macs2CallPeak(val parent: Configurable) extends Macs2 {
   override def beforeGraph(): Unit = {
     if (name.isEmpty) throw new IllegalArgumentException("Name is not defined")
     if (outputdir == null) throw new IllegalArgumentException("Outputdir is not defined")
-    outputNarrow = new File(outputdir + name.get + ".narrowPeak")
-    outputBroad = new File(outputdir + name.get + ".broadPeak")
-    outputXls = new File(outputdir + name.get + ".xls")
-    outputBdg = new File(outputdir + name.get + ".bdg")
-    outputR = new File(outputdir + name.get + ".r")
-    outputGapped = new File(outputdir + name.get + ".gappedPeak")
+    outputNarrow = new File(outputdir, name.get + ".narrowPeak")
+    outputBroad = new File(outputdir, name.get + ".broadPeak")
+    outputXls = new File(outputdir, name.get + ".xls")
+    outputBdg = new File(outputdir, name.get + ".bdg")
+    outputR = new File(outputdir, name.get + ".r")
+    outputGapped = new File(outputdir, name.get + ".gappedPeak")
   }
 
   /** Returns command to execute */

carp/src/main/scala/nl/lumc/sasc/biopet/pipelines/carp/Carp.scala

@@ -14,11 +14,8 @@
   */
 package nl.lumc.sasc.biopet.pipelines.carp
 
-import java.io.File
-
 import nl.lumc.sasc.biopet.core._
 import nl.lumc.sasc.biopet.core.report.ReportBuilderExtension
-import nl.lumc.sasc.biopet.extensions.macs2.Macs2CallPeak
 import nl.lumc.sasc.biopet.extensions.picard.BuildBamIndex
 import nl.lumc.sasc.biopet.extensions.samtools.SamtoolsView
 import nl.lumc.sasc.biopet.pipelines.bammetrics.BamMetrics
@@ -37,7 +34,7 @@ class Carp(val parent: Configurable) extends QScript with MultisampleMappingTrai
   qscript =>
   def this() = this(null)
 
-  override def defaults = super.defaults ++ Map(
+  override def defaults: Map[String, Any] = super.defaults ++ Map(
     "mapping" -> Map(
       "skip_markduplicates" -> false,
       "aligner" -> "bwa-mem"
@@ -46,7 +43,7 @@ class Carp(val parent: Configurable) extends QScript with MultisampleMappingTrai
     "samtoolsview" -> Map("q" -> 10)
   )
 
-  override def fixedValues = super.fixedValues ++ Map(
+  override def fixedValues: Map[String, Any] = super.fixedValues ++ Map(
     "samtoolsview" -> Map(
       "h" -> true,
       "b" -> true
@@ -63,7 +60,8 @@ class Carp(val parent: Configurable) extends QScript with MultisampleMappingTrai
 
     val controls: List[String] = config("control", default = Nil)
 
-    override def summarySettings = super.summarySettings ++ Map("controls" -> controls)
+    override def summarySettings: Map[String, Any] =
+      super.summarySettings ++ Map("controls" -> controls)
 
     override def addJobs(): Unit = {
       super.addJobs()
@@ -90,13 +88,6 @@ class Carp(val parent: Configurable) extends QScript with MultisampleMappingTrai
       buildBamIndex.output =
         swapExt(preProcessBam.get.getParentFile, preProcessBam.get, ".bam", ".bai")
       add(buildBamIndex)
-
-      val macs2 = new Macs2CallPeak(qscript)
-      macs2.treatment = preProcessBam.get
-      macs2.name = Some(sampleId)
-      macs2.outputdir = sampleDir + File.separator + "macs2" + File.separator + sampleId + File.separator
-      macs2.fileformat = if (paired) Some("BAMPE") else Some("BAM")
-      add(macs2)
     }
   }
 
@@ -116,29 +107,22 @@ class Carp(val parent: Configurable) extends QScript with MultisampleMappingTrai
     p
   }
 
-  override def init() = {
+  override def init(): Unit = {
     super.init()
     // ensure that no samples are called 'control' since that is our reserved keyword
     require(!sampleIds.contains("control"),
             "No sample should be named 'control' since it is a reserved for the Carp pipeline")
+    peakCalling.controls = samples.map(x => x._1 -> x._2.controls)
+    peakCalling.outputDir = new File(outputDir, "peak_calling")
   }
 
+  lazy val peakCalling = new PeakCalling(this)
+
   override def addMultiSampleJobs(): Unit = {
     super.addMultiSampleJobs()
-    for ((sampleId, sample) <- samples) {
-      for (controlId <- sample.controls) {
-        if (!samples.contains(controlId))
-          throw new IllegalStateException(
-            "For sample: " + sampleId + " this control: " + controlId + " does not exist")
-        val macs2 = new Macs2CallPeak(this)
-        macs2.treatment = sample.preProcessBam.get
-        macs2.control = samples(controlId).preProcessBam.get
-        macs2.name = Some(sampleId + "_VS_" + controlId)
-        macs2.fileformat = if (paired) Some("BAMPE") else Some("BAM")
-        macs2.outputdir = sample.sampleDir + File.separator + "macs2" + File.separator + macs2.name.get + File.separator
-        add(macs2)
-      }
-    }
+    peakCalling.inputBams = samples.map(x => x._1 -> x._2.preProcessBam.get)
+    peakCalling.paired = paired
+    add(peakCalling)
   }
 }
 

carp/src/main/scala/nl/lumc/sasc/biopet/pipelines/carp/PeakCalling.scala

+package nl.lumc.sasc.biopet.pipelines.carp
+
+import java.io.File
+
+import nl.lumc.sasc.biopet.core.{PipelineCommand, Reference}
+import nl.lumc.sasc.biopet.core.summary.SummaryQScript
+import nl.lumc.sasc.biopet.extensions.macs2.Macs2CallPeak
+import nl.lumc.sasc.biopet.utils.{BamUtils, ConfigUtils}
+import nl.lumc.sasc.biopet.utils.config.Configurable
+import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb
+import org.broadinstitute.gatk.queue.QScript
+
+import scala.concurrent.Await
+import scala.concurrent.duration.Duration
+import scala.concurrent.ExecutionContext.Implicits.global
+
+/**
+  * Created by pjvanthof on 29/05/2017.
+  */
+class PeakCalling(val parent: Configurable) extends QScript with SummaryQScript with Reference {
+
+  def this() = this(null)
+
+  @Input(doc = "Bam files (should be deduped bams)", shortName = "BAM", required = true)
+  protected[carp] var inputBamsArg: List[File] = Nil
+
+  var inputBams: Map[String, File] = Map()
+
+  /** Must return a map with used settings for this pipeline */
+  override def summarySettings: Map[String, Any] = Map()
+
+  /** File to put in the summary for thie pipeline */
+  override def summaryFiles: Map[String, File] = Map()
+
+  var paired: Boolean = config("paired", default = false)
+
+  var controls: Map[String, List[String]] = Map()
+
+  def sampleDir(sampleName: String): File = new File(outputDir, sampleName)
+
+  /** Init for pipeline */
+  override def init(): Unit = {
+    if (controls.isEmpty) controls = config("controls", default = Map()).asMap.map(x => x._1 -> ConfigUtils.any2stringList(x._2))
+    if (inputBamsArg.nonEmpty) {
+      inputBams = BamUtils.sampleBamMap(inputBamsArg)
+
+      val db = SummaryDb.openSqliteSummary(summaryDbFile)
+      for (sampleName <- inputBams.keys) {
+        if (Await.result(db.getSampleId(summaryRunId, sampleName), Duration.Inf).isEmpty) {
+          db.createSample(sampleName, summaryRunId)
+        }
+      }
+    }
+  }
+
+  /** Pipeline itself */
+  override def biopetScript(): Unit = {
+    for ((sampleName, bamFile) <- inputBams) {
+      val macs2 = new Macs2CallPeak(this)
+      macs2.treatment = bamFile
+      macs2.name = Some(sampleName)
+      macs2.outputdir = new File(sampleDir(sampleName), "single_sample_calls")
+      macs2.fileformat = if (paired) Some("BAMPE") else Some("BAM")
+      add(macs2)
+
+      for (control <- controls.getOrElse(sampleName, Nil)) {
+        val macs2 = new Macs2CallPeak(this)
+        macs2.treatment = bamFile
+        macs2.control = inputBams.getOrElse(control, throw new IllegalStateException(s"Control '$control' is not found"))
+        macs2.name = Some(sampleName + "_VS_" + control)
+        macs2.fileformat = if (paired) Some("BAMPE") else Some("BAM")
+        macs2.outputdir = new File(sampleDir(sampleName), sampleName + "_VS_" + control)
+        add(macs2)
+      }
+
+    }
+  }
+}
+
+object PeakCalling extends PipelineCommand
\ No newline at end of file