How to extract proxy SNPs using 1000G vcf for a list of SNPs

Question

[janhuang.cn]

I would like to find the proxy SNPs (r2>0.8) for a list of SNPs. I have a long list, so I do not want to manually check and download from LDlink. I think plink and/or vcftools can do the job. But I came across some problems.

I download the 1000G phase3 vcf (gzziped) from the ftp folder: ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/

Let say my list has three SNPs: rs157595 (chr19), rs157594 (chr19), rs7557280 (chr2). This list is saved in one txt file (SNPneedProxy.txt), with one column.

My expected output will be three txt files, each file contains the proxy SNPs of each provided SNP and the r2.

1) How can I read the .gz file directly in plink?

If I unzip the file, the below script using --vcf can load the file.

plink --vcf ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf

If I do NOT unzip the file, the below script using --data --gen does not work.

plink --data --gen ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf
Error: Failed to open plink.sample.

2) I expect the below script to return a file with all proxy SNPs, but it does not. Did I misunderstand the function of this line?

plink --vcf ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf --r2 --ld-window-r2 0.8 --ld-snp-list SNPneedProxy

3) If I use vcftools, the below script neither return a list of proxy SNPs.

vcftools --gzvcf ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz --remove-indels --snps SNPneedProxy.txt --hap-r2 --min-r2 0.8 --out results

I wonder if I use the wrong function and can anyone suggests a way to do this? Thank you.

Answer 1

[chrchang523]

1. Just use —vcf on the gzipped vcf. Why are you trying to use —data/—gen, which are for a totally different file format?!

[janhuang.cn]

Error: Failed to open
ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.

I got an error when using --vcf for gzipped vcf

[chrchang523]

The error message shows that you didn’t type the right filename (no .gz at the end). Try again.

[janhuang.cn]

But I got the similar error message no matter I put .gz or not.

plink --vcf ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz --r2 --ld-window-r2 0.8 --ld-snp rs157595

Error: Failed to open
ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.

plink --vcf ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf --r2 --ld-window-r2 0.8 --ld-snp rs157595

Error: Failed to open
ALL.chr19.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.

[chrchang523]

Check where your file is located and what its exact name is. "Failed to open" means you got the filename wrong.