Question: Tools to filter the heterozygous SNPs with biased allele depth (AD) values in vcf file.
gravatar for rimgubaev
8 months ago by
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 8 months ago by harold.smith.tarheel4.6k • written 8 months ago by rimgubaev180
gravatar for Pierre Lindenbaum
8 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum130k wrote:

using vcffilterjdk:

 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 8 months ago • written 8 months ago by Pierre Lindenbaum130k

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 REPLYlink written 7 weeks ago by michaeldonaldson10
gravatar for harold.smith.tarheel
8 months ago by
United States
harold.smith.tarheel4.6k wrote:

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

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


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