How can I filter a vcf file by DP and GQ using the R package VariantAnnotation?

Question

[r.jaffe]

Hi guys, need some help filtering my vcf files with the VariantAnnotation package (v1.12.9) from Bioconductor. I have a vcf file (v4.2) containing 28238 SNPs for 160 individuals, and I'm running R version 3.1.3

I opened my file and ran:

hist(geno(vcf)$DP) ##check the distribution of Read Depth
hist(geno(vcf)$GQ) ##check the distribution of Phred-Scaled Genotype Quality

I then tried to filter out all samples with a DP<10 and a GQ<20 and create a new vcf file:

vcf2 <- vcf[geno(vcf)$DP>10 && geno(vcf)$GQ>20]

The problem is that this code does not seem to be working, because both files have exactly the same dimensions:

dim(vcf) ##28238   160
dim(vcf2) ##28238   160

So my questions are:

1. How can I filter vcf files by DP and GQ?
2. Is it better to eliminate SNPs that have low DP & GQ, or to re-code genotypes with low DP and GQ as missing genotypes (./.)?

Any help would be much appreciated, thanks!

Rodolfo

[Martin Morgan]

To fix your code, use & (vector and) rather than && (scalar and) as in response to your other version of this question.

Answer 1

[NB]

Hi,

Not sure how to do it in R but I usually use vcftools to filter my variants based on DP and GQ

vcftools --vcf myfile.vcf --minGQ 20 --minDP 10 --recode --out new_myfile

You can also use GATK's --filterExpression to do the same.

[r.jaffe]

Thank you sukmb!

Yes I used vcftools and bcftools already, but was wandering if I could do it in R.

Answer 2

[Jeremy Leipzig]

when you look at this:

str(geno(vcf)$DP>10)

it probably returns a logical matrix

logi [1:28238, 1:160]

if you use that to subset a data.frame how is R supposed to behave? subset rows, columns, or cells? How do you subset cells?

Maybe you meant to filter by combined depth across samples, which is in info

[r.jaffe]

Thank you Jeremy. Yes I get a logical matrix, but I still don`t know how to subset the rows of my vcf...

[Jeremy Leipzig]

I think you need to be clear about what you are asking for - samples, snps, or sample/snp incidents