BammetricsReport.scala 26.1 KB
Newer Older
Peter van 't Hof's avatar
Peter van 't Hof committed
1
/**
2 3 4 5 6 7 8 9 10 11 12 13 14
  * 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 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.
  */
15 16
package nl.lumc.sasc.biopet.pipelines.bammetrics

17
import java.io.{File, PrintWriter}
18

Peter van 't Hof's avatar
Peter van 't Hof committed
19
import nl.lumc.sasc.biopet.utils.config.Configurable
Vorderman's avatar
Vorderman committed
20 21 22 23 24 25
import nl.lumc.sasc.biopet.core.report.{
  ReportBuilder,
  ReportBuilderExtension,
  ReportPage,
  ReportSection
}
Peter van 't Hof's avatar
Peter van 't Hof committed
26
import nl.lumc.sasc.biopet.utils.ConfigUtils
27
import nl.lumc.sasc.biopet.utils.rscript.{LinePlot, StackedBarPlot}
Peter van 't Hof's avatar
Peter van 't Hof committed
28
import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb
Peter van 't Hof's avatar
Peter van 't Hof committed
29 30
import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb.Implicts._
import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb._
31

Vorderman's avatar
Vorderman committed
32
import scala.collection.mutable.ArrayBuffer
33
import scala.concurrent.{Await, Future}
Peter van 't Hof's avatar
Peter van 't Hof committed
34
import scala.concurrent.duration.Duration
35

Peter van 't Hof's avatar
Peter van 't Hof committed
36
class BammetricsReport(val parent: Configurable) extends ReportBuilderExtension {
37
  def builder = BammetricsReport
38
}
39

40
/**
41 42 43 44
  * Object to create a report for [[BamMetrics]]
  *
  * Created by pjvan_thof on 3/30/15.
  */
