Question: filterVcf: error with param argument
0
gravatar for pau.puigdevall
4.2 years ago by
Spain
pau.puigdevall0 wrote:

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:

...

 

ADD COMMENTlink modified 4.2 years ago by Martin Morgan ♦♦ 23k • written 4.2 years ago by pau.puigdevall0
Answer: filterVcf: error with param argument
0
gravatar for Martin Morgan
4.2 years ago by
Martin Morgan ♦♦ 23k
United States
Martin Morgan ♦♦ 23k wrote:

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)
ADD COMMENTlink written 4.2 years ago by Martin Morgan ♦♦ 23k

Thanks for the report. Now fixed in VariantAnnotation 1.14.10 and 1.15.25.

Valerie

ADD REPLYlink written 4.2 years ago by Valerie Obenchain6.7k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 328 users visited in the last hour