GearsSingle.scala 8.46 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
/**
 * 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
 *
11
 * A dual licensing mode is applied. The source code within this project is freely available for non-commercial use under an AGPL
12
13
14
15
16
 * license; For commercial users or users who do not want to follow the AGPL
 * license, please contact us to obtain a separate license.
 */
package nl.lumc.sasc.biopet.pipelines.gears

Peter van 't Hof's avatar
Peter van 't Hof committed
17
import nl.lumc.sasc.biopet.core.summary.SummaryQScript
Peter van 't Hof's avatar
Peter van 't Hof committed
18
import nl.lumc.sasc.biopet.core.BiopetQScript.InputFile
19
20
import nl.lumc.sasc.biopet.core.{PipelineCommand, SampleLibraryTag}
import nl.lumc.sasc.biopet.extensions.{Gzip, Zcat}
Peter van 't Hof's avatar
Peter van 't Hof committed
21
import nl.lumc.sasc.biopet.pipelines.flexiprep.Flexiprep
Peter van 't Hof's avatar
Peter van 't Hof committed
22
import nl.lumc.sasc.biopet.utils.Logging
Peter van 't Hof's avatar
Peter van 't Hof committed
23
import nl.lumc.sasc.biopet.utils.config.Configurable
24
25
import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb
import nl.lumc.sasc.biopet.utils.summary.db.SummaryDb.{SampleId, SampleName}
26
import org.broadinstitute.gatk.queue.QScript
27
28
29
30
import sun.security.provider.JavaKeyStore.DualFormatJKS

import scala.concurrent.Await
import scala.concurrent.duration.Duration
31

32
/**
Wai Yi Leung's avatar
Wai Yi Leung committed
33
 * Created by wyleung
34
 */
