How to Calculate Genomic Inflation Factor and λgc for GWAS

You have conducted your genome-wide association study (GWAS) and have tested each genetic variant for an association with your trait of interest. Now it is time to investigate if there are any systematic biases that may be present in your association results. A common way to do this is to calculate the genomic inflation factor, also known as lambda gc (λgc). By definition, λgc is defined as the median of the resulting chi-squared test statistics divided by the expected median of the chi-squared distribution. The median of a chi-squared distribution with one degree of freedom is 0.4549364. A λgc value can be calculated from z-scores, chi-square statistics, or p-values, depending on the output you have from the association analysis. Follow these simple steps to calculate lambda GC using R programming language.

(1) Convert your output to chi-squared values
(2) Calculate lambda gc (λgc)
If analysis results your data follows the normal chi-squared distribution, the expected λgc value is 1. If the λgc value is greater than 1, then this may be evidence for some systematic bias that needs to be corrected in your analysis.

Is there a GWAS on that?

A great online resource to find whether a genome-wide association study (GWAS) has been published on a certain trait or disease is the National Human Genome Research Institute (NHGRI) webpage where they maintain a curated catalog of published GWAS.  Here trained curators are constantly scanning PubMed publications and other genomic resources looking for association studies linking a genomic position (usually a tagging SNP) to a disease trait of interest.  Details listed include study size, population, locus, risk allele, odds ratio, p-value, and other pertinent statistics.  Recently, the NHGRI in collaboration with the European Bioinformatics Institute (EBI) released an interactive version of the GWAS catalog called the GWAS Diagram Browser.  This provides a great way to visualize and filter many of the genome-wide significant findings from genome-wide associations studies.  Highlights include filtering by disease, time series views, and some useful downloads.

Another noteworthy resource is the HuGE Navigator GWAS Interagator.  This is a search tool similar to the NHGRI GWAS catalog, but more focused on a search terms.  Handy links are provided to other resources.  Of particular interest are links to visualize the variants in the UCSC Browser.

I am sure other GWAS resources exist as well, but these were the two main ones that first came to mind for me.  If you know other great GWAS resources capable of linking a genomic marker with a disease, please share in the comments below.

Estimate Combined Percent Familial Risk Explained by GWAS Loci

Most current genome-wide association studies (GWAS) include a calculation of the percent familial risk the discovered loci explain.  This calculation indicates how much of the familial risk can be accounted for by the known loci and is usually used as evidence there are additional yet undetected loci that remain to be discovered.  Looking through references, it can be a bit difficult to find exactly how this calculation is performed.  One reference I found that includes a formula for the calculation is by Cox et al. 2007 (PMID:17293864), but I am sure there are plenty others that also include a formula.   While the Cox et al. formula for calculating familial relative risk due to each locus is arranged differently than the ones below, the two are equivalent.  I just find this arrangement less cumbersome to use.  The overall equation is to compare the cumulative risk of the known loci (sum of log lambda k) to the estimated risk of a first degree relative (log lambda 0).

where:

      p is the risk allele frequency for locus k

      r is the per allele odds ratio for locus k.

To make calculations easy, I made a simple R script that does all the calculations automatically.  The input for the script is a file with 3 columns:

(1) Annotation for the SNP - this can be anything, for example: RS number, chromosomal coordinates, etc.
(2) Risk allele frequency - this is the frequency of the risk allele (range: 0-1) equal to p in the above equation.
(3) Per allele odds ratio - odds ratio for every one unit increase in the number of risk alleles.

Note, the risk allele frequency is the frequency of the risk allele and not the minor allele frequency.  The program also needs an estimate of the familial relative risk (lambda 0).  This can usually be done by looking for previous familial studies for the disease.

Here is the R script:

It can be run from the command line by the example command:

Rscript familial_risk_snps.R snp_lst.txt 4

where:
      familial_risk_snps.R is the name of the script.
      snp_lst.txt is the input file with three columns described above.
      4 is the estimate of the familial relative risk of the disease.

Create Triangle Plot from Inferred Genetic Ancestry

I previously posted on how to infer ancestry for a group of study participants using SNP genotypes.  Today, I want to visually plot some of the output in R.  Two informative plots that can be generated from the output are a standard plot of the two eigenvectors with percent ancestry overlaid and a triangle (or ternary) plot with each axis representing percentage of one of the three ancestral populations (ex: European, African, and Asian).

Here is some simple code to plot this in R.  There is no base package to plot the triangle plot, so the plotrix package will need to first be installed.  The ancestry.txt file is the output file from SNPWEIGHTS, but other output could be formatted to work as well.


The output should look similar to the plots below.

Find SNP Overlap from Different SNP Arrays

If you have .bed files for array manifests it is relatively easy to use UNIX to compare the overlap between SNP positions in the two files.  Here is a quick one line of code to do so:


This essentially finds all rows that are the same in the two files.  It requires all unique rows in each file and will only work for one base pair long items.

How to Infer Ancestry from SNP Genotypes

