Added IndelRealigner and BaseRecalibration

protected/biopet-gatk-pipelines/src/main/scala/nl/lumc/sasc/biopet/pipelines/gatk/ShivaGatk.scala

@@ -2,7 +2,8 @@ package nl.lumc.sasc.biopet.pipelines.gatk
 
 import nl.lumc.sasc.biopet.core.PipelineCommand
 import nl.lumc.sasc.biopet.core.config.Configurable
-import nl.lumc.sasc.biopet.pipelines.shiva.{ShivaTrait, Shiva}
+import nl.lumc.sasc.biopet.extensions.gatk._
+import nl.lumc.sasc.biopet.pipelines.shiva.ShivaTrait
 import org.broadinstitute.gatk.queue.QScript
 
 /**
@@ -15,21 +16,74 @@ class ShivaGatk(val root: Configurable) extends QScript with ShivaTrait {
   class Sample(sampleId: String) extends super.Sample(sampleId) {
     override def makeLibrary(id: String) = new this.Library(id)
     class Library(libId: String) extends super.Library(libId) {
-      override def preProcess(input:File): Option[File] = {
-        //TODO: add preproces
-        None
+      override def preProcess(input: File): Option[File] = {
+        val useIndelRealigner: Boolean = config("use_indel_realign", default = true)
+        val useBaseRecalibration: Boolean = config("use_base_recalibration", default = true)
+
+        if (!useIndelRealigner && !useBaseRecalibration) None
+        else {
+          val indelRealignFile = useIndelRealigner match {
+            case true  => addIndelRealign(input, libDir, useBaseRecalibration || libraries.size > 1)
+            case false => input
+          }
+
+          useBaseRecalibration match {
+            case true  => Some(addBaseRecalibrator(indelRealignFile, libDir, libraries.size > 1))
+            case false => Some(indelRealignFile)
+          }
+        }
       }
     }
 
     override def doublePreProcess(input: List[File], isIntermediate: Boolean = false): Option[File] = {
       if (input.size <= 1) super.doublePreProcess(input)
-      else {
-        super.doublePreProcess(input)
-        //TODO: Add indel realignment
+      else super.doublePreProcess(input, true).collect {
+        case file => {
+          config("use_indel_realign", default = true).asBoolean match {
+            case true  => addIndelRealign(file, sampleDir, false)
+            case false => file
+          }
+        }
       }
     }
   }
 
+  def addIndelRealign(inputBam: File, dir: File, isIntermediate: Boolean): File = {
+    val realignerTargetCreator = RealignerTargetCreator(this, inputBam, dir)
+    realignerTargetCreator.isIntermediate = true
+    add(realignerTargetCreator)
+
+    val indelRealigner = IndelRealigner(this, inputBam, realignerTargetCreator.out, dir)
+    indelRealigner.isIntermediate = isIntermediate
+    add(indelRealigner)
+
+    return indelRealigner.o
+  }
+
+  def addBaseRecalibrator(inputBam: File, dir: File, isIntermediate: Boolean): File = {
+    val baseRecalibrator = BaseRecalibrator(this, inputBam, swapExt(dir, inputBam, ".bam", ".baserecal"))
+
+    if (baseRecalibrator.knownSites.isEmpty) {
+      logger.warn("No Known site found, skipping base recalibration, file: " + inputBam)
+      return inputBam
+    }
+    add(baseRecalibrator)
+
+    if (config("use_analyze_covariates", default = false).asBoolean) {
+      val baseRecalibratorAfter = BaseRecalibrator(this, inputBam, swapExt(dir, inputBam, ".bam", ".baserecal.after"))
+      baseRecalibratorAfter.BQSR = baseRecalibrator.o
+      add(baseRecalibratorAfter)
+
+      add(AnalyzeCovariates(this, baseRecalibrator.o, baseRecalibratorAfter.o, swapExt(dir, inputBam, ".bam", ".baserecal.pdf")))
+    }
+
+    val printReads = PrintReads(this, inputBam, swapExt(dir, inputBam, ".bam", ".baserecal.bam"))
+    printReads.BQSR = baseRecalibrator.o
+    printReads.isIntermediate = isIntermediate
+    add(printReads)
+
+    return printReads.o
+  }
 }
 
 object ShivaGatk extends PipelineCommand
\ No newline at end of file

public/shiva/src/main/scala/nl/lumc/sasc/biopet/pipelines/shiva/ShivaTrait.scala

@@ -4,9 +4,9 @@ import java.io.File
 
 import htsjdk.samtools.SamReaderFactory
 import nl.lumc.sasc.biopet.core.summary.SummaryQScript
-import nl.lumc.sasc.biopet.core.{MultiSampleQScript, PipelineCommand}
+import nl.lumc.sasc.biopet.core.MultiSampleQScript
 import nl.lumc.sasc.biopet.extensions.Ln
-import nl.lumc.sasc.biopet.extensions.picard.{AddOrReplaceReadGroups, SamToFastq, MarkDuplicates}
+import nl.lumc.sasc.biopet.extensions.picard.{ AddOrReplaceReadGroups, SamToFastq, MarkDuplicates }
 import nl.lumc.sasc.biopet.pipelines.mapping.Mapping
 import scala.collection.JavaConversions._
 
@@ -62,7 +62,7 @@ trait ShivaTrait extends MultiSampleQScript with SummaryQScript {
             mapping.input_R1 = config("R1")
             mapping.input_R2 = config("R2")
           })
-          case (false, true) =>  config("bam_to_fastq", default = false).asBoolean match {
+          case (false, true) => config("bam_to_fastq", default = false).asBoolean match {
             case true => {
               val samToFastq = SamToFastq(qscript, config("bam"),
                 new File(libDir, sampleId + "-" + libId + ".R1.fastq"),
@@ -78,7 +78,7 @@ trait ShivaTrait extends MultiSampleQScript with SummaryQScript {
               val inputSam = SamReaderFactory.makeDefault.open(config("bam"))
               val readGroups = inputSam.getFileHeader.getReadGroups
 
-              val readGroupOke = readGroups.forall( readGroup => {
+              val readGroupOke = readGroups.forall(readGroup => {
                 if (readGroup.getSample != sampleId) logger.warn("Sample ID readgroup in bam file is not the same")
                 if (readGroup.getLibrary != libId) logger.warn("Library ID readgroup in bam file is not the same")
                 readGroup.getSample == sampleId && readGroup.getLibrary == libId