MultisampleMappingTest.scala 9.89 KB
Newer Older
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.
 */
Peter van 't Hof's avatar
Peter van 't Hof committed
15 16
package nl.lumc.sasc.biopet.pipelines.mapping

Peter van 't Hof's avatar
Peter van 't Hof committed
17
import java.io.{ File, FileOutputStream }
Peter van 't Hof's avatar
Peter van 't Hof committed
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.core.BiopetCommandLineFunction
Peter van 't Hof's avatar
Peter van 't Hof committed
21
import nl.lumc.sasc.biopet.extensions.centrifuge.Centrifuge
Peter van 't Hof's avatar
Peter van 't Hof committed
22
import nl.lumc.sasc.biopet.extensions.picard.{ MarkDuplicates, MergeSamFiles }
Peter van 't Hof's avatar
Peter van 't Hof committed
23
import nl.lumc.sasc.biopet.extensions.sambamba.SambambaMarkdup
Peter van 't Hof's avatar
Peter van 't Hof committed
24
import nl.lumc.sasc.biopet.utils.{ ConfigUtils, Logging }
Peter van 't Hof's avatar
Peter van 't Hof committed
25
import nl.lumc.sasc.biopet.utils.config.Config
Peter van 't Hof's avatar
Peter van 't Hof committed
26
import org.apache.commons.io.FileUtils
Peter van 't Hof's avatar
Peter van 't Hof committed
27 28 29
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
30
import org.testng.annotations.{ AfterClass, DataProvider, Test }
Peter van 't Hof's avatar
Peter van 't Hof committed
31 32

/**
33 34
 * Created by pjvanthof on 15/05/16.
 */
Peter van 't Hof's avatar
Peter van 't Hof committed
35
trait MultisampleMappingTestTrait extends TestNGSuite with Matchers {
Peter van 't Hof's avatar
Peter van 't Hof committed
36 37 38 39

  val outputDir = MultisampleMappingTestTrait.outputDir
  val configMap = MultisampleMappingTestTrait.config(outputDir)

Peter van 't Hof's avatar
Peter van 't Hof committed
40 41 42
  def initPipeline(map: Map[String, Any]): MultisampleMapping = {
    new MultisampleMapping() {
      override def configNamespace = "multisamplemapping"
Peter van 't Hof's avatar
Peter van 't Hof committed
43
      override def globalConfig = new Config(ConfigUtils.mergeMaps(map, configMap))
Peter van 't Hof's avatar
Peter van 't Hof committed
44 45 46 47 48 49 50 51
      qSettings = new QSettings
      qSettings.runName = "test"
    }
  }

  def mergeStrategies = MultisampleMapping.MergeStrategy.values
  def bamToFastq = false
  def correctReadgroups = false
Peter van 't Hof's avatar
Peter van 't Hof committed
52
  def unmappedToGears = false
Peter van 't Hof's avatar
Peter van 't Hof committed
53 54
  def sample1 = Array(true, false)
  def sample2 = Array(true, false)
55 56
  def sample3 = false
  def sample4 = false
Peter van 't Hof's avatar
Peter van 't Hof committed
57 58 59 60

  @DataProvider(name = "mappingOptions")
  def mappingOptions = {
    for (
61
      merge <- mergeStrategies.toArray; s1 <- sample1; s2 <- sample2
Peter van 't Hof's avatar
Peter van 't Hof committed
62
    ) yield Array(merge, s1, s2)
Peter van 't Hof's avatar
Peter van 't Hof committed
63 64 65 66 67
  }

