BammetricsReport.scala 21.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 131 132 133 134 135 136 137 138 139 140 141 142 143 144
    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


    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()))),
145
        List(
Ruben Vorderman's avatar
Ruben Vorderman committed
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
          "Summary" -> ReportSection(
            "/nl/lumc/sasc/biopet/pipelines/bammetrics/alignmentSummary.ssp"),
          "Mapping Quality" -> ReportSection(
            "/nl/lumc/sasc/biopet/pipelines/bammetrics/mappingQuality.ssp",
            Map("showPlot" -> true)),
          "Clipping" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/clipping.ssp",
                                      Map("showPlot" -> true))
        ) ++
          (if (insertsizeMetrics)
             List(
               "Insert Size" -> ReportSection(
                 "/nl/lumc/sasc/biopet/pipelines/bammetrics/insertSize.ssp",
                 Map("showPlot" -> true)))
           else Nil) ++ (if (wgsExecuted)
                           List(
                             "Whole genome coverage" -> ReportSection(
                               "/nl/lumc/sasc/biopet/pipelines/bammetrics/wgsHistogram.ssp",
                               Map("showPlot" -> true)))
                         else Nil) ++
          (if (rnaExecuted)
             List(
               "Rna coverage" -> ReportSection(
                 "/nl/lumc/sasc/biopet/pipelines/bammetrics/rnaHistogram.ssp",
                 Map("showPlot" -> true)))
           else Nil),
        Map("metricsTag" -> metricsTag)
      )
    }
174 175
  }

176
  /**
177
    * Generates the lines for alignmentSummaryPlot
178 179 180
    *
    * @param summary Summary class
    * @param sampleId Default it selects all sampples, when sample is giving it limits to selected sample
181
    *                     * @param libraryLevel Default false, when set true plot will be based on library stats instead of sample stats
182
    */
183
  def alignmentSummaryPlotLines(summary: SummaryDb,
184
                                sampleId: Option[Int] = None,
185
                                libraryLevel: Boolean = false): Seq[String] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
186 187 188 189 190 191 192 193
    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) {
194 195
      summary
        .getStatsForLibraries(runId,
196 197 198 199
                              "bammetrics",
                              "bamstats",
                              sampleId = sampleId,
                              keyValues = statsPaths)
200 201 202 203
        .map(x => (x._1._1, Some(x._1._2)) -> x._2)
    } else
      summary
        .getStatsForSamples(runId,
204 205 206 207
                            "bammetrics",
                            "bamstats",
                            sample = sampleId.map(SampleId),
                            keyValues = statsPaths)
208
        .map(x => (x._1, None) -> x._2)
209
    val summaryPlotLines = ArrayBuffer[String]()
Peter van 't Hof's avatar
Peter van 't Hof committed
210

Peter van 't Hof's avatar
Peter van 't Hof committed
211
    for (((s, l), result) <- results) {
Peter van 't Hof's avatar
Peter van 't Hof committed
212
      val sampleName: String = summary.getSampleName(s).map(_.get)
213 214
      val libName: Option[String] =
        l.flatMap(x => Await.result(summary.getLibraryName(x), Duration.Inf))
215
      val sb = new StringBuffer()
216 217
      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
218 219 220 221
      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
222
      sb.append((mapped - duplicates - secondary) + "\t")
223
      sb.append(duplicates + "\t")
Peter van 't Hof's avatar
Peter van 't Hof committed
224
      sb.append((total - mapped) + "\t")
225
      sb.append(secondary)
226
      summaryPlotLines += sb.toString
227
    }
228 229 230 231 232 233 234 235 236 237 238 239 240 241
    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],
242
                           libraryLevel: Boolean = false): Unit = {
243 244 245 246 247
    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")
248

249 250 251
    for (line <- summaryPlotLines) {
      tsvWriter.println(line)
    }
252 253 254 255 256 257
    tsvWriter.close()

    val plot = new StackedBarPlot(null)
    plot.input = tsvFile
    plot.output = pngFile
    plot.ylabel = Some("Reads")
258
    plot.width = Some(200 + (summaryPlotLines.size * 10))
Peter van 't Hof's avatar
WIP  
Peter van 't Hof committed
259
    plot.title = Some("Aligned_reads")
260 261
    plot.runLocal()
  }
262

