GearsSingle.scala 7.61 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
18
import nl.lumc.sasc.biopet.core.BiopetQScript.InputFile
Peter van 't Hof's avatar
Peter van 't Hof committed
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 org.broadinstitute.gatk.queue.QScript

26
/**
Wai Yi Leung's avatar
Wai Yi Leung committed
27
 * Created by wyleung
28
 */
Peter van 't Hof's avatar
Peter van 't Hof committed
29
class GearsSingle(val parent: Configurable) extends QScript with SummaryQScript with SampleLibraryTag {
30
  def this() = this(null)
31

Wai Yi Leung's avatar
Wai Yi Leung committed
32
  @Input(doc = "R1 reads in FastQ format", shortName = "R1", required = false)
33
  var fastqR1: List[File] = Nil
34

Wai Yi Leung's avatar
Wai Yi Leung committed
35
  @Input(doc = "R2 reads in FastQ format", shortName = "R2", required = false)
36
  var fastqR2: List[File] = Nil
37

Wai Yi Leung's avatar
Wai Yi Leung committed
38
  @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
39
  var bamFile: Option[File] = None
40

Peter van 't Hof's avatar
Peter van 't Hof committed
41 42
  @Argument(required = false)
  var outputName: String = _
43

Peter van 't Hof's avatar
Peter van 't Hof committed
44 45
  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
46 47
  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
48
  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
49
  lazy val seqCount = if (config("gears_use_seq_count", default = false)) Some(new GearsSeqCount(this)) else None
50

Peter van 't Hof's avatar
Peter van 't Hof committed
51 52
  /** Executed before running the script */
  def init(): Unit = {
53 54 55
    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
56
    if (sampleId == null || sampleId == None) Logging.addError("Missing sample ID on GearsSingle module")
Peter van 't Hof's avatar
Peter van 't Hof committed
57 58

    if (outputName == null) {
rhpvorderman's avatar
rhpvorderman committed
59
      outputName = sampleId.getOrElse("noName") + libId.map("-" + _).getOrElse("")
60
    }
61

62
    if (fastqR1.nonEmpty) {
63 64
      fastqR1.foreach(inputFiles :+= InputFile(_))
      fastqR2.foreach(inputFiles :+= InputFile(_))
65
    } else bamFile.foreach(inputFiles :+= InputFile(_))
Peter van 't Hof's avatar
Peter van 't Hof committed
66
  }
67

Wai Yi Leung's avatar
Wai Yi Leung committed
68
  override def reportClass = {
69
    val gears = new GearsSingleReport(this)
Wai Yi Leung's avatar
Wai Yi Leung committed
70
    gears.outputDir = new File(outputDir, "report")
71
    gears.summaryDbFile = summaryDbFile
72 73
    sampleId.foreach(gears.args += "sampleId" -> _)
    libId.foreach(gears.args += "libId" -> _)
Wai Yi Leung's avatar
Wai Yi Leung committed
74 75
    Some(gears)
  }
76

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

79 80 81
  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")
82
      add(Zcat(this, r1) | new Gzip(this) > outputFile)
83 84 85 86 87
      outputFile
    }

    val read2: Option[File] = if (r2.size <= 1) r2.headOption else {
      val outputFile = new File(outputDir, "merged.R2.fq.gz")
88
      add(Zcat(this, r2) | new Gzip(this) > outputFile)
89 90 91
      Some(outputFile)
    }

92 93 94 95 96 97 98 99 100
    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))
101 102
  }

103 104 105 106
  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
