8
$\begingroup$

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?

$\endgroup$
5
$\begingroup$

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 
$\endgroup$
  • 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 Jul 5 '17 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 Jul 5 '17 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 Jul 5 '17 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 Jul 5 '17 at 21:45
4
$\begingroup$

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!

$\endgroup$
2
$\begingroup$

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
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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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