45
object BammetricsReport extends ReportBuilder {
46

47
  /** Name of report */
48 49
  val reportName = "Bam Metrics"

50 51
  def pipelineName = "bammetrics"

52
  /** Root page for single BamMetrcis report */
Peter van 't Hof's avatar
Peter van 't Hof committed
53 54
  def indexPage: Future[ReportPage] =
    bamMetricsPage(summary, sampleId, libId).map { bamMetricsPage =>
55 56 57 58 59 60 61 62 63 64 65 66 67
      ReportPage(
        bamMetricsPage.subPages ::: List(
          "Versions" -> Future(
            ReportPage(List(),
                       List("Executables" -> ReportSection(
                         "/nl/lumc/sasc/biopet/core/report/executables.ssp")),
                       Map())),
          "Files" -> filesPage(sampleId, libId)
        ),
        List(
          "Report" -> ReportSection(
            "/nl/lumc/sasc/biopet/pipelines/bammetrics/bamMetricsFront.ssp")
        ) ::: bamMetricsPage.sections,
Peter van 't Hof's avatar
Peter van 't Hof committed
68 69 70
        Map()
      )
    }
71

Vorderman's avatar
Vorderman committed
72 73 74 75 76
  /** Generates values for bamMetricsPage */
  def bamMetricsPageValues(summary: SummaryDb,
                           sampleId: Option[Int],
                           libId: Option[Int],
                           metricsTag: String = "bammetrics"): Map[String, Any] = {
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
    val wgsExecuted = summary.getStatsSize(runId,
      metricsTag,
      "wgs",
      sample = sampleId.map(SampleId),
      library = libId.map(LibraryId)) >= 1
    val rnaExecuted = summary.getStatsSize(runId,
      metricsTag,
      "rna",
      sample = sampleId.map(SampleId),
      library = libId.map(LibraryId)) >= 1

    val insertsizeMetrics = summary
      .getStatKeys(
        runId,
        metricsTag,
        "CollectInsertSizeMetrics",
        sampleId.map(SampleId).getOrElse(NoSample),
        libId.map(LibraryId).getOrElse(NoLibrary),
        Map("metrics" -> List("metrics"))
      )
      .exists(_._2.isDefined)

    val targetSettings = summary.getSettingKeys(
      runId,
      metricsTag,
      NoModule,
      sample = sampleId.map(SampleId).getOrElse(NoSample),
      library = libId.map(LibraryId).getOrElse(NoLibrary),
      Map("amplicon_name" -> List("amplicon_name"), "roi_name" -> List("roi_name"))
    )
    val targets = (
      targetSettings("amplicon_name"),
      targetSettings("roi_name")
    ) match {
      case (Some(amplicon: String), Some(roi: List[_])) => amplicon :: roi.map(_.toString)
      case (_, Some(roi: List[_])) => roi.map(_.toString)
      case _ => Nil
    }

    val covstatsPlotValuesArray = ArrayBuffer[(String, Map[String,Any])]()
    for (t <- targets) {
      covstatsPlotValuesArray += Tuple2(t,BammetricsReportPage.covstatsPlotValues(summary, runId, sampleId, libId, Some(t)))
    }
    val covstatsPlotValuesList = covstatsPlotValuesArray.toList

    val alignmentSummaryReportValues = BammetricsReportPage.alignmentSummaryValues(
      summary,
      runId,
      samples,
      libraries,
      sampleId,
      libId
    )
    val mappingQualityReportValues = BammetricsReportPage.mappingQualityValues(
      summary,
      runId,
      samples,
      libraries,
      sampleId,
      libId,
      showPlot = true
    )
    val clippingReportValues = BammetricsReportPage.clippingValues(
      summary,
      runId,
      samples,
      libraries,
      sampleId,
      libId,
      showPlot= true
    )
    val insertSizeReportValues = BammetricsReportPage.insertSizeValues(
      summary,
      runId,
      samples,
      libraries,
      sampleId,
      libId,
      showPlot= true
    )
    val wgsHistogramReportValues = BammetricsReportPage.wgsHistogramValues(
      summary,runId, samples, libraries, sampleId, libId, showPlot = true
    )
    val rnaHistogramReportValues = BammetricsReportPage.rnaHistogramValues(
      summary,runId, samples, libraries, sampleId, libId, showPlot = true
    )
    Map("wgsExecuted" -> wgsExecuted,
      "rnaExecuted" -> rnaExecuted,
      "insertsizeMetrics" -> insertsizeMetrics,
      "targetSettings" -> targetSettings,
      "covstatsPlotValuesList" -> covstatsPlotValuesList,
      "alignmentSummaryReportValues" -> alignmentSummaryReportValues,
      "clippingReportValues" -> clippingReportValues,
      "insertSizeReportValues" -> insertSizeReportValues,
      "wgsHistogramReportValues" -> wgsHistogramReportValues,
      "rnaHistogramReportValues" -> rnaHistogramReportValues
    )

Vorderman's avatar
Vorderman committed
175 176
  }

177
  /** Generates a page with alignment stats */
Peter van 't Hof's avatar
Peter van 't Hof committed
178 179 180
  def bamMetricsPage(summary: SummaryDb,
                     sampleId: Option[Int],
                     libId: Option[Int],
Vorderman's avatar
Vorderman committed
181 182
                     metricsTag: String = "bammetrics"): Future[ReportPage] = {

183
    val wgsExecuted: Boolean = summary.getStatsSize(runId,
184 185 186 187
                                           metricsTag,
                                           "wgs",
                                           sample = sampleId.map(SampleId),
                                           library = libId.map(LibraryId)) >= 1
188
    val rnaExecuted: Boolean = summary.getStatsSize(runId,
189 190 191 192 193
                                           metricsTag,
                                           "rna",
                                           sample = sampleId.map(SampleId),
                                           library = libId.map(LibraryId)) >= 1

194
    val insertsizeMetrics: Boolean = summary
195 196 197 198 199 200 201 202
      .getStatKeys(
        runId,
        metricsTag,
        "CollectInsertSizeMetrics",
        sampleId.map(SampleId).getOrElse(NoSample),
        libId.map(LibraryId).getOrElse(NoLibrary),
        Map("metrics" -> List("metrics"))
      )
Peter van 't Hof's avatar
Peter van 't Hof committed
203
      .exists(_._2.isDefined)
204

205
    val targetSettings: Map[String, Option[Any]] = summary.getSettingKeys(
206 207 208 209 210 211 212
      runId,
      metricsTag,
      NoModule,
      sample = sampleId.map(SampleId).getOrElse(NoSample),
      library = libId.map(LibraryId).getOrElse(NoLibrary),
      Map("amplicon_name" -> List("amplicon_name"), "roi_name" -> List("roi_name"))
    )
213
    val targets = (
Peter van 't Hof's avatar
Peter van 't Hof committed
214 215
      targetSettings("amplicon_name"),
      targetSettings("roi_name")
Peter van 't Hof's avatar
Peter van 't Hof committed
216
    ) match {
217 218 219 220
      case (Some(amplicon: String), Some(roi: List[_])) => amplicon :: roi.map(_.toString)
      case (_, Some(roi: List[_])) => roi.map(_.toString)
      case _ => Nil
    }
221

Vorderman's avatar
Vorderman committed
222 223 224 225
    val covstatsPlotValuesArray = ArrayBuffer[(String, Map[String,Any])]()
      for (t <- targets) {
        covstatsPlotValuesArray += Tuple2(t,BammetricsReportPage.covstatsPlotValues(summary, runId, sampleId, libId, Some(t)))
      }
226
    val covstatsPlotValuesList: List[(String, Map[String, Any])] = covstatsPlotValuesArray.toList
Vorderman's avatar
Vorderman committed
227

228
    val alignmentSummaryReportValues: Map[String, Any] = BammetricsReportPage.alignmentSummaryValues(
Vorderman's avatar
Vorderman committed
229 230 231 232 233 234 235
      summary,
      runId,
      samples,
      libraries,
      sampleId,
      libId
    )
236
    val mappingQualityReportValues: Map[String, Any] = BammetricsReportPage.mappingQualityValues(
Vorderman's avatar
Vorderman committed
237 238 239 240 241 242 243 244
      summary,
      runId,
      samples,
      libraries,
      sampleId,
      libId,
      showPlot = true
    )
245
    val clippingReportValues: Map[String, Any] = BammetricsReportPage.clippingValues(
Vorderman's avatar
Vorderman committed
246 247 248 249 250 251 252 253
      summary,
      runId,
      samples,
      libraries,
      sampleId,
      libId,
      showPlot= true
    )
254
    val insertSizeReportValues: Map[String, Any] = BammetricsReportPage.insertSizeValues(
Vorderman's avatar
Vorderman committed
255 256 257 258 259 260 261 262
      summary,
      runId,
      samples,
      libraries,
      sampleId,
      libId,
      showPlot= true
    )
263
    val wgsHistogramReportValues: Map[String, Any] = BammetricsReportPage.wgsHistogramValues(
Vorderman's avatar
Vorderman committed
264 265
      summary,runId, samples, libraries, sampleId, libId, showPlot = true
    )
266
    val rnaHistogramReportValues: Map[String, Any] = BammetricsReportPage.rnaHistogramValues(
Vorderman's avatar
Vorderman committed
267 268
      summary,runId, samples, libraries, sampleId, libId, showPlot = true
    )
Vorderman's avatar
Vorderman committed
269 270 271 272 273 274 275 276 277 278 279 280
    Future {
      ReportPage(
        if (targets.isEmpty) List()
        else
          List(
            "Targets" -> Future.successful(
              ReportPage(
                List(),
                covstatsPlotValuesList.map(covstats =>
                  covstats._1 -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/covstatsPlot.ssp",
                    covstats._2)),
                Map()))),
281
        List(
Vorderman's avatar
Vorderman committed
282
          "Summary" -> ReportSection(
Vorderman's avatar
Vorderman committed
283
            "/nl/lumc/sasc/biopet/pipelines/bammetrics/alignmentSummary.ssp",alignmentSummaryReportValues),
Vorderman's avatar
Vorderman committed
284 285
          "Mapping Quality" -> ReportSection(
            "/nl/lumc/sasc/biopet/pipelines/bammetrics/mappingQuality.ssp",
Vorderman's avatar
Vorderman committed
286
            mappingQualityReportValues),
Vorderman's avatar
Vorderman committed
287
          "Clipping" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/clipping.ssp",
Vorderman's avatar
Vorderman committed
288
            clippingReportValues)
Vorderman's avatar
Vorderman committed
289 290 291 292 293
        ) ++
          (if (insertsizeMetrics)
             List(
               "Insert Size" -> ReportSection(
                 "/nl/lumc/sasc/biopet/pipelines/bammetrics/insertSize.ssp",
Vorderman's avatar
Vorderman committed
294
                 insertSizeReportValues))
Vorderman's avatar
Vorderman committed
295 296 297 298
           else Nil) ++ (if (wgsExecuted)
                           List(
                             "Whole genome coverage" -> ReportSection(
                               "/nl/lumc/sasc/biopet/pipelines/bammetrics/wgsHistogram.ssp",
Vorderman's avatar
Vorderman committed
299
                               wgsHistogramReportValues))
Vorderman's avatar
Vorderman committed
300 301 302 303 304
                         else Nil) ++
          (if (rnaExecuted)
             List(
               "Rna coverage" -> ReportSection(
                 "/nl/lumc/sasc/biopet/pipelines/bammetrics/rnaHistogram.ssp",
Vorderman's avatar
Vorderman committed
305
                 rnaHistogramReportValues))
Vorderman's avatar
Vorderman committed
306 307 308 309
           else Nil),
        Map("metricsTag" -> metricsTag)
      )
    }
