Shiva.scala 9.87 KB
Newer Older
bow's avatar
bow committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
/**
 * 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 that are
 * not part of GATK Queue 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.
 */
Peter van 't Hof's avatar
Peter van 't Hof committed
16
17
package nl.lumc.sasc.biopet.pipelines.shiva

18
import nl.lumc.sasc.biopet.core.{ PipelineCommand, Reference }
Peter van 't Hof's avatar
Peter van 't Hof committed
19
import nl.lumc.sasc.biopet.core.report.ReportBuilderExtension
Peter van 't Hof's avatar
Peter van 't Hof committed
20
import nl.lumc.sasc.biopet.extensions.gatk._
21
import nl.lumc.sasc.biopet.pipelines.bammetrics.TargetRegions
Peter van 't Hof's avatar
Peter van 't Hof committed
22
import nl.lumc.sasc.biopet.pipelines.kopisu.Kopisu
Peter van 't Hof's avatar
Peter van 't Hof committed
23
import nl.lumc.sasc.biopet.pipelines.mapping.MultisampleMappingTrait
Peter van 't Hof's avatar
Peter van 't Hof committed
24
import nl.lumc.sasc.biopet.pipelines.toucan.Toucan
Peter van 't Hof's avatar
Peter van 't Hof committed
25
import nl.lumc.sasc.biopet.utils.config.Configurable
26
import org.broadinstitute.gatk.queue.QScript
Peter van 't Hof's avatar
Peter van 't Hof committed
27

Peter van 't Hof's avatar
Peter van 't Hof committed
28
/**
Peter van 't Hof's avatar
Peter van 't Hof committed
29
30
 * This is a trait for the Shiva pipeline
 *
Peter van 't Hof's avatar
Peter van 't Hof committed
31
32
 * Created by pjvan_thof on 2/26/15.
 */
33
34
class Shiva(val root: Configurable) extends QScript with MultisampleMappingTrait with Reference with TargetRegions { qscript =>

Peter van 't Hof's avatar
Peter van 't Hof committed
35
  def this() = this(null)
Peter van 't Hof's avatar
Peter van 't Hof committed
36

Peter van 't Hof's avatar
Peter van 't Hof committed
37
  override def reportClass: Option[ReportBuilderExtension] = {
38
39
40
41
42
43
    val shiva = new ShivaReport(this)
    shiva.outputDir = new File(outputDir, "report")
    shiva.summaryFile = summaryFile
    Some(shiva)
  }

44
45
46
47
48
49
  override def defaults = Map(
    "haplotypecaller" -> Map("stand_call_conf" -> 30, "stand_emit_conf" -> 0),
    "genotypegvcfs" -> Map("stand_call_conf" -> 30, "stand_emit_conf" -> 0),
    "unifiedgenotyper" -> Map("stand_call_conf" -> 30, "stand_emit_conf" -> 0)
  )

Sander Bollen's avatar
Sander Bollen committed
50
  /** Method to make the variantcalling namespace of shiva */
51
  def makeVariantcalling(multisample: Boolean = false): ShivaVariantcalling with QScript = {
Peter van 't Hof's avatar
Peter van 't Hof committed
52
    if (multisample) new ShivaVariantcalling(qscript) {
Peter van 't Hof's avatar
Peter van 't Hof committed
53
      override def namePrefix = "multisample"
Sander Bollen's avatar
Sander Bollen committed
54
      override def configNamespace: String = "shivavariantcalling"
Peter van 't Hof's avatar
Peter van 't Hof committed
55
56
57
      override def configPath: List[String] = super.configPath ::: "multisample" :: Nil
    }
    else new ShivaVariantcalling(qscript) {
Sander Bollen's avatar
Sander Bollen committed
58
      override def configNamespace = "shivavariantcalling"
Peter van 't Hof's avatar
Peter van 't Hof committed
59
60
61
    }
  }

Peter van 't Hof's avatar
Peter van 't Hof committed
62
  /** Method to make a sample */
63
  override def makeSample(id: String) = new this.Sample(id)
64

Peter van 't Hof's avatar
Peter van 't Hof committed
65
  /** Class that will generate jobs for a sample */
66
  class Sample(sampleId: String) extends super.Sample(sampleId) {
Peter van 't Hof's avatar
Peter van 't Hof committed
67
    /** Method to make a library */
68
    override def makeLibrary(id: String) = new this.Library(id)
Peter van 't Hof's avatar
Peter van 't Hof committed
69

70
    /** Sample specific settings */
71
72
    override def summarySettings = super.summarySettings ++
      Map("single_sample_variantcalling" -> variantcalling.isDefined, "use_indel_realigner" -> useIndelRealigner)
73

Peter van 't Hof's avatar
Peter van 't Hof committed
74
    /** Class to generate jobs for a library */
75
    class Library(libId: String) extends super.Library(libId) {
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