Peter van 't Hof's avatar
Peter van 't Hof committed
35
class GearsSingle(val parent: Configurable) extends QScript with SummaryQScript with SampleLibraryTag {
36
  def this() = this(null)
37

Wai Yi Leung's avatar
Wai Yi Leung committed
38
  @Input(doc = "R1 reads in FastQ format", shortName = "R1", required = false)
Peter van 't Hof's avatar
Peter van 't Hof committed
39
  var fastqR1: List[File] = Nil
40

Wai Yi Leung's avatar
Wai Yi Leung committed
41
  @Input(doc = "R2 reads in FastQ format", shortName = "R2", required = false)
Peter van 't Hof's avatar
Peter van 't Hof committed
42
  var fastqR2: List[File] = Nil
43

Wai Yi Leung's avatar
Wai Yi Leung committed
44
  @Input(doc = "All unmapped reads will be extracted from this bam for analysis", shortName = "bam", required = false)
Peter van 't Hof's avatar
Peter van 't Hof committed
45
  var bamFile: Option[File] = None
46

Peter van 't Hof's avatar
Peter van 't Hof committed
47
48
  @Argument(required = false)
  var outputName: String = _
49

Peter van 't Hof's avatar
Peter van 't Hof committed
50
51
  lazy val krakenScript = if (config("gears_use_kraken", default = false)) Some(new GearsKraken(this)) else None
  lazy val centrifugeScript = if (config("gears_use_centrifuge", default = true)) Some(new GearsCentrifuge(this)) else None
Peter van 't Hof's avatar
Peter van 't Hof committed
52
53
  lazy val qiimeRatx = if (config("gears_use_qiime_rtax", default = false)) Some(new GearsQiimeRtax(this)) else None
  lazy val qiimeClosed = if (config("gears_use_qiime_closed", default = false)) Some(new GearsQiimeClosed(this)) else None
Peter van 't Hof's avatar
Peter van 't Hof committed
54
  lazy val qiimeOpen = if (config("gears_use_qiime_open", default = false)) Some(new GearsQiimeOpen(this)) else None
Peter van 't Hof's avatar
Peter van 't Hof committed
55
  lazy val seqCount = if (config("gears_use_seq_count", default = false)) Some(new GearsSeqCount(this)) else None
Peter van 't Hof's avatar
Peter van 't Hof committed
56

Peter van 't Hof's avatar
Peter van 't Hof committed
57
58
  /** Executed before running the script */
  def init(): Unit = {
Peter van 't Hof's avatar
Peter van 't Hof committed
59
60
61
    if (fastqR1.isEmpty && !bamFile.isDefined) Logging.addError("Please specify fastq-file(s) or bam file")
    if (fastqR1.nonEmpty == bamFile.isDefined) Logging.addError("Provide either a bam file or a R1/R2 file")
    if (fastqR2.nonEmpty && fastqR1.size != fastqR2.size) Logging.addError("R1 and R2 has not the same number of files")
Peter van 't Hof's avatar
Peter van 't Hof committed
62
    if (sampleId == null || sampleId == None) Logging.addError("Missing sample ID on GearsSingle module")
Peter van 't Hof's avatar
Peter van 't Hof committed
63

64
65
66
67
68
69
70
    if (!skipFlexiprep) {
      val db = SummaryDb.openSqliteSummary(summaryDbFile)
      val sample = Await.result(db.getSamples(runId = summaryRunId, name = sampleId).map(_.headOption), Duration.Inf)
      val sId = sample.map(_.id).getOrElse(Await.result(db.createSample(sampleId.getOrElse("noSample"), summaryRunId), Duration.Inf))
      val library = Await.result(db.getLibraries(runId = summaryRunId, name = libId, sampleId = Some(sId)).map(_.headOption), Duration.Inf)
      val lId = library.map(_.id).getOrElse(Await.result(db.createLibrary(libId.getOrElse("noLib"), summaryRunId, sId), Duration.Inf))
    }
Peter van 't Hof's avatar
Peter van 't Hof committed
71
    if (outputName == null) {
Ruben Vorderman's avatar
Ruben Vorderman committed
72
      outputName = sampleId.getOrElse("noName") + libId.map("-" + _).getOrElse("")
73
    }
Peter van 't Hof's avatar
Peter van 't Hof committed
74

Peter van 't Hof's avatar
Peter van 't Hof committed
75
    if (fastqR1.nonEmpty) {
Peter van 't Hof's avatar
Peter van 't Hof committed
76
77
      fastqR1.foreach(inputFiles :+= InputFile(_))
      fastqR2.foreach(inputFiles :+= InputFile(_))
Peter van 't Hof's avatar
Peter van 't Hof committed
78
    } else bamFile.foreach(inputFiles :+= InputFile(_))
Peter van 't Hof's avatar
Peter van 't Hof committed
79
  }
80

Wai Yi Leung's avatar
Wai Yi Leung committed
81
  override def reportClass = {
82
    val gears = new GearsSingleReport(this)
Wai Yi Leung's avatar
Wai Yi Leung committed
83
    gears.outputDir = new File(outputDir, "report")
84
    gears.summaryDbFile = summaryDbFile
85
86
    sampleId.foreach(gears.args += "sampleId" -> _)
    libId.foreach(gears.args += "libId" -> _)
Wai Yi Leung's avatar
Wai Yi Leung committed
87
88
    Some(gears)
  }
89

Peter van 't Hof's avatar
Peter van 't Hof committed
90
91
  protected var skipFlexiprep: Boolean = config("skip_flexiprep", default = false)

Peter van 't Hof's avatar
Peter van 't Hof committed
92
93
94
  protected def executeFlexiprep(r1: List[File], r2: List[File]): (File, Option[File]) = {
    val read1: File = if (r1.size == 1) r1.head else {
      val outputFile = new File(outputDir, "merged.R1.fq.gz")
Peter van 't Hof's avatar
Peter van 't Hof committed
95
      add(Zcat(this, r1) | new Gzip(this) > outputFile)
Peter van 't Hof's avatar
Peter van 't Hof committed
96
97
98
99
100
      outputFile
    }

    val read2: Option[File] = if (r2.size <= 1) r2.headOption else {
      val outputFile = new File(outputDir, "merged.R2.fq.gz")
Peter van 't Hof's avatar
Peter van 't Hof committed
101
      add(Zcat(this, r2) | new Gzip(this) > outputFile)
Peter van 't Hof's avatar
Peter van 't Hof committed
102
103
104
      Some(outputFile)
    }

105
106
107
108
109
110
111
112
113
    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))
114
115
  }

116
117
118
119
  lazy protected val flexiprep: Option[Flexiprep] = if (!skipFlexiprep) {
    Some(new Flexiprep(this))
  } else None

