Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-01-17T19:59:36.020Z
Answers.Text: fyi, plink can directly generate bgzipped VCF files (--recode vcf bgz), and it'll use multiple threads to do the compression.
Answers.UserDisplayName: Christopher Chang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-01-17T21:14:42.387Z
Answers.Text: 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.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-01-18T18:31:13.420Z
Answers.Text: 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.
Answers.UserDisplayName: Christopher Chang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-01-12T01:55:04.370Z
Answers.LastActivityDate: 2018-02-01T04:07:39.073Z
Answers.LastEditDate: 2018-02-01T04:07:39.073Z
Answers.UserDisplayName: gringer
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-05T21:07:57.100Z
Answers.Text: 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.
Answers.UserDisplayName: eric_kernfeld
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-05T12:01:23.483Z
Answers.LastActivityDate: 2017-11-05T12:01:23.483Z
Answers.UserDisplayName: juod
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-31T18:10:33.753Z
Comments.Text: 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?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-10-31T18:34:09.510Z
Comments.Text: That's a good point. I didn't realize rsID's could go stale like that.
Comments.UserDisplayName: eric_kernfeld
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-11-05T11:50:00.780Z
Comments.Text: @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.
Comments.UserDisplayName: juod
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-11-05T12:31:34.960Z
Comments.Text: @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.
Comments.UserDisplayName: terdon
Tags: <snp><gwas><plink>
Title: PLINK clump behavior on missing SNPs?
ThreadText: # 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**: @username_1 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**: @username_2 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](https://imputationserver.sph.umich.edu/index.html) 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
  
  ```
  
  > **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.  
  
  ---
  
  