ExtDataFig1c.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sys\n",
    "sys.path.append('..')\n",
    "from inference.simulation import simBetaPoisson\n",
    "from scipy import special\n",
    "from scipy.stats import poisson,norm\n",
    "from scipy.special import j_roots\n",
    "from scipy.special import beta as beta_fun  \n",
    "from joblib import delayed,Parallel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dBP_dict(at, alpha, bet, lam):\n",
    "    at_ = np.copy(at)\n",
    "    at_.shape = (len(at_), 1)\n",
    "    np.repeat(at_, 40, axis = 1)\n",
    "    def fun(at_, m):\n",
    "        if(max(m) < 1e6):\n",
    "            return(poisson.pmf(at_,m))\n",
    "        else:\n",
    "            return(norm.pdf(at_,loc=m,scale=sqrt(m)))\n",
    "    x,w = j_roots(40,alpha = bet - 1, beta = alpha - 1)\n",
    "    gs = np.sum(w*fun(at_, m = lam*(1+x)/2), axis=1)\n",
    "    prob = 1/beta_fun(alpha, bet)*2**(-alpha-bet+1)*gs\n",
    "    return {k: v for k, v in zip(at.tolist(), prob.tolist())}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getGOF(rna,gene):\n",
    "    from collections import Counter\n",
    "    if np.sum(rna) == 0:\n",
    "        return np.nan\n",
    "    try:\n",
    "        d = dBP_dict(rna, gene[0], gene[1], gene[2])\n",
    "        for key in d:\n",
    "            d[key] *= len(rna)\n",
    "        expected = d\n",
    "        observed = Counter(rna)\n",
    "        chi_sq = sum((((observed[k]-expected[k])**2)/expected[k]) for k in expected)\n",
    "        return chi_sq\n",
    "    except (ValueError, AssertionError):\n",
    "        return np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MCMC(gene):\n",
    "    p = np.array(params.loc[gene][0])\n",
    "    realGOF = getGOF(np.array(UMI.loc[gene][~np.isnan(UMI.loc[gene])]), p)\n",
    "    larger = np.zeros(1000)\n",
    "    for i in np.arange(1000):\n",
    "        rvs = simBetaPoisson(p)\n",
    "        if realGOF < getGOF(rvs, p):\n",
    "                 larger[i] = 1\n",
    "    return np.mean(larger)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "UMI = pd.read_csv('data/SS3_c57_UMIs_concat.csv', index_col=0)\n",
    "params = pd.read_pickle('data/SS3_c57_UMIs_concat_ML.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params = params[params[1]][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params = pd.DataFrame(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=50)]: Done  28 tasks      | elapsed:    2.4s\n",
      "[Parallel(n_jobs=50)]: Done 188 tasks      | elapsed:    9.9s\n",
      "[Parallel(n_jobs=50)]: Done 412 tasks      | elapsed:   21.4s\n",
      "[Parallel(n_jobs=50)]: Done 700 tasks      | elapsed:   36.2s\n",
      "[Parallel(n_jobs=50)]: Done 1052 tasks      | elapsed:   54.8s\n",
      "[Parallel(n_jobs=50)]: Done 1468 tasks      | elapsed:  1.3min\n",
      "[Parallel(n_jobs=50)]: Done 1948 tasks      | elapsed:  1.7min\n",
      "[Parallel(n_jobs=50)]: Done 2492 tasks      | elapsed:  2.2min\n",
      "[Parallel(n_jobs=50)]: Done 3100 tasks      | elapsed:  2.7min\n",
      "[Parallel(n_jobs=50)]: Done 3772 tasks      | elapsed:  3.3min\n",
      "[Parallel(n_jobs=50)]: Done 4508 tasks      | elapsed:  4.0min\n",
      "[Parallel(n_jobs=50)]: Done 5308 tasks      | elapsed:  4.7min\n",
      "[Parallel(n_jobs=50)]: Done 6172 tasks      | elapsed:  5.5min\n",
      "[Parallel(n_jobs=50)]: Done 7100 tasks      | elapsed:  6.3min\n",
      "[Parallel(n_jobs=50)]: Done 7382 out of 7382 | elapsed:  6.5min finished\n"
     ]
    }
   ],
   "source": [
    "GOF_C57 = Parallel(n_jobs=50, verbose = 3)(delayed(MCMC)(gene) for gene in params.index.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "gof_series = pd.Series(GOF_C57, index = params.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "me = np.log10(pd.Series([p[2]*(p[0]/(p[0]+p[1])) for p in params[0]], index=params.index))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuUHVWZ9/Hvz5DQhGsIAZEgDSNC\nlJDQNDcRhpuOZBTIDAojCCKKIGrAV5E4g7B4ZQRFIbwCIwJDGKOAGUmCxAsEkCALNAkNCSSEgCAd\nQugECEQSuT3vH7U7HJM+neqcW3X377PWWV21q86u51Rfnt5Vu/ZWRGBmZlY072p0AGZmZl1xgjIz\ns0JygjIzs0JygjIzs0JygjIzs0JygjIzs0JygjIzs0JygjIzs0JygjIzs0LaqNEB1MI222wTzc3N\njQ7DzMxKzJ49e1lEDMu7f59MUM3NzcyaNavRYZiZWQlJz/Rkf1/iMzOzQnKCMjOzQnKCMjOzQuqT\n96DMzGrhjTfeoL29ndWrVzc6lEJrampi+PDhDBw4sKJ6nKDMzHJqb29n8803p7m5GUmNDqeQIoLl\ny5fT3t7OzjvvXFFdvsRnZpbT6tWrGTp0qJNTNyQxdOjQqrQynaDMzHrAyWn9qnWOnKDMzKyQfA/K\nzGwDTXhpQlXrGzdk3Hr3Wbp0KWeffTYPPPAAQ4YMYdCgQZxzzjmMHTu24uMfcsghXHrppbS2tv5d\n+cyZMzn99NMZOHAgt99+O+PGjWPy5Mm0tbXx3HPPMWbMmIqP3RW3oMy6I1X+MquSiOCYY47h4IMP\n5qmnnmL27NncdNNNtLe31/S4kyZNYvz48bS1tbHDDjswefJkANra2pg+fXrNjusEZWbWS9x1110M\nGjSI008/fU3ZTjvtxFe+8hUg68RxyimnMHLkSPbaay/uvvvubstXrVrF8ccfz4gRIxg7diyrVq1a\n55jXXnstt9xyC+eddx4nnHACTz/9NHvssQevv/463/72t7n55psZPXo0N998c9U/ry/xmZn1Eo8+\n+igtLS1lt1955ZVIYu7cuSxYsICPfvSjLFy4sGz51VdfzeDBg5k/fz6PPPJIl3V//vOf57777uPj\nH/84xx57LE8//TQAgwYN4sILL2TWrFn86Ec/qsnndQvKrNZ8idBq5Mwzz2TUqFHss88+ANx3332c\neOKJAOy+++7stNNOLFy4sGz5vffeu6Z8zz33ZM8992zMBynDCcrMrJf44Ac/yJw5c9asX3nllcyY\nMYOOjo4GRlU7TlBmZr3EYYcdxurVq7n66qvXlL322mtrlg866CAmTZoEwMKFC/nLX/7CbrvtVrb8\n4IMP5mc/+xkA8+bN45FHHulRPJtvvjmvvvpqpR+rLN+DMjPbQHm6hVeTJKZMmcLZZ5/N9773PYYN\nG8amm27KJZdcAsCXvvQlzjjjDEaOHMlGG23EDTfcwMYbb1y2/IwzzuCUU05hxIgRjBgxgr333rtH\n8Rx66KFcfPHFjB49mvHjx3PcccdV9/NGRFUrLILW1tbwhIVWFUW4B9QHf0d7q/nz5zNixIhGh9Er\ndHWuJM2OiNYyb1mHL/GZmVkhOUGZmVkhOUGZmVkhOUGZmVkh1SxBSbpe0guS5pWUbS3pDklPpK9D\nUrkkXSFpkaRHJLWUvOfktP8Tkk6uVbxmZlYstWxB3QB8bK2yc4EZEbErMCOtAxwJ7JpepwFXQ5bQ\ngPOB/YB9gfM7k5qZmfVtNUtQEXEv8OJaxUcDE9PyROCYkvIbI/MAsJWk7YF/Au6IiBcj4iXgDtZN\nemZmjVGN0e57OKzVgAEDGD16NKNGjaKlpYX777+/RyFfcMEFXHrppeuUd3R0sN9++7HXXnsxc+ZM\nxowZw8svv8zLL7/MVVdd1aNjVEu970FtFxFL0vLzwHZpeQfg2ZL92lNZufJ1SDpN0ixJs/rqsB9m\nZptssgltbW08/PDDfPe732X8+PFVqXfGjBmMHDmShx56iIMOOojp06ez1VZb9asEtUZkTwhX7QnE\niLgmIlojonXYsGHVqtbMrLBeeeUVhgzJ7nqsXLmSww8/nJaWFkaOHMnUqVPX7HfRRRfx/ve/nw9/\n+MM8/vjj69TT1tbGOeecw9SpUxk9ejSrVq2iubmZZcuWce655/Lkk08yevRovvGNb9Tts0H9hzpa\nKmn7iFiSLuG9kMoXAzuW7Dc8lS0GDlmr/J46xGlmVkirVq1i9OjRrF69miVLlnDXXXcB0NTUxK23\n3soWW2zBsmXL2H///TnqqKOYM2cON910E21tbbz55pu0tLSsM6TR6NGjy06dcfHFFzNv3jza2trq\n9hk71bsFNQ3o7Il3MjC1pPyk1Jtvf2BFuhT4W+CjkoakzhEfTWVmZv1S5yW+BQsW8Jvf/IaTTjqJ\niCAi+Na3vsWee+7JEUccweLFi1m6dCkzZ85k7NixDB48mC222IKjjjqq0R8ht5q1oCT9nKz1s42k\ndrLeeBcDt0g6FXgG+FTafTowBlgEvAacAhARL0r6v8Cf0n4XRsTaHS/MzPqlAw44gGXLltHR0cH0\n6dPp6Ohg9uzZDBw4kObmZlavXt3oECtSy158/xYR20fEwIgYHhHXRcTyiDg8InaNiCM6k03qvXdm\nRPxDRIyMiFkl9VwfEe9Lr/+uVbxmZr3NggULeOuttxg6dCgrVqxg2223ZeDAgdx9990888wzABx8\n8MFMmTKFVatW8eqrr3Lbbbf16Bi1nlKjO55uw8xsQzVgpPnOe1DZ4YOJEycyYMAATjjhBD7xiU8w\ncuRIWltb2X333QFoaWnhuOOOY9SoUWy77bZrZt/Na+jQoRx44IHsscceHHnkkXz/+9+v+mcqx9Nt\nmHXH021YCU+3kV81pttwC8qsGxNevLziOsZtfVYVIjHrfzxYrJmZFZITlJlZD/TF2yLVVq1z5ARl\nZpZTU1MTy5cvd5LqRkSwfPlympqaKq7L96DMaqzS+1jjqhSHVW748OG0t7fj8T6719TUxPDhwyuu\nxwnKzCyngQMHsvPOOzc6jH7Dl/jMzKyQnKDMzKyQnKDMzKyQnKDMzKyQnKDMzKyQnKDMzKyQnKDM\nzKyQnKDMzKyQnKDMzKyQnKCsb5Mqe5lZwzhBmZlZITlBmZlZITlBmZlZIXk0c7Oiq/RemOcusl7K\nLSgzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyskJygzMyukhiQo\nSWdLelTSPEk/l9QkaWdJD0paJOlmSYPSvhun9UVpe3MjYjYzs/qqe4KStAPwVaA1IvYABgDHA5cA\nl0XE+4CXgFPTW04FXkrll6X9zMysj2vUWHwbAZtIegMYDCwBDgM+nbZPBC4ArgaOTssAk4EfSVKE\nBxiz9Zvw4uWNDsHMNlDdW1ARsRi4FPgLWWJaAcwGXo6IN9Nu7cAOaXkH4Nn03jfT/kPrGbOZmdVf\nIy7xDSFrFe0MvAfYFPhYFeo9TdIsSbM6Ojoqrc7MzBqsEZ0kjgD+HBEdEfEG8EvgQGArSZ2XHIcD\ni9PyYmBHgLR9S2D52pVGxDUR0RoRrcOGDav1ZzAzsxprRIL6C7C/pMGSBBwOPAbcDRyb9jkZmJqW\np6V10va7fP/JzKzv67aThKQDgBOBg4DtgVXAPOB24KcRsaKnB4yIByVNBuYAbwIPAdekOm+S9J1U\ndl16y3XA/0haBLxI1uPPzMz6uLIJStKvgefIWjIXAS8ATcD7gUOBqZJ+GBHTenrQiDgfOH+t4qeA\nfbvYdzXwyZ4ew8zMerfuWlCfiYhla5WtJGv5zAF+IGmbmkVmZmb9Wtl7UJ3JSdKmkt6Vlt8v6ShJ\nA0v3MTMzq7Y8nSTuBZrSCBC/Az4D3FDLoMzMzPIkKEXEa8C/AFdFxCeBD9Y2LDMz6+9yJajUm+8E\nsp52kI2fZ2ZmVjN5EtRZwHjg1oh4VNIuZM8smZmZ1cx6B4uNiN8Dv5c0OK0/RTYauZnVQaUD3o6r\nUhxm9bbeFpSkAyQ9BixI66MkXVXzyMzMrF/Lc4nvcuCfSOPfRcTDwMG1DMrMzCzXWHwR8exaRW/V\nIBYzM7M18kxY+KykDwGRHtAdB8yvbVhmZtbf5WlBnQ6cSTZx4GJgdFo3MzOrmTy9+JaRPQNlZmZW\nN+tNUJKGAV8Amkv3j4jP1S4sMzPr7/Lcg5oKzATuxJ0jzMysTvIkqMER8c2aR2JmZlYiTyeJX0ka\nU/NIzMzMSuRJUOPIktRqSa9IelXSK7UOzMzM+rc8vfg2r0cgZmZmpfKMxSdJJ0o6L63vKGnf2odm\nZmb9WZ5LfFcBBwCfTusrgStrFpGZmRn5evHtFxEtkh4CiIiXJA2qcVxmZtbP5WlBvSFpABCw5sHd\nt2salZmZ9Xt5EtQVwK3AtpIuAu4D/rOmUZmZWb+XpxffJEmzgcMBAcdEhEczNzOzmsozFt/WwAvA\nz0vKBkbEG7UMzMyqRKq8jojK6zDroTyX+OYAHcBC4Im0/LSkOZL2rmVwZmbWf+VJUHcAYyJim4gY\nChwJ/Ar4ElkXdLPakSp7mVmvlSdB7R8Rv+1ciYjfAQdExAPAxjWLzMzM+rU8z0EtkfRN4Ka0fhyw\nNHU9d3dzMzOriTwtqE8Dw4EpZN3Nd0xlA4BP1S40MzPrz/JO+f6VMpsXVTccMzOzTNkWlKSfSBpZ\nZtumkj4n6YQNOaikrSRNlrRA0nxJB0jaWtIdkp5IX4ekfSXpCkmLJD0iqWVDjmlmZr1Ld5f4rgTO\nSwnkF5KuknS9pJnA/cDmwOQNPO4E4DcRsTswCpgPnAvMiIhdgRlpHbJeg7um12nA1Rt4TDMz60XK\nXuKLiDbgU5I2A1qB7YFVwPyIeHxDDyhpS+Bg4LPpOK8Dr0s6Gjgk7TYRuAf4JnA0cGNEBPBAan1t\nHxFLNjQG6z0mvHh5o0MwswbJcw9qJVmyqJadyR72/W9Jo4DZZLP2bleSdJ4HtkvLOwDPlry/PZX9\nXYKSdBpZC4v3vve9VQzXzMwaIU8vvmrbCGgBro6IvYC/8s7lPABSa6lHY6tExDUR0RoRrcOGData\nsGZm1hiNSFDtQHtEPJjWJ5MlrKWStgdIX19I2xeTdW3vNDyVmZlZH5Znyvcue/JtqIh4HnhW0m6p\n6HDgMWAacHIqOxmYmpanASel3nz7Ayt8/8nMrO/LM5LEVZI2Bm4AJkXEiioc9yvApDQz71PAKWTJ\n8hZJpwLP8M5DwNOBMWTPXL2W9jUzsz4uTyeJgyTtCnwOmC3pj8B/R8QdG3rQ1EOwtYtNh3exbwBn\nbuixzMysd8p1DyoingD+g6zb9z8CV6SHbP+llsGZmVn/lece1J6SLiN7mPYw4BMRMSItX1bj+MzM\nrJ/Kcw/q/wHXAt+KiFWdhRHxnKT/qFlkZmbWr+VJUP8MrIqItwAkvQtoiojXIuJ/ahqdmZn1W3nu\nQd0JbFKyPjiVmZmZ1UyeFlRTGu4IyIY+kjS4hjGZWRVVYzzDcVWIw6yn8rSg/lo6xYWkvckGjTUz\nM6uZPC2os4BfSHoOEPBusmnfzczMaibPg7p/krQ70Dk00eMR8UZtwzIzs/4uTwsKYB+gOe3fIomI\nuLFmUZmZWb+33gQl6X+AfwDagLdScQBOUGZmVjN5WlCtwAfSmHhmZmZ1kacX3zyyjhFmZmZ1k6cF\ntQ3wWBrF/G+dhRFxVM2iMjOzfi9Pgrqg1kGYmZmtLU83899L2gnYNSLuTKNIDKh9aGZm1p/lmW7j\nC8Bk4MepaAdgSi2DMjMzy9NJ4kzgQOAVWDN54ba1DMrMzCxPgvpbRLzeuSJpI7LnoMzMzGomT4L6\nvaRvAZtI+gjwC+C22oZlZmb9XZ4EdS7QAcwFvghMBzyTrpmZ1VSeXnxvAz9JLzMzs7rIMxbfn+ni\nnlNE7FKTiMzMzMg/Fl+nJuCTwNa1CcfMzCyz3ntQEbG85LU4Ii4H/rkOsVlvJ1X+MrN+K88lvpaS\n1XeRtajyziNlZma2QfIkmh+ULL8JPA18qibRmJmZJXl68R1aj0DMzMxK5bnE97XutkfED6sXjpmZ\nWSZvL759gGlp/RPAH4EnahWU9Q0TXry80SGYWS+WJ0ENB1oi4lUASRcAt0fEibUMzMzM+rc8Qx1t\nB7xesv56KquIpAGSHpL0q7S+s6QHJS2SdLOkQal847S+KG1vrvTYZmZWfHkS1I3AHyVdkFpPDwIT\nq3DsccD8kvVLgMsi4n3AS8CpqfxU4KVUflnaz8zM+rg8D+peBJxCljReAk6JiP+s5KCShpM97Htt\nWhdwGNnEiJAlwGPS8tG8kxAnA4en/c3MrA/L04ICGAy8EhETgHZJO1d43MuBc4C30/pQ4OWIeDOt\nt5PN3Ev6+ixA2r4i7f93JJ0maZakWR0dHRWGZ2ZmjZZnyvfzgW8C41PRQOCnG3pASR8HXoiI2Rta\nR1ci4pqIaI2I1mHDhlWzajMza4A8vfjGAnsBcwAi4jlJm1dwzAOBoySNIRt8dgtgArCVpI1SK2k4\nsDjtvxjYkazlthGwJbC8guObmVkvkOcS3+sREaQpNyRtWskBI2J8RAyPiGbgeOCuiDgBuBs4Nu12\nMjA1LU9L66Ttd6V4zMysD8uToG6R9GOyFs4XgDupzeSF3wS+JmkR2T2m61L5dcDQVP41shl+zcys\nj8szFt+lkj4CvALsBnw7Iu6oxsEj4h7gnrT8FLBvF/usJpuDyszM+pFuE5SkAcCdacDYqiQlMzOz\nPLq9xBcRbwFvS9qyTvGYmZkB+XrxrQTmSroD+GtnYUR8tWZRmVmxVPpsvPs12QbIk6B+mV5mZmZ1\nUzZBSXpvRPwlIqox7p6ZmVmPdHcPakrngqT/rUMsZmZma3SXoEovOu9S60DMzMxKdZegosyymZlZ\nzXXXSWKUpFfIWlKbpGXSekTEFjWPzszM+q2yCSoiBtQzEDMzs1J554MyMzOrKycoK0+q7GVmVoE8\nD+paPzXhxcsbHYKZ9WNuQZmZWSE5QZmZWSH5Ep+ZrVell3vHVSkO61/cgjIzs0JygjIzs0JygjIz\ns0JygjIzs0JygjIzs0JygjIzs0JygjIzs0JygjIzs0Lyg7p92ISXJjQ6BDOzDeYWlJmZFZITlJmZ\nFZITlJmZFZLvQZlZ7VU6gWVEdeKwXsUtKDMzKyQnKDMzK6S6JyhJO0q6W9Jjkh6VNC6Vby3pDklP\npK9DUrkkXSFpkaRHJLXUO2YzM6u/RrSg3gT+T0R8ANgfOFPSB4BzgRkRsSswI60DHAnsml6nAVfX\nP2QzM6u3uieoiFgSEXPS8qvAfGAH4GhgYtptInBMWj4auDEyDwBbSdq+zmGbmVmdNfQelKRmYC/g\nQWC7iFiSNj0PbJeWdwCeLXlbeypbu67TJM2SNKujo6NmMZuZWX00LEFJ2gz4X+CsiHildFtEBNCj\nfqURcU1EtEZE67Bhw6oYqZmZNUJDEpSkgWTJaVJE/DIVL+28dJe+vpDKFwM7lrx9eCozM7M+rO4P\n6koScB0wPyJ+WLJpGnAycHH6OrWk/MuSbgL2A1aUXAo0s15gwouXV/T+cVWKw3qXRowkcSDwGWCu\npLZU9i2yxHSLpFOBZ4BPpW3TgTHAIuA14JT6hmtmZo1Q9wQVEfcB5cY9ObyL/QM4s6ZBmZlZ4Xgk\nCTMzKyQnKDMzKyQnKDMzKyQnKDMzKyQnKDMzKyRPWFhQE16a0OgQzMwayi0oMzMrJCcoMzMrJF/i\nM7PiU7ln+3sgejT+tBWAW1BmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZI7sVXRqUPyo4b4inW\nzKql0gkPwZMe9kZuQZmZWSE5QZmZWSE5QZmZWSE5QZmZWSE5QZmZWSG5F1+NeLoMs4KpdDw/j+VX\nd25BmZlZITlBmZlZIfkSn5n1C5U+7OsHfevPLSgzMyskJygzMyskJygzMyskJygzMyskJygzMysk\n9+IzM8vDD/rWnVtQZmZWSG5BmZnl4Oeo6q/XtKAkfUzS45IWSTq30fGYmVlt9YoWlKQBwJXAR4B2\n4E+SpkXEY42NzMwsn0oHkB43pP+1wXpFggL2BRZFxFMAkm4CjgacoMysXyjCDAn1TpK9JUHtADxb\nst4O7Fe6g6TTgNPS6kpJj9cptjy2AZY1OohuOL7KOL7KOL7K1C2+szhrQ95WGt9OPXljb0lQ6xUR\n1wDXNDqOrkiaFRGtjY6jHMdXGcdXGcdXmb4cX2/pJLEY2LFkfXgqMzOzPqq3JKg/AbtK2lnSIOB4\nYFqDYzIzsxrqFZf4IuJNSV8GfgsMAK6PiEcbHFZPFPLSYwnHVxnHVxnHV5k+G5/Cw2+YmVkB9ZZL\nfGZm1s84QZmZWSE5QdWApE9KelTS25LKdq+U9LSkuZLaJM0qYHwNGV5K0taS7pD0RPo6pMx+b6Vz\n1yap5p1m1nc+JG0s6ea0/UFJzbWOqYfxfVZSR8k5+3wdY7te0guS5pXZLklXpNgfkdRSr9hyxneI\npBUl5+7bdY5vR0l3S3os/e6u88RsI89hzvh6fg4jwq8qv4ARwG7APUBrN/s9DWxTxPjIOqM8CewC\nDAIeBj5Qp/i+B5ybls8FLimz38o6nrP1ng/gS8B/peXjgZsLFt9ngR/V++ctHftgoAWYV2b7GODX\ngID9gQcLFt8hwK8ace7S8bcHWtLy5sDCLr6/DTuHOePr8Tl0C6oGImJ+RBRpJIu/kzO+NcNLRcTr\nQOfwUvVwNDAxLU8EjqnTcbuT53yUxj0ZOFyqdBKhqsbXMBFxL/BiN7scDdwYmQeArSRtX5/ocsXX\nUBGxJCLmpOVXgflkI+yUatg5zBlfjzlBNVYAv5M0Ow3VVCRdDS9V8Q9cTttFxJK0/DywXZn9miTN\nkvSApFonsTznY80+EfEmsAIYWuO41jl2Uu779a/p8s9kSTt2sb1RGvnzltcBkh6W9GtJH2xUEOnS\n8V7Ag2ttKsQ57CY+6OE57BXPQRWRpDuBd3ex6d8jYmrOaj4cEYslbQvcIWlB+k+uKPHVTHfxla5E\nREgq9yzETun87QLcJWluRDxZ7Vj7kNuAn0fE3yR9kay1d1iDY+ot5pD9vK2UNAaYAuxa7yAkbQb8\nL3BWRLxS7+Ovz3ri6/E5dILaQBFxRBXqWJy+viDpVrLLNFVJUFWIr6bDS3UXn6SlkraPiCXpEsUL\nZeroPH9PSbqH7L+2WiWoPOejc592SRsBWwLLaxTP2tYbX0SUxnIt2b2+oij0cGalf2wjYrqkqyRt\nExF1G0RW0kCyP/6TIuKXXezS0HO4vvg25Bz6El+DSNpU0uady8BHgS57EDVII4eXmgacnJZPBtZp\n8UkaImnjtLwNcCC1nX4lz/kojftY4K5Id4frYL3xrXU/4iiy+wRFMQ04KfVE2x9YUXKZt+Ekvbvz\nfqKkfcn+dtbrnw/Ssa8D5kfED8vs1rBzmCe+DTqH9erl0Z9ewFiy679/A5YCv03l7wGmp+VdyHpa\nPQw8SnbprTDxpfUxZL1xnqxzfEOBGcATwJ3A1qm8Fbg2LX8ImJvO31zg1DrEtc75AC4EjkrLTcAv\ngEXAH4Fd6vxzt774vpt+1h4G7gZ2r2NsPweWAG+kn71TgdOB09N2kU1K+mT6fpbt/dqg+L5ccu4e\nAD5U5/g+THbP+hGgLb3GFOUc5oyvx+fQQx2ZmVkh+RKfmZkVkhOUmZkVkhOUmZkVkhOUmZkVkhOU\nmZkVkhOUmZkVkhOUmZkVkhOUFYKklY2OwUDS/VWqZxNJv5c0IK2/W9JNkp5MgyNPl/T+tK3LedEk\n7VYyd1CbpFcknSVpkKR703BS1of5QV0rBEkrI2KzRsfRG6XhYxQRbzc6lk6SzgQ2iogJKb77gYkR\n8V9p+yhgi4iYKelpslEPyo7JlhLdYmC/iHhG0vlk04tMqvmHsYZxC8oKRdLXJM1Lr7NKys9TNlvs\nfZJ+LunrG1h/s6QFkm6QtFDSJElHSPqDshl89y3Z90RJf0z/vf+4pDUwJbUCHlWaJiXVO1/ST1L5\n7yRtUiaGdeqVtI+yaTCa0jiNj0raoyTeSan+yZIGp/LHJd1INobjjmXq3VTS7cqmOJgn6bgUQ7ny\nlSVxrvO96MHnPIF3xlA8FHijMzkBRMTDETGzB9+6w4EnI+KZtD4lHcP6snqOJ+WXX+VewEpgb7Ix\nxDYFNiMbt2svYB+ysb2ayGbrfAL4+gYepxl4ExhJ9g/abOB6snHMjgampP1GkE1PMTCtXwWclJY7\nxwbchCw5DC2pd3TadgtwYhfH767e7wCXko2nNr4k3gAOTOvXA19P5W8D+3dXL/CvwE9Kjr9l+lqu\nfGX6Wu57sd7PSTaj7/Ml618FLuvme/JnsqkYZgOnldnneuDLJesDgI5G/9z6VduXr+FakXwYuDUi\n/gog6ZfAQWSJZGpErAZWS7qt8w3K5oL6d7I/sMcqGxn+KuB14J7o+hLQnyNibnr/o8CMiAhJc8n+\nAEP2H/vewJ/SAMyb8M60H1+VNDYt70g2p83zqd62VD67pK5S3dV7Idmo5KvJ/qh3ejYi/pCWf5q2\nTQaeiWzm1O7q/RnwA0mXkE233dlqmVumvFO578W0HJ9zG+DlLj57Od3Oi6ZsdPajgPGdZRHxlqTX\nJW0e2Qyu1gc5QVmvFhFPAadKmpyK/gWYHBG3SboZ6CpB/a1k+e2S9bd553dCZPdMxpe+UdIhwBHA\nARHxmrJ5qJq6qPctsiSxti7rTYaStVYGpjr/2vkx19qvc/2vJWVl65XUQjay9HckzYiICyNiYVfl\nXcTUlfV9zlW8c04ga30dW66yWP+8aEcCcyJi6Vpv3ZgsmVsf5XtQViQzgWPSPZZNyaYFmQn8AfhE\nuj+zGfDxbuoYzjvTXr9VQSwzgGPTf/VI2lrSTmSTEL6UktPuwP5Vqhfgx8B5ZEn1kpL3vFfSAWn5\n08B9eeuV9B7gtYj4KfB9oCVt77K8RLnvxXpFxEvAAEmdSeouYOPO+3Xp+HtKOkj55kX7N7LpMCh5\n/1BgWUS8kScm653cgrLCiIg5km4gm0sJsrmfHgKQNI1srpmlZJenVpSppp0sSbVRwT9gEfGYpP8A\nfifpXWTzBJ0J/AY4XdJ84HEoqNAxAAABHklEQVSyeW0qrlfSP5J1JPhZ6oxxv6TDgKfScc6UdD3Z\npIxXA9vmjHdL4PuS3k5lZ6S3jCxT3llfl98LSc05P+rvyC4T3pkun44FLpf0TbJWz9PAWcB2wK3p\nsuRGwM8i4jedlaSk9RHgi2vVfyhwe85YrJdyN3PrFSRtFhErJQ0mu/xzWvojOhS4iOyP2LXAFcCP\nyP4I3lfmHlSvkRLCryJijwaH0iPp8uHZEfGZGtX/S+DciFhYi/qtGNyCst7iGkkfILu3MTEi5gBE\nxHKyWTtLnVLv4OzvpX8e7pY0ICIqudS6jtRpYoqTU9/nFpSZmRWSO0mYmVkhOUGZmVkhOUGZmVkh\nOUGZmVkhOUGZmVkhOUGZmVkhOUGZmVkh/X8GKa/fLjicxgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist([me[gof_series > 0.05], me[gof_series < 0.05]],bins=20, stacked=True, color=['lightgreen', 'red'], label=['Good fit','Bad fit'])\n",
    "plt.xlabel(r'log$_{10}$ mean expression (C57)')\n",
    "plt.ylabel('Frequency (genes)')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig('figures/ExtDataFig10c.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6604\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(gof_series > 0.05))"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}