{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# BayesA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulating Genotypes and Phenotypes" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "using(Distributions)\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nObs = 100\n", "nMarkers = 1000\n", "X = sample([0,1,2],(nObs,nMarkers))\n", "α = randn(nMarkers)\n", "a = X*α\n", "stdGen = std(a)\n", "a = a/stdGen\n", "y = a + randn(nObs)\n", "saveAlpha = α\n", "nothing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "####Centering Genotype Covariates" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "meanXCols = mean(X,1)\n", "X = X - ones(nObs,1)*meanXCols;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Priors" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "seed = 10 # set the seed for the random number generator\n", "chainLength = 2000 # number of iterations\n", "dfEffectVar = 4 # hyper parameter (degrees of freedom) for locus effect variance \n", "nuRes = 4 # hyper parameter (degrees of freedom) for residual variance\n", "varGenotypic = 1 # used to derive hyper parameter (scale) for locus effect variance \n", "varResidual = 1 # used to derive hyper parameter (scale) for locus effect variance \n", "scaleVar = varGenotypic*(dfEffectVar-2)/dfEffectVar # scale factor for locus effects\n", "scaleRes = varResidual*(nuRes-2)/nuRes # scale factor for residual variance\n", "nothing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Function for Sampling Marker Effects" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "get_column (generic function with 1 method)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function get_column(X,nRows,j)\n", " indx = 1 + (j-1)*nRows\n", " ptr = pointer(X,indx)\n", " pointer_to_array(ptr,nRows)\n", "end" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "xpx = [(X[:,i]'X[:,i])[1]::Float64 for i=1:nMarkers]\n", "xArray = Array(Array{Float64,1},nMarkers)\n", "for i=1:nMarkers\n", " xArray[i] = get_column(X,nObs,i)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computing the adjusted right-hand-side efficiently\n", "\n", "We want to compute:\n", "$$\n", "rhs = \\mathbf{X}_j'(\\mathbf{y}_{corr} + \\mathbf{X}_j\\mathbf{\\alpha}_j)\n", "$$\n", "This is more efficiently obtained as\n", "$$\n", "rhs = \\mathbf{X}_j'\\mathbf{y}_{corr} + \\mathbf{X}_j'\\mathbf{X}_j\\mathbf{\\alpha}_j,\n", "$$\n", "using the diagonals of $\\mathbf{X}'\\mathbf{X}$ that have already been computed (line 4 of the function below)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "sampleEffects! (generic function with 1 method)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function sampleEffects!(nMarkers,xArray,xpx,yCorr,α,meanAlpha,vare,locusEffectVarl)\n", " nObs = size(X,1)\n", " for j=1:nMarkers\n", " rhs::Float64 = dot(xArray[j],yCorr) + xpx[j]*α[j] \n", " lhs::Float64 = xpx[j] + vare/locusEffectVar[j]\n", " invLhs::Float64 = 1.0/lhs\n", " mean::Float64 = invLhs*rhs\n", " oldAlpha::Float64 = α[j] \n", " α[j] = mean + randn()*sqrt(invLhs*vare)\n", " BLAS.axpy!(oldAlpha-α[j],xArray[j],yCorr) \n", " end\n", " nothing\n", "end" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1000" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "length(xpx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Function for BayesA\n", "\n", "The intercept is sampled first and the sampleEffects! function is called to sample the marker effects." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "BayesC0! (generic function with 1 method)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "chi1=Chisq(nObs+nuRes)\n", "chi2=Chisq(dfEffectVar+1)\n", "\n", "function BayesC0!(numIter,nMarkers,X,xpx,yCorr,mu,meanMu,α,meanAlpha,vare,locusEffectVar)\n", " for i=1:numIter\n", " # sample residula variance\n", " vare = (dot(yCorr,yCorr)+nuRes*scaleRes)/rand(chi1)\n", " \n", " # sample intercept\n", " yCorr = yCorr+mu\n", " rhs = sum(yCorr) \n", " invLhs = 1.0/(nObs) \n", " mean = rhs*invLhs \n", " mu = mean + randn()*sqrt(invLhs*vare)\n", " yCorr = yCorr - mu \n", " meanMu = meanMu + (mu - meanMu)/i\n", " \n", " # sample effects\n", " sampleEffects!(nMarkers,xArray,xpx,yCorr,α,meanAlpha,vare,locusEffectVar)\n", " meanAlpha = meanAlpha + (α - meanAlpha)/i\n", " \n", " #sameple locus effect variance\n", " \n", " for j=1:nMarkers\n", " locusEffectVar[j] = (scaleVar*dfEffectVar + α[j]^2)/rand(chi2)\n", " end\n", "\n", " if (i%1000)==0\n", " yhat = meanMu+X*meanAlpha\n", " resCorr = cor(a,yhat)\n", " println (\"Correlation of between true and predicted breeding value: \", resCorr)\n", " end\n", " end\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##Run BayesA" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Correlation of between true and predicted breeding value: 0.7886358253499116\n", "Correlation of between true and predicted breeding value: 0.7890681933545429\n", "elapsed time: 1.045476748 seconds (243489164 bytes allocated, 13.68% gc time)\n" ] } ], "source": [ "meanMu = 0\n", "meanAlpha = zeros(nMarkers) \n", "\n", "#initial valus\n", "vare = 1\n", "varEffects = 1.0\n", "mu = mean(y)\n", "yCorr = y - mu\n", "alpha = fill(0.0,nMarkers)\n", "locusEffectVar = fill(varEffects,nMarkers)\n", "\n", "#run it\n", "@time BayesC0!(chainLength,nMarkers,X,xpx,yCorr,mu,meanMu,alpha,meanAlpha,vare,locusEffectVar)" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.3.7", "language": "julia", "name": "julia 0.3" }, "language_info": { "name": "julia", "version": "0.3.7" } }, "nbformat": 4, "nbformat_minor": 0 }