This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Filtering VCF files based on VAF giving incorrect results

I'm trying to filter variants that have VAF<0.2. I used picard FilterVcf, bcftools filter, and tried to use GATK VariantFiltration, and neither of these worked. picard and bcftools only filter a portion of variants but I still get a long list of variants with VAF<0.2 labeled as "PASS". I haven't figured out how to use GATK VariantFiltration for VAF. Any help to figure out the problem would be greatly appreciated. Here are my codes:

java -jar picard.jar FilterVcf  \
    I=input.vcf
    MIN_AB=0.2 \
    O=output.vcf   

bcftools filter --include '(FMT/AD[0:1])/(FMT/AD[0:0]+FMT/AD[0:1]) >=0.25 && (FMT/DP[0:0] >=10)' \
        input.vcf \
        --output output.vcf \
        --soft-filter Filter1 \
        --mode +
gatk picard vcf bcftools
• 1,817 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Arton • Arton

1 answer

using vcfilterjdk: https://jvarkit.readthedocs.io/en/latest/VcfFilterJdk/

 java -jar ${JVARKIT_DIST}/jvarkit.jar vcffilterjdk -e 'return !variant.getGenotypes().stream().filter(G->G.isHet() && G.hasAD() && G.getAD().length==2).anyMatch(G->{double n1=G.getAD()[0];double n2=G.getAD()[1];if(n1+n2<=0) return false; double f=(n2/(n1+n2));return (f<0.2 || f>0.8);});'  --filter Filter1 input.vcf
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thank you. I can't try a completely new tool but I will keep your command for future use. Do you think my commands were correct? For example what does mean that MIN_AB is double?

• 0 views • link
written by Arton • Arton

Do you think my commands were correct?

I'm not familiar enough with this bcftools syntax.

but what if there is no AD, what if both are '0', what if there is a NO-CALL at this position ? what if the variant is multiallelic ?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this