Subtype analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import and set options\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import matplotlib.style\n",
    "import matplotlib as mpl\n",
    "mpl.style.use('classic')\n",
    "import rpy2\n",
    "%load_ext rpy2.ipython\n",
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import glob\n",
    "import os\n",
    "import scipy.stats as stats\n",
    "import matplotlib as mpl\n",
    "import utils\n",
    "from matplotlib_venn import venn2, venn2_circles\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['font.family'] = 'Arial'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cancer subtype analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_table('CHASMplus/results/subtype_specificity.txt')\n",
    "df['mutation type'] = 'missense'\n",
    "#df['gene'] = df['gene'] + ' (M)'\n",
    "#df2 = pd.read_table('CHASMplus/results/subtype_specificity_tsg.txt')\n",
    "#df2['mutation type'] = 'inactivating'\n",
    "#df2['gene'] = df2['gene'] + ' (I)'\n",
    "#df = pd.concat([df, df2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['significant'] = 'no'\n",
    "df.loc[df['q-value']<=0.1, 'significant'] = 'yes'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEUCAYAAABkhkJAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XdYFVf6B/DvS+9IL4KIoChFRbBG\nbNhii4kxvZieTTZlN5uyLTHZ32aTbLKbnt1siinGFhN77wVLLIAigiAdlN7rvff8/pi5eiWIoMDM\nXN7P8/Dcy7kzc9+Z4d6Xc+bMOSSEAGOMMaY0C6UDYIwxxgBOSIwxxlSCExJjjDFV4ITEGGNMFTgh\nMcYYUwVOSIxpGBHxZ5iZDf5jZppARPFEtI6ILhBRIxHlEtEGIpqudGw9gYiyiUgQ0SKTslEAjrZa\nTsg/k3o4RMZumJXSATB2LUS0GMCrAAiADkAlgAAAgQBmE9HvhBDvKxdhjyiC9HmtAwAiigVwCPxP\nJTMj/MfMVI2I5gF4DVIy+hsANyGEFwB/AOvlxd4mIm+FQuwRQoixQogAIcQqucgJ/PllZob/oJna\n/Ul+XCeEeFUIUQsAQogLAB4E8AqABQAajSsQUSwR7SGiBiIqIaKvicjL5PXFcrPWf4jofiI6JzcD\n7iWicNM3J6J7iegMETXJzWavEZGlyesRJk2J9USUTkR/JCJqa2eIaIb83ikmZc/JZaXG9eT3FUR0\nRP79UpOd3By322R9IdciTXkQ0TIiqiOiQiJ6pcNHnDGlCCH4h39U+QPAE4CQf+Z2cJ1wSM1aAkA1\ngAb5eTIAW3mZxXJZlvxYZfI+R022tcikvAxSc6EA8Ln8uj2AArmsAUCpyfIvXyU+WwC18jJ+ctlq\nk/WGyGVL5N9flX/Pln9fBGAsgBKTdfIB/F5ezlhWDqAeQJNJ2TSlzyn/8E97P1xDYmrWz+R5lvGJ\nXAPJb/Xze/nlVwE4AHgfQB8AHpBqE1EA7mi1/f4A5gkhXAH8VS4bSURucu+1v8tltwkhPAAEQ0oE\njxJREIAhkJoOiyE1JXoC+A2A7QD0be2QEKIJwC7516ny43iTRca3em1DG9s4BGChye8BQoh/tVos\nEYAXAB8AOXLZrLZiYkwtOCExrdCZPHcF0LfVj4v82iT58X4AuQDSAcTKZZNbbTNNCGG8DvWzSbkz\ngEGQkg0AfExE+ZA6EbhCup41Ud52FQBvAIeI6E0A5wHMEUK8286+bJQfpxLREHn9dLlsPBFFyPtU\nBOBkO9tpz3+FEHVCiEoAh+Uyl/ZWYExpnJCYmuVAamoCgDBjoRDiFSEECSEIl//7N3KXHz1wOVk5\ny2X+rZYtNXleb/LcwmQ7xvWM27IxlgnpetZMSF2vhwP4I4CtAIqI6N529muT/BgPIE5+/gmk5HYT\ngGly2UYhxPUOx19i8tx4fc2yrQUZUwtOSEy1hBBlkGolgHTtpC2tv2Qvyo+3mSQtJ/n5zFbLmta6\nWn/xXzR57mGyLWf5+VtyjIchJZZAAA8AWAUpmX1FRG3WSIQQeQBOQUpwj8vFuwHsBxACqXYHXK5J\ntbmZdl4DAMM1XmdMdTghMbX7P/lxPhG9TUROAEBEQUS0FtL9SKYOyI/PEpETETkDOEFEZUR0Tyfe\nNxtSZwEAeIUkkQDKiSiPiAYR0UIiqgBwFgCEEN/h8nUnG7TfRGZMNjGQOiCcBrBHLhsBqTPC9nbW\nv5RMiciFiPieQqZ5nJCYqgkhNkPqcCAAvAQpIZRA6uQwT16sBMBO+flbkL7MJ0FqkrsA6XpQPaTm\ntI6+rx6Xk+GLkJrTTgKwBnBaCJEOYAekm3T7AsgmomIAJ+R1dgsh8nF1prWfA3LT3F6Tsr1CiLp2\n1s82eZ4P4J3294gx9eOExFRPCPF/kDoRrICUYFwgdcPeBeA5ACFCiP3yskmQmtD2QKpFNAFYC2Cy\n3ATYmff9L4BHIDWv2UBKfB9B7uEmhKiQ4/rGJK5sSD38brvG5g9BqhkBwD758SSkxAe00buuVWwF\nkBKm8TpYdQd2iTFVo+u/ZsoYY4x1Ha4hMcYYUwVOSIwxxlSBExJjjDFV4ITEGGNMFTghMQaeedVc\n8XnVFj5ZzOwQUX+TmVP7X2PZQUS0DVcO5Gr25Ok52pq2or11FsnrZHdfZF2HiGYA2KJ0HKzj+O5u\n1msRkR8u32PU25RAmjrDLO9fIqLbIE3r0XqsQ6ZinJBYb2aL3pmMIIRYeO2lNI1HNtcgbrLTGCJ6\nhojOyjOY5hHRh/J4bcbXA4hoBRFdlJfJJaL3iMhWft3YnNVIRKFEtE2e6TSbiB5v9V4+8myrJfLM\no4mtx4Oja8/OukR+v3fk96omoo+usm/XnH3VdOZUkzLjDLB72thsBBEdlPc3hYjmG48DTOZYApAl\nx7pR3taXrWI7Kpe/KP9ubBIcTUSb5XgzWo/yLR/DpURUSUS18jGIbmv/W633PBFlyufwLBH9ttXr\nxuMwlYg+l7dfTkQfEJG1cR9NzvVDJM0cW0BEA9tqsuvI+TZZdgoRJcnb/oWIRpm8ZjznLxPR3+Xt\nVRDRm0RkT0QfE1GVXP63Vuc3lIjWy+9fSUSriWiAyevGZsMtRDSTiJLlGI4T0U3GZQB8La8S1Prv\nhamY0jME8k/HfwD8DZdn/6yCNAmcALDZZJnjcpkO0tA0xuX/LL/eX/69BUAepNlL9SbLRcrLOUIa\nNFRAGjm60mSZu+RlOjI76xK5rAlAs7z8I23sW4dmX4XJzKkmZYvlsj2t9lHI79cgv7eQ93UkpEFZ\ni0yWKwLwL0iT+BlniLWStxcgHwMDgAC5zLjeBflYNpiU3WqyT2fkskaTY1gDYFA75/k1k+NeKj8K\nAH9q4zhkyftm+v5PtzoOBjnGSkhNWARpaCUBYHEnzvcik2Pa2GqZXACWrc55MaS/w1qT5dLkc2Ba\ndoe8no98PI3vUSM/L4Q04rppDPnyfpvO9lsAaazBhbj8t6+Tl12o9OeXf679wzUkjSAid0iDfALA\n74Q0y2kspA9cHEkX530hfXiPAfAVQrjj8qCbo1tt0grSdAfukGY+bZDLjXPxPAhpDqIySF+efQC8\nJ7/2sPzYmdlZCUAEpFlMl7exi52efbWDdgNwg/RllwipVeAvQhr4dKzJcmOFEL+HNO5dJaTjMkV+\n7TY5/r3i1wOmFsrbdgewTS57TX58UN6vBEjTsbsB+AcAJwAvtxUsEfUB8AqkpDBKPg7DICWUV4jI\nodUqDZASpiek62EAML31ZgF8J5/DWCF/s7fSkfNt5ABpyvQ+uDxmXyCAyDaWi4R0bNLksn6Q5o7y\ngDxKOi7Pkvs7SMfyR0jHyg3A9wD8ADzVatt95RhcIU37AUh/P5FCiFUAjDMI5wtpRt1VbewzUxlO\nSNoxBtI1j0YAHwKAEOIkgFAALkKIdCHEBSHEXEhftCFE9DtIA40C0pdgax8KIZqFNHJ1ilxmbP6b\nJD+uFkJkyM8XA/AWQkxvtUxHZmc9IoQ4J4SoF22PYn29s69ey9tCiEYhDYT6oVw27GoLC2mK8RXy\nr7fLj8Yv3aVtrPK+EKJcCNEAaaRxABhK0jQZk4y/Q/ryzQPwmFzW+vgYjQFgJz9fQ9JMtZshfVad\ncfn4Gn0rhCiWj6lxNHNn/NpSef9K2ngN6Nj5NqoH8B/5+WaT8tbXbfYKIc4KIZpxeRT0vUKIU/Jx\nTmwVrzGGqZDOfTaA2XJZ6+PVAGlSQwD4yaS8rX1nGsEJSTuMM5hWCCEuTb4mhMgx/Z2I/gypB9Vh\nAM+brN/WuW5rxlTjcsb3uzRCthCittUXWmdmZy1q4/0vEdc/+yrQfucc04n2CuVHv2ts7xv5cT5J\nPfHGQ2oe+rGD2ydIU50bj48TLh8fT7ms9fExMq5jgSunaDdORNiRWW/bOtftHn907HwblRn/5oQQ\njSblrSdLLDd53iw/mm6vSX40XkMyxtAHl/fbTS5rvd/lJjW91rP9Mo3ik6cdxg+3JxFd6hlG0iRx\n04moDxHNgjQlQTOAcCFEEIDP29lmezOmGt/v0hcBEbnLF8dHyheiOzM7awOuQXRs9lVj8rU1WdW1\nnc0GmTz3lR8rjG95lTgOQaqxeQH4N6Qv2o1CiMo2Fu9v8tyY6IT8Hsbj86HJ8bGHdK3FDm0zrlNt\nXKfVcW3d3NneOTR1rePfkfNt1NHZaNtqam2v+dW477832W9HABZCiMGtlr2031dpguRpDDSIE5J2\nHIb0H6U1gD8AABFFAfgOUk0iHJfb8JsBFMhf4nfLZZ0913vkx1uJaIj8/AUAXwFYKn8JdGZ21na/\nIKjjs68ak8k4eT0PAHPa2fRfSZpR1RXAM3LZQfnR9Mu89ayr38qPd8qPbTXXAcDzcu80G1y+bnFC\nCFGPy8fnTiIaQNKoAV8CqCWiD9vaGKSmrXo5nqflfZwBoIaI0uRrTNfjWl/Qe+TH9s53dzMer0eJ\nyFs+ppsBVJHcu7ETjOfWkYgsiGfU1QROSBohhCjH5Q4KfyeiKkht8LYAdgghEiAlLeByD7JiXG6X\n7+x9GUsgXYh2AZBCRJUA/iS/9ob82CWzs8o6OvuqcWbYB4goE0Amrt5kp4N0ja0U0n/fMZA6B/xD\nfr0EUm8uQOp4YFr7+A6Xv8SrcOUMr6YCIPXuqgQwV17ndfm1byFdC/EBcA5SLeQeOd6VbW2s1bWu\nj+XjvhlSs9auq9TSusISXPt8d7cPIR2jcEjHtBTABEjHdE0nt5UtP3pC+ifm6a4JkXUnTkgaIoR4\nFdIMqemQLnznA/gA8sV3IcQ+SLWAXEgf4tOQmr4MACKJyKcT71UPaTbUbyF9SVhDmtH0XiHE9/Iy\nXTk7a0dnX30D0pdnDaRrM58B+PNVNquH1FPOWCM6DamTxDH5PZsh9VwsgvSFb5ytFUKIXFzuGba6\n1bUSU49DGp7GAkAGpOOzXt5GA6Rk/QOkL0VrAEcAzBZCHGhza5I/Q+qFdw5SE1+uvN/PtLPODenI\n+e5uQogiAHGQZss1XhfaCSBeCHGuk5s7BOALSOdUhyuvMzGV4hljGWsDEU2C1GUcAKYIIXa3et34\nwZkshNjTg6ExZra4hsSYCSKaJXe13iEXJeHy9RXGWDfihMTYlfIg9exqgHRT7vweuqDPWK/HTXaM\nMcZUgWtIjDHGVIETEmOMMVXghMQYY0wV+O5lxtpARJbuLi6B4cEDhnr0cXV1cXRycrS3d3Cws3Ow\nt7V1tLWxcbC1tnawtLR0sLSwsDAYDEII6M9mZ/epytVf1AtDk96gb9Ib9I16g76pSddcV1FfWVxc\nU3q2oKrorBCi6tpRMNa7cKcG1msRkfNNw4bHhQYE9Pfx8PBzc3Hxc3Zw9HWws/Wzt7Xz7efr4+Xn\n6WVtbdXx/9t2HDlSVn7awqOt14QQaGhpRHldRU1tU93FFoPuYrOu+WKjrulifXP9xaqG6osltWV5\nuRUFh5pamsrb2gZj5owTEusViMghLjp6wvBBYdGBPj6D3Vxcwnzc3QcN7h/sZm9re+0NdNC2Q4eL\nK89Yel/v+i36FhTXlFZX1Fem1zbVpZTWlZ/JKss9nF2WmyCE0F17C4xpFyckZpaC/fsGTRgxYn54\ncPBQd1fXwZ593AZFDhjgaW93tUG2u8bWhEMXq1KtOjxEU0c0NDeI/MqivNqm2tSqxpozRVUXU9KK\nM7ZUN9QUdOX7MKY0TkjMLBCR3Zy4uFtjh4Tf5O/lFTOoX9DQAX37Olw5a0L323TgYFFtms215lu6\nIQZhQElNWX1JbWl6ZUPV4YySrK2nClM3cA2KaR0nJKZZoYGBPpNjR94XHhw8NsDbZ/SIIYMDbK1t\nrr1iN1q/b39Bwzm7vj35nk26ZmSV5aSX1VUcyCrN2Xky/9RqeUZWxjSFExLTFCKyu3Pa9PtGhkfM\njAgJmRwxYIB7T9eC2rNuz978xkyHAKXeX6fXIassN6e0rnx/dlnunmO5icvkkbwZUz1OSEz1iMhi\n2ugxN08YMWL+wMDAqaMiIvt3pudbT/p51+7cliynfkrHAQB6gx455XmFJbXlBzJLstYdz0taZjrd\nPWNqwwmJqdbQgQODZowZ+1RY//7Txw0dNszF0VE9VaGrWL1jV44+xzno2kv2rGZdM9KKM05kl+X9\nlJh/+uOK+kq+D4qpDickpjpTRo4aOzkm9vFRkRFzIwaEtHlPj1r9uG1ntiHPpb/ScVyNEAJZZbn5\nuRX5a89eOPd5enFmstIxMWbECYmpAhHRvAkTb58wYsT9ccOjp/l7eXVv/+xusmLLtvNU4DZA6Tg6\norS2rDazNHtH2sXM75MKTv/E02wwpXFCYooiItsH58x9YnRExJ2TYmPHOjuov1muPcs2bT1vWeSu\niYRk1NjSKM5ezDiSXZb70/7Mwx+1M107Y92KExJTBBFZ3Dvz5icmx458In7UqGFWlpZKh9Qlvl+/\nOdOm2DNE6Tiuh0EYcKYoPTn1Qtonh7OPf8EdIFhP44TEetyCKfELJsXGPj9z7LjxXTlsjxp8u3ZT\nhl2pV6jScdyIZl2zSC48s+/shXP/PJGXvFHpeFjvwQmJ9ZhZN42fMDl25Eszx42d4ebsos5+2zdo\nyc8bzzmUew9UOo6uUNNY23yq8Mz6lKK0N1MvpJ9QOh5m/jghsW43IXpE+LTRY/40fcyYW/t6ezso\nHU93+mr1+nNOlb5mkZCMimtKq04Xpa5Iyj/9Rl5FIY+fx7oNJyTWbYjI5umFd7x6y8RJjw8JDvZS\nOp6e8L9Va8+5VvubVUIyyinPz08pOvvNzrR9b/LoD6w7cEJi3WLBlPiZ08eMWTxz7LjRlmbSYaEj\n/rt8zTm3ur5mmZAAqeNDSlFa0qmCM68cy03conQ8zLyYZTs+U06gr6/LPTNmvvvKokX3+Xl62Ssd\nT0/TG8y7Y5oFWSDKf8iwvq6+q6YNnvTVjrS9L3M3cdZVuIbEusxd02csmD0+bvGkmJhINQ142pM+\nWbr6nFdjP7OtIZkyCAOS8lOOnCo884fE/NMHlI6HaR/XkNgN83Z3t180Z+6HLz246D5vNzdNjrDQ\nVQx6864hmbIgC0QHRo3u28dv3ZSwuE93px9YzHMysRvBNSR2Q2aOHRc9b+LET+fGTRjTW2tFpj5Y\nsjLDTx+s6fuQrofeoMfxvOQ9pwtTnz9dmJqkdDxMmyyUDoBp16K58x5//p57N86bMJGTkUxv0Csd\ngiIsLSwxKih60szwKdsmhI57mYj4u4V1Gv/RsE4jIvuXHnjwixfuve/jyJCQbp2uW2sMBtGrM7O/\nq6/37Mipb82LmvGju6ObS1dum4j+Q0SbWpW9RUTriCiAiNYQUR0R5crlNvIy1kT0GREVE1E9EW0n\norCujI11DU5IrFNmjB077OOXXt752zvufMTNxcVa6XjURq/v3QkJAKwtrTEhdNytsyOmbQnzCR3U\nhZv+AcBUInIzKVsIYBmAnwFUAIgBcC+AOQD+IS/zWwDTAcwGMBRADYAlXRgX6yLcqYF12P2zZt/9\n/N33vhsVGuqvdCxqZTDzbt8dRUSIDowa62LvvCE2aPhvjuUk7uyCze4HcAHAfABfE9FIAL4AqgEM\nADBGCKEHcJaIngawjYheBtAfQAOAbCFECRH9BkCv6AmpNVxDYh3y6PxbX3z6jjs/52TUPmHgz5Sp\nEM/+A+MHTVgeFzLm8Rvdljxf03IAd8hFdwJYByAYQB8AVURUS0S1ADYBsAEQBOBjAC4AColoN4C7\nAZy60XhY1+MPD2sXEVn89s673n/+7nveDPD2dlI6HrXTGwy9vsmuNR8XL8/4wRM+jA+Le7sLOjv8\nACBebra7HVJznRWAcwCGm/wMg1QLyhNCnAMwCFIiywTwZwCHiKjX3bitdtxkx66KiOxeeXDRt4/e\neutCW2sbpcPRBANfQ2qTi52z7fQhk1+yt7bv72zn9FBNY+11jYUnhEgkogwALwBwBbAFQDyAQABl\nQogKACCi8QCeA3A/ET0AoFkIsRzAz0S0GEAepOtJR25031jX4YTE2hQVGurz1jPPLr17+ox4Cwuu\nSHeUwaDnhHQV1pbWmDxo/B3Odk5+IV79H8wsyc66zk0tA/AnAEuFEM1EtA3AeQBLieiPABwAfAEg\nSQjRSESuAP5CROUA0gHcD6BWfs5UhL9p2K9MHzM24rm77tl878ybORl1koGvIbWLiDAyKDouPmzC\n+gi/sMjr3MxyAHaQEhPkjgzzAOgBJABYD6kDxKPy8p8A+Fr+OQvgFgBzjLUpph5cQ2JXiB85KvyR\nW+avnDBiRLjSsWiRMAhOSB0w2GdghCVZrIjwG3x3StHZ5E6u7gegCMAeY4EQIgvA3LYWlqdif0X+\nYSrGHx52yeTYkYMfmjdvBSej62fghNRhA71DwicMHLsi0n/IsI4sT0Q+RHQ7gHcBfCXXjJgZ4Q8P\nAwBMjIkZ9PC8eSsnxcRebzMKAyAEd2rojIFeAwZPCB27Isp/SHQHFneB1OxWD+Dt7o2MKYGb7Bgm\njogJeWjevFWTY0dGKR2L1ukNht4zG2EXCfUKDgPEDxF+YQtTitJOX205ufu2cw+GxnoY15B6ubjo\n6AEPzpm7auqo0UOVjsUc8I2x1yfUa8DgcQNG/TDQe8AApWNhyuEPTy82OjKq34Oz566aPmZMR5pL\n2DXo9XqQ4FGur9cQ30FR40PGrBzgGdRX6ViYMvjD00t5uLo63j1jxvczx40boXQs5qJZp4MlWfE1\npBsQ5T8k5qYBo1f6u/p6Kh0L63mckHohIrJ8auEd3y+YEh+ndCzmRKfXwcrCkq8h3aDowKhxYweM\n/IaI+Bp3L8MnvJchInru7ns+fnjeLfN5Ur2updPpYUmW/E9eFxjTP2ZWdUPNvwA8q3QsrOfwh6eX\neXT+rS8+futtj1lb8f8iXa1Zp4M18YHtCpYWlogLHfNEV4wSzrSDE1IvcveMmfMXzZ37FxdHR25W\n6gY6nQ6WFlY8aWEXcbJ1tBkTHPt/0YFR3LTcS3BC6iVuHnfT0HtmzvwgyNeP7+PoJk0tzTobSx4V\nvSv5ufp4jQoa8Wk/9wCeh6sX4ITUCxCR/Zy4CZ9Hhw3up3Qs5qy5Rae3tuQKUlcL8wmNHBU04msi\n4oNr5ri9uxd44b77P5gbFzda6TjMXUtLi96aW+y6xej+I6ZXN9a8D+BppWNh3YdrSGbu/lmzH7zv\n5lkP8jQS3a+xpdlgY8EJqTtYWlgiLmTMY3GhY55SOhbWffhbyoxNjIkZdHv81Dc8+/ThCxs9oLm5\nxWDDNaRu42jrYD0ueNTrwwMixysdC+senJDMFBFZzRkf90nMkCF83aiHNLW0GKwsuBW8O/m4eHlG\n+A1+i68nmSdOSGbqt3fc+ffb46dOVTqO3sQgDAZuGu1+wwMib5o8aPxipeNgXY8/PWZo4dSps+6d\nefNTfPNrzzIYeMK4nmBpYYmR/aKfGNo3YqTSsbCuxQnJzBCR3Ywx4/4W4OPjpHQsvY0QQigdQ2/h\n4+LlEeEX9jYR8U3eZoQTkpl59q67F08fM4ZH8FaAXm/gGlIPig6Mmjxp4E1/VDoO1nU4IZmRaaPH\nDJ8/cdKjfB1DGQaDgWtIPcjKwgox/YY9M9h3IE8uaSb4m8tMEBFNGz36zYH9+nkoHUtvZTAIg9Ix\n9Db+rr7eQ/3D/0nEEyOaAz6JZmLRnLlPzp80eabScfRmBsE1JCXE9Bs2fULo2N8rHQe7cZyQzEBU\naKjPnLi4Fxzs7HiCIwUJwTUkJVhbWmNE4NDfDfAMGqB0LOzGcEIyA7dMnPT30ZFRIUrH0dvpdVxD\nUkqgW1//SP8h3MFB4zghaVz8yFGjZ9100x1Kx8EAPTfZKSrcd9Dt3MFB2zghadzUUaNf6MdzHKmC\nQc8JSUnezl59wn0Hvax0HOz68a38GjZnfFz86088OVfpOJjEYBCip+/SLK0tx9rkTcgqy4WNpQ2G\nB0Ti5oh4mM7LpDPo8K+dn2FY3wjMCJ/yq22U11Xg71v/3eb27465DbFBw5FccAY/JW6AhYUFFgyf\ngwi/wZeWWXJ4OeJCRiPEK7jrd7CTwn3DbonyDx93qvBMgtKxsM7jhKRhk2Jjf+fj4WGndBxM0tMD\nNegMOnx1aCl8nL3wzMTHUNtUixUn1gAA5g293OFyx9m9uFhTctXt9HFwxWuzXryibHvqHqQVZyDC\nfzAMwoBVJ9dh/tCbISCw/PgavDH7ZRARCiqL0NjSqIpkBADujm6Og31D/wDgNqVjYZ3HTXYadevk\nybdMHzN2utJxsMv0+p7NSLnlBSitLcddsbfCx8ULIV7BmBkejxN5yZeWKay6gCPZJ+Dt7HnV7ViQ\nBVzsnC/91DTW4nD2cdwVcyvsre1Q11SP+uZ6DAuIwPC+kahvrkdtUx0AYGvqbkwfMqm7d7VTIvwG\nz4oOjJqhdBys8zghaRAR0aQRsc+4u7jwEPwq0tP3IXk7e+LRm+6DrZXtFeU6g84YD1YcX4PZkdPg\nYOPQ4e1uOL0NUf7hGODZHwDgaOsAG0trFFQWIa+yEDaWNnC0dUBeRSFa9C2XllMLFztn2zDv0N8R\nEd8GoTGckDRoYfzUe2aMHTtZ6TjYlfS6nh3KzsnWEYO8L/f2NwgDDmYeQbBHEABgT/pBONo4ILbf\n8A5vM6+iAOeKz19R67EgC8yOnI6P936Jz/Z9jXlDZ8KCLLA1dZfqakdGkf5Dpsb0G8bNdhrD15A0\nhojowz+8+KiTgwP/M6EyBoUH+16XvAUFVRfw/OQnUFJTit3nDuJ3k5/o1DYSzv+CMJ8Q+Lp4X1E+\nPmQ0YvsNBxFga2WL3PJ86A169HMLwIrja6TrTX6DceuwWbBQwSg+Djb2lmHeoc8S0U88Crt2KP+X\nwzpl9vjxN0+KjY1TOg72awaDMgM1CCGwJmkTDp4/ivtG3g4fZy+sOLEGU8MmwN3RrcPbMQgDThWm\nIuYqNSo7a9tLzYNbU3dhZvgUJBWk4GJNMV6Z9izyKwqQVJDSJfvUFSL9B8dFB0QtUDoO1nFcQ9KY\niSNiHnJzduE5YFRIr+/5f8QNwoCVx9fiRF4y7h91ByL9h6C8vhJZZbkoqCzCljO7AAAt+hbklRcg\ntyIfj930QJvbyi7LQ7O+GeEgqIv4AAAbfElEQVS+g9p9z+yyXABAkHsgTuQlY4Bnf9hY2SDEKxi5\n5fmIDojq2p28TrZWtjTQe8AdAH5UOhbWMZyQNGT88Oihb/32GR5AVaWUmH1i/amtOJGfjEVj7kK4\nXxgAwNXOGX+c/twVy313dBX6ewRiyqCrV65zyvPQ19UPdtbt30mwJXUXbg6PBwAQ6FJ3d4Mw9HjX\n92vp7x44rZ97wIDc8vzzSsfCro2b7DQkfuTIx/v7+/NMsCrV0012OeV52JdxCDOGTEaAmz+qG2tQ\n3ViDuuZ6eDp5XPFjZWkFB2t7uNq7AJBqTNWNNVds70J18a+uHbV2vjQbVmSJIPdAAECAmz9SL6Tj\nQnUxzl44h37uAd2zs9fJx8W7T4Rf2GNKx8E6hmtIGuHm4uL0vz//lUdlUDGD3tCj3YyN12s2pezA\nppQdV7z2zvzXYGlx9Zbdk/mnseL4z3jvtjculdU21cLXxafd99yauhuzIy7f/hYdEIXUC+fw0Z7/\nIco/HMMDIq9nV7pV3z7+c4nor0IIndKxsPaR2qrYrG2Pzp//4quPPv4OzwarXu99sSwzkEJ51HWV\n0el1+DFx/aKj2Se+UToW1j7+dtMAIqIRg4fcwslI3fSGnq0hsY6xsrRCiGcw97bTAP6G04Cbhg0f\nMzZq6Cil42Dt4+n51CvYo9+UUK9gdXT/Y1fFCUkDJsXELPByc+NhglSOa0jq5enk7jjEd9DDSsfB\n2scJSeWIiIL79o1XOg52baKHOzWwzvF39Z1LRLbXXpIphROSyk0ZOXLamMgongVTA7iGpG6hXsEh\nMYHD7lY6DnZ1nJBULm549C2uTk58njTAYBB8nlTM0sIS/T0CxysdB7s6/gCpGBFZBfft++spPpkq\nCU5IqtfH3nUskQpGf2Vt4hOjYrNuGj9/XNTQwddekqmB3iC4yU7lgj36DR7kHTJR6ThY2zghqdi4\noUNn2tvxDOVawTUk9bO3sbcY6DWAZ5NVKf4AqRQRkZ+n10il42AdpxcGHoVdAzyc3McpHQNrGyck\nlYoMCR0aFRo6ROk4WMcJPbjJTgP8XXxG+Lv6DlA6DvZrnJBUKi56+Ew/T0++GVZDuJedNng5ezqG\n+YTyUEIqxB8glQoN7BejdAys4wwGA0jw50kr/Fx9uNlOhfgDpEJEZOXj7s4JSUNa9DpYkCV/njTC\nw8FtNBE5Kh0HuxJ/gFRo/PDhk0cMHsxt3Bqi0+lhZcEJSSv6uQf4RQdE3aJ0HOxK/AFSobFRQyc6\nO/A/b1rSrNPBkqy5l51GWFpYwt/VN0LpONiVOCGpUKCvb7TSMbDO0el0sCJLnoFZQ1zsnQcpHQO7\nEicklSEi6uPEHxStadHpYMU1JE1xsLbnz5nKcEJSmSA/v9DQwMAgpeNgndPc0mKwtrTiz5OGuDu6\nhTraOvgoHQe7jD9AKhM7JHxCgLc333+kMc0tLTprSz5tWuLt7Okw0GsAzzWmIpyQVGZQUNBAIr7h\nX2uaW3R6G05ImmJBFvB39ePBi1WEE5LK+Li5hyodA+u85pZmvbUVJySt6WPvwteRVIQTksq4ODkO\nVDoG1nmNLc0GGysbpcNgnWRvwx0b1IQTkoo42tt79PXy5htiNai5uUXYWHANSWvcHFxDba1t3ZWO\ng0k4IanIuKHDJg/s189J6ThY5+kNBr2FBX+ctMbH2cs5zDuUOzaoBH+CVCQkIKCftRXfW6lFQgiD\n0jGwzrO0sIS3s2ew0nEwCSckFfF2d/dTOgZ2ffQGg1A6BnZ9HGzsPZSOgUk4IamIq5OTr9IxsOtj\nMBj0SsfAro+dla2n0jEwCSckFXGyt+eEpFEGg+AakkbZWNlwDUklOCGpiK2NjbfSMbDro+cakmZZ\nW1hxQlIJTkgqQUQWNlbW/MHQKgGuIWmUlSUnJLXghKQeHj4e7nw/hEbpDXpOSBplbWntTkTcvVUF\nOCGpRFRoaKiPu4e90nGw66PX8zUkrXK1d3G3srDiUb9VgBOSSrg5uwS4OPIssVplENztW6tc7Jyt\nA938eZBVFeCEpBLOjo52fKe/dhn0Br4xVqNsrWzg5tAnUOk4GCck1XCyt7dVOgZ2/fi+WG1ztXPm\njg0qwAlJJezt7OyUjoFdPwNnJE2zsbLh67cqwAlJJextbbmGpGEGHstO04jIUukYGCck1bC1seEa\nkobp9Xqe5lfDLMiCE5IKcEJSCVtra64haRgPHaRtFmTB9yGpACcklbCxsuYakoYZuJOdpnGTnTrw\nfwXqwf8caNiQ0CDnvAsl55WOg10f+wq4KR0D44SkGjqDvknpGNj1Gx8dzSO1a1hNY+0qpWNg/F+5\narS0tDQrHQNjvZVer+c2VxXghKQSzS0tXENiTCHcbV8dOCGpRBMnJMYUo9frdUrHwDghqUZjUxMn\nJMYUUtfYUKN0DIwTkmo0cEJiTDEV1TVVSsfAOCGpRm1DAyckxhRgMBhwsaysWuk4GCck1aipr6vR\n6fVKh8FYr1NdV4eSyoo8peNgnJBUIyMvL7W4vJx7+jDWwypqqltOZWTkKh0H44SkGkWlpbkFxcUV\nSsfBWG9TWllZA6BM6TgYJyTVEEI01jTUlygdB2O9TX1DY7UQolHpOBgnJFVpam7mhMRYD2vR67hD\ng0pwQlKRxubmYqVjYKy3qW9sLFI6BibhhKQidQ0NXENirIdV19ZyDzuV4ISkIpU1NZyQGOthJZWV\n+UrHwCSckFTkQlkZN9kx1oOEEMguKuQu3yrBCUlFUrPOp9bW1ysdBmO9RllVlUjJzDypdBxMwglJ\nRY6mpBxOy8kpVzoOxnqLjLy84rScnDNKx8EknJBURAhRV1ZVmaV0HIz1FtV1tflCCJ4cUyU4IalM\nVW0tJyTGekhtQwN3aFARTkgqU1xRcU7pGBjrLUr486YqnJBUJjXrfJrBwGOsMtbd9Ho9TmWcS1Y6\nDnYZJySVSUhO3plVWNCgdByMmbuzOdnV248c2aR0HOwyTkgqU1xenp9VWHhe6TgYM3d5Fy+m1jU0\n8CjfKsIJSYVKKyu5GYGxblZaUcndvVWGE5IKpefkHFc6BsbMXWZBforSMbArcUJSoYNJSZvyiy82\nKR0HY+aqtLJSfzTl9F6l42BX4oSkQinnM1NTzp9PVToOxszVqYyMrORz504oHQe7EicklSouK+dm\nO8a6ycXyskQhBN9foTKckFQq5XzmcSGE0mEwZnaEEDibnX1Y6TjYr3FCUqldx375OS0np1bpOBgz\nN5n5+XX7Tp5YrXQc7Nc4IalUYUnJhcz8PO7+zVgXS8vJPpGek5OtdBzs1zghqVhmfsEBpWNgzNxk\n5hccUjoG1jZOSCq27+TxlYUlJdz9m7EuUlZZqT+YlLhe6ThY2zghqdjhU6eOHz+bekzpOBgzF7+k\nnklOSE46qHQcrG2ckFQuPTd3t9IxMGYuci9c2C+4+6pqcUJSuf0nTizLKSri0b8Zu0GllRUtB5MS\nuXedinFCUrljqWfOJJ1L53smGLtBCcnJh3cePbpP6TjY1XFC0oDUrKxdSsfAmJYJIZB87hzPfaRy\nnJA0YNexX5am5+RUKR0HY1qVlJ5+ceuhhP8pHQdrn5XSAbBrS8nMzHrn2ee3DQoKWqh0LKZWbt+O\nF97/V5uvHf76G3y4fBl+2LrlivLXHnscj86/9VfL6/R6vPn1V1izezfqmxoxOSYWbzz5G3i5uQEA\nfjmTguffexe19fV4+cFFuGfmzZfWfeOL/yE8OBi3x0/twr1j5uR0Zub2rMJCnoxP5TghacS+E8eX\nzomLW+Di6KiaWu3cCRMwKSbm0u8GIfDw64sR6OuDvt7eOJeXiz8//Ahumzzl0jLODg5tbutfS7/H\ntsOH8J8//QmuTs547b+f4bl3/4kf/v4mAOCvn32Ke2bcjGGDBmLR64sxc+w4uLu6ori8HAcTE/Hn\nhx7u3p1lmlVbXy8OJJ78Uek42LWp5suNtW/jwQPr9h4/pqrODfa2tvB2d7/0syXhIApKivH2M88B\nAM7l5WHYwEFXLGNvZ9fmtgwGA15//EmMiohEWFAQHp43H0dOn7r0ekZ+Pm6+aRzGD4+Gq6Mjci5c\nAAB8smolnrhtASwtLbt/h5kmHUxKOrnx4AG+GVYDOCFphBBCHEtNXW0wqHPE/Nr6evz7h6V44b77\n0cfZGcXl5aisqcGAgIAOrf/KoocQP2oUAKCkogLLtm7BuGHDL73e18sLpzIykF98EVW1tfDz9MSF\nsjIcPnUKt0yc2C37xLRPCIFTmRmbeKoJbeCEpCFfrVv7n0OnTmUpHUdbvt+8CTbW1rh7xkwAwLm8\nXFhZWuLd775F7AP3Yfpvn8LK7duvuZ23v1mCEffdg1/OpODVRx+7VP7Koofwhw/ex/hHHsaTC26H\nr4cHPlqxHE/efjvXjthVnUxLK9x0YP9HSsfBOoYTkoYIIeqT0tPWKh1Ha0IILN2yGQ/NnQdrK+my\n5Lm8PABAeHAwvnv9b7hr+gz88eMPsWH//na3dcfUadjw/gcYN3QY7v3rn1FTXwcAuHncTUhetgJJ\ny1bgD/c/gMKSEhw7cwa3TJiIt75ZgtGLHsBTb/8DDY2N3buzTFMOn0r+MT03t1jpOFjHcELSmB2/\nHP1vZn5etdJxmDqVkYGcoiLcNiX+UtmDs+fgxNJleGjeLRgSHIyH592Ce2bejO82bWx3W8F9+2LY\nwEH44IU/oL6xEZsPJlx6zd7WFq5OTgCAj1Ysx1MLF+JkWhrW7d2DnZ/+B3q9Hl+vX9c9O8k0J+V8\nZumWQwnvKx0H6zhOSBpz9PTpswnJye1/q/ew3cePITosDL4eHpfKiAhuzs5XLBca2A8Xykp/tb7B\nYMDWQ4dQUlFxqczezg6BPj4or/717Vf5xRdxMu0s5sZNwPGzqRgeFgYnBweMHx6NE2lnu3DPmJYl\nJCWtOXH2rCqbuFnbOCFp0M5fjr6fd+GCamaTPXn2LMZERl1R9sYX/8Oixa9dUZZyPhMhAYG/Wt/C\nwgJ//c+n+HnP5XFkq+vqkF1YiIGB/X61/AfLluHpO+6EhYUFCARhkMbK1Ov14GEzGQCcL8iv3n70\nyAdKx8E6hxOSBu04cuTormPHVNONNS0nG4OCgq4omzZ6NHYfP4av1q1FdlEhlmxYj9U7d+DJBQsA\nAM0tLSguL4derwcALJozFx+vXIHdx47hbHY2nv3nOwj274vJsbFXbDenqAinMzMwZ3wcAGDowIE4\nfPoUzpw/jy2HEhAdFtYDe8zUbn9i4vqEpKTTSsfBOocTkkbtOnb0XzlFRTVKxwEAJZWV6ON0ZfPc\n2Kih+Pill7F0y2ZM/c2T+HbDenz80isYFREJADiemoqY++9FYanUhPfEbQvw8Lxb8PJHH2De75+H\nlZUlvnptMSwsrvwT/WD5Mvz2jjtBRACA0ZGRuGXiRNz+8ouwt7XDQ3Pn9cAeMzUrKC6u33X0KPes\n0yDiqUG0a/ETT3776C3z71c6DsbU5Ov165b/9bNP71Y6DtZ5XEPSsG2HDv3jbHY2j8/FmOxsdnbp\n1kMJrysdB7s+nJA0LCE5KXXH0SMrlI6DMTUQQmBLwsFvDyQmcldLjeKEpHFbEhL+fjz1TL7ScTCm\ntAOJiekrdmzn2pGGcULSuMT0tMJ1+/Z90tTSrHQojCmmqaUZ2w4f+jTvwgVV3TTOOocTkhn4at3a\nf67ft4+nZma91uaDCQeWbFj/sdJxsBvDCckMCCH0u3755dXcC+roBs5YTyqprGjeffyXt4QQeqVj\nYTeGE5KZWLdv7941e/Ys5W78rLdZt3fvz6t37lTVcFrs+nBCMiM/7d71p/2JJ7mHEes1jqakZK3f\nt+8PSsfBugYnJDOSkZdXsenggXdr6up4MjJm9mrq6gxr9+555xj3MjUbnJDMzNLNm7/6afeuDUrH\nwVh3W71r57pvN274r9JxsK7DCcnMCCHExgMHnjuacjpT6VgY6y7HzpzJ3nTw4O8EXzQ1K5yQzFBC\nclL2qh07/lJcXs7TpzKzU1lTo1u9e+ffEpKTspWOhXUtTkhmatnWLcu/37zpa4OBLycx87J829YV\n323c+JXScbCuxwnJjP37h6W/X79/X8K1l2RMG3YePZq0etfO55WOg3UPTkhmTAjRuGH//mdTs85f\nVDoWxm5UWnZ28aqd259KzcoqVToW1j04IZm5zQkHjy/buvWt2vp6brtjmlVeXd3y3eZNr27Yv59r\n/GaME1Iv8PX6dR8s2bB+GXdIYlqk1+vx9fp1Xy5Zv467eJs5Tki9gBBCvP3NksdW7dixR+lYGOus\nVTt37Hz/h6XPKR0H636ckHoJIUTDiu3bHtx74vhppWNhrKMOnUpO/3HnjseEEDy/Si/ACakXOXL6\nVO7K7dsfSzmfWah0LIxdS1ZhYeWqHTueP3zqVJbSsbCewQmpl1m7d8/h7zdt+n1+cXGt0rEwdjUX\nSksbvl639qWV27dtVjoW1nM4IfVC323auOKrdWvfrKmr4/ljmOpUVFe3fP7zT69/tW7t/5SOhfUs\nTki91P9+/umtz9f8/GWLTqd0KIxd0tDYKP7z0+r3P//5p7eVjoX1PCulA2DKEEIIInrKxsrK5anb\nF95laWmpdEisl2vR6fDZ6h+//HTVypc/wQqlw2EK4ITUiwkh9ET0gLWVle3jt952q4UFV5iZMgwG\nA75Yu2b1v39Y+iSP4N17cULq5YQQLUR0j42V9Y8PzZs3m4iUDon1MkIILN28edubX315nxCCr2v2\nYpyQGIQQjd7u7gvJgn5cNGfuLE5KrKcYDAZ8s3HDli/XrLlDCMHTpfRynJAYAKC4vLyBiG4HwEmJ\n9Qi9Xo8v1q75+f++/OIeTkYM4ITETAghGojodr1ev+KhufPmckcH1l2aW1rw359WL33n228eEkK0\nKB0PUwdOSOwKclJaUFtf/78nblvwgL2dHVeVWJdqaGwUn6xa+eUHy5c9ydeMmCnuVqUxRNSfiAQR\nLSCiDCKqIqJ1ROQpvz6WiA4QUR0RZRPR0519DyFEy3tLv3/o/eXL3i6trOQxxFiXqamrM3y4YvkH\nHyxf9jgnI9YaJyTt+iOAewHMBTAawItENATALgD7AEQDeA3AO0S0sLMbF0KIT1au+ONHK5b/Iaeo\nqLoL42a9VHFFRePHq1b+/eOVK37PXbtZW4j/LrSFiPoDyAIwTwixXi77F4DhABIB3CSEGG2y/FsA\n4oUQI6/3Pe+defPt98+e82FkSIjfjcTOeq/0nJzibzdu+OOSDeu/UjoWpl5cQ9KuTJPn1QCsAQwB\ncKTVcgkABt/IGy3dsvnHz35cdWdCclL6jWyH9U4JSUlnP1q54g5ORuxauFODdrW+tkMA2uo6a4ku\nOM9r9+7ZHxcdPTv3woXPF8ZPncw98Ni1GAwG/Lhz556f9+x6ZP/Jk+eVjoepHyck85IKIL5V2VgA\naV2x8f0nT2YQ0cz8ixffe2DO3Ee93dzsumK7zPxU1dbqv1q39pt/Lf3+KSFEk9LxMG3gJjvz8imA\nKCJ6k4gGEdEDAJ4G8HFXvYEQovn9ZT888++l3z92Mi0tt6u2y8xHem5O2bvff/fie99/9wgnI9YZ\n3KlBY0w6NQwUQmTIZYsBTBVCjCeiyQDeBRAJIBfAe0KI/3RHLJNjYyPmT5z82a2TJ8fxwKxMr9dj\n3f59h7YkHHx+44EDR5WOh2kPJyR2Q4jI7oX77n///lmzH/JwdbVROh6mjILi4rplW7d8+cHyZS9x\nrYhdL05IrEvcd/Ose2eNH/+XuOHRg3kcvN5DCIHdx46d3nBg/ysrt2/bqHQ8TNs4IbEuMyAgwP2u\nadPfWTAl/h4fDw97peNh3au8urpl+baty1Zu3/Z8Rl5ehdLxMO3jhMS63MKp0+bNHDfu9emjxwzn\n2pL5EUJg/8mT6VsOJfzj240bligdDzMfnJBYt/Ds4+bwwOzZby6Ij384yNfPWel4WNfIzM+vXL9/\n33fLt275a35xcZXS8TDzwgmJdat5EyZOjB816o3Z4+Mm2NlwnwetqqmrM6zbv2/LzqNHFm87fPgX\npeNh5okTEut2RGRx/6zZv5k2eszTk2JihnAznnYIIbDr2C8pWw8dem/Z1i1LeFBU1p04IbEe4+Hq\n6nj/rNl/iR816oHosMH+SsfD2ncq49zFbYcPf//txg2vlVVV1SkdDzN/nJBYjxs6cGDfOePj/jRl\n5Kg7B/fv76F0POxKyefSL+49cWLV1kOH/pmYzqNxsJ7DCYkpZvzw4YOnjhr98rhhw+aGBw/gxKSw\nxPT0CwcST67csH//P09nZuQrHQ/rfTghMcWNGDw4eNroMc/EDgm/bUxUVBBfY+pZielpRftOnFi5\n6eDBdzkRMSVxQmKq0cfZ2fnuGTOfGxkesXBSTMxQG2trpUMyWzq9HgeTEjMS09PXrN+3999ns7ML\nlY6JMU5ITHWIyOqeGTMfHjN06D0TR8SMd3dx4cmXukhZZaXuQFLigWNnzvy4ZMP6L4UQbc2hxZgi\nOCEx1SIimhwbO2XiiJjbBwb2mzYmKiqEa02dZzAYcPRMSt6Z8+e37D1+/Iedvxzdo3RMjLWFExLT\nBCKyWRg/9e7RkVGzhwwInhIVEurB15quTgiB1KysqrSc7ITE9PRtP+3e9UVFdXWt0nEx1h5OSExz\nAn19veaMj3skKiR0Wlj//qMH9evnyMlJSkLJGefKMvPyDp7KzNy348iRH7IKC4qUjouxjuKExDQt\nJCAgcOKImAXhAwaM8vP0HBU7JDzE0b73DDSu0+uRmJZ28XxBwYGkc+l7tyQkLL1YXlaudFyMXQ9O\nSMxsEJFV/MhRM0ZHRk4K8vMbFRoYOGJgYD8nc6o9Nbe0IDnjXGlRaempopKSxOSMcyf3HD++prKm\npkbp2Bi7UZyQmNnq5+vXd9zQoTcPCQ6O8HZzj/BwdQ2PCAnp6+rkpHRoHVbX0IBTmRkXLpaVJReU\nlCQlpqUd35xwcKMQgq8HMbPDCYn1GkRkEzNkyKjhAweNDPL3H+TZp89AJ3v7QA9XV7/+/n2dXRwd\nFYutrqEBaTk5FWVVVdk1dXU5ZdVVWfkXL2adOX/+zOHTpxKEEA2KBcdYD+GExHo9InKJDgsbOjCw\nX1iAj4+fh6urv7Ojo5+jnZ2/tZWVh7WVlaOzo6NDHydnR2cHB0snBwd0tPu5Xq9HVW0tyqqqmkoq\nK6oam5qrmltayhqbm0rrGhtLaurqSotKSwtSzmemHDl9+jDXfFhvxgmJsWsgIgsATgBcwoL6+7q7\nuvi6Obv06ePk5GBvZ+egNxh0er1e36LT6Vp0ukuPTS3NupYWXUtxRUVuWk52LoByIYRO4d1hTLU4\nITHGGFMFC6UDYIwxxgBOSIwxxlSCExJjjDFV4ITEGGNMFTghMcYYUwVOSIwxxlTh/wHWEgDB0fQ3\nDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111e228d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "blue_pal = sns.cubehelix_palette(9) # sns.cubehelix_palette(9, start=.5, rot=-.75) # sns.cubehelix_palette(9) #sns.color_palette(\"Reds\", 9)\n",
    "counts = df['significant'].value_counts()\n",
    "with sns.plotting_context('paper', font_scale=1.25):\n",
    "    pie_colors = [ blue_pal[5], blue_pal[1],]\n",
    "    labels = ['yes', 'no', ]\n",
    "    patches, texts, autotexts = plt.pie(counts[labels], colors=pie_colors, labels=labels, \n",
    "                                        autopct='%.1f%%', textprops={'fontsize':14})\n",
    "    [ _.set_fontsize(14) for _ in texts ]\n",
    "    plt.title('Genes with\\ncancer subtype enrichment', \n",
    "              fontdict={'verticalalignment': 'bottom', 'fontweight': 'bold', 'fontsize':16})\n",
    "    _ = plt.axis('equal')\n",
    "    plt.tight_layout()\n",
    "    plt.savefig('CHASMplus/figures/figure_3c.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def format_dataframe(df, ctype):\n",
    "    # only pick significant results for that cancer type\n",
    "    cancer_type_df = df[(df['cancer type']==ctype) & (df['q-value']<=0.1)].copy()\n",
    "\n",
    "    # split the information about subtypes\n",
    "    split_df = cancer_type_df['subtype information'].str.split(';', expand=True)\n",
    "\n",
    "    # extract out the percentages\n",
    "    percent_df = split_df.applymap(lambda x: x.split('%')[0]).astype(float)\n",
    "    subtypes = split_df.iloc[0,:].apply(lambda x: x.split(',')[-1][:-1]).values\n",
    "    n = split_df.iloc[0,:].apply(lambda x: x.split(',')[0].split('/')[-1]).values\n",
    "    subtypes = subtypes + ' (n=' + n + ')'\n",
    "    percent_df.columns = subtypes\n",
    "\n",
    "    # add the information about the percentage of samples mutated\n",
    "    cancer_type_df = pd.concat([cancer_type_df, percent_df], axis=1)\n",
    "\n",
    "    # set the index to be the gene names\n",
    "    cancer_type_df = cancer_type_df.set_index('gene', drop=False)\n",
    "    \n",
    "    return cancer_type_df[subtypes].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_heatmap(mydf, ctype):\n",
    "    cmap = mpl.colors.ListedColormap(['#ffffff']+sns.cubehelix_palette(100, start=.5, rot=-.75).as_hex()[0:100])\n",
    "\n",
    "    type2color = {'missense': \"#3498db\", 'inactivating': \"#e74c3c\"}\n",
    "    \n",
    "    with sns.plotting_context('paper'):\n",
    "        ax = sns.heatmap(mydf, \n",
    "                    annot=True, linecolor='gray', \n",
    "                    linewidth=.5, square=True,\n",
    "                    fmt='.3g',\n",
    "                    vmin=0, vmax=50.,\n",
    "                    xticklabels=1,\n",
    "                    #row_colors=cancer_type_df['mutation type'].apply(lambda x: type2color[x]),\n",
    "                    #col_cluster=False,\n",
    "                    #cbar_kws={'fraction': .02},\n",
    "                    #cbar=False,\n",
    "                    cmap=cmap)\n",
    "        #cg.ax_row_dendrogram.set_visible(False)\n",
    "        \n",
    "\n",
    "\n",
    "        # make sure tick labels are rotated in the appropriate manner\n",
    "        _ = plt.setp(ax.yaxis.get_majorticklabels(), rotation=0)\n",
    "        _ = plt.setp(ax.xaxis.get_majorticklabels(), rotation=90)\n",
    "\n",
    "        # remove x/y labels\n",
    "        ax.set_xlabel('')\n",
    "        ax.set_ylabel('')\n",
    "        \n",
    "        # remove ticks\n",
    "        ax.xaxis.set_ticks_position('none')\n",
    "        ax.yaxis.set_ticks_position('none')\n",
    "\n",
    "        # add cancer type as title\n",
    "        ax.set_title(ctype)\n",
    "\n",
    "        # fudge figure size to be compatible\n",
    "        vertical_fudge = 3\n",
    "        horizontal_fudge = 3\n",
    "        if plot_df.shape[1] <= 2: horizontal_fudge -= 4\n",
    "        if plot_df.shape[0] <= 2: vertical_fudge -= 4\n",
    "        \n",
    "        base_genes = 5.\n",
    "        base_subtypes = 5.\n",
    "        nsubtypes, ngenes = plot_df.shape\n",
    "        myfig = plt.gcf()\n",
    "        myfig.set_size_inches(ngenes/base_genes*4, nsubtypes/base_subtypes*5)\n",
    "    \n",
    "    return myfig, ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1a4ff550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#cancer_type = 'READ'\n",
    "for cancer_type in df['cancer type'].unique():\n",
    "    plot_df = format_dataframe(df, cancer_type)\n",
    "    fig, ax = create_heatmap(plot_df, cancer_type)\n",
    "    if plot_df.shape[1] > 2:\n",
    "        plt.tight_layout()\n",
    "    plt.savefig('CHASMplus/figures/subtype_analysis/{0}.pdf'.format(cancer_type))\n",
    "    plt.clf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>gene</th>\n",
       "      <th>PIK3CA</th>\n",
       "      <th>CTNNB1</th>\n",
       "      <th>PPP2R1A</th>\n",
       "      <th>PTEN</th>\n",
       "      <th>FGFR2</th>\n",
       "      <th>MAP3K1</th>\n",
       "      <th>TP53</th>\n",
       "      <th>NFE2L2</th>\n",
       "      <th>CTCF</th>\n",
       "      <th>KRAS</th>\n",
       "      <th>PIK3R2</th>\n",
       "      <th>CHD3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>CN_HIGH (n=146)</th>\n",
       "      <td>28.08</td>\n",
       "      <td>2.05</td>\n",
       "      <td>30.14</td>\n",
       "      <td>10.27</td>\n",
       "      <td>2.74</td>\n",
       "      <td>0.00</td>\n",
       "      <td>54.79</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>4.79</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CN_LOW (n=118)</th>\n",
       "      <td>48.31</td>\n",
       "      <td>43.22</td>\n",
       "      <td>1.69</td>\n",
       "      <td>32.20</td>\n",
       "      <td>11.02</td>\n",
       "      <td>0.85</td>\n",
       "      <td>4.24</td>\n",
       "      <td>5.08</td>\n",
       "      <td>1.69</td>\n",
       "      <td>16.95</td>\n",
       "      <td>0.85</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSI (n=109)</th>\n",
       "      <td>49.54</td>\n",
       "      <td>12.84</td>\n",
       "      <td>2.75</td>\n",
       "      <td>38.53</td>\n",
       "      <td>11.93</td>\n",
       "      <td>2.75</td>\n",
       "      <td>8.26</td>\n",
       "      <td>7.34</td>\n",
       "      <td>6.42</td>\n",
       "      <td>33.94</td>\n",
       "      <td>5.50</td>\n",
       "      <td>0.92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>POLE (n=6)</th>\n",
       "      <td>83.33</td>\n",
       "      <td>50.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>66.67</td>\n",
       "      <td>0.00</td>\n",
       "      <td>16.67</td>\n",
       "      <td>16.67</td>\n",
       "      <td>0.00</td>\n",
       "      <td>16.67</td>\n",
       "      <td>33.33</td>\n",
       "      <td>0.00</td>\n",
       "      <td>16.67</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "gene             PIK3CA  CTNNB1  PPP2R1A   PTEN  FGFR2  MAP3K1   TP53  NFE2L2  \\\n",
       "CN_HIGH (n=146)   28.08    2.05    30.14  10.27   2.74    0.00  54.79    0.00   \n",
       "CN_LOW (n=118)    48.31   43.22     1.69  32.20  11.02    0.85   4.24    5.08   \n",
       "MSI (n=109)       49.54   12.84     2.75  38.53  11.93    2.75   8.26    7.34   \n",
       "POLE (n=6)        83.33   50.00     0.00  66.67   0.00   16.67  16.67    0.00   \n",
       "\n",
       "gene              CTCF   KRAS  PIK3R2   CHD3  \n",
       "CN_HIGH (n=146)   0.00   4.79    0.00   0.00  \n",
       "CN_LOW (n=118)    1.69  16.95    0.85   0.00  \n",
       "MSI (n=109)       6.42  33.94    5.50   0.92  \n",
       "POLE (n=6)       16.67  33.33    0.00  16.67  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "format_dataframe(df, 'UCEC')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}