BammetricsReport.scala 22.8 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
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

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

Ruben Vorderman's avatar
Ruben Vorderman committed
72 73 74 75 76 77 78 79
  /** Generates values for bamMetricsPage */
  def bamMetricsPageValues(summary: SummaryDb,
                           sampleId: Option[Int],
                           libId: Option[Int],
                           metricsTag: String = "bammetrics"): Map[String, Any] = {
    Map("" -> "")
  }

80
  /** Generates a page with alignment stats */
Peter van 't Hof's avatar
Peter van 't Hof committed
81 82 83
  def bamMetricsPage(summary: SummaryDb,
                     sampleId: Option[Int],
                     libId: Option[Int],
Ruben Vorderman's avatar
Ruben Vorderman committed
84 85
                     metricsTag: String = "bammetrics"): Future[ReportPage] = {

86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
    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"))
      )
Peter van 't Hof's avatar
Peter van 't Hof committed
106
      .exists(_._2.isDefined)
107

108 109 110 111 112 113 114 115
    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"))
    )
116
    val targets = (
Peter van 't Hof's avatar
Peter van 't Hof committed
117 118
      targetSettings("amplicon_name"),
      targetSettings("roi_name")
Peter van 't Hof's avatar
Peter van 't Hof committed
119
    ) match {
120 121 122 123
      case (Some(amplicon: String), Some(roi: List[_])) => amplicon :: roi.map(_.toString)
      case (_, Some(roi: List[_])) => roi.map(_.toString)
      case _ => Nil
    }
124

Ruben Vorderman's avatar
Ruben Vorderman committed
125 126 127 128 129 130
    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

Ruben Vorderman's avatar
Ruben Vorderman committed
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
    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
    )
Ruben Vorderman's avatar
Ruben Vorderman committed
172 173 174 175 176 177 178 179 180 181 182 183
    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()))),
184
        List(
Ruben Vorderman's avatar
Ruben Vorderman committed
185
          "Summary" -> ReportSection(
Ruben Vorderman's avatar
Ruben Vorderman committed
186
            "/nl/lumc/sasc/biopet/pipelines/bammetrics/alignmentSummary.ssp",alignmentSummaryReportValues),
Ruben Vorderman's avatar
Ruben Vorderman committed
187 188
          "Mapping Quality" -> ReportSection(
            "/nl/lumc/sasc/biopet/pipelines/bammetrics/mappingQuality.ssp",
Ruben Vorderman's avatar
Ruben Vorderman committed
189
            mappingQualityReportValues),
Ruben Vorderman's avatar
Ruben Vorderman committed
190
          "Clipping" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/clipping.ssp",
Ruben Vorderman's avatar
Ruben Vorderman committed
191
            clippingReportValues)
Ruben Vorderman's avatar
Ruben Vorderman committed
192 193 194 195 196
        ) ++
          (if (insertsizeMetrics)
             List(
               "Insert Size" -> ReportSection(
                 "/nl/lumc/sasc/biopet/pipelines/bammetrics/insertSize.ssp",
Ruben Vorderman's avatar
Ruben Vorderman committed
197
                 insertSizeReportValues))
Ruben Vorderman's avatar
Ruben Vorderman committed
198 199 200 201
           else Nil) ++ (if (wgsExecuted)
                           List(
                             "Whole genome coverage" -> ReportSection(
                               "/nl/lumc/sasc/biopet/pipelines/bammetrics/wgsHistogram.ssp",
Ruben Vorderman's avatar
Ruben Vorderman committed
202
                               wgsHistogramReportValues))
Ruben Vorderman's avatar
Ruben Vorderman committed
203 204 205 206 207
                         else Nil) ++
          (if (rnaExecuted)
             List(
               "Rna coverage" -> ReportSection(
                 "/nl/lumc/sasc/biopet/pipelines/bammetrics/rnaHistogram.ssp",
Ruben Vorderman's avatar
Ruben Vorderman committed
208
                 rnaHistogramReportValues))
Ruben Vorderman's avatar
Ruben Vorderman committed
209 210 211 212
           else Nil),
        Map("metricsTag" -> metricsTag)
      )
    }