107 108
  /** Method to add jobs */
  def biopetScript(): Unit = {
Peter van 't Hof's avatar
Peter van 't Hof committed
109
    val (r1, r2): (File, Option[File]) = (fastqR1, fastqR2, bamFile) match {
110
      case (r1, _, _) if r1.nonEmpty => executeFlexiprep(r1, fastqR2)
111 112
      case (_, _, Some(bam)) =>
        val extract = new ExtractUnmappedReads(this)
113
        extract.outputDir = outputDir
114
        extract.bamFile = bam
115
        extract.outputName = outputName
Peter van 't Hof's avatar
Peter van 't Hof committed
116
        add(extract)
117
        executeFlexiprep(extract.fastqUnmappedR1 :: Nil, extract.fastqUnmappedR2.toList)
Peter van 't Hof's avatar
Peter van 't Hof committed
118
      case _ => throw new IllegalArgumentException("Missing input files")
119 120
    }

121 122 123
    lazy val combinedFastq = {
      r2 match {
        case Some(r2) =>
124 125 126 127 128 129
          val combineReads = new CombineReads(this)
          combineReads.outputDir = new File(outputDir, "combine_reads")
          combineReads.fastqR1 = r1
          combineReads.fastqR2 = r2
          add(combineReads)
          combineReads.combinedFastq
130 131 132 133
        case _ => r1
      }
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
134
    krakenScript foreach { kraken =>
135
      kraken.outputDir = new File(outputDir, "kraken")
136 137
      kraken.fastqR1 = r1
      kraken.fastqR2 = r2
138
      kraken.outputName = outputName
Peter van 't Hof's avatar
Peter van 't Hof committed
139
      add(kraken)
140 141
    }

142 143 144 145 146 147
    centrifugeScript foreach { centrifuge =>
      centrifuge.outputDir = new File(outputDir, "centrifuge")
      centrifuge.fastqR1 = r1
      centrifuge.fastqR2 = r2
      centrifuge.outputName = outputName
      add(centrifuge)
148
      outputFiles += "centrifuge_output" -> centrifuge.centrifugeOutput
149 150
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
151
    qiimeRatx foreach { qiimeRatx =>
Peter van 't Hof's avatar
Peter van 't Hof committed
152
      qiimeRatx.outputDir = new File(outputDir, "qiime_rtax")
153 154
      qiimeRatx.fastqR1 = r1
      qiimeRatx.fastqR2 = r2
Peter van 't Hof's avatar
Peter van 't Hof committed
155
      add(qiimeRatx)
156
    }
157

Peter van 't Hof's avatar
Peter van 't Hof committed
158
    qiimeClosed foreach { qiimeClosed =>
Peter van 't Hof's avatar
Peter van 't Hof committed
159
      qiimeClosed.outputDir = new File(outputDir, "qiime_closed")
160
      qiimeClosed.fastqInput = combinedFastq
Peter van 't Hof's avatar
Peter van 't Hof committed
161
      add(qiimeClosed)
162
      outputFiles += "qiime_closed_otu_table" -> qiimeClosed.otuTable
163 164
    }

165 166 167 168
    qiimeOpen foreach { qiimeOpen =>
      qiimeOpen.outputDir = new File(outputDir, "qiime_open")
      qiimeOpen.fastqInput = combinedFastq
      add(qiimeOpen)
169
      outputFiles += "qiime_open_otu_table" -> qiimeOpen.otuTable
170 171
    }

Peter van 't Hof's avatar
Peter van 't Hof committed
172 173 174 175
    seqCount.foreach { seqCount =>
      seqCount.fastqInput = combinedFastq
      seqCount.outputDir = new File(outputDir, "seq_count")
      add(seqCount)
176
      outputFiles += "seq_count_count_file" -> seqCount.countFile
Peter van 't Hof's avatar
Peter van 't Hof committed
177 178
    }

Wai Yi Leung's avatar
Wai Yi Leung committed
179
    addSummaryJobs()
180
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
181

182
  /** Pipeline settings shown in the summary file */
Peter van 't Hof's avatar
Peter van 't Hof committed
183 184 185 186
  def summarySettings: Map[String, Any] = Map(
    "skip_flexiprep" -> skipFlexiprep,
    "gears_use_kraken" -> krakenScript.isDefined,
    "gear_use_qiime_rtax" -> qiimeRatx.isDefined,
187 188
    "gear_use_qiime_closed" -> qiimeClosed.isDefined,
    "gear_use_qiime_open" -> qiimeOpen.isDefined
Peter van 't Hof's avatar
Peter van 't Hof committed
189
  )
Peter van 't Hof's avatar
Peter van 't Hof committed
190

191
  /** Statistics shown in the summary file */
192 193
  def summaryFiles: Map[String, File] = Map.empty ++
    (if (bamFile.isDefined) Map("input_bam" -> bamFile.get) else Map()) ++
194 195
    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
196
    outputFiles
197 198 199
}

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