Merge pull request #13 from biopet/fix-BIOPET-525

Adding unassigned reads to krona plot

biopet-extensions/src/main/scala/nl/lumc/sasc/biopet/extensions/centrifuge/Centrifuge.scala

@@ -2,16 +2,19 @@ package nl.lumc.sasc.biopet.extensions.centrifuge
 
 import java.io.File
 
+import nl.lumc.sasc.biopet.core.summary.Summarizable
 import nl.lumc.sasc.biopet.core.{ BiopetCommandLineFunction, Version }
 import nl.lumc.sasc.biopet.utils.config.Configurable
+import nl.lumc.sasc.biopet.utils.tryToParseNumber
 import org.broadinstitute.gatk.utils.commandline.{ Input, Output }
 
+import scala.io.Source
 import scala.util.matching.Regex
 
 /**
  * Created by pjvanthof on 19/09/16.
  */
-class Centrifuge(val root: Configurable) extends BiopetCommandLineFunction with Version {
+class Centrifuge(val root: Configurable) extends BiopetCommandLineFunction with Version with Summarizable {
   @Input(doc = "Input: FastQ or FastA", required = true)
   var inputR1: File = _
 
@@ -129,4 +132,18 @@ class Centrifuge(val root: Configurable) extends BiopetCommandLineFunction with
     }) +
     (if (outputAsStsout) "" else required("-S", output)) +
     optional("--report-file", report)
+
+  /** Must return files to store into summary */
+  override def summaryFiles: Map[String, File] = metFile.map("metrics" -> _).toMap
+
+  /** Must returns stats to store into summary */
+  override def summaryStats: Any = {
+    metFile.map { file =>
+      val reader = Source.fromFile(file)
+      val header = reader.getLines().next().split("\t")
+      val values = reader.getLines().next().split("\t").map(tryToParseNumber(_, true).get)
+      reader.close()
+      Map("metrics" -> header.zip(values).toMap)
+    }.getOrElse(Map())
+  }
 }

gears/src/main/resources/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp

@@ -10,17 +10,23 @@
 <%@ var libId: Option[String] = None %>
 <%@ var args: Map[String, Any] %>
 <%@ var outputDir: File %>
+<%@ var centrifugeTag: Option[String] = None %>
 <%@ var summaryStatsTag: String = "krakenreport" %>
 <%@ var summaryModuleTag: String = "gearskraken" %>
 
 <%
