filterVcf: crash with isSNV() filter argument
2
0
Entering edit mode
@paupuigdevall-8554
Last seen 8.8 years ago
Spain


Hello,

I would like to filter one vcf by selecting only SNVs and in addition, those that only match genomic coordinates of a provided GRanges. I know that the filterVcf() function of Variant Annotation can deal with those requirements by using the isSNV() function on filters and by selecting on param argument the ScanVcfParam option. However, when just selecting the SNV filter I got some message errors and I do not know how to proceed. I've been following the Bioconductor tutorial <http://bioconductor.org/help/course-materials/2014/BioC2014/Lawrence_Tutorial.R>, but still I don't know what is the error with extractROWS.

I use the chr7 default vcf in the extdata folder of the package in order to show you my problem.

 

> library(VariantAnnotation)

> fl <- system.file("extdata", "chr7-sub.vcf.gz", package="VariantAnnotation")
> destination.file <- "~/destination.vcf"

> filterVcf(fl,
+           "hg19",
+           destination.file,
+           filters=FilterRules(list(onlySNV=isSNV)),
+           verbose=TRUE)
starting filter
filtering 3791 records
Error in extractROWS(x, eval(filter, x)) :
  error in evaluating the argument 'i' in selecting a method for function 'extractROWS': Error in rule(envir) :
  call to standardGeneric("isSNV") apparently not from the body of that generic function


May I pass the isSNV function in another way?

I would appreciate you opinion.

Thank you very much,

Pau.

 

 

 

VariantAnnotation filterVcf • 1.5k views
ADD COMMENT
2
Entering edit mode
@michael-lawrence-3846
Last seen 2.4 years ago
United States

This is fixed in release (IRanges 2.2.6) and devel (S4Vectors). Will take some time to propagate.

ADD COMMENT
0
Entering edit mode
@paupuigdevall-8554
Last seen 8.8 years ago
Spain

Thanks Michael!! With the update now seems that the isSNV argument is working perfectly!

However, I'm just stacked again with this function. My goal by using this function is also to filter both with the SNVs and also by specifying genomic coordinates with the param argument using the ScanVcfParam. However, when I pass a GRanges to this function, like is done in the tutorial mentioned in the previous message, it raises an error of the writeVcf function. However, 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", destination.file, param=ScanVcfParam(which=exons_by_gene2))

Thanks in advance.

Pau

ADD COMMENT
1
Entering edit mode
ADD REPLY
0
Entering edit mode

Thanks Martin!! And it is even more efficient that I expected.

ADD REPLY
0
Entering edit mode

It would have been better to submit this as a new question. Hopefully someone more familiar with filterVcf can help out here.

ADD REPLY

Login before adding your answer.

Traffic: 418 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6