Hello,

I am new to the VariantAnnotation package. I am trying to filter a 4-species vcf with the below code. I am having 2 problems: 1) I am not recovering the one site in my test file where all 4 species have RGQ>=30 (along with with variant sites that are not filtered), and 2) the header is of the vcf file is being printed multiple times (~55 times). I will provide any further information that would be helpful.

many thanks,

Stephen

code:

library("VariantAnnotation")

ref <- "ref"
x <- readVcf("test.vcf.gz", genome=ref)

##prefilter
##only those that have RGQ
isRGQ <- function(x) {
    grepl("RGQ", x, fixed=TRUE)
}

##filters
RGQ_filter <- function(x, value=30) {
    ## Filter on RGQ > value; default is 30

    test <- apply(geno(x)$RGQ>=value, 1, sum)==4
    out <- is.na(test) & test  
    out[is.na(out)] <- TRUE
    return(out)
}

prefilters <- FilterRules(list(RGQ_prefilter=isRGQ))
filters <- FilterRules(list(RGQ_filter=RGQ_filter))

file.gz <- "test.vcf.gz"
file.gz.tbi <- "test.vcf.gz.tbi"
destination.file <- "test2.vcf.gz"
size=1
tabix.file <- TabixFile(file.gz, yieldSize=size)
destination.file <- tempfile()

filterVcf(file=tabix.file, genome=ref, destination=destination.file,
          filters=filters, verbose=TRUE)

This is now acting as expected. I stopped using the tabix.file <- TabixFile(file.gz, yieldSize=size), and use the actual file name.

filter has been modified to work as expected:

RGQ_filter <- function(x, value=30) {
    ## Filter on RGQ > value; default is 30
    ##test <- apply(geno(x)$RGQ>=value, 1, sum)==4
    ##out <- is.na(test) & test  
    out <- apply(geno(x)$RGQ>=value, 1, sum)==4
    out[is.na(out)] <- TRUE
    return(out)
}

HTH,

Stephen