310 311
  }

312
  /**
Vorderman's avatar
Vorderman committed
313
    * Generates the lines for alignmentSummaryPlot
314 315 316
    *
    * @param summary Summary class
    * @param sampleId Default it selects all sampples, when sample is giving it limits to selected sample
Vorderman's avatar
Vorderman committed
317
    *                     * @param libraryLevel Default false, when set true plot will be based on library stats instead of sample stats
318
    */
Vorderman's avatar
Vorderman committed
319
  def alignmentSummaryPlotLines(summary: SummaryDb,
320
                                sampleId: Option[Int] = None,
Vorderman's avatar
Vorderman committed
321
                                libraryLevel: Boolean = false): Seq[String] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
322 323 324 325 326 327 328 329
    val statsPaths = Map(
      "Mapped" -> List("flagstats", "Mapped"),
      "Duplicates" -> List("flagstats", "Duplicates"),
      "All" -> List("flagstats", "All"),
      "NotPrimaryAlignment" -> List("flagstats", "NotPrimaryAlignment")
    )

    val results: Map[(Int, Option[Int]), Map[String, Option[Any]]] = if (libraryLevel) {
330 331
      summary
        .getStatsForLibraries(runId,
332 333 334 335
                              "bammetrics",
                              "bamstats",
                              sampleId = sampleId,
                              keyValues = statsPaths)
336 337 338 339
        .map(x => (x._1._1, Some(x._1._2)) -> x._2)
    } else
      summary
        .getStatsForSamples(runId,
340 341 342 343
                            "bammetrics",
                            "bamstats",
                            sample = sampleId.map(SampleId),
                            keyValues = statsPaths)
