3_GWAS_interpretation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GWAS interpretation\n",
    "## In this notebook, we will explore and interpret our GWAS results.\n",
    "### Topics we will cover:\n",
    "Reading and interpreting Manhattan and QQ plots  \n",
    "Multiple testing correction  \n",
    "The effects of population structure correction   \n",
    "Why we take minor allele frequency (MAF) cutoffs in our SNPs data  \n",
    "Use caution when interpreting the relationship between p-value and effect size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.  Initial setup steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1a. Prepare environment\n",
    "Loading packages and functions into R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(lattice) #lattice graphics package (for QQ plots)\n",
    "options(repr.plot.width=12, repr.plot.height=5)## this sets a larger size for figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1b. Define input variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results from population structure corrected GWAS\n",
    "results_file = \"./results/subset_flowering_time_16_gwas.csv\"\n",
    "\n",
    "# results from non-population structure corrected GWAS\n",
    "results_file_nc = \"./results/subset_flowering_time_16_gwas_nc.csv\"\n",
    "\n",
    "# results from GWAS using a minor allele frequency cutoff of 1% rather than 10%\n",
    "# this GWAS has already been run for you\n",
    "### DO NOT change this file during independent work ###\n",
    "results_file_maf1 = './prerun_results/subset_flowering_time_16_gwas_maf1.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.  Load GWAS results\n",
    "Load the results files we generated in the GWAS notebook (2_GWAS.ipynb).\n",
    "Each row in these files represents one SNP and the columns are: chromosome (chr), position (pos), p-value (pvalue), minor allele frequency (maf), and effect size (GVE).  The first two files were generated in the previous notebook, but the GWAS in the third file has already been done for you (results_file_maf1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results <- read.csv(file=results_file, stringsAsFactors=FALSE)\n",
    "results_nc <- read.csv(file=results_file_nc, stringsAsFactors=FALSE)\n",
    "results_maf1<- read.csv(file=results_file_maf1, stringsAsFactors=FALSE)\n",
    "\n",
    "# Make sure what we've loaded is as expected\n",
    "dim(results)\n",
    "head(results)\n",
    "\n",
    "dim(results_nc)\n",
    "head(results_nc)\n",
    "\n",
    "dim(results_maf1) ### to save space, I only outputted the chromosome, position, and p-value for this GWAS, so it only has 3 columns\n",
    "head(results_maf1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Manhattan plot\n",
    "We will start by viewing the p-values for every individual SNP across the genome in a plot know as a __*Manhattan plot*__. Each point in this plot represents one SNP.  The x-axis is the genomic position of the SNP.  The y-axis is the -log10 of the p-value for the SNP-phenotype association test (the GWAS). Manhattan plots use the -log10 so the smallest p-values become more visible.  The point colors alternate to denote the different chromosomes. Can you see \"peaks\" of SNPs that rise above the background?  What does this plot suggest about the architecture of the trait (simple versus complex)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################################\n",
    "### plotting -log10 p-value across the genome by chromosome\n",
    "###\n",
    "### results = results from notebook 2 (a dataframe with size nsnps x 5)\n",
    "### (function uses the column headers, so don't change those!)\n",
    "### chr.lengths = vector of the lengths of each chromosome in bp\n",
    "### title = stsring with plot title\n",
    "###########################################\n",
    "\n",
    "at.chr.lengths <- c(30427671,19698289,23459830,18585056,26975502) # the lengths of the 5 chromosomes in A. thaliana\n",
    "\n",
    "manhattan.plot <- function(results, chr.lengths, title){\n",
    "  # the plot will change colors to denote new chromsomes.  This section lays the groundwork for this\n",
    "  chr.add <- c(0,cumsum(chr.lengths))[1:5]\n",
    "  max.bp <- sum(chr.lengths)\n",
    "  chr.colors <- rep(c(\"blue\", \"dodgerblue\"),ceiling(length(chr.lengths/2)))\n",
    "  chr.mids <- chr.add + (chr.lengths/2)\n",
    "  results.s <- split(results, results$chr)\n",
    "  \n",
    "  # make the actual plot \n",
    "  # generate an empty plot with the right dimensions\n",
    "  plot(results$pos,-log10(results$pvalue), xlim=c(0,max.bp), type=\"n\", xlab=\"Chromosome\", ylab=\"-log10 p-value\", xaxt=\"n\", main=title)\n",
    "  # set up the x axis\n",
    "  axis(1,at=chr.mids,labels=c(1:5))\n",
    "  # does plotting by chromosome\n",
    "  for(up.chr in 1:length(results.s)){\n",
    "    #print(up.chr)\n",
    "    up.c <- results.s[[up.chr]]\n",
    "    up.add <- chr.add[up.chr]\n",
    "    up.c$pos.plot <- up.c$pos + up.add\n",
    "    points(up.c$pos.plot, -log10(up.c$pvalue), col=chr.colors[up.chr])\n",
    "  }\n",
    " # adds Bonferroni cutoff\n",
    "  bon.cutoff = 0.05/nrow(results)\n",
    "  abline(h=-log10(bon.cutoff), col=\"red\", lty=2)\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "manhattan.plot(results=results, chr.lengths=at.chr.lengths, title=\"With K correction\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Multiple testing\n",
    "The red dotted line in the Manhattan plot marks the Bonferroni significance cutoff. Since we perform so many tests in a GWAS (one per SNP), we need to apply a method to correct the p-value threshold so we don't identify too many false positive SNPs. The easiest method is the Bonferroni correction, which is done by dividing a typical significance threshold (usually 0.05) by the number of tests done (i.e. the number of SNPs tested). SNPs above the dotted line in this Manhattan plot would be significant above the Bonferroni cutoff.  There are none in the case of our subset flowering time GWAS!\n",
    "\n",
    "However, the Bonferroni method is an extremely conservative strategy for multiple testing, since it implies that all tests are independent (which we know is untrue - there is LD between SNPs). There are many other methods! Two less conservative types of multiple testing correction that are also applied to GWAS are (1) FDR (false discovery rate) and (2) permutation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. The effects of population structure correction\n",
    "To begin to understand how population structure affects GWAS, generate the following two plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qq.correct <- qqmath(~ -log10(results$pvalue), distribution=function(x){-log10(qunif(1-x))},xlab=\"expected -log10 pvalues\", ylab=\"observed -log10 pvalues\", main=\"QQ plot of p-values with K correction\",panel = function(x, ...) {\n",
    "  panel.abline(a=0,b=1)\n",
    "  panel.qqmath(x, ...)})\n",
    "\n",
    "qq.nc <- qqmath(~ -log10(results_nc$pvalue), distribution=function(x){-log10(qunif(1-x))},xlab=\"expected -log10 pvalues\", ylab=\"observed -log10 pvalues\", main=\"QQ plot of p-values without K correction\",panel = function(x, ...) {\n",
    "  panel.abline(a=0,b=1)\n",
    "  panel.qqmath(x, ...)})\n",
    "\n",
    "\n",
    "print(qq.correct, position = c(0, 0, 0.5, 1), more = TRUE)\n",
    "print(qq.nc, position = c(0.5, 0, 1, 1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are *quantile-quantile (QQ) plots*, which let us assess whether our GWAS p-values are inflated.  We have generated one for the GWAS with population structure correction and one for the GWAS without.\n",
    "\n",
    "If the observed p-values follow a uniform (non-inflated) distribution, the blue dots in these plots will fall along the black diagonal line.  If the GWAS has produced significant associations, we __do__ expect the blue dots to be above the line at the very right of the plot - this is due to the significant p-values that are a result of the GWAS. This is exactly what we observe in the QQ plot of the population-structure corrected GWAS (the first plot).\n",
    "\n",
    "Now compare this plot to the QQ plot for the uncorrected GWAS (the second plot). Here, the blue dots rise high above the solid line and this elevation begins at much lower p-values (i.e. more to the left).  This is p-value inflation due to population structure!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's plot the Manhattan plot for the uncorrected data. How does this differ from the Manhattan plot for corrected data? Are inflated p-values obvious?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "manhattan.plot(results=results_nc, chr.lengths=at.chr.lengths, title=\"Without K correction\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Effect of minor allele frequency cutoff\n",
    "\n",
    "## (Don't run this section in independent work - requires pre-run results)\n",
    "\n",
    "Why do we take a minor allele frequency (MAF) cutoff?  First, low frequency SNPs can have complicated LD patterns with many other SNPs. As a results, they are more likely to be false positives due to confounding, which can be seen in p-value inflation.  Second, we have less power when testing low frequency SNPs given the low sample size of the minor allele.  Testing these SNPs therefore adds to our multiple testing correction problem without giving us much useful information in return.\n",
    "\n",
    "Let's look at a Manhattan plot to see what happens to p-values if we lower the cutoff from 10% to 1%.  The Bonferroni cutoff is much higher because we test many more SNPs.  There are also many more high p-value SNPs that don't seem to be associated with a peak - these are likely false positives due to confounding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "manhattan.plot(results=results_maf1, chr.lengths=at.chr.lengths,title=\"Minor allele frequency cutoff of 1%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The QQ plot also reflects this inflation (compare with the 10% MAF cutoff in our original GWAS on the left)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "qq.maf1 <- qqmath(~ -log10(results_maf1$pvalue), distribution=function(x){-log10(qunif(1-x))},xlab=\"expected -log10 pvalues\", ylab=\"observed -log10 pvalues\", main=\"QQ plot of p-values with MAF cutoff of 1%\",panel = function(x, ...) {\n",
    "  panel.abline(a=0,b=1)\n",
    "  panel.qqmath(x, ...)})\n",
    "\n",
    "qq.maf10 <- qqmath(~ -log10(results$pvalue), distribution=function(x){-log10(qunif(1-x))},xlab=\"expected -log10 pvalues\", ylab=\"observed -log10 pvalues\", main=\"QQ plot of p-values with MAF cutoff of 10%\",panel = function(x, ...) {\n",
    "  panel.abline(a=0,b=1)\n",
    "  panel.qqmath(x, ...)})\n",
    "\n",
    "print(qq.maf1, position = c(0, 0, 0.5, 1), more = TRUE)\n",
    "print(qq.maf10, position = c(0.5, 0, 1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. P-value and effect size are not the same thing.\n",
    "\n",
    "One mistake people often make when interpreting GWAS results is to assume that the SNP with the highest p-value explains more phenotypic variation than the others.  There is a relationship between the effect size of a SNP and its p-value, but this relationship is __not__ absolute.  Additionally, a SNP can be significantly-associated with a phenotype, yet not contribute very much to the *overall* phenotypic variation among the accessions tested.  It is therefore worth considering both the p-value *and* the effect size of a SNP you are interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot(-log10(results$pvalue), abs(results$GVE), xlab=\"-log10 pvalue\", ylab=\"effect size\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "In these three notebooks, we have explored a phenotype, used that phenotype for GWAS, and interpreted the resulting output.  For the rest of the class, you now have the chance to run these notebooks yourself with different phenotypes."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}