Peter van 't Hof's avatar
Peter van 't Hof committed
120
121
  /** Method to add jobs */
  def biopetScript(): Unit = {
Peter van 't Hof's avatar
Peter van 't Hof committed
122
    val (r1, r2): (File, Option[File]) = (fastqR1, fastqR2, bamFile) match {
Peter van 't Hof's avatar
Peter van 't Hof committed
123
      case (r1, _, _) if r1.nonEmpty => executeFlexiprep(r1, fastqR2)
Peter van 't Hof's avatar
Peter van 't Hof committed
124
125
      case (_, _, Some(bam)) =>
        val extract = new ExtractUnmappedReads(this)
Peter van 't Hof's avatar
Peter van 't Hof committed
126
        extract.outputDir = outputDir
Peter van 't Hof's avatar
Peter van 't Hof committed
127
        extract.bamFile = bam
Peter van 't Hof's avatar
Peter van 't Hof committed
128
        extract.outputName = outputName
Peter van 't Hof's avatar
Peter van 't Hof committed
129
        add(extract)
Peter van 't Hof's avatar
Peter van 't Hof committed
130
        executeFlexiprep(extract.fastqUnmappedR1 :: Nil, extract.fastqUnmappedR2.toList)
Peter van 't Hof's avatar
Peter van 't Hof committed
131
      case _ => throw new IllegalArgumentException("Missing input files")
Peter van 't Hof's avatar
Peter van 't Hof committed
132
133
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
134
135
136
    lazy val combinedFastq = {
      r2 match {
        case Some(r2) =>
137
138
139
140
141
142
          val combineReads = new CombineReads(this)
          combineReads.outputDir = new File(outputDir, "combine_reads")
          combineReads.fastqR1 = r1
          combineReads.fastqR2 = r2
          add(combineReads)
          combineReads.combinedFastq
Peter van 't Hof's avatar
Peter van 't Hof committed
143
144
145
146
        case _ => r1
      }
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
147
    krakenScript foreach { kraken =>
Peter van 't Hof's avatar
Peter van 't Hof committed
148
      kraken.outputDir = new File(outputDir, "kraken")
149
150
      kraken.fastqR1 = r1
      kraken.fastqR2 = r2
Peter van 't Hof's avatar
Peter van 't Hof committed
151
      kraken.outputName = outputName
Peter van 't Hof's avatar
Peter van 't Hof committed
152
      add(kraken)
Peter van 't Hof's avatar
Peter van 't Hof committed
153
154
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
155
156
157
158
159
160
    centrifugeScript foreach { centrifuge =>
      centrifuge.outputDir = new File(outputDir, "centrifuge")
      centrifuge.fastqR1 = r1
      centrifuge.fastqR2 = r2
      centrifuge.outputName = outputName
      add(centrifuge)
161
      outputFiles += "centrifuge_output" -> centrifuge.centrifugeOutput
Peter van 't Hof's avatar
Peter van 't Hof committed
162
163
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
164
    qiimeRatx foreach { qiimeRatx =>
Peter van 't Hof's avatar
Peter van 't Hof committed
165
      qiimeRatx.outputDir = new File(outputDir, "qiime_rtax")
166
167
      qiimeRatx.fastqR1 = r1
      qiimeRatx.fastqR2 = r2
Peter van 't Hof's avatar
Peter van 't Hof committed
168
      add(qiimeRatx)
Peter van 't Hof's avatar
Peter van 't Hof committed
169
    }
170

Peter van 't Hof's avatar
Peter van 't Hof committed
171
    qiimeClosed foreach { qiimeClosed =>
Peter van 't Hof's avatar
Peter van 't Hof committed
172
      qiimeClosed.outputDir = new File(outputDir, "qiime_closed")
Peter van 't Hof's avatar
Peter van 't Hof committed
173
      qiimeClosed.fastqInput = combinedFastq
Peter van 't Hof's avatar
Peter van 't Hof committed
174
      add(qiimeClosed)
175
      outputFiles += "qiime_closed_otu_table" -> qiimeClosed.otuTable
176
177
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
178
179
180
181
    qiimeOpen foreach { qiimeOpen =>
      qiimeOpen.outputDir = new File(outputDir, "qiime_open")
      qiimeOpen.fastqInput = combinedFastq
      add(qiimeOpen)
182
      outputFiles += "qiime_open_otu_table" -> qiimeOpen.otuTable
Peter van 't Hof's avatar
Peter van 't Hof committed
183
184
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
185
186
187
188
    seqCount.foreach { seqCount =>
      seqCount.fastqInput = combinedFastq
      seqCount.outputDir = new File(outputDir, "seq_count")
      add(seqCount)
189
      outputFiles += "seq_count_count_file" -> seqCount.countFile
Peter van 't Hof's avatar
Peter van 't Hof committed
190
191
    }

Wai Yi Leung's avatar
Wai Yi Leung committed
192
    addSummaryJobs()
193
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
194

195
  /** Pipeline settings shown in the summary file */
Peter van 't Hof's avatar
Peter van 't Hof committed
196
197
198
199
  def summarySettings: Map[String, Any] = Map(
    "skip_flexiprep" -> skipFlexiprep,
    "gears_use_kraken" -> krakenScript.isDefined,
    "gear_use_qiime_rtax" -> qiimeRatx.isDefined,
Peter van 't Hof's avatar
Peter van 't Hof committed
200
201
    "gear_use_qiime_closed" -> qiimeClosed.isDefined,
    "gear_use_qiime_open" -> qiimeOpen.isDefined
Peter van 't Hof's avatar
Peter van 't Hof committed
202
  )
Peter van 't Hof's avatar
Peter van 't Hof committed
203

204
  /** Statistics shown in the summary file */
Wai Yi Leung's avatar
Wai Yi Leung committed
205
206
  def summaryFiles: Map[String, File] = Map.empty ++
    (if (bamFile.isDefined) Map("input_bam" -> bamFile.get) else Map()) ++
Peter van 't Hof's avatar
Peter van 't Hof committed
207
208
    fastqR1.zipWithIndex.map(x => s"input_R1_${x._2}" -> x._1) ++
    fastqR2.zipWithIndex.map(x => s"input_R2_${x._2}" -> x._1) ++
Wai Yi Leung's avatar
Wai Yi Leung committed
209
    outputFiles
210
211
212
}

/** This object give a default main method to the pipelines */
213
object GearsSingle extends PipelineCommand