      lazy val useIndelRealigner: Boolean = config("use_indel_realigner", default = true)
      lazy val useBaseRecalibration: Boolean = {
        val c: Boolean = config("use_base_recalibration", default = true)
        val br = new BaseRecalibrator(qscript)
        if (c && br.knownSites.isEmpty)
          logger.warn("No Known site found, skipping base recalibration, file: " + inputBam)
        c && br.knownSites.nonEmpty
      }

      override def preProcessBam = if (useIndelRealigner && useBaseRecalibration)
        bamFile.map(swapExt(libDir, _, ".bam", ".realign.baserecal.bam"))
      else if (useIndelRealigner) bamFile.map(swapExt(libDir, _, ".bam", ".realign.bam"))
      else if (useBaseRecalibration) bamFile.map(swapExt(libDir, _, ".bam", ".baserecal.bam"))
      else bamFile

92
      /** Library specific settings */
93
94
95
96
      override def summarySettings = Map(
        "library_variantcalling" -> variantcalling.isDefined,
        "use_indel_realigner" -> useIndelRealigner,
        "use_base_recalibration" -> useBaseRecalibration)
97

98
      lazy val variantcalling = if (config("library_variantcalling", default = false).asBoolean &&
Peter van 't Hof's avatar
Peter van 't Hof committed
99
100
101
102
        (bamFile.isDefined || preProcessBam.isDefined)) {
        Some(makeVariantcalling(multisample = false))
      } else None

Peter van 't Hof's avatar
Peter van 't Hof committed
103
      /** This will add jobs for this library */
Peter van 't Hof's avatar
Peter van 't Hof committed
104
      override def addJobs() = {
105
        super.addJobs()
Peter van 't Hof's avatar
Peter van 't Hof committed
106

107
108
109
110
111
112
113
114
115
        if (useIndelRealigner && useBaseRecalibration) {
          val file = addIndelRealign(bamFile.get, libDir, isIntermediate = true)
          addBaseRecalibrator(file, libDir, libraries.size > 1)
        } else if (useIndelRealigner) {
          addIndelRealign(bamFile.get, libDir, libraries.size > 1)
        } else if (useBaseRecalibration) {
          addBaseRecalibrator(bamFile.get, libDir, libraries.size > 1)
        }

Peter van 't Hof's avatar
Peter van 't Hof committed
116
        variantcalling.foreach(vc => {
Peter van 't Hof's avatar
Peter van 't Hof committed
117
118
          vc.sampleId = Some(sampleId)
          vc.libId = Some(libId)
Peter van 't Hof's avatar
Peter van 't Hof committed
119
          vc.outputDir = new File(libDir, "variantcalling")
120
121
          if (preProcessBam.isDefined) vc.inputBams = Map(sampleId -> preProcessBam.get)
          else vc.inputBams = Map(sampleId -> bamFile.get)
Peter van 't Hof's avatar
Peter van 't Hof committed
122
          add(vc)
Peter van 't Hof's avatar
Peter van 't Hof committed
123
        })
Peter van 't Hof's avatar
Peter van 't Hof committed
124
125
126
      }
    }

Peter van 't Hof's avatar
fix bug    
Peter van 't Hof committed
127
    lazy val variantcalling = if (config("single_sample_variantcalling", default = false).asBoolean) {
Peter van 't Hof's avatar
Peter van 't Hof committed
128
      Some(makeVariantcalling(multisample = false))
Peter van 't Hof's avatar
Peter van 't Hof committed
129
130
    } else None

131
132
133
134
135
136
137
138
    override def keepMergedFiles: Boolean = config("keep_merged_files", default = !useIndelRealigner)