-    val summaries = if (sampleId.isEmpty && libId.isEmpty) {
+    val summaries = if (sampleId.isEmpty) {
         summary.getSampleValues(summaryModuleTag, "stats", summaryStatsTag).map(x => x._1 -> x._2.get.asInstanceOf[Map[String, Any]])
     } else summary.getValue(sampleId, libId, summaryModuleTag, "stats", summaryStatsTag).map(sampleId.get -> _.asInstanceOf[Map[String, Any]]).toList.toMap
 
+    val totalReads = if (sampleId.isEmpty) {
+        centrifugeTag.map {tag => summary.getSampleValues(summaryModuleTag, "stats", tag, "metrics", "Read").map(x => x._1 -> x._2.getOrElse(0L).toString.toLong) }
+    } else centrifugeTag.flatMap(tag => summary.getValue(sampleId, libId, summaryModuleTag, "stats", tag, "metrics", "Read"))
+      .map(value => Map(sampleId.get -> value.toString.toLong))
+
     val tempFile = File.createTempFile("krona.", ".xml")
     tempFile.deleteOnExit()
-    GearsKraken.convertKrakenSummariesToKronaXml(summaries, tempFile)
+    GearsKraken.convertKrakenSummariesToKronaXml(summaries, tempFile, totalReads)
 
     val output = ReportBuilder.renderTemplate("/nl/lumc/sasc/biopet/core/report/krona.ssp",
     args ++ Map("kronaXml" -> tempFile))

gears/src/main/scala/nl/lumc/sasc/biopet/pipelines/gears/GearsCentrifuge.scala

@@ -38,6 +38,7 @@ class GearsCentrifuge(val root: Configurable) extends QScript with SummaryQScrip
     val centrifugeCmd = centrifuge | new Gzip(this) > centrifugeOutput
     centrifugeCmd.threadsCorrection = -1
     add(centrifugeCmd)
+    addSummarizable(centrifuge, "centrifuge")
 
     makeKreport("centrifuge", unique = false)
     makeKreport("centrifuge_unique", unique = true)

gears/src/main/scala/nl/lumc/sasc/biopet/pipelines/gears/GearsKraken.scala

@@ -120,7 +120,7 @@ object GearsKraken {
     convertKrakenSummariesToKronaXml(summaries, outputFile)
   }
 
-  def convertKrakenSummariesToKronaXml(summaries: Map[String, Map[String, Any]], outputFile: File): Unit = {
+  def convertKrakenSummariesToKronaXml(summaries: Map[String, Map[String, Any]], outputFile: File, totalReads: Option[Map[String, Long]] = None): Unit = {
 
     val samples = summaries.keys.toList.sorted
 
@@ -156,7 +156,7 @@ object GearsKraken {
             if (k == "root") {
               val unclassified = summaries(sample)("unclassified").asInstanceOf[Map[String, Any]]("size").asInstanceOf[Long]
               <val>
-                { getValue(sample, (path ::: k :: Nil).tail, "size").getOrElse(0).toString.toLong + unclassified }
+                { totalReads.flatMap(_.get(sample)).getOrElse(getValue(sample, (path ::: k :: Nil).tail, "size").getOrElse(0).toString.toLong + unclassified) }
               </val>
             } else {
               <val>

gears/src/main/scala/nl/lumc/sasc/biopet/pipelines/gears/GearsReport.scala

@@ -50,7 +50,7 @@ object GearsReport extends MultisampleReportBuilder {
       )), Map())), List(
         "Unique mappings" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp",
           Map("summaryStatsTag" -> "centrifuge_unique_report")
-        )), Map("summaryModuleTag" -> "gearscentrifuge")))
+        )), Map("summaryModuleTag" -> "gearscentrifuge", "centrifugeTag" -> Some("centrifuge"))))
       else Nil) ::: (if (krakenExecuted) List("Kraken analysis" -> ReportPage(List(), List(
         "Krona plot" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp"
         )), Map()))
@@ -92,9 +92,9 @@ object GearsReport extends MultisampleReportBuilder {
     )), Map())), List(
       "Unique mappings" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp",
         Map("summaryStatsTag" -> "centrifuge_unique_report")
-      )), Map("summaryModuleTag" -> "gearscentrifuge")))
-    else Nil) ::: (if (krakenExecuted) List("Kraken" -> ReportPage(List(), List(
-      "Kraken analysis" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp"
+      )), Map("summaryModuleTag" -> "gearscentrifuge", "centrifugeTag" -> Some("centrifuge"))))
+    else Nil) ::: (if (krakenExecuted) List("Kraken analysis" -> ReportPage(List(), List(
+      "Krona plot" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp"
       )), Map()))
     else Nil) ::: (if (qiimeClosesOtuTable.isDefined) List("Qiime closed reference analysis" -> ReportPage(List(), List(
       "Krona plot" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/qiimeKrona.ssp"
@@ -126,9 +126,9 @@ object GearsReport extends MultisampleReportBuilder {
       )), Map())), List(
         "Unique mappings" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp",
           Map("summaryStatsTag" -> "centrifuge_unique_report")
-        )), Map("summaryModuleTag" -> "gearscentrifuge")))
-      else Nil) ::: (if (krakenExecuted) List("Kraken" -> ReportPage(List(), List(
-        "Kraken analysis" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp"
+        )), Map("summaryModuleTag" -> "gearscentrifuge", "centrifugeTag" -> Some("centrifuge"))))
+      else Nil) ::: (if (krakenExecuted) List("Kraken analysis" -> ReportPage(List(), List(
+        "Krona plot" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp"
         )), Map()))
       else Nil) ::: (if (qiimeClosesOtuTable.isDefined) List("Qiime closed reference analysis" -> ReportPage(List(), List(
         "Krona plot" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/qiimeKrona.ssp"

gears/src/main/scala/nl/lumc/sasc/biopet/pipelines/gears/GearsSingle.scala

@@ -94,18 +94,21 @@ class GearsSingle(val root: Configurable) extends QScript with SummaryQScript wi
       Some(outputFile)
     }
 
-    if (!skipFlexiprep) {
-      val flexiprep = new Flexiprep(this)
-      flexiprep.inputR1 = read1
-      flexiprep.inputR2 = read2
-      flexiprep.sampleId = if (sampleId.isEmpty) Some("noSampleName") else sampleId
-      flexiprep.libId = if (libId.isEmpty) Some("noLibName") else libId
-      flexiprep.outputDir = new File(outputDir, "flexiprep")
-      add(flexiprep)
-      (flexiprep.fastqR1Qc, flexiprep.fastqR2Qc)
-    } else (read1, read2)
+    flexiprep.map { f =>
+      f.inputR1 = read1
+      f.inputR2 = read2
+      f.sampleId = Some(sampleId.getOrElse("noSampleName"))
+      f.libId = Some(libId.getOrElse("noLibName"))
+      f.outputDir = new File(outputDir, "flexiprep")
+      add(f)
+      (f.fastqR1Qc, f.fastqR2Qc)
+    }.getOrElse((read1, read2))
   }
 
+  lazy protected val flexiprep: Option[Flexiprep] = if (!skipFlexiprep) {
+    Some(new Flexiprep(this))
+  } else None
+
   /** Method to add jobs */
   def biopetScript(): Unit = {
     val (r1, r2): (File, Option[File]) = (fastqR1, fastqR2, bamFile) match {

gears/src/main/scala/nl/lumc/sasc/biopet/pipelines/gears/GearsSingleReport.scala

@@ -41,7 +41,7 @@ object GearsSingleReport extends ReportBuilder {
         else Nil) ++
         (if (centrifugeExecuted) List("Centrifuge analysis" ->
           ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp",
-            Map("summaryModuleTag" -> "gearscentrifuge", "summaryStatsTag" -> "centrifuge_unique_report")))
+            Map("summaryModuleTag" -> "gearscentrifuge", "summaryStatsTag" -> "centrifuge_unique_report", "centrifugeTag" -> Some("centrifuge"))))
         else Nil),
       pageArgs
     )

