We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

Allele Depth (AD) / Allele Balance (AB) Filtering in GATK 4


I am trying to filter my GATK 4.0.3 - HaplotypeCaller generated multi-sample VCF for allele depth (AD) annotation at sample genotype-level (so available in "FORMAT" fields of each sample).

I think prior to GATK 4, this annotation was available as "Allele Balance" (AB) ratios (generated by AlleleBalanceBySample), but it is not available anymore in GATK 4. So I tried to filter genotypes based on AD field, that is exactly the same thing but indicated in "X,Y" format, so in an array format of integers. This array format makes it difficult to filter based on depth of alternative allele divided by depth of all alleles at a specific site.

Can you please recommend any solution to this problem? If I could turn this array into a ratio, I could easily filter genotypes using VariantFiltration or other tools such as vcflib/vcffilter. I also tried the below code (following https://gatkforums.broadinstitute.org/gatk/discussion/1255/what-are-jexl-expressions-and-how-can-i-use-them-with-the-gatk):

gatk VariantFiltration -R $ref -V $vcf -O $output --genotype-filter-expression 'vc.getGenotype("Sample1").getAD().1 / vc.getGenotype("Sample1").getAD().0 > 0.33' --set-filtered-genotype-to-no-call --genotype-filter-name 'ABfilter'

This worked, but strangely it filters the variant for all samples if only one of the sample have allele depths that are not in balance (defined by the filter). If it worked only for Sample1, I was planning to write a quick loop for all the samples for instance. I tried the same with GATK 3.8, but still it filters whole variant for all the samples if it is filtered in just one sample.


  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭


    Can you post some before and example VCF records? I would not expect this behavior.


  • psb21psb21 BejaMember
    edited October 2018

    Hello, @Sheila,

    I'm having the same issue in the last GATK release. I make a filter for a single sample (last part of the command), but it is being applied to all of them:

    gatk VariantFiltration \
    --filter-name LowQD --filter-expression "QD < 2.0" \
    --filter-name LowMQ --filter-expression "MQ < 40.0" \
    --filter-name HighFS --filter-expression "FS > 60.0" \
    --filter-name LowMQRankSum --filter-expression "MQRankSum < -12.5" \
    --filter-name HighSOR --filter-expression "SOR > 3.0" \
    --filter-name LowReadPosRankSum --filter-expression "ReadPosRankSum < -8.0" \
    --set-filtered-genotype-to-no-call true \
    --genotype-filter-expression "GQ <= 20" --genotype-filter-name LowGQ \
    --genotype-filter-expression "vc.getGenotype('H1253.1').getAD().1 / vc.getGenotype('H1253.1').getDP() < 0.30 && vc.getGenotype('H1253.1').getDP() < 30" 
    --genotype-filter-name "H1253.1_LowCov"

    Which results in somethings like this in the samples fields:

    chr1    116253409       .       C       T       80.71   PASS    ANN=T|intron_variant|MODIFIER|CASQ2|ENSG00000118729|Transcript|ENST00000261448|protein_coding||8/10|ENST00000261448.5:c.839-5496G>A|||||||rs543183086||-1||SNV|HGNC|1513|YES||1||||chr1:g.116253409C>T|0.0022|0.0076|AFR||||||||4.096|0.145370||rs543183086|0.00000e+00|-5.19000005722046|0|-1.5789999961853,T|intron_variant|MODIFIER|CASQ2|ENSG00000118729|Transcript|ENST00000456138|protein_coding||6/8|ENST00000456138.2:c.626-5496G>A|||||||rs543183086||-1||SNV|HGNC|1513|||1||||chr1:g.116253409C>T|0.0022|0.0076|AFR||||||||4.096|0.145370||rs543183086|0.00000e+00|-5.19000005722046|0|-1.5789999961853;BaseQRankSum=2.08;ClippingRankSum=0;DP=39727;ExcessHet=3.0178;FS=0;InbreedingCoeff=-0.0019;MLEAC=1;MLEAF=0.0005721;MQ=60;MQRankSum=0;QD=6.21;ReadPosRankSum=-2.004;SOR=0.283      GT:AD:DP:FT:GQ:PL       ./.:22,0:22:H1253.1_LowCov:60:0,60,677  ./.:21,0:21:H1253.1_LowCov:54:0,54,810  ./.:24,0:24:H1253.1_LowCov:63:0,63,945 ....


    Post edited by shlee on
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi Pedro (@psb21),

    I think perhaps you need a || or && between these two variant contexts:

    "vc.getGenotype('H1253.1').getAD().1 / vc.getGenotype('H1253.1').getDP()

    Here are some documents you may find helpful:

    I have a vague recollection you might want to try removing the quotes around your JEXL expressions. Let me know if these work for you.

Sign In or Register to comment.