Skip to content

R GWAS V: Results

Jump to bottom
Bioinformatics Core Facility Gothenburg edited this page Sep 19, 2017 · 9 revisions

GWAS Results

Now we will use the two result files and the frequency file (gwas-data.clean.frq) to investigate the results and visualize them.

To pass this part of the course you should for both GWASes present:

a) A table with the top 10 SNPs with the following columns: 1. SNP - Name of SNP
2. Chr - Chromosome
3. Pos - Position
4. A1 - Minor allele
5. A2 - Major allele
6. MAF - Minor allele frequency
7. BETA, OR - Effect or OR (depending on phenotype)
8. SE - Standard Error of the effect, for the quantitative trait
9. P - P-value
10. N - Number of individuals

b) A manhattan plot
c) A qq-plot
d) A short summary of what you have found: Were there any significant results? Are the results the same in both analyses and what do you think is the reason for that? Put this together in one document (word, powerpoint or similar) so it is easy to read.

Merge results with marker information.

For this part you need the files diseasetrait.assoc, quanttrait.qassoc and gwas-data.clean.frq files.

Not all information we want is in the result files from PLINK. Most notably, the minor allele frequency (MAF) and the number of individuals (N) that had a non-missing genotype for the different SNPs are missing. But this information is in the gwas-data.clean.frq file.

  1. Read the gwas-data.clean.frq file to R using the read.table() function and store it in a data frame called info.
  2. Look at the data frame using the head() function.

The A1, A2 and MAF columns we can use as they are. The NCHROBS columns tells how many alleles that were non missing. If we divide that column by 2 we will get the number of individuals analyzed for each SNP. The SNP column we need to keep as an ID column so we can merge the SNP information with the correct result. The CHR column we have no use of so this one we will remove.

  1. Create a new column N in your data frame that is NCHROBS/2. To add one extra column in the end of a data frame you can just set it in the way you usually do with a variable (write your own code instead of the dots)
info$N <- …
  1. Remove the two columns CHR and NCHROBS from your data frame.

Your info data frame should now have five columns: SNP, A1, A2, MAF and N.

Disease Trait GWAS Results

  1. Read the PLINK results from the disease trait association analysis and store it in a data frame called dis.
  2. The results from PLINK contains a bit more information than we need for the final result table. Remove the columns you don’t need: A1 and A2 (since they are included in the info data frame) F_A, F_U and CHISQ.

We are now going to merge the results with the SNP information using the merge() function. This function takes a lot of arguments. If you look at the help file you will see them all. We are going to use

  • x= Your result data frame
  • y= Your SNP information data frame
  • by.x= The ID column name (written within citation marks) in your result data frame, i.e the SNP column.
  • by.y= The ID column name (written within citation marks) in your SNP information data frame, i.e the SNP column.
  • all.x=T Since we want to use all the results even if they are not included in the information file. If there are top hits without information we will solve that later.
  • all.y=F Since we don’t need frequency and N for SNPs without result.
  1. Merge the results with the SNP information in to a new data frame called dismerge.
  2. Make sure you have all the columns you need for the final result table using head()

Quantitative Trait GWAS Results

Now we will do the same for the quantitative trait GWAS

  1. Read in the quantitative results to a data frame called quant.
  2. Remove the columns you don’t need for the final table (keep the SNP column since it’s the ID column)
  3. Merge the quant results with the SNP information to a quantmerge data frame
  4. Check that you have all you need for the final table.

Sort and clean results

Right now the SNPs are sorted in alphabetic order. This was done by the merge() function. Before they were sorted by position. Nothing of this is really relevant for us when we want to know which SNPs are the most significant. Instead we want to sort the results based on the p-value and put the SNP with the smallest p-value in the top.

For this we are going to use the function order(). This function takes a vector and returns a numerical vector with the order of the first vector. Note that it doesn’t sort the input vector.

Disease Trait GWAS Results

  1. Find out the order of the p-values in the dismerge data frame by giving dismerge$P as argument to the order() function. Store the order in a vector called disord.
  2. Reorder the rows in your results using the disord vector.
  3. You can use the same technique to reorder the columns of your result data frame so that they match the order that your final reported table should have. In this case you have to define the order vector yourself.

In your results there are some SNPs that have Beta or OR equal to NA. This happens when there is no variation (or very little variation) in the genotypes in the population. The same SNPs will also have SE and P equal to NA. SNPs without variation is of no interest to us so we will filter them out. One way is to remove all SNPs with MAF=0 but this will leave a few SNPs with very low MAF but still no result. Instead we can check that the reason for that we don’t have a result is really due to small variations.

  1. The following code will produce a summary of the MAF for SNPs where the Odds Ratio is NA. If the maximum value of these MAF is still very low we can exclude all SNPs where the result is NA:\
summary(disresults$MAF[is.na(disresults$OR)])

Hopefully you can conclude that the reason for a NA result is low variance in the genotypes. To find the rows with NA-results we can use a combination of two functions: is.na() and which(). The is.na() function takes a numerical vector and returns a boolean vector with TRUE on the positions where the input vector had NA and FALSE on all other positions. The which() function takes a boolean vector and returns the positions for TRUE as numbers. If we use them together we can get the row numbers where we have NA results:\

ix<-which(is.na(disresults$OR))
  1. Find the row numbers where OR=NA and remove the corresponding rows from the data frame.
  2. When you are pleased with the result data frame, export the 10 top rows (i.e the 10 most significant SNPs and their results) to a file and save it for the report.

Quantitative Trait GWAS Results

Now we will do the same for the quantitative GWAS results.

  1. Sort your quantmerge data frame with the smallest p-value on top using the order() function.
  2. Order the columns so they match the final table.
  3. Check the MAF for SNPs with NA result and remove the SNPs without result from the result data frame.
  4. Export the top 10 SNP results to a file and save for the report.

Next: Plots

Previous: Data analysis

Clone this wiki locally