Selecting sites from VCF which have an alt AD > 10

Question

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
3   187451609   rs1880101   A   G   39794.03    .   AC=2;AF=1.00;AN=2;BaseQRankSum=1.859;ClippingRankSum=0.000;DB;DP=995;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=0.000;QD=24.56;ReadPosRankSum=0.406;SOR=8.234  GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0
4   1803279 .   T   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0.000;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0.00;MQ=59.97;MQRankSum=0.000;ReadPosRankSum=-2.523;SOR=6.357   GT:AD:DP:GQ:PL  0/0:211,23:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0.000;DB;DP=215;ExcessHet=3.0103;FS=1.110;MLEAC=1;MLEAF=0.500;MQ=59.97;MQRankSum=0.000;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809  GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803671 .   C   A   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-0.880;ClippingRankSum=0.000;DP=450;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.953;SOR=0.572    GT:AD:DP:GQ:PL  0/0:445,2:447:99:0,1272,15958
4   1803681 .   T   C   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.654;ClippingRankSum=0.000;DP=483;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.422;SOR=0.664    GT:AD:DP:GQ:PL  0/0:479,2:481:99:0,1408,18538
4   1803703 .   A   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.704;ClippingRankSum=0.000;DP=458;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=0.299;SOR=0.497 GT:AD:DP:GQ:PL  0/0:454,2:456:99:0,1325,18095
4   1803704 rs2234909   T   C   6676.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0.000;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=14.71;ReadPosRankSum=0.324;SOR=0.849   GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=8.083;ClippingRankSum=0.000;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=16.51;ReadPosRankSum=0.180;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
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

A dropbox link for file

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?

cross-posted: gatkforums.broadinstitute.org/gatk/discussion/comment/… — Pierre, Commented Jul 3, 2017 at 16:23

Matt Bashton · Accepted Answer · 2020-12-03 20:45:58Z

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[0: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
4   1803279 .   T   G   0   PASS    AC=0;AF=0;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0;MQ=59.97;MQRankSum=0;ReadPosRankSum=-2.523;SOR=6.357 GT:AD:DP:GQ:PL  0/0/0:211,3,34:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0;DB;DP=215;ExcessHet=3.0103;FS=1.11;MLEAC=1;MLEAF=0.5;MQ=59.97;MQRankSum=0;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809   GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803704 rs2234909   T   C   6676.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=14.71;ReadPosRankSum=0.324;SOR=0.849  GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=8.083;ClippingRankSum=0;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=16.51;ReadPosRankSum=0.18;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
4   1805296 rs3135883   G   A   3876.03 PASS    AC=2;AF=1;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=60;QD=29.22;SOR=9.401 GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0

Edited (Dec 2020): bcftools filter syntax has been updated, answer now reflects changes found in version v1.11 (not sure at what point these were introduced). The old query 'AD[1-] > 10' would produce the error The FORMAT tag AD can have multiple subfields, run as AD[sample:subfield]. A discussion of these changes can be found here.

So just revisiting this and the query now errors in current versions of bcftools: The FORMAT tag AD can have multiple subfields, run as AD[sample:subfield] — Matt Bashton, Commented Dec 3, 2020 at 18:43
Okay I think the answer is now: 'AD[0:1-] > 10' just testing... — Matt Bashton, Commented Dec 3, 2020 at 19:04

Pierre · Accepted Answer · 2017-07-05 16:24:32Z

5

using vcfilterjs

and the following script:

function accept(vc)
    {
    var i,j;
    for(i=0;i< vc.getNSamples();++i)
        {
        var genotype = vc.getGenotype(i);
        if(!genotype.hasAD()) continue;
        var ad = genotype.getAD();
        /* loop over AD starting from '1' =  first ALT */
        for(j=1  ;j< ad.length ;++j)
            {
            if(ad[j]>10) return true;
            }
        }
    return false;
    }
accept(variant);

usage:

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

edited Jul 5, 2017 at 16:24

answered Jul 3, 2017 at 15:41

Pierre

1,5768 silver badges11 bronze badges

1

$\begingroup$ 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. $\endgroup$
– Matt Bashton
Commented Jul 5, 2017 at 16:36
$\begingroup$ @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 $\endgroup$
– Pierre
Commented Jul 5, 2017 at 17:06
1

$\begingroup$ @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 ` $\endgroup$
– Pierre
Commented Jul 5, 2017 at 17:32
1

$\begingroup$ Nice I'll have to give the jdk version a go too. I have no issues your braces or indentation style ¯\_(ツ)_/¯ $\endgroup$
– Matt Bashton
Commented Jul 5, 2017 at 21:45

Add a comment |

Cotton Seed · Accepted Answer · 2017-07-09 19:15:01Z

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!

Stack Exchange Network

Selecting sites from VCF which have an alt AD > 10

3 Answers 3

Your Answer

Not the answer you're looking for? Browse other questions tagged
variant-calling
vcf
bioawk
or ask your own question.

Hot Network Questions

Selecting sites from VCF which have an alt AD > 10

3 Answers 3

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged variant-callingvcfbioawk or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
variant-calling
vcf
bioawk
or ask your own question.