  @Test(dataProvider = "mappingOptions")
  def testMultisampleMapping(merge: MultisampleMapping.MergeStrategy.Value, sample1: Boolean, sample2: Boolean): Unit = {
    val map: Map[String, Any] = {
Peter van 't Hof's avatar
Peter van 't Hof committed
68
      var m: Map[String, Any] = configMap
Peter van 't Hof's avatar
Peter van 't Hof committed
69 70
      if (sample1) m = ConfigUtils.mergeMaps(MultisampleMappingTestTrait.sample1, m)
      if (sample2) m = ConfigUtils.mergeMaps(MultisampleMappingTestTrait.sample2, m)
71 72
      if (sample3) m = ConfigUtils.mergeMaps(MultisampleMappingTestTrait.sample3, m)
      if (sample4) m = ConfigUtils.mergeMaps(MultisampleMappingTestTrait.sample4, m)
Peter van 't Hof's avatar
Peter van 't Hof committed
73 74 75 76 77 78
      m ++ Map(
        "merge_strategy" -> merge.toString,
        "bam_to_fastq" -> bamToFastq,
        "correct_readgroups" -> correctReadgroups,
        "unmapped_to_gears" -> unmappedToGears
      )
Peter van 't Hof's avatar
Peter van 't Hof committed
79 80
    }

81
    if (!sample1 && !sample2 && !sample3 && !sample4) { // When no samples
Peter van 't Hof's avatar
Peter van 't Hof committed
82
      intercept[IllegalStateException] {
Peter van 't Hof's avatar
Peter van 't Hof committed
83 84
        initPipeline(map).script()
      }
Peter van 't Hof's avatar
Peter van 't Hof committed
85
      Logging.errors.clear()
86 87 88 89
    } else if (sample4 && !bamToFastq && !correctReadgroups) {
      intercept[IllegalStateException] {
        initPipeline(map).script()
      }
Peter van 't Hof's avatar
Peter van 't Hof committed
90 91 92 93
    } else {
      val pipeline = initPipeline(map)
      pipeline.script()

94
      val numberFastqLibs = (if (sample1) 1 else 0) + (if (sample2) 2 else 0) + (if (sample3 && bamToFastq) 1 else 0) + (if (sample4 && bamToFastq) 1 else 0)
Peter van 't Hof's avatar
Peter van 't Hof committed
95
      val numberSamples = (if (sample1) 1 else 0) + (if (sample2) 1 else 0) + (if (sample3) 1 else 0) + (if (sample4) 1 else 0)
Peter van 't Hof's avatar
Peter van 't Hof committed
96

Peter van 't Hof's avatar
Peter van 't Hof committed
97 98
      val pipesJobs = pipeline.functions.filter(_.isInstanceOf[BiopetCommandLineFunction])
        .flatMap(_.asInstanceOf[BiopetCommandLineFunction].pipesJobs)
99

Peter van 't Hof's avatar
Peter van 't Hof committed
100 101 102 103
      if (merge == MultisampleMapping.MergeStrategy.PreProcessMarkDuplicates) {
        ""
      }

Peter van 't Hof's avatar
Peter van 't Hof committed
104
      import MultisampleMapping.MergeStrategy
105
      pipeline.functions.count(_.isInstanceOf[MarkDuplicates]) shouldBe (numberFastqLibs +
Peter van 't Hof's avatar
Peter van 't Hof committed
106
        (if (merge == MergeStrategy.MarkDuplicates || merge == MergeStrategy.PreProcessMarkDuplicates) numberSamples else 0))
Peter van 't Hof's avatar
Peter van 't Hof committed
107 108
      pipeline.functions.count(_.isInstanceOf[MergeSamFiles]) shouldBe (
        (if (sample2 && (merge == MergeStrategy.MergeSam || merge == MergeStrategy.PreProcessMergeSam)) 1 else 0))
Peter van 't Hof's avatar
Peter van 't Hof committed
109 110
      pipeline.functions.count(_.isInstanceOf[SambambaMarkdup]) shouldBe
        (if (merge == MergeStrategy.PreProcessSambambaMarkdup) numberSamples else 0)
111 112 113 114 115 116 117 118
      pipeline.samples.foreach {
        case (sampleName, sample) =>
          if (merge == MergeStrategy.None) sample.bamFile shouldBe None
          sample.summaryStats shouldBe Map()
          sample.libraries.foreach {
            case (libraryId, library) =>
              library.summaryStats shouldBe Map()
          }
Peter van 't Hof's avatar
Peter van 't Hof committed
119
      }
120

121
      pipesJobs.count(_.isInstanceOf[Centrifuge]) shouldBe (if (unmappedToGears) (numberFastqLibs + numberSamples) else 0)
122

Peter van 't Hof's avatar
Peter van 't Hof committed
123
      pipeline.summarySettings.get("merge_strategy") shouldBe Some(merge.toString)
Peter van 't Hof's avatar
Peter van 't Hof committed
124 125
    }
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
126 127 128 129 130

  // remove temporary run directory all tests in the class have been run
  @AfterClass def removeTempOutputDir() = {
    FileUtils.deleteDirectory(outputDir)
  }
Peter van 't Hof's avatar
Peter van 't Hof committed
131 132
}

Peter van 't Hof's avatar
Peter van 't Hof committed
133 134 135 136
class MultisampleMappingTest extends MultisampleMappingTestTrait {
  override def sample1 = Array(true)
}

137 138 139 140 141 142
class MultisampleMappingNoSamplesTest extends MultisampleMappingTestTrait {
  override def sample1 = Array(false)
  override def sample2 = Array(false)
  override def mergeStrategies = MultisampleMapping.MergeStrategy.values.filter(_ == MultisampleMapping.MergeStrategy.PreProcessMarkDuplicates)
}

143 144 145 146 147 148 149
class MultisampleMappingGearsTest extends MultisampleMappingTestTrait {
  override def sample1 = Array(true)
  override def sample2 = Array(false)
  override def unmappedToGears = true
  override def mergeStrategies = MultisampleMapping.MergeStrategy.values.filter(_ == MultisampleMapping.MergeStrategy.PreProcessMarkDuplicates)
}

150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
class MultisampleMappingBamTest extends MultisampleMappingTestTrait {
  override def sample1 = Array(false)
  override def sample2 = Array(false)
  override def sample3 = true
  override def mergeStrategies = MultisampleMapping.MergeStrategy.values.filter(_ == MultisampleMapping.MergeStrategy.PreProcessMarkDuplicates)
}

class MultisampleMappingWrongBamTest extends MultisampleMappingTestTrait {
  override def sample1 = Array(false)
  override def sample2 = Array(false)
  override def sample4 = true
  override def mergeStrategies = MultisampleMapping.MergeStrategy.values.filter(_ == MultisampleMapping.MergeStrategy.PreProcessMarkDuplicates)
}

class MultisampleMappingCorrectBamTest extends MultisampleMappingTestTrait {
  override def sample1 = Array(false)
  override def sample2 = Array(false)
  override def correctReadgroups = true
  override def sample4 = true
  override def mergeStrategies = MultisampleMapping.MergeStrategy.values.filter(_ == MultisampleMapping.MergeStrategy.PreProcessMarkDuplicates)
}

class MultisampleMappingBamToFastqTest extends MultisampleMappingTestTrait {
  override def sample1 = Array(false)
  override def sample2 = Array(false)
  override def bamToFastq = true
  override def sample3 = true
  override def sample4 = true
  override def mergeStrategies = MultisampleMapping.MergeStrategy.values.filter(_ == MultisampleMapping.MergeStrategy.PreProcessMarkDuplicates)
}

Peter van 't Hof's avatar
Peter van 't Hof committed
181
object MultisampleMappingTestTrait {
Peter van 't Hof's avatar
Peter van 't Hof committed
182
  def outputDir = Files.createTempDir()
183 184 185

  val inputDir = Files.createTempDir()

Peter van 't Hof's avatar
Peter van 't Hof committed
186
  def inputTouch(name: String): File = {
187
    val file = new File(inputDir, name).getAbsoluteFile
Peter van 't Hof's avatar
Peter van 't Hof committed
188 189 190 191 192 193
    Files.touch(file)
    file
  }

  private def copyFile(name: String): Unit = {
    val is = getClass.getResourceAsStream("/" + name)
Peter van 't Hof's avatar
Peter van 't Hof committed
194 195 196
    val out = new File(inputDir, name)
    out.deleteOnExit()
    val os = new FileOutputStream(out)
Peter van 't Hof's avatar
Peter van 't Hof committed
197 198 199 200 201 202 203
    org.apache.commons.io.IOUtils.copy(is, os)
    os.close()
  }

  copyFile("ref.fa")
  copyFile("ref.dict")
  copyFile("ref.fa.fai")
204
  copyFile("empty.sam")
Peter van 't Hof's avatar
Peter van 't Hof committed
205

Peter van 't Hof's avatar
Peter van 't Hof committed
206
  def config(outputDir: File) = Map(
207
    "skip_write_dependencies" -> true,
Peter van 't Hof's avatar
Peter van 't Hof committed
208 209 210 211 212
    "name_prefix" -> "test",
    "cache" -> true,
    "dir" -> "test",
    "vep_script" -> "test",
    "output_dir" -> outputDir,
213
    "reference_fasta" -> (inputDir + File.separator + "ref.fa"),
Peter van 't Hof's avatar
Peter van 't Hof committed
214 215 216 217 218 219 220
    "fastqc" -> Map("exe" -> "test"),
    "input_alleles" -> "test",
    "fastqc" -> Map("exe" -> "test"),
    "seqtk" -> Map("exe" -> "test"),
    "sickle" -> Map("exe" -> "test"),
    "cutadapt" -> Map("exe" -> "test"),
    "bwa" -> Map("exe" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
221
    "sambamba" -> Map("exe" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
222
    "samtools" -> Map("exe" -> "test"),
223
    "igvtools" -> Map("exe" -> "test", "igvtools_jar" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
224
    "wigtobigwig" -> Map("exe" -> "test"),
225 226
    "kraken" -> Map("exe" -> "test", "db" -> "test"),
    "krakenreport" -> Map("exe" -> "test", "db" -> "test"),
Peter van 't Hof's avatar
Peter van 't Hof committed
227 228 229
    "centrifuge" -> Map("exe" -> "test"),
    "centrifugekreport" -> Map("exe" -> "test"),
    "centrifuge_index" -> "test",
Peter van 't Hof's avatar
Peter van 't Hof committed
230 231 232 233 234 235 236 237 238 239 240 241 242
    "md5sum" -> Map("exe" -> "test")
  )

  val sample1 = Map(
    "samples" -> Map("sample1" -> Map("libraries" -> Map(
      "lib1" -> Map(
        "R1" -> inputTouch("1_1_R1.fq"),
        "R2" -> inputTouch("1_1_R2.fq")
      )
    )
    )))

  val sample2 = Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
243
    "samples" -> Map("sample2" -> Map("libraries" -> Map(
Peter van 't Hof's avatar
Peter van 't Hof committed
244 245 246 247 248 249 250 251 252 253
      "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")
      )
    )
    )))
254 255 256 257

  val sample3 = Map(
    "samples" -> Map("sample3" -> Map("libraries" -> Map(
      "lib1" -> Map(
258
        "bam" -> (inputDir + File.separator + "empty.sam")
259 260 261 262 263 264 265
      )
    )
    )))

  val sample4 = Map(
    "samples" -> Map("sample4" -> Map("libraries" -> Map(
      "lib1" -> Map(
266
        "bam" -> (inputDir + File.separator + "empty.sam")
267 268 269
      )
    )
    )))
Peter van 't Hof's avatar
Peter van 't Hof committed
270
}