ShivaReport.scala 7.46 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.{ Summary, SummaryValue }
24
25

/**
26
27
 * With this extension the report is executed within a pipeline
 *
28
29
 * Created by pjvan_thof on 3/30/15.
 */
30
class ShivaReport(val root: Configurable) extends ReportBuilderExtension {
31
  def builder = ShivaReport
32
33
}

34
/** Object for report generation for Shiva pipeline */
35
36
37
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
38
trait ShivaReportTrait extends MultisampleMappingReportTrait {
39

40
  def variantcallingExecuted = summary.getValue("shiva", "settings", "multisample_variantcalling") match {
Peter van 't Hof's avatar
Peter van 't Hof committed
41
    case Some(true) => true
Peter van 't Hof's avatar
Peter van 't Hof committed
42
    case _          => false
43
44
  }

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

Peter van 't Hof's avatar
Peter van 't Hof committed
47
48
  override def pipelineName = "shiva"

49
50
51
52
  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)

53
  /** Root page for the shiva report */
54
  override def indexPage = {
55
    val regions = regionsPage
56
57
    val oldPage = super.indexPage

58
    oldPage.copy(subPages = oldPage.subPages ++ regionsPage)
59
  }
60

61
62
  /** Generate a page with all target coverage stats */
  def regionsPage: Option[(String, ReportPage)] = {
63
64
65
66
67
68
    val roi = summary.getValue("shiva", "settings", "regions_of_interest")
    val amplicon = summary.getValue("shiva", "settings", "amplicon_bed")

    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
69
      ReportPage(
70
        List(),
Peter van 't Hof's avatar
Peter van 't Hof committed
71
72
73
        List("Coverage" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/covstatsMultiTable.ssp")),
        Map("target" -> name)
      )
74
75
76
    }

    amplicon match {
77
      case Some(x: String) => regionPages += (x + " (Amplicon)") -> createPage(x, amplicon = true)
78
79
80
81
      case _               =>
    }

    roi match {
82
83
      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))
84
85
86
      case _                =>
    }

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

101
  /** Files page, can be used general or at sample level */
102
103
  override def filesPage: ReportPage = {
    val vcfFilesSection = if (variantcallingExecuted) List("VCF files" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/shiva/outputVcfFiles.ssp",
Peter van 't Hof's avatar
Peter van 't Hof committed
104
105
      Map("sampleId" -> None)))
    else Nil
106
107
    val oldPage = super.filesPage
    oldPage.copy(sections = oldPage.sections ++ vcfFilesSection)
108
109
  }

110
111
112
113
114
  /** Single sample page */
  override def samplePage(sampleId: String, args: Map[String, Any]): ReportPage = {
    val variantcallingSection = if (variantcallingExecuted) List("Variantcalling" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/shiva/sampleVariants.ssp")) else Nil
    val oldPage = super.samplePage(sampleId, args)
    oldPage.copy(sections = variantcallingSection ++ oldPage.sections)
115
116
  }

117
  /** Name of the report */
118
  def reportName = "Shiva Report"
119

120
121
  /**
   * Generate a stackbar plot for found variants
Peter van 't Hof's avatar
Peter van 't Hof committed
122
   *
123
124
125
126
127
128
   * @param outputDir OutputDir for the tsv and png file
   * @param prefix Prefix of the tsv and png file
   * @param summary Summary class
   * @param libraryLevel Default false, when set true plot will be based on library stats instead of sample stats
   * @param sampleId Default it selects all sampples, when sample is giving it limits to selected sample
   */
129
  def variantSummaryPlot(outputDir: File,
130
131
132
                         prefix: String,
                         summary: Summary,
                         libraryLevel: Boolean = false,
133
134
135
                         sampleId: Option[String] = None,
                         caller: String = "final",
                         target: Option[String] = None): Unit = {
136
137
138
139
140
141
142
    val tsvFile = new File(outputDir, prefix + ".tsv")
    val pngFile = new File(outputDir, prefix + ".png")
    val tsvWriter = new PrintWriter(tsvFile)
    if (libraryLevel) tsvWriter.print("Library") else tsvWriter.print("Sample")
    tsvWriter.println("\tHomVar\tHet\tHomRef\tNoCall")

    def getLine(summary: Summary, sample: String, lib: Option[String] = None): String = {
143
144
      val path = target match {
        case Some(t) => List("shivavariantcalling", "stats", s"multisample-vcfstats-$caller-$t", "genotype")
Peter van 't Hof's avatar
Peter van 't Hof committed
145
        case _       => List("shivavariantcalling", "stats", s"multisample-vcfstats-$caller", "genotype")
146
147
148
149
150
      }
      val homVar = new SummaryValue(path :+ "HomVar", summary, Some(sample), lib).value.getOrElse(0).toString.toLong
      val homRef = new SummaryValue(path :+ "HomRef", summary, Some(sample), lib).value.getOrElse(0).toString.toLong
      val noCall = new SummaryValue(path :+ "NoCall", summary, Some(sample), lib).value.getOrElse(0).toString.toLong
      val het = new SummaryValue(path :+ "Het", summary, Some(sample), lib).value.getOrElse(0).toString.toLong
151
152
153
154
155
156
157
158
159
160
      val sb = new StringBuffer()
      if (lib.isDefined) sb.append(sample + "-" + lib.get + "\t") else sb.append(sample + "\t")
      sb.append(homVar + "\t")
      sb.append(het + "\t")
      sb.append(homRef + "\t")
      sb.append(noCall)
      sb.toString
    }

    if (libraryLevel) {
161
      for (
162
        sample <- summary.samples if sampleId.isEmpty || sample == sampleId.get;
163
164
        lib <- summary.libraries(sample)
      ) {
165
166
167
        tsvWriter.println(getLine(summary, sample, Some(lib)))
      }
    } else {
168
      for (sample <- summary.samples if sampleId.isEmpty || sample == sampleId.get) {
169
170
171
172
173
174
175
176
177
178
        tsvWriter.println(getLine(summary, sample))
      }
    }

    tsvWriter.close()

    val plot = new StackedBarPlot(null)
    plot.input = tsvFile
    plot.output = pngFile
    plot.ylabel = Some("VCF records")
179
180
    if (libraryLevel) {
      plot.width = Some(200 + (summary.libraries.filter(s => sampleId.getOrElse(s._1) == s._1).foldLeft(0)(_ + _._2.size) * 10))
181
    } else plot.width = Some(200 + (summary.samples.count(s => sampleId.getOrElse(s) == s) * 10))
182
183
    plot.runLocal()
  }
184
}