Tools to filter the heterozygous SNPs with biased allele depth (AD) values in vcf file.

Question: Tools to filter the heterozygous SNPs with biased allele depth (AD) values in vcf file.

8 months ago by

rimgubaev • 180

Russia/Moscow/Skoltech

rimgubaev • 180 wrote:

Hello everyone! There is the following problem in my vcf file produced by GATK: a large proportion of heterozygous sites has a bias in allele depth values (AD). For example, AD is 43,7 which means that the 43 reads support one allele while 7 support another. In order to avoid possible false calls in GWAS, I want to switch such SNPs to no-calls (./.).

I know that it is possible to write a heavy-duty R or Python script substituting the heterozygous calls (0/1) for no-calls (./.) in case of significant AD bias, let's say when minor allele to major allele ratio is greater than 0.3. However, I wonder if there are ready-to-go solutions. I tried to find such a filter in GATK but failed.

In case if you faced a similar problem or know about the ready solution please share!

heterozygous snps ad allele depth gatk vcf • 397 views

ADD COMMENT • link •

modified 8 months ago by harold.smith.tarheel • 4.6k • written 8 months ago by rimgubaev • 180

8 months ago by

Pierre Lindenbaum ♦ 130k

France/Nantes/Institut du Thorax - INSERM UMR1087

Pierre Lindenbaum ♦ 130k wrote:

using vcffilterjdk: http://lindenb.github.io/jvarkit/VcfFilterJdk.html

 java -jar dist/vcffilterjdk.jar -rc -e 'return new VariantContextBuilder(variant).genotypes( variant.getGenotypes().stream().map(G->{ if(!G.hasAD() || !G.isHet()) return G; final int ad[]=G.getAD(); if(ad.length!=2) return G; final double treshold = 0.3; final double R= ad[0]; final double A= ad[1]; final double sum = R+A; if(sum == 0 || A/sum < treshold || (1.0 - A/sum) < treshold) { return GenotypeBuilder.createMissing(G.getSampleName(),G.getPloidy()); } return G; }).collect(Collectors.toList()) ).make();'  input.vcf

ADD COMMENT • link modified 8 months ago • written 8 months ago by Pierre Lindenbaum ♦ 130k

Hi Pierre,

Can you clarify how the above script relates to allele balance (i.e., the proportion of reads corresponding to the reference/alternate allele), across heterozygotes? I would like to change any heterozygous call with AB<0.25 & AB>0.75 to a no-call (./.). The vcf file I am working with does not have the AB field but it does have AD.

Thank you!

ADD REPLY • link written 7 weeks ago by michaeldonaldson • 10

8 months ago by

harold.smith.tarheel • 4.6k

United States

harold.smith.tarheel • 4.6k wrote:

VCFtools '--freq' flag will output allele frequencies; you can then parse those values using 'awk' to replace the genotype call.

ADD COMMENT • link written 8 months ago by harold.smith.tarheel • 4.6k

Please log in to add an answer.

Similar posts • Search »