BammetricsReport.scala 5.52 KB
Newer Older
1
2
package nl.lumc.sasc.biopet.pipelines.bammetrics

3
import java.io.{ PrintWriter, File }
4

5
import nl.lumc.sasc.biopet.core.report.{ ReportBuilder, ReportPage, ReportSection }
6
import nl.lumc.sasc.biopet.core.summary.{ SummaryValue, Summary }
Peter van 't Hof's avatar
Peter van 't Hof committed
7
import nl.lumc.sasc.biopet.extensions.rscript.{XYPlot, StackedBarPlot}
8
9
10
11
12

/**
 * Created by pjvan_thof on 3/30/15.
 */
object BammetricsReport extends ReportBuilder {
13
14
  // FIXME: Not yet finished

15
16
17
18
19
20
21
22
23
24
25
26
27
28
  val reportName = "Bam Metrics"

  def indexPage = ReportPage(Map(
    "Bam Metrics" -> bamMetricsPage
  ), List(
    "Report" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/bamMetricsFront.ssp")
  ),
    Map()
  )

  def bamMetricsPage = ReportPage(
    Map(),
    List(
      "Summary" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/alignmentSummary.ssp"),
Peter van 't Hof's avatar
Peter van 't Hof committed
29
      "Bam Stats" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/bamStats.ssp"),
30
      "Insert Size" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/insertSize.ssp"),
31
32
      "RNA (optional)" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/rna.ssp"),
      "Target (optional)" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/target.ssp"),
33
      "GC Bias" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/bammetrics/gcBias.ssp")
34
    ),
35
36
37
    Map()
  )

38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
  def alignmentSummaryPlot(outputDir: File,
                           prefix: String,
                           summary: Summary,
                           libraryLevel: Boolean = false,
                           sampleId: Option[String] = None): Unit = {
    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")

    def getLine(summary: Summary, sample: String, lib: Option[String] = None): String = {
      val mapped = new SummaryValue(List("bammetrics", "stats", "biopet_flagstat", "Mapped"),
        summary, Some(sample), lib).value.getOrElse(0).toString.toLong
      val duplicates = new SummaryValue(List("bammetrics", "stats", "biopet_flagstat", "Duplicates"),
        summary, Some(sample), lib).value.getOrElse(0).toString.toLong
      val total = new SummaryValue(List("bammetrics", "stats", "biopet_flagstat", "All"),
        summary, Some(sample), lib).value.getOrElse(0).toString.toLong
      val secondary = new SummaryValue(List("bammetrics", "stats", "biopet_flagstat", "NotPrimaryAlignment"),
        summary, Some(sample), lib).value.getOrElse(0).toString.toLong
      val sb = new StringBuffer()
      if (lib.isDefined) sb.append(sample + "-" + lib.get + "\t") else sb.append(sample + "\t")
      sb.append((mapped - duplicates) + "\t")
      sb.append(duplicates + "\t")
      sb.append((total - mapped - secondary) + "\t")
      sb.append(secondary)
      sb.toString
    }

    if (libraryLevel) {
68
69
70
71
      for (
        sample <- summary.samples if (sampleId.isEmpty || sample == sampleId.get);
        lib <- summary.libraries(sample)
      ) {
72
73
74
75
76
77
78
79
80
81
82
83
84
85
        tsvWriter.println(getLine(summary, sample, Some(lib)))
      }
    } else {
      for (sample <- summary.samples if (sampleId.isEmpty || sample == sampleId.get)) {
        tsvWriter.println(getLine(summary, sample))
      }
    }

    tsvWriter.close()

    val plot = new StackedBarPlot(null)
    plot.input = tsvFile
    plot.output = pngFile
    plot.ylabel = Some("Reads")
Peter van 't Hof's avatar
Peter van 't Hof committed
86
    plot.width = Some(1200)
Peter van 't Hof's avatar
Peter van 't Hof committed
87
    plot.title = Some("Aligned reads")
88
89
90
    plot.runLocal()
  }

Peter van 't Hof's avatar
Peter van 't Hof committed
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
  def insertSizePlot(outputDir: File,
                           prefix: String,
                           summary: Summary,
                           libraryLevel: Boolean = false,
                           sampleId: Option[String] = None): Unit = {
    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 {
      sampleId match {
        case Some(sample) => tsvWriter.println("\t" + sample)
        case _ => tsvWriter.println(summary.samples.mkString("Sample\t","\t", ""))
      }
    }

    var map: Map[Int, Map[String, Int]] = Map()

    def fill(sample: String, lib: Option[String]): Unit = {
      new SummaryValue(List("bammetrics", "stats", "CollectInsertSizeMetrics", "histogram", "content"),
        summary, Some(sample), lib).value.getOrElse(List(List("insert_size","All_Reads.fr_count"))) match {
        case l:List[_] => {
          l.tail.foreach(_ match {
            case l:List[_] => {
              val insertSize = l.head.toString.toInt
              val count = l.tail.head.toString.toInt
              val old = map.getOrElse(insertSize, Map())
              map += insertSize -> (old + ((sample + lib.getOrElse("")) -> count))
            }
            case _ => throw new IllegalStateException("Must be a list")
          })
        }
        case _ => throw new IllegalStateException("Must be a list")
      }
    }

    summary.samples.foreach(fill(_, None))

    for ((insertSize, counts) <- map) {
      tsvWriter.print(insertSize)
      for (sample <- summary.samples) {
        tsvWriter.print("\t" + counts.getOrElse(sample, "0"))
      }
      tsvWriter.println()
    }

    tsvWriter.close()

    val plot = new XYPlot(null)
    plot.input = tsvFile
    plot.output = pngFile
    plot.ylabel = Some("Reads")
    plot.xlabel = Some("Insertsize")
    plot.width = Some(1200)
    plot.title = Some("Insertsize")
    plot.runLocal()

  }
150
}