213 214
  }

215
  /**
216
    * Generates the lines for alignmentSummaryPlot
217 218 219
    *
    * @param summary Summary class
    * @param sampleId Default it selects all sampples, when sample is giving it limits to selected sample
220
    *                     * @param libraryLevel Default false, when set true plot will be based on library stats instead of sample stats
221
    */
222
  def alignmentSummaryPlotLines(summary: SummaryDb,
223
                                sampleId: Option[Int] = None,
224
                                libraryLevel: Boolean = false): Seq[String] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
225 226 227 228 229 230 231 232
    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) {
233 234
      summary
        .getStatsForLibraries(runId,
235 236 237 238
                              "bammetrics",
                              "bamstats",
                              sampleId = sampleId,
                              keyValues = statsPaths)
239 240 241 242
        .map(x => (x._1._1, Some(x._1._2)) -> x._2)
    } else
      summary
        .getStatsForSamples(runId,
243 244 245 246
                            "bammetrics",
                            "bamstats",
                            sample = sampleId.map(SampleId),
                            keyValues = statsPaths)
247
        .map(x => (x._1, None) -> x._2)
248
    val summaryPlotLines = ArrayBuffer[String]()
Peter van 't Hof's avatar
Peter van 't Hof committed
249

Peter van 't Hof's avatar
Peter van 't Hof committed
250
    for (((s, l), result) <- results) {
Peter van 't Hof's avatar
Peter van 't Hof committed
251
      val sampleName: String = summary.getSampleName(s).map(_.get)
252 253
      val libName: Option[String] =
        l.flatMap(x => Await.result(summary.getLibraryName(x), Duration.Inf))
254
      val sb = new StringBuffer()
255 256
      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
257 258 259 260
      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
261
      sb.append((mapped - duplicates - secondary) + "\t")
262
      sb.append(duplicates + "\t")
Peter van 't Hof's avatar
Peter van 't Hof committed
263
      sb.append((total - mapped) + "\t")
264
      sb.append(secondary)
265
      summaryPlotLines += sb.toString
266
    }
267 268 269 270 271 272 273 274 275 276 277 278 279 280
    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],
281
                           libraryLevel: Boolean = false): Unit = {
282 283 284 285 286
    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")
287

288 289 290
    for (line <- summaryPlotLines) {
      tsvWriter.println(line)
    }
291 292 293 294 295 296
    tsvWriter.close()

    val plot = new StackedBarPlot(null)
    plot.input = tsvFile
    plot.output = pngFile
    plot.ylabel = Some("Reads")
297
    plot.width = Some(200 + (summaryPlotLines.size * 10))
Peter van 't Hof's avatar
WIP  
Peter van 't Hof committed
298
    plot.title = Some("Aligned_reads")
299 300
    plot.runLocal()
  }
301

302
  /**
Ruben Vorderman's avatar
Ruben Vorderman committed
303
    * This is a generic method to create a Map that can be turned into a plot
304 305 306 307 308 309 310 311 312
    * @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
    */
313 314 315 316 317 318 319
  def summaryForPlot(summary: SummaryDb,
                     statsPaths: Map[String, List[String]],
                     yKeyList: List[String],
                     xKeyList: List[String],
                     pipeline: PipelineQuery,
                     module: ModuleQuery,
                     libraryLevel: Boolean = false,
320
                     sampleId: Option[Int] = None,
321
                     libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
322
    val results: Map[(Int, Option[Int]), Map[String, Option[Array[Any]]]] = if (libraryLevel) {
323 324 325 326 327 328 329 330
      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,
331 332 333 334
                            pipeline,
                            module,
                            sample = sampleId.map(SampleId),
                            keyValues = statsPaths)
335
        .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
336 337
    val tables: Array[Map[String, Array[Any]]] = results.map {
      case ((sample, library), map) =>
338 339
        val sampleName = Await
          .result(summary.getSampleName(sample), Duration.Inf)
Peter van 't Hof's avatar
Peter van 't Hof committed
340
          .getOrElse(throw new IllegalStateException("Sample must be there"))
341 342
        val libraryName =
          library.flatMap(l => Await.result(summary.getLibraryName(l), Duration.Inf))
Peter van 't Hof's avatar
Peter van 't Hof committed
343 344
        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
345
        Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
346 347 348 349
          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
350
        )
Peter van 't Hof's avatar
Peter van 't Hof committed
351
    }.toArray
