I am working with the filterVcf() function of the VariantAnnotation package.

My goal is filtering variants by selecting SNVs that matches some genomic coordinates by passing the functions isSNV() and ScanVcfParam() as respective filters and param arguments. I have followed the tutorial of Lawrence (http://bioconductor.org/help/course-materials/2014/BioC2014/Lawrence_Tutorial.R).

However, when I pass a GRanges to this function, like is done in the tutorial, it raises an error of the writeVcf function which I don't know how to proceed. On the other hand, when I do a test using the function with readVcf() and specifying the same param, it seems that everything goes smoothly. May I be doing something wrong?

>library(VariantAnnotation)
>library(TxDb.Hsapiens.UCSC.hg19.knownGene)

> txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
> exons_by_gene <- exonsBy(txdb,"gene")
> exons_by_gene <- resize(exons_by_gene, width(exons_by_gene)+101, fix="start", use.names=TRUE, ignore.strand=FALSE)
> exons_by_gene <- resize(exons_by_gene, width(exons_by_gene)+101, fix="end", use.names=TRUE, ignore.strand=FALSE)
> exons_by_gene2 <- unlist(exons_by_gene)
> exons_by_gene2 <- exons_by_gene2[seqnames(exons_by_gene2)=="chr22"]
> seqlevelsStyle(exons_by_gene2) <- "NCBI"
>
> fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
> destination.file <-"~/destination.vcf"
>
> filterVcf(fl,
+           "hg19",
+           destination.file,
+           filters=FilterRules(list(onlySNV=isSNV)),
+           verbose=TRUE,
+           param=ScanVcfParam(which=exons_by_gene2))
starting filter
found header lines for 3 ‘fixed’ fields: ALT, QUAL, FILTER
found header lines for 22 ‘info’ fields: LDAF, AVGPOST, ..., VT, SNPSOURCE
found header lines for 3 ‘geno’ fields: GT, DS, GL
filtering 3594 records
Error in writeVcf(vcfChunk, filtered) :
  error in evaluating the argument 'filename' in selecting a method for function 'writeVcf': Error: object 'filtered' not found

> test <- readVcf(fl, "hg19", param=ScanVcfParam(which=exons_by_gene2))
found header lines for 3 ‘fixed’ fields: ALT, QUAL, FILTER
found header lines for 22 ‘info’ fields: LDAF, AVGPOST, ..., VT, SNPSOURCE
found header lines for 3 ‘geno’ fields: GT, DS, GL

> rowRanges(test)
GRanges object with 3594 ranges and 5 metadata columns:

> test2 <- readVcf(fl, "hg19")

found header lines for 3 ‘fixed’ fields: ALT, QUAL, FILTER
found header lines for 22 ‘info’ fields: LDAF, AVGPOST, ..., VT, SNPSOURCE
found header lines for 3 ‘geno’ fields: GT, DS, GL

> rowRanges(test2)
GRanges object with 10376 ranges and 5 metadata columns:

...

This seems to be a bug in filterVcf. A work-around (and maybe more efficient for a large file and many ranges...) is to iterate through the entire file and use the range criterion as a second filter. Here's a factory for making a range-based filter

withinRange <- function(rng)
    function(x) x %within% rng

and then it's use

filters <- FilterRules(list(
    onlySNV=isSNV,
    withinRange=withinRange(exons_by_gene2)))
dest <- filterVcf(fl, "hg19", tempfile(), filters=filters, verbose=TRUE)