263
  /**
Ruben Vorderman's avatar
Ruben Vorderman committed
264
    * This is a generic method to create a Map that can be turned into a plot
265 266 267 268 269 270 271 272 273
    * @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
    */
274 275 276 277 278 279 280
  def summaryForPlot(summary: SummaryDb,
                     statsPaths: Map[String, List[String]],
                     yKeyList: List[String],
                     xKeyList: List[String],
                     pipeline: PipelineQuery,
                     module: ModuleQuery,
                     libraryLevel: Boolean = false,
281
                     sampleId: Option[Int] = None,
282
                     libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
283
    val results: Map[(Int, Option[Int]), Map[String, Option[Array[Any]]]] = if (libraryLevel) {
284 285 286 287 288 289 290 291
      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,
292 293 294 295
                            pipeline,
                            module,
                            sample = sampleId.map(SampleId),
                            keyValues = statsPaths)
296
        .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
297 298
    val tables: Array[Map[String, Array[Any]]] = results.map {
      case ((sample, library), map) =>
299 300
        val sampleName = Await
          .result(summary.getSampleName(sample), Duration.Inf)
Peter van 't Hof's avatar
Peter van 't Hof committed
301
          .getOrElse(throw new IllegalStateException("Sample must be there"))
302 303
        val libraryName =
          library.flatMap(l => Await.result(summary.getLibraryName(l), Duration.Inf))
Peter van 't Hof's avatar
Peter van 't Hof committed
304 305
        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
306
        Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
307 308 309 310
          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
311
        )
Peter van 't Hof's avatar
Peter van 't Hof committed
312
    }.toArray
313 314 315 316 317 318 319
    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
320
    * @param plotTables Tables to be written
321 322 323 324 325 326 327 328 329
    * @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
330
                           plotTables: Array[Map[String, Array[Any]]],
331
                           yKeyList: List[String],
332
                           xKeyList: List[String],
333 334 335 336 337 338
                           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
339

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

342 343 344 345 346 347
    LinePlot(
      tsvFile,
      pngFile,
      xlabel = xlabel,
      ylabel = ylabel,
      title = title,
Ruben Vorderman's avatar
Ruben Vorderman committed
348
      hideLegend = plotTables.size > 40,
349 350 351
      /* changed from results.size. Original results in summaryForPlot*/
      removeZero = removeZero
    ).runLocal()
Peter van 't Hof's avatar
Peter van 't Hof committed
352 353
  }

354 355 356
  def insertSizePlotTables(summary: SummaryDb,
                           libraryLevel: Boolean = false,
                           sampleId: Option[Int] = None,
Ruben Vorderman's avatar
Ruben Vorderman committed
357 358 359 360 361
                           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
362 363 364 365 366 367 368 369 370
    summaryForPlot(summary,
                   statsPaths,
                   "insert_size" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "CollectInsertSizeMetrics",
                   libraryLevel,
                   sampleId,
                   libraryId)
Ruben Vorderman's avatar
Ruben Vorderman committed
371
  }
372

