VcfFilter.scala 3.72 KB
Newer Older
Peter van 't Hof's avatar
Peter van 't Hof committed
1
package nl.lumc.sasc.biopet.tools
Peter van 't Hof's avatar
Peter van 't Hof committed
2
3
4
5
6
7

import htsjdk.variant.variantcontext.writer.AsyncVariantContextWriter
import htsjdk.variant.variantcontext.writer.VariantContextWriterBuilder
import htsjdk.variant.vcf.VCFFileReader
import java.io.File
import nl.lumc.sasc.biopet.core.BiopetJavaCommandLineFunction
Peter van 't Hof's avatar
Peter van 't Hof committed
8
import nl.lumc.sasc.biopet.core.ToolCommand
Peter van 't Hof's avatar
Peter van 't Hof committed
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import nl.lumc.sasc.biopet.core.config.Configurable
import org.broadinstitute.gatk.utils.commandline.{ Output, Input }
import scala.collection.JavaConversions._

class VcfFilter(val root: Configurable) extends BiopetJavaCommandLineFunction {
  javaMainClass = getClass.getName

  @Input(doc = "Input vcf", shortName = "I", required = true)
  var inputVcf: File = _
    
  @Output(doc = "Output vcf", shortName = "o", required = false)
  var outputVcf: File = _
  
  var minSampleDepth: Option[Int] = _
  var minTotalDepth: Option[Int] = _
  var minAlternateDepth: Option[Int] = _
  var minSamplesPass: Option[Int] = _
Peter van 't Hof's avatar
Peter van 't Hof committed
26
  var filterRefCalls: Boolean = _
Peter van 't Hof's avatar
Peter van 't Hof committed
27
28
29
30
31
32
33
34
35
  
  override val defaultVmem = "8G"
  memoryLimit = Option(4.0)
  
  override def afterGraph {
    minSampleDepth = config("min_sample_depth")
    minTotalDepth = config("min_total_depth")
    minAlternateDepth = config("min_alternate_depth")
    minSamplesPass = config("min_samples_pass")
Peter van 't Hof's avatar
Peter van 't Hof committed
36
    filterRefCalls = config("filter_ref_calls")
Peter van 't Hof's avatar
Peter van 't Hof committed
37
38
39
40
41
42
43
44
  }
  
  override def commandLine = super.commandLine + 
    required("-I", inputVcf) + 
    required("-o", outputVcf) +
    optional("-minSampleDepth", minSampleDepth) +
    optional("-minTotalDepth", minTotalDepth) +
    optional("-minAlternateDepth", minAlternateDepth) + 
Peter van 't Hof's avatar
Peter van 't Hof committed
45
46
    optional("-minSamplesPass", minSamplesPass) +
    conditional(filterRefCalls, "-filterRefCalls")
Peter van 't Hof's avatar
Peter van 't Hof committed
47
48
}

Peter van 't Hof's avatar
Peter van 't Hof committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
object VcfFilter extends ToolCommand {
  case class Args (inputVcf:File = null, outputVcf:File = null, minSampleDepth: Int = -1, minTotalDepth: Int = -1,
                   minAlternateDepth: Int = -1, minSamplesPass: Int = 0, filterRefCalls: Boolean = false) extends AbstractArgs

  class OptParser extends AbstractOptParser {
    opt[File]('I', "inputVcf") required() maxOccurs(1) valueName("<file>") action { (x, c) =>
      c.copy(inputVcf = x) }
    opt[File]('o', "outputVcf") required() maxOccurs(1) valueName("<file>") action { (x, c) =>
      c.copy(outputVcf = x) } text("output file, default to stdout")
    opt[Int]("minSampleDepth") unbounded() action { (x, c) =>
      c.copy(minSampleDepth = x ) }
    opt[Int]("minTotalDepth") unbounded() action { (x, c) =>
      c.copy(minTotalDepth = x ) }
    opt[Int]("minAlternateDepth") unbounded() action { (x, c) =>
      c.copy(minAlternateDepth = x) }
    opt[Int]("minSamplesPass") unbounded() action { (x, c) =>
      c.copy(minSamplesPass = x) }
    opt[Boolean]("filterRefCalls") unbounded() action { (x, c) => 
      c.copy(filterRefCalls = x) }
  }
  
Peter van 't Hof's avatar
Peter van 't Hof committed
70
71
72
73
  /**
   * @param args the command line arguments
   */
  def main(args: Array[String]): Unit = {
Peter van 't Hof's avatar
Peter van 't Hof committed
74
75
    val argsParser = new OptParser
    val commandArgs: Args = argsParser.parse(args, Args()) getOrElse sys.exit(1)
Peter van 't Hof's avatar
Peter van 't Hof committed
76
    
Peter van 't Hof's avatar
Peter van 't Hof committed
77
78
    val reader = new VCFFileReader(commandArgs.inputVcf, false)
    val writer = new AsyncVariantContextWriter(new VariantContextWriterBuilder().setOutputFile(commandArgs.outputVcf).build)
Peter van 't Hof's avatar
Peter van 't Hof committed
79
    writer.writeHeader(reader.getFileHeader)
Peter van 't Hof's avatar
Peter van 't Hof committed
80
81
    for (record <- reader) {
      val genotypes = for (genotype <- record.getGenotypes) yield {
Peter van 't Hof's avatar
Peter van 't Hof committed
82
83
84
        genotype.getDP >= commandArgs.minSampleDepth && 
        List(genotype.getAD:_*).tail.count(_ >= commandArgs.minAlternateDepth) > 0 &&
        !(commandArgs.filterRefCalls && genotype.isHomRef)
Peter van 't Hof's avatar
Peter van 't Hof committed
85
86
      }
      
Peter van 't Hof's avatar
Peter van 't Hof committed
87
      if (record.getAttributeAsInt("DP", -1) >= commandArgs.minTotalDepth && genotypes.count(_ == true) >= commandArgs.minSamplesPass)
Peter van 't Hof's avatar
Peter van 't Hof committed
88
89
90
91
92
93
        writer.add(record)
    }
    reader.close
    writer.close
  }
}