0
$\begingroup$

I am working with the public WGS files from the HGDP from here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/

I want to convert the VCFs to plink binary format (bed, bim, fam). However, there is a strange error:

user@desktop:/media/luks8tb1/data/genomics/HGDP$ plink --vcf hgdp_chr21_chrpos.vcf_hg37.vcf --recode --snps-only 'just-acgt' --biallelic-only --allow-no-sex --make-bed --out hgdp_chr21_chrpos_chrpos_hg37
PLINK v1.90b6.17 64-bit (28 Apr 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to hgdp_chr21_chrpos_chrpos_hg37.log.
Options in effect:
  --allow-no-sex
  --biallelic-only
  --make-bed
  --out hgdp_chr21_chrpos_chrpos_hg37
  --recode
  --snps-only just-acgt
  --vcf hgdp_chr21_chrpos.vcf_hg37.vcf

64259 MB RAM detected; reserving 32129 MB for main workspace.
--vcf: 35k variants complete.
Error: ALT allele duplicates REF allele on line 37074 of .vcf file.

I.e., it gets 35k variants into the VCF before it meets a single invalid line, and then gives up. Is there no way to simply skip these lines? I read through the plink1.9 docs but didn't see any way to deal with this problem.

The conversion to plink binary worked fine on the original files, this allele issue was caused by my use of CrossMap to move the data to hg37 from the original hg38 (see this question, I used it on the VCFs).

This 2016 biostars thread has the same issue. One person suggests using:

cat infile.vcf | grep '^#' > infile_filtered.vcf
cat infile.vcf | grep -v '^#' | awk '{if ($4 != $5) print $0}' >> infile_filtered.vcf

I.e., the first prints the comments in the top of the file, and the second prints all the non-comment rows where the REF and ALT columns aren't identical. Another person gives an awk solution:

awk '($0 ~/^#/ || $4!=$5)' input.vcf > out.vcf

Is this really the best solution?

In my case, I tried their code, but it does not work. Looking at the offending line in the VCF, I see:

user@desktop:/media/luks8tb1/data/genomics/HGDP$ sed -n 37074p hgdp_chr21_chrpos_hg37_filtered.vcf
chr21   10370717    rs183658892 T   T,A 131621  PASS

So that's why the sed/awk solutions don't work. The REF allele is listed among others in the ALT, it's not the only one.

$\endgroup$
2
  • 1
    $\begingroup$ Rather than just excluding it, I'd probably try and investigate why your liftover tool has maybe bugged out, given you the same allele in the ref and alt fields. What happens if you look at the same position on your data before lifting it over, does the vcf field look like then? $\endgroup$
    – user438383
    Aug 28 '20 at 7:26
  • $\begingroup$ Good question. It looks like this chr21 9892675 rs183658892 C T,A so no REF=ALT issue, but it has mismatching alleles. I ended up extracting all the allele metadata in the lifted VCF and finding all the problem variants by loading the metadata into memory in R. There were only 6 REF=ALT problems, all in multi-allele variants. I didn't want these anyway, so I removed them with bcftools view -k -m2 -M2 -v snps before lifting and conversion to plink. I will post my answer in a bit when I am sure no more issues. $\endgroup$ Aug 28 '20 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.