rbioinformaticsbioconductorgenomevcf-variant-call-format

How to save Variant Call Format (VCF) file to disk in R using VariantAnnotation Package

I've searched the web for this without much luck. More or less you always get to the example from the VariantAnnotation Package. And since this example works fine on my computer I have no idea why the VCF I created does not.

The problem: I want to determine the number and location of SNPs in selected genes. I have a large VCF file (over 5GB) that has info on all SNPs on all chromosomes for several mice strains. Obviously my computer freezes if I try to do anything on the whole genome scale, so I first determined genomic locations of genes of interest on chromosome 1. I then used the VariantAnnotation Package to get only the data relating to my genes of interest out of the VCF file:

library(VariantAnnotation)
param<-ScanVcfParam(
  info=c("AC1","AF1","DP","DP4","INDEL","MDV","MQ","MSD","PV0","PV1","PV2","PV3","PV4","QD"), 
  geno=c("DP","GL","GQ","GT","PL","SP","FI"),
  samples=strain, 
  fixed="FILTER",
  which=gnrng
  )

The code above is taken out of a function I wrote which takes strain as an argument. gnrng refers to a GRanges object containing genomic locations of my genes of interest. 

vcf<-readVcf(file, "mm10",param)

This works fine and I get my vcf (dim: 21783 1) but when I try to save it won't work

file.vcf<-tempfile()
writeVcf(vcf, file.vcf)
Error in .pasteCollapse(ALT, ",") : 'x' must be a CharacterList

I even tried in parallel, doing the example from the package first and then substituting for my VCF file:

#This is the example:
out1.vcf<-tempfile()
in1<-readVcf(fl,"hg19")
writeVcf(in1,out1.vcf)

This works just fine, but if I only substitute in1 for my vcf I get the same error.

I hope I made myself clear... And any help will be greatly appreciated!! Thanks in advance!

Solution

  • Thanks for reporting this bug. The problem is fixed in version 1.9.47 (devel branch). The fix will be available in the release branch after April 14.

    The problem was that you selectively imported 'FILTER' from the 'fixed' field but not 'ALT'. writeVcf() was throwing an error because there was no ALT value to write out. If you don't have access to the version with the fix, a work around would be to import the ALT field.

    ScanVcfParam(fixed = c("ALT", "FILTER"))

    You can see what values were imorted with the fixed() accessor:

    fixed(vcf)

    Please report and bugs or problems on the Bioconductor mailing list Martin referenced. More Bioc users will see the question and you'll get help more quickly.

    Valerie