# Selecting sites from VCF which have an alt AD > 10

I have high-depth variant calling created using the HaplotypeCaller with --output_mode EMIT_ALL_SITES I'm interested in finding all sites (regardless of genotype call heterozygous or homozygous) where at least one of the alternative alleles have an AD value (Allelic Depth) greater than 10, I.e. are supported by more than 10 reads. Also ideally I want back more than just the first alternative allele. Note that I don't want back lines of VCF were we only see an AD count for the ref allele only.

So in the example VCF snippet below I'm wanting to select lines: 6,7,8,12,13 and 14, which have GT:AD values 1/1:1,988:989 0/1:116,92 0/1:220,234 0/1:62,611 1/1:0,109 respectively.

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  12908_DIAG
3   187446740   .   T   .   Infinity    .   AN=2;DP=1095;MQ=60.00   GT:AD:DP    0/0:1095:1095
3   187446741   .   C   .   Infinity    .   AN=2;DP=1117;MQ=60.00   GT:AD:DP    0/0:1117:1117
3   187446752   .   A   .   Infinity    .   AN=2;DP=1297;MQ=60.00   GT:AD:DP    0/0:1297:1297
3   187446763   .   C   .   Infinity    .   AN=2;DP=1494;MQ=60.00   GT:AD:DP    0/0:1494:1494
3   187451574   .   C   .   Infinity    .   AN=2;DP=1493;MQ=60.00   GT:AD:DP    0/0:1493:1493
4   1805296 rs3135883   G   A   3876.03 .   AC=2;AF=1.00;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.22;SOR=9.401    GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0


I'd initially considered using GATK's SelectVariants but I'm not sure JEXL has the ability to select out what I want specifically other than a blanket AD > 10 which will give me both ref and alt alleles with AD > 10. Perhaps there is a bioawk solution or something more elaborate with coreutils which could successfully return sites with an alt AD count > 10?

• – Pierre Jul 3 '17 at 16:23

using vcfilterjs

and the following script:

function accept(vc)
{
var i,j;
for(i=0;i< vc.getNSamples();++i)
{
var genotype = vc.getGenotype(i);
/* loop over AD starting from '1' =  first ALT */
{
}
}
return false;
}
accept(variant);


usage:

java -jar dist/vcffilterjs.jar -f script.js Test.vcf

• I think this system you've created looks to be very flexible I'd be very interested in seeing some more walkthroughs and documented examples, on Github because I'm sure it can solve a lot of complex variant filtration issues. – Matt Bashton Jul 5 '17 at 16:36
• @MatthewBashton thanks ! and there is a much more faster alternative version since today (not javascript , but 100% java based ) lindenb.github.io/jvarkit/VcfFilterJdk.html – Pierre Jul 5 '17 at 17:06
• @DanielStandage not braces, just lambdas using vcffilterjdk :  java -jar dist/vcffilterjdk.jar -e 'return variant.getGenotypes().stream().filter(G->G.hasAD() && java.util.Arrays.stream(G.getAD()).skip(1).filter(AD->AD>10).findAny().isPresent()).findAny().isPresent();' input.vcf  – Pierre Jul 5 '17 at 17:32
• Nice I'll have to give the jdk version a go too. I have no issues your braces or indentation style ¯\_(ツ)_/¯ – Matt Bashton Jul 5 '17 at 21:45

You can do this in Hail:

from hail import *
hc = HailContext()
(hc.import_vcf('test.vcf')
.filter_variants_expr('gs.exists(g => g.ad[1:].exists(d => d > 10))')
.export_vcf('filtered.vcf'))


This works with any number of samples and will keep the variants where at least one sample has a genotype with an alternate allele support by more than 10 reads.

To verify we got the expected 6 variants:

>>> hc.import_vcf('filtered.vcf').count()
(1L, 6L)


count returns the number of samples (1) and number of variants (6).

Take a look at the getting started page or tutorials if you want to try it out!

This now works with the development version of Bcftools v1.5 (commit 4f134df). Thanks to Petr Danecek for adding the feature. I expect this feature to make its way into the next release of Bcftools:

git clone git://github.com/samtools/htslib.git
git clone git://github.com/samtools/bcftools.git
(cd bcftools; make)

bgzip Test.vcf
./bcftools/bcftools index Test.vcf.gz
./bcftools/bcftools filter -i 'AD[1-] > 10' Test.vcf.gz


Output without header (I have modified the second line to be tri-allelic to demonstrate the filtering works):

3   187451609   rs1880101   A   G   39794   PASS    AC=2;AF=1;AN=2;BaseQRankSum=1.859;ClippingRankSum=0;DB;DP=995;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=60;MQRankSum=0;QD=24.56;ReadPosRankSum=0.406;SOR=8.234   GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0