{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hierarchical Bayesian modelling of gene expression time series across irregularly sampled replicates and clusters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook accompanies the paper\n", "\n", "**Hierarchical Bayesian modelling of gene expression time series\n", "across irregularly sampled replicates and clusters** \n", "_James Hensman, Neil D. Lawrence and Magnus Rattray_ \n", "\n", "\n", "The code in this notebook uses the **development** version of GPy . \n", "\n", "You can download this notebook at \n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#config\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "from matplotlib import pyplot as plt\n", "import numpy as np\n", "import GPy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simple Hierarchical sampling" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we'll illustrate samples from the model by constructing GPy kernels and drawing random samples. We'll create a grid of time-points X on which to visualize the samples, and two kernels (covariance functions), which will form the upper and lower part of the hierarchy." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "T = np.linspace(0, 10, 400).reshape(400,1)\n", "\n", "#construce kernel (covariance function) objects\n", "kern_upper = GPy.kern.Matern32(input_dim=1, variance=1.0, lengthscale=3.0)\n", "kern_lower = GPy.kern.Matern32(input_dim=1, variance=0.1, lengthscale=5.0)\n", "\n", "#compute the covariance matrices\n", "K_upper = kern_upper.K(T)\n", "K_lower = kern_lower.K(T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll draw the functions representing the underlying function describing the $n^\\textrm{th}$ gene $g_n(t)$, and the replications of realisations $f_{nr}(t)$. The key concept is that the draw $g_n(t)$ becomes the mean of the replicate-functions $f_{nr}(t)$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ " /home/maxz/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/figure.py:1653: UserWarning:This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect.\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "g = np.random.multivariate_normal(np.zeros(400), K_upper)\n", "for r in range(5):\n", " f = np.random.multivariate_normal(g, K_lower)\n", " plt.plot(T, f, 'b', linewidth=1, label='$f_{nr}(t)$' if r==0 else None)\n", "plt.plot(T, g, 'r', linewidth=2.5, label='$g_n(t)$')\n", "plt.legend()\n", "plt.xlabel('time')\n", "plt.ylabel('Simulated gene expression')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hierarchical covariance functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the functions $g_n(t)$ and $f_{nr}(t)$ are related linearly, they are _jointly distributed_ as a Gaussian process, which has two inputs. The first input is time, and the second indexes which function we are interested in. We'll use 0...R to represent $f_{n1}(t)\\ldots f_{nR}(t)$.\n", "\n", "GPy can construct these covariances from their underlying covariances using the hierarchical object." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "

hierarchy.	Value	Constraint
upper.variance	1.0	+ve
upper.lengthscale	2.0	+ve
lower.variance	0.1	+ve
lower.lengthscale	4.0	+ve

" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#construct a hierarchical GPy kernel. \n", "kern_upper = GPy.kern.Matern32(input_dim=1, variance=1.0, lengthscale=2.0, active_dims=[0], name='upper')\n", "kern_lower = GPy.kern.Matern32(input_dim=1, variance=0.1, lengthscale=4.0, active_dims=[0], name='lower')\n", "k_hierarchy = GPy.kern.Hierarchical(kernels=[kern_upper, kern_lower])\n", "\n", "#construct a 'grid' on which to examine samples\n", "#the first column contains time (T), thrice. the second columns contains\n", "#zeros, ones and twos to represent three replications. \n", "X = np.r_[np.c_[T, np.zeros_like(T)], np.c_[T, np.ones_like(T)], np.c_[T, np.ones_like(T)*2]]\n", "\n", "#print to show the parameters of the hierarchical kernel are as those of both separate kernels\n", "k_hierarchy\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because of the way we have structured the matrix X to contain blocks of time, the covariance contains this block structure. The diagonal (strong) block contain the covaraince of points on the _same_ replicate, and the off-diagonal blocks contain covariance between points of different replicates. \n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.matshow(k_hierarchy.K(X), cmap=plt.cm.hot); plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's now possible to draw samples from the replicates simultaneously. Notice that here we don't have a draw for the underlying function $g_n(t)$: this has been marginalised out. We can perform inferences for $g$, but since it's an unobserved latent process, it's not required for computation. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sample = np.random.multivariate_normal(np.zeros(X.shape[0]), k_hierarchy.K(X))\n", "for r in range(3):\n", " index = X[:,1]==r\n", " plt.plot(X[index,0], sample[index], 'b', linewidth=1)\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Drosophila data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we model the data in \n", "\n", "**Gene expression divergence recapitulates the developmental hourglass model** \n", "_Alex T. Kalinka,\tKarolina M. Varga,\tDave T. Gerrard,\t Stephan Preibisch,\t David L. Corcoran,\t Julia Jarrells,\tUwe Ohler,\tCasey M. Bergman and Pavel Tomancak\t_ \n", "Nature 468, 811–814 (09 December 2010)\n", "\n", "The data are available at http://publications.mpi-cbg.de/publications-sites/4240/Supplemental_Information_Summary.html#Raw_data. These are a time series experiment across several species, with several replicates. We have downloaded these data, and filtered for only the Melanogaster rows. This results in 57 measurements of 3590 genes, spread (unevenly) across 8 replicates. \n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import urllib\n", "urllib.urlretrieve('http://staffwww.dcs.sheffield.ac.uk/people/J.Hensman/data/kalinka09_mel.csv', 'kalinka_data.csv')\n", "urllib.urlretrieve('http://staffwww.dcs.sheffield.ac.uk/people/J.Hensman/data/kalinka09_mel_pdata.csv', 'kalinka_pdata.csv')\n", "expression = np.loadtxt('kalinka_data.csv', delimiter=',', usecols=range(1, 57))\n", "gene_names = np.loadtxt('kalinka_data.csv', delimiter=',', usecols=[0], dtype=np.str)\n", "replicates, times = np.loadtxt('kalinka_pdata.csv', delimiter=',').T\n", "\n", "#normalize data row-wise\n", "expression -= expression.mean(1)[:,np.newaxis]\n", "expression /= expression.std(1)[:,np.newaxis]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inference in a single gene" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's a simple plot of the data in a single gene. Each sup-plot represents one replicate. Note how some replicates do not contain all the time points. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(10,1.85))\n", "gene_index = 12 # hoppefully a typical gene\n", "for r in range(1,9):\n", " plt.subplot(1,8,r)\n", " index = replicates==r\n", " plt.plot(times[index], expression[gene_index, index], 'kx', mew=1.2, ms=4)\n", " plt.xlabel('replicate %i'%r)\n", " \n", " if r==1:plt.ylabel('gene expression')\n", "GPy.plotting.matplot_dep.base_plots.align_subplots(1,8, xlim=(0,22))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To construct a GP model, we need to for the matrix X as for the sampling above, with one column representing the times of observation, and the other containing the replicate index.\n" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " :0: FutureWarning:IPython widgets are experimental and may change in the future.\n" ] } ], "source": [ "X = np.hstack((times.reshape(-1,1), replicates.reshape(-1,1)))\n", "Y = expression[gene_index].reshape(-1,1) #must be a column vector for GPy\n", "m = GPy.models.GPRegression(X=X, Y=Y, kernel=k_hierarchy)\n", "m.likelihood.variance = 0.01\n", "m.optimize('bfgs', messages=1)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(10,1.85))\n", "#to plot the mean function g_n(t), predict using only the 'upper' kernel\n", "Xplot = np.linspace(0,22, 100)[:,None]\n", "mu, var = m.predict(Xplot, kern=kern_upper)\n", "ax = plt.subplot(1,9,1)\n", "GPy.plotting.matplot_dep.base_plots.gpplot(Xplot, mu, mu - 2*np.sqrt(var), mu + 2*np.sqrt(var), ax=ax, edgecol='r', fillcol='r')\n", "plt.ylabel('gene expression')\n", "plt.xlabel('underlying\\nfunction')\n", "\n", "#plot each of the functions f_{nr}(t)\n", "for r in range(1,9):\n", " ax = plt.subplot(1,9,r+1)\n", " m.plot(fixed_inputs=[(1,r)],ax=ax, which_data_rows=X[:,1]==r)\n", " plt.xlabel('replicate %i'%r)\n", " plt.plot(Xplot, mu, 'r--', linewidth=1)\n", " \n", "GPy.plotting.matplot_dep.base_plots.align_subplots(1,9, xlim=(0,22), ylim=(-2,2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above, we've plotted the hierarchical GP posterior, including our posterior over the underlying function. For each replicate, the shaded area represents a 95 confidence interval in another point _in that replicate_. We also show the mean of the underlying function as red dashes: it's clear that there is a strong batch effect which our model copes with well. Note that this would be hard to normalize-out because of the disorganised replicate structure. \n", "\n", "To demonstrate what we've gained from the hierarchical model, let's compare with two simpler methods. First, fitting a GP to all the data without considering the replicate structure, and second, averaging at each time point before fitting a GP." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(10,3))\n", "#lump all the data together into a single GP\n", "m_lump = GPy.models.GPRegression(X[:,0:1], Y, kernel=GPy.kern.Matern32(1))\n", "m_lump.optimize()\n", "m_lump.plot(ax=plt.subplot(1,3,1))\n", "plt.plot(Xplot, mu, 'r--', linewidth=1)\n", "plt.xlabel('Replicate structure ignored')\n", "\n", "#average each time point\n", "X_unique = np.unique(X[:,0])\n", "Y_u = []\n", "for x in X_unique:\n", " Y_u.append(Y[X[:,0]==x].mean())\n", "Y_u = np.array(Y_u).reshape(-1,1)\n", "m_ave = GPy.models.GPRegression(X_unique[:,None], Y_u, kernel=GPy.kern.Matern32(1))\n", "m_ave.optimize()\n", "m_ave.plot(ax=plt.subplot(1,3,2))\n", "plt.plot(Xplot, mu, 'r--', linewidth=1)\n", "plt.xlabel('Each time point\\n averaged independently')\n", "\n", "#plot the hierarchical result\n", "ax = plt.subplot(1,3,3)\n", "GPy.plotting.matplot_dep.base_plots.gpplot(Xplot, mu, mu - 2*np.sqrt(var), mu + 2*np.sqrt(var), ax=ax, edgecol='r', fillcol='r')\n", "GPy.plotting.matplot_dep.base_plots.align_subplots(1,3, xlim=(0,22), ylim=(-2,2))\n", "plt.xlabel('Underlying function from\\n the hierarchical model')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Deeper Hierarchies" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's straightforward to add another layer of the hierarchy to ther model. Here, we'll take some genes that have been clustered together and build a Hierarchical GP with three layers. At the top, is an underlying function which governs the behaviour of all the genes. Beneath that is another GP which models gene-specific deviations from the cluster, and below that a GP to model replicate variation as above. \n", "\n", "We'll construct a design matrix X which contains three columns. The first is time, the next indexes the gene and the third indexes the replicate. We'll arrange the observed gene expression values into a single column vector Y similarly. " ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [], "source": [ "k_cluster = GPy.kern.Matern32(input_dim=1, active_dims=[0], name='cluster', variance=1.5, lengthscale=2.5)\n", "k_gene = GPy.kern.Matern32(input_dim=1, active_dims=[0], name='gene', variance=0.1, lengthscale=2.5)\n", "k_replicate = GPy.kern.Matern32(input_dim=1, active_dims=[0], name='replicate', variance=0.1, lengthscale=2.5)\n", "\n", "k_hierarchy = GPy.kern.Hierarchical([k_cluster, k_gene, k_replicate])" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "clustered_gene_names = np.array(['ac', 'bib', 'yellow-e3', 'Tom', 'CG13333', 'tld'], dtype=np.str)\n", "index = [i for i,gn in enumerate(gene_names) if gn in clustered_gene_names]\n", "Y_cluster = expression[index]\n", "\n", "X = []\n", "for i in range(len(index)):\n", " X.append(np.hstack((times.reshape(-1,1), np.ones((times.size,1))*i, replicates.reshape(-1,1) + i*8)))\n", "X = np.vstack(X)\n", "Y = np.vstack([y.reshape(-1,1) for y in Y_cluster])" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [], "source": [ "m = GPy.models.GPRegression(X, Y, kernel=k_hierarchy)\n", "m.optimize(messages=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's it for inference. All that remains is to plot the GP predictions. Here we'll follow the convention above of plotting the replicates of a gene in rows, with the first column showing the underlying function for that gene. \n", "\n", "In the bottom row we'll plot a prediction for a new 'unseen' gene. Note that the errorbars are larger since they comprise the gene-variance, replicate-varaince and the noise. The very bottom-left plot shows the posterior of the top-level underlying function (green). " ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "f, ax_array = plt.subplots(len(clustered_gene_names)+1, 9, figsize=(10,10))\n", "#create the kernel for making predictions of the gene-specific functions (turn off replicate variance)\n", "k_pred = k_hierarchy.copy()\n", "k_pred.replicate.variance=0\n", "for g, gn in enumerate(clustered_gene_names):\n", " Xpred = np.hstack((Xplot, Xplot*0+g, Xplot*0))\n", " mu, var = m.predict(Xpred, kern=k_pred)\n", " ax = ax_array[g,0]\n", " GPy.plotting.matplot_dep.base_plots.gpplot(Xplot, mu, mu - 2*np.sqrt(var), mu + 2*np.sqrt(var), ax=ax, edgecol='r', fillcol='r')\n", " ax.set_ylabel(gn)\n", " for r in range(8):\n", " ax = ax_array[g,r+1]\n", " m.plot(fixed_inputs=[(1,g),(2,r+g*8)], ax=ax, which_data_rows=np.logical_and(X[:,1]==g, X[:,2]==((r+1)+g*8)))\n", "\n", " \n", "#create the kernel for making predictions of the overall function (turn off the replicate- and gene-variance)\n", "k_pred = k_hierarchy.copy()\n", "k_pred.replicate.variance = 0.0\n", "k_pred.gene.variance = 0.0\n", "mu, var = m.predict(Xpred, kern=k_pred)\n", "ax = ax_array[-1,0]\n", "GPy.plotting.matplot_dep.base_plots.gpplot(Xplot, mu, mu - 2*np.sqrt(var), mu + 2*np.sqrt(var), ax=ax, edgecol='g', fillcol='g')\n", "\n", "#predict the gene expression for a completely new gene (give it index 100)\n", "Xpred = np.hstack((Xplot, Xplot*0+100, Xplot*0 + 100))\n", "mu, var = m.predict(Xpred)\n", "[GPy.plotting.matplot_dep.base_plots.gpplot(Xplot, mu, mu - 2*np.sqrt(var), mu + 2*np.sqrt(var), ax=ax_array[-1,i], edgecol='y', fillcol='y') for i in range(1,9)]\n", "\n", "#set some labels\n", "ax_array[-1,0].set_ylabel('Unseen gene')\n", "ax_array[-1,0].set_xlabel('Underlying\\nfunctions')\n", "[ax_array[-1,i].set_xlabel('replicate %i'%(i)) for i in range(1,9)]\n", "GPy.plotting.matplot_dep.base_plots.align_subplot_array(ax_array, xlim=(0,22), ylim=(-2,2))\n", "plt.tight_layout()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }