ShivaReport.scala 6.65 KB
Newer Older
Peter van 't Hof's avatar
Peter van 't Hof committed
1
2
3
4
5
6
7
8
9
10
/**
 * 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
 *
11
 * A dual licensing mode is applied. The source code within this project is freely available for non-commercial use under an AGPL
Peter van 't Hof's avatar
Peter van 't Hof committed
12
13
14
 * license; For commercial users or users who do not want to follow the AGPL
 * license, please contact us to obtain a separate license.
 */
15
16
package nl.lumc.sasc.biopet.pipelines.shiva

Peter van 't Hof's avatar
Peter van 't Hof committed
17
import java.io.{ File, PrintWriter }
18

19
import nl.lumc.sasc.biopet.core.report._
20
import nl.lumc.sasc.biopet.pipelines.mapping.MultisampleMappingReportTrait
Peter van 't Hof's avatar
Peter van 't Hof committed
21
import nl.lumc.sasc.biopet.utils.config.Configurable
22
import nl.lumc.sasc.biopet.utils.rscript.StackedBarPlot
Peter van 't Hof's avatar
Peter van 't Hof committed
23
import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb
Peter van 't Hof's avatar
Peter van 't Hof committed
24
25
import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb.Implicts._
import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb.{ NoModule, NoSample, SampleId }
Peter van 't Hof's avatar
Peter van 't Hof committed
26

Peter van 't Hof's avatar
Peter van 't Hof committed
27
import scala.concurrent.{ Await, Future }
Peter van 't Hof's avatar
Peter van 't Hof committed
28
import scala.concurrent.duration.Duration
29
30

/**
31
32
 * With this extension the report is executed within a pipeline
 *
33
34
 * Created by pjvan_thof on 3/30/15.
 */
Peter van 't Hof's avatar
Peter van 't Hof committed
35
class ShivaReport(val parent: Configurable) extends ReportBuilderExtension {
36
  def builder = ShivaReport
37
38
}

39
/** Object for report generation for Shiva pipeline */
40
41
42
object ShivaReport extends ShivaReportTrait

/** Trait for report generation for Shiva pipeline, this can be extended */
Peter van 't Hof's avatar
Peter van 't Hof committed
43
trait ShivaReportTrait extends MultisampleMappingReportTrait {
44

Peter van 't Hof's avatar
Peter van 't Hof committed
45
  def variantcallingExecuted = summary.getSettingKeys(runId, "shiva", NoModule, keyValues = Map("multisample_variantcalling" -> List("multisample_variantcalling"))).get("multisample_variantcalling")
Peter van 't Hof's avatar
Peter van 't Hof committed
46
    .flatten match {
Peter van 't Hof's avatar
Peter van 't Hof committed
47
48
49
      case Some(true) => true
      case _          => false
    }
50

51
52
  override def frontSection = ReportSection("/nl/lumc/sasc/biopet/pipelines/shiva/shivaFront.ssp")

Peter van 't Hof's avatar
Peter van 't Hof committed
53
54
  override def pipelineName = "shiva"

55
56
57
58
  override def additionalSections = super.additionalSections ++ (if (variantcallingExecuted) List("Variantcalling" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/shiva/sampleVariants.ssp",
    Map("showPlot" -> true, "showTable" -> false)))
  else Nil)

59
  /** Root page for the shiva report */
Peter van 't Hof's avatar
Peter van 't Hof committed
60
  override def indexPage: Future[ReportPage] = Future {
61
    val regions = regionsPage
Peter van 't Hof's avatar
Peter van 't Hof committed
62
    val oldPage = Await.result(super.indexPage, Duration.Inf)
63

64
    oldPage.copy(subPages = oldPage.subPages ++ regionsPage)
65
  }
66

67
  /** Generate a page with all target coverage stats */
Peter van 't Hof's avatar
Peter van 't Hof committed
68
  def regionsPage: Option[(String, Future[ReportPage])] = {
69
70
71
    val shivaSettings = Await.result(summary.getSetting(runId, "shiva"), Duration.Inf).get
    val roi = shivaSettings.get("regions_of_interest")
    val amplicon = shivaSettings.get("amplicon_bed")
72
73
74
75

    var regionPages: Map[String, ReportPage] = Map()

    def createPage(name: String, amplicon: Boolean = false): ReportPage = {
Peter van 't Hof's avatar
Peter van 't Hof committed
76
      ReportPage(
77
        List(),
Peter van 't Hof's avatar
Peter van 't Hof committed
78
79
80
        List("Coverage" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/covstatsMultiTable.ssp")),
        Map("target" -> name)
      )
81
82
83
    }

    amplicon match {
84
      case Some(x: String) => regionPages += (x + " (Amplicon)") -> createPage(x, amplicon = true)
85
86
87
88
      case _               =>
    }

    roi match {
89
90
      case Some(x: String)  => regionPages += x -> createPage(x, amplicon = false)
      case Some(x: List[_]) => x.foreach(x => regionPages += x.toString -> createPage(x.toString, amplicon = false))
91
92
93
      case _                =>
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
94
95
    if (regionPages.nonEmpty) Some("Regions" -> Future(ReportPage(
      regionPages.map(p => p._1 -> Future(ReportPage(Nil,
96
        List(
Peter van 't Hof's avatar
Peter van 't Hof committed
97
          "Variants" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/shiva/sampleVariants.ssp", Map("showPlot" -> true)),
98
99
100
          "Coverage" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/covstatsMultiTable.ssp")
        ),
        Map("target" -> Some(p._1.stripSuffix(" (Amplicon)")))
Peter van 't Hof's avatar
Peter van 't Hof committed
101
      ))).toList.sortBy(_._1),
102
      List(),
Peter van 't Hof's avatar
Peter van 't Hof committed
103
      Map()))
Peter van 't Hof's avatar
Peter van 't Hof committed
104
105
    )
    else None
106
107
  }

108
  /** Single sample page */
Peter van 't Hof's avatar
Peter van 't Hof committed
109
  override def samplePage(sampleId: Int, args: Map[String, Any]): Future[ReportPage] = Future {
110
    val variantcallingSection = if (variantcallingExecuted) List("Variantcalling" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/shiva/sampleVariants.ssp")) else Nil
Peter van 't Hof's avatar
Peter van 't Hof committed
111
    val oldPage: ReportPage = super.samplePage(sampleId, args)
112
    oldPage.copy(sections = variantcallingSection ++ oldPage.sections)
113
114
  }

115
  /** Name of the report */
116
  def reportName = "Shiva Report"
117

118
119
  /**
   * Generate a stackbar plot for found variants
Peter van 't Hof's avatar
Peter van 't Hof committed
120
   *
121
122
123
124
125
   * @param outputDir OutputDir for the tsv and png file
   * @param prefix Prefix of the tsv and png file
   * @param summary Summary class
   * @param sampleId Default it selects all sampples, when sample is giving it limits to selected sample
   */
126
  def variantSummaryPlot(outputDir: File,
127
                         prefix: String,
Peter van 't Hof's avatar
Peter van 't Hof committed
128
129
                         summary: SummaryDb,
                         sampleId: Option[Int] = None,
130
131
                         caller: String = "final",
                         target: Option[String] = None): Unit = {
132
133
134
    val tsvFile = new File(outputDir, prefix + ".tsv")
    val pngFile = new File(outputDir, prefix + ".png")
    val tsvWriter = new PrintWriter(tsvFile)
Peter van 't Hof's avatar
Peter van 't Hof committed
135
136
137
138
139
140
141
    tsvWriter.print("Sample")
    val field = List("HomVar", "Het", "HomRef", "NoCall")
    tsvWriter.println(s"\t${field.mkString("\t")}")

    val samples = Await.result(summary.getSamples(runId = runId, sampleId = sampleId), Duration.Inf)
    val statsPaths = {
      (for (sample <- Await.result(summary.getSamples(runId = runId), Duration.Inf)) yield {
Peter van 't Hof's avatar
Peter van 't Hof committed
142
        field.map(f => s"${sample.name};$f" -> List("total", "genotype", "general", sample.name, f)).toMap
Peter van 't Hof's avatar
Peter van 't Hof committed
143
144
145
146
147
148
      }).fold(Map())(_ ++ _)
    }

    val moduleName = target match {
      case Some(t) => s"multisample-vcfstats-$caller-$t"
      case _       => s"multisample-vcfstats-$caller"
149
150
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
151
    val results = summary.getStatKeys(runId, "shivavariantcalling", moduleName, sampleId.map(SampleId).getOrElse(NoSample), keyValues = statsPaths)
Peter van 't Hof's avatar
Peter van 't Hof committed
152
153

    for (sample <- samples if sampleId.isEmpty || sample.id == sampleId.get) {
Peter van 't Hof's avatar
Peter van 't Hof committed
154
      tsvWriter.println(sample.name + "\t" + field.map(f => results.get(s"${sample.name};$f").getOrElse(Some("0")).get).mkString("\t"))
155
156
157
158
159
160
161
162
    }

    tsvWriter.close()

    val plot = new StackedBarPlot(null)
    plot.input = tsvFile
    plot.output = pngFile
    plot.ylabel = Some("VCF records")
Peter van 't Hof's avatar
Peter van 't Hof committed
163
    plot.width = Some(200 + (samples.count(s => sampleId.getOrElse(s) == s) * 10))
164
165
    plot.runLocal()
  }
166
}