MuTect2.scala 3.7 KB
Newer Older
akaljuvee's avatar
-  
akaljuvee committed
1
2
package nl.lumc.sasc.biopet.pipelines.shiva.variantcallers.somatic

Peter van 't Hof's avatar
Peter van 't Hof committed
3
import nl.lumc.sasc.biopet.extensions.bcftools.BcftoolsReheader
4
import nl.lumc.sasc.biopet.extensions.gatk.CombineVariants
Peter van 't Hof's avatar
Peter van 't Hof committed
5
import nl.lumc.sasc.biopet.extensions._
6
import nl.lumc.sasc.biopet.extensions.gatk.gather.BqsrGather
pjvan_thof's avatar
pjvan_thof committed
7
import nl.lumc.sasc.biopet.utils.{IoUtils, Logging}
Peter van 't Hof's avatar
Peter van 't Hof committed
8
9
import nl.lumc.sasc.biopet.utils.config.Configurable
import org.broadinstitute.gatk.queue.extensions.gatk.TaggedFile
akaljuvee's avatar
-  
akaljuvee committed
10
11
12
13
14
15
16

class MuTect2(val parent: Configurable) extends SomaticVariantCaller {

  def name = "mutect2"

  override val mergeVcfResults: Boolean = false

Peter van 't Hof's avatar
Peter van 't Hof committed
17
  def defaultPrio: Int = -1
akaljuvee's avatar
-  
akaljuvee committed
18

Peter van 't Hof's avatar
Peter van 't Hof committed
19
  lazy val runConEst: Boolean = config("run_contest", default = false)
akaljuvee's avatar
-  
akaljuvee committed
20
21
22

  def biopetScript(): Unit = {

pjvan_thof's avatar
Typo    
pjvan_thof committed
23
    if (tnPairs.isEmpty) Logging.addError("No tumor-normal found in config")
Peter van 't Hof's avatar
Peter van 't Hof committed
24
    val outputFiles = for (pair <- tnPairs) yield {
akaljuvee's avatar
-  
akaljuvee committed
25

Peter van 't Hof's avatar
Peter van 't Hof committed
26
27
28
29
30
31
32
33
34
      val bqsrFile =
        if (inputBqsrFiles.contains(pair.tumorSample) &&
            inputBqsrFiles.contains(pair.normalSample)) {
          val gather = new BqsrGather()
          gather.inputBqsrFiles =
            List(inputBqsrFiles(pair.tumorSample), inputBqsrFiles(pair.normalSample))
          gather.outputBqsrFile =
            new File(outputDir, s"${pair.tumorSample}-${pair.normalSample}.bqsr")
          add(gather)
akaljuvee's avatar
-  
akaljuvee committed
35

Peter van 't Hof's avatar
Peter van 't Hof committed
36
37
38
39
40
41
42
          Some(gather.outputBqsrFile)
        } else None

      val outputFile =
        new File(outputDir, s"${pair.tumorSample}-${pair.normalSample}.$name.vcf.gz")

      val muTect2 = new gatk.MuTect2(this)
pjvan_thof's avatar
pjvan_thof committed
43
44
      inputBams.get(pair.tumorSample).foreach(muTect2.input_file :+= TaggedFile(_, "tumor"))
      inputBams.get(pair.normalSample).foreach(muTect2.input_file :+= TaggedFile(_, "normal"))
Peter van 't Hof's avatar
Peter van 't Hof committed
45
      muTect2.BQSR = bqsrFile
akaljuvee's avatar
-  
akaljuvee committed
46

Peter van 't Hof's avatar
Peter van 't Hof committed
47
48
      if (runConEst) {
        val namePrefix = outputFile.getAbsolutePath.stripSuffix(".vcf.gz")
pjvan_thof's avatar
pjvan_thof committed
49
        val contEst = new gatk.ContEst(this)
Peter van 't Hof's avatar
Peter van 't Hof committed
50
51
        inputBams.get(pair.tumorSample).foreach(contEst.input_file :+= TaggedFile(_, "eval"))
        inputBams.get(pair.normalSample).foreach(contEst.input_file :+= TaggedFile(_, "genotype"))
pjvan_thof's avatar
pjvan_thof committed
52
        contEst.output = new File(s"$namePrefix.contamination.txt")
Peter van 't Hof's avatar
Peter van 't Hof committed
53
        contEst.BQSR = bqsrFile
Peter van 't Hof's avatar
Peter van 't Hof committed
54
        add(contEst)
akaljuvee's avatar
-  
akaljuvee committed
55

Peter van 't Hof's avatar
Peter van 't Hof committed
56
57
        val contaminationPerSample: File = new File(s"$namePrefix.contamination.short.txt")
        val awk: Awk = Awk(this, "BEGIN{OFS=\"\\t\"}{if($1 != \"name\") print $1,$4;}")
pjvan_thof's avatar
pjvan_thof committed
58
        awk.input = contEst.output
Peter van 't Hof's avatar
Peter van 't Hof committed
59
60
61
62
        add(awk > contaminationPerSample)

        muTect2.contaminationFile = Some(contaminationPerSample)
      }
akaljuvee's avatar
-  
akaljuvee committed
63

Peter van 't Hof's avatar
Peter van 't Hof committed
64
65
      if (inputBqsrFiles.contains(pair.tumorSample) && inputBqsrFiles
            .contains(pair.normalSample)) {
Peter van 't Hof's avatar
Peter van 't Hof committed
66
67
68
69
70
71
72
73
74
        val gather = new BqsrGather()
        gather.inputBqsrFiles =
          List(inputBqsrFiles(pair.tumorSample), inputBqsrFiles(pair.normalSample))
        gather.outputBqsrFile = new File(swapExt(outputFile, "vcf.gz", "bqsr.merge"))
        add(gather)

        muTect2.BQSR = Some(gather.outputBqsrFile)
      }

pjvan_thof's avatar
pjvan_thof committed
75
      outputDir.mkdirs()
Peter van 't Hof's avatar
Peter van 't Hof committed
76
77
78
79
80
81
      val renameFile = new File(outputDir, s".rename.${pair.tumorSample}-${pair.normalSample}.txt")
      IoUtils.writeLinesToFile(renameFile,
                               List(
                                 s"TUMOR ${pair.tumorSample}",
                                 s"NORMAL ${pair.normalSample}"
                               ))
Peter van 't Hof's avatar
Peter van 't Hof committed
82

Peter van 't Hof's avatar
Peter van 't Hof committed
83
84
85
      val pipe = muTect2 | BcftoolsReheader(this, renameFile) | new Bgzip(this) > outputFile
      pipe.threadsCorrection = -2
      add(pipe)
Peter van 't Hof's avatar
Peter van 't Hof committed
86
87
88
89
      add(Tabix(this, outputFile))

      outputFile
    }
Peter van 't Hof's avatar
Peter van 't Hof committed
90
91
92
93
94
95
96

    if (outputFiles.size > 1) {
      add(CombineVariants(this, outputFiles, outputFile))
    } else if (outputFiles.nonEmpty) {
      add(Ln(this, outputFiles.head, outputFile))
      add(Ln(this, outputFiles.head + ".tbi", outputFile + ".tbi"))
    }
akaljuvee's avatar
-  
akaljuvee committed
97
  }
pjvan_thof's avatar
pjvan_thof committed
98
}