Repository navigation
R GWAS III: Quality control
We will try to look at the ped-file by reading part of it in to R. Since it contains two columns per SNP (for very many SNPs) plus some individual information, it is heavy for the computer to open the file as it is. So we use PLINK to extract a small part of it.
- Run:
plink --file gwas-data --chr 1 --from-bp 3696000 --to-bp 3697000 --recode --out tmp - Open the tmp.ped and the tmp.map file in a text editor of your choice (not Word...) or Excel, and have a look.
The first 6 columns of tmp.ped contains information about the individual and the following columns are alleles for the SNPs.
The actual ped file you downloaded is quite large. To be able to work a bit faster with it we can transform it into a binary file. Then it is no longer possible to open and look at it, but the future calculations will run much faster.
At the same time we will include some quality control (QC) steps so that we have clean, high quality data to work with. In the QC we will filter out
- SNPs with more than 10% missing genotypes
- individuals with more than 10% missing genotypes
- SNPs that are not in Hardy Weinberg Equilibrium using the limit p ≤ 0.001.
To do this we will use the following commands (all at the same time):
-
--filefollowed by the name of your data (in this case gwas-data) tells PLINK what input data to use -
--missingwill calculate missing rates for SNPs and individuals -
--mind 0.1will filter out individuals with more than 10% missing genotypes -
--geno 0.1will filter out SNPs with 10% missing genotypes -
--hwe 0.001will do a Hardy Weinberg Equilibrium test and filter out SNPs with p ≤0.001 -
--make-bedcreates a binary (.bed) file (and a .bim with SNP info, and a .fam with sample info) -
--outfollowed by the name of the output data will print the new filtered binary data into a file.
- Open a command prompt in the folder where your data is stored.
- Use the following PLINK command to do all of the steps mentioned above at the same time:
plink --file gwas-data --missing --mind 0.1 --geno 0.1 --hwe 0.001 --make-bed --out gwas-data.clean
If you already got the binary data, you just tell PLINK that your input is binary, by using --bfile instead of --file:
plink --bfile gwas-data --missing --mind 0.1 --geno 0.1 --hwe 0.001 --make-bed --out gwas-data.clean
Q1: How many individuals passed this QC, and how many were removed?
Q2: How many markers passed QC?
One thing that we are often interested in and that we sometimes filter the results on is allele frequency. It is a good idea to not filter your genotype data on allele frequency. But when we look at the results later on it is good to know the allele frequencies, so we are going to calculate them. That is done using the PLINK argument --freq. Since we created a binary file we will now change the --file argument to --bfile instead.
- Run PLINK with the following command to create a file with allele frequency data:
plink --bfile gwas-data.clean --freq --out gwas-data.clean
You can open the gwas-data.clean.frq file you got in a text editor and look at it. You will see 6 columns for chromosome, SNP, minor allele (A1), major allele (A2), minor allele frequency (MAF) and number of non-missing alleles in your data.