{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Inferring Transcription Factor Activities by Combining Binding Information with Gene Expression Profiles\n", "\n", "## A Model and Analysis in Progress\n", "\n", "### 30th May 2014\n", "\n", "### M. Arif Rahman and Neil D. Lawrence\n", "\n", "#### modified 30th March 2015 with Sura Al-Zaki\n", "\n", "In this notebook we design a covariance function for reconstructing transcription factor activities given gene expression profiles and a connectivity matrix (binding data) between genes and transcription factors. Our modelling framework builds on ideas in [Sanguinetti et al (2006)](http://bioinformatics.oxfordjournals.org/content/22/14/1753.short) who used a linear-Gaussian statespace modelling framework to infer the transcription factor activity of a group of genes. \n", "\n", "We note that the linear Gaussian model is equivalent to a Gaussian process with a particular covariance function. We therefore build a model directly from the Gaussian process perspective to achieve the same effect. We introduce a computational trick, based on judicious application of singluar value decomposition, to enable us to efficiently fit the Gaussian process in a reduced 'TF activity' space. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import scipy as sp\n", "import pylab as pb\n", "import GPy\n", "import scipy.io\n", "import matplotlib.pyplot as plt\n", "from IPython.display import display\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we load in the classic [Spellman et al (1998)](http://www.molbiolcell.org/content/9/12/3273.full) Yeast Cell Cycle data set. The cdc15 time series data has 23 time points. We can load this gene expression data in with GPy." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Acquiring resource: spellman_yeast\n", "\n", "Details of data: \n", "Two colour spotted cDNA array data set of a series of experiments to identify which genes in Yeast are cell cycle regulated.\n", "\n", "Please cite:\n", "Paul T. Spellman, Gavin Sherlock, Michael Q. Zhang, Vishwanath R. Iyer, Kirk Anders, Michael B. Eisen, Patrick O. Brown, David Botstein, and Bruce Futcher 'Comprehensive Identification of Cell Cycle-regulated Genes of the Yeast Saccharomyces cerevisiae by Microarray Hybridization.' Molecular Biology of the Cell 9, 3273-3297\n", "\n", "After downloading the data will take up 2510955 bytes of space.\n", "\n", "Data will be stored in /home/maxz/tmp/GPy-datasets/spellman_yeast.\n", "\n", "Do you wish to proceed with the download? [yes/no]\n", "yes\n", "combined.txt\n", "Downloading http://genome-www.stanford.edu/cellcycle/data/rawdata/combined.txt -> /home/maxz/tmp/GPy-datasets/spellman_yeast/combined.txt\n", "[==============================] 2.395/2.395MB\n", "Time series of synchronized yeast cells from the CDC-15 experiment of Spellman et al (1998). Two colour spotted cDNA array data set of a series of experiments to identify which genes in Yeast are cell cycle regulated.\n" ] } ], "source": [ "data = GPy.util.datasets.spellman_yeast_cdc15()\n", "Y = data['Y'].fillna(0) # Replace missing values with zero following Sanguinetti et al.\n", "t = data['t']\n", "print data['info'], data['details']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More details on the data are given in [this notebook](../datasets/spellmanYeastData.ipynb).\n", "\n", "We can make a simple helper function to plot genes from the data set (which are provided as a `pandas` array)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEWCAYAAACZnQc8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFDlJREFUeJzt3T93G9eZx/Hfs6s6gqQ3QDJJH9FKpQZnSbtxF9HSGzBp\nJz1la1WILnxMLdOvJecFxFLoKmooUUHDZlf/+jX/vIA1BbvfPFvcO+QQHBAEZsDBXHw/5+AQmL93\nCOKH4TMzd8zdBQBIw7/U3QAAQHUIdQBICKEOAAkh1AEgIYQ6ACSEUAeAhBDqqJSZrZjZrpn908wO\nzexbM7tcU1semtmrguHPzew/4/N/Fjz+z8x+Fce/Lxj/pGib4nJne4YtmtnrON9Wwfh5MzssWFb+\n9/ijmS2X/X1gOhDqqIyZPZS0LmlVUkvSgqQ5Sdt1tMfdv5A0lw9EM1uU9IG7/zE36bxCe7PHFXf/\nJVuMpMXcuA8Utulefplm9khhez03vCVpS9LXcd43kp72NPM7SSe+IMzsrqS7kpbjfJ9Jemhmt4b/\nLWDaEOqoRAywVUnz7v6Du//i7m/d/SNJu2Y2E6dbzO2BbmV7vGY2F/doV+Me/o9mdj23/ML54rjn\nZrbep2nLCoH4q/j6kaRPe6bpxvYePXrGH+bGvZX0vUKwZ64rF+Y5i5JeZ78PhS+8+Vy770ralWQ9\nv8d1SYvu/jKuc1vSF5Ju99lG4AihjqrckPTG3Q96R7j7HXc/iIH1RCFUW5L2FPZUM9fD5H5V0gtJ\nD6WjoDtrvnVJfy1qlLtvSnol6d9jiB66+w89k9npOYvHx7Z8qLAHnq1jw90/l9TtWfffFII9c0Mh\nxGVmc5JWFMJaPdOc+j26+3fufmdAOwFdqrsBSMa8QthKOgqtH3PjP1MIxxfu/o847POeenLX3f8c\nnz/ScXDfPmu+uCd7ls8UwvS9Qomk167ZiVzfdfffZpsi6XXP+Bfu/pcB68za9rMkmdmSpMeSluKo\nJwr/RbzvmWVO0qkaO3BehDqqsqtcecDd92KwS9J/KOxhX5O01BPk+bJFfng+RX89YL4zufu+mT1W\nqJW/K5hkUbkvpJ52uEIQv4mvr0jaNrOFc3yZZHv2TxW2/9/c/Z2ZrUh65e7/MLMrPbPsSrpasJzL\nku64++NB68R0o/yCqmxLms/Xwd39IJYRPoiDdiX9zd2vZg+dLE/08+OI8+V1dXqvOLOXtTU+emvq\n+fFvFfay5wuWU2RbIcB/n/tCWZS0Er+ksnLMYTzu8Frh9zjbs5w7CgdPgTMR6qiEu3cV6sPbZnbL\nzFrxdL3nCnu7rhCGi2a2EMc/UqybD/D0rPnieq73nz1Mdsa4K3G5R48By9pV+O/h7BWGkotLehwP\nBM/F/16ys1pmdFwOmolfGtnv8Xlue5cUjhv0OxgMHCldfsmdLvZrd/+y7PLQXO6+YWZdhdP9nirs\ndX6jEKiz7v6zmX2iUC+fk/RcxzVm6XRJxeNyuwPmW4/rentW8wqWn3ndO8DMFt39ZZ/p9yR9Kenz\novbm3FDYo9/NT+Pu/5pbz34c9ktugo1Yw8+2d1fS3fPW8THdrEx/6ma2oPCv6b6ZPZH06Dx1RgDA\neJQtv8zpuLa5p5Pn7gIALlip8ou7588Vnlefc4UBABejkgOlZjavcOVc0eliAIALUtV56gvufq9o\nhJlxE1QAGIG7D7ra+ZTSe+pmtuLuG/F50dV6cvdkHw8ePKi9DWxffOzvh1Nc9vfT27bU3zu279Rj\nVKVCPfZ4tx47XzrUEFf5AZXqdqWNDWl/P/zsdgfPAySoVKi7+wsPV/j9Jv7sd14vMD7drnT/vvT1\n19LMTPh5/z7BjqnEFaUltdvtupswVo3Yvp2dEOSteCFoqxVe7+ycOVsjtq0Etm86lbr46FwrMPNx\nrwMAUmNm8joOlAIAJgehjuny7NnpWnu3G4YDCSDUMV1u3jx5EDU7yHrzZr3tAipCTR3TJwvy1dVw\n+mP+ICswIUatqRPqmE4HB9LsbDivfWam7tYAp3CgFDgvLlRCwgh1TBcuVELiKL9gujx7Fg6K5mvo\n3W64UOnjj+trF9CDmjoAJISaOoD0cZ3BQIQ6gObgOoOBKL8AaJYpuc6AmjqA6TEF1xlQUwcwHcZx\nncE4a/UXfByAUAfQHOO6zmCctfoLPg5A+QVAc4zzOoNx1upHWDY1dQAoa5y1+iGXTU0dAMoYZ59A\nF9jfEKEOAOPsE+iC+xui/AIA46zVj7jsWmvqZvbQ3b/oM45Qx/DoeAtTrraaupmtSLpVdjmVo4+I\nZuNycGAkpUPd3R9L2qugLdUiFJqt1TquPR4cHNckE7wcHKhSVeWXLXf/qM+4+sovU9JHRNKm4HJw\noAinNBZptUKgz86GnwR6s3DbOWBoaYc6odBcTbztHMdxMAEupPzy4MGDo9ftdlvtdrv0OgfKh0Kr\ndfo1JlsTz37hbw4ldDoddTqdo9dfffVVPac0mtmSpMeS7rr7XwrG11NTb2IooPk4joOK0PcLMCk4\nuIsKcKAUmAQcx0HNCHWgKk08uIvkUH4BqsJxHFSImjoAJISaOgCAUAdQMy7aqhShDqBedL5XKWrq\nAOrHRVuncKAUQLNx0dYJHCgF0FxctFUZQh1Avbhoq1KUXwDUi4u2ClFTB4CEUFMHABDqAJASQh0A\nEkKoA0BCCHUASAihDgAJIdSbip7tmo33D2NCqDcVPds1G+8fxoSLj5qMnu2ajfcPZ+CK0mlFz3bN\nxvuHPrii9CJNSj2Unu2ajfcPY1A61M3slpktmNlyFQ1qhEmoh9KzXbPx/mFMSpVfzGxe0qy7b8ZQ\nf+Xub3umSbP8Unc9lJ7tmo33DwPUUlM3s3VJW+7+0swWJM27+0bPNGmGukQ9FMDY1FVTb0k6zL2+\nVnJ5zUE9FMAEquJA6dDfJI1HPRTAhLpUcv6upKvx+RVJPxVNtLa2dvS83W6r3W6XXG3NdnZO1tBb\nrfCaeiiAEXU6HXU6ndLLKVtTvy7phrt/Z2arkp67+7ueadKtqQOThgOwyailpp6d6RIPknZ7Ax3A\nBZuE021RK64oBVJT9+m2qATdBAA4xum2jUc3AQCCcZxuOyldY2AgQh1IybhOtx22Vs+XQG0IdSAl\nZ51uW0a2nPv3Q2kn++LoV6vngG1tqKkDOL9havUcsC2FmjqA8Rq2Vt9qhUCfnQ0/CfQLQagDGGyU\nWj39I9WC8guAwYa9UjX/JdBqnX6NgThPHcDkoLuC0gh1AEgIB0oBAIQ6AKSEUAeAhBDqAJAQQh0A\nEkKoA0BCCHUASAihDgAJIdQBICGEOgAkhFAHgIQQ6gCQEEIdABJSSaib2cMqlgMAKKd0qJvZiqRb\nFbQF48Kd3YGpUTrU3f2xpL0K2oJx4c7uwNSgpj4NWq3je0oeHHBbsabhPy0MgVCfFtzZvbn4TwtD\nuDRoAjNbLhh86O6b513J2tra0fN2u612u33eWVGV3ju7s6feHPn/tFZXef8S1el01Ol0Si+nknuU\nmtmWu3/UZxz3KK3bpNzZnZsRl3NwEP7T2t+XZmbKL4/3Y6LVdo9SM1uSdMPMPi27LIzJzs7JAM/2\n/HZ2LrYdlBFG1/ufVm+NfRS8H2ly97E+wipwLn//u/v79yeHvX8fhqfi/Xv3P/3JfX8//OzdXpyW\n/c6y31Xv6yqWzfsxcWJ2Dp25lZRfzkL5ZQiTUiYZt6rLCKkbd5mE92Mi1VZ+QYWm4dTDcZQRUvfx\nx6f/BlqtagKd9yM57KlPolT3nKblP5Gm4P2YaOypp2IS9pzGdbHLpBywRcD7kST21CfJpOw5TUo7\ngCk26p46oT5JJum84SzIudgFqAWhjuqlWtsHGoCaOqo1CbV9AEMj1HFavoY+M3N8miXBDkw8yi84\nbZJq+8CUoqYOAAmhpg4AINQBICWEOgAkhFAHgIQQ6gCQEEIdABJCqANAQgh1AEgIoQ4ACak/1Md1\nQwYAmEL1h/rNmyc7i8o6k7p5s952AUADTUbfL9yQAQBOaH6HXtyQAQCO1Nahl5ktx8f6yAvhhgwA\nUIlSoW5mC5JeuPt3kubi6+FwQwYAqEzZPfU5SYvx+V58PZydnZM19FYrvN7ZKdk0AJg+ldXUzWxL\n0l13f9cznJtkAMCQRq2pX6po5fOSXvcGemZtbe3oebvdVrvdrmK1AJCMTqejTqdTejkD99TNbLlg\n8KG7b+amWXX3jT7zs6cOAEOq7ZRGM1tx98fx+YK7b/eMJ9QBYEi1nNJoZouS1s3sRzM7lER6A0CN\nJufiIwDAkdouPgIATA5CHQASQqgDQEIIdQBICKEOAAkh1AEgIYQ6ACSEUAeAhBDqAJAQQh0AEkKo\nA0BCCHUASAihDgAJIdQBICGEOgAkhFAHgIQQ6gCQEEIdABJCqANAQgh1AEgIoQ4ACSHUASAhpUPd\nzJbMbMHMvq2iQQCA0ZUKdTNbkLTg7tuS5szsd9U0CwAwCnP3ahZk9srdbxQM96rWAQDTwszk7jbs\nfJcqWPFlSSuSvim7LABAOVXuqW9J+szd93uGs6cOAEMa2566mS0XDD50900zm5fk7v5W0htJS5I2\neideW1s7et5ut9Vut4dtJwAkrdPpqNPplF5OqT11M1uV9Mbdt+PZL1vu/kPPNOypA8CQRt1TLxvq\nlyXdji/n3P1ewTSEOgAMqZZQP9cKCHUAGNqooc4VpQCQEEIdABJCqANAQgh1AEgIoQ4ACSHUASAh\nhDoAJIRQB4CEEOoAkBBCHQASQqgDQEIIdQBICKEOAAkh1AEgIYQ6ACSEUAeAhBDqAJAQQh0AEkKo\nA0BCCHUASAihDgAJIdQBICGVhbqZrVa1LADAaCoJdTNblPRhFcsCAIyuqj11r2g5AIASSoe6mV13\n9+0qGgMAKKeKPfWrFSwDAFCBS4MmMLPlgsGH7r553r30tbW1o+ftdlvtdnuYNmKSPXsm3bwptVrH\nw7pdaWdH+vjj+toFNEyn01Gn0ym9HHMfvRxuZrfi02uSViQtu/vbnmm8zDow4bpd6f596euvQ7D3\nvgYwEjOTu9uw85Uqv7j7prtvKhwovSwOmE6fVisE+P370sEBgQ7UrNSe+rlWwJ76dDg4kGZnpf19\naWam7tYAjVfLnjogKZRcNjZCoG9shNcAakGoo5x8DX1m5rgUQ7ADtaD8gnI4+wUYi1HLL4Q6AEwg\nauoAAEIdAFJCqANAQgh1AEgIoQ4ACSHUASAhhDoAJIRQB4CEEOoAkBBCHQASQqgDQEIIdQBICKEO\nAAkh1AEgIYQ6ACSEUAeAhBDqAJAQQh0AEkKoA0BCCHUASEjpUDezh/HncvnmNE+n06m7CWOV8val\nvG0S2zetqthTXzaz/5G0W8GyGif1P6yUty/lbZPYvml1qYJlLLv7ZgXLAQCUVMWe+lUzWzCz1QqW\nBQAowdy9mgWZrUt67u7bPcOrWQEATBl3t2HnGVh+6XMA9NDdN81sRdJPsfzyk6Q5SSdCfZRGAQBG\nMzDU3f27M0bvSnoVn1+T9LyKRgEARlOqph5LLYtmdkvS/7r7u2qahYuWnZqae30rHitZPmtYUxRs\n36lTcZu8fUCm9IFSd9+Mjz/nh6f6AUkxDGIZ7Vbu9bx09KUtM7teNKyGpo6kd/uiE6fiNnz7luNj\nPTcsiS/lPtuWzGfQzJZiu7/NDSv13o3litImf0DOIZkwyLj7Y0l7uUG3Jb2Pz/ckLcZh3Z5hjVCw\nfVI4Ffe37v4yvi7a5olnZguSXsQy6Vz84F+Xmv+lXLRtcVQSn8G4PQux3XP93qdht29c3QQ08gNy\nTkmEwQAtSYe519f6DGuy3lNxm7p9czr+m9uLr+/o9BdwE7+Ue7dtNj5P4jPo7tvu/sf48qq7v1V4\n70rtUFVx8VGRpn5AzuNq/Iadd/cNpbutSZ+1lJ0AYGYf5vYAG7fNPScyzEv6XtIHCmejZRr5pVyw\nbX+Nz5P5DJrZZUkrkr6Jgy6r5A7VuEJdauAH5DxSCYMBupKuxuctHQdENuyKToZGo8S65GHPqbj5\nbW7c9sV/0V+7+1szkxL6m8xt2zsprc+gu/8sacPMtszsTRxcalvGVX5p9Aekn3jAJjvglkQY9PG9\nwrYp/nzeM2xWzT59dU/Si/j8mqT/VvO3b8Hd78XnRV/KTf47Pdq2lD6DZjafq4+/kbSkCt67cYV6\n0z8g/aQYBjKzJUk3zOxTSYq1vexATtfd3xUNq63BQyrYvlOn4jZ8+1ZiGSJrfzJfygXbltJncEEn\nA3xXFbx3lXUTcGrB4V/cPUlzAy5gapTcXsJsdhpnqtuKyWdmi5KeKNRcr0pacveXRX+TTfs7PWPb\nkvgMxnr67fhyLv/fiEq8d2MLdQDAxePORwCQEEIdABJCqANAQgh1AEgIoQ4ACSHUASAhhDomhpld\nzs5Bjl2Nlr7vbb7L1jOmGXpd+bYCk2Scfb8Aw7qm0EvdZuyXpZTYZ8h/DZpuxHUdtXWEeYGxIdQx\nSb7Q8eX7JumGwiXRXyp0Rzov6ZGkDxWuMFyOHVg9UuiZUNmw+Py2u38Zr0z8ot8yJP06t67PFfrW\nuBGnO1S4ii+7VP2Vu9/ItfUP7v5DURvMbC4uwxX671iOHTgBY0OoY5KsS7oSb2qeL224u9+Owz5z\n94/i8ztm9vs4/oaZtRTumfubgmX3XYZC/yGZWXf/JF7C/VohvM9q6w/xzkpFbbgl6ZW734v9llyV\nRKhjrAh1TJJ+XY5mXZL+rOM7GHUVOkGaV7hrzJM4/L10VHrZisO8zzJ+jsvIeyGFLlFjF7bnaesH\nkmZ72yDpsaR7ZrYV19nvCwKoDAdKMakG9SmdjX8t6Y2733b32wodQEmh9PKyeNahXZOOOpgq8qpf\nGyR97+4fKYT6SkXtAfoi1DFJDiXNx7JIvqc5z/088Tx3/8otM3uleN/KAv2WcRZX2HOfj3vbi7l5\nsrb+4Yw2vJL0NA77RNLTAesDSqOXRgBICHvqAJAQQh0AEkKoA0BCCHUASAihDgAJIdQBICGEOgAk\nhFAHgIT8P0Sa6CwsJriOAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def plot_gene(gene_name='YER124C'):\n", " plt.plot(data['t'], data['Y'][gene_name], 'rx')\n", " plt.title('Gene: ' + gene_name)\n", " plt.xlabel('time/minutes')\n", "plot_gene('YER124C')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our second data set is from ChiP-chip experiments performed on yeast by [Lee et al (2002)](http://www.cs.gsu.edu/~wkim/index_files/ref/TR.pdf). These give us the binding information between transcription factors and genes. In this notebook we are going to try and combine this binding information with the gene expression information to infer transcription factor activities." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Acquiring resource: lee_yeast_ChIP\n", "\n", "Details of data: \n", "Binding location analysis for 106 regulators in yeast. The data consists of p-values for binding of regulators to genes derived from ChIP-chip experiments.\n", "\n", "Please cite:\n", "Tong Ihn Lee, Nicola J. Rinaldi, Francois Robert, Duncan T. Odom, Ziv Bar-Joseph, Georg K. Gerber, Nancy M. Hannett, Christopher T. Harbison, Craig M. Thompson, Itamar Simon, Julia Zeitlinger, Ezra G. Jennings, Heather L. Murray, D. Benjamin Gordon, Bing Ren, John J. Wyrick, Jean-Bosco Tagne, Thomas L. Volkert, Ernest Fraenkel, David K. Gifford, Richard A. Young 'Transcriptional Regulatory Networks in Saccharomyces cerevisiae' Science 298 (5594) pg 799--804. DOI: 10.1126/science.1075090\n", "\n", "After downloading the data will take up 1674161 bytes of space.\n", "\n", "Data will be stored in /home/maxz/tmp/GPy-datasets/lee_yeast_ChIP.\n", "\n", "Do you wish to proceed with the download? [yes/no]\n", "yes\n", "binding_by_gene.tsv\n", "Downloading http://jura.wi.mit.edu/young_public/regulatory_network/binding_by_gene.tsv -> /home/maxz/tmp/GPy-datasets/lee_yeast_ChIP/binding_by_gene.tsv\n", "[==============================] 9.047/9.047MB\n" ] } ], "source": [ "data = GPy.util.datasets.lee_yeast_ChIP()\n", "# set S to find relationships where p-value is less than 1e-3\n", "S = data['Y'].T<1e-3" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(23, 6178)\n" ] } ], "source": [ "print Y.shape" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(113, 6270)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More details on the data are give in [this notebook](../datasets/leeYeastChip.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Matching the Data\n", "\n", "The two data sets have slightly different overlapping sets of genes. Fortunately, with `pandas` it is easy to do the matching." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "There are n= 6099 overlapping genes.\n" ] } ], "source": [ "intersect = list(set(S.columns) & set(Y.columns))\n", "# this is the new value for n\n", "print 'There are n=', len(intersect), 'overlapping genes.'\n", "# Now reset S and Y to have the overlap of genes\n", "Y = Y[intersect].T # now an n by T matrix\n", "S = S[intersect].T # now an n by q matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Model for Transcription Factor Activities\n", "\n", "We are working with *log* expression levels in a matrix $\\mathbf{Y} \\in \\Re^{n\\times T}$ and we will assume a linear (additive) model giving the relationship between the expression level of the gene and the corresponding transcription factor activity which are unobserved, but we represent by a matrix $\\mathbf{F} \\in \\Re^{q\\times T}$. Our basic assumption is as follows. Transcription factors are in time series, so they are likely to be temporally smooth. Further we assume that the transcription factors are potentially correlated with one another (to account for transcription factors that operate in unison). \n", "\n", "#### Correlation Between Transcription Factors \n", "If there are $q$ transcription factors then correlation between different transcription factors is encoded in a covariance matrix, $\\boldsymbol{\\Sigma}$ which is $q\\times q$ in dimensionality. \n", "\n", "#### Temporal Smoothness\n", "\n", "Further we assume that the log of the transcription factors' activities is temporally smooth, and drawn from an underlying Gaussian process with covariance $\\mathbf{K}_t$. \n", "\n", "#### Intrinsic Coregionalization Model\n", "\n", "We assume that the joint process across all $q$ transcription factor activities and across all time points is well represented by an *intrinsic model of coregionalization* where the covariance is given by the Kronecker product of these terms.\n", "$$\\mathbf{K}_f = \\mathbf{K}_t \\otimes \\boldsymbol{\\Sigma}$$\n", "\n", "This is known as an intrinsic coregionalization model [Wackernagel, (2003)](http://books.google.co.uk/books/about/Multivariate_Geostatistics.html?id=Rhr7bgLWxx4C). See [Alvarez et al (2012)](http://www.nowpublishers.com/articles/foundations-and-trends-in-machine-learning/MAL-036) for a machine learning orientated review of these methods. The matrix $\\boldsymbol{\\Sigma}$ is known as the coregionalization matrix.\n", "\n", "### Relation to Gene Expressions\n", "\n", "We now assume that the $j$th gene's expression is given by the product of the transcription factors that bind to that gene. Because we are working in log space, that implies a log linear relationship. At the $i$th time point, the log of the $j$th gene's expression, $\\mathbf{y}_{:,i}$ is linearly related to the log of the transcription factor activities at the corresponding time point, $\\mathbf{f}_{:, i}$. This relationship is given by the binding information from $\\mathbf{S}$. We then assume that there is some corrupting Gaussian noise to give us the final observation.\n", "\n", "$$\\mathbf{y}_{:, i} = \\mathbf{S}\\mathbf{f}_{:, i} + \\boldsymbol{\\epsilon}_i$$\n", "\n", "where the Gaussian noise is sampled from \n", "\n", "$$\\boldsymbol{\\epsilon}_i \\sim \\mathcal{N}(\\mathbf{0}, \\sigma^2 \\mathbf{I}).$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gaussian Process Model of Gene Expression\n", "\n", "We consider a vector operator which takes all the separate time series in $\\mathbf{Y}$ and stacks the time series to form a new vector $n\\times T$ length vector $\\mathbf{y}$. A similar operation is applied to form a $q \\times T$ length vector $\\mathbf{f}$. Using Kronecker products we can now represent the relationship between $\\mathbf{y}$ and $\\mathbf{f}$ as follows: \n", "$$\n", "\\mathbf{y} = \\left[\\mathbf{I} \\otimes \\mathbf{S}\\right] \\mathbf{f} + \\boldsymbol{\\epsilon}.\n", "$$\n", "[Standard properties of multivariate Gaussian distributions](../background/gaussianProperties.ipynb) tell us that \n", "\n", "$$\\mathbf{y} \\sim \\mathcal{N}(\\mathbf{0}, \\mathbf{K}),$$\n", "\n", "where\n", "$$\n", "\\mathbf{K} = \\mathbf{K}_t \\otimes \\mathbf{S} \\boldsymbol{\\Sigma} \\mathbf{S}^\\top + \\sigma^2 \\mathbf{I}.\n", "$$\n", "and the log likelihood of the model can be written as:\n", "$$\n", "L = -\\frac{1}{2} \\log |\\mathbf{K}| - \\frac{1}{2} \\mathbf{y}^\\top \\mathbf{K}^{-1} \\mathbf{y}.\n", "$$\n", "The covariance matrix in this log likelihood is of size $n$ by $T$ where $n$ is number of genes and $T$ is number of time points. For many experiments this covariance matrix is prohibitively large. It will require $\\mathcal{O}(n^3T^3)$ operations to invert and $\\mathcal{O}(n^2T^2)$ storage. We need to look for ways of representing the same covariance with alternative approximations.\n", "\n", "In the worst case, because the vector $\\mathbf{y}$ contains $T\\times n$ points ($T$ time points for each of $n$ genes) we are faced with $O(T^3n^3)$ computational complexity. We are going to use a rotation trick to help. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Improving the Computational Complexity\n", "\n", "However, we can get a drastic reduction in the size of the covariance function by considering the singular value decomposition of $\\mathbf{S}$. \n", "\n", "The matrix $\\mathbf{S}$ is $n$ by $q$ matrix, where $q$ is the number of transcription factors. It contains a 1 if a given transcription factor binds to a given gene, and zero otherwise.\n", "\n", "\n", "#### Rotating the Basis of a Multivariate Gaussian\n", "\n", "For any multivariate Gaussian you can rotate the data set and compute a new roated covariance which is valid for the rotated data set. Mathematically this works by first inserting $\\mathbf{R}\\mathbf{R}^\\top$ into the likelihood at three points as follows:\n", "\n", "$$L = -\\frac{1}{2} \\log |\\mathbf{K}\\mathbf{R}^\\top\\mathbf{R}| - \\frac{1}{2} \\mathbf{y}^\\top\\mathbf{R}^\\top\\mathbf{R} \\mathbf{K}^{-1}\\mathbf{R}^\\top\\mathbf{R} \\mathbf{y} + \\text{const}$$\n", "\n", "The rules of determinants and a transformation of the data allows us to rewrite the likelihood as\n", "\n", "$$L = -\\frac{1}{2} \\log |\\mathbf{R}^\\top\\mathbf{K}\\mathbf{R}| - \\frac{1}{2} \\hat{\\mathbf{y}}^\\top \\left[\\mathbf{R}^\\top\\mathbf{K}\\mathbf{R}\\right]^{-1}\\hat{\\mathbf{y}} + \\text{const}$$\n", "\n", "where we have introduced a rotated version of our data set: $\\hat{\\mathbf{y}}=\\mathbf{R} \\mathbf{y}$. Geometrically what this says is that if we want to maintain the same likelihood, then when we rotate our data set by $\\mathbf{R}$ we need to rotate either side of the covariance matrix by $\\mathbf{R}$, which makes perfect sense when we recall the properties of the multivariate Gaussian. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### A Kronecker Rotation\n", "\n", "In this notebook we are using a particular structure of covariance which involves a Kronecker product. The rotation we consider will be a Kronecker rotation (see [Stegle et al, 2011](http://papers.nips.cc/paper/4281-efficient-inference-in-matrix-variate-gaussian-models-with-iid-observation-noise.pdf)). We are going to try and take advantage of the fact that the matrix $\\mathbf{S}$ is square meaning that $\\mathbf{S}\\boldsymbol{\\Sigma}\\mathbf{S}^\\top$ is not full rank (it has rank of most $q$, but is size $n\\times n$, and we expect number of transcription factors $q$ to be less than number of genes $n$). \n", "\n", "When ranks are involved, it is always a good idea to look at singular value decompositions (SVDs). The SVD of $\\mathbf{S}$ is given by:\n", "$$\\mathbf{S} = \\mathbf{Q} \\boldsymbol{\\Lambda} \\mathbf{V}^\\top$$\n", "where $\\mathbf{V}^\\top \\mathbf{V} = \\mathbf{I}$, $\\boldsymbol{\\Lambda}$ is a diagonal matrix of positive values, $\\mathbf{Q}$ is a matrix of size $n\\times q$: it matches the dimensionality of $\\mathbf{S}$, but we have $\\mathbf{Q}^\\top \\mathbf{Q} = \\mathbf{I}$. Note that because it is not square, $\\mathbf{Q}$ is not in itself a rotation matrix. However it could be seen as the first $q$ columns of an $n$ dimensional rotation matrix (assuming $n$ is larger than $q$, i.e. there are more genes than transcription factors). \n", "\n", "If we call the $n-q$ missing columns of this rotation matrix $\\mathbf{U}$ then we have a valid rotation matrix $\\mathbf{R}=\\begin{bmatrix} \\mathbf{Q}& \\mathbf{U}\\end{bmatrix}$. Although this rotation matrix is only rotating across the $n$ dimensions of the genes, not the additional dimensions across time. In other words we are choosing $\\mathbf{K}_t$ to be unrotated. To represent this properly for our covariance we need to set $\\mathbf{R} = \\mathbf{I} \\otimes \\begin{bmatrix} \\mathbf{Q}& \\mathbf{U}\\end{bmatrix}$. This gives us a structure that when applied to a covariance of the form $\\mathbf{K}_t\\otimes \\mathbf{K}_n$ it will rotate $\\mathbf{K}_n$ whilst leaving $\\mathbf{K}_t$ untouched.\n", "\n", "When we apply this rotation matrix to $\\mathbf{K}$ we have to consider two terms, the rotation of $\\mathbf{K}_t \\otimes \\mathbf{S}\\boldsymbol{\\Sigma}\\mathbf{S}^\\top$, and the rotation of $\\sigma^2 \\mathbf{I}$.\n", "\n", "Rotating the latter is easy, because it is just the identity multiplied by a scalar so it remains unchanged\n", "$$\n", "\\mathbf{R}^\\top\\mathbf{I}\\sigma^2 \\mathbf{R}= \\mathbf{I}\\sigma^2\n", "$$\n", "The former is slightly more involved, for that term we have\n", "$$\n", "\\left[\\mathbf{I}\\otimes \\begin{bmatrix}\\mathbf{Q} & \\mathbf{U}\\end{bmatrix}^\\top \\right]\\mathbf{K}_t \\otimes \\mathbf{S}\\boldsymbol{\\Sigma}\\mathbf{S}^\\top\\left[ \\mathbf{I} \\otimes \\begin{bmatrix}\\mathbf{Q} & \\mathbf{U}\\end{bmatrix}\\right]=\\mathbf{K}_t \\otimes \\begin{bmatrix}\\mathbf{Q} & \\mathbf{U}\\end{bmatrix}^\\top \\mathbf{S} \\boldsymbol{\\Sigma}\\mathbf{S}^\\top \\begin{bmatrix}\\mathbf{Q} & \\mathbf{U}\\end{bmatrix}.\n", "$$ \n", "Since $\\mathbf{S} = \\mathbf{Q}\\boldsymbol{\\Lambda}\\mathbf{V}^\\top$ then we have\n", "$$\n", "\\begin{bmatrix}\\mathbf{Q} & \\mathbf{U}\\end{bmatrix}^\\top \\mathbf{X}\\boldsymbol{\\Sigma}\\mathbf{X}^\\top\\begin{bmatrix}\\mathbf{Q} & \\mathbf{U}\\end{bmatrix} = \\begin{bmatrix}\\boldsymbol{\\Lambda} \\mathbf{V}^\\top \\boldsymbol{\\Sigma}\\mathbf{V} \\boldsymbol{\\Lambda} &\\mathbf{0} \\\\ \\mathbf{0} & \\mathbf{0}\\end{bmatrix}.\n", "$$\n", "This prompts us to split our vector $\\hat{\\mathbf{y}}$ into a $q$ dimensional vector $\\hat{\\mathbf{y}}_q =\\mathbf{Q}^\\top \\mathbf{y}$ and an $n-q$ dimensional vector $\\hat{\\mathbf{y}}_u = \\mathbf{U}^\\top \\mathbf{y}$. The Gaussian likelihood can be written as\n", "$$\n", "L = L_u + L_q + \\text{const}\n", "$$\n", "where\n", "$$\n", "L_q = -\\frac{1}{2} \\log |\\mathbf{K}_t\\otimes\\boldsymbol{\\Lambda}\\mathbf{V}^\\top\\boldsymbol{\\Sigma}\\mathbf{V}\\boldsymbol{\\Lambda}+\\sigma^2\\mathbf{I}| - \\frac{1}{2} \\hat{\\mathbf{y}}_q^\\top \\left[\\mathbf{K}_t\\otimes \\boldsymbol{\\Lambda}\\mathbf{V}^\\top\\boldsymbol{\\Sigma}\\mathbf{V}\\boldsymbol{\\Lambda}+\\sigma^2\\mathbf{I}\\right]^{-1} \\hat{\\mathbf{y}}_q\n", "$$\n", "and\n", "$$\n", "L_u = -\\frac{T(n-q)}{2} \\log \\sigma^2 -\\frac{1}{2\\sigma^2} \\hat{\\mathbf{y}}_u^\\top \\hat{\\mathbf{y}}_u\n", "$$\n", "Strictly speaking we should fit these models jointly, but for the purposes of illustration we will firstly use a simple procedure. Firstly, we fit the noise variance $\\sigma^2$ on $\\hat{\\mathbf{y}}_u$ alone using $L_u$. Once this is done, fix the value of $\\sigma^2$ in $L_q$ and optimize with respect to the other parameters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inferring the Transcription Factors\n", "\n", "The model for $\\hat{y}_q$ assumes that we have some latent function $\\mathbf{f}$ which is sampled according to the covariance $\\mathbf{K}_f = \\mathbf{K}_t \\otimes \\mathbf{\\Sigma}$ which is then multiplied by $\\mathbf{I}\\otimes \\mathbf{V}\\boldsymbol{\\Lambda}$ and corrupted with Gaussian noise. To recover predictions for $\\mathbf{f}$ we first define $\\mathbf{g} = \\mathbf{I}\\otimes \\mathbf{V}\\boldsymbol{\\Lambda} \\mathbf{f}$, then we have that posterior predictions for $\\mathbf{g}$ may be had through,\n", "$$\n", "g|\\hat{y} \\sim \\mathcal{N}(\\boldsymbol{\\mu}_g, \\boldsymbol{\\Sigma}_g)\n", "$$\n", "where\n", "$$\n", "\\mu_g = \\mathbf{K}_*\\otimes\\boldsymbol{\\Lambda}\\mathbf{V}^\\top \\boldsymbol{\\Sigma}\\mathbf{V}\\boldsymbol{\\Lambda} \\left[\\mathbf{K}_t\\otimes\\boldsymbol{\\Lambda}\\mathbf{V}^\\top \\boldsymbol{\\Sigma}\\mathbf{V}\\boldsymbol{\\Lambda} + \\sigma^2 \\mathbf{I}\\right]^{-1} \\hat{\\mathbf{y}}_q\n", "$$\n", "and\n", "$$\n", "\\mathbf{C}_g = \\mathbf{K}_{*,*}\\otimes\\boldsymbol{\\Lambda}\\mathbf{V}^\\top \\boldsymbol{\\Sigma}\\mathbf{V}\\boldsymbol{\\Lambda} - \\mathbf{K}_*\\otimes\\boldsymbol{\\Lambda}\\mathbf{V}^\\top \\boldsymbol{\\Sigma}\\mathbf{V}\\boldsymbol{\\Lambda} \\left[\\mathbf{K}_t\\otimes\\boldsymbol{\\Lambda}\\mathbf{V}^\\top \\boldsymbol{\\Sigma}\\mathbf{V}\\boldsymbol{\\Lambda} + \\sigma^2 \\mathbf{I}\\right]^{-1} \\mathbf{K}_*\\otimes\\boldsymbol{\\Lambda}\\mathbf{V}^\\top \\boldsymbol{\\Sigma}\\mathbf{V}\\boldsymbol{\\Lambda}\n", "$$\n", "We can then easily find the values for the latent transcription factors,\n", "$$\n", "\\mathbf{f} = \\mathbf{I}\\otimes \\mathbf{V}\\boldsymbol{\\Lambda}^{-1}\\mathbf{g} .\n", "$$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# step 1, find the SVD of S.\n", "n, q = S.shape\n", "T = Y.shape[1]\n", "R, Lambda, V = scipy.linalg.svd(S)\n", "# Extract first q columns for Q\n", "Q = R[:, :q]\n", "# remaining columns for U\n", "U = R[:, q:]\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(6099, 113)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S.shape" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(6099, 113)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q.shape\n", "#U.shape\n", "#Y.shape" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sigma2 found as 0.178461463478\n" ] } ], "source": [ "# Find sigma2 by looking at variance of y_u\n", "Y_u = np.dot(U.T, Y)\n", "sigma2 = 1./(T*(n-q))*(Y_u*Y_u).sum()\n", "print \"sigma2 found as\", sigma2\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Prepare the data for processing in GPy\n", "Y_q = np.dot(Q.T, Y) # project data onto the principal subspace of X \n", "\n", "# Generate the input associated with each Y, the TF and the time point.\n", "x0, x1 = np.asarray(np.meshgrid(t.flatten(),np.arange(q)))\n", "X = np.hstack([x0.flatten()[:, None], x1.flatten()[:, None]])\n", "y = Y_q.flatten()[:, None]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "\n", "kern = GPy.kern.RBF(1, active_dims=[0])*GPy.kern.Coregionalize(1,q,rank=5, active_dims=[1])\n", "m = GPy.models.GPRegression(X, y, kern)\n", "m.mul.rbf.lengthscale = 50\n", "m.Gaussian_noise.variance = sigma2\n", "#m.Gaussian_noise.variance.constrain_fixed()\n", "#()#kern, t, y_q)\n", "#m.sigma2 = constrained_fixed(sigma2)\n", "#m.optimize(messages=True)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "

\n", "Model: GP regression
\n", "Objective: 3113.60415866
\n", "Number of Parameters: 681
\n", "Number of Optimization Parameters: 681
\n", "Updates: True
\n", "

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
GP_regression. valueconstraintspriors
mul.rbf.variance 1.0 +ve
mul.rbf.lengthscale 50.0 +ve
mul.coregion.W (113, 5)
mul.coregion.kappa (113,) +ve
Gaussian_noise.variance0.178461463478 +ve
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(m)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "_ = m.optimize()\n", "#m.plot(fixed_inputs=[(1, 1)]) # this would plot ACE2." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "

\n", "Model: GP regression
\n", "Objective: 1150.97215161
\n", "Number of Parameters: 681
\n", "Number of Optimization Parameters: 681
\n", "Updates: True
\n", "

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
GP_regression. valueconstraintspriors
mul.rbf.variance 0.0225223175544 +ve
mul.rbf.lengthscale 5.76167717687 +ve
mul.coregion.W (113, 5)
mul.coregion.kappa (113,) +ve
Gaussian_noise.variance0.0817952261337 +ve
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the current design the model is switching off the temporal correlation. The next step in the analysis will be to reimplement the same model as described by [Sanguinetti et al (2006)](http://bioinformatics.oxfordjournals.org/content/22/14/1753.short) and recover their results. That will involve using an Ornstein Uhlbeck covariance and joint maximisation of the likelihoood of $L_u$ and $L_q$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Things to do:\n", "1. Code the model such that $L_q$ and $L_u$ are jointly optimized.\n", "2. Recover the true latent functions associated with the transcription factor activities. \n", "3. Implement the same model with the OU covariance function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This work is funded by a Commonwealth Scholarship to Arif Rahman and by the EU FP7 RADIANT project." ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "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.11" } }, "nbformat": 4, "nbformat_minor": 0 }