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

Peter van 't Hof's avatar
Peter van 't Hof committed
18
import java.io.{ File, FileOutputStream }
19

20
import com.google.common.io.Files
Peter van 't Hof's avatar
Peter van 't Hof committed
21
import nl.lumc.sasc.biopet.extensions.gatk.{ BaseRecalibrator, IndelRealigner, PrintReads, RealignerTargetCreator }
22
import nl.lumc.sasc.biopet.extensions.picard.MarkDuplicates
23
import nl.lumc.sasc.biopet.extensions.tools.VcfStats
24
import nl.lumc.sasc.biopet.utils.ConfigUtils
25
import nl.lumc.sasc.biopet.utils.config.Config
26
27
28
import org.broadinstitute.gatk.queue.QSettings
import org.scalatest.Matchers
import org.scalatest.testng.TestNGSuite
Peter van 't Hof's avatar
Peter van 't Hof committed
29
import org.testng.annotations.{ DataProvider, Test }
30
31

/**
Peter van 't Hof's avatar
Peter van 't Hof committed
32
33
34
35
 * Class for testing shiva
 *
 * Created by pjvan_thof on 3/2/15.
 */
36
trait ShivaTestTrait extends TestNGSuite with Matchers {
37
38
  def initPipeline(map: Map[String, Any]): Shiva = {
    new Shiva() {
Sander Bollen's avatar
Sander Bollen committed
39
      override def configNamespace = "shiva"
40
      override def globalConfig = new Config(ConfigUtils.mergeMaps(map, ShivaTest.config))
41
42
43
44
45
46
47
      qSettings = new QSettings
      qSettings.runName = "test"
    }
  }

