ShivaTest.scala 9.32 KB
Newer Older
bow's avatar
bow committed
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
bow's avatar
bow committed
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
16
package nl.lumc.sasc.biopet.pipelines.shiva

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

19
import com.google.common.io.Files
Peter van 't Hof's avatar
Peter van 't Hof committed
20
import nl.lumc.sasc.biopet.extensions.gatk.{ BaseRecalibrator, IndelRealigner, PrintReads, RealignerTargetCreator }
21
import nl.lumc.sasc.biopet.extensions.picard.MarkDuplicates
22
import nl.lumc.sasc.biopet.extensions.tools.VcfStats
Peter van 't Hof's avatar
Peter van 't Hof committed
23
import nl.lumc.sasc.biopet.utils.{ ConfigUtils, Logging }
24
import nl.lumc.sasc.biopet.utils.config.Config
25
26
27
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
28
import org.testng.annotations.{ DataProvider, Test }
29
30

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

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

53
54
55
56
57
58
59
  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
60
61
  def dbsnp = true
  def svCalling = false
62
  def cnvCalling = false
Peter van 't Hof's avatar
Peter van 't Hof committed
63
  def annotation = false
64

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

83
84
    }

85
    if (!sample1 && !sample2) { // When no samples
86
87
88
      intercept[IllegalArgumentException] {
        initPipeline(map).script()
      }
Peter van 't Hof's avatar
Peter van 't Hof committed
89
      Logging.errors.clear()
90
91
92
93
    } 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
      pipeline.summarySettings.get("annotation") shouldBe Some(annotation)
      pipeline.summarySettings.get("sv_calling") shouldBe Some(svCalling)
107
      pipeline.summarySettings.get("cnv_calling") shouldBe Some(cnvCalling)
Peter van 't Hof's avatar
Peter van 't Hof committed
108

Peter van 't Hof's avatar
Peter van 't Hof committed
109
110
111
112
113
114
115
116
117
118
      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
119
120
      }

121
122
123
124
      pipeline.functions.count(_.isInstanceOf[VcfStats]) shouldBe (
        (if (multisampleCalling) 2 else 0) +
        (if (sampleCalling) numberSamples * 2 else 0) +
        (if (libraryCalling) numberLibs * 2 else 0))
125
126
127
128
    }
  }
}

129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
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
150
151
152
153
154
155
156
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
}
157
158
159
160
161
162
163
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
164
165
166
167
168
169
170
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
}
171

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

182
183
184
185
186
187
188
189
190
191
192
  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")

193
194
  val config = Map(
    "name_prefix" -> "test",
Peter van 't Hof's avatar
Peter van 't Hof committed
195
196
197
    "cache" -> true,
    "dir" -> "test",
    "vep_script" -> "test",
198
    "output_dir" -> outputDir,
Peter van 't Hof's avatar
Peter van 't Hof committed
199
    "reference_fasta" -> (outputDir + File.separator + "ref.fa"),
200
201
202
203
204
205
206
207
208
209
210
211
212
    "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"),
213
    "igvtools" -> Map("exe" -> "test", "igvtools_jar" -> "test"),
214
215
216
    "wigtobigwig" -> Map("exe" -> "test"),
    "md5sum" -> Map("exe" -> "test"),
    "bgzip" -> Map("exe" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
217
218
219
220
221
222
223
224
225
226
227
    "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",
Peter van 't Hof's avatar
Peter van 't Hof committed
228
229
230
231
232
233
      "translocations_only" -> false),
    "freec" -> Map(
      "exe" -> "test",
      "chrFiles" -> "test",
      "chrLenFile" -> "test"
    )
234
235
236
237
238
  )

  val sample1 = Map(
    "samples" -> Map("sample1" -> Map("libraries" -> Map(
      "lib1" -> Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
239
240
        "R1" -> inputTouch("1_1_R1.fq"),
        "R2" -> inputTouch("1_1_R2.fq")
241
242
243
244
245
      )
    )
    )))

  val sample2 = Map(
246
    "samples" -> Map("sample3" -> Map("libraries" -> Map(
247
      "lib1" -> Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
248
249
        "R1" -> inputTouch("2_1_R1.fq"),
        "R2" -> inputTouch("2_1_R2.fq")
250
251
      ),
      "lib2" -> Map(
252
253
        "R1" -> inputTouch("2_2_R1.fq"),
        "R2" -> inputTouch("2_2_R2.fq")
254
255
256
257
      )
    )
    )))
}