344
        .map(x => (x._1, None) -> x._2)
Vorderman's avatar
Vorderman committed
345
    val summaryPlotLines = ArrayBuffer[String]()
Peter van 't Hof's avatar
Peter van 't Hof committed
346

Peter van 't Hof's avatar
Peter van 't Hof committed
347
    for (((s, l), result) <- results) {
Peter van 't Hof's avatar
Peter van 't Hof committed
348
      val sampleName: String = summary.getSampleName(s).map(_.get)
349 350
      val libName: Option[String] =
        l.flatMap(x => Await.result(summary.getLibraryName(x), Duration.Inf))
351
      val sb = new StringBuffer()
352 353
      if (libName.isDefined) sb.append(sampleName + "-" + libName.get + "\t")
      else sb.append(sampleName + "\t")
Peter van 't Hof's avatar
Peter van 't Hof committed
354 355 356 357
      val mapped = ConfigUtils.any2long(result("Mapped"))
      val duplicates = ConfigUtils.any2long(result("Duplicates"))
      val total = ConfigUtils.any2long(result("All"))
      val secondary = ConfigUtils.any2long(result("NotPrimaryAlignment"))
Peter van 't Hof's avatar
Peter van 't Hof committed
358
      sb.append((mapped - duplicates - secondary) + "\t")
359
      sb.append(duplicates + "\t")
Peter van 't Hof's avatar
Peter van 't Hof committed
360
      sb.append((total - mapped) + "\t")
361
      sb.append(secondary)
Vorderman's avatar
Vorderman committed
362
      summaryPlotLines += sb.toString
363
    }