Self-reported ancestry is poor metric to use when attempting to statistically adjust for the effects of ancestry since several individuals falsely report their ancestry or are simply unaware of their true ancestry.  Worse yet, sometimes you don't even have information collected on an individual's ancestry.  As many of you know, if you have SNP genotyping data you can rather precisely estimate the ancestry of an individual.  Classically, to do this one needed to combine genotypes from their study sample with genotypes from a reference panel (eg: HapMap or 1000 Genomes), find the intersection of SNPs in each dataset, and then run a clustering program to see which samples clustered with the reference ancestral populations.  Not a ton of work, but a minor annoyance at best.  Luckily, a relatively new program was just released that, in essence, does a lot of this ground work for you.  It is called SNPWEIGHTS and can be downloaded here.  Essentially, the program takes SNP genotypes as input, finds the intersection of the sample genotypes with reference genotypes, weights them based on pre-configured parameters to construct the first couple of principle components (aka. eigenvectors) and then calculates an individual's percentage ancestry for each ancestral population.  Here is how to run the program.

First, make sure you have Python installed on your system and that your genotyping data is in EIGENSTRAT format.  A brief tutorial to convert to EIGENSTRAT format using the convertf tool is here.

Next, download the SNPWEIGHTS package here and a reference panel. I usually use the European, West African and East Asian ancestral populations, but there are other options on the SNPWEIGHTS webpage as well.

Then, create a parameter file with directories of input files, your input population (designated "AA", "CO", and "EA"), and a output file.  An example of one is below:


Finally, run the program using the command inferancestry.py --par par.SNPWEIGHTS.  For the program to run correctly, make sure the inferancestry.info and snpwt.co files are in the same directory as your inferancestry.py file.

For more details, see the SNPWEIGHTS paper or the README file included in the SNPWEIGHTS zip folder.  For code on generating eigenvector plots with overlaid ancestry percentages and triangle plots according to percent ancestry, see this post.

What Are SNP Ambiguity Codes and What Do They Mean?

That's a good question.  In fact one that I had myself.  Here's what I found:

Apparently single nucleotide polymorphism (SNP) ambiguity codes were constructed by the International Union of Pure and Applied Chemistry (IUPAC) to denote nucleotide changes in SNPs.  Here is a table of the meaning of each code taken from the ENSEMBLE SNPView website.

IUPAC Code	Mnemonic	Meaning	Complement
A	Adenine	A	T
C	Cytosine	C	G
G	Guanine	G	C
T/U	Thymidine	T	A
K	Keto	G or T	M
M	Amino	A or C	K
S	Strong	C or G	S
W	Weak	A or T	W
R	Purine	A or G	Y
Y	Pyrimidine	C or T	R
B	not A	C, G, or T	V
D	not C	A, G, or T	H
H	not G	A, C, or T	D
V	not T or U	A, C, or G	B
N	any	G, A, T, or C	N

How to get Phased Imputed Genotypes for Haplotype Analysis

I wanted to compare haplotypes from a set of cases genotyped on one platform to a set of reference controls genotyped on another platform for a small chromosomal region.  I looked into a few methods to do this and found that phasing the genotypes using SHAPEIT and then imputing off the 1000 Genomes reference panel using IMPUTE2 was the optimal approach to have a set of overlapping genotypes for a haplotype analysis.  Phasing from the .ped and .map files was easy using SHAPEIT.  The command I used with my .ped and .map files was:


Other file types can also be used as direct input to SHAPEIT (ex: .bed/.bim/.fam and .gen/.sample).  Genetic recombination maps can be downloaded here for build36 and build37.  Also the backslashes (\) are not necessary.  They just help organize the code by telling the computer to keep reading the input on the next line.

Once you have your .haps and .sample file generated, you are then ready to impute using IMPUTE2.  I used the command below:


Again the backslashes are just for a cleaner visualization of the code.  It is important when using the 1000 Genomes as your reference panel to have the SNP coordinates in hg19 coordinates.  LiftOver can help you convert from one build to another.  If you need the 1000 Genomes reference panel, it can be downloaded here.  The genetic map file is from the above SHAPEIT download.  Finally, if you want phased imputed results, the -phase command is essential to include.

You should now have your phased imputed genotypes to do whatever you wish with.  I chose to convert the phased imputed genotypes into a .fasta file with entries for each haplotype (see Python script).  Then I created parsimony trees with MEGA and visualized them in HapView.

Download One Thousand Genome Data for Haploview

Haploview has a built-in portal to download HapMap data, but Haploview development hasn't kept pace with developing a way to automatically download 1000G SNP data.  Searching for a way to visualize the higher density SNP coverage of the 1000G project, I found it was not all too difficult to do this manually.  It involves a couple of extra steps.

First, determine the genomic coordinates of region you are interested in.  This needs to be in hg19 coordinates.  If you have hg18 coordinates, liftOver is a useful tool to convert coordinates from one human genome build to another (liftOver format is chr:start-end, for example: chr8:1000-50000).

Next, go to this 1000G website, and plug in your genomic coordinates of interest.  Here the coordinates should not include chr in the chromosome name (for example: 8:1000-50000).  Then on the next page select ancestral populations you are interested in (you can select multiple populations by holding down Ctrl).  Give the website a few seconds to generate the files.  Eventually a link to a marker information file (.info) and linkage pedigree file (.ped) will appear.  Right click on each of these files and save them to your computer.

Now, fire up Haploview and select Open new data.  Go to the Linkage Format tab and browse for your .ped file in the Data File field and your .info file in the Locus Information File field (the .info file field is usually automatically generated after selecting your .ped file if your .ped and .info files have the same prefix).  Haploview will load the files and you should be ready to visualize the LD structure.  Enjoy using 1000G data in Haploview!