Question: Tools to filter the heterozygous SNPs with biased allele depth (AD) values in vcf file.
0
gravatar for rimgubaev
4 months ago by
rimgubaev180
Russia/Moscow/Skoltech
rimgubaev180 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!

ADD COMMENTlink modified 4 months ago by harold.smith.tarheel4.5k • written 4 months ago by rimgubaev180
2
gravatar for Pierre Lindenbaum
4 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum128k 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 COMMENTlink modified 4 months ago • written 4 months ago by Pierre Lindenbaum128k
0
gravatar for harold.smith.tarheel
4 months ago by
United States
harold.smith.tarheel4.5k wrote:

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

ADD COMMENTlink written 4 months ago by harold.smith.tarheel4.5k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1267 users visited in the last hour