ShivaVariantcallingTest.scala 6.59 KB
Newer Older
bow's avatar
bow committed
1
/**
Peter van 't Hof's avatar
Peter van 't Hof committed
2 3 4 5
 * Due to the license issue with GATK, this part of Biopet can only be used inside the
 * LUMC. Please refer to https://git.lumc.nl/biopet/biopet/wikis/home for instructions
 * on how to use this protected part of biopet or contact us at sasc@lumc.nl
 */
Peter van 't Hof's avatar
Peter van 't Hof committed
6
package nl.lumc.sasc.biopet.pipelines.shiva
7

Peter van 't Hof's avatar
Peter van 't Hof committed
8
import java.io.{ File, FileOutputStream }
Peter van 't Hof's avatar
Peter van 't Hof committed
9 10

import com.google.common.io.Files
Peter van 't Hof's avatar
Peter van 't Hof committed
11
import nl.lumc.sasc.biopet.core.BiopetPipe
Peter van 't Hof's avatar
Peter van 't Hof committed
12
import nl.lumc.sasc.biopet.extensions.Freebayes
Peter van 't Hof's avatar
Peter van 't Hof committed
13
import nl.lumc.sasc.biopet.extensions.bcftools.{ BcftoolsCall, BcftoolsMerge }
14
import nl.lumc.sasc.biopet.utils.config.Config
Peter van 't Hof's avatar
Peter van 't Hof committed
15
import nl.lumc.sasc.biopet.extensions.gatk.CombineVariants
Peter van 't Hof's avatar
Peter van 't Hof committed
16 17
import nl.lumc.sasc.biopet.extensions.gatk.broad.{ HaplotypeCaller, UnifiedGenotyper }
import nl.lumc.sasc.biopet.extensions.tools.{ MpileupToVcf, VcfFilter, VcfStats }
18
import nl.lumc.sasc.biopet.pipelines.shiva.ShivaVariantcalling
Peter van 't Hof's avatar
Peter van 't Hof committed
19
import nl.lumc.sasc.biopet.utils.ConfigUtils
20
import org.apache.commons.io.FileUtils
Peter van 't Hof's avatar
Peter van 't Hof committed
21 22 23
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
24
import org.testng.annotations.{ AfterClass, DataProvider, Test }
Peter van 't Hof's avatar
Peter van 't Hof committed
25

Peter van 't Hof's avatar
Peter van 't Hof committed
26 27
import scala.collection.mutable.ListBuffer

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

  @DataProvider(name = "shivaVariantcallingOptions")
Peter van 't Hof's avatar
Peter van 't Hof committed
44
  def shivaVariantcallingOptions = {
45
    val bool = Array(true, false)
46

Sander van der Zeeuw's avatar
Sander van der Zeeuw committed
47
    (for (
48
      bams <- 0 to 2;
Sander van der Zeeuw's avatar
Sander van der Zeeuw committed
49 50
      raw <- bool;
      bcftools <- bool;
51 52 53 54 55
      bcftools_singlesample <- bool;
      haplotypeCallerGvcf <- bool;
      haplotypeCallerAllele <- bool;
      unifiedGenotyperAllele <- bool;
      unifiedGenotyper <- bool;
Peter van 't Hof's avatar
Peter van 't Hof committed
56 57 58 59 60 61
      haplotypeCaller <- bool;
      freebayes <- bool;
      varscanCnsSinglesample <- bool
    ) yield Array[Any](bams, raw, bcftools, bcftools_singlesample, unifiedGenotyper,
      haplotypeCaller, haplotypeCallerGvcf, haplotypeCallerAllele, unifiedGenotyperAllele,
      freebayes, varscanCnsSinglesample)
Peter van 't Hof's avatar
Peter van 't Hof committed
62
    ).toArray
Peter van 't Hof's avatar
Peter van 't Hof committed
63 64 65
  }

  @Test(dataProvider = "shivaVariantcallingOptions")