373
  /**
374 375 376 377
    * 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
378
    * @param insertSizePlotTables Plot map generated by insertSizePlotTables function.
379
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
380
  def insertSizePlot(outputDir: File,
Peter van 't Hof's avatar
Peter van 't Hof committed
381
                     prefix: String,
382
                     insertSizePlotTables: Array[Map[String, Array[Any]]]): Unit = {
Ruben Vorderman's avatar
Ruben Vorderman committed
383 384 385 386 387 388 389 390
    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
391
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
392

393 394 395 396
  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
397 398 399
    val statsPaths = Map(
      "mapping_quality" -> List("mapping_quality", "histogram", "values"),
      "count" -> List("mapping_quality", "histogram", "counts")
400
    )
Ruben Vorderman's avatar
Ruben Vorderman committed
401 402 403 404 405 406 407 408 409
    summaryForPlot(summary,
                   statsPaths,
                   "mapping_quality" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "bamstats",
                   libraryLevel,
                   sampleId,
                   libraryId)
410 411
  }

412 413 414 415 416 417 418 419 420 421 422 423
  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
424

425 426 427
  def clippingPlotTables(summary: SummaryDb,
                         libraryLevel: Boolean = false,
                         sampleId: Option[Int] = None,
Ruben Vorderman's avatar
Ruben Vorderman committed
428
                         libraryId: Option[Int] = None): Array[Map[String, Array[Any]]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
429 430 431
    val statsPaths = Map(
      "clipping" -> List("clipping", "histogram", "values"),
      "count" -> List("clipping", "histogram", "counts")
432
    )
Ruben Vorderman's avatar
Ruben Vorderman committed
433 434 435 436 437 438 439 440 441
    summaryForPlot(summary,
                   statsPaths,
                   "clipping" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "bamstats",
                   libraryLevel,
                   sampleId,
                   libraryId)
442

Ruben Vorderman's avatar
Ruben Vorderman committed
443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464
  }
  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")
465
    )
Ruben Vorderman's avatar
Ruben Vorderman committed
466 467 468 469 470 471 472 473 474
    summaryForPlot(summary,
                   statsPaths,
                   "coverage" :: Nil,
                   "count" :: "high_quality_coverage_count" :: Nil,
                   "bammetrics",
                   "wgs",
                   libraryLevel,
                   sampleId,
                   libraryId)
475 476
  }

477
  /**
478 479 480 481
    * 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
482
    * @param wgsHistogramPlotTables Plot map generated by wgsHistogramPlotTables function.
483
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
484
  def wgsHistogramPlot(outputDir: File,
Peter van 't Hof's avatar
Peter van 't Hof committed
485
                       prefix: String,
Ruben Vorderman's avatar
Ruben Vorderman committed
486 487 488 489 490 491 492 493 494
                       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
495
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
496

Ruben Vorderman's avatar
Ruben Vorderman committed
497 498 499 500 501 502 503 504
  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")
    )
505 506 507 508 509 510 511 512 513
    summaryForPlot(summary,
                   statsPaths,
                   "position" :: Nil,
                   "count" :: Nil,
                   "bammetrics",
                   "rna",
                   libraryLevel,
                   sampleId,
                   libraryId)
Ruben Vorderman's avatar
Ruben Vorderman committed
514
  }
515

Peter van 't Hof's avatar
Peter van 't Hof committed
516
  /**
517 518 519 520
    * 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
521
    * @param rnaHistogramPlotTables Plot map generated by rnaHistogramPlotTables function.
522
    */
Peter van 't Hof's avatar
Peter van 't Hof committed
523 524
  def rnaHistogramPlot(outputDir: File,
                       prefix: String,
Ruben Vorderman's avatar
Ruben Vorderman committed
525
                       rnaHistogramPlotTables: Array[Map[String, Array[Any]]],
Peter van 't Hof's avatar
Peter van 't Hof committed
526
                       libraryLevel: Boolean = false,
Peter van 't Hof's avatar
Peter van 't Hof committed
527 528
                       sampleId: Option[Int] = None,
                       libraryId: Option[Int] = None): Unit = {
Ruben Vorderman's avatar
Ruben Vorderman committed
529
    writePlotFromSummary(outputDir,
530 531 532 533 534 535 536
                         prefix,
                         rnaHistogramPlotTables,
                         "position" :: Nil,
                         "count" :: Nil,
                         "Relative position",
                         "Coverage",
                         "RNA coverage")
Peter van 't Hof's avatar
Peter van 't Hof committed
537 538
  }

Peter van 't Hof's avatar
Peter van 't Hof committed
539
  def mergeTables(tables: Array[Map[String, Array[Any]]],
540 541
                  mergeColumn: String,
                  defaultValue: Any = 0): Map[String, Array[Any]] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
542 543
    val keys = tables.flatMap(x => x(mergeColumn)).distinct
    (for (table <- tables; (columnKey, columnValues) <- table if columnKey != mergeColumn) yield {
544 545
      columnKey -> keys.map(x =>
        table(mergeColumn).zip(columnValues).toMap.getOrElse(x, defaultValue))
Peter van 't Hof's avatar
Peter van 't Hof committed
546 547 548 549
    }).toMap + (mergeColumn -> keys)
  }

  def writeTableToTsv(tsvFile: File, table: Map[String, Array[Any]], firstColumn: String): Unit = {
550
    require(table.map(_._2.length).toList.distinct.size == 1,
551
            "Not all values has the same number or rows")
Peter van 't Hof's avatar
Peter van 't Hof committed
552 553 554
    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
555 556 557
    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
558 559 560
    }
    writer.close()
  }
Ruben Vorderman's avatar
Ruben Vorderman committed
561
}