ShivaTest.scala 10.1 KB
Newer Older
bow's avatar
bow committed
1
/**
2
3
4
5
6
7
8
9
10
11
12
13
14
  * 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 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.
  */
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, IOException}
18

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

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

  @DataProvider(name = "shivaOptions")
pjvan_thof's avatar
pjvan_thof committed
52
  def shivaOptions: Array[Array[Any]] = {
53
54
55
    for (s1 <- sample1; s2 <- sample2;
         realign <- realignProvider; baseRecalibration <- baseRecalibrationProvider)
      yield Array("", s1, s2, realign, baseRecalibration)
56
57
  }

58
59
60
61
62
  def sample1 = Array(false, true)
  def sample2 = Array(false, true)
  def realignProvider = Array(false, true)
  def baseRecalibrationProvider = Array(false, true)
  def multisampleCalling: Boolean = true
Peter van 't Hof's avatar
Peter van 't Hof committed
63
64
  def dbsnp = true
  def svCalling = false
65
  def cnvCalling = false
Peter van 't Hof's avatar
Peter van 't Hof committed
66
  def annotation = false
Peter van 't Hof's avatar
Peter van 't Hof committed
67
  def usePrintReads = true
68

Peter van 't Hof's avatar
Peter van 't Hof committed
69
70
  private var dirs: List[File] = Nil

71
  @Test(dataProvider = "shivaOptions")
72
73
74
75
76
  def testShiva(f: String,
                sample1: Boolean,
                sample2: Boolean,
                realign: Boolean,
                baseRecalibration: Boolean): Unit = {
Peter van 't Hof's avatar
Peter van 't Hof committed
77
    val outputDir = ShivaTest.outputDir
Peter van 't Hof's avatar
Peter van 't Hof committed
78
    dirs :+= outputDir
79
    val map = {
Peter van 't Hof's avatar
Peter van 't Hof committed
80
      var m: Map[String, Any] = ShivaTest.config(outputDir)
Peter van 't Hof's avatar
Peter van 't Hof committed
81
82
      if (sample1) m = ConfigUtils.mergeMaps(ShivaTest.sample1, m)
      if (sample2) m = ConfigUtils.mergeMaps(ShivaTest.sample2, m)
Peter van 't Hof's avatar
Peter van 't Hof committed
83
      if (dbsnp) m = ConfigUtils.mergeMaps(Map("dbsnp_vcf" -> "test.vcf.gz"), m)
84
85
86
87
88
89
90
91
92
93
94
95
      ConfigUtils.mergeMaps(
        Map(
          "multisample_variantcalling" -> multisampleCalling,
          "use_indel_realigner" -> realign,
          "use_base_recalibration" -> baseRecalibration,
          "sv_calling" -> svCalling,
          "cnv_calling" -> cnvCalling,
          "annotation" -> annotation,
          "use_printreads" -> usePrintReads
        ),
        m
      )
96

97
98
    }

99
    if (!sample1 && !sample2) { // When no samples
100
101
102
      intercept[IllegalArgumentException] {
        initPipeline(map).script()
      }
Peter van 't Hof's avatar
Peter van 't Hof committed
103
      Logging.errors.clear()
104
105
106
107
    } else {
      val pipeline = initPipeline(map)
      pipeline.script()

108
109
110
      val numberLibs = (if (sample1) 1 else 0) + (if (sample2) 2 else 0)
      val numberSamples = (if (sample1) 1 else 0) + (if (sample2) 1 else 0)

Peter van 't Hof's avatar
Peter van 't Hof committed
111
      pipeline.functions.count(_.isInstanceOf[MarkDuplicates]) shouldBe (numberLibs + numberSamples)
112

113
      // Gatk preprocess
114
115
116
117
118
119
      pipeline.functions.count(_.isInstanceOf[IndelRealigner]) shouldBe (if (realign) numberSamples
                                                                         else 0)
      pipeline.functions.count(_.isInstanceOf[RealignerTargetCreator]) shouldBe (if (realign)
                                                                                   numberSamples
                                                                                 else 0)
      pipeline.functions.count(_.isInstanceOf[BaseRecalibrator]) shouldBe (if (dbsnp && baseRecalibration)
pjvan_thof's avatar
pjvan_thof committed
120
                                                                             numberLibs * 2
121
122
123
124
                                                                           else 0)
      pipeline.functions.count(_.isInstanceOf[PrintReads]) shouldBe (if (dbsnp && baseRecalibration && usePrintReads)
                                                                       numberLibs
                                                                     else 0)
125

Peter van 't Hof's avatar
Peter van 't Hof committed
126
127
      pipeline.summarySettings.get("annotation") shouldBe Some(annotation)
      pipeline.summarySettings.get("sv_calling") shouldBe Some(svCalling)
128
      pipeline.summarySettings.get("cnv_calling") shouldBe Some(cnvCalling)
Peter van 't Hof's avatar
Peter van 't Hof committed
129

Peter van 't Hof's avatar
Peter van 't Hof committed
130
      pipeline.samples foreach {
Peter van 't Hof's avatar
Peter van 't Hof committed
131
        case (_, sample) =>
Peter van 't Hof's avatar
Peter van 't Hof committed
132
          sample.summarySettings.get("single_sample_variantcalling") shouldBe None
Peter van 't Hof's avatar
Peter van 't Hof committed
133
134
          sample.summarySettings.get("use_indel_realigner") shouldBe Some(realign)
          sample.libraries.foreach {
Peter van 't Hof's avatar
Peter van 't Hof committed
135
            case (_, lib) =>
Peter van 't Hof's avatar
Peter van 't Hof committed
136
              lib.summarySettings.get("library_variantcalling") shouldBe None
Peter van 't Hof's avatar
Peter van 't Hof committed
137
              lib.summarySettings.get("use_indel_realigner") shouldBe None // Should not exist anymore
138
139
              lib.summarySettings.get("use_base_recalibration") shouldBe Some(
                baseRecalibration && dbsnp)
Peter van 't Hof's avatar
Peter van 't Hof committed
140
          }
Peter van 't Hof's avatar
Peter van 't Hof committed
141
142
      }

pjvan_thof's avatar
pjvan_thof committed
143
144
      pipeline.functions.count(_.isInstanceOf[VcfStats]) shouldBe (if (multisampleCalling) 2
                                                                   else 0)
145
146
    }
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
147
148

  // remove temporary run directory all tests in the class have been run
Peter van 't Hof's avatar
Peter van 't Hof committed
149
  @AfterClass def removeTempOutputDir(): Unit = {
Peter van 't Hof's avatar
Peter van 't Hof committed
150
151
152
153
    dirs.filter(_.exists()).foreach { dir =>
      try {
        FileUtils.deleteDirectory(dir)
      } catch {
Peter van 't Hof's avatar
Peter van 't Hof committed
154
155
        case e: IOException if e.getMessage.startsWith("Unable to delete directory") =>
          Logging.logger.error(e.getMessage)
Peter van 't Hof's avatar
Peter van 't Hof committed
156
157
      }
    }
Peter van 't Hof's avatar
Peter van 't Hof committed
158
  }
159
160
}