352 353 354 355 356 357 358
    tables
  }

  /**
    * This is a generic method to create plots
    * @param outputDir Outputdir of the plot
    * @param prefix Files will start with this name
Ruben Vorderman's avatar
Ruben Vorderman committed
359
    * @param plotTables Tables to be written
360 361 362 363 364 365 366 367 368
    * @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,
Ruben Vorderman's avatar
Ruben Vorderman committed
369
                           plotTables: Array[Map[String, Array[Any]]],
370
                           yKeyList: List[String],
371
                           xKeyList: List[String],
372 373 374 375 376 377
                           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
378

Ruben Vorderman's avatar
Ruben Vorderman committed
379
    writeTableToTsv(tsvFile, mergeTables(plotTables, yKeyList.head), yKeyList.head)
Peter van 't Hof's avatar
Peter van 't Hof committed
380

381 382 383 384 385 386
    LinePlot(
      tsvFile,
      pngFile,
      xlabel = xlabel,
      ylabel = ylabel,
      title = title,
Ruben Vorderman's avatar
Ruben Vorderman committed
387
      hideLegend = plotTables.size > 40,
388 389 390
      /* changed from results.size. Original results in summaryForPlot*/
      removeZero = removeZero
    ).runLocal()
Peter van 't Hof's avatar
Peter van 't Hof committed
391 392
  }

393 394 395
  def insertSizePlotTables(summary: SummaryDb,
                           libraryLevel: Boolean = false,
                           sampleId: Option[Int] = None,
Ruben Vorderman's avatar
Ruben Vorderman committed
396 397 398 399 400
                           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")
    )
Ruben Vorderman's avatar
Ruben Vorderman committed
401 402 403 404 405 406 407 408 409
    summaryForPlot(summary,
                   statsPaths,
                   "insert_size" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "CollectInsertSizeMetrics",
                   libraryLevel,
                   sampleId,
                   libraryId)
Ruben Vorderman's avatar
Ruben Vorderman committed
410
  }
411

