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.