  @DataProvider(name = "shivaOptions")
  def shivaOptions = {
48
    for (
49
50
51
      s1 <- sample1; s2 <- sample2;
      realign <- realignProvider; baseRecalibration <- baseRecalibrationProvider
    ) yield Array("", s1, s2, realign, baseRecalibration)
52
53
  }

54
55
56
57
58
59
60
  def sample1 = Array(false, true)
  def sample2 = Array(false, true)
  def realignProvider = Array(false, true)
  def baseRecalibrationProvider = Array(false, true)
  def multisampleCalling: Boolean = true
  def sampleCalling = false
  def libraryCalling = false
Peter van 't Hof's avatar
Peter van 't Hof committed
61
62
  def dbsnp = true
  def svCalling = false
63
  def cnvCalling = false
Peter van 't Hof's avatar
Peter van 't Hof committed
64
  def annotation = false
65

66
  @Test(dataProvider = "shivaOptions")
67
68
  def testShiva(f: String, sample1: Boolean, sample2: Boolean,
                realign: Boolean, baseRecalibration: Boolean): Unit = {
69
70
    val map = {
      var m: Map[String, Any] = ShivaTest.config
Peter van 't Hof's avatar
Peter van 't Hof committed
71
72
      if (sample1) m = ConfigUtils.mergeMaps(ShivaTest.sample1, m)
      if (sample2) m = ConfigUtils.mergeMaps(ShivaTest.sample2, m)
73
      if (dbsnp) m = ConfigUtils.mergeMaps(Map("dbsnp" -> "test.vcf.gz"), m)
74
75
76
77
      ConfigUtils.mergeMaps(Map(
        "multisample_variantcalling" -> multisampleCalling,
        "single_sample_variantcalling" -> sampleCalling,
        "library_variantcalling" -> libraryCalling,
78
        "use_indel_realigner" -> realign,
Peter van 't Hof's avatar
Peter van 't Hof committed
79
80
        "use_base_recalibration" -> baseRecalibration,
        "sv_calling" -> svCalling,
81
        "cnv_calling" -> cnvCalling,
Peter van 't Hof's avatar
Peter van 't Hof committed
82
        "annotation" -> annotation), m)
83

84
85
    }

86
    if (!sample1 && !sample2) { // When no samples
87
88
89
90
91
92
93
      intercept[IllegalArgumentException] {
        initPipeline(map).script()
      }
    } else {
      val pipeline = initPipeline(map)
      pipeline.script()

94
95
96
97
      val numberLibs = (if (sample1) 1 else 0) + (if (sample2) 2 else 0)
      val numberSamples = (if (sample1) 1 else 0) + (if (sample2) 1 else 0)

      pipeline.functions.count(_.isInstanceOf[MarkDuplicates]) shouldBe (numberLibs + (if (sample2) 1 else 0))
98

99
100
101
      // Gatk preprocess
      pipeline.functions.count(_.isInstanceOf[IndelRealigner]) shouldBe (numberLibs * (if (realign) 1 else 0) + (if (sample2 && realign) 1 else 0))
      pipeline.functions.count(_.isInstanceOf[RealignerTargetCreator]) shouldBe (numberLibs * (if (realign) 1 else 0) + (if (sample2 && realign) 1 else 0))
Peter van 't Hof's avatar
Peter van 't Hof committed
102
      pipeline.functions.count(_.isInstanceOf[BaseRecalibrator]) shouldBe (if (dbsnp && baseRecalibration) (numberLibs * 2) else 0)
103
      pipeline.functions.count(_.isInstanceOf[PrintReads]) shouldBe (if (dbsnp && baseRecalibration) numberLibs else 0)
104

Peter van 't Hof's avatar
Peter van 't Hof committed
105
106
107
      pipeline.summarySettings.get("annotation") shouldBe Some(annotation)
      pipeline.summarySettings.get("sv_calling") shouldBe Some(svCalling)

Peter van 't Hof's avatar
Peter van 't Hof committed
108
109
110
111
112
113
114
115
116
117
      pipeline.samples foreach {
        case (sampleId, sample) =>
          sample.summarySettings.get("single_sample_variantcalling") shouldBe Some(sampleCalling)
          sample.summarySettings.get("use_indel_realigner") shouldBe Some(realign)
          sample.libraries.foreach {
            case (libId, lib) =>
              lib.summarySettings.get("library_variantcalling") shouldBe Some(libraryCalling)
              lib.summarySettings.get("use_indel_realigner") shouldBe Some(realign)
              lib.summarySettings.get("use_base_recalibration") shouldBe Some(baseRecalibration && dbsnp)
          }
Peter van 't Hof's avatar
Peter van 't Hof committed
118
119
      }

120
121
122
123
      pipeline.functions.count(_.isInstanceOf[VcfStats]) shouldBe (
        (if (multisampleCalling) 2 else 0) +
        (if (sampleCalling) numberSamples * 2 else 0) +
        (if (libraryCalling) numberLibs * 2 else 0))
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
class ShivaDefaultTest extends ShivaTestTrait
class ShivaNoDbsnpTest extends ShivaTestTrait {
  override def sample1 = Array(true)
  override def sample2 = Array(false)
  override def realignProvider = Array(true)
  override def dbsnp = false
}
class ShivaLibraryCallingTest extends ShivaTestTrait {
  override def sample1 = Array(true, false)
  override def sample2 = Array(false, true)
  override def realignProvider = Array(false)
  override def baseRecalibrationProvider = Array(false)
  override def libraryCalling = true
}
class ShivaSampleCallingTest extends ShivaTestTrait {
  override def sample1 = Array(true, false)
  override def sample2 = Array(false, true)
  override def realignProvider = Array(false)
  override def baseRecalibrationProvider = Array(false)
  override def sampleCalling = true
}
Peter van 't Hof's avatar
Peter van 't Hof committed
149
150
151
152
153
154
155
class ShivaWithSvCallingTest extends ShivaTestTrait {
  override def sample1 = Array(true)
  override def sample2 = Array(false)
  override def realignProvider = Array(false)
  override def baseRecalibrationProvider = Array(false)
  override def svCalling = true
}
156
157
158
159
160
161
162
class ShivaWithCnvCallingTest extends ShivaTestTrait {
  override def sample1 = Array(true)
  override def sample2 = Array(false)
  override def realignProvider = Array(false)
  override def baseRecalibrationProvider = Array(false)
  override def cnvCalling = true
}
Peter van 't Hof's avatar
Peter van 't Hof committed
163
164
165
166
167
168
169
class ShivaWithAnnotationTest extends ShivaTestTrait {
  override def sample1 = Array(true)
  override def sample2 = Array(false)
  override def realignProvider = Array(false)
  override def baseRecalibrationProvider = Array(false)
  override def annotation = true
}
170

171
172
object ShivaTest {
  val outputDir = Files.createTempDir()
Peter van 't Hof's avatar
Peter van 't Hof committed
173
174
175
176
177
178
  new File(outputDir, "input").mkdirs()
  def inputTouch(name: String): String = {
    val file = new File(outputDir, "input" + File.separator + name)
    Files.touch(file)
    file.getAbsolutePath
  }
179

180
181
182
183
184
185
186
187
188
189
190
  private def copyFile(name: String): Unit = {
    val is = getClass.getResourceAsStream("/" + name)
    val os = new FileOutputStream(new File(outputDir, name))
    org.apache.commons.io.IOUtils.copy(is, os)
    os.close()
  }

  copyFile("ref.fa")
  copyFile("ref.dict")
  copyFile("ref.fa.fai")

191
192
  val config = Map(
    "name_prefix" -> "test",
Peter van 't Hof's avatar
Peter van 't Hof committed
193
194
195
    "cache" -> true,
    "dir" -> "test",
    "vep_script" -> "test",
196
    "output_dir" -> outputDir,
Peter van 't Hof's avatar
Peter van 't Hof committed
197
    "reference_fasta" -> (outputDir + File.separator + "ref.fa"),
198
199
200
201
202
203
204
205
206
207
208
209
210
211
    "gatk_jar" -> "test",
    "samtools" -> Map("exe" -> "test"),
    "bcftools" -> Map("exe" -> "test"),
    "fastqc" -> Map("exe" -> "test"),
    "input_alleles" -> "test",
    "variantcallers" -> "raw",
    "fastqc" -> Map("exe" -> "test"),
    "seqtk" -> Map("exe" -> "test"),
    "sickle" -> Map("exe" -> "test"),
    "cutadapt" -> Map("exe" -> "test"),
    "bwa" -> Map("exe" -> "test"),
    "samtools" -> Map("exe" -> "test"),
    "macs2" -> Map("exe" -> "test"),
    "igvtools" -> Map("exe" -> "test"),
212
213
214
    "wigtobigwig" -> Map("exe" -> "test"),
    "md5sum" -> Map("exe" -> "test"),
    "bgzip" -> Map("exe" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
215
216
217
218
219
220
221
222
223
224
225
226
    "tabix" -> Map("exe" -> "test"),
    "breakdancerconfig" -> Map("exe" -> "test"),
    "breakdancercaller" -> Map("exe" -> "test"),
    "pindelconfig" -> Map("exe" -> "test"),
    "pindelcaller" -> Map("exe" -> "test"),
    "pindelvcf" -> Map("exe" -> "test"),
    "clever" -> Map("exe" -> "test"),
    "delly" -> Map("exe" -> "test"),
    "pysvtools" -> Map(
      "exe" -> "test",
      "exclusion_regions" -> "test",
      "translocations_only" -> false)
227
228
229
230
231
  )

  val sample1 = Map(
    "samples" -> Map("sample1" -> Map("libraries" -> Map(
      "lib1" -> Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
232
233
        "R1" -> inputTouch("1_1_R1.fq"),
        "R2" -> inputTouch("1_1_R2.fq")
234
235
236
237
238
      )
    )
    )))

  val sample2 = Map(
239
    "samples" -> Map("sample3" -> Map("libraries" -> Map(
240
      "lib1" -> Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
241
242
        "R1" -> inputTouch("2_1_R1.fq"),
        "R2" -> inputTouch("2_1_R2.fq")
243
244
      ),
      "lib2" -> Map(
245
246
        "R1" -> inputTouch("2_2_R1.fq"),
        "R2" -> inputTouch("2_2_R2.fq")
247
248
249
250
      )
    )
    )))
}