mapping/src/main/scala/nl/lumc/sasc/biopet/pipelines/mapping/MultisampleMappingReport.scala

@@ -65,7 +65,7 @@ trait MultisampleMappingReportTrait extends MultisampleReportBuilder {
         )), Map())), List(
           "Unique mappings" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp",
             Map("summaryStatsTag" -> "centrifuge_unique_report")
-          )), Map("summaryModuleTag" -> "gearscentrifuge")))
+          )), Map("summaryModuleTag" -> "gearscentrifuge", "centrifugeTag" -> Some("centrifuge"))))
         else Nil) ++
         List("Reference" -> ReportPage(List(), List(
           "Reference" -> ReportSection("/nl/lumc/sasc/biopet/core/report/reference.ssp", Map("pipeline" -> pipelineName))
@@ -128,7 +128,7 @@ trait MultisampleMappingReportTrait extends MultisampleReportBuilder {
       "Libraries" -> generateLibraryPage(args),
       "Alignment" -> BammetricsReport.bamMetricsPage(summary, Some(sampleId), None)) ++
       (if (centrifugeExecuted) List("Centriguge analysis" -> ReportPage(List("Non-unique" -> ReportPage(List(), List("All mappings" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp",
-        Map("summaryStatsTag" -> "centrifuge_report")
+        Map("summaryStatsTag" -> "centrifuge_report", "centrifugeTag" -> Some("centrifuge"))
       )), Map())), List(
         "Unique mappings" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp",
           Map("summaryStatsTag" -> "centrifuge_unique_report")
@@ -162,7 +162,7 @@ trait MultisampleMappingReportTrait extends MultisampleReportBuilder {
         )), Map())), List(
           "Unique mappings" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp",
             Map("summaryStatsTag" -> "centrifuge_unique_report")
-          )), Map("summaryModuleTag" -> "gearscentrifuge")))
+          )), Map("summaryModuleTag" -> "gearscentrifuge", "centrifugeTag" -> Some("centrifuge"))))
         else Nil) ::: (if (krakenExecuted) List("Dustbin analysis" -> ReportPage(List(), List(
           "Krona Plot" -> ReportSection("/nl/lumc/sasc/biopet/pipelines/gears/krakenKrona.ssp"
           )), Map()))