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

Hi,

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.

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @maegsul
    Hi,

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

    Thanks,
    Sheila

  • psb21psb21 BejaMember

    Hello, @Sheila,

    I'm having the same issue in the last GATK 4.0.9.0 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 ....

    Best,
    Pedro

Sign In or Register to comment.