6
$\begingroup$

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.

$\endgroup$
4
  • 1
    $\begingroup$ Are you sure they're missing? I mean, are you sure it isn't just a question of wrong/outdated/changed rsID? I've never done anything like this, but isn't there a way of using actual SNPs (genomic coordinates, reference allele and alternate allele) somehow? Perhaps by feeding it a vcf file? $\endgroup$
    – terdon
    Oct 31 '17 at 18:10
  • $\begingroup$ That's a good point. I didn't realize rsID's could go stale like that. $\endgroup$ Oct 31 '17 at 18:34
  • $\begingroup$ @terdon On the contrary. rsIDs are designed to be a stable alternative to just using position and allele, so that when reference assembly changes, SNPs can be tracked by rsIDs. $\endgroup$
    – juod
    Nov 5 '17 at 11:50
  • 1
    $\begingroup$ @juod yes but, alas, they can change. Certainly in cases where rsIDs are merged for example. In those cases, one of the merged IDs is rendered obsolete. $\endgroup$
    – terdon
    Nov 5 '17 at 12:31
5
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Unfortunately, I do not have even an individual-level dataset. But, I am accepting this answer because so many SNPs are missing, and I doubt different dbSNP versions can account for so many missing SNPs. $\endgroup$ Nov 5 '17 at 21:07
3
$\begingroup$

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
$\endgroup$
3
  • $\begingroup$ fyi, plink can directly generate bgzipped VCF files (--recode vcf bgz), and it'll use multiple threads to do the compression. $\endgroup$ Jan 17 '18 at 19:59
  • $\begingroup$ Yes, but the plink-generated bgz-VCF files didn't work properly with the Michigan Imputation Server last time I tried, and it's a real pain testing that out due to the queue length. $\endgroup$
    – gringer
    Jan 17 '18 at 21:14
  • $\begingroup$ Thanks for mentioning this, I'll investigate. I noticed that the recommended pipelines now include a "bcftools sort" step after VCF file generation even though that VCF would be position-sorted, so there may be a poorly documented assumption about how multiple variants sharing the same bp coordinate must be sorted, and/or an assumption about chromosome order. $\endgroup$ Jan 18 '18 at 18:31

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.