Vorderman's avatar
Vorderman committed
364 365 366 367 368 369 370 371 372 373 374 375 376 377
    summaryPlotLines
  }

  /**
    * Generate a stackbar plot for alignment stats
    *
    * @param outputDir OutputDir for the tsv and png file
    * @param prefix Prefix of the tsv and png file
    * @param summaryPlotLines A sequence of strings written to the summary tsv
    * @param libraryLevel Default false, when set true plot will be based on library stats instead of sample stats
    */
  def alignmentSummaryPlot(outputDir: File,
                           prefix: String,
                           summaryPlotLines: Seq[String],
378
                           libraryLevel: Boolean = false): Unit = {
Vorderman's avatar
Vorderman committed
379 380 381 382 383
    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("\tMapped\tDuplicates\tUnmapped\tSecondary")
384

Vorderman's avatar
Vorderman committed
385 386 387
    for (line <- summaryPlotLines) {
      tsvWriter.println(line)
    }
388 389 390 391 392 393
    tsvWriter.close()

    val plot = new StackedBarPlot(null)
    plot.input = tsvFile
    plot.output = pngFile
    plot.ylabel = Some("Reads")
Vorderman's avatar
Vorderman committed
394
    plot.width = Some(200 + (summaryPlotLines.size * 10))
Peter van 't Hof's avatar
WIP  
Peter van 't Hof committed
395
    plot.title = Some("Aligned_reads")
396 397
    plot.runLocal()
  }
398

399
  /**
Vorderman's avatar
Vorderman committed
400
    * This is a generic method to create a Map that can be turned into a plot
401 402 403 404 405 406 407 408 409
    * @param libraryLevel If enabled the plots will show data per library
    * @param sampleId If set only this sample is shown
    * @param libraryId If set onlt this library is shown
    * @param statsPaths Paths in summary where the tables can be found
    * @param yKeyList Keys to search from, first has prio over second one
    * @param xKeyList Keys to search from, first has prio over second one
    * @param pipeline Query for the pipeline
    * @param module Query for the module
    */
410 411 412 413 414 415 416
  def summaryForPlot(summary: SummaryDb,
                     statsPaths: Map[String, List[String]],
                     yKeyList: List[String],
                     xKeyList: List[String],
                     pipeline: PipelineQuery,
                     module: ModuleQuery,
                     libraryLevel: Boolean = false,
417
                     sampleId: Option[Int] = None,
418
                     libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
419
    val results: Map[(Int, Option[Int]), Map[String, Option[Array[Any]]]] = if (libraryLevel) {
420 421 422 423 424 425 426 427
      summary
        .getStatsForLibraries(runId, pipeline, module, sampleId = sampleId, keyValues = statsPaths)
        .map(x =>
          (x._1._1, Some(x._1._2)) -> x._2.map(x =>
            x._1 -> x._2.map(ConfigUtils.any2list(_).toArray)))
    } else
      summary
        .getStatsForSamples(runId,
428 429 430 431
                            pipeline,
                            module,
                            sample = sampleId.map(SampleId),
                            keyValues = statsPaths)
432
        .map(x => (x._1, None) -> x._2.map(x => x._1 -> x._2.map(ConfigUtils.any2list(_).toArray)))
Peter van 't Hof's avatar
Peter van 't Hof committed
433 434
    val tables: Array[Map[String, Array[Any]]] = results.map {
      case ((sample, library), map) =>
435 436
        val sampleName = Await
          .result(summary.getSampleName(sample), Duration.Inf)
Peter van 't Hof's avatar
Peter van 't Hof committed
437
          .getOrElse(throw new IllegalStateException("Sample must be there"))
438 439
        val libraryName =
          library.flatMap(l => Await.result(summary.getLibraryName(l), Duration.Inf))
Peter van 't Hof's avatar
Peter van 't Hof committed
440 441
        val yKey = yKeyList.find(x => map.contains(x) && map(x).isDefined).getOrElse("none")
        val xKey = xKeyList.find(x => map.contains(x) && map(x).isDefined).getOrElse("none")
Peter van 't Hof's avatar
Peter van 't Hof committed
442
        Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
443 444 445 446
          yKeyList.head -> map.getOrElse(yKey, None).getOrElse(Array()),
          (sampleName + libraryName.map("-" + _).getOrElse("")) -> map
            .getOrElse(xKey, None)
            .getOrElse(Array())
Peter van 't Hof's avatar
Peter van 't Hof committed
447
        )
Peter van 't Hof's avatar
Peter van 't Hof committed
448
    }.toArray
