Reverse DE.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import pandas as pd\n",
    "\n",
    "rcParams['axes.spines.right'] = False\n",
    "rcParams['axes.spines.top'] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "counts = pd.read_csv('~/Google Drive/Research Projects/Finished/Power-analysis/Public data/GSE60768/GSE60768_umi_counts.csv', index_col=0)\n",
    "\n",
    "counts = counts.astype(int)\n",
    "\n",
    "sample_info = pd.DataFrame({'condition': counts.columns.str.split('_').str.get(1)}, index=counts.columns)\n",
    "\n",
    "sample_info['total_counts'] = counts.sum()\n",
    "sample_info['num_genes'] = (counts > 0).sum()\n",
    "sample_info['ERCC'] = counts[counts.index.str.startswith('ERCC-')].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ESC    48\n",
       "NSC    48\n",
       "2i     48\n",
       "Name: condition, dtype: int64"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_info.condition.value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have two types of cells, you might want to find what molecular features distinguish them from each other. A popular assay for this is RNA-sequencing. If you measure the RNA from different genes in the two cell types, you can identify which RNAs are more abundant in one cell type or the other. This is known as _differential expression_ (DE) analysis, and we usually say that genes are _upregulated_ or _downregulated_ depending if they are more or less abundant. (I'd argue \"enriched\" or \"depleted\" would be better terms, because \"regulated\" suggests some causality you're not measuring.)\n",
    "\n",
    "Abstracting away many details about normalisation and data noise, say $ x^g $ is the gene expression, and $ y $ be and indicator of cell type such that $ y = 1 $ for one, and $ y = -1 $ for the other. In differential expression analysis, for every gene $ g $ we investigate the relation\n",
    "$$\n",
    "x^g = \\beta_0^g + \\beta_1^g \\cdot y + \\varepsilon\n",
    "$$\n",
    "\n",
    "with regards to the data, and ask the question of whether $ \\beta_1^g $ is different from zero in a meaningful way.\n",
    "\n",
    "To make the example more concrete, let's consider the data from Velten et al, where the authors studied mES cells ($ N = 96 $) and NS cells ($ N = 48 $). Say that $ y = 1 $ for mESC, and $ y = -1 $ for NSC. For example, if $ \\beta_1^g $ is postive the gene is more abundant in mESCs, and the magnitude of $ \\beta_1^g $ is the _effect size_.\n",
    "\n",
    "For this simple example, let's investigate 200 genes from the data (selected by having high variance) with expression on a log scale. For the sake of simplicity, let's assume normal distributed noise $ \\varepsilon \\sim \\mathcal{N}(0, \\sigma^2_g) $.\n",
    "\n",
    "The model described above can be implemented in Stan in the following way\n",
    "\n",
    "    data {\n",
    "      int<lower=0> N;\n",
    "      int<lower=0> G;\n",
    "      matrix[N, G] X;\n",
    "      vector[N] y;\n",
    "    }\n",
    "    parameters {\n",
    "      vector[G] beta0;\n",
    "      vector[G] beta;\n",
    "      real<lower=0> sigma[G];\n",
    "    }\n",
    "    model {\n",
    "\n",
    "      beta ~ normal(0, 1.);\n",
    "      beta0 ~ normal(0, 1.);\n",
    "\n",
    "      for (i in 1:G) {\n",
    "        col(X, i) ~ normal(beta0[i] + y * beta[i], sigma[i]);\n",
    "      }\n",
    "\n",
    "    }\n",
    "\n",
    "(To keep it simple, we collect all the genes in a matrix and analyse them all at once).\n",
    "\n",
    "Running the model, we obtain samples from the posterior distribution of the Effect size of each gene ($ \\beta_1^g $). We plot the mean of this, with 95% confidence intervals (CI)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "with open('DE.pkl', 'rb') as fh:\n",
    "    model = pickle.load(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = 200\n",
    "idx = counts.pipe(np.log1p).var(1).sort_values(ascending=False).head(p).index\n",
    "tmp = np.log1p(counts.T[idx])\n",
    "y = (sample_info['condition'] != 'NSC')\n",
    "n = tmp.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = {\n",
    "    'N': n,\n",
    "    'G': p,\n",
    "    'X': tmp.values,\n",
    "    'y': 2 * (y.astype(int).values - 0.5)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "s = model.sampling(data=data, iter=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "beta = s.extract('beta')['beta']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAADgCAYAAAAwjrGWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+cHlV9L/DPdyGGF4mQkATIDyIREPcHSiVitavZ+KsY\ntBBJtywRJe4FIrKtt/SC9PHVinVttEWvBClRgwEvu5qCKLeiqN2scV+tt0TEkE2wUooSAiQEQiAK\nm5Dv/WPnWWdn58yc+fXMeZ79vF+vfWXz7DwzZ86cOfOdM+fMEVUFEREREREVo6nsBBARERERNTIG\n3EREREREBWLATURERERUIAbcREREREQFYsBNRERERFQgBtxERERERAViwE1ERLkSkQ0i8umy00FE\n5AoG3EREjhCRR0XkdyLyvIjsE5F/E5HVItLkW2aDiIyIyAu+n18Y1tchIoe9ZZ4XkV+KyKra7dFY\nGnbWcptERK5hwE1E5Jb3qeorAbwKwBoA1wBYH1jmc6o63ffz+oj17VLV6QCOAfA/AXxFRE4vJOVE\nRBSKATcRkYNU9TlVvRvAnwH4kIi0ZVyfquo9AJ4B8Lrq5yLyWhH5oYg847WAd/r+tkxEtnut44+L\nyF95n18iIkP+9YuIisipgc+mAfgegHm+1vh5InK2iGwRkf0i8pSIfD7LvhERuY4BNxGRw1T1PwDs\nBPDWLOsRkSYR+RMAswE87H02DcAPAfQBOB7AhQBuEpEW72vrAVzutbi3ARhImPYDAN4Dr5Xd+9kF\n4IsAvqiqxwA4BcDGLPtGROQ6BtxERO7bBeA43///yuvjXf25NeK780RkH4DfAbgLwF+q6s+9v70X\nwKOq+jVVPeR9fieAP/X+fhBAi4gco6rPqur9Oe3PQQCnishsVX1BVX+a03qJiJzEgJuIyH3zMdoV\npOofVXWG7+dDEd/dpaozMNqH+wYAb/f97VUA3uQP3gGsBHCi9/cLACwD8GsR+bGIvDmn/ekG8BoA\nD4nIfSLy3pzWS0TkpCPLTgAREZmJyBsxGnAPxS0bRVVfEpFrAPxSRM5X1W8DeAzAj1X1XYbv3Afg\nPBGZAuBKjHb9OAnAAQBH+9J4Ytj3q6sJWe+vAHR5b195P4A7RGSW1wWFiKjhsIWbiMhBInKM1/L7\nDQD/R1UfzLpOVR0BcD2Av/E++hcArxGRi0VkivfzRhFpFpFXiMhKETlWVQ8C2A/gsPe9XwBoFZEz\nReQoAJ+M2OxTAGaJyLG+ffuAiMxR1cMA9nkfHw79NhFRA2DATUTklv8rIs9jtPW5AuDzAILvzr46\n8B7upxOs/xYAC0Xkfar6PIB3Y3Sw5C4ATwL4LICp3rIXA3hURPYDWI3R7iZQ1f8E8CkAPwLwK0S0\nvqvqQwD6ATzidVuZB+AcAMMi8gJGB1BeqKq/S7APRER1RVQnPO0jIiIiIqKcsIWbiIiIiKhADLiJ\niIiIiArEgJuIiIiIqEAMuImIiIiICsSAm4iIiIioQA058c0555yj3//+98tOBhERERE1NrFZqCFb\nuJ9+OskraYmIiIiIitOQATcRERERkSsYcBMRERERFYgBNxERERFRgRhwExEREREVqNSAW0RuEZHd\nIrLN8PcOEXlORB7wfv6m1mkkIiIiIsqi7NcCbgBwI4DbIpb5iaq+tzbJSa+jowMAMDg4WGo6iIiI\niMgtpbZwq+pmAM+UmQYiIiIioiLVQx/ut4jIVhH5noi0lp0YIiIiIqIkyu5SEud+AAtV9QURWQbg\n2wBOC1tQRC4DcBkALFy4sHYpJCIiIiKK4HQLt6ruV9UXvN/vATBFRGYblv2yqi5W1cVz5sypaTqJ\niIiIiEycDrhF5EQREe/3szGa3r3lpoqIiIiIyF6pXUpEpB9AB4DZIrITwN8CmAIAqnozgBUAPiIi\nhwD8DsCFqqolJZeIiIiIKLFSA25V7Yr5+40YfW0gEREREVFdcrpLCRERERFRvWPATURERERUIAbc\nREREREQFYsBNRERERFQgBtw56O/vx/DwMDZv3oy2tjb09/eXnSQiIiIicoTrM006r7+/H5VKBRs3\nbkR7ezuGhobQ3d0NAOjqinwJCxERERFNAtKIr7VevHixbtmypSbbamtrw9q1a7F06dKxzzZt2oSe\nnh5s27atJmkgIiIiolKI1UIMuLM54ogj8OKLL2LKlCljnx08eBBHHXUUXn755ZqkgYiIiIhKYRVw\nsw93Rs3NzRgaGhr32dDQEJqbm0tKERERERG5hAF3RpVKBd3d3di0aRMOHjyITZs2obu7G5VKpeyk\nEREREZEDOGgyo+rAyM7OTuzduxctLS3o7e3lgEkiIiIiAsA+3Lnp6OgAAAwODtZ0u0RERERUGvbh\nJiIiIiIqGwNuIiIiIqICMeAmIiIiIioQA24iIiIiogIx4CYiIiIiKhADbiLKXUdHx9ibe4iIiCY7\nBtxERERERAViwE1EREREVCAG3EREREREBWLATURERERUoFIDbhG5RUR2i8g2w99FRG4QkYdFZKuI\nvKHWaUyDA8aIiIiIqOrIkre/AcCNAG4z/P09AE7zft4E4J+8f50zODgYu0w1CLdZlojIFusWIiK3\nldrCraqbATwTsch5AG7TUT8FMENE5tYmdURERERE2bneh3s+gMd8/9/pfUZEkwC7ZxERUSNwPeC2\nJiKXicgWEdmyZ8+espNTOH8g4lpQ4lp6qLb6+/sxPDyMzZs3o62tDf39/WUniYiIqFSuB9yPAzjJ\n9/8F3mcTqOqXVXWxqi6eM2dOTRJHROP19/ejUqlg48aNeOmll7B27VpUKhUG3URENKm5HnDfDeCD\n3ttK/hDAc6r6RNmJIqJwvb29WL9+PZYuXYopU6Zg6dKlWL9+PXp7e8tOGhERUWlKfUuJiPQD6AAw\nW0R2AvhbAFMAQFVvBnAPgGUAHgbwWwCrykkpEdnYsWMH2tvbx33W3t6OHTt2lJQiIiKi8pUacKtq\nV8zfFcBHa5QcIsqoubkZQ0NDWLp06dhnQ0NDaG5uLjFVRERE5XK9SwkR1ZFKpYLu7m5s2rQJBw8e\nxKZNm9Dd3Y1KpVJ20hoWB6nmhwO+f8+VvHAlHURZlT3xDVnixBZUD7q6Rh9adXZ2Yu/evWhpaUFv\nb+/Y55Qv/yDV9vZ2DA0Nobu7GwCY50REDmHAPUkwYKda6erqwrp16wCwvBXNP0gVwNgg1Z6eHgbc\ndcamjvYvwzqdqL6wSwkRUZ3iIFUKYhcMIjcx4K4R9rMkSobnTLzqIFU/DlIll7k8aRtRkRhw1wAn\nAyFKhueMHQ5SJZo8XLlBcSUd9YZ9uHNWbZXbu3cv2traUKlU2M+SKCGeM3Y4SJWo9th/ntJgC3eO\nTK1yefez9D9qP+mkk7B161ZnHruzGwDlgX2T7XV1daG1tRVve9vbsG3bNgbblJkr9bgr6SDKAwPu\nHJmmtT722GNz62fpD+pvvfVWNDU14c4773TisTu7AVBe2DeZao3B3ShX6vEi08EuEVQGBtw5MrXK\n7du3L7d+lv6g/rOf/Sw2bNgwIcDv7e3Na5dSp82F9FD9Yt9kqiVXgkwXuFKPu5IOytdkvtlhH+4c\nmaa1bmlpQaVSyaWfpT+od+2xu2vpofJk7duYR99k9rMkW2WPGQgb+xO13aTLJ+FKPe5KOojyEtvC\nLSKvEZF/FZFt3v9fJyKfKD5p9SeqVS7Yz3LdunWp7vL8j9pde+zuWnrKNJnv4vPCvsluasSyXWZw\nZ9u67u/yctVVV4Uun0e3mFrV43Fp5fWEGo6qRv4A+DGAswH83PfZtrjvlflz1llnaVn6+vp09uzZ\nKiLa2tqqfX19Y39bsmSJLlmyZMLvNqrL9/X16aJFi3RgYEBvu+02XbhwoQ4MDOjIyIgODAzookWL\nxm0zbNtF8actLj1p1GIf8lJPaXVZlnycTMeglvvaiPna2tqqAwMD4z4bGBjQ1tZWJ7btr1tNyy9Y\nsCCX+rfIetx/fZw7d27kNopMR9YyXMY5EBVb1FM6GrH+UMvY1Cbgvs/71x9wP2C7gTJ+ygy4Vc0F\nKo+AW3V8gV+wYIHOnDkztvDXqpAXWSnU04laT2l1GQNuOwy4sym6sSBKU1OTjoyMjPtsZGREm5qa\nxv7vD7JNy0+fPj23m4Yi6nGbm4ZgWou6ntRbwF1m+cw7HY1Yf2iOAff3AJwC4H7v/ysAfM92A2X8\nNFrAHVbpJFlXnpWWTbqLOqHq6UStp7S6wOacyWudlE0Rx6rWwtJaVguiTfDpD7JNy4tIbOCeRN7H\n0+amISytRZSrLOsso5yU+QQm73TUUz2RgFVsavOWko8CWAfgtSLyOICPAfhIfp1aKIqpf9/u3bsz\nfX8yjr4PY9MftZ76rNZTWolcUtaYAZs38vj7M5uWf9WrXpWqz3Ot6gx/P/ky+2dn6ede1vXUlQGk\nWdMx2V+9GRtwq+ojqvpOAHMAvFZV21X10cJTRgDMr0bas2dPpu/z1Ur5yvOiVVbQ3AjBeiPsA00u\nXV1d6O3tRWdnJ6ZOnYqenp4Jb+TxB9krVqzARRddhOXLl49b/jOf+Uxhr9LM47yyuWko+rWfWQPm\nsq6nrgwgzZIONv5ZvBZQRF4G8A8ArlUd7VMiIver6huKTtxk19HRge3bt4feUe7du9dqHUnuSPka\nNffwmPxennlRdL7yuBWvkfK4q6sL69atAxC+P2GvyZw3bx7OPPPMCcvn8frZIlSD7PXr12PFihUY\nHh7G8uXLsX///pqlNevrH8u6nvrzrr29HUNDQ+ju7s4U6KdJX5Z0lP3qTRfYvId7GKMt4T8QkT9T\n1WcASLHJoqpZs2aFvtt71qxZsd/t6OjAcccdF/r9PO6MbU5Y0zKNdLGsZ2mOA48dxWEZyV8wKA9r\ncY4L3MuU5KahKFm7RJjm2ii6pTmPeQnKTocr3WLKZNOH+5CqXg3gqwB+IiJnAdBik0VVc+bMCX30\nNmfOnEzf54x9lIeiu3AMDg46Fzg0mkbuhuPKvrmSjrIF+8kff/zxNd1+1q4ZZc6A68q8BGnT4Uq3\nmDLZtHALAKjqN0VkGEAfgIWFporGHH/88bj88ssn3FFWWzHSfn+yPMJJohatcvXS8ldP3TeKUq/p\nzkst93+y57WLGvFYZO2a4UpLcz0qoltMvbEJuP9H9RdV3SYibwVwXnFJoqCwx4S2Abfp+1WT8UI3\nGfe50cQdwyKnvp7M5afIfCUKk+d5lkfA7HK3HZfxZiUi4BaRt6vqAIBXicirAn9+IY+Ni8g5AL4I\n4AgAX1XVNYG/dwD4DoD/9j76lqp+Ko9tkz0XW37LCnoma7BVxH7n1b8/uPzu3bvHRsP7W1IAFFq5\nN3rZ8L9loJb5WkuNeuzo9xgwlyePvK/nejaqhXsJgAEA7wv5mwL4VpYNi8gRAL4E4F0AdgK4T0Tu\nVtXtgUV/oqrvzbItV9RjAak3LpyMk7UV0IW8B4A9e/Zg48aNdTUavqi8y2u9HR0dGB4eTpWvWdPg\nSrkid7GMUD0wBtyq+rfev6sK2vbZAB5W1UcAQES+gdGuKsGAm+rMZK70imwF9F9UeIEx27t3b8OM\nhnfpOEfla/Am8/DhwzUfEGfL1RvipMfYtHyS9biUFy6VdaIixL6lRET+QkSOkVFfFZH7ReTdOWx7\nPoDHfP/f6X0W9BYR2Soi3xOR1oh0XiYiW0Rki+2kMGXJY7Yl/9sbJsObHIJ5ZjvTZlXStwSkfauA\nSxMNTdY3I1Rfpek32UbDh8laHkz5Om/evAkTWuzbty/xOZqF7b65MPlGluOQZ11vmxdFzg7YaNeu\nenlrU9Zj2mjHrVZsXgv4YVXdD+DdAGYBuBjAmuiv5OZ+AAtV9XUA1gL4tmlBVf2yqi5W1cW2r8wr\ngwsVvuuClVZYntX6gm6rHt416lognvVmKqiIV2G6lmdlCMvXlStX4sCBAxNuMm+//fbQ2XBt87Go\n/C7yhrjMMpJm2zZ5kfR6ZZOOyRCs1WIK87TlrYwYhPXnKOvXAgJYBuA2VR0WkTwmvnkcwEm+/y/w\nPhvjBfrV3+8RkZtEZLaqPp3D9gsTVZlwtqXkwvLs9ttvR2dnZ6b12jwGT/qoPOvECFkezZf13SzC\nBjiuXLky0zpdfhWmiwOQbYXl64wZM/DQQw9lmg03jbTlNe0NcSN2d7DJC9P1qrOzE+vWrXO+e1tR\naYtaX9yg7bJfEFB2DFJUOXG5HFbZtHD/TER+gNGA+14ReSWAwzls+z4Ap4nIIhF5BYALAdztX0BE\nTqwG9yJytpfe4mrxGqiHFlCTvFsibZnyLMsF3V8pmlrNbZYJCpsYYeXKlTh8eOIpE9eSf/7552PX\nrl1W+Z3lKUCa/czLnj17rFtH4/jLZ29vL+bMmRM6OUPRrU+1aN0qS9jEJaauJjaz4aaRpay7NvlG\nmS1/NnlRRN3rgiJb2cPqtLK6FYap5xik3tkE3N0APg7gjar6WwCvAJB5IKWqHgJwJYB7AewAsNFr\nPV8tIqu9xVYA2CYivwBwA4ALVbWuZ7l0rcK3VWZQZsqzLBd0m0DPNhgMBnrnnnsuOjs7MXXqVPT0\n9GDGjBlWrW/+loc77rgDfX19uOuuu6zyO+zxsD+tURf2PIPepEwD8YIX9LibPdvymffj1GC+hqWj\n0buMBbuafOITn8Dy5cvHWp/zriPiynpV2I1PrWYK9JeLshoq4tjkRZa6t+ibCVe7Kbg+aLteY5BG\nEBtwq+phVb1fVfd5/9+rqlvz2Liq3qOqr1HVU1S11/vsZlW92fv9RlVtVdXXq+ofquq/5bHdMtlW\n+LWupONa5dIEZXm19JlajbP01bcJ9GyWCQvgvvvd745rXbUJtjs6OrB9+/ax7dkGFVVZWqJsg97q\n/uZZLm1aR21aNG3LZ9GDWk2tW3/+53+ey/rzFjyeLS0t6OjoSHScjz/+ePT29o7dZN5yyy3WN4ph\naYirJ2zKeliZueqqq/B3f/d349La09NTaHcjl8efdHV1xeaF6Xrl8jipsrk+aLvM6emzqvenhzYt\n3JQjm0ouj0o6yQXTplUuSVBm2ofqOpOeNGF5ZttqbGIT6NksYwrg0rQO+7eXNIBO0xJVPQ7Tpk2z\n+m6wnCTp8mJiGojnv6Db3HzYls+iH6cmPU/KZKpnHn744cT1T7WryaxZs3D77bdb3yia6omobdmU\n9agyE+wWU2S/1aQ3zlW1ar2NywvT9crVVz7mLc1xKGLQdl46Ojqwbt26Ccf0qKOOSjR7dXCdeZTV\nuPU0wgsnGHCXIK6SM1XSTU12hyusYL744ou4/PLLQ5e36XMWDD77+/vH7tjDgmZTIHr11VcnvsCG\n5VnWCt8m0LNZJs8+jv7tJQ2gkz4F8AfPN910Ey655BJjl4DqsfWXk6RdXkyCraNhN1PBPO7v70dP\nT8+4Lgu2/YhtH6eabljjbmRr3Z85C1M988ILL6TuYpT0hiPNDatNWXel77Er6ciiljcoaZie0pQl\nrE5L+hSl6BuuvK+nVf50570PLr1yNy2b93B/3eYzyk/WVrikBdOmz5k/GPz617+Oj3/84/jKV75i\nvNM07cO+fftyaxHOwibQs1kmSWAcF6z5t7d9+3asXLnSOoBO+hTAHzxffPHFWLNmDS699NLQLgEf\n+9jHsHXr1nHlJEs/2rC0Vyv/SqWCPXv2jFven8fVm8m1a9eOC/SnTp1qdcNh8zg1S8uvzU1aWnlf\nwEzn6P79+1MHiXE35sEynyYgtSnrRYz7CBN3TGqVjnqXtqtAmWOLqsLKgOs3KfWoEQZ72jSZjpts\nxpuS/axikkNA9kENSQumTZ8zfzB4xRVXYMOGDZEBvWkfDhw4MCFtO3fuxIsvvliT/uqmN1kEA71q\nGuKCQds+jrYXBv/2rr/++kTdaJK0WgRvsrq6urBjxw5MmzZtXJeAJ598EkcddRTuvPNOtLS0JOry\nkqRb0e7du41dm5YuXTqWx8FA/8knn8TRRx+NXbt24eDBg7jgggsi88umS1ealt/q/jz00ENW6fAr\nql9i3A2e6Rw95phjEgXNfnE35sEyHxeQmvImrqwXMe4jjbLT4erAQr8sXQXKHPAdVMtxV2X0Zc7y\n3u880toQgz1VNfQHwLUAngdwCMB+7+d5jL6W7+9N33Ph56yzzlLXLVmyRJcsWRL6t76+Pl20aJEO\nDAzoyMiIDgwM6KJFi7Svr89q3a2trTowMDDus4GBAW1tbZ2wndmzZ6uI6Ny5c0O351+mtbVVm5ub\nVUR0ZGRk3LpGRka0qakpdh9OPvnkcWnr6+vThQsXjltu7ty52tzcbMyzYJr8+WLK1yVLlmhzc/OE\nNM2dO1fnz58f+nk1DabvmvKor69vXDqWLFmis2fPDj0ms2fPnpDupL+H5VHc76b0BI+tvyz5j2mw\njPX19ekpp5wyLg9M5XDBggWh+T1z5kxjufXncTV9YWUseNxM51jU35qamkLLt6nci4h1OoLbjTrX\nw849U5qD+2Yq6/7zypR//vPhtttuizw/w/atmu7p06dHlvmo/W9ubo485+LytbpuU/5FHX9TnqY9\n3/zpWLBggc6cOTP2mCZJn219mLQuiUpPXutVtb9ehaUp6pxMkoaobdjkjc35lnSdaeqMJPtjc6yy\nLGPKF1Na4/IlSV1pGyvlyCo2jV/A8eA67KfeA27V6Eo0js0JGVymUqnoscceO257pgtyVGAUtw/B\ndZ5yyimxF2V/niUNfP3fNQWZxx57bGwwbPpudZ/jLkhJLgxhv0ddtNNcYE0XiOCxDQaf1UDaf5MW\nFpQtWrTIGLiaArG4G7ngcTBdqMNuYoKi/mZar6mc2JQf00Ur6U1J2AU8bN+C5TXshqj6eVhQahs0\nR5WxuDLvT5s/DVdeeWXktoPnXNaA00ZcsBF3cxR1zqUNylTj63rb+iBuW/5lTHWR7c1EUNQNblya\nkjRkpGGbN3HpqLK9rkeV7bQ3KMH9sSkDWZaJypewtNrktU1MkbSBMie5BdzLARzr+/8MAOfbbqCM\nn0YIuG2XMYk7seOClahlXvnKV2a6a/WnzfaiXF2X6QQOC1CCF56obZmCSttW/eB+2lY8NoGL/0Jt\nCm6DF1ib1r2wZYJBQdQNUVxQZroxM+WlTYDlT1/chTrq/In6m03Lr//ztDdTquZgw6Z12CRY1rME\nZVHlPq6M2Qb9/u3602rKG/85l7YVPEndGhfQ2wTSNnWAn2364oIvm/23yZew4+Ovi5I8CbHdh7iy\nHpf3ebR62hyHuGtLVZLAMKpsx50XtvuT5LinWSau/ohaj+32VNM/IclZbgH3AyGf/dx2A2X8MOCO\n/75Nq0LUMknv1E2fJ6lso07gqGAtLjD0t1CmbdWPqixsLspRlZZNq64pKE26repF1d+CFXchTdLV\nIqxbUXUf4m7kgulL0voaFHduJWnFS3szpWo+nkluRMPS7c+bLEFZkhvcqG5YtjeK/rTapDvJk6eo\n4xCVn3Fl0iaQDjtPgjf2Nl3jgmxuSuL23yZf4o5PlqdNcV3DTN+J6qaTpCtDFJvjYFsGkl7rTOu0\necIcdxyjfrdZ3nadWVq4bbanGn8O1EhuAffWkM8etN1AGT8MuOO/n6WFO8njurgTPEllG3UC2wR9\nplYYf8tl2lb9qMrCv682fUqDv/v3zeYCG1X5p0lD0pbM6vaCfa/juirFtUoF05fkEX2Q7cXD5ubF\n9gYnyeNQ001JVKufqfUxS1BmClzi0hfcZ9vuIf602gS7SZ48mX6PYhP02wTSwfMkyVOHrOlLsv9x\n2/UfH9PvYXkQ170kyXgFm7EbSQK9KDbHwbYuStJ1Jqpsmxoy4m7YbOtW0/4nXSaq/gje+ISlzbaM\nNloL9y0APg/gFO/n8wA22G6gjB+XAu60FZvtMmm/bxPoZhmUFpYG0+9J+icnDQCCrQHBx9rBPqtR\nlVySCstmv21/T9rCXcTgviyVqk3lbzOgLSwNSS7UpnVFbcO29Tpuf5I++k7T6mcayJq2L3TYuVFN\nX5LuPFFlMhgY2wzG9R+rPFq4TWVBNfkNrimQDra6JgmUo+TRh9svbrs2Ldxpym6S9Nmck0m6Mpjk\n0d/aL68WblNDRlia/PtSRF9/mzoj6poelTbbAcWN1od7GoA1ALYAuA/AZwBMs91AGT8MuO2+bxOs\npAmMTGmIOnmXJAjckrQUJg0M0g7CMX1uG8REtTQnfTSfZXBo1PGJOo55tp5kWd72nLHdRpIbsKh+\n8mm6vyS9mYgq63HnVZJuSKrJW1ZtA+OkT7yy1Bn+9ZjyNWkXrqigKpgO277xcYFDkqdCSc4x07bi\n+nBn6ZNtkz6bczLuJjNO0v7WNvWPbdn2H0/T28PitpXkfE3z1DNsf/wvXgimO81T6+o5E7dvwTxL\n218/o3wC7rEFHQ+y/T/1FnAnrRSzbtsmHUnTmnT9we9GnXhxrVVVYSddkm4xS1IEHmH76U+PzeC7\n4Bti4gY4Jn0K4F9/kpuP4PFJW17zKld5njO22zCVybAuRlFdTdIO8Eyyn1FlPS4oS3qjmbTFLMlT\nkKQ3GnH7luYJTNL9tHlCZlOu4gZ/mxRxjpmY6iL/72nGH9imz+acjGuYiGPbTSFrWQ0uHxXEJrlR\nDP4tyc247VNP01ONvMZlBc+ZPMpugaxi0yPj3tMtIm8B8FUA0wEsFJHXA7hcVa9I+epvolA2M15G\n6erqwrp16wAAg4ODY593d3dj/fr1aG9vx9DQUOTEE8cffzwuv/xydHZ2Yu/evWhpaYmduMTEP4EK\ngLFJGZYvX4677roLS5cuRX9//9gU6f70xe1bdQIC/36a9uHEE0/EXXfdhZ6eHgwNDY2lB5g8s95F\n5VeU6kQuwfJz5JFHhh7bzs7Ose/u2bMHGzduxNKlS8cmbSgy76uTrISVddO5UZV0SvbqREH+8yRq\n+uqw88q0vKmsm8Ttm/84ABibqKunpyd2BsCw/TzqqKPGtjc8PIy9e/fi6aefHqtTbI5zWLnq7u4O\nLVe2aa2VqLqoo6MDw8PDhZZ123PyyCOPxKWXXopHHnkktnwGmSb22r59+9j//RP2RNXdfnFlO3jN\n+PSnP40uQY+DAAAb6klEQVR3vOMd6OzsxLZt26zSHsZU/8yYMSO2Hqvua7Wst7W1oVKpjMsjf7qT\nzh5bnXhvMlyXbGaa/AKAP8bohDdQ1V8AeFuRiaJiuTr7mGnGy+OOOw7AaAWVNGBKOuV59TutljM2\nRrGZOtt2ivSkqvswa9assZkjy571LkqaY1sL/hlW/eVn3759sRcVfxCbNe9tZrFLU9arTOde1EUv\n6fTVweXXrVtXk3oojxv5YH1gmjl26tSpVjPPhpWr3t5e/OY3v8lt+upaznzo559ttIh6xvacrM6e\nCyDx9OpxM6ACxdTdSYNVW6YZkZ977jnj9vzl56qrrpowE+jChQtDZx62yTu/sPISds5UuXqtsGET\ncENVHwt89HIBaaFJrqiKOq8AOijuxsVm6mxTBfvMM8/kkkZ/sFENyHp6ekIDsnquyKps9yHJvoaV\nn7gAdXBwEC0tLWPLVPP+0ksvTRwMm4I7U9CdpqwXHSSVJXgcqmymhI46v01Tir/00kuhgfT27dsn\nlLewG5a8pq8OTpd+/vnnY9euXdi8eTNOOukkbN26tbBpwU0BcV71LpDunEzC5ga5iOA4abBqK+xm\nvLe317i96dOnj5WflpaWsUaballfv349Dh06NJZH/vUkbVwIKy/nnnsu9uzZU/ObxaLZBNyPed1K\nVESmiMhfAUh+u03jNEJwk7e0FbXpwmjTkm9zHNIeK1PFM3369NCKqsp/gU267eDywYtQV1cX1q5d\ni1mzZuV681ELLp0zNgFq8PifeOKJ+O1vf4vXvva1ifLeFNxlfQpSNTg4iO3bt084966//vpxj9CT\nrtOVY2Vq3atUKqnXGdUFJ2nLfxFp9be+3nHHHWPd1m699VY0NTXhzjvvHNdamXfQ7c+DSqUyLniK\n21awZd4mbYODg7jhhhtyu2m0eVqUJDi2faqc5UlY3DbCyqVpe/7uOaYbi127do3l0fbt27Fy5Ups\n2rQJK1aswEUXXYTly5eneqpcqVTw3e9+16qBwdWn9SaxfbgBrAbwRQDzATwO4AcAPlpkoshdRV9E\nk/bfdFlY/89qxePvz7py5Urcfvvt4/py9vb25pIGU3/Hem+5LJtNX/+o459E0v7VacX1hS5SkdtL\n2t/cRlH9TvNKq6l/bVtbGzZs2BDaRxzAhH66WfuN+5/O+Os3/75mWd4vr/PNv76oa1HUmAkbYf2i\n89iHqPWG7WPY9h566KEJXUSCZb25uXlcHgXrw3nz5uHMM89MfB03jX0K9iuvR8YWbhH5rPfrUlVd\nqaonqOrxqvoBVc23pieqc6YWPdMjfv/n119//YRHfUkudFF3+UU83vXvq0stmbVm030jj+5MeT4q\nn6yytDqHKbILTh5p9be++oPvqMGAwW5LebR8hz2dWb9+fWiDwuDgIJqamqyXD1NU90HTttLWrWHd\nxKr5nWUfotYbtR9R3XNsn7pEpTvJdaKofuwuiGrhXiYiHwdwLYB/rlF6KKO0b2NI+50y1tloim5Z\njGqpMW3P9eNWi/S5kgdFPKXIUk+4rhb7lufbjIrgb331t1AmeVtFHm9HSTpg1RRspRk0Wgtp6+4s\nb86pyvuNPMD41vHqk9cVK1ZgeHgYy5cvx/79+3N5QhQl7Rud6qEuiwq4vw/gWQDTRWQ/AAGg1X9V\n9ZgapG9Sy1qAXC6ALqeNKEqty67rwV2cRq3H0nZ/q8X+BLsKVIOna665Bpdccgk2bNgwrsuG6W0V\nWQNdU9cb0yDQqO4Lwa4Shw8frptzIMh0I5J2zEQe6w1257nuuuvGBdn+LiJ5MB3PYFed6667Djfe\neCP279+fW1enskQF3J9Q1f8lIt9R1fNqliKqG65eCF1Wr3mWJd1Jv1v08iZJ+j4WmY4wZQd39Vpu\n42Q95q4wdWcL6187f/58XHDBBdi3b99Ya2Vvb2+iwNi07WA6TO8bN3URCesX3d3djXPPPTfxO6+T\nqOUTn8HBQbS1tVm34iZJU5axBcHW8eA7wPMcTxX1DvPgzWJ1HomkffpdFBVw/zuANwDYX9TGReQc\njA7IPALAV1V1TeDv4v19GYDfArhEVe8vKj1UnEa5sFF+XAriwi4ARVbstueDS3lUS7W46ar1MU8r\nj2DQZuKsJIGxbVqSTHhUTScwcdBo2oF0/nRmyce8r19ZB1yaZOl+VqvB2UD8wMhqeR0eHh57JWF1\nOdcmgkoiKuB+hYhcBOAtIvL+4B9V9VtZNiwiRwD4EoB3AdgJ4D4RuVtV/c8+3gPgNO/nTQD+yfuX\n6ojNha2WgUWWfsuTNQCyUc95E3YBKKpiz/Imhizq+fgUIckxr+ZdVLeGes7fIt7k4l93kn7OYct/\n4AMfKHS+gihFnK95v1GlKkv3s1rO+Gg7MDLrpFWuiXoP92oAbwUwA8D7Aj/vzWHbZwN4WFUfUdUR\nAN8AEOy6ch6A27y56n8KYIaIzM1h2+SxGT2cZIRxmLAZuZKMPE+SpqxpdVkj71uRbN7rW8vBWkne\n3BAmzXuKi5K0TLpUhpMe8+BkMlHvB7blUn7k/SaXPOU1IVAaWc9Xk6LeqJL2bSFp37yTpgybjmdL\nS8u4z0xvaPL36XehHrRlbOFW1SEAQyKyRVXXF7Dt+QD8M1juxMTW67Bl5gN4ooD05I7dKEbV28jz\nesDWeDu23QaiBmsB+eZlllabslrH8+ZC3ZhkgF6lUknUrYHnXr5Mfbvzmq/AZHBwEEcccUSm61ct\nynoeA0rzGJwd1Z/fz/Z4msYAhPXpr4d60Bhwi8jVqvo5VV0vIn+qqv/s+9tnVPWva5NEOyJyGYDL\nAGDBggXYtm3bhGVOP/10TJkyBU899VToLG15/v3WW2/FTTfdFFog3v72t499/5577sGDDz6IZ599\nFieddBIOHDiAffv24fTTT0d3dzeWLVuWaPsAMDIyErn/IyMjOHjw4IRlisqfU089NfTCduqpp45L\nQ1j6Dhw4AADYtm2bcf3VZQ4ePDju7/7vHjhwAEcffTQAWKX/xhtvHPuuaf/864/a/yLSZ7N+P//6\na338y/z7Jz/5ydBuA6tXr8YZZ5wx9t0PfvCDEyr2VatW4Yorrog8vknLJwAcd9xxkcG9af9U1fja\nr9WrV2PFihWYMmUKvvnNb2LPnj2RxzdYdqt/B+Lrj6zH5+abb8aaNWvG1Y2rVq3CoUOHcPHFF2de\nv6l8qypEZOz7H/zgB7Fq1Sp87WtfG1dHL126FNdee+2E9D322GPGx+BRxz9N/gXrBiD6/Dat378e\n//7X4vj769Ck9VswfWeccQauu+66sWDwtNNOGzuHg+enqfwHj79fVPrirl9Rx+eee+6ZEAd8+MMf\nxmOPPYZly5YZj4/N9acqbBsrV67ESy+9ZF1+brzxxrH1r127FieffDK+9rWvYdWqVThw4MCE9GUt\n32ecccZYnfXss8/itNNOw3XXXYeurq5x18dp06bhkksuGVuupaUFV199Nb7whS8Yu4P546u06Uv6\n97a2tgnLh1LV0B8A94f9Hvb/ND8A3gzgXt//rwVwbWCZdQC6fP//JYC5ces+66yztGytra06MDAw\n7rOBgQFtbW0d+39fX58uWrRIBwYG9LbbbtOFCxfqwMCAjoyM6MDAgC5atEj7+voSbXfJkiW6ZMkS\n49/7+vp09uzZKiLa2tqaeP1p+PfTZt+q+2CbVtM+x+VFFkny0SZ9WdJq890i88JlTU1NOjIyMu6z\nkZERbWpqmrBsmnMjTd4nPR/86xER6/3Jmu6i2NSNWSSpD8KOuSl9M2fODP189uzZuaTbJK+6waa+\nKapcJF1vEemL+m7U3+LO16jvxpX1PK4Bpm1Uy2WWvC+6nrA9JsHlktTrNWIV90YNmhTD72H/T+M+\nAKeJyCKMThl/IYCLAsvcDeBKEfkGRrubPKeqddGdxKYbhe2Uu3k9IilrVH6aATkuPzqvl7cbUHxX\nEWD82wtqMbV5lgFqSd9r7KKyupjFvTovboDevn37CnmzBLkt7nyNqitqUdYbeWZGE5t63UVRgybV\n8HvY/xNT1UMArgRwL4AdADaq6rCIrBaR1d5i9wB4BMDDAL4C4Iqs260Vm0EeNlPu5nliJhm8GDVd\neBpJBuQM5jDFb5EDKvIaBOrSgKlaybtcxbGdljittIOO0w5QCxvYlHR/yh5sVOYAOBtRA7rSTuVN\n9S3t+Wpb1rOck6ZtzJo1K/O5XuY1KmrbRdfrhTE1fQN4GaPv4H4ewCHv9+r/D9o2oZfx40KXEpvH\nxv5HQXk9Zo16RJPkMUwRj5KSrNM2rWGPhNM+srfl0uOsMroIZVFGV4a4PKrlI1TTdpMsn+WYF31u\nlJ2GPLoFZelCUIS8tscuJeHfzdp1MYpNWWpubs50PoRtY+7cuTp//vzUXddcqwPDlnPs2mfXldp2\nwXr6cSHgVo0vELXuw50kqC874E7aB96fZyeffHKhfUSL7oNqy4XgKakyAu647dZLwJ3Hd10qu3lf\nLLP0jU9yMa/XgLvMbRVx05CE6XgmKTN5b7u6ztmzZ2c+J5OMRYhbr4vlrdY3iCkw4HZBXIHwnygL\nFizQmTNnZroIRW2vFpVL2rSlSaupQslrYFmWtNWCK8FTEgy47dJU1HddejqTZh/yalCwWacLF/la\nt+I5FMCESpO+qPq6Vo1QUd/N65wMbiPtel2sA104F2NYxaZRgyapBmym3LUV977PImcTy5tNWk39\n3qdNm1bogApX8pHvN6ek6nWwkY1GOx/KGJzdiGNKomYUdaHMFHVOpl2vi2XAxTSlETVokupI2Cxo\nlUplwiAJl2cTC4pLq2mwyIwZMwofUOFCPro++IzcU7eDjSw02vlQxAy9k1FUUO1CmSnqnGzkc71e\nsYW7QUTdxbscVGdhmq3qc5/7HIDyW6CLVtbsa1S/XHk6U4RGOx9caH1tBFEtvS6UmaLOyXo41xul\n5doWA+4GUXblnKUrTFpxFUot3qlcpnqoUMk9Ye+ebgSNdj40cveftNKU16ig2pUyU9Q52ajner1i\nwN0gJmvlPNkrlMm+/0R+ac4HV88bF1pfG0FcUG1bZlwtJ1Q/GHA3CFbORESNw5XW10bAhonGVG/H\nkgF3g2DlTJRdvVXgVHu17jbHQJFssHy4jwF3A2HlTJRM3Ks0i8RzlGhy4LlOAANuIqqRMoNbU3pq\n/Z7jvDTCBdz1fXA9fURUX/gebpqgGhht3rwZbW1tE97lndbg4CAvYpOU7XviixJW9lx9z3FHR8fY\nW3+IiKgxMOCmccoOjKgxuRjclv0qTUqmqIYAokbGhi53MOCmcVwMjKj+uRjcujDLHNlhQwA1OgbG\njY8BN43jYmBE9c/F4JZTH9cPNgQQUb3joMmC2dyx+pcp+w7XtQl0ys6PKC6nzTUuvid+Mr9Ks97K\nLhsCKCtOakNlY8A9CUVVLi4GRlT/XA1u+SrN+pCkIYDHkYhcxICbxnE1MKL6x+CW0prMDQE8V4ga\nAwNumiBpYFT2+5VNaZwsF6rJsp80ebEhgIjqHQNuyqSeJw8hovrBJyREVM/4lhLKhG8PICIiIorG\ngJsy4dsDiIiIiKKVEnCLyHEi8kMR+ZX370zDco+KyIMi8oCIbKl1Oimei+9XJqpXnE2RiKgxldXC\n/XEA/6qqpwH4V+//JktV9UxVXVybpNW3Ws9WxclDiPLB2RSJiBpXWYMmzwPQ4f1+K4BBANeUlBbK\ngG8PIMqHfzwEgLHxED09PTyfiIjqXFkB9wmq+oT3+5MATjAspwB+JCIvA1inql+uSeooEb49gCg7\njocgImpchQXcIvIjACeG/GlcXwNVVRFRw2raVfVxETkewA9F5CFV3WzY3mUALgOAhQsXZkg5EU0m\nrtwkJplNkYiI6kthfbhV9Z2q2hby8x0AT4nIXADw/t1tWMfj3r+7AdwF4OyI7X1ZVRer6uI5c+bk\nv0NERAXieAgiosZVVpeSuwF8CMAa79/vBBcQkWkAmlT1ee/3dwP4VE1TSUS5cqU12UUcDxGP5YeI\n6lVZAfcaABtFpBvArwF0AoCIzAPwVVVdhtF+3XeJSDWdfar6/ZLSS0RUOI6HICpXR0cHAJ5/lL9S\nAm5V3QvgHSGf7wKwzPv9EQCvr3HSiIiIiIhyVVYLNzmOd/dERERE+eDU7kREREREBWLATURERERU\nIAbcREREREQFYsBNREREk15/fz+Gh4exefNmtLW1ob+/v+wkUQPhoEkiIiKa1Pr7+1GpVLBx40a0\nt7djaGgI3d3dAMB34VMuRNU0q3r9Wrx4sW7ZsqXsZBARJcb3ABPVXltbG9auXYulS5eOfbZp0yb0\n9PRg27ZtJaaM6oDYLMQuJURERDSp7dixA+3t7eM+a29vx44dO0pKETUaBtxERA4ZHBxk6zZRjTU3\nN2NoaGjcZ0NDQ2hubi4pRdRoGHATERHRpFapVNDd3Y1Nmzbh4MGD2LRpE7q7u1GpVMpOGjUIDpok\nIiKiSa06MLKzsxN79+5FS0sLent7OWCScsNBk0RERETgoGVKhYMmiYiIiIjKxoCbiIiIiKhADLiJ\niIiIiArEgJuIiIiIqEAMuImIiIiICsTXAhIRERGBbyeh4rCFm4iIiIioQAy4iYiIiIgKxICbiIiI\niKhADLiJiIiIiApUSsAtIn8qIsMiclhEFkcsd46I/FJEHhaRj9cyjUREREREeSirhXsbgPcD2Gxa\nQESOAPAlAO8B0AKgS0RaapM8IiIiIqJ8lPJaQFXdAQAiErXY2QAeVtVHvGW/AeA8ANsLTyARERER\nUU5cfg/3fACP+f6/E8CbTAuLyGUALvP++4KI/LLAtJnMBvB0CdutV8yv5JhnyTC/kmOeJcP8So55\nlgzzK7la5tn3VfWcuIUKC7hF5EcATgz5U0VVv5P39lT1ywC+nPd6kxCRLapq7JNO4zG/kmOeJcP8\nSo55lgzzKznmWTLMr+RczLPCAm5VfWfGVTwO4CTf/xd4nxERERER1Q2XXwt4H4DTRGSRiLwCwIUA\n7i45TUREREREiZT1WsDlIrITwJsBfFdE7vU+nyci9wCAqh4CcCWAewHsALBRVYfLSG8CpXZpqUPM\nr+SYZ8kwv5JjniXD/EqOeZYM8ys55/JMVLXsNBARERERNSyXu5QQEREREdU9Btw54IyY8UTkJBHZ\nJCLbvVlG/8L7/JMi8riIPOD9LCs7ra4QkUdF5EEvX7Z4nx0nIj8UkV95/84sO52uEJHTfeXoARHZ\nLyIfYxn7PRG5RUR2i8g232fGMiUi13r12i9F5I/LSXW5DHn2DyLykIhsFZG7RGSG9/nJIvI7X1m7\nubyUl8OQX8ZzkGXMmGff9OXXoyLygPc5y5g5nnC6LmOXkoy8GTH/E8C7MPqu8PsAdKkqJ+jxEZG5\nAOaq6v0i8koAPwNwPoBOAC+o6j+WmkAHicijABar6tO+zz4H4BlVXePd3M1U1WvKSqOrvPPycYy+\nu38VWMYAACLyNgAvALhNVdu8z0LLlIzO7NuP0UnI5gH4EYDXqOrLJSW/FIY8ezeAAVU9JCKfBQAv\nz04G8C/V5SYjQ359EiHnIMvYqLA8C/z9egDPqeqnWMYi44lL4HBdxhbu7MZmxFTVEQDVGTHJR1Wf\nUNX7vd+fx+hA2PnlpqounQfgVu/3WzFaydBE7wDwX6r667IT4hJV3QzgmcDHpjJ1HoBvqOpLqvrf\nAB7GaH03qYTlmar+wBvYDwA/xehrawnGMmbCMoboPBMRwWjDVH9NE+WwiHjC6bqMAXd2YTNiMpCM\n4N2h/wGA/+d91OM9mr2FXSTGUQA/EpGfyehMqgBwgqo+4f3+JIATykma8y7E+AsUy5iZqUyxbrPz\nYQDf8/1/kfeo/8ci8tayEuWgsHOQZSzeWwE8paq/8n3GMuYJxBNO12UMuKmmRGQ6gDsBfExV9wP4\nJwCvBnAmgCcAXF9i8lzTrqpnAngPgI96jx3H6Gh/MPYJC5DR9/b/CYB/9j5iGbPEMpWMiFQAHAJw\nu/fREwAWeuftXwLoE5FjykqfQ3gOpteF8Y0HLGOekHhijIt1GQPu7DgjpiURmYLRk+N2Vf0WAKjq\nU6r6sqoeBvAVTMLHiSaq+rj3724Ad2E0b57y+q9V+7HtLi+FznoPgPtV9SmAZcyCqUyxbosgIpcA\neC+Ald7FHd4j673e7z8D8F8AXlNaIh0RcQ6yjEUQkSMBvB/AN6ufsYyNCosn4HhdxoA7O86IacHr\nh7YewA5V/bzv87m+xZYD2Bb87mQkItO8wSAQkWkA3o3RvLkbwIe8xT4E4DvlpNBp41qEWMZimcrU\n3QAuFJGpIrIIwGkA/qOE9DlHRM4BcDWAP1HV3/o+n+MN2IWIvBqjefZIOal0R8Q5yDIW7Z0AHlLV\nndUPWMbM8QQcr8uOrPUGG403Sr06I+YRAG6pgxkxy/BHAC4G8GD19UYA/hpAl4icidFHP48CuLyc\n5DnnBAB3jdYrOBJAn6p+X0TuA7BRRLoB/Bqjg2nI492cvAvjy9HnWMZGiUg/gA4As2V0tt+/BbAG\nIWVKVYdFZCOA7RjtNvHRyfb2CMCYZ9cCmArgh945+lNVXQ3gbQA+JSIHARwGsFpVbQcQNgRDfnWE\nnYMsY6PC8kxV12PiWBSAZQwwxxNO12V8LSARERERUYHYpYSIiIiIqEAMuImIiIiICsSAm4iIiIio\nQAy4iYiIiIgKxICbiIiIiKhADLiJiBwnIieISJ+IPCIiPxORfxeR5WWni4iI7DDgJiJymDfJw7cB\nbFbVV6vqWRh9P++CclNGRES2GHATEbnt7QBGVPXm6geq+mtVXQsAInKEiPyDiNwnIltF5HLv8w4R\nGRSRO0TkIRG53QveISJniciPvdbyewMzAcJbZoOI3CAi/+a1rK/wrfdffMvd6E1zDhF5VET+XkQe\nEJEtIvIGb/3/JSKri8wkIiKXMeAmInJbK4D7I/7eDeA5VX0jgDcCuNSbvhgA/gDAxwC0AHg1gD8S\nkSkA1gJY4bWW3wKg17DuuQDaAbwXo7O42fiNqp4J4CcANgBYAeAPAVxn+X0ioobDqd2JiOqIiHwJ\no0HwiBdkvxvA66ot0ACOBXAagBEA/6GqO73vPQDgZAD7ALTh99OSHwHgCcPmvq2qhwFsF5ETLJN4\nt/fvgwCmq+rzAJ4XkZdEZIaq7rPfWyKixsCAm4jIbcMALqj+R1U/KiKzAWzxPhIAPap6r/9LItIB\n4CXfRy9jtM4XAMOq+maLbfu/L96/hzD+6ehRhu8cDnz/MHjNIaJJil1KiIjcNgDgKBH5iO+zo32/\n3wvgI15XEYjIa0RkWsT6fglgjoi82Vt+ioi0JkjPrwG0iMhUEZkB4B0JvktENCmxtYGIyGGqqiJy\nPoAviMjVAPYAOADgGm+Rr2K0q8j93qDIPQDOj1jfiNf95AYRORaj14H/jdGWdJv0PCYiGwFsA/Df\nAH6easeIiCYRUdWy00BERERE1LDYpYSIiIiIqEAMuImIiIiICsSAm4iIiIioQAy4iYiIiIgKxICb\niIiIiKhADLiJiIiIiArEgJuIiIiIqEAMuImIiIiICvT/AffwdIo+XMvHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x117ec3518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(12, 3)\n",
    "plt.errorbar(np.arange(beta.shape[1]), beta.mean(0), yerr=2*beta.std(0), fmt='o', color='k', markerfacecolor='w');\n",
    "plt.axhline(0, lw=1, ls='--', color='lightgrey', zorder=-1)\n",
    "plt.xlabel('Gene num')\n",
    "plt.ylabel('Effect size');\n",
    "plt.title('DE Results');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several of the 200 genes have effect sizes such that the CI is far away from 0. A handy way to quantify the uncertainty of the effect sizes is to invsetigate the probability of the effect size being 0, let's call this a _P-value_. A simple way to do that in this setting is\n",
    "$$\n",
    " P = \\min( p(\\beta_1^g < 0 | y, x^g), p(-\\beta_1^g < 0 | y, x^g) ).\n",
    "$$\n",
    "\n",
    "In other words, we just count how many of the posterior samples are on the wrong side of 0 for the effect size. Comparing the effect size with the P-value is known as a _volcano plot_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pval = np.sum((beta < 0), 0) / beta.shape[0]\n",
    "pval = np.minimum(pval, 1 - pval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAADgCAYAAABYS05fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX98VOWV/z9naEggk8kMP/IDhpAhWDY/dFWgC5oqoVbQ\nfmGnalcSBKNZrJVkv+5X10Wu31LWjYpd+1VS0GLBGMsMrUIjVbetmth22rpbqgLRqFgiCyKVQjME\ncciv8/1j5o4zyczkJrkzk5k579freTH3uc9zz7n3DifPnOc85yFmhiAIgqAfhngrIAiCkGyIYRUE\nQdAZMayCIAg6I4ZVEARBZ8SwCoIg6IwYVkEQBJ0RwyoIw4CIComIiegL8dZFGLuIYRXiAhF9SESf\nEVEXEXUS0e+I6HYiMgS0aSSibiI6G1D2h7neIiLq97XpIqL3iOiW2N2RIHyOGFYhnixj5iwAMwE8\nBOBfAWwf0OZhZjYGlL+NcL3jzGwEYPJd60kiKomK5oIQATGsQtxhZjcz7wVwI4CbiahslNdjZm4G\n8FcAgwwrEd1IRPsG1P0zEe31ff4aEb1JRGeI6CgRfSecLN/I+6qA4+8Q0Y8Cjhf4RuOdRLSfiBYF\nnKsmosO+EXYHEa0czX0LYwcxrMKYgZn/G8AxAF8ezXWIyEBEXwdgBnAwRJOfAZhDRBcE1FUBcPg+\nfwpgta//1wB8i4jsI9BjOoAXAfw7gEkA7gawm4imElEmgM0ArvGN2i8D8NZwZQhjEzGswljjOLxG\nSOVu32hPLU9H6DuNiDoB/AXABgCrmPm9gY2Y+RyA5wFUAoDPwP4NgL2+868x80Fm7mfmAwCcAK4c\nwb3cBOAlZn7Jd62XAewDcK3vfD+AMiKawMwfM/PbI5AhjEHEsApjjekATgcc/wczmwPKzRH6Hve1\nmcTMFzPzLgAgoicCJr/W+9o64DOs8I5Wm30GF0T0d0TUSkQnicgN4HYAU0ZwLzMBfCPwDwOAcgD5\nzPwpvK6P2wF8TEQvEtHfjECGMAYRwyqMGYhoPryG1aXndZn59oDJrwd81S8DmEpEF8NrYB0BXRzw\njl5nMHM2gCcAUJjLfwpgYsBxXsDnowCeGfCHIZOZH/Lp9Qtm/iqAfADvAnhylLcqjBHEsApxh4hM\nRPS/AOwC8CNmDuUX1RVm7gHwLIDvwut6eDngdBaA08zsIaIvwTuiDcdbAFYQURoRzQNwQ8C5HwFY\nRkRLiGgcEWX4wsKsRJRLRH/v87WeB3AWXteAkASIYRXiyc+IqAvekZ0C4HsABsae3jMgjvUvOsp3\nALgKwLPM3BtQfweAf/Pp9m0AP4lwjf8LoAjeCISNCBj5MvNRAH8PYD2Ak/De57/A+//OAOD/wOtT\nPg2vD/dbutyVEHdIEl0LgiDoi4xYBUEQdEYMqyAIgs6IYRUEQdAZMayCIAg6I4ZVEARBZ5Iyp+TS\npUv55z//ebzVEAQh+Qi3UCSIpByx/uUveoY6CoIgDI+kNKyCIAjxJKkMKxEtI6Jtbrc73qoIgpDC\nJJVhZeafMfNt2dnZ8VZFEIQUJqkM60hxOp0oKyvDuHHjUFZWBqfTGTd5sdYlFXE6nZgxYwaysrJg\nMBiQlZWFiRMnwmAwYNKkSVi6dKn/HWRmZvrbZWdn+/8dP348xo8fD5PJBIPBgIkTJ4a8ntpnxowZ\nyMnJ8bex2WwR363T6YTNZvNfLyMjw9934sSJQfJUHVRZ06ZNw+TJk0PqPXnyZNTV1QXdfyhdtHwP\nE/G7GjOdmTnpyty5c1krDoeDbTYbt7S0cHd3N7e0tLDNZmOHw6H5GsMhkrxY65KKOBwOzsnJYavV\nGvScCwoKuKmpiVtaWjgvL4/tdjs3NTVxQUHBoPehKArn5eVxZmYmK4rCOTk5g9pZrVbOyclhRVF4\n2rRpnJWVNahNQUFByHfrcDgGtc3Ly2NFUbipqYmnTJkS8rzdbmebzcZ2u53z8vLC6p2enh5RFy3f\nw0T8ruqksyYbFHcjGI0yHMNaWlrKLS0tQXUtLS1cWlqq+RrDIZK8WOuSipSWlnJRUVHE59zS0sIW\niyXi+2hpaWGz2RzxekVFRUFttb7bob4j4eRZLJagf8PpbTKZIuqi5XuYiN9VnXTWZIOSKrsVES0D\nsGz27NlrDh06pKnPuHHj4PF4kJaW5q/r6elBRkYG+vr6dNcxkjwAMdUlFRk3bhyAyM+5p6cH6enp\nIKKw7Twej79NpOup59LT03H+/HlN73ao70g4eaqMSLK06KLl/0Ss/9/ogU46p14cK49g8qq4uBgu\nV3DCepfLheLiYr3VG1JerHVJRYqLi2Gz2SI+Z5fLBbPZHPF9uFwuZGdnR7yezWYLaqv13Q71HQkn\nz2w2B/0bTu+srKyIumj5HibidzWmOmsd2iZSER+rEA7xsYqPVXysMTCszN4HXlpaygaDgUtLS6P+\n5YgkL9a6pCIOh4OtVisbjUYmIjYajTxhwgQmIrZYLLxkyRL/O5g4caK/nclk8v+blpbGaWlpnJWV\nxUTEEyZMCHk9tY/VauWpU6f62xQWFkZ8tw6HgwsLC/3XS09P9/edMGFCkDxVB1VWfn4+T5o0KaTe\nkyZN4tra2qD7D6WLlu9hIn5XddBZfKyCIAg6Iz5WQRCEeJBUhlUQBGEskFSGVXIFCIIwFkgqwyqu\nAEEQxgJJZVgFQRDGAkllWMUVIAjCWGDMG1YiKiaiJ4joOSL6VqS24goQBGEsEFXDSkQ7iOgTImob\nUL+UiN4jog+IaF2kazBzOzPfDuAfAFweTX2F1ERNJaemDVRTytXV1YVNMbd06dKgdHxLly6N4x0I\nY41oj1gbAQR944hoHIAtAK4BUAKgkohKiOhCInphQMnx9VkO4EUAL0USJq4AYbg4nU4oigK73Y7C\nwkLs3r0bHo8HDQ0NeO6552C32/3HiqLA6XRi6dKl2L9/P5qbm3H+/Hk0Nzdj//79YlyFz9G6RGuk\nBUAhgLaA44UAfhFwfC+AezVe60Ut7Ya7pFVIXdRUclpSyqnH4dLumUymmOgsxJWxsaSViAoBvMDM\nZb7jGwAsZeZ/9B2vAvB3zFwbpv8iANcBSAdwgJm3hGl3G4DbAKCgoGDukSNH9L0RISlRU8mpKfUi\npZRTj5k5ZNq99PR09Pf3x/wehJiSHEtamfk1Zv4nZv5mOKPqa7eNmecx87ypU6fGUkUhgVFTyWlJ\nKaceh0u7l5WVFROdhbFPPAzrRwBmBBxbfXWCEHMURUFNTQ3sdjtqamrQ2tqKnp4etLa2oqqqCna7\n3X9cU1MDRVGwcOFCVFVVDWq7cOHCeN+OMFbQ6jMYacFgH+sXABwGYAMwHsB+AKV6yhQfqzAc1FRy\natpANaVcbW1t2BRzS5YsCUrHt2TJkjjegRBD4u9jJSIngEUApgD4M4ANzLydiK4F8CiAcQB2MHO9\nnnLnzZvH+/bt0/OSgiAIgEYf6xeiqQEzV4apfwlDhE4JgiAkKmN+8krQRiLu8R5vhvPMAhcRWCwW\nGAwGZGVlYcaMGfKshcFo9RkkUkk1H2si7j8Ub4bzzNS2iqIM6lNQUMA5OTnyrFMHbXNLWhsmUkk1\nw5qIe7zHm+E8s6EWERQVFcmzTh3iP3kVL1Jt8ioR93iPN8N5ZloWEQCQZ50aJMcCAWFoEnGP93gz\nnGc21CICm80mz1oIRuvQNpFKqrkCxMc6fMTHKowQ8bGmEom4x3u8CXxmFouFiYjz8/M5MzOTiYjN\nZrN/EYDRaAyqV+usVqs869RCDKsgDEWokWtBQQE3NTVxS0sL5+Xl8be+9S0ZnQoqYlgFYSiGig5o\naWlhi8Xi/ywRACmPJhuUVFEBRLQMwLLZs2evOXToULzVERKAoaIDAtMBSgSAgFSMCmDZ80oYgsDV\nVpMmTQIzo7i4OGj1lMvlwsSJE1FWVoaNGzfCbDb76/Pz8yUCQBiSqOYKEISxhLoNy/bt21FeXg6X\ny4WamhqsWLEC69atQ29vL6xWK6qrq7F161ZYrVZUVVVh3rx5/tSAXV1d+PrXvx7vWxHGOlp9BolU\nxMcqhCKSP7WlpYWNRiMXFRUFTU6p9RaLhRVFER+rkHo+VpVUW3klaCOSP9Xj8SA9PT3kliuqX1Vd\nfQWIjzWFST0fqyBEoqCgIOxqK5fLhczMzJD+VnVllayyErQiPlYhJXA6nTh79iyqq6vR2Njo97FW\nVVVhwYIFQX7V6upqv7919erV6O7uxooVK1BdXQ2Px4P7778/3rcjjHHEFSCkBGVlZWhoaMCJEydQ\nX1+P9vZ22Gw2nD17FufPn8fWrVtRWenNy97a2orly5cjNzcXx48fh8fjQWZmJsxmMx5++GF/OyEl\n0eQKSCrDKnGsQjgi+VeBYJ9poN9VMoQJA0g9HytLHKsQhkiZqWw226B6tf3EiRNlpwBh2CSVYRWE\ncCiKgurq6qAtq6urq3Hq1Cm43e6g+ptvvhl2ux2rV6/G1q1bsXfvXhgMBtx5551iXAVtaI3LSqQi\ncaxCKBwOBxcWFgZlq1IUJSjLVVFREWdkZPD06dMHxbNKDKvAEscqk1dCZML5XYeKZxWfa0qTej5W\nQRgO4fyu2dnZslOAMCqSyrAS0TIi2uZ2u+OtijAGGbjddUVFBWpqaoL8q5WVlZg/f37I+q6uLiiK\nEu/bEBIBrT6DRCriYxUGEm4rliVLlrDZbPbvvKAoChcUFHBtbW2Q3zUzM1MSXAvMGm1Q3I1gNIoY\nVmEg4RKwmEwmzs3NDTlRJfuHCSGQyStBUIm0QGDatGno7u7Go48+isrKSv8EFhGhuLgYiqLIaitB\nRdPkleQKEFICdaKqoqLCX6cuBGhoaMDy5cuxdu1aAEBeXh5KSkrQ1tYWL3WFBCepDGvAktZ4qyKM\nMdRRp9PpDEpyXV9fj/Lycpw7dw6vvPKKP9HKo48+Gm+VhQRGXAFCyjB58mT09PTg7NmzyMzMxIQJ\nE/DYY48hLy8PdXV1aGtrQ2trK2699VZ0dHTEW11hbCKuAEFQcTqdMBgMsFgseP755/2j1tWrV8Pj\n8WDz5s0AgPLycvzP//xPnLUVEp2kimMVhIGosasrV65Eb28vVq1ahYqKCqSlpaGiogJNTU0wGo3+\nySmXywWj0Sg5AYTRoTV8IBEKgGUAts2ePVuPsAohwQkXuxoYOtXd3c0GgyHoPBFJiJUQDoljFVKb\nSJsHBh5bLJagBQJqP0m4IoRADKuQujgcDjYajX6DqY4+B45Q8/LyWFEU/7G6WEBtJwgD0GSDkmry\nSsKtBMDrV1UUBXv37g0KrQK8MapGoxHp6emYNm0abrjhBjQ3N+PBBx+EzWZDWloaKisr0draKglX\nhJGj1QIDKAdwi+/zVAA2rX1jXWTEmtqEcwEUFRX5facOh4OnT58e5H8tKCjgpqYmWcYqREKTDdIU\nx0pEGwDMAzCHmb9IRNMAPMvMl0fP5I8ciWNNbSLlWd25c6c/AmDGjBlIT09HR0cHCgoK0Nvbi+PH\njyM7OxtbtmyRZaxCKHSNY/06gEsAvAEAzHyciLJGqFjUEFeAAIRfvlpSUhJkLNUdWEPlDxCjKowG\nrXGs3ewd2jIAEFFm9FQKhojsRPQkEf2YiK6O1JZlM8EgBuYfTZXYTEVRBuVTrampgaIoQc8kOzsb\nGzduDOqr5g8QhFGhxV8A4G4APwBwGMAaAL8HUKeh3w4AnwBoG1C/FMB7AD4AsE6jDhYA27W0FR+r\nthjOZCZwHys1KiDUMxkYFZCXl8e1tbXxVl8Yu2izmZobAl8F8F0A/wHgqxr7XAHg0kDDCmAcgD8B\nmAVgPID9AEoAXAjghQElJ6DfIwAuHUKeLBDwoSWGM9WIlJN1YKLrVPkDJAwbfQ3rSAuAwgGGdSGA\nXwQc3wvg3gj9CcAmAFdplSkjVvbHagaSyrGZkeJaiWhQoutU/gMkRESTDdLkYyWiLiI64yseIuoj\nojNa+oZgOoCjAcfHfHXhqANwFYAbiOj2CDreRkT7iGjfyZMnR6ha8hBuo7xU9B8GxrV6PB40NDTg\nrrvuQl1dHVwuF2bNmuX3vwLeRCzt7e1x1lpIaLRaYA4eQdoBPKSxfSGCR6w3APhhwPEqAN8frh6R\nSiqOWAN9ilarladOncq5ubkp62NVcTgcQUtW1ftvaWlhs9nMOTk57HA4gkapMmIVIhBdVwCANzW2\nG2hYh+UKGElJNcMaOCnT1NTEBQUF/s9FRUVMRFxYWJgyRjXwj4zZbA6anFL/uKhukcLCQmb+3E2i\nLhRIlWclDBv9DCuA6wLKDQAeAvB7jX0HGtYvwBtdYMPnk1elWhXWUlLNsAZOyqT6pFVtbS3n5+eH\nHamrz0L9V/U5t7S0sNFoZKPRKEZViISuhvWpgPIkAAUBM/YR+jkBfAygB15fao2v/loA78MbHaBo\nVVZrSTXDGjhRlcqTVg6Hg81mc8Q/LOqzsNlsQZmsCgoK/MeCEAH9DGuilVQzrDJi9aKOQCP9YVHT\nBCqKwrm5uUxEXFRUxJMmTUpJH7QwbEZvWAE0ANgcrmgVEuuSaoY1nI811Sat1AmqcH9Y1LSAROSP\nWbXZbNzU1MREJAsDBC3oYlhvjlS0Col1STXDyjw4KqCwsHDQTHiyo45EB66uUo2p0WhkRVGC+qhZ\nr4qKilJiVC+MGk02SHZpFZIGg8GA7OxsuN1umM1muN1umEwmuN1uvPrqq7Db7ejq6oLZbMbKlSvR\n0NDgz3r19NNPo7q6Gn19ffG+DWFsoym7ldYFAlOJ6D+I6CUialHL6PQTBP1wOp3Iy8vDnj17cP78\neezevRu5ublYtGgRsrOzcd1116Gurs5/7rnnnvMvEJg5cyasVmtKLp4QooSWYS2AXwKoAdAO4Ep4\nk6ts0josjnVJRVdAqhPOt6rGsU6fPj3kOTU6IFX80MKo0c8VQER/ZOa5RHSAmS/y1f2BmedHz+SP\nHHEFpB7hkltnZGTA4/EgIyMj6Ge+6gIAgJKSEiiKIjlYBS3o5wqANw4VAD4moq8R0SUAJo1ILUHQ\nGafTCZPJFDY3gsvlgs1mG3QuKysLO3fuRFtbmxhVQVe0GtZ/J6JsAHfBm5v1hwD+OWpaCcIwqK+v\nx9q1awclt165ciXsdjtuvvlmuN3uoHNVVVVYvHgx6uvr462+kIRo3Zrlv5jZDcANoGKoxoIQS9rb\n2/Hmm2+itLQUdXV1aG9vR3FxMU6cOIHm5mZs2rQJALBmzRocPnzYHxXwve99DxkZGXHWXkhKtDhi\n4V1+qk5gWbQ6cGNdIImuUxKLxRJy4spisQTVDcy9mior0gRd0WSLNLkCmPmLAO4DUArgj0T0AhHd\nFBVLPwpY9rxKSTo7O1FVVTXop35nZ2dQO5fLhezsbNTX1wftgyUIuqPVAvPno8IpAJoA9A23b6yK\nhFulFuqKq8A9rhRFYZPJNGhpLxH5V2FJeJUwAnQNtzLBuwX2CgBFAH4K4CfM/MfomfyRI+FWqUW4\nUKv09HSUlJT4fa52ux27du0CAGRkZKCtrS1eKguJi67hVvsBXAzg35j5i8z8r2PVqAqpR7htaLKz\ns9HQ0ODfjuWZZ56B2+1GV1eXuACEqKI1KmAW+4a2RJTHzCeiqJMgDIuKigpUVVXB4XCgvLwcLpcL\n1dXVuOmmm3DrrbfiyJEjyMzMRG9vL6ZMmYKHH35Y4laFqKLJsHKwv+AleLe0FoS443Q6sXfvXtTU\n1PhDrUwmE/r7+3HZZZehoaEh3ioKKciws1sR0ZvMfEmU9NEF8bGmDmVlZWhoaEBFxefh1a2trViz\nZg36+vrQ0dERR+2EJERXH2sgT46gT0wgomVEtM3tdsdbFSFGtLe3o7y8PKiuvLwcHR0dOHLkiH9L\na0GIJRENKxFlENGdRPR9IvomEX2BmbfGSrnhwhLHmnKEm7iy2WyYNWsW1q9fj7KyMowbNw5lZWVi\naIWYMNSI9WkA8wAcBHANgEeirpEgDANFUVBdXR20OKC6uhpdXV3YsGEDjhw5EhRupSiKGFch6gw1\neVXCzBcCABFtB/Df0VdJELThdDqxfv16HD16FMuXL0dvb6//3JQpU/Dee+9h1qxZaG9vh8vlQk1N\nDaqqqlBfXy9RAUJUGWrEqqYLBDP3RmooCLHE6XRi3bp12LFjB86fP4+9e/ciJycHRqMR27ZtQ1NT\nE7Zv345rrrkGaWlpqKiowPbt29Hc3Iz29vZ4qy8kORGjAoioD8Cn6iGACQDO+T4zM5uiruEIkKiA\n5CdSNIC6qqq1tRV1dXX+FVZq4uvi4mJZdSWMFE1RAUm1mSARLQOwbPbs2WsOHToUb3WEKBJpxwAA\n6Ovr8x+rOwe0trbiuuuuw9atW8UVIIyUqIVbjVkkKiB1iBQNoG4KqB4HJr6+6aabxKgKUUfrklZB\nGFOo0QCNjY1By1jPnTuHDRs2oLW1FZWVlTh79qz/5/8jjzwiRlWICWJYhYRENZCBuQDGjx+PiRMn\norq6GiaTCZWVlfjGN77hz7sqRlWIFUnlChBSi8rKSnR0dKCkpAR79+7FqVOncPToUfT19WHPnj34\n9a9/7Y8GkL2thJiiNXFrIhT4tmaxWCwMYFD56KOPmJl5w4YNcl7Oy3k5P5Lz+iW6TjQk3Cq1CBd6\ndeutt6Kjo2NQ2JUgjILUiwoQkh+n0zlo7b+iKFi1atWgZa3nzp3DfffdJ3tbCTFHJq+EhEE1otu3\nb/dHAtTU1OCLX/wizpw5g6985SvIzMyE2WzGww8/jLy8PFx//fXYsmWLTFwJMSWpRqySNjC5qa+v\nx/bt21FRURG0TPX111/H888/71/aajB4v9bl5eVwu91iVIWYk1SGlWWBQNLidDpx5MgRXHXVVUHp\n/8rLy9HV1RVkbBsbG7F+/Xq4XC7/YgFBiCXiChDGPKoLYO/evUEuAADIy8tDQUFBUPvy8nIcOXIE\nK1euxCOPSKZLIQ4MN6RpLBf4wq1mz57NQvJQWlrKLS0tQXUtLS1cVFTEVquVp0yZwg6HI+ic0Wjk\n2traWKsqJD/abJHWholU5s6dq8cDTHkcDgeXlpaywWDg0tLSIOMVSwwGA3d3dwfVdXd3MxFxTk4O\nK4rCRUVF3N3dzS0tLVxQUMBTp06Ni65C0qPJBiWVj1XQD/Xnd0NDAzweDxoaGuKWfT9cwpWSkhLs\n2rULzc3NOHz4MDIyMrBmzRqcOXMGjz32WMz1FAQ/Wi1wIhSIK0A3wv38Li0tjbkuDoeDbTYbt7S0\n+EelNpuNHQ4Hd3d3s8FgYIvF4h9ZE1HMdRRShuRwBQCYBWA7gOe09hFXwOgJ/Pkd6BIwGo1xcQk4\nHI4g46nq0NLSwtnZ2X5/aryMv5AyxN+wAtgB4BMAbQPqlwJ4D8AHANZpvJYY1hiijlgjjRZjjcPh\n4Pz8/CBd8vPz2W63+/WNl25CyjAmDOsVAC4NNKwAxgH4k28kOh7AfgAlAC4E8MKAkhPQb0jDKq4A\n/VANalFR0ZhxCTAzE9GgCTV1IiueE2xCyhB/w8peY1c4wLAuBPCLgON7Adyr4ToyYo0BgT/7rVYr\nExF3d3cPihCIpR8zULbFYmFFUYLOy89/IYaMWcN6A4AfBhyvAvD9CP0nA3jCN8oNa4AB3AZgH4B9\nBQUFuj/NVCDUz36z2cyKogyqz8/Pj8no0OFwcE5ODhcVFbHBYOCioiK/TvF2TQgpSXIY1pEUGbGO\njFCRAIqicHZ2dtzcAVarlQsKCoKMekFBAWdmZsY9vlZISTTZoKjnYyWiQgAvMHOZ73ghgO8w8xLf\n8b0AwMwP6iVT8rGOjHA7n6anp+P8+fMhd0RVd0CNFhMmTMD06dPR0dGB4uJiKIqCvLw8LF++HF1d\nXVGVLQghGLP5WP8A4AIishHReAArAOyNgx7CAMIF4pvN5pD10U5w4nQ6kZ2djSeffDJokcKxY8fw\n6aefRlW2IIwKrUPbkRQATgAfA+gBcAxAja/+WgDvw+s3VfSWK66AkREutKq2tjaoXlEUNpvNUf8p\nbrVa/b5VVY6aI6CwsDAqMgVhCMaGjzUeRQzryAmXH0CtJ6JBsaTRmDxyOBycm5s7SE5TUxMTkfhV\nhXghhlXQF3X1UywmsiJltJLRqhBHxLAK+qG6CdSlrtGOa42U0UpGq0Ic0WSDJLuVoAl1W5Ti4mJs\n3LhxUOarvLw8XTNfZWdnh5wwM5lMstWKMPbRaoETqciIVX8CR6pmszmq7gCHw8EZGRkhfaySuUqI\nM5psUNTjWOOBxLHqz6RJk7B7925UVFSEjXfVI67V6XRi3bp1aGxsxLFjx7Bx40YcPnwYM2fOxMqV\nK9Hc3Iy2trbR3o4gjJQxG8cqJCCdnZ2oqalBa2tr2HhXPeJa6+vr0djYiIqKCqxatQoffPABXn31\nVYwbNw47duyAoiijliEI0UZGrIImysrKYLfbsXXrVnR2diIvLw87d+70b+6nbtw3Wv9npNVfANDf\n3z+q6wvCKJERq6AfiqLgmWeewR133IHCwkLceuutqKurQ0ZGBq677jpcf/31ukwqhRsNz5o1CyUl\nJaO+viDEAhmxCpqZPHky+vv74Xa7YTab4Xa7YTKZYDAYcOrUKV1kBPpY1dFwdXU1PB4PHn30UYkI\nEOJN6o1YiWgZEW1zu93xViUp6e7uxp49e9Df34/Tp0+jr68Pe/bsQXd396iv7XQ6UVZWhptuugkG\ngwE33ngj0tPTsXz5cvT394tRFRILreEDiVQk3Co6qEmvA1GD9kfDWNr+RRCGQBYICPoyc+bMkP7P\nmTNnjuq699xzDwwGA6666ipccsklOHHiBLZv3476+vpRXVcQ4oUYVkEzDzzwAKqrq9Ha2oqenh60\ntraiuroaDzzwwIiv6XQ60dPTEzI1YHt7u47aC0Ls+EK8FYgVPT09OHbsGDweT7xViSsZGRmwWq1B\n4UxaUX2cdXV1aG9vR3FxMR566CEA3nAstU5RFM3+0Pr6ejidTlRUVAAAKioqsH37dqxZsybq+V4F\nIWpo9RkkUgnlYz18+DCfPHmS+/v7R+RYSQb6+/v55MmTfPjwYd2u6XA4Qm6dotU/KslWhARDfKyB\neDweTJ5Z8sA8AAAPOklEQVQ8GUSaoiWSEiLC5MmTdR21r1+/3r9SKi0tDRUVFWhsbMRdd92lqX9u\nbm5Yv61EAQiJSsoYVgApbVRV9H4GR44cQXl5eVBdeXk5Tpw4MWS2q7q6Onz22WeD/LarVq0ald9W\nEOJNUhlWPeNYe3t7cfDgQRw5ckQHzbz8+c9/RlVVFWbNmoW5c+di4cKF+OlPf6rb9eNBZmZm2JVS\nQ83q79y5E3v27MFDDz3kX8W1Zs0adHZ2ymhVSGy0+gwSqYTysb7zzjuD6k6dOsV79+7lX/3qV9zb\n2+uv37NnD0+fPp3nzJnDU6dO5cWLF/OxY8eGcL1Epr+/nxcsWMCPP/64v+7DDz/kzZs3j+q6IyHU\nsxgpVqt1UHo/k8nEmZmZbDQa/Vu5hPKXRisuVhCiiCYbFHcjGI2ixbBu3ryZs7Oz+eqrr+aLL76Y\nZ8+ezW1tbXzgwAHOycnh3/3ud8zs/Y++YcMGnj9/Pvf39/PRo0f5m9/8Js+aNYsvueQS3rx5c5BR\nDscrr7zCV1xxRchzvb29fPfdd/O8efP4wgsv5CeeeIKZmVtbW/nKK6/k66+/nufMmcNVVVX+ybd9\n+/bxFVdcwZdeeilfffXVfPz4cWZmfuyxx7i4uJgvvPBCvvHGG0PK09OwOhwOzszM9G/6N2nSJDab\nzYMmtKZPnx5kXB0OB5tMppB5XS0Wi276CYLOiGENJNCY/Pa3v2Wr1codHR3+uh07dvAFF1zAa9eu\n5Y0bNwb17evr4wsuuIBffvllnjlzJt9zzz38zjvv8GuvvcaXX34533HHHRHeg5fHHnuM77zzzpDn\nfvCDH/D999/PzMwej4fnzp3Lhw8f5tbWVjaZTHz06FHu6+vjBQsW8G9+8xvu7u7mhQsX8ieffMLM\nzLt27eJbbrmFmZnz8/PZ4/EwM/Nf//rXkPL0NKzMzLW1tf4NBi0WCxcVFUVMhK1GEtjtds7Lywsy\nwHl5eVxbW6urfoKgI5psUMrEsQby1FNP4c4770RhYaG/rrq6Go8++ijeffddXH755UHtDQYDioqK\n8Mwzz2DRokXYtGkTAG8mppdeegk2mw3r1q3DjBkzNOuwdu1auFwujB8/HjNnzsSBAwfw3HPPAQDc\nbjcOHTqE8ePH40tf+hKsVisA4OKLL8aHH34Is9mMtrY2fPWrXwUA9PX1IT8/HwBw0UUXYeXKlbDb\n7bDb7SN+RsOhoaEBl112Gex2O7q6uuB2u0NOaKkB/4GRBHV1dbj++uvR2dmJrKwsrF69Gg0NDTHR\nWxCiRVJNXmnlzJkzyMnJCaojIuTk5GDWrFnYvXt30LmTJ0/i9ddfx+nTp3HNNdcEnTOZTFiwYAHe\nfPPNiDJLS0vxxhtv+I+3bNmCV199FSdPngQzo6GhAW+99RbeeustdHR04OqrrwYAfx5SwJurtLe3\nF8yM0tJSf/uDBw/il7/8JQDgxRdfxNq1a/HGG29g/vz56O3tHf4DGgGVlZXo6uqC2WyGzWbzT2ip\nyVUyMjJgNBpRV1cXFEnQ0NCA06dP4/z58+jq6hKjKiQFSWVYtUYFLFmyBI2NjUFJkw8fPox9+/bh\n29/+Nt59913ccsst+NWvfoWf/OQnWLx4MWprazFnzhzs378/6Fp9fX1oa2sLGv2GYvHixfB4PHj8\n8cf9defOnfPr8/jjj6OnpwcA8P777+PTTz8Ne605c+bg5MmT+P3vfw/Au6rs7bffRn9/P44ePYqK\nigps2rQJbrcbZ8+ejaiXnsycORNXXnklTp06herqatx3331Bmw42Nzfj2WefBeAd7QeGY7lcLmRm\nZsZMV0GIKlp9BolUhvKxejweXrRoEVdUVHBjYyNv2rSJp0+fzlu3bmVm5tOnT/OGDRt4wYIFvHTp\nUv7xj3/M/f39fOjQIZ4yZQo3NzdzX18fd3Z28h133MGLFy/W5Jw5fvw433jjjVxYWMjz58/nRYsW\n8a5du7ivr4/vvfdeLisr49LSUl60aBF3dnZya2srf+1rX/P3X7t2LT/11FPMzPzmm2/yl7/8Zb7o\noou4pKSEt23bxt3d3Xz55Zf7r/Pggw+G1ENvH6uKw+Hg3Nxczs7O5szMzLCTU6Wlpf4VWk1NTf7P\nVqs1KnoJgo7I5FUgA42Jx+Php556ilesWMG33347v/7665qe6muvvcYXXXQRWywWNhqNXFlZyadP\nn9bUd6wQLcPK/Hn4VW5ubtjlqgaDgZm9RtZoNLLRaOScnBxZwiokAmJYA9HTmPT39/OJEyf4zJkz\nul0zlkTTsKq5VXNzcyOOWJk/N7IWi0WMqpAoaLJBSeVjjRVEhNzcXGRlZcVblTFHZWUl6uvrkZaW\nhq6uLlRWVgYtV62pqfHvtOpyuZCVlYUtW7bISishqUiqPa+IaBmAZbNnz15z6NChoHNqSjshds/C\nZrPhs88+g9FoREdHB0wmE9auXYsNGzbA5XKhqqoK3d3duu2XJQgxIPX2vGLmnzHzbdnZ2eHOx1ij\nsUcsn8EDDzwAZkZPTw9eeeUVbN68Gbt27UJ6ejpWrVqF/v5+fP/734+ZPoIQK1JmgUBGRgZOnTqV\n0qkDmRmnTp1CRkZGTOSpP+/vueceLF++HJ9++ilMJhMAwGw2DyshtiAkEknlClAJtf217CDgZTQ7\nCAiCoM0VkFQj1gAf66BzaWlpsNlssVdKEISUI6V8rIIgCLEgqQyrIAjCWCCpDKueOwgIgiCMlKSc\nvCKikwCG2lNlCoC/xEAdkT125Ivs1JOvt+y/MPPSoRolpWHVAhHtY+Z5Ijt15Ivs+JCK955UrgBB\nEISxgBhWQRAEnUllw7pNZKecfJGdevLjIjtlfayCIAjRIpVHrIIgCFEhqQ0rEU0iopeJ6JDvX0uI\nNnOI6K2AcoaI7vSd+w4RfRRw7lo9ZfvafUhEB33X3zfc/qO47xlE1EpE7xDR20T0vwPODfu+iWgp\nEb1HRB8Q0boQ54mINvvOHyCiS7X21UH2Sp/Mg0T0OyL624BzIZ+/zvIXEZE74Hl+W2tfHWT/S4Dc\nNiLqI6JJvnOjunci2kFEnxBRW5jz0XznQ8mO6jsfEq0ZsROxAHgYwDrf53UANg3RfhyAEwBm+o6/\nA+DuaMoG8CGAKaPVfbh9AeQDuNT3OQvA+wBKRnLfvuf2JwCzAIwHsF+9VkCbawH8J7xJLBYA+C+t\nfXWQfRkAi+/zNarsSM9fZ/mLALwwkr6jlT2g/TIALTre+xUALgXQFuZ8VN65RtlRe+daSlKPWAH8\nPYCnfZ+fBmAfov1XAPyJmYdaXBAN2Xr2H7IvM3/MzG/4PncBaAcwfZg6qnwJwAfMfJiZuwHs8ukw\nUKcm9vI6ADMR5WvsOyrZzPw7Zv6r7/B1ANbh3uBo5Eep70j6VwJwRjg/LJj51wBOR2gSrXc+pOwo\nv/MhSXbDmsvMH/s+nwCQO0T7FRj8xavz/aTYMZyf48OQzQBeIaI/EtFto9B9JLIBAERUCOASAP8V\nUD2c+54O4GjA8TEMNtLh2mjpO1rZgdTAO4pSCff89ZZ/me95/icRlQ6z72hlg4gmAlgKYHdA9Wjv\nfaT6jfa+h4ve73xIEj5tIBG9AiAvxCkl8ICZmYjChkAQ0XgAywHcG1D9OID74X0R9wN4BMCtOssu\nZ+aPiCgHwMtE9K7vr3HE/jretxHe/2x3MvMZX3XE+05UiKgC3v9k5QHVQz5/HXgDQAEzn/X5q5sB\nXKCzjKFYBuC3zBw4yovFvceVeL3zhDeszHxVuHNE9Gciymfmj30/QT6JcKlrALzBzH8OuLb/MxE9\nCeAFvWUz80e+fz8hop/C+zPp1wAi9tdDNhGlwWtUdzLzHq33HYKPAMwIOLb66rS0SdPQd7SyQUQX\nAfghgGuY2b/JVoTnr5v8gD9YYOaXiGgrEU3RqvtoZAcw6NeYDvc+Uv1G+841EcV3PiTJ7grYC+Bm\n3+ebATwfoe0g/5PPKKl8HUDIGciRyiaiTCLKUj8DuDpAxnB0H4lsArAdQDszf2/AueHe9x8AXEBE\nNt/If4VPh4E6rfbNFC8A4Pa5K7T0HZVsIioAsAfAKmZ+P6A+0vPXU36e73mDiL4E7/+7U1r6jla2\nT2Y2gCsR8D3Q6d6HIlrvfEii/M6HJpozY/EuACYDeBXAIQCvAJjkq58G4KWAdpnwftGzB/R/BsBB\nAAfgffH5esqGd1Z0v6+8DUAZqr+Ossvh/al/AMBbvnLtSO8b3hng9+Gd7VV8dbcDuN33mQBs8Z0/\nCGBepL7DfM9Dyf4hgL8G3Oe+oZ6/zvJrfdffD+9EymWxunffcTWAXQP6jfre4R2IfAygB14/aU0M\n3/lQsqP6zocqsvJKEARBZ5LdFSAIghBzxLAKgiDojBhWQRAEnRHDKgiCoDNiWAVBEHRGDKswpiFv\nNqbA7GPrfPVfJm9WrreIaAIRfdd3/N0RyFg/zPb/RkRhF2gIgoRbCWMaIjrLzMYQ9U8AcDHzj3zH\nbnjjdfv0kiEII0VGrELCQUT/COAfANxPRDuJaC8AI4A/EtGNRDSViHYT0R985XJfPyMRPUXeXJwH\niOh6InoIwATfyHfnADnjiKiRvHlMDxLRP/vqG4noBiKaFzCSPki+nAxEVEREPydvko/fENHfxPQB\nCXEn4XMFCEnPBCJ6K+D4QWb+IRGVw5vj9DnAP+q82PfZAeD/MbPLt7TxFwCKAfxfeJdVXuhrZ2Hm\n3URUq/YdwMUApjNzma+9OfAkM+/ztYHPBfFz36lt8K4AOkREfwdgK4DFOjwLIUEQwyqMdT4LY/Qi\ncRWAEt/yfAAwkTeL11XwrksHAPDn+TrDcRjALCJqAPAigF+GakREN8KbdPlqn5zLADwbID99mPoL\nCY4YViEZMQBYwMyewMoAQ6cJZv4rebf0WALvOvR/wID0iURUBu+OC1cwcx8RGQB0juCPgZBEiI9V\nSEZ+CaBOPSAi1ci9DGBtQL2awLuHvCkUgyBvaj8DM+8GcB+8o9LA82Z4k4GsZuaTgD9FYAcRfcPX\nhihgvyUhNRDDKox11IkltTykoc8/AZjnm6B6B97RJgD8OwCLbzJqP4AKX/02AAcGTl7Bm9X+NZ+P\n90cIToIOeLcTmQngSVU/X/1KADU+GW9jmNuOCImPhFsJgiDojIxYBUEQdEYMqyAIgs6IYRUEQdAZ\nMayCIAg6I4ZVEARBZ8SwCoIg6IwYVkEQBJ0RwyoIgqAz/x9ZtBiWe0I57gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e3181d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(5, 3)\n",
    "plt.yscale('log')\n",
    "plt.scatter(beta.mean(0), np.maximum(pval, 1./beta.shape[0]), c='w', edgecolor='k', label='Genes');\n",
    "plt.axhline(0.05, lw=1, ls='--', c='k');\n",
    "plt.gca().invert_yaxis();\n",
    "plt.ylabel('P-value')\n",
    "plt.xlabel('Effect size');\n",
    "plt.legend(scatterpoints=3);\n",
    "plt.title('DE P-values');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case we drew 2,000 samples, which puts a limit to the smallest P-value we can observe as 1 / 2,000, causing the plateau in the figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>beta</th>\n",
       "      <th>pval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000060461</th>\n",
       "      <td>1.296872</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000091955</th>\n",
       "      <td>0.463096</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000028936</th>\n",
       "      <td>0.383841</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000063856</th>\n",
       "      <td>0.427333</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000067288</th>\n",
       "      <td>0.365123</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000033307</th>\n",
       "      <td>0.390917</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000028691</th>\n",
       "      <td>0.462913</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000063524</th>\n",
       "      <td>0.521403</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000079435</th>\n",
       "      <td>0.369893</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000049422</th>\n",
       "      <td>0.864710</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        beta  pval\n",
       "ENSMUSG00000060461  1.296872   0.0\n",
       "ENSMUSG00000091955  0.463096   0.0\n",
       "ENSMUSG00000028936  0.383841   0.0\n",
       "ENSMUSG00000063856  0.427333   0.0\n",
       "ENSMUSG00000067288  0.365123   0.0\n",
       "ENSMUSG00000033307  0.390917   0.0\n",
       "ENSMUSG00000028691  0.462913   0.0\n",
       "ENSMUSG00000063524  0.521403   0.0\n",
       "ENSMUSG00000079435  0.369893   0.0\n",
       "ENSMUSG00000049422  0.864710   0.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({'pval': pval, 'beta': beta.mean(0)}, index=idx).sort_values('pval').head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The reason I'm writing about this, is that I had a conversation with [Tomás](https://twitter.com/tomsgoms) about this in relation to our notion of cell types.\n",
    "\n",
    "**It's kind of backwards!**\n",
    "\n",
    "We had the cell types, and then investigated which genes were expressed in the cells. In essence, from a machine learning perspective, we are assessing if the cell type label can predict the gene expression. But what we _want_ to do is investigate how gene expression predicts cell type!\n",
    "\n",
    "So can we do it the other way around? Keeping the notation like above, we want to investigate\n",
    "$$\n",
    "y = \\beta_0 + \\sum_{g=1}^G \\beta^g \\cdot x^g + \\varepsilon.\n",
    "$$\n",
    "\n",
    "Now, if $ \\beta^g $ is positive, the gene will be a predictor for mESC identity, and the magnitude of this will inform about how important it is for determining the cell type. (I think we can still call this _effect size_ in a meaningful way.)\n",
    "\n",
    "Let's refer to this as _reverse differential expression_, and implement it in Stan in this way:\n",
    "\n",
    "    data {\n",
    "      int<lower=0> N;\n",
    "      int<lower=0> G;\n",
    "      matrix[N, G] X;\n",
    "      real y[N];\n",
    "    }\n",
    "    parameters {\n",
    "      real beta0;\n",
    "      vector[G] beta;\n",
    "      real<lower=0> sigma;\n",
    "    }\n",
    "    model {\n",
    "\n",
    "      beta ~ normal(0, 1.);\n",
    "      beta0 ~ normal(0, 1.);\n",
    "\n",
    "      y ~ normal(beta0 + X * beta, sigma);\n",
    "    }\n",
    "\n",
    "After sampling, we can plot the effect sizes of the genes like above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "with open('LR.pkl', 'rb') as fh:\n",
    "    model = pickle.load(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = 200\n",
    "idx = counts.pipe(np.log1p).var(1).sort_values(ascending=False).head(p).index\n",
    "tmp = np.log1p(counts.T[idx])\n",
    "y = (sample_info['condition'] != 'NSC')\n",
    "n = tmp.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = {\n",
    "    'N': n,\n",
    "    'G': p,\n",
    "    'X': tmp.values,\n",
    "    'y': 2 * (y.astype(int).values - 0.5)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "s = model.sampling(data=data, iter=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "beta = s.extract('beta')['beta']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10e2a3cf8>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAADgCAYAAADBheobAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUHVWdL/DvDxITkhDyJkBAImCeC1HQUScjiYggKBBt\nMoRckJDLw0fucHVGjI0XmEU7AYd5hfHyMIAwdEMEg9xhLqg3RM1VR5GLEJIwgxmchEcSGmIgkXRC\nfvePU9VUn9Sus3e9q873s1ZWuk+fc2rXrr13/WrX3rtEVUFERERERO4OKDoBRERERERVxWCaiIiI\niCgmBtNERERERDExmCYiIiIiionBNBERERFRTAymiYiIiIhiYjBNRESlIyLPi8jHik4HEVErDKaJ\nqG15AdsfROQNEXlZRO4UkRFFpysuETlaRNTbnzdEZIuI/LOInNr0vuB++/9uMnznNSKyx3vPdhH5\nmYh8KJ89GpCGf8pzm0REthhME1G7+5SqjgBwAoD3AliS58alIe22eJS3T+8B8EMAK0Xkoqb3fEpV\nRwT+fTHi++7zvm8cgMcAfDfl9BIRVRaDaSIiAKr6MoBH0QiqAQAiMkRE/lpE/tPr5b1ZRA7y/rZe\nRD4ZeO8gEdkmIu/zfv+g14u7XUR+IyKzA+9dLSJdIvJ/AewC8C4RuUhENorI6yLyHyKyIPD+i73t\nvSYij4rIO233SVX/HsA1AK5PGrSr6l4A9wA4QkTGB9L3SRF5MtBzfXzgb1eKyAvefj0rIqd4r98p\nItcF3jdbRDY3b1NETgfwNQB/6vWO/8Z73ZhfRER5YjBNRARARCYB+ASA5wIvLwXwbjQC7GMBHAHg\nf3h/6wEwP/De0wC8oqpPiMgRAB4GcB2AMQD+HMADwQAUwAUALgVwMIBtAP4BwCdU9WAAHwbwpJeu\ns9EIJj8NYDyAn3rbdvE9ABMATHH83AAi8g4AFwLoBfCa99p7AdwO4DIAYwHcAuAh70JkCoAvAni/\nt1+nAXjeZZuq+giAb8DrHVfV94jIcBjyi4gobwymiajdPSgirwPYBGArgKuBxvALNILd/66qr6rq\n62gEded5n+sGcJaIDPN+Px9vB7n/BcC/qOq/qOo+Vf0hgMcBnBHY7p2q+ozX27sXwD4AM0XkIFV9\nSVWf8d53OYC/UtX13nu/AeAE295pz4ve/2Oa9nt74N8lEZ+fJyLbAfwBwCUAOry0AI08ukVV/1VV\n31LV7wDYDeCDAN4CMATAdBEZrKrPq+pvHdIdxZRfRES5YjBNRO3uHK93czaAqWiMCwYavcDDAPza\nDzgBPOK9DlV9DsB6AJ/yAuqz0AiwAeCdAM4NBqsAZgE4LLDdTf4PqroTwJ+iETi/JCIPi8jUwHf9\nfeB7XgUgaPSS2/Lf+2rTfo8K/Lst4vMrVHUUgEMBrAVwYuBv7wTw5aZ9PRLA4V4eXYHGMJOtInKv\niBzukO5QLfKLiChXDKaJiACo6o8B3Angr72XXkGjJ3ZGIOA8xJuI5/OHepwNYJ0XPAKNQPnupmB1\nuKouDW6yafuPquqpaATcGwDcFviuy5q+6yBV/ZnD7s1Fo9f9WYfP7EdVX0GjJ/oaEfEvDDYB6GpK\n3zBV7fE+062qs9AIuhXA9d7ndqJxseKbGLXpkLSY8ouIKFcMpomI3vZ3AE4Vkfeo6j40ArS/FZEJ\nACAiR4jIaYH33wvg4wA+h7d7pQHgn9DosT5NRA4UkaHeBLtJYRsVkUNF5GxvLPBuAG+gMYwBAG4G\nsEREZnjvPUREzrXZGe97v4jG0JUl3j4loqrPojFR8yveS7cBuFxE/shbmWS4iJwpIgeLyBQR+aiI\nDAHwJhoXJ34angRwhoiMEZGJaPRgm2wBcLQ/gbJFfhER5YrBNBGRR1W3AbgLb08yvBKNCYm/EJEd\nAH6EwCQ+VX0JwM/RmAB3X+D1TWj0Vn8NjcmFmwD8Bcxt7gEAvoTG2OZXAZyMRoAOVV2JRm/uvV4a\n1qIxUTLKdhHZCeBpNMZpn6uqtze953/JwHWmV7b4zqBvArhURCao6uNojKO+CY1Jic8BuMh73xA0\nJnG+AuBlNCZB+ksP3g3gN2hMSPwBAvkXwl+Kr1dEnkBEfhER5U1U97t7RkREREREFtgzTUREREQU\nE4NpIiIiIqKYGEwTEREREcXEYJqIiIiIKCYG00REREREMQ0qOgEuTj/9dH3kkUeKTgYRERER1Z/Y\nvKlSPdOvvPJK0UkgIiIiIupXqWCaiIiIiKhMGEwTEREREcXEYJqIiIiIKCYG00REREREMTGYJiIi\nIiKKicE0ERFRScyePRuzZ88uOhlE5IDBNBERERFRTAymiYiIiIhiYjBNRERERBQTg2kiIqKMcAw0\nUf0xmCYiIiKizNX14pLBNBFRhLo2/kRElA4G00REREREMTGYJiIiIiKKicF0SngrOD9lz+uypy+J\n4L7VeT+JiIhsMZim0mKwRkRpYXtCRFlhME1UMJ7kB6pbfhS1P3XLRyKismIwTeSAAQoRpY3tykDM\nD6oaBtNEJZLHSYQnqnrh8SQiKhaD6Rqp0km1Smklylod6kMd9oGIKA4G0xnjCaY6eKzam+vxZ3mp\nBx7HcEnyhXlK7YbBNBGVFk/K1cdjSFlguaIyYTBdAqZGgY0FUTVVse5WMc1ERGXAYJqIiIicleUC\nrCzpqBPmqRsG0xVUpUJepbQSERERuWIwTUSVwAuzdDAfiYjSxWA6AzxZEWWLdYyywHKVH+Y11QmD\n6YrLavIiGzoiiqtd2o922U8iisZgmmqDJzai9NjUp7LVuazSU7b9JDMeKyoCg2mqLTaqbphf9cVj\nW288vlQnVSzPhQbTInK7iGwVkbVFpoPsFVXIy74Wd1nSUVXMPyIiqqqie6bvBHB6wWkgqj0Gq+Xj\nekzqcAzbcZ+JTOpcvuu8b2EKDaZV9ScAXi0yDZSOdqs4ceSRRzwO6WA+VgePVfaYx8Vh3ldD0T3T\nRESpKftwIHobj0n2sshjHjei/ZU+mBaRS0XkcRF5fNu2bUUnZ4CqNCrBdKaZ5ixu2VYlT8kOjye1\ni7KU9bKkg6idlD6YVtVbVfUkVT1p/PjxhaaFjZQd5lP1JTmGZTn+ZUlHHTAvqU7yLs+sP/VX+mCa\niKJVtaGuarrzxDyiLLBcpaOM+VjGNNmoarp9RS+N1wPg5wCmiMhmEVlUZHqIiKgeqn5ypv3leUxZ\nfsjFoCI3rqrzi9y+Db8yrV69ulTfRdnisXqbTV4wv8ohj+PAY90a84jKgG13fjjMg9pOVXscsppI\nWnXMi/iYd9XBY9U+XI+1zftZfrLFYLoiWBGykXUjxONGlD7WKy57l4V233+Kj8F0QcpYacuYJqIs\npdkDxPpDaWOZKh/bY1LmY2dKWx32rSgMpomIiCwxkCCiZgymqSWePCgNLEfMAyIKV/a2oaj0lT1f\nfAymiYiIaq4qQUkZMe+yUad8ZTBN1Mbq1JjZasd9pvJhOSRXwTLD8lMuDKaJiIgokaST2oiqjME0\nERERFYpBN1UZg+k2xsaLqJpYdyktLEuUBMtPA4NpagvtWOHbcZ+JiMqM7XI9DSo6AURV0dPTg2ee\neQa9vb2YOXMm9u3bhwkTJmD16tVFJ41qxD/RslwNxHwhVywzlBf2TBNZ2Lp1Kzo7O7FixQrs3r0b\ny5Ytw/bt27F169aik0ZEVFtp9eSyR5iyxGCayMK2bduwfPlyzJkzB4MHD8acOXNwzz33YNu2bUUn\njahWyhL0lCUdRFXRznWGwzyIYB7C4evt7cWsWbMGfGbWrFno7e3NLE28RUlEdVLGNq2Maaqbdshj\n9kxT27MZwjF27FisWbNmwOfWrFmDsWPH5p1cIsqAf0H9k5/8BDNnzuQQLnLSzr2yZBFMi8i7ReT/\niMha7/fjReSq7JNGlA+bIRzjx4/HokWL8Nhjj2HPnj147LHHsGDBAowfP77AlOeLwQbVRXPgk9ec\niKoEXKzrRG5seqZvA7AEwB4AUNWnAJyXZaLqgg1SNdgM4ZgwYQK6urowb948DBkyBIsXL8aoUaMG\nDAVxUdRJNW6Z5ATMgZrzsaenp+gkUQJx5kTUtX1nXSdyZxNMD1PVXza9tjeLxNQJGyR7RffW2A7h\nmD9/PmbMmIGPfOQjWLt2bexAuihJyiQnYL4tLB87OztZtyvMdU5Endt31nXKSl0vQAG7YPoVETkG\ngAKAiHQAeCnTVFWYX1g2bdoUq0EqOrBsR+0yhCPJSbKICZhlFZaPy5cvZ7BRYa5zIuoccLKuUxbq\nfAEK2AXTXwBwC4CpIvICgCsAfC7TVFVUsLDs2rWLDVJFpD2Eo6ySnCTzmoBZhZ6LKgcbeeev7faK\n7kRwvaCuchlopWqTravQZlC9L0ABi2BaVTeq6scAjAcwVVVnqerzmaesgoKFZdq0aZk0SO3ScOS9\nn1UfwmEjyUkyj9775p6Lc845By+++GIuZcClvJUp2DClO+z1vHuGqtQT5XpBXaYykLYq3amrUhlr\nd3W+AAXsVvN4S0SWAtilqq97rz2RecpKLuxkFSwsnZ2dqTdIdW84gnn65S9/ubb7GZTnRUOSk2Ta\nvfdh+x28GL3//vvR3d2NlStXZl4GXOtVWD4uWrQo92DDlO7nnnsu9PWXX345156hqvVEuVxQVyng\ndJW0rufZplWtjLlKKy/LcEeqzheggN0wj2e89/1ARMZ4r0l2SSq/np6e0JPVqFGj+gvL/Pnz0dXV\nhUsuuSS1oQNFNxxZVshgYDB9+nTcc889oftZp575vC+Okp4kg8FGZ2cntm3bFus4mOpP8GK0q6sr\nt7LeXK9efvllDBs2DBs2bAjdt7B87Orqyv1uhqk9eOONN0Jf3759e6yeobh1Lqonqur1uExDw7LI\ny7h1PY02zWV/6tzbmdb5Ic73JFmtyLS9IUOG1PYCFLALpveq6lcAfBvAT0XkRHiTEdtV84neP/lu\n374dCxYs6C8sEydOxK5duzB16tRUhg4U2XBkfQs+GBisX7/euJ+mRqG58i9evLiQk7XLiaCIi6M0\nhrMkbeRNgfLIkSP7L0ajykDagvXKD/Rvu+22yH1rzsf58+enni6XdPtmzZqFHTt2hL4+fPhw554h\n22Pt0hM1YsSI2N9ZJlkMDXPd5zQCpTSD46Rtmuv24vR2lqVctUpHWucH1+9JulqRaXu7d+8uzQVo\nFmyCaQEAVb0PwJ8CuAPAu7JMVNkFT/TNJ9+LL74Yc+fOTVxYirxNUsQt+GBgYBpvPnLkyNBKunnz\n5v0C/fvvvz/3YSKuJ4Kq9qokbeRNgfKOHTv6ey6ay0BPTw+mTZsGAKmfAIP1Ks8e8aRM7UHwoiT4\nuog49wxFHetWw7JMPVGDBg1qmcd53bUpKrBKa0x7GoFSmsFx0jbNdXuuw23KMlTSJh1pnR9cvyfp\nakVR26vz3CSbYPq/+j+o6loAfwLgv2WWopIKNn6HHHKI8eR73XXXYeXKlRg7dmzqPX953CYxbTvq\nFnyr2+I2goGBaby5qcdNVQek58EHH0R3d3fuAZHriaBMY8jyvLUaFSjv2bMHn/nMZ7Bu3br+uzx3\n3303vvrVr7bsLY4reELOs0c8KVMgMWLEiNDXJ02a5NwzFHWsWw3Lau6JuvDCC/Hmm29aDTfJ465N\nUYGVabtxxrSnESilGRwnbdNst+e3Vxs2bOhvM2zKdBF3A1t1TpnSkVavu+v3hB2DzZs3480336zc\nBO08GYNpEfmo9+M7ReTT/j8ApwB4I5fUlURz4/f5z38e559/fionX5eKlsdtEtO2TbfgbW+LtxIM\nDDo6OnD++efv18NvqqQ7d+4ccAyKCohcTzzNwdBVV12FuXPnore3N9desjxurQYFL5aaA+W77roL\nQ4cOxdSpU3HjjTdi3rx5+PznP48777zTaQy9y8VBcPzrsGHDEp8I8po0ZBq3e+yxxxrbCVPPkGmM\nZFTvt82wLH97U6dOxeDBg/HAAw9g+vTpLfM4j7s2eQRWLu17nDHtaQRKcYPjsH1LOjHTdnvBOyHB\nNqNVB1bedwNtOqdM6bDNy7h3iEzHpPkY9PT04KqrrsJDDz0Ue4K2SxkoyzAcV1E90yd7/38q5N8n\nM05XqTQ3ftdddx0WLVqEuXPnJjr5xqloWd8miRqHGXYLPq3b4s2BwYMPPojDDz98wH6aKmlw4idg\nHiaS1ZWxX/ldx6Q27/Ptt99uNXQm7cYm6a1V14sAf3Juq0DZL+vNF0tA9Bh602oWUenyt3XkkUcm\n2re8Jw2Z2gOXdiJsW1dccQWeeuop9Pb2DpgHEnaXyGZITrCM2ax0lEfvVtaBlWv7HmdMu2vg4pqv\npu8fMmRI6L4BSNThY7O9qAnqreTda2rTOWVKh80k16iJ+/4d4xdffNGp9775GHz96183ttFhkkzO\nLcswnDiMwbSqXu39vzDk38X5JbF4YY3f1VdfjR07dux38nW5CktS0bISdlUadQs+qlcqGPQdeeSR\neOqppyIDwFYBQHMlDd42Dp7wzznnnP47B1nPGg5W/m9961u46KKLnAIxf5/Hjh1rdYLIorFxDSri\nXgSE7XdUoOyz6R21Wc3C5mSbdN+KmjSURNhqJkOHDsUDDzxgnAdiGpZlGpITLGP+xdTixYuNJ9u0\ne7fCViPIOrAyrRJjCpqbx7TbtB+ugYvrhbCpzf39739vLJ9JOnxM+7N7926rOyGtOhqyWtLQVN5s\nOqeaj0PwXNnV1YXx48cb89I0cb/5jrFL733zMdi4caPzRWfcMlD0imVJ2Kwz/WciMlIavi0iT4jI\nx9PYuIicLiLPishzIvLVNL4zC1GNbpKrMNuKlucSMsHGxuYWvKlnPjhj/zvf+Q4OOOCA/pOzawAY\nbKj8xiV427j5hP/ggw+io6Mjl1nDwcp/wQUXYOnSpf3LIboEYrYBbRaNTZygwvUiIMm2TSdA0xh6\n0+u2PY5J9i1Ob2erNevj7EOr749Kc1dX14CeqOA8EH+ZtGCPdXBYlulOQ3MHwfz587Fs2TLj3JK0\ne7f8nvY0hyS0Ylolpvmi2zSm3bb9cFnGLs7FYthQnSyf8BsWiNlMULdZJSaLJQ2jVr+Iat/CjoPr\nudKUL0nvGAePgc2wrLRUdVI+YDcB8WJV3QHg4wDGArgAwNKkGxaRAwH8I4BPAJgOYL6ITE/6vWny\nT0JhtzqDD2mIexVmW9H83gDbCX5JhgEEGxubW/CmnvngjP3rr7/e6TZRkO1kneaJn8uWLctl1nBz\n5Z8/fz7Wr1+P4cOHOwVitgFtWsFaUFRQ4RqI2aTHdts+0wnQdjWLuCuBuOxb3KE+pvLdPHSp1feY\nRN3JMKXZdmnK5gvYww8/3HinIU4HQVq9W8097WkNSTAJy9dgcNN80R02pj3OhZzr0CDXbeTxhF8T\nmwnqNqvEBPe/1byBJEvB+atfRLVvYcfB9Vxpypc05w0153dac3vq9mAX66XxAJwB4C5VfQbpPLTl\nAwCe8x5X3gfgXgBnp/C9qWh+sETzrc40HtJgU9Gae2Cj1lYOW14pznrQLrfgTYFOcDJNkoqd5mSd\nLNhOigTcJiOagg3XxsbmBGs6hgBafjZp42fbUxTW+2YazxtczSLJSiC2+xZsK0y9jqag0VS+9+7d\nm8oJzPT9weUkm9M8efJkq2E1YSsX2XYQZHm3qFVPe1pDEsKYhn01t4H+RTeA0O3GuUjNeqWOrJ/w\nG8VmgnrUOaFVoJxk+FxUPtq0b8HPu54rTfmSxkRqX3COS9TdjCTrl/txiqlNr8KDXQZZvOfXIvID\nAJMBLBGRgwHsS2HbRwDYFPh9M4A/an6TiFwK4FIAmDRpEtauXbvfF02ZMgWDBw/Gli1bQhuOOH+/\n5ppr+hsmALjuuutwyimnoKOjA/feey8AYOfOnQAaY4kBoK+vD2vXru1/3f952LBh6Onpwdq1a/Hq\nq69iypQp2Lt3L8aMGYOLLroI8+bNQ29vL6ZMmYKDDjoIw4cP7//sli1bcP/99/enw28gzzrrLCxZ\nsgQrVqzArFmzsGbNGixYsAB/+MMf8L3vfQ9z5sxBT09P/3rQs2bNwrXXXoubbroJGzZsGJCGPXv2\nYPDgwejr68OePXv6tz169GisWbOmf9tAo0KOHj16wH52dHTglltuQV9fH2699VYsXLhwwGf9XgzT\n96gqRARbtmzZL+9aTdYxfafr8fH/Htx//z1R6Rs5ciQWLVqE5cuXDzgO/vKJfvp6enrw9a9/vb+8\n+Xnvf/+ECRNw3nnnoaOjA6+99hqOO+64/rIQPD4jR47EwoULcccddwzY3vjx40PTZyo/HR0dA/bf\nv+1+9NFH44477sDChQuxcePGlp817f/IkSMjy/+iRYsiy8/OnTtDj8+WLVsGlPtrr70Wc+fOxY4d\nO/rrz5gxY3Duueeio6MDfX19eOihh1ruf9jxN+V1cN+Agb2OADBo0CBccskl2LhxIyZOnNh/V2na\ntGno6+vrP+5R5fuNN97AVVdd1V8eJk6c2F+P/XTs3r07dv0JLicZluYLL7wQd91114DttRo+06p8\nfvSjH8XRRx89oIzt3LlzQP284YYb8PTTT+O1117rryeTJk2yqp8333wznn766f3ahqgApVX+AcDd\nd9/df4eyue6Glc9gnfHz1Q9uwtqrsO9v1fYC+7dvrYLjuO27//1jxozpf7//kCK/zEyfPh0jRowY\ncO5qPj6m/DWdf4KfD2sfx44di+OOOw4/+9nPcOqpp2Lr1q2h+xMc/mGqP9u2bcOKFSv6P9u83Ou+\nffswatSo0PLpX0CG5ePOnTtbtm/Bz5vOlWPHjg3Nv+HDh2PJkiX9MYSfLwCs2i7b4+OfH/bt29ff\ni96cT1/60pdC8zjs/BusJ2Fxit+mN7efCxcuxBlnnDGgjGcR/wXNnDlzv9dCqWrkPzR6r98HYJT3\n+1gAx7f6nMX3dgD4duD3CwDcFPWZE088UfNywAEHaF9f34DX+vr6VET6fz/55JP15JNPbvnztGnT\ndPLkybpq1Srt6+vTVatW6WGHHabTpk1r+VkRCU3HiBEjdNWqVf2vdXd36zHHHDPg/TNmzOh/T3d3\nd2Qa0kp3d3e3jhs3TkVEDzvsMF21apXeddddetRRR8Xa/3Hjxg3YT1XVVatW6cEHH+y8P83pmzRp\nko4ePVpFRGfMmKHTpk2zOp5R3+l/TzDvkuy//3Or7bmWHxFp+T1Rn3VJT1g5mjx5cqz8NpWHcePG\nOe9/3GPb/J6otqJV/bHZH9d9tsmvVvkStt+26XApn63KiW09iapvxxxzTOz8c01TVL6Gfc8RRxzh\n9Lq/3STHJ+6+Rb3f9Zg3v7+7u9u6DXA9VqNHj26ZL8Hj1nyu7Ozs1EMOOWTAOa3VsXJp32zPFUnz\nOO75zReVT8FYw6bsmeKU4GfDzvGTJ0/W7u5uzZlVTNtymIeq7lPVJ1R1u/d7r6o+ZReqR3oBwJGB\n3yd5r5VCmmPCkkwasxlGEJzcEpwsEOyRiTMhwXWyRvOtG39ozGc/+1ns27evf2kelzHgpuEPcR5A\nEUxf0kmRQWG3im3Hnsfhcms6ySOdbYc5tJoAlfSJWkFRvW9ZjMGzyWvXJ3YG18e2ua2ZZFy67XKS\nzfkStt+2Q5GymMnfamWgqInAu3fvdhp6E5Um0wOqbMbMt1qlwvaZAs1DEG3WEw6rG67te5LhYEFJ\nH1dtYjPk0Be1YlDwXBl82q/tw4lch4EG0918roybx0A2j7s35ZPr8JTg95g+q00PY/PPG11dXYn3\nIws2Y6az8isAx4nIZBF5B4DzADxUYHoGSHNMWNYnw2ChDqY7eJKPO27ZpUKGrcftj6nctGkTjj/+\n+Mgx4GGiGnzXxiKYviSTIm25jD3Piqn82EzWSetRvWnO0Ha9OMjjqaGuT+yMmsgXFtAkuSAw1Z+J\nEyc650sWKyEE2Txt0XQRbJoIDACbNm3C0qVLE6+4ZHpAVXBd86gx861WqWjeZ1P7ZuoYMQXggDkQ\nc21Dw97v2lmU5sW1TfpcVwwydULZPJzI364/DCZOuv1zZVQeJ3nqcJIFCkz55Nr5GPwe02dN502/\nXpdNYcG0qu4F8EUAjwJYD2CFNiY3lkLzoPskJw6XiUy2vQfBk2GwUAfXcA2uB53H7GuboClOL31a\nV9i2Ez3SfihKkTOUk/TWuAZPeayb7npxkMdTQ01thWnN9uHDh7ecyGezz7YXBK3unLjkS5YPjbJZ\nT9x0EdyqjKWx4pIpiA2uax61UofLPkfVDdugzr9LtGnTJqs2N267l2Qio837k3JdMSg4eS8qYIy7\nSlAcNhd1NttPukCBKZ9cOx+D3xOMU2zunvl5XjY260zfbfNaHKr6L6r6blU9RlVL13ef1onDZqH8\nqNnErU6GzTN3g2u4+utBmwpsmj10NieGIteRDKYvyVqlgNuJJ+v1bFuJ21tj+qyJy7rpwaUlXcS5\nOHBZhzeuVsMigquKuK7Raxv4ugZDWT9N1ZXNeuKmQDLttfnDhuFEbdtlpQ6bfY5Kt03HSPB8YlPe\nkqxm4frI8bw7F1xXDAou9xrMa5uHE2W1bK3NRZ1Nz36wwyM4hCXOA7eC+WRaYSWq7AfbHz9OaXX3\nbNGiRejs7LTKs7zZ9EzPCP7irQ99YjbJqSebhfKb10126bGNegpjqwKb5gnU5sRQZC9tMH1XXnll\n6O1Ym+EPrieerG+Px5FFgO+yLFqSpSXjXhwk7ZVxPRGaxs3HuUvUKvCt8mN4fTbriUflXVp1zDT3\nw7TcWJI7L3HaBpteQNc1oZPM63F95HjY0Ku4F9e2XC4cTb2mNg8n8vPLZtlalzpqGl7hs+2QCnZ4\nJH2wS3PZ9deaj3Nxbnv3rKury3kITV6MwbSILBGR1wEcLyI7vH+vA9gK4Pu5pbAmWi2Un2Td5Dhr\n9TYX9jSGNtiko8heWpuJHjbHocihKmnJIsC3WTc97pjCJNv2JemViXsiDBs3n3Q+RlhdrfJjeMOe\ncGrq4TddBJvGJMdhmvvR/LjvsHXN0xqG0yqPzjzzTOu1i23KW5I7hqa2xGVy5Zlnnpn6HaMk5zRT\nJ1SrhxMsGSwcAAATX0lEQVT5wwRtHjTmUkej7kQD9hdvNhP/XO4Su5bdJMchi/NGmozBtKr+laoe\nDOCbqjrS+3ewqo5V1SU5prFWWq2bHJT2qgMmpsqf1soWQUX30raa6FH2oSppSvtJYEUeW9eHI7j2\nyiQNVoPlyh9jbTuuNiiPSZ55cn2okOkiOM0yFrX2d1gZO/bYYzMt92F59PDDD+930RHkWt6yWPXG\ndnJlZ2cnHn744VTvqqR5p8b1TpipbUn6oDGbO9Gt2Ez8S/MucR3umNmyGebxSxE5xP9FREaJyDkZ\npqnWTJXQ1OuRR49t0ts9rvLupU17fHMeQ1Wa09zT0+P8mTgNVtLGr8ge+FbbTtIrk/Zj0ydOnIhd\nu3Zh6tSpqSwfl+YkzzzZXqS0ughOU1T9NpWxLMt9nNUvXMtb3sO+gqLOP2k+4jvNc1pUfpnalqQd\nZr4knRZRE//Sekx4UJXvmLmyCaavVtXf+794601fnV2S6i3NdZPTksbtnrjSXjmjWRbjm7MeqhKW\n5iuuuGK/tXWT7KdJnRu/JL0yeT02vRWXSZ5VeAxvGXvUi54w3CxOHqW1hnSaw75MwVrU+SeLR3yn\nISq/TG1Lmh1mSS7eTENYTI8JT6KM9TsrNsF02HtsHkNOIdJcNzktedzuCZPHLaCk45vDVn/IejhD\n2JqiQ4cOjVybO60guM6Nn81yTKYTWxoBVhr122WSZ9GTXG0UOSHZpOihaM3i5lEaa0gnYTPxfuvW\nrbEeetRKHuXKlF+mMepFdpi12gfTXK6knShlrN9ZsQmmHxeRvxGRY7x/fwPg11knrM7KNhktzQfU\nuMijFzRJcOi6XKGtVr3xzWnu6upq+YCZ5s/09PRg8eLFzrfs6t74xV3dpiwBlsskz6LbFRtl6wX2\nlSkvw/Io69Uv0mITrMV56FErRZarqGdUlKlcBWXViVLW+p0Fm2B6MYA+APcBuBfAmwC+kGWiKF9p\nPqDGRR69oEmCwyyCfZve+OY02wzDCX7Gn1C6bNmyREsw1b3xK7rnLo6yBPVpKWJ/sh5aljbTEmFV\nOuatJiPaPPQIyHbJwTSVoa1wkVUnStHHIU8tg2lV3amqXwVwsqq+X1W/pqo7c0gb5aiIyp9HL2iS\n4DCLYN8mQG9O8+TJk50eh5tkQmk7NX5VVbUTdSt57k9VVxeo0hJhYVq19a0eepTHkoNVuLDKSpad\nKHVrr0xsnoD4YRFZh8YjvyEi7xGRb2WeMqq9PHpBkwSHWQT7NgF6c5p3795tXFs37DPr1q1LZQmm\nujd+VH2uwVCdJ9iWWZy2PusL+7QvrFavXo3Vq1enkra8sRMlOZthHn8L4DQAvQCgqr8B8JEsE0Xt\nIa8KHDc4LHK5qOalwJYuXWr9UJ66j3smAuIFQ0VPsK1ywJVE3LY+7yUH63xh1arssRMlGZtgGqq6\nqemltzJIC+UozjrGWShzBQ47Adx4441Yt25d7O+MG6C7PGClncY91w1vO9uLEwy1utDMM9htt8C6\nbG190RdWVC82wfQmEfkwABWRwSLy5/CGfFA1hfXodHZ28sQdIu2ximn2xpt65gDwll0FVXU8b1Hi\nBEO80CQf7+BRmmyC6cvRWL3jCAAvADgBXM2j0uI8UYvSk1YPTVTPXBl6gdjL6qbdbjsnFScY4thQ\n8vHCitJkfPiKiFyvqlcCmKOqC3JME2WMt7fCVe2Wa5mPY7CXddasWVizZg0WLBjYjFQtv7NW5uNZ\nRn4wtHz58gFlzGbI1C233AKgUQZnz54dOw0sw9U1YcIEXHbZZZg3bx56e3sxffp0XlhRbFE902eI\niABYkldiKB9lvr3VbuMIkyjzcWQvq7syH88yYi8zJVWGO3hVwXNztKhg+hEArwE4XkR2iMjrwf9z\nSh9loMpP1KK3lfk2JXtZ3ZX5eJYVgyEiKgPjMA8AV6nqX4jI91X17NxSVCFVvUoLu73V1dXVf+uT\nqqHMtyn9XtY5c+b0v8Ze1mhlPp5VbeuoNR5bouSigumfA3gfAPZCO6pC49Q8bhAAg+kKSnP8Z5ri\njmdtd2U9nlQvNueoKpzH2kGVjkMwrVVKdxqigul3iMj5AD4sIp9u/qOqfi+7ZNVHuxUoIqDcvaxU\nDmwbKU/tHOhR9qKC6csBLAAwCsCnmv6mABhMtwlTw8MGiaLk1ctaVDlk+c8G85XIDetM8YzBtKqu\nAbBGRB5X1eU5pomIMsSGt7U651Gd942ywTKTvbTymMeqGFHrTH9FVW9Q1eUicq6qfjfwt2+o6tfy\nSWJ7YUUgIiIX7X7eqOr+VzXdtL+oYR7nAbjB+3kJgO8G/nY6AAbTbYiVvxzKdhzKlh6qh3YvV+2+\n/82qmh9VTTfZiwqmxfBz2O9toSwVoizpoPbDspcv5jcRxZVH+8E2qiEqmFbDz2G/E1EINjTU7lgH\niNy0S52p035GBdPv8Z50KAAOCjz1UAAMzTxlVCllqBTNaShDmqg4PP71wuNJRGUVtZrHgXkmhIio\nytol2Euyn1xms77qcAzrsA95YD7tL6pnmkok78LLytIeynKcGWRVU5mPT5nTRhQHy3R5MZgmolLh\nQ1jcZNFTTETtjW2Dm0KCaRE5F8A1AKYB+ICqPl5EOmywQKWHeVluWfcO8/hTlbC8ko9loThVyfui\neqbXAvg0gFsK2j4REVEmqhIAxGGzb3Xef8peFctPIcG0qq4HAJG2XK6aaiJY4atY+auGeUzUGucf\nZCtpPrb7cajr/pd+zLSIXArgUgA46qijCk4NEVE11PWkRUTVqt9VSmtcmQXTIvIjABND/tSpqt+3\n/R5VvRXArQBw0kknVfphMa4Fqh0KYJR2uZ1Yh32g+knSXrFMZ495nC3mb76qnt+ZBdOq+rGsvpvq\npchl/6pegckOjzOlhWUpe1yhhqqm9MM8KLk6NC512AcilmOidLTLncssMF/SV9TSeHMBLAMwHsDD\nIvKkqp5WRFqIyirNBi+LxpMNcr1U6XiWvW60Cwa0RA1FreaxEsDKIrbd7tiwUbtrxzrQjvtMRJQX\nDvMoAZ7o3JQ9v8qePspe1mWAZYyyUMZyy7JOVcBgOiVcqeNtZdy3MqSJEx8pCZYZqjOWb6oyBtM1\nVZaGqSzpICIiooF4jk4Hg+mKY0XIBvM1PuYdlR0n5BJRmhhME5UIT8jJMP+IKC9R7Q3bovbCYJr6\nsfJTVaRVVouc68De0fRltf/tnq9EFO2AohNARERERFRV7JkmIiKitsS7DpQGBtNENcWTRHXwWFEr\nXKOZqLwYTBMRERG1OV58xcdgmoiIiHLFwI3qhME0ERFRSTDIJKoeBtNERBEY3NQDjyMRZYXBNBER\nUUYYxBPVH9eZJiIiIiKKiT3TRC2UpWcpz3SUZZ+JiIjKjsF0jTAAIiof1ktqhWWEqNoYTBMRlQgD\nKyKiamEwTYVgwEBERER1wGC6jTGgrR8eUyIionxxNQ8iIiIiopgYTBMRERERxcRgmoiIiIgoJgbT\nREREREQxMZgmIiIiIoqJwTQRERERUUwMpomIiIiIYmIwTUREREQUE4NpIiIiIqKYCgmmReSbIrJB\nRJ4SkZUiMqqIdBARERERJVFUz/QPAcxU1eMB/BuAJQWlg4iIiIgotkFFbFRVfxD49RcAOopIBxER\nVdvq1auLTgIRtblCgukmFwO4r+hEEFUVgwkiIqLiZBZMi8iPAEwM+VOnqn7fe08ngL0A7on4nksB\nXAoARx11VAYpJSIiIiKKJ7NgWlU/FvV3EbkIwCcBnKKqGvE9twK4FQBOOukk4/uIiIiIiPJWyDAP\nETkdwFcAnKyqu4pIAxERERFRUkWNmb4JwBAAPxQRAPiFql5eUFqIKCMcz01ERHVX1GoexxaxXSIi\nIiKiNPEJiEREREREMTGYJiIiIiKKicE0EREREVFMDKaJiIiIiGJiME1EREREFBODaSIiIiKimCTi\n4YOlIyLbAPyugE2PA/BKAdutMuaZG+aXO+aZG+aXO+aZG+aXO+aZm7zz6xVVPb3VmyoVTBdFRB5X\n1ZOKTkeVMM/cML/cMc/cML/cMc/cML/cMc/clDW/OMyDiIiIiCgmBtNERERERDExmLZza9EJqCDm\nmRvmlzvmmRvmlzvmmRvmlzvmmZtS5hfHTBMRERERxcSeaSIiIiKimBhMtyAip4vIsyLynIh8tej0\nlI2IHCkij4nIOhF5RkT+zHv9GhF5QUSe9P6dUXRay0REnheRp728edx7bYyI/FBE/t37f3TR6SwD\nEZkSKEdPisgOEbmCZWwgEbldRLaKyNrAa8YyJSJLvHbtWRE5rZhUF8eQX98UkQ0i8pSIrBSRUd7r\nR4vIHwJl7ebiUl4cQ54Z6yHLWGh+3RfIq+dF5Env9bYvYxHxROnbMQ7ziCAiBwL4NwCnAtgM4FcA\n5qvqukITViIichiAw1T1CRE5GMCvAZwDYB6AN1T1rwtNYEmJyPMATlLVVwKv3QDgVVVd6l24jVbV\nK4tKYxl5dfIFAH8EYCFYxvqJyEcAvAHgLlWd6b0WWqZEZDqAHgAfAHA4gB8BeLeqvlVQ8nNnyK+P\nA1ilqntF5HoA8PLraAD/7L+vXRny7BqE1EOWsfD8avr7jQB+r6p/yTIWGU9chJK3Y+yZjvYBAM+p\n6kZV7QNwL4CzC05TqajqS6r6hPfz6wDWAzii2FRV1tkAvuP9/B00GhEa6BQAv1XVIh7eVGqq+hMA\nrza9bCpTZwO4V1V3q+p/AHgOjfaubYTll6r+QFX3er/+AsCk3BNWYoYyZsIyFpFfIiJodDr15Jqo\nEouIJ0rfjjGYjnYEgE2B3zeDgaKRd2X9XgD/6r202LtdejuHLOxHAfxIRH4tIpd6rx2qqi95P78M\n4NBiklZq52HgyYdlLJqpTLFta+1iAP878Ptk7/b7j0XkT4pKVEmF1UOWsWh/AmCLqv574DWWMU9T\nPFH6dozBNKVCREYAeADAFaq6A8D/BPAuACcAeAnAjQUmr4xmqeoJAD4B4Ave7cB+2hh/xTFYASLy\nDgBnAfiu9xLLmAOWKXsi0glgL4B7vJdeAnCUV2e/BKBbREYWlb6SYT2MZz4GdgywjHlC4ol+ZW3H\nGExHewHAkYHfJ3mvUYCIDEaj4N+jqt8DAFXdoqpvqeo+ALehzW7vtaKqL3j/bwWwEo382eKNGfPH\njm0tLoWl9AkAT6jqFoBlzJKpTLFtMxCRiwB8EsAC78QN7zZyr/fzrwH8FsC7C0tkiUTUQ5YxAxEZ\nBODTAO7zX2MZawiLJ1CBdozBdLRfAThORCZ7vWLnAXio4DSVijfuazmA9ar6N4HXDwu8bS6Atc2f\nbVciMtybXAERGQ7g42jkz0MAPuu97bMAvl9MCktrQE8Oy5gVU5l6CMB5IjJERCYDOA7ALwtIX6mI\nyOkAvgLgLFXdFXh9vDf5FSLyLjTya2MxqSyXiHrIMmb2MQAbVHWz/wLLmDmeQAXasUFFbLQqvBnd\nXwTwKIADAdyuqs8UnKyy+WMAFwB42l/iB8DXAMwXkRPQuB3zPIDLikleKR0KYGWj3cAgAN2q+oiI\n/ArAChFZBOB3aExOIfRfdJyKgeXoBpaxt4lID4DZAMaJyGYAVwNYipAyparPiMgKAOvQGM7whXZa\nZQEw5tcSAEMA/NCrn79Q1csBfATAX4rIHgD7AFyuqrYT8WrDkGezw+ohy1h4fqnqcuw/9wNgGQPM\n8UTp2zEujUdEREREFBOHeRARERERxcRgmoiIiIgoJgbTREREREQxMZgmIiIiIoqJwTQRERERUUwM\npomICiQih4pIt4hs9B4v/3MRmVt0uoiIyA6DaSKigngPKXgQwE9U9V2qeiIaa9BOKjZlRERki8E0\nEVFxPgqgT1Vv9l9Q1d+p6jIAEJEDReSbIvIrEXlKRC7zXp8tIqtF5H4R2SAi93iBOUTkRBH5sdfL\n/WjTE+rgvedOEfkHEfmZ1yPeEfjefw687ybv8doQkedF5K9E5EkReVxE3ud9/29F5PIsM4mIqMwY\nTBMRFWcGgCci/r4IwO9V9f0A3g/gEu+xuQDwXgBXAJgO4F0A/lhEBgNYBqDD6+W+HUCX4bsPAzAL\nwCfReMKYjf9U1RMA/BTAnQA6AHwQwLWWnyciqh0+TpyIqCRE5B/RCHD7vAD64wCO93uOARwC4DgA\nfQB+qaqbvc89CeBoANsBzMTbj8M+EMBLhs09qKr7AKwTkUMtk/iQ9//TAEao6usAXheR3SIySlW3\n2+8tEVE9MJgmIirOMwA+4/+iql8QkXEAHvdeEgCLVfXR4IdEZDaA3YGX3kKjPRcAz6jqhyy2Hfy8\neP/vxcA7lkMNn9nX9Pl94PmEiNoUh3kQERVnFYChIvK5wGvDAj8/CuBz3vANiMi7RWR4xPc9C2C8\niHzIe/9gEZnhkJ7fAZguIkNEZBSAUxw+S0TUltiTQERUEFVVETkHwN+KyFcAbAOwE8CV3lu+jcbw\njSe8CYbbAJwT8X193pCQfxCRQ9Bo4/8OjR5wm/RsEpEVANYC+A8A/y/WjhERtRFR1aLTQERERERU\nSRzmQUREREQUE4NpIiIiIqKYGEwTEREREcXEYJqIiIiIKCYG00REREREMTGYJiIiIiKKicE0ERER\nEVFMDKaJiIiIiGL6/9ZzPQF0+J/zAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1187a18d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(12, 3)\n",
    "plt.errorbar(np.arange(beta.shape[1]), beta.mean(0), yerr=2*beta.std(0), fmt='o', color='k', markerfacecolor='w');\n",
    "plt.axhline(0, lw=1, ls='--', color='lightgrey', zorder=-1)\n",
    "plt.xlabel('Gene num')\n",
    "plt.ylabel('Effect size');\n",
    "plt.title('Reverse DE Results')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are not exactly stellar. All effect sizes are quite small, and very uncertain! The P-values illustrate this as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pval = np.sum((beta < 0), 0) / beta.shape[0]\n",
    "pval = np.minimum(pval, 1 - pval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10dea0a20>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAADgCAYAAABYS05fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX94HGW59793Sg+BbpNNzM83a7t1i7RNoVhiLYeCrQJF\nAVsbfpyiB4WCh0J6LOr7ImdfrR5PLkTPOVoQgWIRQrupyo/KEaSIiS/0IFwURSiNCLZwWsqPUGmB\nQmib3O8fO7PMbmZmZ5PZ3dnN93Ndc2V35pmZe57Mfud+7ud57hFVBSGEEP+oKLYBhBBSblBYCSHE\nZyishBDiMxRWQgjxGQorIYT4DIWVEEJ8hsJKSAARkS+KyOZi20FGBoWVpBCRF0TkXRF5W0ReEZFb\nRSRUbLtGiohERUSN63lbRF4VkV+JyKkZ5azXbS4/cjjmt0TkoFFmr4g8IiInFOaKSKlAYSWZnKWq\nIQDHAfgIgKsKeXJJ4vd9GTauaRaA3wC4W0S+mFHmLFUNWZYOl+P9zDhePYDNAO4SEfHZZlLCUFiJ\nLar6CoBNSAosAEBEDheRfxeR/zG8vxtF5AhjW5+InGkpe5iI9IvIbOP7XMO72ysifxKR+ZayvxOR\nThH5bwDvAPiQ0RTeLiJvicgOEfmcpfxFxvneEJFNIjLZ6zWp6moA3wJwzWgFXFUPArgNQBOAD2Ru\nF5EbROTfM9b9UkS+Ynz+uoj81bjGbSLyWbvzWDzvwyzrficiF1u+29aJ8aD6gYi8JiJvisjTIjJz\nNNdNskNhJbaISATApwA8b1n9XQAfRlJspwJoAfBNY1s3gKWWsgsBvK6qfxCRFgD3Avg3ALUAvgbg\nThGpt5T/RwBfAjARQD+AawF8SlUnAvh7AE8adi0C8C8AliDpMT5snDsX7gLQAODoHPdLQ0QOB/BF\nADtV9XWbIt0AzjO9WRGpAXAagA3G9r8COAlANYBvA1gnIs0jsMOtTk4DcDKS/7dqAOcC2JPrOUiO\nqCoXLlBVAHgBwNsA3gKgAH6LZDMaAATAfgAxS/kTAOwwPk819jvS+L4ewDeNz1cCuD3jXJsAfMH4\n/DsA/2rZNgHAXgDtAI7I2O/XAJZZvlcg6eVOtrmeqHEdh2WsrzTWn5hx3XstyyUOdfQtAAeMMq8B\n6AFwvENZAfA/AE42vl8CoMel/p8EsMj4/EUAm52uw6izi7PVCYBPAPgLgLkAKop9j42VhR4ryWSx\nJr3E+QCmAagz1tcDOBLAE0Zzfi+A+431UNXnAfQBOEtEjgTwGQAJY9/JAM4x9zP2nQfA6p3tND+o\n6n4A5wG4FMDLInKviEyzHGu15Th/Q1LAWnK4RrPs3zKuO2xZbnbZ/+dGmQZV/YSqPgEAIvKMpfPr\nJE2q3Aa878mfj+QDB0b5C0TkScu1zMT79Z0LjnWiqj0AfgTgegCvicgaEakawTlIDlBYiS2q+v8A\n3ArAjBG+DuBdAK0W8anWZCeOiRkOWARgmyG2QFI0b88Qrgmq+l3rKTPOv0lVT0VSfP8M4GbLsf4p\n41hHqOojOVzeZ5H0Np/NYZ+sqGqrvt/59bCxuhvA2UbM82MA7gQA4/vNADoAfEBVwwC2IimImew3\n/h5pWddk+exaJ6p6raoeD2AGkiGB/+3LBRNHKKzEjR8COFVEZqnqEJJC8AMRaQAAEWkRkYWW8huQ\njOktx/veKgCsQ9KTXSgi40SkUkTmG3HcYYhIo4gsEpEJAN5Dspk+ZGy+EcBVItJqlK0WkXO8XIxx\n3A4AqwBcZVxTXlHVPyL5UPoJgE2qutfYNAHJh0m/YduFSHqsdsfoB/ASgM8b9XcRgJiliGOdiMhH\nReRjIjIeSYEewPt1SfIEhZU4Yvygu/B+B9WVSHZmPSoibwJ4EJYOIFV9GcDvkexs+pll/U4kvdh/\nQVJIdiLpNTndfxUAvgJgN5LN2o8jKdZQ1bsBXANgg2HDViQ72dzYKyL7ATwN4NMAzlHVWzLK/Jek\nj2O9O8sxcyEB4BRYHjaqug3AfyBZX68COAbAf7sc4xIk62wPgFYAKQ89S51UIflAfAPAi8b+3/fj\noogzkgwDEUII8Qt6rIQQ4jMUVkII8RkKKyGE+AyFlRBCfIbCSgghPnNY9iKlx+mnn673339/sc0g\nhJQfnrKYlaXH+vrrdvkwCCGkMJSlsBJCSDGhsBJCiM9QWAkhxGcorIQQ4jMUVkII8ZtiZ9r2cwFw\nFoA1NTU1imRKtrTlpZdeUlXVVatWcTu3czu3j2S7Jy0qy+xWbW1tumXLlmKbQQgpP8buOFZCCCkm\nFFZCCPEZCishhPgMhZUQQnyGwkoIIT5DYSWEEJ+hsBJCiM9QWAkhxGcorIQQ4jMUVkII8ZmyElYR\nOUtE1uzbt6/YphBCxjBlJayq+l+q+qXq6upim0IIGcOUlbASQkgQoLASQojPUFgJIcRnKKyEEOIz\nFFZCCPEZCishhPgMhZUQQnymrISVEwQIIUGgrISVEwQIIUGgrISVEEKCAIWVEEJ8hsJKCCE+Q2El\nhBCfobASQojPUFgJIcRnKKyEEOIzZSWsnCBACAkCZSWsnCBACAkCZSWshBASBCishBDiMxRWQgjx\nGQorIYT4DIWVEEJ8pqyElcOtCCFBoKyElcOtCCFBoKyElRBCgkDghVVEPiQia0XkjmLbQgghXsir\nsIrILSLymohszVh/uog8KyLPi8jX3Y6hqttVdVk+7SSEED85LM/HvxXAjwB0mStEZByA6wGcCmAX\ngMdF5B4A4wBcnbH/Rar6Wp5tJIQQX8mrsKrqQyISzVg9B8DzqrodAERkA4BFqno1gDPzaQ8hhBSC\nYsRYWwDstHzfZayzRUQ+ICI3AviIiFzlUu5LIrJFRLb09/f7Zy0hhORIvkMBo0ZV9wC41EO5NQDW\nAEBbW5vm2y5CCHGiGB7rSwA+aPkeMdYRQkhZUAxhfRzAUSIyRUT+DsA/ALinCHYQQkheyPdwq24A\nvwdwtIjsEpFlqnoIQAeATQD6APxcVZ/Jpx2EEFJI8j0qYKnD+vsA3JfPcxNCSLEI/MwrQggpNSis\nhBDiMxRWQgjxGQorIYT4DIWVEEJ8hsJKCCE+Q2ElhBCfobASQojPeBZWEZknIhcan+tFZEr+zCKE\nkNLFk7CKyCoAVwIw0/aNB7AuX0YRQkgp49Vj/SyAzwDYDwCquhvAxHwZNRbp7u7GzJkzMW7cOMyc\nORPd3d3FNokQMkK85go4oKoqIgoAIjIhjzaNObq7uxGPx7F27VrMmzcPmzdvxrJlydd8LV1qm26B\nEBJgvHqsPxeRmwCEReQSAA8CuDl/Zo0tOjs7sXbtWixYsADjx4/HggULsHbtWnR2dhbbNELICBBV\nb8n2ReRUAKcBEACbVPU3+TRsNLS1temWLVuKbYZnxo0bh4GBAYwfPz617uDBg6isrMTg4GARLSOE\nZCBeCnlOG2gIaWDFtJSZPn06Nm/ejAULFqTWbd68GdOnTy+iVYSQkeJ1VMBbIvKmsQyIyKCIvJlv\n48YK8Xgcy5YtQ29vLw4ePIje3l4sW7YM8Xi82KYRQkaAJ49VVVMjAEREACwCMDdfRo01zA6qFStW\noK+vD9OnT0dnZyc7rggpUXKeeaVJNgJYmAd7xixLly7F1q1bMTg4iK1bt1JUCfGRQg9n9OSxisgS\ny9cKAG0ABvJiESGE+EgxhjN69VjPsiwLAbyFZDiA5BlOHCBkdBRjOKPXGOuFebOAOMKJA4SMnr6+\nPsybNy9t3bx589DX15e3c7oKq4hcB8BxoKuq/rPvFpEU1ictgNSTtr29HQDFlRAvFGM4YzaPtXRG\n2ZchTk/affv2pYZiUVwJcccczpjZ8svrzEZVLbvl+OOP13KgtbVVe3p60tb19PSk1re2thbJMkJK\ni0Qioa2trVpRUaGtra2aSCRGeihPGuRpSquI1COZNnAGgEqLKH8if5I/ckptSqsTTjHWzs5OnH32\n2ZzySkjh8TSl1euogPUA+gBMAfBtAC8AeHxEZhHPLF26FJ2dnWhvb0dlZSVWrFiRmjjgJUbEEQWE\nFAkvbi2AJ4y/T1nWPe7VLS70Ui6hAJNEIqFTpkzRnp4ePXDggPb09OiUKVNcmzMj2YcQkhVvmump\nEPCo8XcTgDMAfATAX72epNBLuQmrau4xIrf4LCFkxPgaYz0TwMMAPgjgOgBVAL6tqvf47kL7QLnE\nWEcDUxESkhd8jbE+pqr7VHWrqi5Q1eODKqokiTl2zwpTEZJSpCT7Cry4tQD+AuABAMsA1Hh1h4u1\nFCoU4OMQDt9hjJWUAwG8j/2LsWpSXOcA+E8A2wH8CsDnve5b6KUQwhrAf/gwgiz8hHghgH0F/sVY\nrYhInSGwn1PVcf75zv6Rzxhrd3c3Ojs78eKLL+Kee+5JmybX29uLFStWYOvWrXk5NyFjjQD2FfgX\nYxWRKhH5goj8GsAjAF5G0oMdU5gD9q+77jq88847BU/sQEi5kS1+WrJ9BV7cWgA7APwAwAleXeFi\nLvkKBVibJQFsohBSUngJpwUw5OZfjBXG21yNz01eD16sJV/CWlFRoQcOHFDVQP7DCSkpvDonAesr\n8E9Y03YA/pDrPoVeCuGxqib/4bFYTEUkCP/wvBGwG5uUCVZHxeTAgQNaUVFRJIs84UmDcn7nFTwG\nb8uRzLepNjU1YWhoCOvXry/b91RZ48oDAwO47rrrEI/HS2MsIQk0JRs/9YJXBTYXAJfluk+hFiRf\nHbNm6tSpPj2chlMu3pvX62AsmeSLEg2njT4UgGSKwJUAfgTgnwAc5vXAxVzKMVeAn+RyQ5doc40U\niVwdjxJ0VHwR1p8BWGeI6kYAq70euJgLhdWdXLxQeqzEKyXqgeaKL8L6tOXzYSiBjistY2Ed6dM9\nc79cvNAx8mMhPpDtIVyC3qkdvgjrH9y+B3UpR2EdqcDZ7RcOhz17oYlEQiORiIZCIRURjUajpfqD\nGNMUQtTcHth292FjY6NGIpFSu598EdZBAG8ay1sADlk+v+n1JIVeylFYR9okt9svHo9rc3NzVpGm\nt1oeFOr/6HaPOm2LxWKldk+NXlhLdSlHYR1pJ5LTfubYWzcPhvHV8qBQ/0c3AXe7f0vsnqKwlhN+\neqxe414cEVCajCam7ve5Ozo6tLW1VUOhkON9WGL3FIW1nPAzxmrul+2Y9FhLj9HG1PNlS1dXl06a\nNMn2Xiuxe4rCWm74NSrAq3DmMzZXJj3EgWM0MfV82BKPx1P/50gkog0NDSoiGovFtKurqxTj9hRW\nvyk3MfDSROzo6NCamhoVEa2pqdGOjo5Rn5edYvljNDF1vxER2/9zMWzxEQqrn5SjGBTLY2WIIX8E\nqW5rampsbampqSm4LT5CYfWTIN2wflGIGKudl89OMX+x1nE0GtWGhoZAOABl+n+msPpJmd4keR0V\n4CTckUik7B5SxcKsY2ssMxwOa21tbdGb2uXojCiF1V/K9CZxZbT5Z53qLBqNll1YpViYHUSZ9dnc\n3DwsE7/Vq41EInkX3nIMnymF1V/K9CZxxetwGSeyTXEs4Q6MwGDWX66x8kmTJhWkV74M/88UVr8p\nlZvETzs7Ojo0HA47DvCuqalxzeVqbaJavxN/8DIBIJvwem155XJflcpvZQRQWMcK1ps4Eonk7Fm6\nYf4o3X68Tsfv6OjQpqamNFuampp8GbJFkiQSiawTALIJr5e4eS4ttjJv3VFYxwKZN3EsFvM1Fmz+\nKN28Hqfjj8W4dDHo6OhwnQAQiUQ0FouleY+ZHmskEnH1MJnDNwWFdSyQeRP7PXrBPL6bF+J0/HId\nSRFEnJreiURiWAtm0qRJWldXl4qxNjQ02LZyzHn+FRUVGgqFtKurK+2cY/T/TmHNhVKNCWXexH57\nC5kdWOaoAGteVrvjJxIJra6uLmfPJfAkEgnHQfrV1dWpez0ajXqaBjtp0qS03wU9VgqrK6UcE7Ib\nEpUtxuqUgcjpoWKWNwU1cwB6Y2Nj2j4dHR1aXV2tIjIsxpr54yT5wbynvXiPdmWcxDEWizHGSmH1\nRik/Ye1u4oaGBo1Go7ZCaVe+qalJ4/F41h+BNSxgCnEsFtNIJJJ2fKunE4/HNRwO25Yl+cMa+852\nb9uVccs5wFEBFFZPlHpMKJeb2MsPzemh4sUTcTt+KdVpqWPe017+Z0FKNVgCUFi9Usoea654eYi4\nCWA2ER9ppvgy9nDyhludWe9pazmncceZM7OqqqpsO734fxmDwgrgLABrpk6dmlNN5TsmFCTRcJum\naqYHHE2T3S025zSGtcxjcr6RS7KV0dSpl5DPGGbsCau5BGlUQNBEw8s01Xg87uqduNWV3fU2NjZq\nfX19ak77aMZIjlWcpqW69dKP9J4u9dBYnqGwBoEgiob5g3N7D5HbjzQcDrt2dllHEcRisbRt8Xhc\na2pqmEZQ/YmNZw7sF5FR25VtQsEYh8IaBIIsGm62WW3M5nU7/egyj+90HLtxlOX+Q861JWP3v+rq\n6tLGxsa0Y2RmtRqJXXax1YaGBoZmklBYg0AQPVYTr9NUs12D04Micz83rytI4ZJCkOt9YVfe7+nL\nbnZFo9ERH7PMoLAGgaDFWL3YlhkLzeZ1ex2eNZI0gkHq+POTbHVqN4nDOj548eLFWlVVNexdZHYP\nOWtoJjMMk6tddsctt/9NFiisQSHIN2DmUBy7JNZunqbdFFe343vxsKz7mOLR2trq2PlVirh5rE4P\nvIULF2o4HLad0WaOuDDHoFoF2fqwzDaRxKsnHWSHIc9QWIk/OA0gz4zvZfthjXSwulUYyiWfq1td\nZBPdqqoq2+3hcHjYLLrm5uZUnWW+itp8RU5m4hUvghnkEFeeobAS/8j0ukfa4WQ9TiQSycljMv8G\noeNvNGRrmrtNJ00kEioijtvj8XjaerPOMl9F7Raf9dLCCnKnbJ6hsJL8ka8XDXoZqRAKhVxDKkEP\nvWSb0ZT50LJO4qiurk7zWDPDJk6p/TLDMKP9/9FjpbCSPDDaH5bT/m5xWLO529XVpfF4PJVByxrf\nDXrsz8nTN3vdE4lEWn5Uu0kcEydOTDX5vUwaMGOuXV1dafFuJ+/WC0Gv5zxCYSX5YzQ/rEQioaFQ\naJhHaTZn7d44unjx4lQHjdsspKB7Um7NeNXh00mdJnHU1tY6xlpbWlrSJnMsXLhQjzjiiGEC7TWr\nmRNBbhnkEQoryS8j+WG5CbIpgOYLDK0vIGxsbNTly5erqru3bI4gMOfTm695DoVCw+wrhjA4CWUo\nFFLV4U10tya7Wyw2s36dslVZE16PEWEcLRRWEjycRDEcDqdm97glclF179yxy3pv95rnQjdlrR1W\nVjEz8zKYCU7sJlXYvek2FotlrSfrulAoNFY7m/yGwkqCh1UUrR5jKBTSurq6tNwBmR6l2evt5H05\nNY3NMID1s18hAy9er9UDz8yz0NTUpFVVVY6Cb4ZAMpvwhx9+uO2DpLGx0bYDy23EAMkJCisJHtYY\nol2c1ByCZdcx09TUpNXV1SkvL3Nfp/il6ZVZP/sxXMjJ67W+6qa2ttZ1vK8ZQ86cZZVtUoXplWY+\nnMyxqXblzQfSGOts8hsKKwkephg5jaM0xcEUgcycoKZ42mVgyjZ33m+P1e4Y1pfwdXV1ZX2hYubr\nTuLxeNqMqGwPC+u6UCiUikdnCmhXV1fqHIypjgoKKykO2ZrH2Qa5mx0zXV1dtrld4/G452xP+Yyx\n2oUsrB6mKY7Z8iyYL+izDiEzp646hT2scVTrONdwOKzV1dVM+5c/KKyk8HgVLLfxnK2trRqLxRw9\n0HA4bLstHo/rhAkTNBQKqYhoQ0ODTpgwwVHgOzo6UvkRrIlMvGLNX2Ber1VIzREKbuNyzWFiTvVm\nF2OdNGmS1tXVOY5zNUMmbPbnhbEnrCN9NQvxj1ySeDjNQDIHybs1g4888sg07zQejw8ToClTpjgm\nfrYTsubmZtskNE4kEolhHqX1+s2muZ3XLSKpsEe2BDXxeFyrqqrSHhBdXV0aCoUch2+ZDyg2+33H\nkxaJJgWprGhra9MtW7akrTt48CB27dqFgYGBIlkVDCorKxGJRDB+/Pi8HH/cuHEYGBhIO/7BgwdR\nWVmJwcHBtLLd3d3o7OxEX18fpk+fjng8jqVLl6a2LV++HHfffTcWLFiQ2qe3txdLlizBtddei2ef\nfRbXX3899u3bh4kTJ2Ljxo3Dyra3t+Nvf/vbsPNVV1fjzjvvHFb+oosuwi233IJly5ahs7MzZY/X\n6+3u7sbKlSsxceJEbN++HY2NjVi2bBk2btyIvr4+TJw4EYODg/jyl7+MtWvXIpFIYNeuXbjsssvw\nzjvvpNWDWW8DAwM4/PDDMTQ0lGbrJZdcgh07dniub+IL4qmUVwUupcUuFLB9+3bt7+/XoaGhkTyl\nyoKhoSHt7+/X7du35+0cfs58cpthZT226f05zeSyps+za7KbWBOdeLXZbl6/NWvU4sWLHcetWkdH\nWMeqhsPhVApAsz4zm/bmsKqgzzQrQ8ZeKMBc7IR127ZtY1pUTYaGhnTbtm15O77fA++tY0BDodCw\nMZrWDquurq60/LCmWJnN/Hg8njZQv6WlJdUUzxx47zW5c21tbVpIIzP26yR85oB9uzitGSdtaGhI\nCXFHR4dGo9G0t+g6DVtjPDWvUFit5FNMSo1814WfU0WtwmQnUqaQ2QlMU1OT1tbWpjxQu1fAWEcN\nNDU16fLly7WiosIxLmwXL12+fHnqejPjwtmmnZp15CS+1vqzlrPakvlAoajmFQqrlVzF5N1339Un\nn3xSd+/endN+brzyyiu6dOlSnTJlis6ePVvnzp2rd911l2/H90opPWQyZ2plCpspZG7TO81xnG5v\npTU/h0KhtOm1VtzOYfWErW+mddonGo2mxvN6naxg93JGa34ECmpBoLBasROTl156STdu3KiPPfZY\nWpjghhtu0Lq6Om1tbdWamho9++yzde/evVnq252hoSGdO3eu3nDDDal1L7zwgl577bWjOu5IKCVh\nzRQm65jN1tZWbWlpSXl+TuJkip/XmVlm2ebm5lSowGyGW1PvmbPE7F6V0tjYqLW1ta4zr8x4bLZJ\nBE514VSO5BUKqxWrmAwNDemVV16p4XBYzzjjDP3whz+sbW1tumvXLr3vvvs0Go1qX1+fqqru379f\nL774Yj3nnHNSxzn//PM1Go3qCSecoOvWrfMUu33wwQf15JNPtt126NAh/drXvqZtbW16zDHH6I03\n3qiqqr29vfrxj39c29vb9eijj9bzzz8/da4tW7boySefrLNnz9bTTjst5VmvXr1ap0+frsccc4ye\nd955tucLirB6CRk4dWCZzfe6ujptaWlxnXVlNr3r6+tdB87bdRJNmjRJly9frrFYTCsrK7Wuri51\nDHPev1timcymektLy7DOt9ra2mFz/u3ipIynBgIKqxWrmKxfv15nzZqlr7/+uqomhXbVqlX6yU9+\nUs8880y97bbb0vZ9++23taamRn//+99rQ0ODXnPNNfrcc8/pvffeqzNmzNDvfe972f4Zunr1al25\ncqXttptuukm/853vqKrqwMCAHn/88bp9+3bt7e3Vqqoq3blzpw4ODurcuXP14Ycf1gMHDugJJ5yg\nr732mqqqbtiwQS+88EJVVW1ubtaBgQFVVX3jjTdszxcEYc1FJKwCHI1Gtba2NjUJIBqNakdHh9bX\n1zt6hj09PVpfX287brauri4txmomKrGes7q6WhcvXqxHHnnksGNke+uBFTvhtnqvXuLSYzQHapCg\nsFqxiskpp5yid9xxR9r29957T+vq6vS4447TRx55ZNj+06ZN03PPPVdXrVqVtn779u1aW1ur+/fv\nH7aPlUxhveyyy/TYY4/VtrY2bW9v16OOOkpnzZqls2bN0mg0qps2bdLe3l495ZRTUvtceumlevvt\nt+vTTz+tEydOTJWfOXOmnnrqqaqqunDhQm1vb9fbb79d33rrLVtbgiCsI23WuglyIpHQlpaWVJzT\nOpXVOkrAOgrAHAplDRXYDYGqrq52jNF6ffusdfgXhbFkobBasYrJnDlz9KGHHhpWZurUqXrBBRfo\nV7/61bT1W7du1YaGBp09e7Y++uijw/abMWOG/ulPfxq23opdKKC/v18nT56sS5Ys0fvvv3/YPr29\nvXrGGWekvl9++eX605/+VJ966imdO3eu7XkOHTqkPT09esUVV+i0adP04MGDw8oEQVhHml0q2xRR\ncyRAZnPd6c0E1umsjY2NaU32hoaGtJldbvkN7I5tl6avpqZmFLVGAgCF1YpVTL7xjW+kms4mDz30\nkE6aNElfeOEFjUQi+pWvfEU3b96st9xyi06ePFlvvvlmbW9v15tuuiltv71792o4HNb+/n6nf4Sq\nJsMNc+bM0R//+MepdS+++KJOnjxZb7rpJl20aFHqR/vss8/q22+/7Sis7733nsZisZRnfeDAAd26\ndasODg7qjh07Uuuam5ttwwFBENaReqxuw5fM9zjZ5Xx1y9VqFWSnyQjWjrLMY5jesJlAJRaL6fLl\ny23n8Oeaj4AEDgqrFauY7NmzR1tbW3XJkiW6bt06/eY3v6n19fX6y1/+UlVVd+7cqVdccYXOmTNH\nFy9erA888ICqJsW3ublZH3roIR0aGtJXX31VlyxZohdddJGn/8ju3bv1vPPO02g0qh/96Ed1/vz5\numHDBh0cHNSrrrpKZ86cqa2trTp//nzdu3evo7Cqqv7xj3/Uk046SY899lidMWOGrlmzRg8cOKAn\nnnhi6jhXX321rR1BENaRdsS4CbIpqHZlvGaZsjuuuX8kErGN01ZWVmosFtMJEyakhQ8ikUgqNDGS\nJC8kkFBYrWSKyZtvvqmrV6/Wc889V1euXKnPPPOMp1r9xS9+obFYTOvq6rSqqkovv/xyfffddz3t\nGxSCIKyq/r8zyy2JtlP6PWteVLscp+Y6007rcKvW1tZUnlMzsQpff1L2jD1hhUt2Kz/FZHBwUHfv\n3p21wyqoBEVYR4qTILvNRsrMFZDpITt5rJm5Up2GVdXX13Oc6dhg7AmruXBKqzvlXBduXnDmC/0y\n3z2V+epHVVuBAAAHLklEQVQSM1mK2QnW1tbmmBvVLacqe/7LCgqrlXIWk1xhXdgnuc5FlK35VO3K\ncThV2eJJg8ZMPta+vj5MmzYNIt7SKZYrqoo///nPmD59erFNIaQU8SQgFfm2IihUVlZiz549KMcH\niVdUFXv27EFlZWWxTSGkrDms2AYUikgkgl27dqG/v7/YphQV8w0ChJD8MWaEdfz48ZgyZUqxzSCE\njAHGTCiAEEIKBYWVEEJ8hsJKCCE+U5bDrUSkH8CLllV1AF4vkjmjhbYXB9pePIJs/+uqenq2QmUp\nrJmIyBZVbSu2HSOBthcH2l48St1+gKEAQgjxHQorIYT4zFgR1jXFNmAU0PbiQNuLR6nbPzZirIQQ\nUkjGisdKCCEFoyyFVURqReQ3IvKc8bfGoVxYRO4QkT+LSJ+InFBoW21s8mS7UXaciPxRRH5VSBud\n8GK7iHxQRHpFZJuIPCMiXy6GrRZ7TheRZ0XkeRH5us12EZFrje1PicjsYthphwfbP2fY/LSIPCIi\ns4phpx3ZbLeU+6iIHBKRswtp32gpS2EF8HUAv1XVowD81vhux2oA96vqNACzAPQVyD43vNoOAF9G\nMGw28WL7IQBfVdUZAOYCuFxEZhTQxhQiMg7A9QA+BWAGgKU2tnwKwFHG8iUANxTUSAc82r4DwMdV\n9RgA30FAYpcebTfLXQPggcJaOHrKVVgXAbjN+HwbgMWZBUSkGsDJANYCgKoeUNW9BbPQmay2A4CI\nRACcAeAnBbLLC1ltV9WXVfUPxue3kHwwtBTMwnTmAHheVber6gEAG5C8BiuLAHQZSY4fBRAWkeZC\nG2pDVttV9RFVfcP4+iiAoKQ181LvALACwJ0AXiukcX5QrsLaqKovG59fAdBoU2YKgH4APzWa0z8R\nkQkFs9AZL7YDwA8B/B8AQwWxyhtebQcAiEgUwEcAPJZfsxxpAbDT8n0Xhou8lzLFIFe7lgH4dV4t\n8k5W20WkBcBnEZAWQq6UbNpAEXkQQJPNprj1i6qqiNgNfTgMwGwAK1T1MRFZjWTT9Ru+G5vBaG0X\nkTMBvKaqT4jI/PxYaY8P9W4eJ4SkN7JSVd/010piRUQWICms84ptSw78EMCVqjpUim/9KFlhVdVT\nnLaJyKsi0qyqLxvNNrumxC4Au1TV9JbugHs80zd8sP1EAJ8RkU8DqARQJSLrVPXzeTI5hQ+2Q0TG\nIymq61X1rjyZ6oWXAHzQ8j1irMu1TDHwZJeIHItkuOhTqrqnQLZlw4vtbQA2GKJaB+DTInJIVTcW\nxsTRUa6hgHsAfMH4/AUAv8wsoKqvANgpIkcbqz4JYFthzHPFi+1XqWpEVaMA/gFATyFE1QNZbZfk\nL2UtgD5V/c8C2mbH4wCOEpEpIvJ3SNblPRll7gFwgTE6YC6AfZZwRzHJaruITAJwF4B/VNW/FMFG\nJ7LarqpTVDVq3ON3ALisVEQVQHm+pRXAB5DslX4OwIMAao31/wvAfZZyxwHYAuApABsB1JSK7Zby\n8wH8qth2e7UdyeaoGnX+pLF8uog2fxrAXwD8FUDcWHcpgEuNz4JkD/ZfATwNoK3Y9ZyD7T8B8Ial\nnrcU22avtmeUvRXA2cW2OZeFM68IIcRnyjUUQAghRYPCSgghPkNhJYQQn6GwEkKIz1BYCSHEZyis\nJLCIyKCIPGlZvm6sP8nIjPWkiBwhIt83vn9/BOf4lxzL/6uIOE6SIARgomsSYETkbVUN2ay/EcBm\nVV1nfN+H5JjZQb/OQchooMdKSgoRuRjAuQC+IyLrReQeACEAT4jIeSJSLyJ3isjjxnKisV9IRH5q\n5CZ9SkTaReS7AI4wPN/1GecZJyK3ishWY58rjPW3isjZItJm8aSfNvMiiEhMRO4XkSdE5GERmVbQ\nCiKBoGRzBZAxwREi8qTl+9Wq+hMRmYfkbLM7gJTXeZzxOQHgB6q62ZjSuQnAdCST6+zTZG5SiEiN\nqt4pIh3mvhkcB6BFVWca5cPWjaq6xSgDIwRxv7FpDZKzh54TkY8B+DGAT/hQF6SEoLCSIPOug+i5\ncQqAGZaMSFVGJq1TkJyTDgDQ9/OUOrEdwIdE5DoA98Ih2bKInIdklrTTjPP8PYBfWM5/eI72kzKA\nwkrKjQoAc1V1wLoy19RzqvqGJF9lshDJOeznArgo45gzAXwLwMmqOigiFQD2juBhQMoMxlhJufEA\nkpnnAQAiYorcbwBcbllvvo/roJHGMA0RqQNQoap3Avi/SHql1u1hAN0ALlDVfgDQZF7ZHSJyjlFG\nJEDvmSKFg8JKgozZsWQu3/Wwzz8DaDM6qLYh6W0CwL8BqDE6o/4EYIGxfg2ApzI7r5DMaP87I8a7\nDsBVGdsXAZgM4GbTPmP95wAsM87xDOxfOULKHA63IoQQn6HHSgghPkNhJYQQn6GwEkKIz1BYCSHE\nZyishBDiMxRWQgjxGQorIYT4DIWVEEJ85v8DrbWjnwuzpcUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e318400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(5, 3)\n",
    "plt.yscale('log')\n",
    "plt.scatter(beta.mean(0), np.maximum(pval, 1e-13), c='w', edgecolor='k', label='Genes');\n",
    "plt.axhline(0.05, lw=1, ls='--', c='k');\n",
    "plt.gca().invert_yaxis();\n",
    "plt.ylabel('P-value')\n",
    "plt.xlabel('Effect size');\n",
    "plt.legend(scatterpoints=3);\n",
    "plt.title('Reverse DE P-values')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>beta</th>\n",
       "      <th>pval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000022708</th>\n",
       "      <td>-0.420413</td>\n",
       "      <td>0.1355</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000038900</th>\n",
       "      <td>0.501867</td>\n",
       "      <td>0.1515</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000083084</th>\n",
       "      <td>-0.391816</td>\n",
       "      <td>0.1565</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000046330</th>\n",
       "      <td>-0.386397</td>\n",
       "      <td>0.1600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000003429</th>\n",
       "      <td>-0.329270</td>\n",
       "      <td>0.1725</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000037805</th>\n",
       "      <td>0.425457</td>\n",
       "      <td>0.1815</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000060126</th>\n",
       "      <td>-0.498119</td>\n",
       "      <td>0.1845</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000046364</th>\n",
       "      <td>-0.512343</td>\n",
       "      <td>0.1905</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000021754</th>\n",
       "      <td>0.387441</td>\n",
       "      <td>0.1940</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000064341</th>\n",
       "      <td>-0.252449</td>\n",
       "      <td>0.2055</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        beta    pval\n",
       "ENSMUSG00000022708 -0.420413  0.1355\n",
       "ENSMUSG00000038900  0.501867  0.1515\n",
       "ENSMUSG00000083084 -0.391816  0.1565\n",
       "ENSMUSG00000046330 -0.386397  0.1600\n",
       "ENSMUSG00000003429 -0.329270  0.1725\n",
       "ENSMUSG00000037805  0.425457  0.1815\n",
       "ENSMUSG00000060126 -0.498119  0.1845\n",
       "ENSMUSG00000046364 -0.512343  0.1905\n",
       "ENSMUSG00000021754  0.387441  0.1940\n",
       "ENSMUSG00000064341 -0.252449  0.2055"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({'pval': pval, 'beta': beta.mean(0)}, index=idx).sort_values('pval').head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, negative results are also results.\n",
    "\n",
    "Can we improve this somehow? We can think a little about the expected biology. While biology is complex and intricite, and everything interacts with everything, the results of this way of thinking might not be very actionable. What we expect (or rather hope) is that a small number of key genes determine cell type.\n",
    "\n",
    "In the statistical sense, it means our prior expectation on the effect sizes is that most of the time they are 0. Allen Riddell wrote an excellent post about this concept and the \"Horseshoe prior\" [here](https://www.ariddell.org/horseshoe-prior-with-stan.html). Based on the code in the post, we can make a _sparse_ version of the reverse DE in the following way\n",
    "\n",
    "    data {\n",
    "      int<lower=0> N;\n",
    "      int<lower=0> G;\n",
    "      matrix[N, G] X;\n",
    "      real y[N];\n",
    "    }\n",
    "    parameters {\n",
    "      real beta0;\n",
    "      vector[G] beta;\n",
    "      vector<lower=0>[G] lambda;\n",
    "      real<lower=0> tau;\n",
    "      real<lower=0> sigma;\n",
    "    }\n",
    "    model {\n",
    "      lambda ~ cauchy(0, 1);\n",
    "      tau ~ cauchy(0, 1);\n",
    "      for (i in 1:G) {\n",
    "        beta[i] ~ normal(0, lambda[i] * tau);\n",
    "      }\n",
    "\n",
    "      beta0 ~ normal(0, 1.);\n",
    "\n",
    "      y ~ normal(beta0 + X * beta, sigma);\n",
    "    }\n",
    "\n",
    "Again, we perform the sampling and plot the effect sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with open('SparseLR.pkl', 'rb') as fh:\n",
    "    model = pickle.load(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "s = model.sampling(data=data, iter=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "beta = s.extract('beta')['beta']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAADgCAYAAAAwjrGWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8HHV9//HXmxDhpwiBXLgkhFuhniS1VKPVmso5XoGq\noGJKTPmBPRbwp6i/2gr09Pfz1tOiViuirSCxQCUHUxGlFAU1QTy/egHxAiGgiFjuCQcjESUJ5PP7\nY2fjns3O7szuzl7OeT8fj/M4u7OzM5/5zndmPvPd78woIjAzMzMzs2Ls1u0AzMzMzMymMifcZmZm\nZmYFcsJtZmZmZlYgJ9xmZmZmZgVywm1mZmZmViAn3GZmZmZmBXLCbWZmVkXSJZL+rttxmNnU4ITb\nzKYEScsk/ZekX0p6VNL/k/S8bsdVT5LUbZP0qyTmr0p6VrfjaoWkeyT9RtIWSZuTdXKmpN0qxqlc\n7vLfD1OmNyhpRzLOFkl3SnpT55ZoZwz3dXKeZja1OOE2s74naW/gGuACYD9gPvA+YGsB89q9zZP8\nUETsRSnm+4FVbZ5+QwUs06sj4pnAIcB5wNnsulwfioi9Kv5+v870HkjKaG/gfwOflvS7bY7ZzKww\nTrjNbCo4CiAixiLiqYj4TURcHxE/ApB0WtLi/YmkBfwOSS8tf1nSmyRtSFpQ75Z0RsVng5Luk3S2\npIeAf5U0R9I1SQvuo5K+WW7BlXSQpCslbZL0M0lvz7IAEfEbYA1wdOVwSX+exPYLSddJOiQZ/i+S\n/rFq3C9J+stGcUh6r6TPS/qspMeA0yQ9X9LNkh6T9LCkj1aM/4KkpXqzpB9KGsy4TL+MiKuBPwVO\nlbQky/fqTC8i4lrgUeDZFfE9K/l14NGkBXx5xWfHS7o9Wbf3S/qrZPhpksYrpy8pJP1O1bBnAF8G\nDqpojT+oXnmZmVVzwm1mU8GPgackXSrpOEn71hjnD4GfAnOA9wBfkLRf8tlG4FWUWlDfBPyTpOdU\nfPcASi3nhwCnA+8C7gPmAvsDfwNEknT/B/BDSi3WLwXeKemVjRYgSexWAHdVDDshmfbrknl9ExhL\nPh4D/lSSknH3BV4BXJExjhOAzwOzgMuB84HzI2Jv4AhKyT+S5gP/CfxdUgZ/BVwpaW6jZSqLiO9S\nKq8/zvqdWiTtJuk1lNbhXcmwZwBfBVYD84CTgX+WtCj52irgjKTFfQmwNs88I+Jx4DiSVvbk7wFS\nysvMrBYn3GbW9yLiMWAZEMCngU2Srpa0f8VoG4GPRcT2iPgccCfwJ8n3/zMifpq0oH4DuJ7JyeEO\n4D0RsTVpid4OHAgckkzvmxERwPOAuRHx/ojYFhF3J/GcXCf8v5K0GdiSLMMpFZ+dCfxDRGyIiCeB\nvweOTlq5v5ksbznOk4BvJclglji+FRFfjIgdFcv0O5LmRMSvIuLbyXh/BlwbEdcm434VuBk4vs4y\n1fIApYR90nJX/F1a57sHJWX0G+Aq4C8j4vvJZ68C7omIf42IJ5PhVwJvSD7fDiyStHdE/CIibskZ\nd5q08jIz24UTbjObEpKk9LSIWECpJfMg4GMVo9yfJMVlP0/GIWkV/3bSJWEzpWRyTsW4myLiiYr3\nH6bUwnp90gXlnGT4ISTJYfmPUgt1ZeJf7R8jYhZwKKWEsrJv8iHA+RXTehQQMD9ZlisotYoDvJFS\nS3XWOO6timOYUtecOyTdJOlVFdN6Q9W0llE64chjfhL/pOWu+Du1zncfSMpob+DjwEsqPjsE+MOq\n+FZS+lUC4PWU1ufPJX1D0gtzxp0mrbzMzHbR7gtlzMy6LiLukHQJcEbF4PmSVJF0LwSulrQHpRbR\n/wl8KSK2S/oipcR25ySrpr+FUreSdyX9ktdKuolSEvuziDiyiZj/W9I7gEslXZO0Ot8LjEbE5Slf\nG6OU9J9HqcvMa5PhWeKoXqafACuS7iivAz4vaXYyrX+LiL/Iu0xlKt0tZj4w3mjceiJiq6SzgTsl\nnRgRX0zi+0ZEvDzlOzcBJ0iaCbyNUtePg4HHgadXxHhAre+XJ1NjujXLK+mCYmY2iVu4zazvJRfN\nvUvSguT9wZRafit/5p8HvF3STElvAAaAa4GnAXsAm4AnJR1HqS90vfm9StLvJP2nfwk8RanbyXeB\nLSpdYPk/JM2QtEQZb0+YdNd4gFI/cYBPAedKWpzMd58k9vL43wceAS4GrouIzclHueOQ9GeS5kbE\nDqA8nR3AZ4FXS3plMp09VbqQdEGj5ZG0d9LyewXw2Yi4NUs51BMR24CPAP83GXQNcJSkU5J1O1PS\n8yQNSHqapJWS9omI7cBjyTJBqX/7YklHS9oTeG+d2T4MzJa0T8WypZWXmdkunHCb2VSwhVIL73ck\nPU4p0b6NUit02XeAIyklqKPASRExkbRWv51Sy+cvKHXNuLrB/I4Evgb8CvgW8M8RsS4inqLUp/ho\n4Gf8NhneJ21CNXwYeLekPSLiKuCDlC6EfCxZpuOqxl8NvCz5D0CTcRwLrJf0K0oXBJ6c3O3lXkoX\nWP4NpZOSe4G/pv7x4z8kbUnGHQE+Suli1Erv1uT7cD9SZ3rVPgMslPTqZP29glL/9AeAhyiV2R7J\nuKcA9yTldyal7iZExI+B91Najz+hTut7RNxB6deEu5NuKweRUl45lsHMphFN7tJoZjb1SDoNeHNE\nLOt2LGZmNv24hdvMzMzMrEBOuM3MzMzMCuQuJWZmZmZmBXILt5mZmZlZgZxwm5mZmZkVaEo++ObY\nY4+Nr3zlK90Ow8zMzMymNjUeZYq2cD/ySJ7buZqZmZmZFWdKJtxmZmZmZr3CCbeZmZmZWYGccJuZ\nmZmZFcgJt5mZmZlZgZxwm5mZmZkVyAl3mwwODjI4ONjtMMzMzMysxzjhNjMzMzMrkBNuMzMzM7MC\nOeE2MzMzMyuQE24zMzMzswJ1NeGWdKykOyXdJemcOuM9T9KTkk7qZHxmZmZmZq3qWsItaQbwSeA4\nYBGwQtKilPE+CFzf2QjNzMzMzFrXzRbu5wN3RcTdEbENuAI4ocZ4ZwFXAhs7GZyZmZmZWTt0M+Ge\nD9xb8f6+ZNhOkuYDrwX+pYNxmZmZmZm1Ta9fNPkx4OyI2NFoREmnS7pZ0s2bNm3qQGhmZmZmZo3t\n3sV53w8cXPF+QTKs0lLgCkkAc4DjJT0ZEV+snlhEXARcBLB06dIoJGIzMzMzs5y6mXDfBBwp6TBK\nifbJwBsrR4iIw8qvJV0CXFMr2TYzMzMz61VdS7gj4klJbwOuA2YAn4mI9ZLOTD7/VLdiMzMzMzNr\nl262cBMR1wLXVg2rmWhHxGmdiMnMzMzMrJ16/aJJMzMzM7O+5oTbzMzMzKxATrjNzMzMzArkhLsA\ng4ODDA4OdjsMMzMzM+sBTrjNzMzMzArkhNvMzMzMrEBOuM3MzMzMCuSE28zMzMysQE64zczMzMwK\n5ITbzMzMzKxATrjNzMzMzArkhNvM2s73ojczM/stJ9xmZmZmZgVywm1mZmZmViAn3Nb33H3BzMzM\nepkTbjMzMzOzAu3e7QCmgrGxMdavX8/ExARLlixhx44dzJs3r9thmZmZmVkPcAt3i8bGxhgZGWHN\nmjVs3bqVCy64gM2bN7Nx48Zuh2ZmZmZmPcAJd4tGR0dZtWoVQ0NDzJw5k6GhIS6//HI2bdrU7dDM\nzMzMrAc44W7Rhg0bWLZs2aRhy5YtY2JioksRmZmZmVkvccLdooGBAcbHxycNGx8fZ/bs2V2KyMzM\nzMx6iRPuFo2MjDA8PMy6devYvn0769atY+XKlcydO7fboZmZmZlZD+jqXUokHQucD8wALo6I86o+\nXwmcDQjYArwlIn7Y8UDrWLFiBQDLly9nYmKCRYsWMWvWLN+lxMzMzMyALrZwS5oBfBI4DlgErJC0\nqGq0nwHHRMTvAR8ALupslNmsWLGCxYsX8+IXv5jbbrvNybaZmZmZ7dTNLiXPB+6KiLsjYhtwBXBC\n5QgR8V8R8Yvk7beBBR2O0czMzMysJd1MuOcD91a8vy8ZlmYY+HKhEZmZmZmZtVlfPGlS0hClhHtZ\nnXFOB04HWLhwYYciMzMzMzOrr5st3PcDB1e8X5AMm0TSs4GLgRMiIvXm1hFxUUQsjYilvkOImZmZ\nmfWKbibcNwFHSjpM0tOAk4GrK0eQtBD4AnBKRPy4CzGamZmZmbWka11KIuJJSW8DrqN0W8DPRMR6\nSWcmn38K+L/AbOCfJQE8GRFLuxWzmZmZmVleXe3DHRHXAtdWDftUxes3A2/udFxmZmZmZu3iJ03a\nlDU4OMjg4GC3wzAzM7Npzgl3n3DyaGZmZtafnHCbmZmZmRXICbeZmZlZRv7F2ZrhhNusAN4hm5mZ\nWVnDhFvSUZK+Lum25P2zJf1t8aH1p7GxMdavX8+NN97IkiVLGBsb63ZIZmZmZtZFWVq4Pw2cC2wH\niIgfUXpIjVXZuHEjIyMjrFmzhq1bt3LBBRcwMjIy7ZJut+6aFcvbmJlZf8mScD89Ir5bNezJIoLp\nd5s2bWLVqlUMDQ0xc+ZMhoaGWLVqFaOjo90OzaY5J2hm04+3e7PekSXhfkTSEUAASDoJeLDQqPrU\nxMQEy5YtmzRs2bJlbNiwoUsRmZmZmVm3ZUm43wpcCDxL0v3AO4G3FBpVn5o9ezbj4+OTho2PjzMw\nMNCliGpzq4dNVa7bZmbWixom3BFxd0S8DJgLPCsilkXEPYVH1ofmzp3L8PAw69atY/v27axbt47h\n4WFGRka6HVoqJyhmZlbNx4bOmgrlPRWWoUhZ7lLylKTzgF9HxJZk2C2FR9aH5s2bx+joKMuXL2eP\nPfbgrLPOYnR0lBUrVnQlnqlQ+RstQzvvClNUeU2F9ZCH79RjZmbtMlWOobtnGGc9pcT8ekl/GhGP\nAio2rP61YsUKLrzwQgBuuOGG7gaTU7lC90vcY2NjO+8Ks2zZMsbHxxkeHgYo7CSn38qo09qxTlzG\nZtOb9wE2FWXpw/1kRLwbuBj4pqTnklxAad3Ta2d8WeJpd8yjo6O+K0yPafc66bV63oypsAxTgdfD\nb7kszDovSwu3ACLic5LWA6uBhYVGZZbBhg0bpuRdYbK07nS6BSjr/PplnfRjC1o/xgz9G7dNH66j\n1glZWrjfXH4REbcBfwy8vbCIrBB5WzT6oQVkYGCgL+4Kk0Uny7vIeU2ldWLTQz/s66w/9Epd6pU4\naunl2IqWmnBLekny8hBJryv/AS8FftWR6MzqGBkZ6bu7wkx13VgnU2EH3uoypH2/cnivlVMz8fTa\nMmTVL3H3S5xTxXQv7+m2/PW6lBwDrAVeXeOzAL5QSERTWN6frbKOPx1/Disvc/muMBMTEyxatCj3\nXWFaKbvK7zYznV5Yb+2MIcs6aaXMOlFera7TXtCuOt1NvRJHHv26D7D+0U/1pahY+3kfnZpwR8R7\nkv9v6lw401M/VZp2VfZ2LXPaXWHKt6abmJhgyZIljIyM1Ez6io6vU9PtJf12p55eqMfdmn4z80gb\nfyrU7XYtQ73pZCm/Isqy0/u0dtWTIvbXnT456sVto9eXpxfLrFUNL5qU9A7gX4EtwKeB5wDnRMT1\nBcc2JbRjR1AreeykTmyY7UzkN27cyBNPPDHp1nQrV67kAx/4ALfffnvuaeZVvb527NjBvHnzWp5u\n3gNPK3WvmYShFb1wEO70/Hr9gNIrLUm9UE5pZdHp2Io42elmA0QRcbRz/WQ5Hpf370WdgE7lE9zp\nJstdSv48Is6X9EpgNnAK8G+AE+4KeSt91o1l48aNNe9rvOeeezaVxFXuLA4++GAef/xxNm/e3LbE\nsKhkE7KX2aZNm1izZg1DQ0MADA0Ncfnll7N8+fK2xFFPrftQr1y5sunpTYWdaqeWoci6V6mVxGAq\nrM8sOt3iWEQra6+vq16Mr5Mnnd1Y/lrH41b272mKaoEvcntodf+bNv2s0+3F7aHabhnGKT/k5njg\nsohYjx980zGbNm3a5b7Gb3zjG3nggQd2Pslv48aNu3yv+ml/GzdunLSzuPTSS9ltt9248sor2bp1\nKyeeeGLDaTZSOf2tW7dywQUXsHnzZjZu3Fgznkaa+Q7AxMREzVvTTUxM5F4myHdhR637UF9++eVs\n2rQp0/ebWeY88bVbkfPOUxb16l4r0+3kMvSjvOu/E3W1m9tDEbpRxo3qbS+u92bleTJureNxef+e\ndVvvdFm0sq6yfHdwcJBFixbV3f+Wp5N3f5hnv94PsrRwf0/S9cBhwLmSngnsaMfMJR0LnA/MAC6O\niPOqPlfy+fHAr4HTIqJvHiufpTU57eytcnhl8jg2Nsbq1au56qqrWLZsGe973/v4xCc+wR133LHz\n+8CkLhXlcZ566imuvvpqhoaGWLJkCZdccglDQ0O7TDPLWXutuNNalk844YTcrQL1WhIanfHOnj2b\n8fHxnXFA6dZ0s2fPzrcC6yxz5fqsfA00ney3s/Ukb2tDvbo6NDTUdD3OOo/q8avLolY9r5zvE088\nsbNuw66/apQPGmeccUam6bZD1jpcXRaV1xtk0UzLUrOtVUUPb1Q2eaab9xe8on4hade20Up8edYD\nkHk/1KlflfIuQ9YW0LxPxq3XmFP0fqWdZd3u9VZ97H/ooYd4+tOfvktekvf4lmW6RdW3QkRE3T9K\nreDPAWYl72cDz270vQzTnQH8FDgceBrwQ2BR1TjHA1+m1KL+AuA7Wab93Oc+N7rpmGOOiYGBgTjs\nsMNi7dq1cdlll8XChQtj7dq1sW3bthgZGYl99tknJMWBBx64c/jatWvjwAMPjPnz5+/87uLFi2Pt\n2rU7p135fvXq1TvHq5zuXnvtVXOc3XbbLbZt2xYRMel19TwiItauXRtz5szZuTzHHHNMrF69OubM\nmZMat6Sd0yzPe/HixZPiqTX9yukuXrw4BgYGYs6cOZO+s3r16jjiiCNS5z0wMNAwvupxyvM65phj\nJi1no2WuXJ/V6/aII45IXdZG8663zAsWLIh99923btxZlr/WcjaqqwcccEBT9ThtndSbTq31n1bP\nK6dZWZ8r61512e27776p081ST7IML8+vut7XqsPVZbF27do47LDDYvXq1Znqar1lyFKna9WryvpQ\nbz0XMXzBggWp21at/WPadNO+m7Zu6y1zZRnl3Q4rlyetrjdzrMgTX/WypZVjeXjlNlJdbxcvXhxv\ne9vbau5nai1bvfqWpVzT9ml59z1p86pe1vL+evHixbmOS2nH3GaOP2lx1yvr1atX19zWsxwfstSl\nLK/L+98seUmtelVZFuW46023XrmuXr06Oixb3pt1xHb/AS8Erqt4fy5wbtU4FwIrKt7fCRzYaNq9\nkHBXbphpSXJ1optl401LlNMS68px0l5XJysREdu2bQtJNTfStAR9n332aRhP+bPKZCjtAFNrQ0ub\n9zOf+cxMB8y8B4hay5z2uhxr3iShPO+0Zc6aPDRaP9UnUI3qar1lzbtOskwn7eQty3zTxqkuu3on\nmmkndXlOkCvnVzmvLHFXlt2CBQsyJdJ513WW5C5t/5NlPbQyvPqENe9+Jku5Vu8naiWZafUny4li\ndXJbuTxZEpG8ZVYvvrQEulY5Vg5vlDyVT8CzNP5kaajIcmKV9xiatYGk3rGv0UlK2raedb9Sa5tO\n2z4b7TPnzZuX+YQlb13K8jrLvjhtf1hrWyqf1KVtG43Ktdxo0UE9n3CfRKkbSfn9KcAnqsa5BlhW\n8f7rwNJG0+6FhLtyp5W2MaZVwHoJauXOM0tinTaPtI2lPF51Ra6XoFfG1yhBrXeAzXLwTNtBprWi\npyV9aYle2k4+7deBWvFcdtllsddee+VusU3bueRNaLOeQNVqnaj+bpZfRbKsk7y/rlQmA3m3n3pl\nl/f7eU+Q857gZq3Paa1BWVv2a7W65z1hL3p4o7KpXLZG0827n0ibZtb9WForXt4GkrxllmU/Uf3L\nY7ks0obXS3KyxJRlW8wSdxH7gHrLkzbv8rSrf21MK6csMVVPs9ZxIus2k6WRp5W6lHd7yDv9eid1\neRswKtfh4sWLo4OmV8INnA7cDNy8YMGCuPXWW3f5K6+shx56qNDPX/jCF8Z+++2Xa8PJsgMqd0s4\n9NBD625EeVoAZs2aFZJqVvC8O9dmWqhbabGvLJssB492JXpZ1tV+++0XS5cujQ9+8IOx3377ZWqx\nTdu5pCUeeZetMkHP0jpRb1nzrpO8J1CS4pBDDslcD8vlUq7b9couz0ld3hPkLOPkTZLSTg4PP/zw\nOOqooxrWpSz1qpWyaNfwenGkrbdWyrU83XrbUpbyqtdY0sq+P0uZtXLyVq+Rorzt1TsBb1es9fYB\nzZ58ZJ1Xll/Y8o6Tt/ErbZrNbDOV63Hvvfdu+/pp5biU9Ziet2Ejrf7stttuHcv/ol0JN/BvWYbl\n/ZvqXUry/jSUZWMs//SUpRtBnp+xKn96ypLEZTm7rNdCkOWAVP5O9QacVjZpP5tm2WlXzitvYpzW\nB7f65+W8SWMrLZFZyijLstVrbcg7v2Z+XanVFaiV5Sxr1Ae8Xa1BaeuzUXehww47bOdJddoyrF37\n2xPwZlsT856wFz28/Fmjn6zzdrVopU7mTYbrLU+W8s67zFlOMtJOtOt1kUm7DqFdJ0GtbDPtalkt\na3QNSTPHuEb7lSzTzLrNtHIsa2WfkaW+ZVm39U7qysonYHka4RYv7sMWbuCWqvczgNuzzqDOdHcH\n7qZ095PyRZOLq8b5k6qLJr+bZdq9kHDnvfihmQsKsl4oVzlO9bSq+3hm2UjTWrLLqjeOeolx2gGp\nrHoDrp53eXnSLgxKO8DmPYCnHZBqXfCzevXqXfpGZzk4pZVf3r5yaWXUqOW3vJOuXp7Ki6Ty9IfP\neiFa3gvd8sy33glR2olm3oN7lkSvXh/+WvWnXqtpZR0pz79cNo1aolrpVlYddxHD166d3B+1Vren\n6vqa9wKwPI0UWfqsNmosqVyeLPPOso1l6Z9ca/vJciFvo3pbPgFPiylvQ0XWvt2N6lKWk9pG+4OI\nXRvMsuyjs5R3vW26Vn3OuvxZumXWOz7k6WOft75lqVeNWrgjdr3xQKNGuL7rw520OG8BngQeS/62\nABPAP2SdQd2Zl+5C8uPkbiUjybAzgTOT1wI+mXx+a5b+29FDCXfa61oJcFrCWOtuA3mmWy+OWq+z\n3h2iVktAWfXG0Wg5ax2QmjlpyDK/PAeItGVuVI7l91mShCw7l7TEI2+sWVptal1s10p9yzNOo7vZ\nNDvftBOitANs3oN7eV00uqNM1u2zcppZyqXR+m0mucuy7RU1vNb6ybufybpu0xKptPqTN3GvXp4s\nDSRZ6nqW+GrVwyxd7xqVb+UJeFpMzd4FI2vDU959fZYGkmb2S83WvXoJY72T8bT1n+WXx3bdRaaV\nu/ZkrS+VJ3W1to1667q8TjusbS3cbUmuO/nX6wl32us8B+S0abUaR6OK3GiHUu9EIW050w5IzZ40\nNLPBpyXGzcRQfp81SchzC796tw5sFGujlt+8J3jtfl0vAWp1+tWarSfNnCC3Emv1eqs1r1rjZWnV\nKx+UGi1/u9dzluWv/ixvg0TldNK2w7S7dKTd4aVefam3DWf5frvLO8svfpXLnGfby1s2eZat1X1A\nvflWyrM8rW7rjfYfeRo/0l4XdSzL8jqtXPOMk/WkrlG51ppfh2TKTbM8+Oa7kvaJiF8CSJoFDEbE\nFzN8d9pLe7BEreErVqzgwgsv3Pl5N57MNW/ePM444wyWL1/OxMQEixYt4qCDDuLoo4+uGVOt8WfN\nmlX3ZvTVywnsfF9EGaRNszx8/fr1kx6Us2LFCg444ACWL1/Obbfd1vR8586dy/DwMKtWrZp0o/+5\nc+dOiilP+VVP84ADDuDXv/41z3rWs7jtttvqllf1upo/fz6vf/3r2bx58y7zzfp43HY+RreoBxZl\n1aieVA/PW++rZS27rNtYvfVbva5HR0d3ebBHM9teu9Z/o+nUKoOPfOQjO5eh3n52xowZNR9W8thj\nj6Vun1k0uw3X+n698s5bxjfccANLlizZ+UREmPwQsuplHh4e5pBDDsm17WU5rlXv47NotA+onFan\njqHNHOPqxVQ2Nja2c5qtPhm5nceyRmVchFpldOutt+58X2s99kK+lFeWhPs9EXFV+U1EbJb0HsAJ\nd5t0qlJn1cyBoFcqfjNlWS8xbkXWHXWeHUejxKNSngNSnnIrqr4WtR4aaSahKUtbb+0uo+rk7u1v\nf/vOp62NjIzUHK96/eZd13kaC6qH5/1u1mk1k8QNDAzUTOLKj6NO2z7TYigq0ctbfvXKYMOGDTWT\nuMcff5wLL7xw0jKPjo4C1N326s2rlbpevTxjY2NN7QNarW9Zptvuk9G0JBl+Wz/zqrUPHR4e5uMf\n/3ju+PpNPyxPloR7tya/Z1VaObC3c7qtKGpeWabbrnlXT6fZFowsB79WD8TtTDy6JUviUn5d2eqT\nNQHKMt9O6ERMtR4XPzw8zJ577tlfjzjuoJGRkZpJSLmVvxPbUt6628r0If0kY2BgoO4yt/KLTbOx\nVio3HGRpUCgqhk5PKy1JLp8I5Tnxq3Usq/VrVjPL0K79bztPgvpNlsT5ZkkfpXTxIsBbge8VF9L0\n1otJdiuKOtC0++SlUz9PFZ2UdbNetGPe3ezWUG+63SrX6vlu2rSJNWvWTOoqsGrVKs4666ymY+zF\nk5d2qpXEZUlC+lm9k4w0vfITfadOgorQTLzNJsn15t2P5ddPsTYrS8J9FvB/gM8BAXyVUtJtfaTX\nkrh605wOG16/6ZWEs9en26ys8aT199ywYUMBURWrk+ugqCSuqFa8VvXCSUY/7TN6YX31Y5Jcrd0N\nYc1Ot1c1TLgj4nHgHEnPSF5bn8vzE5X1v7GxMdavX8/ExARLlixhx44d7n7QQe3crtIuKhsYGGjb\nPNrN+5XumCpJ3FTg8s9vKpZZw4Rb0h8BFwN7AQsl/T5wRkT8r6KDs9ZMhQpb9DLkbaHqtzIdGxvb\npc/vypUrux1WYTq9fjo9v0b9PVvVb/W737m8+4fXVe/pt3WSpUvJPwGvBK4GiIgfSnpxoVGZTQOd\n2FmMjo6C2IypAAAOnklEQVTWvD3Y8uXLOxZDP+mF8qh3N5F29fc0a6QXtgWzqSTT3UYi4l5JlYOe\nKiYcmwq8o+4dabcHy3OPV+stzXYVmArbZS8vQy/HZr3JdWZ6yZJw35t0KwlJM4F3AP13hY5ZB/XK\njrTePYhtavJ1GWbF8nZlzciScJ8JnA/MB+4Hrsd3KZnyvEPprKLKu5nbg5mZmVl7pSbckj4YEWcD\nQxExda+y6hNOgPNzmfXG7cGsPVyfi+FyNbNOqNfCfbykc4BzgX/vUDzWp3zQ6l2+PZiZZeH9g1lx\n6iXcXwF+Aewl6TFAlB58IyAiYu8OxNeXvNOyormOmVmzvP8w67x6CfffRsRfS/pSRJzQsYjMcvLB\nw8zMrMTHxN60W53PvpX8f6wTgZiZmZl1S/mpvDfeeCNLlixhbGys2yHZFFKvhftpkt4I/JGk11V/\nGBFfKC4ss6ljOrY2TMdlNrP+VeupvMPDwwC+yNzaol7CfSawEpgFvLrqswCccNuU4iSxN3g9mFmn\n1Xoq76pVqzjrrLOccFtbpCbcETEOjEu6OSJWdTAmMzMzm6a6cdKd9lTeDRv8nD9rj9Q+3JLeDRAR\nqyS9oeqzvy86MDMzs3a54YYb/OuJpSo/lbfS+Pg4AwMDXYrIppp6F02eXPH63KrPji0gFjMzM7OO\nKz+Vd926dWzfvp1169YxPDzMyMhIt0OzKaJeH26lvK713szMzKwvdfupvP71Zeqrl3BHyuta73OR\ntB/wOeBQ4B5geUT8omqcg4HLgP2T+V0UEee3Ml8zMzOzWvxUXitSvS4lvy/pMUlbgGcnr8vvf6/F\n+Z4DfD0ijgS+nryv9iTwrohYBLwAeKukRS3O18zMzMyso+rdpWRGgfM9ARhMXl8K3ACcXTX/B4EH\nk9dbJG0A5gO3FxiXmZmZmVlb1WvhLtL+SUIN8BClbiOpJB0K/AHwnWLDMjMzMzNrr3p9uFsi6WvA\nATU+mnTJb0SEpNQ+4ZL2Aq4E3hkRqY+Zl3Q6cDrAwoULm4rZzMxsqnA/ZLPeUVjCHREvS/tM0sOS\nDoyIByUdCGxMGW8mpWT78kaPko+Ii4CLAJYuXdrSRZ1mZmZmZu3SrS4lVwOnJq9PBb5UPYIkAauA\nDRHx0Q7GZmZmZmbWNt1KuM8DXi7pJ8DLkvdIOkjStck4LwJOAV4i6QfJ3/HdCdfMzMzMrDmFdSmp\nJyImgJfWGP4AcHzyehw/YMfMzMzM+ly3WrjNzMzMzKYFJ9xmZmZmZgVywm1mZmZmVqCu9OE2M7Pm\n+f7KZmb9xS3cZmZmZmYFcsJtZmZmZlYgJ9xmZj1ibGyM9evXc+ONN7JkyRLGxsa6HZKZmbWB+3Cb\nmfWAsbExRkZGWLNmDcuWLWN8fJzh4WEAVqxY0eXozMysFW7hNjPrAaOjo6xatYqhoSFmzpzJ0NAQ\nq1atYnR0tNuhmZlZixQR3Y6h7ZYuXRo333xzt8MwM8tsxowZPPHEE8ycOXPnsO3bt7Pnnnvy1FNP\ndTEys+ljcHAQ8J2ALJdMT0V3C7eZWQ8YGBhgfHx80rDx8XEGBga6FJGZmbWLE24zsx4wMjLC8PAw\n69atY/v27axbt47h4WFGRka6HZqZmbXIF02amfWA8oWRy5cvZ2JigkWLFjE6OuoLJs06yF1JrCju\nw21m1kPch9TMrK+4D7eZmZmZWbc54TYzMzMzK5ATbjMzMzOzAjnhNjMzMzMrkBNuMzMzM7MCOeE2\nMzMzMyuQE24zMzMzswI54TYzMzMzK1BXEm5J+0n6qqSfJP/3rTPuDEnfl3RNJ2M0MzMzM2uHbrVw\nnwN8PSKOBL6evE/zDmBDR6IyMzMzM2uzbiXcJwCXJq8vBU6sNZKkBcCfABd3KC4zMzMzs7bqVsK9\nf0Q8mLx+CNg/ZbyPAe8GdnQkKjMzMzOzNtu9qAlL+hpwQI2PRirfRERIihrffxWwMSK+J2kww/xO\nB04HWLhwYVMxm5mZmZm1W2EJd0S8LO0zSQ9LOjAiHpR0ILCxxmgvAl4j6XhgT2BvSZ+NiD9Lmd9F\nwEUAS5cu3SWBNzMzMzPrhm51KbkaODV5fSrwpeoRIuLciFgQEYcCJwNr05JtMzMzM7Ne1a2E+zzg\n5ZJ+ArwseY+kgyRd26WYzMzMzMzarrAuJfVExATw0hrDHwCOrzH8BuCGwgMzMzMzM2szP2nSzMzM\nzKxATrjNzMzMzArkhNvMzMzMrEBOuM3MzMzMCuSE28zMzMysQE64zczMzMwK5ITbzMzMzKxATrjN\nzMzMzArkhNvMzMzMrEBOuM3MzMzMCtSVR7ubmVltN9xwQ7dDMDOzNnMLt5mZmZlZgZxwm5mZmZkV\nyAm3mZmZmVmBnHCbmZmZmRXICbeZmZmZWYGccJuZmZmZFcgJt5mZmZlZgRQR3Y6h7SRtAn7ehVnP\nAR7pwnz7lcsrP5dZPi6v/Fxm+bi88nOZ5ePyyq+TZfZIRBzbaKQpmXB3i6SbI2Jpt+PoFy6v/Fxm\n+bi88nOZ5ePyys9llo/LK79eLDN3KTEzMzMzK5ATbjMzMzOzAjnhbq+Luh1An3F55ecyy8fllZ/L\nLB+XV34us3xcXvn1XJm5D7eZmZmZWYHcwm1mZmZmViAn3G0g6VhJd0q6S9I53Y6nF0k6WNI6SbdL\nWi/pHcnw90q6X9IPkr/jux1rr5B0j6Rbk3K5ORm2n6SvSvpJ8n/fbsfZKyT9bkU9+oGkxyS903Xs\ntyR9RtJGSbdVDEutU5LOTfZrd0p6ZXei7q6UMvuwpDsk/UjSVZJmJcMPlfSbirr2qe5F3h0p5ZW6\nDbqOpZbZ5yrK6x5JP0iGu46l5xM9vS9zl5IWSZoB/Bh4OXAfcBOwIiJu72pgPUbSgcCBEXGLpGcC\n3wNOBJYDv4qIf+xqgD1I0j3A0oh4pGLYh4BHI+K85ORu34g4u1sx9qpku7wf+EPgTbiOASDpxcCv\ngMsiYkkyrGadkrQIGAOeDxwEfA04KiKe6lL4XZFSZq8A1kbEk5I+CJCU2aHANeXxpqOU8novNbZB\n17GSWmVW9flHgF9GxPtdx+rmE6fRw/syt3C37vnAXRFxd0RsA64ATuhyTD0nIh6MiFuS11uADcD8\n7kbVl04ALk1eX0ppJ2O7einw04joxgOwelZE3Ag8WjU4rU6dAFwREVsj4mfAXZT2d9NKrTKLiOsj\n4snk7beBBR0PrEel1LE0rmPULzNJotQwNdbRoHpYnXyip/dlTrhbNx+4t+L9fTiRrCs5Q/8D4DvJ\noLOSn2Y/4y4SkwTwNUnfk3R6Mmz/iHgwef0QsH93Qut5JzP5AOU6li6tTnnfls2fA1+ueH9Y8lP/\nNyT9cbeC6kG1tkHXscb+GHg4In5SMcx1LFGVT/T0vswJt3WUpL2AK4F3RsRjwL8AhwNHAw8CH+li\neL1mWUQcDRwHvDX52XGnKPUHc5+wKpKeBrwG+PdkkOtYRq5T+UgaAZ4ELk8GPQgsTLbbvwRWS9q7\nW/H1EG+DzVvB5MYD17FEjXxip17clznhbt39wMEV7xckw6yKpJmUNo7LI+ILABHxcEQ8FRE7gE8z\nDX9OTBMR9yf/NwJXUSqbh5P+a+V+bBu7F2HPOg64JSIeBtexDNLqlPdtdUg6DXgVsDI5uJP8ZD2R\nvP4e8FPgqK4F2SPqbIOuY3VI2h14HfC58jDXsZJa+QQ9vi9zwt26m4AjJR2WtKydDFzd5Zh6TtIP\nbRWwISI+WjH8wIrRXgvcVv3d6UjSM5KLQZD0DOAVlMrmauDUZLRTgS91J8KeNqlFyHWsobQ6dTVw\nsqQ9JB0GHAl8twvx9RxJxwLvBl4TEb+uGD43uWAXSYdTKrO7uxNl76izDbqO1fcy4I6IuK88wHUs\nPZ+gx/dlu3d6hlNNcpX624DrgBnAZyJifZfD6kUvAk4Bbi3f3gj4G2CFpKMp/fRzD3BGd8LrOfsD\nV5X2K+wOrI6Ir0i6CVgjaRj4OaWLaSyRnJy8nMn16EOuYyWSxoBBYI6k+4D3AOdRo05FxHpJa4Db\nKXWbeOt0u3sEpJbZucAewFeTbfTbEXEm8GLg/ZK2AzuAMyMi6wWEU0JKeQ3W2gZdx0pqlVlErGLX\na1HAdQzS84me3pf5toBmZmZmZgVylxIzMzMzswI54TYzMzMzK5ATbjMzMzOzAjnhNjMzMzMrkBNu\nMzMzM7MCOeE2M+txkvaXtFrS3ZK+J+lbkl7b7bjMzCwbJ9xmZj0secjDF4EbI+LwiHgupfvzLuhu\nZGZmlpUTbjOz3vYSYFtEfKo8ICJ+HhEXAEiaIenDkm6S9CNJZyTDByXdIOnzku6QdHmSvCPpuZK+\nkbSWX1f1JECScS6R9HFJ/5W0rJ9UMd1rKsb7RPKYcyTdI+kfJP1A0s2SnpNM/6eSziyykMzMepkT\nbjOz3rYYuKXO58PALyPiecDzgL9IHl8M8AfAO4FFwOHAiyTNBC4ATkpayz8DjKZM+0BgGfAqSk9x\ny+K/I+Jo4JvAJcBJwAuA92X8vpnZlONHu5uZ9RFJn6SUBG9LkuxXAM8ut0AD+wBHAtuA70bEfcn3\nfgAcCmwGlvDbx5LPAB5Mmd0XI2IHcLuk/TOGeHXy/1Zgr4jYAmyRtFXSrIjYnH1pzcymBifcZma9\nbT3w+vKbiHirpDnAzckgAWdFxHWVX5I0CGytGPQUpX2+gPUR8cIM8678vpL/TzL519E9U76zo+r7\nO/Axx8ymKXcpMTPrbWuBPSW9pWLY0yteXwe8JekqgqSjJD2jzvTuBOZKemEy/kxJi3PE83NgkaQ9\nJM0CXprju2Zm05JbG8zMelhEhKQTgX+S9G5gE/A4cHYyysWUuorcklwUuQk4sc70tiXdTz4uaR9K\nx4GPUWpJzxLPvZLWALcBPwO+39SCmZlNI4qIbsdgZmZmZjZluUuJmZmZmVmBnHCbmZmZmRXICbeZ\nmZmZWYGccJuZmZmZFcgJt5mZmZlZgZxwm5mZmZkVyAm3mZmZmVmBnHCbmZmZmRXo/wP72zp4jwyU\nMwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11878dd68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(12, 3)\n",
    "plt.errorbar(np.arange(beta.shape[1]), beta.mean(0), yerr=2*beta.std(0), fmt='o', color='k', markerfacecolor='w');\n",
    "plt.axhline(0, lw=1, ls='--', color='lightgrey', zorder=-1)\n",
    "plt.xlabel('Gene num')\n",
    "plt.ylabel('Effect size');\n",
    "plt.title('Sparse Reverse DE Results');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the uncertainty is not very large for most of the genes! A small number of the genes have larger effect sizes, though with pretty large CI's. We can look at the volcano plot to get a clearer summary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pval = np.sum((beta < 0), 0) / beta.shape[0]\n",
    "pval = np.minimum(pval, 1 - pval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAADgCAYAAABYS05fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X18VOWZ//HPlWwgQGACyHOMhIfSkGhdwYfWFI2tFm0R\n9te6GrVFm9baCqyKv9Y6ru1uN1u1vvqrP9r1ocW2rCbY2pYqda0PyW7NCi1qUbFpq6IuoBVEDYhG\nAlz7xzmTTkJCJmEmk5n5vl+v82LO83UOyZV77vs+9zF3R0REkicv3QGIiGQbJVYRkSRTYhURSTIl\nVhGRJFNiFRFJMiVWEZEkU2IVySFmNtXM3Mz+Jt2xZDMl1gxjZlVm9piZtZrZG2b232Z2fLrjOhQz\n+5GZ7TWzt8OYHzKz96c7rsNhZi+Z2btmttvM3gr/Ty41s7y4beKvOzY91cPxTjWzA+E2u83sT2Z2\n8cBdkSSTEmsGMbNRwFpgBTAGmAL8E/BeCs6V7BLNje5eRBDzNmBlko/fqxRc0wJ3HwkcBVwPfIWD\nr+tGdy+Kmz5wiOO9Et6jUeGxvm9ms5McswwAJdbM8j4Ad29w9/3u/q67P+juTwOY2UVhCfa7YYn2\nj2b2kdjOZnaxmbWEJaLNZvaFuHWnmtlWM/uKmf0F+KGZHWFma8MS2Rtm9misRGZmk83sZ2a2w8xe\nNLNliVyAu78L/AQ4Nn65mX02jO1NM/u1mR0VLr/FzG7qsu0vzezK3uIws6+b2T1mdqeZ7QIuMrMT\nzOxxM9tlZq+Z2bfjtj8pLHm+ZWZPmdmpCV5Tq7vfC5wLLDazykT2O8Tx3N3XAG8CByVWMzvXzB7v\nsuwKM7s3/PxxM/t9eI1bzOzrPZ0rLHl/NG7+62Z2Z9x8j/ck/HnbHP48vWhmFxzOdWcVd9eUIRNB\nSWYn8GPgTGB0l/UXAfuAK4ACgl/0VmBMuP7jwHTAgFOAd4DjwnWnhvveAAwFhgHfBG4Nj1UAfDjc\nNw94ArgOGAJMAzYDH+sh7h8B/xJ+HgH8O/BU3PqFwPNAOfA3wLXAY+G6ecAWwML50cC7wOTe4gC+\nDrQDi8JthwHrgE+H64uAk8LPU8J7e1a47enh/Lgerukl4KPdLP8f4ItdrzuB/9tTga3h5zzg78LY\nZ3Wz7XBgNzAzbtkG4Ly4Yx0dHucY4DVgUbhuKuDA33R3HeE9u7O3exL+P+6KxQdMAirS/TsyWCaV\nWDOIu+8Cqgh+Mb4P7DCze81sQtxm24HvuHu7u98N/IkgoeLuv3L3FzzwX8CDBMky5gDwNXd/z4OS\nZTvBL8xR4fEe9eC36HiChPPP7r7X3TeH8Zx3iPCvMrO3CBJCFfDpuHWXAt909xZ33wf8K3BsWGp9\nNLzeWJyfAta5+ysJxrHO3de4+4G4a5phZke4+9vuvj7c7kLgfne/P9z2IeBxgqTSF68QVNN0uu64\n6ceH2HdyeI9eB75G8AfgT103cvd3gF8CNQBmNhN4P3BvuP4/3f2Z8DqeBhoI/pD2VW/35ABQaWbD\n3P1Vd3+2H+fISkqsGSZMPhe5ewlQSVBy+07cJtvC5BfzcrgNZnamma0Pv9a/RfALckTctjvcvS1u\n/lsEJckHw698V4fLjyJMArEJuAaIT/Bd3eTuxQQlpneBWXHrjgJujjvWGwQl4ynhtawmTCLA+cBd\nfYhjS5c4agmqVP5oZhvM7BNxxzqny7GqCP6w9MWUMP5O1x03LT7Evq+E24xx92PdfTWAmd0a1/h1\nTbhtPZ3vyZow4WJmJ5pZU1g90krwh+uIridLQI/3xN33EHwjuhR41cx+ZRneIJlMSqwZzN3/SPB1\nM75Ob4qZWdx8KfCKmQ0FfgbcBEwIk9z9BAms45Bdjr/b3Ze7+zTgbODKsM52C/Bil4Qx0t17Ld25\n+/8A/0CQSIeFi7cAX+hyvGHu/li4vgH4VFiCPTG8jth+vcXR9Zqec/caYDxBtcc9ZjYiPNa/dznW\nCHe/vrdrirGgd8YUoDnRfRLh7pf6Xxu//jVc/BAwzsyOJUiw9XG71BOUXo909whBdY7RvT0EVQsx\nE+M+H/KeuPuv3f10gj8+fyT4tiAosWYUM3u/mS03s5Jw/kiCX6r1cZuNB5aZWYGZnUNQb3k/QR3k\nUGAHsM/MzgTO6OV8nzCzGWGibgX2E3z9+x2w24KGrmFmlm9mlZZgt6/wK+UrwCXholuBr5pZRXje\nSBh7bPvfE3w9/gHwa3d/K1zV5zjM7EIzG+fuB4DYcQ4AdwILzOxj4XEKLWjQK+ntesxsVFjyXU1Q\nP/lMIvfhcLh7O/BTgm8VYwgSbcxI4A13bzOzEwhKtD3ZCJwX/rzMJahqienxnpjZBDNbGP5Reg94\nm+A+CkqsmWY3QYntt2a2hyChbgKWx23zW2AmQSKqAz7l7jvdfTewjKBF/k2CX7Z7eznfTOBhgl+a\ndcC/uXuTu+8HPkHQsv8if016kT5cy7eAL5vZUHf/BUHpcbUFrfebCBrn4tUDHyWuZNbPOOYDz5rZ\n28DNBA0+77r7FoJGtGsI/vhsAf4vh/4duc/MdofbRoFvA137nn7ZOvdjff0Qx+ur2D35aVg3HfMl\n4J/D2K4j+D/vyT8SNGi+SdB1L/7+Huqe5AFXEvyBfIOgDveLSbmqLBBraZUsYGYXAZ9z96p0xyKS\ny1RiFRFJMiVWEZEkU1WAiEiSqcQqIpJkSqwiIkmWlWMyzp8/3x944IF0hyEi2aenBy06ycoS6+uv\nJ7OroIhI32RlYhURSSclVhGRJMuqxGpmC8zs9tbW1nSHIiI5LCv7sc6dO9cff/zx3jcUkQ7t7e1s\n3bqVtra23jfOcoWFhZSUlFBQUNB1VUKNV1nZK+BwNDQ0UFdXR0tLC+Xl5USjUWpqanrfUSTDbd26\nlZEjRzJ16lQ6jzyZW9ydnTt3snXrVsrKyvp1DCXWOA0NDUSjUVauXElVVRXNzc3U1tYCKLlK1mtr\na8v5pApgZowdO5YdO3b0+xhZVcd6uOrq6li5ciXV1dUUFBRQXV3NypUrqaurS3doIgMi15NqzOHe\nh6xKrIfbeNXS0kJVVecR96qqqmhpaUlGeCJZZf/+/axbt47f/OY3vPde8t7A/tprr3H++eczbdo0\n5syZwwc/+EF+8YtfJO34AyGrEqu73+ful0QifRlv+a/Ky8tpbu78Vo3m5mbKy8uTEZ5IxmltbaWh\noYFVq1bx2muvdSzfsGEDM2fO5Atf+AJXXXUVRx11FGvXrj3s87k7ixYtYt68eWzevJknnniC1atX\ns3Xr1sM+9oDyQfCq2GRPc+bM8f6or6/3srIyb2xs9L1793pjY6OXlZV5fX19v44nkkn+8Ic/dJpf\nu3atjxkzxhcsWODnnHOOFxcX+2233eZ79uzxiRMn+s9//vOObR977DEfO3asv/zyy97a2urLli3z\n8ePH++jRo33x4sW+bdu2hGJ4+OGHfd68ed2u27dvn1911VU+d+5cP/roo/3WW291d/empiY/5ZRT\n/JOf/KTPmjXLzz//fD9w4IC7uz/++OM+b948P+644/yMM87wV155xd3db775Zi8vL/ejjz7azz33\n3ITuRyihHJT2JJiKqb+J1T1IrhUVFZ6Xl+cVFRVKqpIz4hPJm2++6WPGjPH169d3LHvhhRd83Lhx\nftNNN/n8+fMP2v9LX/qS19XV+bx583zx4sX+wgsv+LZt2/wrX/mKz5w50/fs2dNrDDfffLNffvnl\n3a677bbb/Bvf+Ia7u7e1tfmcOXN88+bN3tTU5KNGjfItW7b4/v37/aSTTvJHH33U9+7d6x/84Ad9\n+/bt7u6+evVqv/jii93dfdKkSd7W1tZxrb3djzgJ5aBB3yvAzBYBHwdGASvd/cFDbLsAWDBjxox+\nn6+mpkY9ACTn/fKXv+TUU0/lxBNP7Fg2bdo0Fi9ezCOPPEJJycHvWCwpKWHjxo3s3LmTpqYm8vKC\nmsbrr7+eTZs2cffdd3PxxV1fCXZol112Gc3NzQwZMoSjjjqKp59+mnvuuQcIqimee+45hgwZwgkn\nnNAR07HHHstLL71EcXExmzZt4vTTTweCOuFJk4K3mR9zzDFccMEFLFq0iEWLFvX9BvUipXWsZnaH\nmW03s01dls83sz+Z2fNx76rvlruvcffPE7y//Nxetj2sOlYRCbS1tVFUVHTQ8qKiIiZPnsx9991H\nfCNxe3s7q1evZuzYscybN68jqcaceuqpPPNM7y+vraio4Mknn+yY/973vscjjzzCjh07cHdWrFjB\nxo0b2bhxIy+++CJnnBG8aHjo0KEd++Tn57Nv3z7cnYqKio7tn3nmGR58MCiX/epXv+Kyyy7jySef\n5Pjjj2ffvn0kU6obr35E8FbMDmaWD3yP4C2cs4EaM5ttZkeb2dou0/i4Xa8N9xORFDvzzDNZu3Yt\n27Zt61i2e/duVq1axUUXXcR5551HVVUVd9xxBw0NDZx22mmUlZVx9tlns379+qCeMc769euZNWtW\nr+c97bTTaGtr45ZbbulY9s477wDwsY99jFtuuYX29nYA/vznP7Nnz54ejzVr1ix27NjBunXrgCD5\nP/vssxw4cIAtW7ZQXV3NDTfcQGtrK2+//XbiNycBKa0KcPffmNnULotPAJ53980AZrYaWOju3yR4\nlXEn4Tvtrwf+w92f7LpeRJKvtLSUa665huOPP57Pfe5zDBs2jDvuuIMzzzyTk08+mZNPPpk1a9aw\nevVq9u7dS21tLRdccAH5+flce+21XH755Vx33XUMHTqU7373u2zYsIEf/vCHvZ7XzFizZg1XXHEF\nN954I+PGjWPEiBHccMMNnHPOObz00kscd9xxuDvjxo1jzZo1PR5ryJAh3HPPPSxbtozW1lb27dvH\n5Zdfzvve9z4uvPBCWltbcXeWLVtGcXFxMm9f6scKCBPrWnevDOc/Bcx398+F858GTnT3JT3svwxY\nDGwANrr7rT1sdwlwCUBpaemcl19+OclXIpLdYo9xx3vyySe5++672bt3LwsXLuSUU07ptfP866+/\nzvLly/npT3/Kvn37OOuss/j2t7/NtGnTUhl+0nV3P0hwrIBBn1j7Q4OwiPRdD4mk3w4cOIC7k5+f\nn7RjDqTDSazp6BWwDTgybr4kXCYiWaRrA1YuSceVbwBmmlmZmQ0BzgPuTUMcIiIpkeruVg3AOmCW\nmW01s1p33wcsAX4NtAA/cfdnUxmHiCQm1VWDmeJw70OqewV029Pe3e8H7k/luUWkbwoLC9m5cydj\nx47N6VGu3IPxWAsLC/t9jEH/5JWIDIySkhK2bt16WOOQZovYGwT6S4lVRAAoKCjo94j50lnuNtuJ\niKSIEquISJIpsYqIJJkSq4hIkimxiogkmRKriEiSKbGKiCSZEquISJIpsYqIJJkSq4hIkimxiogk\nmRKriEiSKbGKSNZqaGigsrKS/Px8KisraWhoGJDzanQrEclKDQ0NRKNRVq5cSVVVFc3NzdTW1gJQ\nU9PtUNFJk/KXCaaDXiYoIpWVlaxYsYLq6uqOZU1NTSxdupRNmzb197CD4y2t6aDEKiL5+fm0tbVR\nUFDQsay9vZ3CwkL279/f38MmlFizqo7VzBaY2e2tra3pDkVE0qy8vJzm5uZOy5qbm5P6iu+eZFVi\ndff73P2SSCSS7lBEJM2i0Si1tbU0NTXR3t5OU1MTtbW1RKPRlJ9bjVcikpViDVRLly6lpaWF8vJy\n6urqUt5wBapjFRHpi9yrYxURGQyUWEVEkkyJVUQkyZRYRUSSTIlVRCTJlFhFRJJMiVVEJMmUWEVE\nkkyJVUQkyZRYRUSSLKsSq0a3EpHBIKsSq0a3So50vc5CJFtodCvpJJ2vsxDJFhrdSjpJ0essRLKF\nXs0ifZei11mIZAsNGyh9l87XWYhkCyVW6SSdr7MQyRZqvJJO0vk6C5FsoTpWEZHEqY5VRCQdEk6s\nZlZlZheHn8eZWVnqwuofPXklIoNBQonVzL4GfAX4arioALgzVUH1l568EpHBINES698BZwN7ANz9\nFWBkqoISEclkiSbWvR60cjmAmY1IXUgiIpkt0cT6EzO7DSg2s88DDwPfT11YIiKZK6F+rO5+k5md\nDuwCZgHXuftDKY1MRCRDJfyAQJhIlUxFRHqRUGI1s92E9avAEIJeAXvcfVSqAhMRyVju3qeJ4MmD\nRcD1fd031ROwALh99OjRsYa2TtO2bdvc3f1rX/ua1mu91mt9f9YnlIv6/Uirmf3e3f+2XzunmB5p\nFZEUSeiR1kSrAv5P3GweMBdo60dQIiJZL9HGqwVxn/cBLwELkx6NiEgWSLS71cWpDkREJFscMrGa\n2QqCittuufuypEckIpLheiuxZlQLkJktABbMmDEj3aGISA7TQNciIolLaq+AcQTDBs4GCmPL3f20\nfoUmIpLFEh2E5S6gBSgD/omgV8CGFMUkIpLREk2sY919JdDu7v/l7p8FVFoVEelGov1Y28N/XzWz\njwOvAGNSE5KISGZLNLH+i5lFgOXACmAUcEXKouon9QoQkcEgoV4BZjbO3XcMQDxJoV4BIpIiSX39\n9X+b2YNmVmtmow8jKBGRrJdQYnX39wHXAhXAE2a21swuTGlkIiIZKtESK+7+O3e/EjgBeAP4ccqi\nEhHJYAklVjMbZWaLzew/gMeAVwkS7KBiZgvM7PbW1tZ0hyIiOSzRxqsXgTXAT9x9XcqjOkxqvBKR\nFEneI63ANA8zsJlNdPe/9DssEZEsl2jjVXyx9v4UxSIikhUSbryKk1BROFnMbJqZrTSzexLYVnWs\nIpJ2/Ums3090QzO7w8y2m9mmLsvnm9mfzOx5M7v6UMdw983uXpvI+dz9Pne/JBKJJBqiiEjS9fYG\ngULgUmAG8Ayw0t3/rQ/H/xHwXWBV3DHzge8BpwNbgQ1mdi+QD3yzy/6fdfftfTifiEja9dZ49WOC\nAVgeBc4kGI/1HxI9uLv/xsymdll8AvC8u28GMLPVwEJ3/ybwiUSP3ZWZXQJcAlBaWtrfw4iIHLbe\nqgJmu/uF7n4b8Cngw0k45xRgS9z81nBZt8xsrJndCvytmX21p+3c/XZ3n+vuc8eNG5eEMEVE+qe3\nEmtsuEDcfZ/ZgLZbxc67k6A6QkQkI/SWWD9gZrvCzwYMC+eNoBfWqH6ccxtwZNx8SbhMRCQrHDKx\nunt+Cs65AZhpZmUECfU84PwUnEdEJC36090qYWbWAKwDZpnZVjOrdfd9wBLg1wTv0fqJuz+byjhE\nRAZSoo+09ou71/Sw/H70BJeIZKmUllhFRHKREquISJIpsYqIJJkSq4hIkimxiogkmRKriEiSKbGK\n5JiGhgYqKyvJz8+nsrKShoaGdIeUdVLaj1VEBpeGhgai0SgrV66kqqqK5uZmamuD4Y5rarrtdi79\nkNDLBDONXiYo0r3KykpWrFhBdXV1x7KmpiaWLl3Kpk2bDrGnhBIaiUqJVSSH5Ofn09bWRkFBQcey\n9vZ2CgsL2b9/fxojyxgJJVbVsYrkkPLycpqbmzsta25upry8PE0RZSclVpEcEo1Gqa2tpampifb2\ndpqamqitrSUajaY7tKyixiuRHBJroFq6dCktLS2Ul5dTV1enhqskUx2riEjiVMcqIpIOSqwiIkmm\nxCqSIfTEVOZQYpVBQUnj0GJPTK1YsYK2tjZWrFhBNBrVfRqs3D1rJmABcPuMGTNcMkd9fb2XlZV5\nY2Oj79271xsbG72srMzr6+vTHdqgUVFR4Y2NjZ2WNTY2ekVFRZoiylkJ5SL1CpC002OWvdMTU4OG\negVIZmhpaaGqqqrTsqqqKlpaWtIU0eCjJ6YyixKrpJ2SRu/0xFRmUWKVtOtL0sjVRq6amhrq6ur4\nzGc+w5gxY/jIRz5CNlbjZY1EK2MzaZozZ87hV1HLgKqvr/eKigrPy8vzioqKbhuucr2RK9evf5BI\nrCE90Q0zaVJizU653jKe69c/SKhXgGSXXG8Zz/XrHyTUK0CyS643cuX69WcSJVbJGLneMp7r159J\nNB6rZIxcH0s0168/k6iOVXJKQ0MDdXV1HYkpGo0qMUlfqI5VBP7a9zUvL4/ly5d3Gshk+fLlLF26\nNN0hSpZRYpWsE/8Qwfjx47nyyitZsWIFs2fP5q677qK6upqCggKqq6u56667uPPOOw/rQYNcfWhB\nDiHRflmZMKHRrXJe107006dP7+j7mZeX53v37u20/d69ezseSuiPJUuWeHFxcccxotGoOu1nNz0g\nILmnayf6+GR6qA72eXl5fT5XfX29T5o06aAnoaLRqBcXF3tJSckhnySTjKTEKrmna6k0Ppl290jo\npEmTPBqN9rnEWl9f76NHjz5koi4tLfVVq1bp0dPsosQquadrqbS+vt5LS0s7kmmsNBkrSZqZT5w4\n0evr6xMaryB2zLKysl6rFuIfN9Wjp1lDiVVyT3el0iOOOKJTMo0lzMbGRo9EIj5s2DA3Mx81apQv\nWrSo1wFOYnWpkUik2xJrcXGx19fXdyRZd+/0WTKaEqvkpviS59SpU72kpMTNzCdMmNAp4U6ZMsWL\ni4s7LZswYYIXFhZ6SUlJj1UEeXl5XlJS4mPGjOlUGm5sbPSJEyf6kiVL3N1VYs1OSqwi8ZYsWeKR\nSKSj5DpixIiDqg2mT5/uZuaRSMRHjRrlZnbQcUpKSnzUqFGel5fn48aN80gk4mbmxcXFHolEfNWq\nVR3HmTp1qkejUZ80aZLqWLODEqtIVyUlJR1Jz8w66khj3abMzEePHu15eXkeiUR8yJAhHQlxyZIl\nPmLEiINKvhMnTvQvfvGLHVULXddPmjSp2wQtGUmJVaSr+IanWKt+rNtUNBr18ePH+/Tp0z0vL8+n\nT5/ekVyPOeYYj0QiPdarRiKRjhJvd+tHjx49kJcpqaPEKtKdWFepRYsW+cSJEzuSYUlJyUF1pqWl\npT5ixAgfPny4T548uceeAGbWUa/bU08ByQq5l1j15JUkKlZKHT58eEeVQFFRUbelzaKioo6Sak8P\nGRQVFXWUcjXKf1bLvcQam1RilUSYWUdDVGNjY6c615hYaTS2rrvuXKWlpV5YWOjTp0/34uJinzJl\nit5Llb2UWEUOJVb6HDNmjE+YMKEjIcaLdcsaMWJER91rrAtXXl6ejxo1yo844ggfPny4Dx8+3Jcs\nWZLwgwaSkZRYRQ4lVvqMRqMdpdeuLfolJSU+fPjwg5aXlpZ6JBLxESNGeCQS8aKionRfjgwMJVaR\n3sRKl7H61WHDhvn48ePdzLyoqMiHDx/eUVUQr7Gx0UeOHOmFhYVeWFioUmnuSCgH6Q0CIgRjql59\n9dW0trYybNgw3nnnHcyCweJ37drF7NmzO71toL29naFDh1JUVMSBAwd4++230xm+DJyE3iCgd16J\n8Nf3SX35y19m27ZtAIwfP56Ghgaqqqpobm6mtra2Y9vm5mZGjhzJ7t270xazDF4qsYp0MXLkSPLy\n8lizZg3V1dUdy5uamli6dCkrVqygpqaG7du3M3nyZIqLi9m0aVMaI5YBlFCJNasSq5ktABbMmDHj\n888991y6w5EMdeSRR7Jt2zbee+89CgoKOpbHvv5HIhHa2toYOnQoQ4cO5Tvf+Y5eSJg7cu9lgu5+\nn7tfEolE0h2KZLAbb7yRoqIimpubOy1vbm4mEokEjRNmuLuSqnQrqxKrSDLU1NTwoQ99iJqaGpqa\nmmhvb6epqYmamhreffdd2tramDlzJq2trUqq0i0lVpFuPPDAA5xzzjksXLiQoUOHcvbZZ7Nr1y7c\nnc9//vM89dRT6Q5RBrGsqmONUeOViKRI7tWxiogMBlmVWM1sgZnd3tramu5QRCSHZVViVa8AERkM\nsiqxiogMBlnZeGVmO4CX+7DLEcDrKQonlTIx7kyMGRT3QBrMMb/u7vN72ygrE2tfmdnj7j433XH0\nVSbGnYkxg+IeSJkYc1eqChARSTIlVhGRJFNiDdye7gD6KRPjzsSYQXEPpEyMuRPVsYqIJJlKrCIi\nSZaTidXMxpjZQ2b2XPjv6G62KTSz35nZU2b2rJn9Uzpi7RJTInEfaWZNZvaHMO5/SEescfH0GnO4\n3R1mtt3M0jpitJnNN7M/mdnzZnZ1N+vNzP5/uP5pMzsuHXF2iam3mN9vZuvM7D0zuyodMXYngbgv\nCO/xM2b2mJl9IB1x9kuiL8fKpgm4Ebg6/Hw1cEM32xhQFH4uAH4LnJQBcU8Cjgs/jwT+DMwezDGH\n6+YBxwGb0hhrPvACMA0YAjzV9d4BZwH/Ef58nAT8Ns0/E4nEPB44HqgDrkpnvH2M+0PA6PDzmem+\n132ZcrLECiwEfhx+/jGwqOsGHoi9Ia4gnNJdIZ1I3K+6+5Ph591ACzBlwCI8WK8xA7j7b4A3Biqo\nHpwAPO/um919L7CaIP54C4FV4c/HeqDYzCYNdKBxeo3Z3be7+wagPR0B9iCRuB9z9zfD2fVAyQDH\n2G+5mlgnuPur4ee/ABO628jM8s1sI7AdeMjdfztQAfYgobhjzGwq8LcEpe106VPMaTYF2BI3v5WD\n/yglss1AGmzxJKqvcdcSfFPICFn7llYzexiY2M2qaPyMu7uZdVsSdff9wLFmVgz8wswq3T2ldYDJ\niDs8ThHwM+Byd9+V3CgPOldSYhbpjplVEyTWqnTHkqisTazu/tGe1pnZa2Y2yd1fDb/Gbe/lWG+Z\nWRMwH0hpYk1G3GZWQJBU73L3n6co1A7JvNdptg04Mm6+JFzW120G0mCLJ1EJxW1mxwA/AM50950D\nFNthy9WqgHuBxeHnxcAvu25gZuPCkipmNgw4HfjjgEXYvUTiNmAl0OLu3x7A2HrSa8yDyAZgppmV\nmdkQ4DyC+OPdC3wm7B1wEtAaV9WRDonEPBj1GreZlQI/Bz7t7n9OQ4z9l+7Ws3RMwFjgEeA54GFg\nTLh8MnB/+PkY4PfA0wSl1OsyJO4qgka2p4GN4XTWYI45nG8AXiVoYNkK1KYp3rMIelK8AETDZZcC\nl4afDfheuP4ZYO4g+LnoLeaJ4T3dBbwVfh6VAXH/AHgz7uf48XTHnOikJ69ERJIsV6sCRERSRolV\nRCTJlFgMZXT9AAACXUlEQVRFRJJMiVVEJMmUWEVEkkyJVQYtM9tvZhvjpqvD5R8OR+7aaGbDzOxb\n4fy3+nGOa/q4/T+bWY8PRIiABrqWQczM3nb3om6W3wo0u/ud4XwrQf/Y/ck6h8jhUIlVMoqZfQ74\ne+AbZnaXmd0LFAFPmNm54RNzPzOzDeF0crhfkZn9MBzb82kz+6SZXQ8MC0u+d3U5T76Z/cjMNoX7\nXBEu/5GZfcrM5saVpJ+JjYFgZtPN7AEze8LMHjWz9w/oDZJBIWvHCpCsMCwcXSzmm+7+AzOrAta6\n+z3QUeo8NvxcD/w/d28OH4n8NVAO/CPB46dHh9uNdvefmdmS2L5dHAtMcffKcPvi+JXu/ni4DWEV\nxAPhqtsJnhx6zsxOBP4NOC0J90IyiBKrDGbv9pD0DuWjwOxgyAQARoUjfX2U4Hl0APyv43z2ZDMw\nzcxWAL8CHuxuIzM7l2CA7jPC83wI+Gnc+Yf2MX7JAkqskm3yCN700Ba/MC7RJcTd3wxfBfIxgufX\n/x74bJdjVgJfB+a5+34zywPe6scfA8kyqmOVbPMgsDQ2Y2axJPcQcFnc8ti7t9rDYRY7MbMjgDx3\n/xlwLUGpNH59McHAMZ9x9x0AHox7+6KZnRNuY5ZJ72mSpFFilcEs1rAUm65PYJ9lwNywgeoPBKVN\ngH8BRoeNUU8B1eHy24GnuzZeEYxm/59hHe+dwFe7rF8IHAV8PxZfuPwCoDY8x7Mc/GoXyQHqbiUi\nkmQqsYqIJJkSq4hIkimxiogkmRKriEiSKbGKiCSZEquISJIpsYqIJJkSq4hIkv0vc+EIFkMg6VAA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e163a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(5, 3)\n",
    "plt.yscale('log')\n",
    "plt.scatter(beta.mean(0), np.maximum(pval, 1e-13), c='w', edgecolor='k', label='Genes');\n",
    "plt.axhline(0.05, lw=1, ls='--', c='k');\n",
    "plt.gca().invert_yaxis();\n",
    "plt.ylabel('P-value')\n",
    "plt.xlabel('Effect size');\n",
    "plt.legend(scatterpoints=3);\n",
    "plt.title('Sparse Reverse DE P-values');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Three of the genes have particularly small P-values, in order: mt-Nd2, Dppa5a, Ckb.\n",
    "\n",
    "I'm not really expecting very relevant results from this analysis, because the noise models are very crude, and I haven't corrected for any technical factors. But Dppa5a is a well known mESC pluripotency marker, and Ckb is known to be highly abundant in brain (NSC's are neural stem cells). While not very scientific, it's fun that it \"makes sense\".\n",
    "\n",
    "I just wanted to explore Bayesian thinking in differential expression, and give some small Stan examples on how to investigate small conceptual ideas of this.\n",
    "\n",
    "This post is available as a notebook [here](https://github.com/vals/ReverseDE/blob/master/Reverse%20DE.ipynb), with all analysis and code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>beta</th>\n",
       "      <th>pval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000064345</th>\n",
       "      <td>-0.288968</td>\n",
       "      <td>0.0080</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000060461</th>\n",
       "      <td>0.245447</td>\n",
       "      <td>0.0195</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000001270</th>\n",
       "      <td>-0.147287</td>\n",
       "      <td>0.0335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000057666</th>\n",
       "      <td>0.115933</td>\n",
       "      <td>0.1545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000064358</th>\n",
       "      <td>0.077362</td>\n",
       "      <td>0.2100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000008682</th>\n",
       "      <td>-0.084028</td>\n",
       "      <td>0.2115</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000052551</th>\n",
       "      <td>0.078640</td>\n",
       "      <td>0.2435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000060126</th>\n",
       "      <td>-0.061452</td>\n",
       "      <td>0.2490</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000049421</th>\n",
       "      <td>-0.025748</td>\n",
       "      <td>0.2990</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ENSMUSG00000029580</th>\n",
       "      <td>-0.033113</td>\n",
       "      <td>0.3105</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        beta    pval\n",
       "ENSMUSG00000064345 -0.288968  0.0080\n",
       "ENSMUSG00000060461  0.245447  0.0195\n",
       "ENSMUSG00000001270 -0.147287  0.0335\n",
       "ENSMUSG00000057666  0.115933  0.1545\n",
       "ENSMUSG00000064358  0.077362  0.2100\n",
       "ENSMUSG00000008682 -0.084028  0.2115\n",
       "ENSMUSG00000052551  0.078640  0.2435\n",
       "ENSMUSG00000060126 -0.061452  0.2490\n",
       "ENSMUSG00000049421 -0.025748  0.2990\n",
       "ENSMUSG00000029580 -0.033113  0.3105"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({'pval': pval, 'beta': beta.mean(0)}, index=idx).sort_values('pval').head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ENSMUSG00000064345\n",
      "ENSMUSG00000060461\n",
      "ENSMUSG00000001270\n",
      "ENSMUSG00000057666\n",
      "ENSMUSG00000064358\n",
      "ENSMUSG00000008682\n",
      "ENSMUSG00000052551\n",
      "ENSMUSG00000060126\n",
      "ENSMUSG00000049421\n",
      "ENSMUSG00000029580\n"
     ]
    }
   ],
   "source": [
    "for g in _.index:\n",
    "    print(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}