    lazy val useIndelRealigner: Boolean = config("use_indel_realigner", default = true)

    override def preProcessBam = if (useIndelRealigner && libraries.values.flatMap(_.preProcessBam).size > 1) {
      bamFile.map(swapExt(sampleDir, _, ".bam", ".realign.bam"))
    } else bamFile

Peter van 't Hof's avatar
Peter van 't Hof committed
139
    /** This will add sample jobs */
140
141
    override def addJobs(): Unit = {
      super.addJobs()
Peter van 't Hof's avatar
Peter van 't Hof committed
142

143
144
145
146
      if (useIndelRealigner && libraries.values.flatMap(_.preProcessBam).size > 1) {
        addIndelRealign(bamFile.get, sampleDir, false)
      }

147
      preProcessBam.foreach { bam =>
Peter van 't Hof's avatar
Peter van 't Hof committed
148
        variantcalling.foreach(vc => {
Peter van 't Hof's avatar
Peter van 't Hof committed
149
150
          vc.sampleId = Some(sampleId)
          vc.outputDir = new File(sampleDir, "variantcalling")
151
          vc.inputBams = Map(sampleId -> bam)
Peter van 't Hof's avatar
Peter van 't Hof committed
152
          add(vc)
Peter van 't Hof's avatar
Peter van 't Hof committed
153
        })
Peter van 't Hof's avatar
Peter van 't Hof committed
154
      }
Peter van 't Hof's avatar
Peter van 't Hof committed
155
156
157
    }
  }

158
  lazy val multisampleVariantCalling = if (config("multisample_variantcalling", default = true).asBoolean) {
Peter van 't Hof's avatar
Peter van 't Hof committed
159
160
161
    Some(makeVariantcalling(multisample = true))
  } else None

Peter van 't Hof's avatar
Peter van 't Hof committed
162
  lazy val svCalling = if (config("sv_calling", default = false).asBoolean) {
Peter van 't Hof's avatar
Peter van 't Hof committed
163
    Some(new ShivaSvCalling(this))
Peter van 't Hof's avatar
Peter van 't Hof committed
164
165
  } else None

Peter van 't Hof's avatar
Peter van 't Hof committed
166
  lazy val cnvCalling = if (config("cnv_calling", default = false).asBoolean) {
Peter van 't Hof's avatar
Peter van 't Hof committed
167
    Some(new Kopisu(this))
Peter van 't Hof's avatar
Peter van 't Hof committed
168
169
  } else None

170
171
172
173
174
  lazy val annotation = if (multisampleVariantCalling.isDefined &&
    config("annotation", default = false).asBoolean) {
    Some(new Toucan(this))
  } else None

Peter van 't Hof's avatar
Peter van 't Hof committed
175
  /** This will add the mutisample variantcalling */
Peter van 't Hof's avatar
Peter van 't Hof committed
176
  override def addMultiSampleJobs() = {
177
178
    super.addMultiSampleJobs()

179
    multisampleVariantCalling.foreach(vc => {
Peter van 't Hof's avatar
Peter van 't Hof committed
180
      vc.outputDir = new File(outputDir, "variantcalling")
181
      vc.inputBams = samples.flatMap { case (sampleId, sample) => sample.preProcessBam.map(sampleId -> _) }
Peter van 't Hof's avatar
Peter van 't Hof committed
182
      add(vc)
Peter van 't Hof's avatar
Peter van 't Hof committed
183

184
      annotation.foreach { toucan =>
Peter van 't Hof's avatar
Peter van 't Hof committed
185
        toucan.outputDir = new File(outputDir, "annotation")
186
        toucan.inputVcf = vc.finalFile
Peter van 't Hof's avatar
Peter van 't Hof committed
187
        add(toucan)
Peter van 't Hof's avatar
Peter van 't Hof committed
188
      }
Peter van 't Hof's avatar
Peter van 't Hof committed
189
    })
Peter van 't Hof's avatar
Peter van 't Hof committed
190

Peter van 't Hof's avatar
Peter van 't Hof committed
191
    svCalling.foreach { sv =>
Peter van 't Hof's avatar
Peter van 't Hof committed
192
      sv.outputDir = new File(outputDir, "sv_calling")
193
      sv.inputBams = samples.flatMap { case (sampleId, sample) => sample.preProcessBam.map(sampleId -> _) }
Peter van 't Hof's avatar
Peter van 't Hof committed
194
      add(sv)
Peter van 't Hof's avatar
Peter van 't Hof committed
195
196
197
198
199
200
201
    }

    cnvCalling.foreach { cnv =>
      cnv.outputDir = new File(outputDir, "cnv_calling")
      cnv.inputBams = samples.flatMap { case (sampleId, sample) => sample.preProcessBam.map(sampleId -> _) }
      add(cnv)
    }
Peter van 't Hof's avatar
Peter van 't Hof committed
202
203
  }

Peter van 't Hof's avatar
Peter van 't Hof committed
204
  /** Location of summary file */
Peter van 't Hof's avatar
Peter van 't Hof committed
205
  def summaryFile = new File(outputDir, "Shiva.summary.json")
Peter van 't Hof's avatar
Peter van 't Hof committed
206

Peter van 't Hof's avatar
Peter van 't Hof committed
207
  /** Settings of pipeline for summary */
208
  override def summarySettings = super.summarySettings ++ Map(
209
210
211
    "annotation" -> annotation.isDefined,
    "multisample_variantcalling" -> multisampleVariantCalling.isDefined,
    "sv_calling" -> svCalling.isDefined,
212
    "cnv_calling" -> cnvCalling.isDefined,
213
214
215
    "regions_of_interest" -> roiBedFiles.map(_.getName.stripSuffix(".bed")),
    "amplicon_bed" -> ampliconBedFile.map(_.getName.stripSuffix(".bed"))
  )
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251