66 67 68
  def testShivaVariantcalling(bams: Int,
                              raw: Boolean,
                              bcftools: Boolean,
69
                              bcftoolsSinglesample: Boolean,
70 71 72 73 74
                              unifiedGenotyper: Boolean,
                              haplotypeCaller: Boolean,
                              haplotypeCallerGvcf: Boolean,
                              haplotypeCallerAllele: Boolean,
                              unifiedGenotyperAllele: Boolean,
Peter van 't Hof's avatar
Peter van 't Hof committed
75
                              freebayes: Boolean,
76
                              varscanCnsSinglesample: Boolean) = {
Peter van 't Hof's avatar
Peter van 't Hof committed
77 78 79
    val callers: ListBuffer[String] = ListBuffer()
    if (raw) callers.append("raw")
    if (bcftools) callers.append("bcftools")
80
    if (bcftoolsSinglesample) callers.append("bcftools_singlesample")
81 82 83 84 85
    if (unifiedGenotyper) callers.append("unifiedgenotyper")
    if (haplotypeCallerGvcf) callers.append("haplotypecaller_gvcf")
    if (haplotypeCallerAllele) callers.append("haplotypecaller_allele")
    if (unifiedGenotyperAllele) callers.append("unifiedgenotyper_allele")
    if (haplotypeCaller) callers.append("haplotypecaller")
Peter van 't Hof's avatar
Peter van 't Hof committed
86
    if (freebayes) callers.append("freebayes")
87
    if (varscanCnsSinglesample) callers.append("varscan_cns_singlesample")
88
    val map = Map("variantcallers" -> callers.toList)
Peter van 't Hof's avatar
Peter van 't Hof committed
89 90
    val pipeline = initPipeline(map)

91
    pipeline.inputBams = (for (n <- 1 to bams) yield n.toString -> ShivaVariantcallingTest.inputTouch("bam_" + n + ".bam")).toMap
Peter van 't Hof's avatar
Peter van 't Hof committed
92

93 94 95 96 97
    val illegalArgumentException = pipeline.inputBams.isEmpty ||
      (!raw && !bcftools &&
        !haplotypeCaller && !unifiedGenotyper &&
        !haplotypeCallerGvcf && !haplotypeCallerAllele && !unifiedGenotyperAllele &&
        !bcftoolsSinglesample)
Peter van 't Hof's avatar
Peter van 't Hof committed
98 99 100 101 102 103 104 105

    if (illegalArgumentException) intercept[IllegalArgumentException] {
      pipeline.script()
    }

    if (!illegalArgumentException) {
      pipeline.script()

Peter van 't Hof's avatar
Peter van 't Hof committed
106 107
      val pipesJobs = pipeline.functions.filter(_.isInstanceOf[BiopetPipe]).flatMap(_.asInstanceOf[BiopetPipe].pipesJobs)

Peter van 't Hof's avatar
Peter van 't Hof committed
108
      pipeline.functions.count(_.isInstanceOf[CombineVariants]) shouldBe (1 + (if (raw) 1 else 0) + (if (varscanCnsSinglesample) 1 else 0))
Peter van 't Hof's avatar
Peter van 't Hof committed
109 110 111 112
      pipesJobs.count(_.isInstanceOf[BcftoolsCall]) shouldBe (if (bcftools) 1 else 0) + (if (bcftoolsSinglesample) bams else 0)
      pipeline.functions.count(_.isInstanceOf[BcftoolsMerge]) shouldBe (if (bcftoolsSinglesample && bams > 1) 1 else 0)
      pipesJobs.count(_.isInstanceOf[Freebayes]) shouldBe (if (freebayes) 1 else 0)
      pipesJobs.count(_.isInstanceOf[MpileupToVcf]) shouldBe (if (raw) bams else 0)
Peter van 't Hof's avatar
Peter van 't Hof committed
113
      pipeline.functions.count(_.isInstanceOf[VcfFilter]) shouldBe (if (raw) bams else 0)
114 115 116 117 118
      pipeline.functions.count(_.isInstanceOf[HaplotypeCaller]) shouldBe (if (haplotypeCaller) 1 else 0) +
        (if (haplotypeCallerAllele) 1 else 0) + (if (haplotypeCallerGvcf) bams else 0)
      pipeline.functions.count(_.isInstanceOf[UnifiedGenotyper]) shouldBe (if (unifiedGenotyper) 1 else 0) +
        (if (unifiedGenotyperAllele) 1 else 0)
      pipeline.functions.count(_.isInstanceOf[VcfStats]) shouldBe (1 + callers.size)
Peter van 't Hof's avatar
Peter van 't Hof committed
119 120
    }
  }
121 122 123 124

  @AfterClass def removeTempOutputDir() = {
    FileUtils.deleteDirectory(ShivaVariantcallingTest.outputDir)
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
125 126 127 128
}

object ShivaVariantcallingTest {
  val outputDir = Files.createTempDir()
Peter van 't Hof's avatar
Peter van 't Hof committed
129
  new File(outputDir, "input").mkdirs()
130
  def inputTouch(name: String): File = {
Peter van 't Hof's avatar
Peter van 't Hof committed
131 132 133 134
    val file = new File(outputDir, "input" + File.separator + name).getAbsoluteFile
    Files.touch(file)
    file
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
135

Peter van 't Hof's avatar
Peter van 't Hof committed
136 137 138 139 140 141 142 143 144 145 146
  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")

Peter van 't Hof's avatar
Peter van 't Hof committed
147 148 149
  val config = Map(
    "name_prefix" -> "test",
    "output_dir" -> outputDir,
Peter van 't Hof's avatar
Peter van 't Hof committed
150 151 152
    "cache" -> true,
    "dir" -> "test",
    "vep_script" -> "test",
Peter van 't Hof's avatar
Peter van 't Hof committed
153
    "reference_fasta" -> (outputDir + File.separator + "ref.fa"),
Peter van 't Hof's avatar
Peter van 't Hof committed
154 155
    "gatk_jar" -> "test",
    "samtools" -> Map("exe" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
156
    "bcftools" -> Map("exe" -> "test"),
157 158
    "md5sum" -> Map("exe" -> "test"),
    "bgzip" -> Map("exe" -> "test"),
159
    "tabix" -> Map("exe" -> "test"),
160
    "input_alleles" -> "test.vcf.gz"
Peter van 't Hof's avatar
Peter van 't Hof committed
161 162
  )
}