LmmGwas/demo_LmmGWAS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GWAS Analysis using Linear Mixed Effect Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we show you how to use the VarianceComponentModels.jl of the Open Mendel project to do a standard GWAS that accounts for relatedness or population substructure using as a random effect. The data used in this tutorial are from the Mendel version 16.0 option 29 (http://software.genetics.ucla.edu/download?package=1). They are simulated data and are freely available but please acknowledge the Open Mendel project if you use them. Strictly, because the data are only snps on chromosome 19, but with ~140,000 snps you will get an idea of the capabilities of Open Mendel. \n",
    "\n",
    "To use this tutorial you will need to have installed SnpArrays, VarianceComponentModels and MendelPlots from the Open Mendel project. To do so, please open julia in a terminal, then use the julia package manager. Add SnpArrays first, then VarianceComponentModels and finally MendelPlots. \n",
    "\n",
    "] #invokes the package manager\n",
    "\n",
    "add https://github.com/OpenMendel/SnpArrays.jl.git\n",
    "\n",
    "add https://github.com/OpenMendel/VarianceComponentModels.jl.git\n",
    "\n",
    "add  https://github.com/OpenMendel/MendelPlots.jl.git\n",
    "\n",
    "You will also need these registered packages, DataFrames, CSV, Distributions, DelimitedFIles, and LinearAlgebra and can add them using the package manager if you haven't already:\n",
    "\n",
    "add Distributions.jl \n",
    "\n",
    "add DelimitedFiles.jl  \n",
    "\n",
    "add LinearAlgebra.jl\n",
    "\n",
    "add CSV.jl\n",
    "\n",
    "add DataFrames.jl\n",
    "\n",
    "\n",
    "The tutorial has been tested with julia 1.1.0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "versioninfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Required Packages  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# packages from openMendel\n",
    "using SnpArrays,VarianceComponentModels\n",
    "# packages from Julia base\n",
    "using Distributions, DelimitedFiles, LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in the family structure and the trait\n",
    "\n",
    "In this example we will use one of the two simulated traits found in the fam file.  We will also use sex as a covariate.  In the fam file sex is denoted as F or M.  We arbitrarily choose M (male) to be the reference group and so change M to 0 when we define the sex variable.  The effect of sex is the change in effect size from male to female, which is changed from a coding of F to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pedLMM = readdlm(\"SNP_29a.fam\", ','; header = false)\n",
    "Trait1 = convert(Vector{Float64}, pedLMM[:, 7])\n",
    "# Trait2 = convert(Vector{Float64}, pedLMM[:, 8])\n",
    "# Y = [Trait1 Trait2]\n",
    "sex = map(x -> strip(x) == \"F\" ? 1.0 : 0.0,  pedLMM[:, 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check that the data were read in correctly by typing the name of the variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Trait1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in genotypes and calculate GRM\n",
    "\n",
    "We use SnpArrays to read in the binary snp file.  We also use SnpArrays to calculate the genetic relationship matrix (GRM). In this example we exclude any snps with a minor allele frequency (maf) less than 0.05. Using SNP with maf >0.05 helps insure that the GRM is accurate because rare variants can bias the GRM. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snpbinLMM = SnpArray(\"SNP_29a.bed\")\n",
    "ex29agrm = grm(snpbinLMM; method = :GRM, minmaf=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to know the order of the snps in the bed file so we read in the bim file.  Of course we need to exclude those with maf >0.05 to match the set of snps we used in the GRM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# columns are: :chrom, :snpid, :?, :pos, :allele1, :allele2\n",
    "snpLMM = readdlm(\"SNP_29a.bim\"; header = false)\n",
    "snpLMM = snpLMM[maf(snpbinLMM) .> 0.05,:]\n",
    "snpid = map(x -> strip(string(x)), snpLMM[:, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the data for the covariates.  I\n",
    "In this case we have only sex as a covariate but we could have used other covariates as desired. The ones(n) sets of a variable that has value 1.0 for all individuals. This allows for the estimate of the grand mean $\\mu$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, snps = size(snpbinLMM[:,maf(snpbinLMM) .> 0.05])\n",
    "X = [ones(n) sex]\n",
    "p = size(X,2)  # no. covariates\n",
    "n, snps, p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare to fit LmmGWAS\n",
    "\n",
    "First we analyze the data under the null model of no snp effects. The next three commands set up the data.  Then we need to decide which algorithm we wish to use to get our estimates. We have chosen the MM algorithm. Alternatively we could have used Fisher scoring (FS). The next set of commands then implements the optimization.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit null model once to store nessary information for alternative model \n",
    "nulldata    = VarianceComponentVariate(Trait1, X, (2ex29agrm, Matrix{Float64}(I, n, n)))\n",
    "nulldatarot = TwoVarCompVariateRotate(nulldata)\n",
    "nullmodel   = VarianceComponentModel(nulldata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "algorithm = :MM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if algorithm == :MM\n",
    "    logl_null,_,_,Σcov, = mle_mm!(nullmodel, nulldatarot; verbose = true)\n",
    "elseif algorithm == :FS\n",
    "    logl_null,_,_,Σcov, = mle_fs!(nullmodel, nulldatarot; verbose = true)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Heritablity of `Trait1`\n",
    "\n",
    "We now calculate the narrow sense heritability and its standard error.   The equation for heritability(h) is h = $\\sigma^2_a / (\\sigma^2_a + \\sigma^2_e)$.  Note that in this version of the VarianceComponentModels.jl, we allow for only two variance components, the additive genetic variance and the dependent environmental variance.  In future implementations we will allow for more variance components.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h, hse = heritability(nullmodel.Σ, Σcov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The heritability of this simulated trait is rather on the high side ($72\\%$) for a human trait, which explains why we can get away with only 212 individuals in this GWAS. \n",
    "\n",
    "## GWAS of Trait1\n",
    "\n",
    "We now prepare our alternative models in order to conduct our GWAS of Trait1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## fit alternative model with SNPs, push null model info to alternative model \n",
    "T = eltype(sex)\n",
    "altdatarot = TwoVarCompVariateRotate(nulldatarot.Yrot,\n",
    "    zeros(T, n, size(X, 1) + 1), nulldatarot.eigval, nulldatarot.eigvec,\n",
    "    nulldatarot.logdetV2)\n",
    "copyto!(altdatarot.Xrot, nulldatarot.Xrot) # last column ramains zero\n",
    "altmodel = VarianceComponentModel(altdatarot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop over all SNPs to calculate LRT pvalues for LmmGWAS\n",
    "\n",
    "The following routine shows you how you can write some simple julia code to execute a GWAS. This tutorial is set up to run all the snps. Note that if you are running this tutorial on your old laptop, then be prepared to wait a while for this step to finish. If you see the counter (for every 1000 snps processed) progressing, then the program is working so just be patient.  Alternatively you might wish to try out the tutorial for a much smaller example.  Then loop through only the first few snps, for example the first 100. If so, comment out with a # \"@time for snp in 1:snps\" and remove the # on \"@time for snp in 1:testrun\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvalue   = ones(snps)\n",
    "genovec  = zeros(T, n)\n",
    "testrun  = 100\n",
    "\n",
    "snpsidx = vcat(1:size(snpbinLMM,2))[maf(snpbinLMM) .> 0.05]\n",
    "#@time for snp in 1:testrun\n",
    "@time for snp in 1:snps \n",
    "    # append (rotated) genotype vector to covariate matrix\n",
    "    Base.copyto!(genovec, @view(snpbinLMM[:,snpsidx[snp]]), model=ADDITIVE_MODEL, center=true, scale=true, impute=true)\n",
    "    tmp_mat = similar(genovec)\n",
    "    LinearAlgebra.mul!(tmp_mat, transpose(altdatarot.eigvec), genovec)\n",
    "    altdatarot.Xrot[:, end] = tmp_mat\n",
    "    # initialize mean effects to null model fit\n",
    "    fill!(altmodel.B, zero(T))\n",
    "    copyto!(altmodel.B, nullmodel.B)\n",
    "    copyto!(altmodel.Σ[1], nullmodel.Σ[1])\n",
    "    copyto!(altmodel.Σ[2], nullmodel.Σ[2])\n",
    "    # fit alternative model\n",
    "    if algorithm == :MM\n",
    "        logl_alt, vcmodel_mle, Σse, Σcov, Bse, Bcov = mle_mm!(altmodel, altdatarot; verbose = false)\n",
    "    elseif algorithm == :FS\n",
    "        logl_alt, = mle_fs!(altmodel, altdatarot; verbose = false)\n",
    "    end\n",
    "    # LRT statistics and its pvalue\n",
    "    lrt = - 2(logl_null - logl_alt)\n",
    "    pvalue[snp] = ccdf(Chisq(1), lrt)\n",
    "#    println(snp,\": \", hapmap_snpdata.snpid[snp], \n",
    "    if mod(snp, 1000) == 1\n",
    "            println(snp)\n",
    "    end    \n",
    "#        \"\\n\\tMAF: \", @sprintf(\"%0.3f\", maf[snp]), \n",
    "#        \"\\n\\tLRT p: \", @sprintf(\"%0.3f\", pvalue[snp]))\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Output results to file\n",
    "In some situations you may want to save GWAS results for future use, for example for use as part of a meta analysis. In the next set of commands, we show you how to make and save a comma delimited file with the snp id, the position of the snp, the chromosome location of the snp, the minor allele frequency, and the result provided as the negative log of the p-value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DataFrames\n",
    "using CSV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# maf, = summarize(snpbinLMM)\n",
    "plot_frame = DataFrame(snpid = snpLMM[:,2],\n",
    "   AdjBasepairs = snpLMM[:,4], \n",
    "   Chromosome = snpLMM[:,1], \n",
    "   MAF = maf(snpbinLMM)[maf(snpbinLMM).>0.05],\n",
    "   Pvalue = pvalue)\n",
    "\n",
    "CSV.write(\"lmmGWAS_output_pVal.txt\", plot_frame)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manhattan Plot\n",
    "\n",
    "One of the most common ways to display the results of a GWAS is as a plot of negative log base 10 of the pvalues versus chromosomal position, a Manhattan plot.  For your convenience, we have developed a Julia plotting module as part of the Open Mendel project, MendelPlots.  We demonstrate its use below.  The Manhattan plot will look at little different than the typical one because this example only includes markers from chromosome 19. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_frame = CSV.read(\"lmmGWAS_output_pVal.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using MendelPlots\n",
    "\n",
    "manhattan(plot_frame; pvalvar = \"Pvalue\", chrvar = \"Chromosome\", \n",
    "    posvar = \"AdjBasepairs\", outfile = \"lmmGWAS_manhattan.png\", fontsize = 18pt, linecolor = \"red\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    " display(\"image/png\", read(\"lmmGWAS_manhattan.png\")) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "This tutorial demonstrates how with just a little extra Julia coding, an Open Mendel user can use the VarianceComponentModels module to conduct a GWAS that takes into account possible relatedness or population substructure among individuals as a random effect.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.1",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}