This is a test version of Biostars. For the public version, visit https://www.biostars.org.
how to switch REF and ALT alleles in VCF files if the REF is incorrect according to RefSeq?

I have a VCF file, and I believe it is converted from PED file by PLINK, as illustrated in this blog:

There is one comment saying ##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome"> in the VCF file.

For a some variant loci, the REF and ALT had been switched in the VCF file for unknown reason. For example, it should be G at locus 1234 in RefSeq, and the variant is T. But the VCF file records T(REF) and G(ALT).

I only have the VCF file, and do not have the original PED file. Is there any tool or method to check if the REF alleles are correct by RefSeq and switch REF and ALT columns (or just remove this loci) in the VCF if they're wrong?

Thanks!

vcf ref allele alt allele plink
• 11,957 views • link
updated by Shicheng Guo Shicheng Guo • written by miaowzai miaowzai

I guess the solution should be like this way:

  1. check the reference allele with genome reference (hg19.fa or hg38.fa)

  2. if the reference allele is not same with the human genome, then check whether the alternative allele is same with human reference

  3. if the alternative allele is same with human genome, then switch it. For the sample gentoype, 0 -> 1 and 1 -> 0

  4. if neither reference and alternative allele is not same as reference, then .... remove it??

• 0 views • link
updated by Ram Ram • written by Shicheng Guo Shicheng Guo

2 answers

Given a file with the correct reference alleles for each variant ID, you can use plink's --a2-allele flag to fix them in the VCF.

• 0 views • link
written by chrchang523 chrchang523

Thanks for the quick answer. I never used plink. I looked at the manual of --a2-allele in plink documentation and I don't completely understand. Could you elaborate the method? Thank you!

• 0 views • link
written by miaowzai miaowzai

It depends on what your RefSeq file looks like. The basic structure of the command would be

plink --vcf [name of plink-exported VCF with incorrect reference alleles] --a2-allele [name of RefSeq file] [1-based column index of ref alleles] [1-based column index of variant IDs] --recode vcf --real-ref-alleles
• 0 views • link
written by chrchang523 chrchang523

Now I get the best solution with bcftools

bcftools view -T ^remove.txt $i -Ov | bcftools norm -m-any -Ov --check-ref w -f ~/hpc/db/hg19/hg19.fa | bcftools annotate -x ID,INFO,FORMAT -I +'%CHROM:%POS' -Oz -o $i.trim.vcf.gz
• 0 views • link
written by Shicheng Guo Shicheng Guo

That's bcftools, not vcftools. They're two different tools.

• 0 views • link
written by Ram Ram

Log in to answer this question.

More posts like this