161
162
class ShivaDefaultTest extends ShivaTestTrait
class ShivaNoDbsnpTest extends ShivaTestTrait {
Peter van 't Hof's avatar
Peter van 't Hof committed
163
164
165
  override def sample1 = Array(x = true)
  override def sample2 = Array(x = false)
  override def realignProvider = Array(x = true)
166
167
  override def dbsnp = false
}
Peter van 't Hof's avatar
Peter van 't Hof committed
168
class ShivaNoPrintReadsTest extends ShivaTestTrait {
Peter van 't Hof's avatar
Peter van 't Hof committed
169
  override def realignProvider = Array(x = false)
Peter van 't Hof's avatar
Peter van 't Hof committed
170
171
  override def usePrintReads = false
}
Peter van 't Hof's avatar
Peter van 't Hof committed
172
class ShivaWithSvCallingTest extends ShivaTestTrait {
Peter van 't Hof's avatar
Peter van 't Hof committed
173
174
175
176
  override def sample1 = Array(x = true)
  override def sample2 = Array(x = false)
  override def realignProvider = Array(x = false)
  override def baseRecalibrationProvider = Array(x = false)
Peter van 't Hof's avatar
Peter van 't Hof committed
177
178
  override def svCalling = true
}
179
class ShivaWithCnvCallingTest extends ShivaTestTrait {
Peter van 't Hof's avatar
Peter van 't Hof committed
180
181
182
183
  override def sample1 = Array(x = true)
  override def sample2 = Array(x = false)
  override def realignProvider = Array(x = false)
  override def baseRecalibrationProvider = Array(x = false)
184
185
  override def cnvCalling = true
}
Peter van 't Hof's avatar
Peter van 't Hof committed
186
class ShivaWithAnnotationTest extends ShivaTestTrait {
Peter van 't Hof's avatar
Peter van 't Hof committed
187
188
189
190
  override def sample1 = Array(x = true)
  override def sample2 = Array(x = false)
  override def realignProvider = Array(x = false)
  override def baseRecalibrationProvider = Array(x = false)
Peter van 't Hof's avatar
Peter van 't Hof committed
191
192
  override def annotation = true
}
193

