BedtoolsCoverage.scala 3.4 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
 * license; For commercial users or users who do not want to follow the AGPL
 * license, please contact us to obtain a separate license.
 */
15
package nl.lumc.sasc.biopet.extensions.bedtools
Peter van 't Hof's avatar
Peter van 't Hof committed
16

Sander Bollen's avatar
Sander Bollen committed
17
import java.io.{ File, PrintWriter }
Peter van 't Hof's avatar
Peter van 't Hof committed
18

Sander Bollen's avatar
Sander Bollen committed
19
import nl.lumc.sasc.biopet.core.Reference
Peter van 't Hof's avatar
Peter van 't Hof committed
20
import nl.lumc.sasc.biopet.utils.config.Configurable
Peter van 't Hof's avatar
Peter van 't Hof committed
21
22
import org.broadinstitute.gatk.utils.commandline.{ Argument, Input, Output }

Sander Bollen's avatar
Sander Bollen committed
23
import scala.io.Source
Peter van 't Hof's avatar
Peter van 't Hof committed
24

25
/** Extension for bedtools coverage */
Sander Bollen's avatar
Sander Bollen committed
26
class BedtoolsCoverage(val root: Configurable) extends Bedtools with Reference {
27

bow's avatar
bow committed
28
  @Input(doc = "Input file (bed/gff/vcf/bam)")
29
  var input: File = null
bow's avatar
bow committed
30
31

  @Input(doc = "Intersect file (bed/gff/vcf)")
32
  var intersectFile: File = null
bow's avatar
bow committed
33
34

  @Output(doc = "output File")
35
  var output: File = null
bow's avatar
bow committed
36

bow's avatar
bow committed
37
  @Argument(doc = "depth", required = false)
Peter van 't Hof's avatar
Peter van 't Hof committed
38
  var depth: Boolean = false
Peter van 't Hof's avatar
Peter van 't Hof committed
39

Peter van 't Hof's avatar
Peter van 't Hof committed
40
41
  @Argument(doc = "sameStrand", required = false)
  var sameStrand: Boolean = false
Peter van 't Hof's avatar
Peter van 't Hof committed
42

Peter van 't Hof's avatar
Peter van 't Hof committed
43
44
  @Argument(doc = "diffStrand", required = false)
  var diffStrand: Boolean = false
bow's avatar
bow committed
45

Sander Bollen's avatar
Sander Bollen committed
46
  @Argument(doc = "sorted", required = false)
47
48
  var sorted: Boolean = config("sorted", default = false, freeVar = false)

Peter van 't Hof's avatar
Peter van 't Hof committed
49
  override def defaultCoreMemory = 4.0
50

51
  /** Returns command to execute */
bow's avatar
bow committed
52
  def cmdLine = required(executable) + required("coverage") +
53
    required("-a", input) +
bow's avatar
bow committed
54
55
    required("-b", intersectFile) +
    conditional(depth, "-d") +
Peter van 't Hof's avatar
Peter van 't Hof committed
56
57
    conditional(sameStrand, "-s") +
    conditional(diffStrand, "-S") +
58
    conditional(sorted, "-sorted") +
59
    (if (sorted) required("-g", BedtoolsCoverage.getGenomeFile(referenceFai, jobTempDir)) else "") +
60
    (if (outputAsStsout) "" else " > " + required(output))
Sander Bollen's avatar
Sander Bollen committed
61

Peter van 't Hof's avatar
Peter van 't Hof committed
62
63
64
}

object BedtoolsCoverage {
65
  /** Returns defaul bedtools coverage */
66
  def apply(root: Configurable, input: File, intersect: File, output: Option[File] = None,
67
            depth: Boolean = false, sameStrand: Boolean = false, diffStrand: Boolean = false): BedtoolsCoverage = {
Peter van 't Hof's avatar
Peter van 't Hof committed
68
69
70
    val bedtoolsCoverage = new BedtoolsCoverage(root)
    bedtoolsCoverage.input = input
    bedtoolsCoverage.intersectFile = intersect
71
    output.foreach(bedtoolsCoverage.output = _)
Peter van 't Hof's avatar
Peter van 't Hof committed
72
    bedtoolsCoverage.depth = depth
Peter van 't Hof's avatar
Peter van 't Hof committed
73
74
    bedtoolsCoverage.sameStrand = sameStrand
    bedtoolsCoverage.diffStrand = diffStrand
75
    bedtoolsCoverage
Peter van 't Hof's avatar
Peter van 't Hof committed
76
  }
77
78
79
80
81
82
83
84
85

  private var genomeCache: Map[(File, File), File] = Map()

  def getGenomeFile(fai: File, dir: File): File = {
    if (!genomeCache.contains((fai, dir))) genomeCache += (fai, dir) -> createGenomeFile(fai, dir)
    genomeCache((fai, dir))
  }

  /**
Sander Bollen's avatar
Sander Bollen committed
86
87
88
   * Creates the genome file. i.e. the first two columns of the fasta index
   * @return
   */
89
90
91
92
93
94
95
96
97
98
99
  def createGenomeFile(fai: File, dir: File): File = {
    val tmp = File.createTempFile(fai.getName, ".genome", dir)
    tmp.deleteOnExit()
    val writer = new PrintWriter(tmp)
    Source.fromFile(fai).
      getLines().
      map(s => s.split("\t").take(2).mkString("\t")).
      foreach(f => writer.println(f))
    writer.close()
    tmp
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
100
}