Practical One
Quality control in genome-wide association studies
In this practical, we will go through the steps in performing quality control (QC) of genotype data from a simulated genome-wide association study of 1000 cases and 1000 controls, typed for 317,503 autosomal and X chromosome SNPs. We will begin by performing sample QC, including calculation of call rates, heterozygosity, and sex discordance. We will then perform SNP QC, including calculation of call rates and deviation from Hardy-Weinberg equilibrium.
This practical is based on “Data quality control in genetic case-control association studies” (Anderson et al. 2010, Nature Protocols 5: 1564-73).
In order to run this practical, you will require the following resources:
- Computing workstation with Unix or Linux operating system;
- PLINK software (
- Statistical software for graphical presentation of results, such as R (
- Genome-wide association binary ped files and analysis script files (which can be downloaded from
Before starting the practical, you will need to unpack the genome-wide association binary ped files and additional analysis script files. You can do this by typing the following commad at the shell prompt in the directory where you downloaded the file “raw-GWA-data.tgz”:
tar -xvzfraw-GWA-data.tgz
Part A: Sample QC
Identification of individuals with discordant sex information
At the shell prompt, type:
plink --bfile raw-GWA-data --check-sex --out raw-GWA-data
This command will calculate the mean homozygosity rate across X-chromosome markers for each individual in the study. We can produce a list of individuals with discordant sex data by typing:
grep PROBLEM raw-GWA-data.sexcheck > raw-GWA-data.sexprobs
Open the file “raw-GWA-data.sexprobs” to obtain the family IDs (column 1) and individual ID (column 2) for these individuals. Column 3 denotes ascertained sex and column 4 denotes sex according to genotype data. When the homozygosity rate is more than 0.2 but less than 0.8, the genotype data are inconclusive regarding the sex of an individual and these are marked in column 4 by “0”.
In general, any discordancies in sex should be reported to study co-ordinators to double check records for errors. In situations in which discrepancy cannot be resolved, add the family ID (FID) and individual ID (IID) of the samples to a file named “fail-sexcheck-qc.txt” (one individual per line, tab delimited). This file can be usedto exclude these samples from downstream analyses.
Identification of individuals with elevated missing data rates or outlying heterozygosity rate
At the shell prompt, type:
plink --bfile raw-GWA-data --missing --out raw-GWA-data
This command will create the files “raw-GWA-data.imiss” and “raw-GWA-data.lmiss”. The fourth column in the file “raw-GWA-data.imiss” (N_MISS) denotes the number of
missing SNPs and the sixth column (F_MISS) denotes the proportion of missing SNPs per individual.
At the shell prompt type:
plink --bfile raw-GWA-data --het --out raw-GWA-data
This command will create the file “raw-GWA-data.het”, in which the third column denotes the observed number of homozygous genotypes [O(Hom)] and the fifth column denotes the number of non-missing genotypes [N(NM)] per individual.
You can alculate the observed heterozygosity rate per individual using the formula:
Het = (N(NM) − O(Hom))/N(NM)
Create a graph in which the observed heterozygosity rate per individual is plotted on the x axis and the proportion of missing SNPs per individuals is plotted on the y axis. This can be carried out using standard software such as Excel or R. A script for calculating the heterozygosity rate and producing the graph using R is supplied: “imiss-vs-het.Rscript”. At the shell prompt, type:
R CMD BATCH imiss-vs-het.Rscript
This will create the graph “raw-GWA-data.imiss-vs-het.pdf”. Examine the plot to decide reasonable thresholds at which to exclude individuals based on elevated missing genotype rates or extreme heterozygosity. The red lines indicate suggested thresholds: missing data rate of 3% and heterozygosity rate ± 3 SD from the mean. This file can easily be adapted to highlight different thresholds. Add the FID and IID of the samples failing this QC to a file named “fail-imisshet-qc.txt” (one individual per line, tab delimited). This file can be usedto exclude these samples from downstream analyses.
Identification of duplicated or related individuals
To minimize computational complexity, create an “independent” set of SNPs to generate theidentity by descent (IBS) matrix. This can be done by “pruning” the data so that no pair of SNPs (within a given genomic interval) has an r2 value greater than a given threshold, typically chosen to be 0.2. Typically, we would also exclude regions of strong LD, such as the MHC, which are listed in the file “high-LD-regions.txt”.
Please note that the data for this practical have been simulated without LD between SNPs. We have created a “pruned” list of SNPs to use called “raw-GWA-data.prune.in” which you should use for this practical. However, in general, you would use the following command to generate the pruned SNP list:
plink --file raw-GWA-data --exclude high-LD-regions.txt --range \
--indep-pairwise50 5 0.2 --out raw-GWA-data
To generate IBS between each pair of individuals, type the following command at the shell prompt:
plink --bfile raw-GWA-data --extract raw-GWA-data.prune.in --genome--out raw-GWA-data
Please note that this step is computationally demanding and will take some time to run.
We have created a script file to identify all pairs of individuals with IBS > 0.185 which outputs the ID of the individual from the pair with lowest call rate (from the previously created file “raw-GWA-data.imiss”). This script file can easily be adapted to other data sets and thresholds for “related” individuals.
To run the script file, type the following command at the shell prompt:
perl run-IBD-QC.pl raw-GWA-data
The command creates a file called “fail-IBD-QC.txt” which can be used to exclude these samples from downstream analyses.
Removal of all individuals failing sample QC
At the shell prompt, type the following command to concatenate all files listing individuals who have failed previous QC steps:
cat fail-* | sort -k1 | uniq > fail-qc-inds.txt
The file “fail-qc-inds.txt” should now contain a list of unique individuals failing the previous QC steps. To remove them from the data set, type the following command at the shell prompt:
plink --bfile raw-GWA-data --remove fail-qc-inds.txt --make-bed --out clean-inds-GWA-data
The binary ped file “clean-inds-GWA-data” can be used for subsequent SNP QC analyses.
Part B: SNP QC
Identification of all SNPs with an excessive missing data rate
To calculate the missing genotype rate for each SNP, type the following command at the shell prompt:
plink --bfile clean-inds-GWA-data--missing --out clean-inds-GWA-data
The third column in the file “clean-inds-GWA-data.lmiss”(N_MISS) denotes the number of
missing genotypes and the fifth column (F_MISS) denotes the proportion of missing genotypes per SNP.
A script for plotting a histogram of the missing genotype rateusing R is supplied: “lmiss-hist.Rscript”. At the shell prompt, type:
R CMD BATCH lmiss-hist.Rscript
This will create the graph “clean-inds-GWA-data.lmiss.pdf”. Examine the plot to decide a reasonable threshold at which to exclude SNPs based on elevated missing data. The dotted line indicates a suggested threshold of a missing data rate of 5%. This file can easily be adapted to highlight different thresholds.
To remove SNPs with call rate less than 5%, simply add the “--geno 0.05” option to the PLINK command line. We will do this below when creating our final cleaned data set.
Test SNPs for different genotype call rates between cases and controls
To test for differential call rates between cases and controls for each SNP, type the following command at the shell prompt:
plink --bfile clean-inds-GWA-data --test-missing --out clean-inds-GWA-data
The output of this test can be found in the file “clean-inds-GWA-data.missing”. We have created a script to highlight all SNPs with significant differences in case and control call rates (p<10-5) from this output file. This script can be easily adapted to other data sets and thresholds for differential call rates.
To run the script file, type the following command at the shell prompt:
perl run-diffmiss-qc.pl clean-inds-GWA-data
The command creates a file called “fail-diffmiss-qc.txt”, which can be used to exclude these SNPs from downstream association analyses.
Removal of all SNPs failing QC
To remove low-quality SNPs, type the following command at the shell prompt:
plink --bfileclean-inds-GWA-data --exclude fail-diffmiss-qc.txt \
--geno 0.05--hwe 0.00001 --make-bed --out clean-GWA-data
In addition to removing SNPs identified with differential call rates between cases and controls, this command removes SNPs with call rate less than 5% with --geno option and deviation from HWE (p<10-5) with the --hwe option. One additional option is --maf which excludes all SNPs with minor allele frequency less than a specified threshold.
This command will produce cleaned binary ped files for downstream association analyses: “clean-GWA-data.bed”, “clean-GWA-data.bim” and “clean-GWA-data.fam”.