PLINK clump behavior on missing SNPs?

Question

I have a long list of autoimmune-associated SNPs, and I want to boil it down so that I get one SNP representing each LD block. I chose to use PLINK's --clump option for this. I'm roughly following this tutorial (but analyzing my own data).

I don't have the full data from the original GWAS's, so I am using 1000 Genomes data as the reference for LD estimation. But, the 1000 Genomes data don't have all the SNPs I need. I end up a few hundred warnings:

Warning: 'rs12345678' is missing from the main dataset, and is a top variant.

How does PLINK handle this situation? More importantly, how should I handle this situation? I don't want to omit those SNPs, but I also don't want to take them all, for the same reason that I am doing this in the first place.

Answer 1

Most likely, the GWASs that generated your summary statistics used other imputation panels than 1000G, like HRC. Clearly, PLINK can't estimate the LD for SNPs that aren't found in the reference, and most likely just ignores them. I think the easiest solution is to use whatever individual-level dataset you have to estimate the correlations - even if it's small or from a different population. The estimate will be biased, but you don't really need high precision for clumping.

Then again, you can skip PLINK and try searching for SNP prioritization methods based on summary statistics - it's an area of active development, and something might be available to suit your needs.

Answer 2

You can use imputation for guessing at what the missing SNPs might be based on known LD patterns in populations. This procedure will give you an idea of whether the recombinational history of genomic regions is important, rather than specific SNPs. Given that you're looking at LD-based measures, this should be okay (although the logic behind the calculation is a bit circular).

If you don't have a big compute cluster available, the Michigan Imputation Server is a free imputation service that can be used for academic purposes. Data is fed into the service as bgzip-compressed VCF files (output from plink as an uncompressed VCF file, then use bgzip to compress the files), and returned as a larger dataset with imputed SNPs derived from the HRC or 1kg reference panels. There's a bit of initial legwork required: you need to make sure that the variants as defined in the VCF file have the same location and strand as defined in the reference panel, otherwise you'll wait weeks in the queue and end up empty handed. To avoid hassle, I'd recommend removing any complementary SNPs (e.g. A/T, C/G), then flipping the strand for any remaining genotypes where the alleles don't match up.

I've recently done this again. Here's a summary of the things required to get MIS imputation working properly from a PLINK dataset:

1. Retrieve site data for HRC panel
2. Generate list of complementary variants, rename chrX to chr23 to fit PLINK expectation
3. extract out variant names (remove duplicated rs numbers)
4. change 'RS' to 'rs' in plink BIM file
5. filter PLINK file against variants
6. adjust position of variants to fit with HRC panel
7. sort HRC variants by rs number
8. identify variants that need flipping
9. flip variants using PLINK, create combined VCF file
10. create VCF index for bcftools
11. extract per-chromosome VCF files
12. convert chromosome '23' back to chromosome 'X'
13. upload to MIS, remembering to separate out autosomal and X-chromosomal files

And a full example:

## Retrieve site data for HRC panel
wget 'ftp://ngs.sanger.ac.uk/production/hrc/HRC.r1-1/HRC.r1-1.GRCh37.wgs.mac5.sites.vcf.gz'
## Generate list of complementary variants, rename chrX to chr23 to fit PLINK expectation
zgrep -v '^#' HRC.r1-1.GRCh37.wgs.mac5.sites.vcf.gz | cut -f 1-5 | \
  awk '{if(!($4 ~ /[AT]/ && $5 ~ /[AT]/) && !($4 ~ /[CG]/ && $5 ~ /[CG]/)){print}}' | perl -pe 's/^X/23/' > HRC_complementary_variants.tsv
## extract out variant names (remove duplicated rs numbers)
cut -f 3 HRC_complementary_variants.tsv | grep '^rs' | uniq -u > HRC_complementary_variants_namesOnly.txt
## change 'RS' to 'rs' in plink BIM file
perl -i -pe 's/RS/rs/' plinkFileHg19.bim
## filter PLINK file against variants
plink --bfile plinkFileHg19 --extract HRC_complementary_variants_namesOnly.txt --make-bed --out noComp_plinkFileHg19
## adjust position of variants to fit with HRC panel
plink --bfile noComp_plinkFileHg19 --update-chr HRC_complementary_variants.tsv 1 3 \
  --update-map HRC_complementary_variants.tsv 2 3 --make-bed --out repositioned_nc_plinkFileHg19
## sort HRC variants by rs number
sort -k 3,3 HRC_complementary_variants.tsv > rsSorted_HRC_complementary_variants.tsv
## identify variants that need flipping
join -1 3 -2 2 <(pv rsSorted_HRC_complementary_variants.tsv) <(sort -k 2,2 repositioned_nc_plinkFileHg19.bim) | \
  awk '{if(!($4 == $9) && !($4 == $8)){print $1}}' > toFlip_HRC_complementary_variants.txt
## flip variants using PLINK, create combined VCF file
plink --bfile repositioned_nc_plinkFileHg19 --flip toFlip_HRC_complementary_variants.txt --recode vcf bgz --out flipped_rnc_plinkFileHg19
## create VCF index for bcftools
bcftools index flipped_rnc_plinkFileHg19.vcf.gz
## extract per-chromosome VCF files
for x in $(seq 1 23); do echo ${x};
  bcftools filter -r ${x} flipped_rnc_plinkFileHg19.vcf.gz -O z -o chr${x}_flipped_rnc_plinkFileHg19.vcf.gz;
done
## convert chromosome '23' back to chromosome 'X'
zcat chr23_flipped_rnc_plinkFileHg19.vcf.gz | perl -pe 's/ID=23/ID=X/;s/^23\t/X\t/' > chrX_flipped_rnc_plinkFileHg19.vcf
bgzip chrX_flipped_rnc_plinkFileHg19.vcf