194
object ShivaTest {
pjvan_thof's avatar
pjvan_thof committed
195
  def outputDir: File = Files.createTempDir()
Peter van 't Hof's avatar
Peter van 't Hof committed
196

pjvan_thof's avatar
pjvan_thof committed
197
  val inputDir: File = Files.createTempDir()
Peter van 't Hof's avatar
Peter van 't Hof committed
198

Peter van 't Hof's avatar
Peter van 't Hof committed
199
  def inputTouch(name: String): String = {
Peter van 't Hof's avatar
Peter van 't Hof committed
200
    val file = new File(inputDir, name)
Peter van 't Hof's avatar
Peter van 't Hof committed
201
    Files.touch(file)
Peter van 't Hof's avatar
Peter van 't Hof committed
202
    file.deleteOnExit()
Peter van 't Hof's avatar
Peter van 't Hof committed
203
204
    file.getAbsolutePath
  }
205

206
207
  private def copyFile(name: String): Unit = {
    val is = getClass.getResourceAsStream("/" + name)
Peter van 't Hof's avatar
Peter van 't Hof committed
208
    val os = new FileOutputStream(new File(inputDir, name))
209
210
211
212
213
214
215
216
    org.apache.commons.io.IOUtils.copy(is, os)
    os.close()
  }

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

Peter van 't Hof's avatar
Peter van 't Hof committed
217
  def config(outputDir: File) = Map(
218
    "skip_write_dependencies" -> true,
219
    "name_prefix" -> "test",
Peter van 't Hof's avatar
Peter van 't Hof committed
220
221
222
    "cache" -> true,
    "dir" -> "test",
    "vep_script" -> "test",
223
    "output_dir" -> outputDir,
Peter van 't Hof's avatar
Peter van 't Hof committed
224
    "reference_fasta" -> (inputDir + File.separator + "ref.fa"),
225
226
227
228
229
230
    "gatk_jar" -> "test",
    "samtools" -> Map("exe" -> "test"),
    "bcftools" -> Map("exe" -> "test"),
    "fastqc" -> Map("exe" -> "test"),
    "input_alleles" -> "test",
    "variantcallers" -> "raw",
Peter van 't Hof's avatar
Typo    
Peter van 't Hof committed
231
    "rscript" -> Map("exe" -> "test"),
232
233
234
235
236
237
238
    "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"),
239
    "igvtools" -> Map("exe" -> "test", "igvtools_jar" -> "test"),
240
241
242
    "wigtobigwig" -> Map("exe" -> "test"),
    "md5sum" -> Map("exe" -> "test"),
    "bgzip" -> Map("exe" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
243
244
245
246
247
248
249
250
    "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"),
251
252
253
    "pysvtools" -> Map("exe" -> "test",
                       "exclusion_regions" -> "test",
                       "translocations_only" -> false),
Peter van 't Hof's avatar
Peter van 't Hof committed
254
255
256
257
258
    "freec" -> Map(
      "exe" -> "test",
      "chrFiles" -> "test",
      "chrLenFile" -> "test"
    )
259
260
261
  )

  val sample1 = Map(
262
263
264
265
266
267
268
269
    "samples" -> Map(
      "sample1" -> Map(
        "libraries" -> Map(
          "lib1" -> Map(
            "R1" -> inputTouch("1_1_R1.fq"),
            "R2" -> inputTouch("1_1_R2.fq")
          )
        ))))
270
271

  val sample2 = Map(
272
273
274
275
276
277
278
279
280
281
282
283
284
    "samples" -> Map(
      "sample3" -> Map(
        "libraries" -> Map(
          "lib1" -> Map(
            "R1" -> inputTouch("2_1_R1.fq"),
            "R2" -> inputTouch("2_1_R2.fq")
          ),
          "lib2" -> Map(
            "R1" -> inputTouch("2_2_R1.fq"),
            "R2" -> inputTouch("2_2_R2.fq")
          )
        ))))
}