412
  /**
413 414 415 416
    * Generate a line plot for insertsize
    *
    * @param outputDir OutputDir for the tsv and png file
    * @param prefix Prefix of the tsv and png file
Ruben Vorderman's avatar
Ruben Vorderman committed
417
    * @param insertSizePlotTables Plot map generated by insertSizePlotTables function.
418
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
419
  def insertSizePlot(outputDir: File,
Peter van 't Hof's avatar
Peter van 't Hof committed
420
                     prefix: String,
421
                     insertSizePlotTables: Array[Map[String, Array[Any]]]): Unit = {
Ruben Vorderman's avatar
Ruben Vorderman committed
422 423 424 425 426 427 428 429
    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
430
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
431

432 433 434 435
  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
436 437 438
    val statsPaths = Map(
      "mapping_quality" -> List("mapping_quality", "histogram", "values"),
      "count" -> List("mapping_quality", "histogram", "counts")
439
    )
Ruben Vorderman's avatar
Ruben Vorderman committed
440 441 442 443 444 445 446 447 448
    summaryForPlot(summary,
                   statsPaths,
                   "mapping_quality" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "bamstats",
                   libraryLevel,
                   sampleId,
                   libraryId)
449 450
  }

451 452 453 454 455 456 457 458 459 460 461 462
  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")
  }
Ruben Vorderman's avatar
Ruben Vorderman committed
463

464 465 466
  def clippingPlotTables(summary: SummaryDb,
                         libraryLevel: Boolean = false,
                         sampleId: Option[Int] = None,
Ruben Vorderman's avatar
Ruben Vorderman committed
467
                         libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
468 469 470
    val statsPaths = Map(
      "clipping" -> List("clipping", "histogram", "values"),
      "count" -> List("clipping", "histogram", "counts")
471
    )
Ruben Vorderman's avatar
Ruben Vorderman committed
472 473 474 475 476 477 478 479 480
    summaryForPlot(summary,
                   statsPaths,
                   "clipping" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "bamstats",
                   libraryLevel,
                   sampleId,
                   libraryId)
481

Ruben Vorderman's avatar
Ruben Vorderman committed
482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503
  }
  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")
504
    )
Ruben Vorderman's avatar
Ruben Vorderman committed
505 506 507 508 509 510 511 512 513
    summaryForPlot(summary,
                   statsPaths,
                   "coverage" :: Nil,
                   "count" :: "high_quality_coverage_count" :: Nil,
                   "bammetrics",
                   "wgs",
                   libraryLevel,
                   sampleId,
                   libraryId)
514 515
  }

516
  /**
517 518 519 520
    * 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
Ruben Vorderman's avatar
Ruben Vorderman committed
521
    * @param wgsHistogramPlotTables Plot map generated by wgsHistogramPlotTables function.
522
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
523
  def wgsHistogramPlot(outputDir: File,
Peter van 't Hof's avatar
Peter van 't Hof committed
524
                       prefix: String,
Ruben Vorderman's avatar
Ruben Vorderman committed
525 526 527 528 529 530 531 532 533
                       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
534
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
535

Ruben Vorderman's avatar
Ruben Vorderman committed
536 537 538 539 540 541 542 543
  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")
    )
544 545 546 547 548 549 550 551 552
    summaryForPlot(summary,
                   statsPaths,
                   "position" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "rna",
                   libraryLevel,
                   sampleId,
                   libraryId)
Ruben Vorderman's avatar
Ruben Vorderman committed
553
  }
554

Peter van 't Hof's avatar
Peter van 't Hof committed
555
  /**
556 557 558 559
    * 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
Ruben Vorderman's avatar
Ruben Vorderman committed
560
    * @param rnaHistogramPlotTables Plot map generated by rnaHistogramPlotTables function.
561
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
562 563
  def rnaHistogramPlot(outputDir: File,
                       prefix: String,
Ruben Vorderman's avatar
Ruben Vorderman committed
564
                       rnaHistogramPlotTables: Array[Map[String, Array[Any]]],
Peter van 't Hof's avatar
Peter van 't Hof committed
565
                       libraryLevel: Boolean = false,
Peter van 't Hof's avatar
Peter van 't Hof committed
566 567
                       sampleId: Option[Int] = None,
                       libraryId: Option[Int] = None): Unit = {
Ruben Vorderman's avatar
Ruben Vorderman committed
568
    writePlotFromSummary(outputDir,
569 570 571 572 573 574 575
                         prefix,
                         rnaHistogramPlotTables,
                         "position" :: Nil,
                         "count" :: Nil,
                         "Relative position",
                         "Coverage",
                         "RNA coverage")
Peter van 't Hof's avatar
Peter van 't Hof committed
576 577
  }

Peter van 't Hof's avatar
Peter van 't Hof committed
578
  def mergeTables(tables: Array[Map[String, Array[Any]]],
579 580
                  mergeColumn: String,
                  defaultValue: Any = 0): Map[String, Array[Any]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
581 582
    val keys = tables.flatMap(x => x(mergeColumn)).distinct
    (for (table <- tables; (columnKey, columnValues) <- table if columnKey != mergeColumn) yield {
583 584
      columnKey -> keys.map(x =>
        table(mergeColumn).zip(columnValues).toMap.getOrElse(x, defaultValue))
Peter van 't Hof's avatar
Peter van 't Hof committed
585 586 587 588
    }).toMap + (mergeColumn -> keys)
  }

  def writeTableToTsv(tsvFile: File, table: Map[String, Array[Any]], firstColumn: String): Unit = {
589
    require(table.map(_._2.length).toList.distinct.size == 1,
590
            "Not all values has the same number or rows")
Peter van 't Hof's avatar
Peter van 't Hof committed
591 592 593
    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
594 595 596
    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
597 598 599
    }
    writer.close()
  }
Ruben Vorderman's avatar
Ruben Vorderman committed
600
}