449 450 451 452 453 454 455
    tables
  }

  /**
    * This is a generic method to create plots
    * @param outputDir Outputdir of the plot
    * @param prefix Files will start with this name
Vorderman's avatar
Vorderman committed
456
    * @param plotTables Tables to be written
457 458 459 460 461 462 463 464 465
    * @param yKeyList Keys to search from, first has prio over second one
    * @param xKeyList Keys to search from, first has prio over second one
    * @param xlabel X label shown on the plot
    * @param ylabel Y label shown on the plot
    * @param title Title of the plot
    * @param removeZero
    */
  def writePlotFromSummary(outputDir: File,
                           prefix: String,
Vorderman's avatar
Vorderman committed
466
                           plotTables: Array[Map[String, Array[Any]]],
467
                           yKeyList: List[String],
468
                           xKeyList: List[String],
469 470 471 472 473 474
                           xlabel: Option[String] = None,
                           ylabel: Option[String] = None,
                           title: Option[String] = None,
                           removeZero: Boolean = true): Unit = {
    val tsvFile = new File(outputDir, prefix + ".tsv")
    val pngFile = new File(outputDir, prefix + ".png")
Peter van 't Hof's avatar
Peter van 't Hof committed
475

Vorderman's avatar
Vorderman committed
476
    writeTableToTsv(tsvFile, mergeTables(plotTables, yKeyList.head), yKeyList.head)
Peter van 't Hof's avatar
Peter van 't Hof committed
477

478 479 480 481 482 483
    LinePlot(
      tsvFile,
      pngFile,
      xlabel = xlabel,
      ylabel = ylabel,
      title = title,
Vorderman's avatar
Vorderman committed
484
      hideLegend = plotTables.size > 40,
485 486 487
      /* changed from results.size. Original results in summaryForPlot*/
      removeZero = removeZero
    ).runLocal()
Peter van 't Hof's avatar
Peter van 't Hof committed
488 489
  }

490 491 492
  def insertSizePlotTables(summary: SummaryDb,
                           libraryLevel: Boolean = false,
                           sampleId: Option[Int] = None,
Vorderman's avatar
Vorderman committed
493 494 495 496 497
                           libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
    val statsPaths = Map(
      "insert_size" -> List("histogram", "insert_size"),
      "count" -> List("histogram", "All_Reads.fr_count")
    )
Vorderman's avatar
Vorderman committed
498 499 500 501 502 503 504 505 506
    summaryForPlot(summary,
                   statsPaths,
                   "insert_size" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "CollectInsertSizeMetrics",
                   libraryLevel,
                   sampleId,
                   libraryId)
Vorderman's avatar
Vorderman committed
507
  }
508

