ShivaVariantcallingTest.scala 6.48 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 56 57
      bcftools_singlesample <- bool;
      haplotypeCallerGvcf <- bool;
      haplotypeCallerAllele <- bool;
      unifiedGenotyperAllele <- bool;
      unifiedGenotyper <- bool;
      haplotypeCaller <- bool
    ) yield Array[Any](bams, raw, bcftools, bcftools_singlesample, unifiedGenotyper, haplotypeCaller, haplotypeCallerGvcf, haplotypeCallerAllele, unifiedGenotyperAllele)
Peter van 't Hof's avatar
Peter van 't Hof committed
58
    ).toArray
Peter van 't Hof's avatar
Peter van 't Hof committed
59 60 61
  }

  @Test(dataProvider = "shivaVariantcallingOptions")
62 63 64
  def testShivaVariantcalling(bams: Int,
                              raw: Boolean,
                              bcftools: Boolean,
65
                              bcftoolsSinglesample: Boolean,
66 67 68 69 70
                              unifiedGenotyper: Boolean,
                              haplotypeCaller: Boolean,
                              haplotypeCallerGvcf: Boolean,
                              haplotypeCallerAllele: Boolean,
                              unifiedGenotyperAllele: Boolean,
Peter van 't Hof's avatar
Peter van 't Hof committed
71
                              freebayes: Boolean,
72
                              varscanCnsSinglesample: Boolean) = {
Peter van 't Hof's avatar
Peter van 't Hof committed
73 74 75
    val callers: ListBuffer[String] = ListBuffer()
    if (raw) callers.append("raw")
    if (bcftools) callers.append("bcftools")
76
    if (bcftoolsSinglesample) callers.append("bcftools_singlesample")
77 78 79 80 81
    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
82
    if (freebayes) callers.append("freebayes")
83
    if (varscanCnsSinglesample) callers.append("varscan_cns_singlesample")
84
    val map = Map("variantcallers" -> callers.toList)
Peter van 't Hof's avatar
Peter van 't Hof committed
85 86
    val pipeline = initPipeline(map)

87
    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
88

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

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

    if (!illegalArgumentException) {
      pipeline.script()

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

Peter van 't Hof's avatar
Peter van 't Hof committed
104
      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
105 106 107 108
      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
109
      pipeline.functions.count(_.isInstanceOf[VcfFilter]) shouldBe (if (raw) bams else 0)
110 111 112 113 114
      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
115 116
    }
  }
117 118 119 120

  @AfterClass def removeTempOutputDir() = {
    FileUtils.deleteDirectory(ShivaVariantcallingTest.outputDir)
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
121 122 123 124
}

object ShivaVariantcallingTest {
  val outputDir = Files.createTempDir()
Peter van 't Hof's avatar
Peter van 't Hof committed
125
  new File(outputDir, "input").mkdirs()
126
  def inputTouch(name: String): File = {
Peter van 't Hof's avatar
Peter van 't Hof committed
127 128 129 130
    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
131

Peter van 't Hof's avatar
Peter van 't Hof committed
132 133 134 135 136 137 138 139 140 141 142
  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
143 144 145
  val config = Map(
    "name_prefix" -> "test",
    "output_dir" -> outputDir,
Peter van 't Hof's avatar
Peter van 't Hof committed
146 147 148
    "cache" -> true,
    "dir" -> "test",
    "vep_script" -> "test",
Peter van 't Hof's avatar
Peter van 't Hof committed
149
    "reference_fasta" -> (outputDir + File.separator + "ref.fa"),
Peter van 't Hof's avatar
Peter van 't Hof committed
150 151
    "gatk_jar" -> "test",
    "samtools" -> Map("exe" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
152
    "bcftools" -> Map("exe" -> "test"),
153 154
    "md5sum" -> Map("exe" -> "test"),
    "bgzip" -> Map("exe" -> "test"),
155
    "tabix" -> Map("exe" -> "test"),
156
    "input_alleles" -> "test.vcf.gz"
Peter van 't Hof's avatar
Peter van 't Hof committed
157 158
  )
}