<p>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).</p>

<p>If you don't have a big compute cluster available, the <a href="https://imputationserver.sph.umich.edu/index.html" rel="nofollow noreferrer">Michigan Imputation Server</a> 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 <code>plink</code> as an uncompressed VCF file, then use <code>bgzip</code> 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.</p>

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

<ol>
<li>Retrieve site data for HRC panel</li>
<li>Generate list of complementary variants, rename chrX to chr23 to fit PLINK expectation</li>
<li>extract out variant names (remove duplicated rs numbers)</li>
<li>change 'RS' to 'rs' in plink BIM file</li>
<li>filter PLINK file against variants</li>
<li>adjust position of variants to fit with HRC panel</li>
<li>sort HRC variants by rs number</li>
<li>identify variants that need flipping</li>
<li>flip variants using PLINK, create combined VCF file</li>
<li>create VCF index for bcftools</li>
<li>extract per-chromosome VCF files</li>
<li>convert chromosome '23' back to chromosome 'X'</li>
<li>upload to MIS, remembering to separate out autosomal and X-chromosomal files</li>
</ol>

<p>And a full example:</p>

<pre><code>## 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]/ &amp;&amp; $5 ~ /[AT]/) &amp;&amp; !($4 ~ /[CG]/ &amp;&amp; $5 ~ /[CG]/)){print}}' | perl -pe 's/^X/23/' &gt; HRC_complementary_variants.tsv
## extract out variant names (remove duplicated rs numbers)
cut -f 3 HRC_complementary_variants.tsv | grep '^rs' | uniq -u &gt; 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 &gt; rsSorted_HRC_complementary_variants.tsv
## identify variants that need flipping
join -1 3 -2 2 &lt;(pv rsSorted_HRC_complementary_variants.tsv) &lt;(sort -k 2,2 repositioned_nc_plinkFileHg19.bim) | \
  awk '{if(!($4 == $9) &amp;&amp; !($4 == $8)){print $1}}' &gt; 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/' &gt; chrX_flipped_rnc_plinkFileHg19.vcf
bgzip chrX_flipped_rnc_plinkFileHg19.vcf
</code></pre>

CC BY-SA 3.0
2018-01-17T19:59:36.020Z
fyi, plink can directly generate bgzipped VCF files (--recode vcf bgz), and it'll use multiple threads to do the compression.
Christopher Chang
CC BY-SA 3.0
2018-01-17T21:14:42.387Z
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.
gringer
CC BY-SA 3.0
2018-01-18T18:31:13.420Z
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.
Christopher Chang
CC BY-SA 3.0
2018-01-12T01:55:04.370Z
2018-02-01T04:07:39.073Z
2018-02-01T04:07:39.073Z
gringer
<p>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.</p>

<p>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.</p>

CC BY-SA 3.0
2017-11-05T21:07:57.100Z
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.
eric_kernfeld
CC BY-SA 3.0
2017-11-05T12:01:23.483Z
2017-11-05T12:01:23.483Z
juod
<p>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 <a href="https://neurogenetics.qimrberghofer.edu.au/iSECA/LD_clumping_tutorial.html" rel="noreferrer">this</a> tutorial (but analyzing my own data).</p>

<p>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:</p>

<p><code>Warning: 'rs12345678' is missing from the main dataset, and is a top variant.</code></p>

<p>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.</p>

CC BY-SA 3.0
2017-10-31T18:10:33.753Z
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?
terdon
CC BY-SA 3.0
2017-10-31T18:34:09.510Z
That's a good point. I didn't realize rsID's could go stale like that.
eric_kernfeld
CC BY-SA 3.0
2017-11-05T11:50:00.780Z
@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.
juod
CC BY-SA 3.0
2017-11-05T12:31:34.960Z
@juod yes but, alas, they can change. Certainly in cases where rsIDs are [merged](https://www.ncbi.nlm.nih.gov/books/NBK44467/#Build.do_you_have_tables_that_show_those) for example. In those cases, one of the merged IDs is rendered obsolete.
terdon
<snp><gwas><plink>
PLINK clump behavior on missing SNPs?
# PLINK clump behavior on missing SNPs?

(Asked by: username_0 on 2017-10-31)

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](https://neurogenetics.qimrberghofer.edu.au/iSECA/LD_clumping_tutorial.html) 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.

**Tags:** snp, gwas, plink
> **Comments:**
> **username_1**: 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?  
> **username_0**: That's a good point. I didn't realize rsID's could go stale like that.  
> **username_2**: @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.  
> **username_1**: @juod yes but, alas, they can change. Certainly in cases where rsIDs are [merged](https://www.ncbi.nlm.nih.gov/books/NBK44467/#Build.do_you_have_tables_that_show_those) for example. In those cases, one of the merged IDs is rendered obsolete.  

---

## 2 answers

### Answer by username_2 on 2017-11-05. Score: 6 | Accepted answer

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.

> **Comments:**
> **username_0**: 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.  

---

### Answer by username_3 on 2018-01-12. Score: 3

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:

Retrieve site data for HRC panel
Generate list of complementary variants, rename chrX to chr23 to fit PLINK expectation
extract out variant names (remove duplicated rs numbers)
change 'RS' to 'rs' in plink BIM file
filter PLINK file against variants
adjust position of variants to fit with HRC panel
sort HRC variants by rs number
identify variants that need flipping
flip variants using PLINK, create combined VCF file
create VCF index for bcftools
extract per-chromosome VCF files
convert chromosome '23' back to chromosome 'X'
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

> **Comments:**
> **username_4**: fyi, plink can directly generate bgzipped VCF files (--recode vcf bgz), and it'll use multiple threads to do the compression.  
> **username_3**: 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.  
> **username_4**: 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.  

---

