Merge branch 'Fix-BIOPET-479' of github.com:biopet/biopet into Fix-BIOPET-479

basty/src/test/scala/nl/lumc/sasc/biopet/pipelines/basty/BastyTest.scala

@@ -98,7 +98,7 @@ class BastyTest extends TestNGSuite with Matchers {
       val numberLibs = (if (sample1) 1 else 0) + (if (sample2) 2 else 0)
       val numberSamples = (if (sample1) 1 else 0) + (if (sample2) 1 else 0)
 
-      pipeline.functions.count(_.isInstanceOf[MarkDuplicates]) shouldBe (numberLibs + (if (sample2) 1 else 0))
+      pipeline.functions.count(_.isInstanceOf[MarkDuplicates]) shouldBe (numberLibs + numberSamples)
 
       // Gatk preprocess
       pipeline.functions.count(_.isInstanceOf[IndelRealigner]) shouldBe (numberLibs * (if (realign) 1 else 0) + (if (sample2 && realign) 1 else 0))

biopet-core/src/main/resources/nl/lumc/sasc/biopet/tools/plotHeatmap.R

@@ -14,6 +14,7 @@ rownames(heat) <- heat[,1]
 heat<- heat[,-1]
 heat<- as.matrix(heat)
 
+textMargin <- max(sapply(rownames(heat), nchar)) + 4
 colNumber <- 50
 col <- rev(colorRampPalette(brewer.pal(11, "Spectral"))(colNumber))
 for (i in (colNumber+1):(colNumber+round((dist(range(heat)) - dist(range(heat[heat < 1]))) / dist(range(heat[heat < 1])) * colNumber))) {
@@ -22,7 +23,7 @@ for (i in (colNumber+1):(colNumber+round((dist(range(heat)) - dist(range(heat[he
 col[length(col)] <- "#00FF00"
 
 png(file = outputArg, width = 1200, height = 1200)
-heatmap.2(heat, trace = 'none', col = col, Colv=NA, Rowv=NA, dendrogram="none", margins = c(12, 12), na.color="#00FF00")
+heatmap.2(heat, trace = 'none', col = col, Colv=NA, Rowv=NA, dendrogram="none", margins = c(textMargin, textMargin), na.color="#00FF00")
 dev.off()
 
 hc <- hclust(d = dist(heat))
@@ -31,5 +32,5 @@ plot(as.dendrogram(hc), horiz=TRUE, asp=0.02)
 dev.off()
 
 png(file = outputArgClustering, width = 1200, height = 1200)
-heatmap.2(heat, trace = 'none', col = col, Colv="Rowv", dendrogram="row",margins = c(12, 12), na.color="#00FF00")
+heatmap.2(heat, trace = 'none', col = col, Colv="Rowv", dendrogram="row",margins = c(textMargin, textMargin), na.color="#00FF00")
 dev.off()

biopet-core/src/main/scala/nl/lumc/sasc/biopet/core/BiopetCommandLineFunction.scala

@@ -54,7 +54,7 @@ trait BiopetCommandLineFunction extends CommandLineResources { biopetFunction =>
     }
     writer.println("set -eubf")
     writer.println("set -o pipefail")
-    writer.println(this.commandLine)
+    writer.println(lines.mkString("\n"))
     jobDelayTime.foreach(x => writer.println(s"sleep $x"))
     writer.close()
   }
@@ -99,7 +99,13 @@ trait BiopetCommandLineFunction extends CommandLineResources { biopetFunction =>
     beforeGraph()
     internalBeforeGraph()
 
-    if (jobOutputFile != null) this.commandDirectory = this.jobOutputFile.getAbsoluteFile.getParentFile
+    if (jobOutputFile != null) {
+      this.commandDirectory = this.jobOutputFile.getAbsoluteFile.getParentFile
+      this match {
+        case s: ScatterGatherableFunction => s.scatterGatherDirectory = new File(this.commandDirectory, ".scatter")
+        case _                            =>
+      }
+    }
 
     super.freezeFieldValues()
   }

biopet-core/src/main/scala/nl/lumc/sasc/biopet/core/BiopetQScript.scala

@@ -126,17 +126,21 @@ trait BiopetQScript extends Configurable with GatkLogging { qscript: QScript =>
       }
     })
 
-    logger.info("Adding report")
-    this match {
-      case q: MultiSampleQScript if q.onlySamples.nonEmpty && !q.samples.forall(x => q.onlySamples.contains(x._1)) =>
-        logger.info("Write report is skipped because sample flag is used")
-      case _ => reportClass.foreach { report =>
-        for (f <- functions) f match {
-          case w: WriteSummary => report.deps :+= w.jobOutputFile
-          case _               =>
+    val writeHtmlReport: Boolean = config("write_html_report", default = true)
+
+    if (writeHtmlReport) {
+      logger.info("Adding report")
+      this match {
+        case q: MultiSampleQScript if q.onlySamples.nonEmpty && !q.samples.forall(x => q.onlySamples.contains(x._1)) =>
+          logger.info("Write report is skipped because sample flag is used")
+        case _ => reportClass.foreach { report =>
+          for (f <- functions) f match {
+            case w: WriteSummary => report.deps :+= w.jobOutputFile
+            case _               =>
+          }
+          report.jobOutputFile = new File(report.outputDir, ".report.out")
+          add(report)
         }
-        report.jobOutputFile = new File(report.outputDir, ".report.out")
-        add(report)
       }
     }
 

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/bwa/BwaMem.scala

@@ -33,7 +33,7 @@ class BwaMem(val parent: Configurable) extends Bwa with Reference {
   var R2: File = _
 
   @Input(doc = "The reference file for the bam files.", shortName = "R")
-  var reference: File = null
+  var reference: File = _
 
   @Output(doc = "Output file SAM", shortName = "output")
   var output: File = _
@@ -42,7 +42,7 @@ class BwaMem(val parent: Configurable) extends Bwa with Reference {
   var k: Option[Int] = config("k")
   var r: Option[Float] = config("r")
   var S: Boolean = config("S", default = false)
-  var M: Boolean = config("M", default = true)
+  var M: Boolean = config("M", default = false)
   var w: Option[Int] = config("w")
   var d: Option[Int] = config("d")
   var c: Option[Int] = config("c")
@@ -75,7 +75,7 @@ class BwaMem(val parent: Configurable) extends Bwa with Reference {
     if (reference == null) reference = referenceFasta()
   }
 
-  def cmdLine = {
+  def cmdLine: String = {
     required(executable) +
       required("mem") +
       optional("-k", k) +

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/clever/CleverFixVCF.scala

@@ -38,6 +38,8 @@ class CleverFixVCF(val parent: Configurable) extends BiopetJavaCommandLineFuncti
   @Argument(doc = "Samplename")
   var sampleName: String = _
 
+  override def defaultCoreMemory = 4.0
+
   override def cmdLine = super.cmdLine +
     required("-i", input) +
     required("-o", output) +

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

@@ -56,6 +56,8 @@ class HaplotypeCaller(val parent: Configurable) extends CommandLineGATK with Sca
   @Argument(fullName = "maxGGAARExtension", shortName = "maxGGAARExtension", doc = "the maximum extent into the full active region extension that we're willing to go in genotyping our events for GGA mode", required = false, exclusiveOf = "", validation = "")
   var maxGGAARExtension: Option[Int] = config("maxGGAARExtension")
 
+  var useNewAFCalculator: Boolean = config("useNewAFCalculator", default = false)
+
   /** Include at least this many bases around an event for calling indels */
   @Argument(fullName = "paddingAroundIndels", shortName = "paddingAroundIndels", doc = "Include at least this many bases around an event for calling indels", required = false, exclusiveOf = "", validation = "")
   var paddingAroundIndels: Option[Int] = config("paddingAroundIndels")
@@ -441,6 +443,7 @@ class HaplotypeCaller(val parent: Configurable) extends CommandLineGATK with Sca
     optional("-paddingAroundSNPs", paddingAroundSNPs, spaceSeparated = true, escape = true, format = "%s") +
     repeat("-comp", comp, formatPrefix = TaggedFile.formatCommandLineParameter, spaceSeparated = true, escape = true, format = "%s") +
     repeat("-A", annotation, spaceSeparated = true, escape = true, format = "%s") +
+    conditional(useNewAFCalculator, "--useNewAFCalculator") +
     repeat("-XA", excludeAnnotation, spaceSeparated = true, escape = true, format = "%s") +
     repeat("-G", group, spaceSeparated = true, escape = true, format = "%s") +
     conditional(debug, "-debug", escape = true, format = "%s") +

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/sambamba/Sambamba.scala

@@ -25,6 +25,6 @@ abstract class Sambamba extends BiopetCommandLineFunction with Version {
 
   executable = config("exe", default = "sambamba", namespace = "sambamba", freeVar = false)
   def versionCommand = executable
-  def versionRegex = """sambamba v(.*)""".r
+  def versionRegex = """sambamba v?(.*)""".r
   override def versionExitcode = List(0, 1)
 }
\ No newline at end of file

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/sambamba/SambambaMarkdup.scala

@@ -21,7 +21,6 @@ import org.broadinstitute.gatk.utils.commandline.{ Input, Output }
 
 /** Extension for sambemba markdup  */
 class SambambaMarkdup(val parent: Configurable) extends Sambamba {
-  override def defaultThreads = 4
 
   @Input(doc = "Bam File")
   var input: File = _
@@ -32,17 +31,29 @@ class SambambaMarkdup(val parent: Configurable) extends Sambamba {
   var removeDuplicates: Boolean = config("remove_duplicates", default = false)
 
   // @doc: compression_level 6 is average, 0 = no compression, 9 = best
-  val compressionLevel: Option[Int] = config("compression_level", default = 6)
-  val hashTableSize: Option[Int] = config("hash-table-size", default = 262144)
-  val overflowListSize: Option[Int] = config("overflow-list-size", default = 200000)
-  val ioBufferSize: Option[Int] = config("io-buffer-size", default = 128)
+  val compressionLevel: Option[Int] = config("compression_level")
+  val hashTableSize: Option[Int] = config("hash-table-size")
+  val overflowListSize: Option[Int] = config("overflow-list-size")
+  val ioBufferSize: Option[Int] = config("io-buffer-size")
+  val showProgress: Boolean = config("show-progress", default = true)
+
+  override def defaultThreads = 4
+  override def defaultCoreMemory = 4.0
+
+  @Output
+  private var indexOutput: File = _
+
+  override def beforeGraph(): Unit = {
+    indexOutput = new File(output + ".bai")
+  }
 
   /** Returns command to execute */
-  def cmdLine = required(executable) +
+  def cmdLine: String = required(executable) +
     required("markdup") +
     conditional(removeDuplicates, "--remove-duplicates") +
     optional("-t", nCoresRequest) +
     optional("-l", compressionLevel) +
+    conditional(showProgress, "--show-progress") +
     optional("--hash-table-size=", hashTableSize, spaceSeparated = false) +
     optional("--overflow-list-size=", overflowListSize, spaceSeparated = false) +
     optional("--io-buffer-size=", ioBufferSize, spaceSeparated = false) +
@@ -51,10 +62,11 @@ class SambambaMarkdup(val parent: Configurable) extends Sambamba {
 }
 
 object SambambaMarkdup {
-  def apply(root: Configurable, input: File, output: File): SambambaMarkdup = {
+  def apply(root: Configurable, input: File, output: File, isIntermediate: Boolean = false): SambambaMarkdup = {
     val markdup = new SambambaMarkdup(root)
     markdup.input = input
     markdup.output = output
+    markdup.isIntermediate = isIntermediate
     markdup
   }
 

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/sambamba/SambambaMerge.scala

@@ -21,7 +21,6 @@ import org.broadinstitute.gatk.utils.commandline.{ Input, Output }
 
 /** Extension for sambemba merge  */
 class SambambaMerge(val parent: Configurable) extends Sambamba {
-  override def defaultThreads = 4
 
   @Input(doc = "Bam File[s]")
   var input: List[File] = Nil
@@ -30,13 +29,29 @@ class SambambaMerge(val parent: Configurable) extends Sambamba {
   var output: File = _
 
   // @doc: compression_level 6 is average, 0 = no compression, 9 = best
-  val compressionLevel: Option[Int] = config("compression_level", default = 6)
+  val compressionLevel: Option[Int] = config("compression_level")
+  val header: Boolean = config("header", default = false)
+  val showProgress: Boolean = config("show-progress", default = true)
+  val filter: Option[String] = config("filter")
+
+  override def defaultThreads = 4
+  override def defaultCoreMemory = 4.0
+
+  @Output
+  private var indexOutput: File = _
+
+  override def beforeGraph(): Unit = {
+    indexOutput = new File(output + ".bai")
+  }
 
   /** Returns command to execute */
-  def cmdLine = required(executable) +
+  def cmdLine: String = required(executable) +
     required("merge") +
     optional("-t", nCoresRequest) +
     optional("-l", compressionLevel) +
+    optional("-F", filter) +
+    conditional(header, "--header") +
+    conditional(showProgress, "--show-progress") +
     required(output) +
-    repeat("", input)
+    repeat(input)
 }

biopet-tools-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/tools/VcfStatsForSv.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.tools
+
+import java.io.File
+
+import nl.lumc.sasc.biopet.core.summary.Summarizable
+import nl.lumc.sasc.biopet.core.ToolCommandFunction
+import nl.lumc.sasc.biopet.utils.config.Configurable
+import nl.lumc.sasc.biopet.utils.ConfigUtils
+import org.broadinstitute.gatk.utils.commandline.{ Argument, Input, Output }
+
+class VcfStatsForSv(val parent: Configurable) extends ToolCommandFunction with Summarizable {
+  def toolObject = nl.lumc.sasc.biopet.tools.vcfstats.VcfStatsForSv
+
+  mainFunction = false
+
+  @Input(required = true)
+  var inputFile: File = _
+
+  @Argument(required = true)
+  var histogramBinBoundaries: Array[Int] = _
+
+  @Output(required = true)
+  var outputFile: File = _
+
+  override def defaultCoreMemory = 1.0
+
+  override def cmdLine = super.cmdLine +
+    required("-i", inputFile) +
+    required("-o", outputFile) +
+    repeat("--histBinBoundaries", histogramBinBoundaries)
+
+  def summaryStats: Map[String, Any] = ConfigUtils.fileToConfigMap(outputFile)
+
+  def summaryFiles: Map[String, File] = Map.empty
+
+}

biopet-tools/src/main/resources/nl/lumc/sasc/biopet/tools/vcfstats/plotHeatmap.R

@@ -14,6 +14,7 @@ rownames(heat) <- heat[,1]
 heat<- heat[,-1]
 heat<- as.matrix(heat)
 
+textMargin <- max(sapply(rownames(heat), nchar)) + 4
 colNumber <- 50
 col <- rev(colorRampPalette(brewer.pal(11, "Spectral"))(colNumber))
 for (i in (colNumber+1):(colNumber+round((dist(range(heat)) - dist(range(heat[heat < 1]))) / dist(range(heat[heat < 1])) * colNumber))) {
@@ -22,7 +23,7 @@ for (i in (colNumber+1):(colNumber+round((dist(range(heat)) - dist(range(heat[he
 col[length(col)] <- "#00FF00"
 
 png(file = outputArg, width = 1200, height = 1200)
-heatmap.2(heat, trace = 'none', col = col, Colv=NA, Rowv=NA, dendrogram="none", margins = c(12, 12), na.color="#00FF00")
+heatmap.2(heat, trace = 'none', col = col, Colv=NA, Rowv=NA, dendrogram="none", margins = c(textMargin, textMargin), na.color="#00FF00")
 dev.off()
 
 hc <- hclust(d = dist(heat))
@@ -31,5 +32,5 @@ plot(as.dendrogram(hc), horiz=TRUE, asp=0.02)
 dev.off()
 
 png(file = outputArgClustering, width = 1200, height = 1200)
-heatmap.2(heat, trace = 'none', col = col, Colv="Rowv", dendrogram="row",margins = c(12, 12), na.color="#00FF00")
+heatmap.2(heat, trace = 'none', col = col, Colv="Rowv", dendrogram="row",margins = c(textMargin, textMargin), na.color="#00FF00")
 dev.off()

biopet-tools/src/main/scala/nl/lumc/sasc/biopet/tools/vcfstats/VcfStatsForSv.scala

+package nl.lumc.sasc.biopet.tools.vcfstats
+
+import java.io.File
+
+import htsjdk.variant.vcf.VCFFileReader
+import nl.lumc.sasc.biopet.utils.{ ConfigUtils, ToolCommand }
+
+import scala.collection.JavaConversions._
+
+object VcfStatsForSv extends ToolCommand {
+  /** Commandline arguments */
+  case class Args(inputFile: File = null, outputFile: File = null, histBinBoundaries: Array[Int] = Array()) extends AbstractArgs
+
+  /** Parsing commandline arguments */
+  class OptParser extends AbstractOptParser {
+    opt[File]('i', "inputFile") required () maxOccurs 1 valueName "<file>" action { (x, c) =>
+      c.copy(inputFile = x)
+    } validate {
+      x => if (x.exists) success else failure("Input VCF required")
+    } text "Input VCF file (required)"
+
+    opt[File]('o', "outputFile") required () maxOccurs 1 valueName "<file>" action { (x, c) =>
+      c.copy(outputFile = x)
+    } text "Output file (required)"
+
+    opt[Int]("histBinBoundaries") required () unbounded () action { (x, c) =>
+      c.copy(histBinBoundaries = c.histBinBoundaries :+ x)
+    } text "When counting the records, sv-s are divided to different size classes, this parameter should give the boundaries between these classes."
+  }
+
+  protected var cmdArgs: Args = _
+
+  def main(args: Array[String]): Unit = {
+    val argsParser = new OptParser
+    cmdArgs = argsParser.parse(args, Args()) getOrElse (throw new IllegalArgumentException)
+
+    logger.info(s"Parsing file: ${cmdArgs.inputFile}")
+
+    val stats: Map[String, Any] = getVariantCounts(cmdArgs.inputFile, cmdArgs.histBinBoundaries)
+
+    ConfigUtils.mapToYamlFile(stats, cmdArgs.outputFile)
+
+  }
+
+  /** Parses a vcf-file and counts sv-s by type and size. Sv-s are divided to different size classes, the parameter histogramBinBoundaries gives the boundaries between these classes. */
+  def getVariantCounts(vcfFile: File, histogramBinBoundaries: Array[Int]): Map[String, Any] = {
+    val delCounts, insCounts, dupCounts, invCounts = Array.fill(histogramBinBoundaries.size + 1) { 0 }
+    var traCount = 0
+
+    val reader = new VCFFileReader(vcfFile, false)
+    for (record <- reader) {
+      record.getAttributeAsString("SVTYPE", "") match {
+        case "TRA" | "CTX" | "ITX" => traCount += 1
+        case svType => {
+          val size = record.getEnd - record.getStart
+          var i = 0
+          while (i < histogramBinBoundaries.size && size > histogramBinBoundaries(i)) i += 1
+          svType match {
+            case "DEL" => delCounts(i) += 1
+            case "INS" => insCounts(i) += 1
+            case "DUP" => dupCounts(i) += 1
+            case "INV" => invCounts(i) += 1
+            case _     => logger.warn(s"Vcf file contains a record of unknown type: file-$vcfFile, type-$svType")
+          }
+        }
+      }
+    }
+    reader.close()
+
+    Map("DEL" -> delCounts, "INS" -> insCounts, "DUP" -> dupCounts, "INV" -> invCounts, "TRA" -> traCount)
+  }
+
+}

biopet-utils/src/main/resources/nl/lumc/sasc/biopet/utils/rscript/plotXY.R

@@ -13,6 +13,11 @@ parser$add_argument('--llabel', dest='llabel', type='character')
 parser$add_argument('--title', dest='title', type='character')
 parser$add_argument('--hideLegend', dest='hideLegend', type='character', default="false")
 parser$add_argument('--removeZero', dest='removeZero', type='character', default="false")
+parser$add_argument('--xLog10', dest='xLog10', type='character', default="false")
+parser$add_argument('--yLog10', dest='yLog10', type='character', default="false")
+
+parser$add_argument('--xLog10Breaks', dest='xLog10Breaks', nargs='+', type='integer')
+parser$add_argument('--xLog10Labels', dest='xLog10Labels', nargs='+', type='character')
 
 arguments <- parser$parse_args()
 
@@ -31,14 +36,29 @@ DF1 <- melt(DF, id.var="Rank")
 if (arguments$removeZero == "true") DF1 <- DF1[DF1$value > 0, ]
 if (arguments$removeZero == "true") print("Removed 0 values")
 
-ggplot(DF1, aes(x = Rank, y = value, group = variable, color = variable)) +
+plot = ggplot(DF1, aes(x = Rank, y = value, group = variable, color = variable)) +
   xlab(xlab) +
   ylab(arguments$ylabel) +
-  guides(fill=guide_legend(title=arguments$llabel)) +
+  guides(color=guide_legend(title=arguments$llabel)) +
   theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
   ggtitle(arguments$title) +
   theme_bw() +
   theme(legend.position = legendPosition) +
   geom_line()