  /** Adds indel realignment jobs */
  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)

    indelRealigner.out
  }

  /** Adds base recalibration jobs */
  def addBaseRecalibrator(inputBam: File, dir: File, isIntermediate: Boolean): File = {
    val baseRecalibrator = BaseRecalibrator(this, inputBam, swapExt(dir, inputBam, ".bam", ".baserecal"))

    if (baseRecalibrator.knownSites.isEmpty) return inputBam
    add(baseRecalibrator)

    if (config("use_analyze_covariates", default = true).asBoolean) {
      val baseRecalibratorAfter = BaseRecalibrator(this, inputBam, swapExt(dir, inputBam, ".bam", ".baserecal.after"))
      baseRecalibratorAfter.BQSR = Some(baseRecalibrator.out)
      add(baseRecalibratorAfter)

      add(AnalyzeCovariates(this, baseRecalibrator.out, baseRecalibratorAfter.out, swapExt(dir, inputBam, ".bam", ".baserecal.pdf")))
    }

    val printReads = PrintReads(this, inputBam, swapExt(dir, inputBam, ".bam", ".baserecal.bam"))
    printReads.BQSR = Some(baseRecalibrator.out)
    printReads.isIntermediate = isIntermediate
    add(printReads)

    printReads.out
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
252
}
Peter van 't Hof's avatar
Peter van 't Hof committed
253

Peter van 't Hof's avatar
Peter van 't Hof committed
254
/** This object give a default main method to the pipelines */
Peter van 't Hof's avatar
Peter van 't Hof committed
255
object Shiva extends PipelineCommand