Adding gender calling haplotypecaller

shiva/src/main/scala/nl/lumc/sasc/biopet/pipelines/shiva/ShivaVariantcalling.scala

@@ -20,10 +20,10 @@ import nl.lumc.sasc.biopet.core.summary.SummaryQScript
 import nl.lumc.sasc.biopet.extensions.Tabix
 import nl.lumc.sasc.biopet.extensions.gatk.{ CombineVariants, GenotypeConcordance }
 import nl.lumc.sasc.biopet.extensions.tools.VcfStats
-import nl.lumc.sasc.biopet.extensions.vt.{VtDecompose, VtNormalize}
+import nl.lumc.sasc.biopet.extensions.vt.{ VtDecompose, VtNormalize }
 import nl.lumc.sasc.biopet.pipelines.bammetrics.TargetRegions
 import nl.lumc.sasc.biopet.pipelines.shiva.variantcallers.{ VarscanCnsSingleSample, _ }
-import nl.lumc.sasc.biopet.utils.{BamUtils, Logging}
+import nl.lumc.sasc.biopet.utils.{ BamUtils, Logging }
 import nl.lumc.sasc.biopet.utils.config.Configurable
 import org.broadinstitute.gatk.queue.QScript
 import org.broadinstitute.gatk.queue.extensions.gatk.TaggedFile
@@ -56,12 +56,13 @@ class ShivaVariantcalling(val parent: Configurable) extends QScript
     if (inputBamsArg.nonEmpty) inputBams = BamUtils.sampleBamMap(inputBamsArg)
     if (genders == null) genders = {
       val samples: Map[String, Any] = config("genders")
-      samples.map { case (sampleName, gender) =>
-        sampleName -> (gender.toString.toLowerCase match {
-          case "male"   => Gender.Male
-          case "female" => Gender.Female
-          case _ => Gender.Unknown
-        })
+      samples.map {
+        case (sampleName, gender) =>
+          sampleName -> (gender.toString.toLowerCase match {
+            case "male"   => Gender.Male
+            case "female" => Gender.Female
+            case _        => Gender.Unknown
+          })
       }
     }
   }

shiva/src/main/scala/nl/lumc/sasc/biopet/pipelines/shiva/variantcallers/HaplotypeCallerGvcf.scala

@@ -19,7 +19,9 @@
  */
 package nl.lumc.sasc.biopet.pipelines.shiva.variantcallers
 
+import nl.lumc.sasc.biopet.core.MultiSampleQScript.Gender
 import nl.lumc.sasc.biopet.extensions.gatk
+import nl.lumc.sasc.biopet.extensions.gatk.CombineGVCFs
 import nl.lumc.sasc.biopet.utils.config.Configurable
 
 /** Gvcf mode for haplotypecaller */
@@ -34,14 +36,51 @@ class HaplotypeCallerGvcf(val parent: Configurable) extends Variantcaller {
 
   def getGvcfs = gVcfFiles
 
+  val genderAwareCalling: Boolean = config("gender_aware_calling", default = false)
+  val haploidRegions: Option[File] = config("hap̦loid_regions")
+  val haploidRegionsMale: Option[File] = config("haploid_regions")
+  val haploidRegionsFemale: Option[File] = config("haploid_regions")
+
   override def fixedValues = Map("haplotypecaller" -> Map("emitRefConfidence" -> "GVCF"))
 
   def biopetScript() {
     gVcfFiles = for ((sample, inputBam) <- inputBams) yield {
-      val hc = gatk.HaplotypeCaller(this, List(inputBam), new File(outputDir, sample + ".gvcf.vcf.gz"))
-      hc.BQSR = inputBqsrFiles.get(sample)
-      add(hc)
-      sample -> hc.out
+      if (genderAwareCalling) {
+        val finalFile = new File(outputDir, sample + ".gvcf.vcf.gz")
+        val haploidBedFiles: List[File] = genders.getOrElse(sample, Gender.Unknown) match {
+          case Gender.Female => haploidRegions.toList ::: haploidRegionsFemale.toList ::: Nil
+          case Gender.Male   => haploidRegions.toList ::: haploidRegionsMale.toList ::: Nil
+          case _             => haploidRegions.toList
+        }
+
+        val haploidGvcf = if (haploidBedFiles.nonEmpty) {
+          val hc = gatk.HaplotypeCaller(this, List(inputBam), new File(outputDir, sample + ".haploid.gvcf.vcf.gz"))
+          hc.BQSR = inputBqsrFiles.get(sample)
+          hc.intervals = haploidBedFiles
+          add(hc)
+          Some(hc.out)
+        } else None
+
+        val hcDiploid = gatk.HaplotypeCaller(this, List(inputBam), new File(outputDir, sample + ".diploid.gvcf.vcf.gz"))
+        hcDiploid.BQSR = inputBqsrFiles.get(sample)
+        hcDiploid.excludeIntervals = haploidBedFiles
+        add(hcDiploid)
+
+        haploidGvcf match {
+          case Some(file) =>
+            val combine = new CombineGVCFs(this)
+            combine.variant = Seq(hcDiploid.out, file)
+            combine.out = new File(outputDir, sample + ".gvcf.vcf.gz")
+            add(combine)
+            sample -> combine.out
+          case _ => sample -> hcDiploid.out
+        }
+      } else {
+        val hc = gatk.HaplotypeCaller(this, List(inputBam), new File(outputDir, sample + ".gvcf.vcf.gz"))
+        hc.BQSR = inputBqsrFiles.get(sample)
+        add(hc)
+        sample -> hc.out
+      }
     }
 
     val genotypeGVCFs = gatk.GenotypeGVCFs(this, gVcfFiles.values.toList, outputFile)