2_GWAS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running GWAS using limix\n",
    "## Using our phenotype from the first notebook, we are going to run GWAS.\n",
    "\n",
    "### Steps of this notebook:\n",
    "Preparing input data  \n",
    "Running GWAS using the limix library for python  \n",
    "Outputting the results.\n",
    "\n",
    "Limix is freely available [here](https://github.com/limix/limix) and has an extensive [documentation](https://limix.readthedocs.io/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.  Initial setup steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1a. Set up the environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import h5py\n",
    "from limix.qtl import scan\n",
    "from bisect import bisect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1b. Define variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# phenotype file\n",
    "pheno_file = './data/subset_flowering_time_16.csv'\n",
    "\n",
    "# genotype file\n",
    "geno_file = './data/subset_all_chromosomes_binary.hdf5'\n",
    "\n",
    "# kinship matrix file\n",
    "kin_file = './data/kinship_ibs_binary_mac5.h5py'\n",
    "\n",
    "# minor allele frequency threshold - we will talk more about why we do this in the next notebook!\n",
    "# 10% is the standard cutoff we apply to *Arabidopsis thaliana* GWAS, but this may be different in your organism of choice.\n",
    "MAF_thrs = 0.1\n",
    "\n",
    "# output file for results with K correction\n",
    "output_file = './results/subset_flowering_time_16_gwas.csv'\n",
    "\n",
    "# output file for results without K correction\n",
    "output_file_nc = './results/subset_flowering_time_16_gwas_nc.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Load input for GWAS\n",
    "\n",
    "As a reminder, the linear mixed model for GWAS is:\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "Y   = X\\beta + K + \\epsilon\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "\n",
    "We need three variables to run GWAS in limix:\n",
    "1.  Y = Phenotype matrix (Y in model)\n",
    "2.  G = Genotype matrix (X in model)\n",
    "3.  K = K matrix (K in model)\n",
    "\n",
    "limix will use this input to estimate the association p-value and effect size\n",
    "($\n",
    "\\begin{align}\n",
    "\\beta\n",
    "\\end{align}\n",
    "$) for each SNP.\n",
    "\n",
    "\n",
    "\n",
    "The first step to generating these three variables is to load the data.\n",
    "### Please Note!\n",
    "This script is written to handle data from the *Arabidopsis thaliana* 1001 genome project.  If you are running GWAS using another organism, you will more than likely start with data that is formatted differently!  If this is the case, some parts of this code *will not* work for you out of the box.  However, if you understand the __steps__ that this code takes, you will be able to adapt it to suit your specific input data.\n",
    "\n",
    "Therefore, I would like you to __focus on the general approach__ to generating the Y, G, and K matrices rather than trying to understand every line of code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2a. Load phenotypes\n",
    "The phenotype data are stored in a 2-column .csv file.\n",
    "The first column specifies the accession identifier (\"ecotypeid\"), the second column contains the phenotype value.  Our flowering time dataset has 200 accessions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load phenotype data\n",
    "pheno = pd.read_csv(pheno_file, index_col = 0)\n",
    "\n",
    "# encode the index (the accessions column) to UTF8 so it matches encoding of accessions in the genotype data.\n",
    "# this will allow us to compare the two datasets in section 3a. \n",
    "pheno.index = pheno.index.map(lambda x: str(int(x)).encode('UTF8'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove accessions with missing or non numerical data\n",
    "pheno = pheno.dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# does pheno match our expectations?\n",
    "print(pheno.shape)\n",
    "print(pheno.head(n=5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The phenotypes are loaded."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2b. Load genotypes\n",
    "The genotypes we're going to use are a subset of SNPs obtained from whole-genome resequencing of 1,135 *Arabidopsis thaliana* accessions ([1001 genomes](http://1001genomes.org/)).\n",
    "\n",
    "Genotype data is stored as an [hdf5](https://www.h5py.org/) file, which is a composite data type.  This means that one hdf5 file stores multiple related data sets. \n",
    "\n",
    "The genotype hdf5 file we are using here consists of three data sets:\n",
    "1.  'accessions' contains the accession identifiers.\n",
    "2.  'positions' provides the SNP positions.\n",
    "3.  'snps' gives the SNP calls themselves. SNPs are coded as 0 for reference allele and 1 for alternate allele.\n",
    "\n",
    "The 'positions' dataset is also associated with a small file called an attribute which provides information about the chromosome location ('chr_index').  We will use the attribute later when outputting GWAS results.\n",
    "\n",
    "(If you are interested in learning more about how to use hdf5 files, check out https://www.pythonforthelab.com/blog/how-to-use-hdf5-files-in-python/ after class for a more detailed introduction.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load genotype data\n",
    "geno_hdf = h5py.File(geno_file, 'r')\n",
    "# structure of hdf5 file\n",
    "# does geno_hdf match our expectations? Here, \"key\" refers to the three different data sets\n",
    "for key in geno_hdf.keys():\n",
    "    print(key)\n",
    "    print(geno_hdf[key].shape)\n",
    "    print(geno_hdf[key][0:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2c.  Load K matrix\n",
    "The kinship matrix we are using is an IBS (identity by state) matrix calculated with the complete SNP dataset (rather than the SNP subset we are using).  As discussed in lecture, there any many options for calculating a K matrix, but we are using a pre-calculated matrix to save ourselves the long computation time.\n",
    "\n",
    "The input file is also an hdf5 file.  It contains 2 datasets:\n",
    "1. 'accessions' gives the accession identifiers\n",
    "2. 'kinship' is the kinship matrix.  (In this case, it is calculated for all 1135 accessions in the complete 1001 genomes dataset.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "kin_hdf = h5py.File(kin_file, 'r')\n",
    "\n",
    "print(kin_hdf['accessions'].shape)\n",
    "print(kin_hdf['accessions'][0:5])\n",
    "print(kin_hdf['kinship'].shape)\n",
    "print(kin_hdf['kinship'][0:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.  Generating matrices for GWAS\n",
    "Now that the data have been loaded, we need to do some final manipulations to generate the appropriate input matrices.\n",
    "\n",
    "We have __*three main objectives*__ here:\n",
    "\n",
    "1.  Since limix will not work with missing data, we need to make sure that all three matrices include the __same set of accessions__.\n",
    "2.  We also need to make sure that the accessions appear in the __same order__ in all three matrices.\n",
    "3.  We will also remove any SNPs that don't meet our __minor allele frequency__ threshold. (We will return to the reason for this step in the next notebook.)\n",
    "\n",
    "Again, this code will work for the 1001 genomes data out of the box, but will likely need to be modified for other data.  Therefore, focus on the __steps__ rather than on each individual line of code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3a.  Generate phenotype matrix (Y)\n",
    "\n",
    "We now have the phenotypes stored as the variable called 'pheno'.  To make the phenotype matrix (Y), we need to __*remove non-genotyped accessions*__ and __*order the data*__ to match the genotype data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove non-genotyped accessions from phenotype\n",
    "acn_genotyped = [acn for acn in pheno.index if acn in geno_hdf['accessions'][:]]\n",
    "# subset phenotype data\n",
    "pheno = pheno.loc[acn_genotyped]\n",
    "# order genotypes in phenotype according SNP-matrix\n",
    "acn_indices = [np.where(geno_hdf['accessions'][:] == acn)[0][0] for acn in pheno.index]\n",
    "acn_indices.sort()\n",
    "acn_order = geno_hdf['accessions'][acn_indices]\n",
    "pheno = pheno.loc[acn_order]\n",
    "# transform to a numpy phenotype matrix (Y)\n",
    "Y = pheno.to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The phenotype matrix (Y) is now ready."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3b. Generate genotype matrix (G)\n",
    "Now we can finish the genotype matrix (G).  Here we remove:\n",
    "1. accessions that are not phenotyped \n",
    "2. SNPs with a minor allele frequency below our threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# subset SNP matrix for phenotyped genotypes\n",
    "# be patient - this step can take a few minutes\n",
    "G = geno_hdf['snps'][:, acn_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove SNPs with minor allele frequency below set threshold\n",
    "# count allele 1 and 0 for each SNP\n",
    "AC1 = G.sum(axis = 1)\n",
    "AC0 = G.shape[1] - AC1\n",
    "AC = np.vstack((AC0,AC1))\n",
    "# define the minor allele for each position\n",
    "MAC = np.min(AC, axis = 0)\n",
    "# calculate minor allele frequency\n",
    "MAF = MAC/G.shape[1]\n",
    "# select SNPs with MAF above threshold \n",
    "SNP_indices = np.where(MAF >= MAF_thrs)[0]\n",
    "SNPs_MAF = MAF[SNP_indices]\n",
    "G = G[SNP_indices, :]\n",
    "\n",
    "# transpose SNP-matrix into accessions x SNPs matrix\n",
    "G = G.transpose()\n",
    "geno_hdf.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The genotype matrix (G) is now ready."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3c. Generate kinship matrix (K)\n",
    "To finish the kinship array, subset for accessions that are phenotyped and genotyped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# subset kinship matrix for phenotyped and genotyped accessions\n",
    "acn_indices = [np.where(kin_hdf['accessions'][:] == acn)[0][0] for acn in pheno.index]\n",
    "acn_indices.sort()\n",
    "K = kin_hdf['kinship'][acn_indices, :][:, acn_indices]\n",
    "kin_hdf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The kinship matrix (K) is now ready."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Check input variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have all three of our input variables.  If you are running GWAS using something other than the 1001 genomes data, __*these are the matrices that you need to generate to run limix*__.  Let's quickly look at these variables to make sure we understand their format.  The number of accessions should be the same in all three matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4a. Y matrix (phenotypes)\n",
    "*Y is a numpy array.\n",
    "The number of rows = number of accessions.\n",
    "The number of columns = 1.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(type(Y))\n",
    "print(Y.shape)\n",
    "print(Y[0:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4b. G matrix (genotypes)\n",
    "*G is a numpy array.  The number of rows = number of accessions.  The number of columns = the number of snps.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(type(G))\n",
    "print(G.shape)\n",
    "print(G[0:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4c. K matrix (K matrix)\n",
    "*K is a numpy array.  The number of rows and columns = the number of accessions*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(type(K))\n",
    "print(K.shape)\n",
    "print(K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5.  Run GWAS\n",
    "We are going to run two GWAS with our phenotype of interest.  The first one will include a population structure correction (i.e. the model will include the K matrix) and the second one will not (we will compare the results in the next notebook)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5a.  Run limix with K correction\n",
    "\n",
    "Now that we have phenotype (Y), genotype (G), and K (K) matrices, we can run GWAS using the 'scan' function from limix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "r = scan(G, Y, K = K, lik = 'normal', M = None, verbose = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the output generated by limix.  Limix has run the linear mixed model for each SNP, and this output summarizes the range of values calculated across all these tests.\n",
    "Although this output is quite detailed, we are primarily interested in two values from limix:\n",
    "\n",
    "1. __P-value__ - This is the p-value for the test that a SNP is associated with the phenotype.  Limix calculates this as likelihood-ratio test between our full GWAS model (H2 in the output) and a model that doesn't contain the genotype term (H0 in the output).  The range of these values across all tested SNPs is given under the heading \"Likelihood-ratio test p-values\".\n",
    "2. __Effect size__ - This is the effect of a SNP in the linear model (the beta in the GWAS equation).   The range of these values across all tested SNPs is given as \"cand. effsizes\" in the table below \"Hypothesis 2\".\n",
    "\n",
    "In the next section, we will output the p-values and effect sizes for all SNPs tested."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5b.  Output results for limix with K correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will look at the results of these GWAS is the next notebook.  This cell saves the GWAS results as a .csv file that gives chromosome (chr), position (pos), p-value (pvalue), minor allele frequency (maf), and effect size (GVE) for each SNP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save results\n",
    "\n",
    "# get chromosomes, positions, and minor allele frequencies\n",
    "# you would need to recode this section if working with something other than the 1001 genomes data\n",
    "geno_hdf = h5py.File(geno_file, 'r')\n",
    "chr_index = geno_hdf['positions'].attrs['chr_regions']\n",
    "chromosomes = [bisect(chr_index[:, 1], snp_index) + 1 for snp_index in SNP_indices]\n",
    "positions_all = geno_hdf['positions'][:]\n",
    "positions = [positions_all[snp] for snp in SNP_indices]\n",
    "mafs = SNPs_MAF #from section 3b\n",
    "\n",
    "# get p-value and effect size\n",
    "# these variables are output from limix and should work regardless of initial input data format\n",
    "#extract p-values\n",
    "pvalues = r.stats.pv20.tolist()\n",
    "#extract effect sizes\n",
    "effsizes = r.effsizes['h2']['effsize'][r.effsizes['h2']['effect_type'] == 'candidate'].tolist()\n",
    "\n",
    "gwas_results = pd.DataFrame(list(zip(chromosomes, positions, pvalues, mafs, effsizes)), columns = ['chr', 'pos', 'pvalue', 'maf', 'GVE'])\n",
    "gwas_results.to_csv(output_file, index = False)\n",
    "geno_hdf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5c.  Run limix without K correction\n",
    "\n",
    "We will also run limix without population structure correction (i.e. without K in the model) so we can compare results in the next notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "r_nc = scan(G, Y, lik=\"normal\", M=None, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This output is in the same format as that of section 5a."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5d.  Output results for limix without K correction\n",
    "We will look at the results of these GWAS is the next notebook.  This cell saves the GWAS results as a .csv file that gives chromosome (chr), position (pos), p-value (pvalue), minor allele frequency (maf), and effect size (GVE) for each SNP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save results\n",
    "# link chromosome and position to p-values and effect sizes\n",
    "geno_hdf = h5py.File(geno_file, 'r')\n",
    "chr_index = geno_hdf['positions'].attrs['chr_regions']\n",
    "chromosomes = [bisect(chr_index[:, 1], snp_index) + 1 for snp_index in SNP_indices]\n",
    "positions_all = geno_hdf['positions'][:]\n",
    "positions = [positions_all[snp] for snp in SNP_indices]\n",
    "mafs = SNPs_MAF #from section 3b\n",
    "\n",
    "pvalues = r_nc.stats.pv20.tolist()\n",
    "effsizes = r_nc.effsizes['h2']['effsize'][r_nc.effsizes['h2']['effect_type'] == 'candidate'].tolist()\n",
    "\n",
    "gwas_results = pd.DataFrame(list(zip(chromosomes, positions, pvalues, mafs, effsizes)), columns = ['chr', 'pos', 'pvalue', 'maf', 'GVE'])\n",
    "gwas_results.to_csv(output_file_nc, index = False)\n",
    "geno_hdf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### In this notebook, we have learned about input data, run GWAS, and output GWAS results.\n",
    " ### Let's move on to 3_GWAS_interpretation.ipynb, where we will explore these results."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}