text: user1769: 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?
  
  user492: 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.
  
  user1024: 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!
  
  user71: 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 
  
  
original_text: 1769: <p>I have high-depth variant calling created using the HaplotypeCaller with <code>--output_mode EMIT_ALL_SITES</code> I'm interested in finding all sites (regardless of genotype call heterozygous or homozygous) where at least one of the <strong>alternative</strong> alleles have an <code>AD</code> value (Allelic Depth) greater than 10, <em>I</em>.<em>e</em>. 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.</p>
  <p>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 <code>1/1:1,988:989</code> <code>0/1:116,92</code> <code>0/1:220,234</code> <code>0/1:62,611</code> <code>1/1:0,109</code> respectively.</p>
  <pre><code>#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
  </code></pre>
  <p><a href="https://www.dropbox.com/s/gm4pftuk9i1gx6p/test.vcf?dl=1" rel="nofollow noreferrer">A dropbox link for file</a></p>
  <p>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 &gt; 10 which will give me both ref and alt alleles with AD &gt; 10.  Perhaps there is a bioawk solution or something more elaborate with coreutils which could successfully return sites with an alt AD count &gt; 10?</p>
   492: <p>This now works with the development version of Bcftools v1.5 (commit 4f134df). Thanks to Petr Danecek for <a href="https://github.com/samtools/bcftools/issues/639" rel="nofollow noreferrer">adding the feature</a>. I expect this feature to make its way into the next release of Bcftools:</p>
  <pre><code>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-] &gt; 10' Test.vcf.gz 
  </code></pre>
  <p>Output without header (I have modified the second line to be tri-allelic to demonstrate the filtering works):</p>
  <pre><code>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
  </code></pre>
  <p><strong>Edited</strong> (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 <code>'AD[1-] &gt; 10'</code> would produce the error <code>The FORMAT tag AD can have multiple subfields, run as AD[sample:subfield]</code>. A discussion of these changes can be found <a href="https://github.com/samtools/bcftools/issues/757" rel="nofollow noreferrer">here</a>.</p>
   1024: <p>You can do this in <a href="https://hail.is" rel="nofollow noreferrer">Hail</a>:</p>
  
  <pre><code>from hail import *
  hc = HailContext()
  (hc.import_vcf('test.vcf')
   .filter_variants_expr('gs.exists(g =&gt; g.ad[1:].exists(d =&gt; d &gt; 10))')
   .export_vcf('filtered.vcf'))
  </code></pre>
  
  <p>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.</p>
  
  <p>To verify we got the expected 6 variants:</p>
  
  <pre><code>&gt;&gt;&gt; hc.import_vcf('filtered.vcf').count()
  (1L, 6L)
  </code></pre>
  
  <p><a href="https://hail.is/hail/hail.VariantDataset.html#hail.VariantDataset.count" rel="nofollow noreferrer"><code>count</code></a> returns the number of samples (1) and number of variants (6).</p>
  
  <p>Take a look at the <a href="https://hail.is/hail/getting_started.html" rel="nofollow noreferrer">getting started page</a> or <a href="https://hail.is/hail/tutorials-landing.html" rel="nofollow noreferrer">tutorials</a> if you want to try it out!</p>
   71: <p>using <a href="http://lindenb.github.io/jvarkit/VCFFilterJS.html" rel="nofollow noreferrer">vcfilterjs</a></p>
  
  <p>and the following script:</p>
  
  <pre><code>function accept(vc)
      {
      var i,j;
      for(i=0;i&lt; 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&lt; ad.length ;++j)
              {
              if(ad[j]&gt;10) return true;
              }
          }
      return false;
      }
  accept(variant);
  </code></pre>
  
  <p>usage:</p>
  
  <pre><code>java -jar dist/vcffilterjs.jar -f script.js Test.vcf 
  </code></pre>
  