multisample vcf filter with bcftools, condition true for ALL samples

Question

[Axzd]

Hello, I have a vcf file with 11 samples.

I would like to keep sites where, for all samples, this condition is true:

MAX(AD)/SUM(AD) <= 0.6 & MAX(AD)/SUM(AD) >=0.4

In other words, I want to keep sites where, for all samples, the allele frequencies to support the SNP are between 0.4 et 0.6.

I am using this command

bcftools view -i 'MAX(AD)/SUM(AD) <= 0.6 & MAX(AD)/SUM(AD) >=0.4' biallelicSNP.vcf.gz > biallelicSNP.freq04_06.vcf

However, this results in a very low number of sites (compared to the initial number). Therefore, I am wondering if bcftools is not doing here something else than I think it does. What's your opinion? Thank you

EDIT: I think I was wrong, and the proper way of doing it is through VAF tag, so I did

bcftools +fill-tags biallelicSNP.vcf.gz -- -t FORMAT/VAF > biallelicSNP.VAF.vcf
bcftools view -i 'VAF>= 0.4 & VAF <= 0.6' biallelicSNP.VAF.vcf > biallelicSNP.freq04_06.vcf 

I would still apreciate some opinion, especially on why the first method didn't work.

Thanks!

EDIT2: actually it doesn't work as it doesn't check the condition is satisfied for ALL samples.

Answer 1

[Pierre Lindenbaum]

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

 java -jar dist/jvarkit.jar vcffilterjdk -e 'return variant.getGenotypes().stream().filter(G->G.hasAD() && G.getAD().length==2).allMatch(G->{final int[] ad=G.getAD();double sum=ad[0]+ad[1];int max=Math.max(ad[0],ad[1]);if(sum==0) return false; double f = max/sum;return f<=0.6 && f>=0.4;});' in.vcf

view raw README.md hosted with ❤ by GitHub

biostars admin: I cannot post code anymore (lang blablabla is not supported)

[Axzd]

Thanks Pierre, I am going to use it for comparison

[Pierre Lindenbaum]

let's try...

java -jar dist/jvarkit.jar vcffilterjdk -e 'return variant.getGenotypes().stream().filter(G->G.hasAD() && G.getAD().length==2).allMatch(G->{final int[] ad=G.getAD();double sum=ad[0]+ad[1];int max=Math.max(ad[0],ad[1]);if(sum==0) return false; double f = max/sum;return f<=0.6 && f>=0.4;});' in.vcf

[Pierre Lindenbaum]

ok. bablablablablabalabala weird bug

[Axzd]

Both commands work equally, at least for me

[GenoMax]

It is the same command. Just posted two different ways.