509
  /**
510 511 512 513
    * Generate a line plot for insertsize
    *
    * @param outputDir OutputDir for the tsv and png file
    * @param prefix Prefix of the tsv and png file
Vorderman's avatar
Vorderman committed
514
    * @param insertSizePlotTables Plot map generated by insertSizePlotTables function.
515
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
516
  def insertSizePlot(outputDir: File,
Peter van 't Hof's avatar
Peter van 't Hof committed
517
                     prefix: String,
518
                     insertSizePlotTables: Array[Map[String, Array[Any]]]): Unit = {
Vorderman's avatar
Vorderman committed
519 520 521 522 523 524 525 526
    writePlotFromSummary(outputDir,
                         prefix,
                         insertSizePlotTables,
                         "insert_size" :: Nil,
                         "count" :: Nil,
                         "Insert size",
                         "Reads",
                         "Insert size")
Peter van 't Hof's avatar
Peter van 't Hof committed
527
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
528

529 530 531 532
  def mappingQualityPlotTables(summary: SummaryDb,
                               libraryLevel: Boolean = false,
                               sampleId: Option[Int] = None,
                               libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
533 534 535
    val statsPaths = Map(
      "mapping_quality" -> List("mapping_quality", "histogram", "values"),
      "count" -> List("mapping_quality", "histogram", "counts")
536
    )
Vorderman's avatar
Vorderman committed
537 538 539 540 541 542 543 544 545
    summaryForPlot(summary,
                   statsPaths,
                   "mapping_quality" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "bamstats",
                   libraryLevel,
                   sampleId,
                   libraryId)
546 547
  }

548 549 550 551 552 553 554 555 556 557 558 559
  def mappingQualityPlot(outputDir: File,
                         prefix: String,
                         mappingQualityPlotTables: Array[Map[String, Array[Any]]]): Unit = {
    writePlotFromSummary(outputDir,
                         prefix,
                         mappingQualityPlotTables,
                         "mapping_quality" :: Nil,
                         "count" :: Nil,
                         "Mapping Quality",
                         "Reads",
                         "Mapping Quality")
  }
Vorderman's avatar
Vorderman committed
560

561 562 563
  def clippingPlotTables(summary: SummaryDb,
                         libraryLevel: Boolean = false,
                         sampleId: Option[Int] = None,
Vorderman's avatar
Vorderman committed
564
                         libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
565 566 567
    val statsPaths = Map(
      "clipping" -> List("clipping", "histogram", "values"),
      "count" -> List("clipping", "histogram", "counts")
568
    )
Vorderman's avatar
Vorderman committed
569 570 571 572 573 574 575 576 577
    summaryForPlot(summary,
                   statsPaths,
                   "clipping" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "bamstats",
                   libraryLevel,
                   sampleId,
                   libraryId)
578

Vorderman's avatar
Vorderman committed
579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600
  }
  def clippingPlot(outputDir: File,
                   prefix: String,
                   clippingPlotTables: Array[Map[String, Array[Any]]]): Unit = {
    writePlotFromSummary(outputDir,
                         prefix,
                         clippingPlotTables,
                         "clipping" :: Nil,
                         "count" :: Nil,
                         "Clipping",
                         "Reads",
                         "Clipping")
  }

  def wgsHistogramPlotTables(summary: SummaryDb,
                             libraryLevel: Boolean = false,
                             sampleId: Option[Int] = None,
                             libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
    val statsPaths = Map(
      "coverage" -> List("histogram", "coverage"),
      "count" -> List("histogram", "count"),
      "high_quality_coverage_count" -> List("histogram", "high_quality_coverage_count")
601
    )
Vorderman's avatar
Vorderman committed
602 603 604 605 606 607 608 609 610
    summaryForPlot(summary,
                   statsPaths,
                   "coverage" :: Nil,
                   "count" :: "high_quality_coverage_count" :: Nil,
                   "bammetrics",
                   "wgs",
                   libraryLevel,
                   sampleId,
                   libraryId)
611 612
  }

613
  /**
614 615 616 617
    * Generate a line plot for wgs coverage
    *
    * @param outputDir OutputDir for the tsv and png file
    * @param prefix Prefix of the tsv and png file
Vorderman's avatar
Vorderman committed
618
    * @param wgsHistogramPlotTables Plot map generated by wgsHistogramPlotTables function.
619
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
620
  def wgsHistogramPlot(outputDir: File,
Peter van 't Hof's avatar
Peter van 't Hof committed
621
                       prefix: String,
Vorderman's avatar
Vorderman committed
622 623 624 625 626 627 628 629 630
                       wgsHistogramPlotTables: Array[Map[String, Array[Any]]]): Unit = {
    writePlotFromSummary(outputDir,
                         prefix,
                         wgsHistogramPlotTables,
                         "coverage" :: Nil,
                         "count" :: "high_quality_coverage_count" :: Nil,
                         "Coverage",
                         "Bases",
                         "Whole genome coverage")
Peter van 't Hof's avatar
Peter van 't Hof committed
631
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
632

Vorderman's avatar
Vorderman committed
633 634 635 636 637 638 639 640
  def rnaHistogramPlotTables(summary: SummaryDb,
                             libraryLevel: Boolean = false,
                             sampleId: Option[Int] = None,
                             libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
    val statsPaths = Map(
      "position" -> List("histogram", "normalized_position"),
      "count" -> List("histogram", "All_Reads.normalized_coverage")
    )
Vorderman's avatar
Vorderman committed
641 642 643 644 645 646 647 648 649
    summaryForPlot(summary,
                   statsPaths,
                   "position" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "rna",
                   libraryLevel,
                   sampleId,
                   libraryId)
Vorderman's avatar
Vorderman committed
650
  }
Vorderman's avatar
Vorderman committed
651

Peter van 't Hof's avatar
Peter van 't Hof committed
652
  /**
653 654 655 656
    * Generate a line plot for rna coverage
    *
    * @param outputDir OutputDir for the tsv and png file
    * @param prefix Prefix of the tsv and png file
Vorderman's avatar
Vorderman committed
657
    * @param rnaHistogramPlotTables Plot map generated by rnaHistogramPlotTables function.
Vorderman's avatar
Vorderman committed
658
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
659 660
  def rnaHistogramPlot(outputDir: File,
                       prefix: String,
Vorderman's avatar
Vorderman committed
661
                       rnaHistogramPlotTables: Array[Map[String, Array[Any]]],
Peter van 't Hof's avatar
Peter van 't Hof committed
662
                       libraryLevel: Boolean = false,
Peter van 't Hof's avatar
Peter van 't Hof committed
663 664
                       sampleId: Option[Int] = None,
                       libraryId: Option[Int] = None): Unit = {
Vorderman's avatar
Vorderman committed
665
    writePlotFromSummary(outputDir,
Vorderman's avatar
Vorderman committed
666 667 668 669 670 671 672
                         prefix,
                         rnaHistogramPlotTables,
                         "position" :: Nil,
                         "count" :: Nil,
                         "Relative position",
                         "Coverage",
                         "RNA coverage")
Peter van 't Hof's avatar
Peter van 't Hof committed
673 674
  }

Peter van 't Hof's avatar
Peter van 't Hof committed
675
  def mergeTables(tables: Array[Map[String, Array[Any]]],
676 677
                  mergeColumn: String,
                  defaultValue: Any = 0): Map[String, Array[Any]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
678 679
    val keys = tables.flatMap(x => x(mergeColumn)).distinct
    (for (table <- tables; (columnKey, columnValues) <- table if columnKey != mergeColumn) yield {
680 681
      columnKey -> keys.map(x =>
        table(mergeColumn).zip(columnValues).toMap.getOrElse(x, defaultValue))
Peter van 't Hof's avatar
Peter van 't Hof committed
682 683 684 685
    }).toMap + (mergeColumn -> keys)
  }

  def writeTableToTsv(tsvFile: File, table: Map[String, Array[Any]], firstColumn: String): Unit = {
686
    require(table.map(_._2.length).toList.distinct.size == 1,
687
            "Not all values has the same number or rows")
Peter van 't Hof's avatar
Peter van 't Hof committed
688 689 690
    val keys = table.keys.filterNot(_ == firstColumn).toList.sorted
    val writer = new PrintWriter(tsvFile)
    writer.println((firstColumn :: keys).mkString("\t"))
Peter van 't Hof's avatar
Peter van 't Hof committed
691 692 693
    table(firstColumn).zipWithIndex.foreach {
      case (c, i) =>
        writer.println((c :: keys.map(x => table(x)(i))).mkString("\t"))
Peter van 't Hof's avatar
Peter van 't Hof committed
694 695 696
    }
    writer.close()
  }
Vorderman's avatar
Vorderman committed
697
}