gene-popularity.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/pkerp/.virtualenvs/stuff/lib/python2.7/site-packages/IPython/core/magics/pylab.py:161: UserWarning: pylab import has clobbered these variables: ['f']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n",
      "  \"\\n`%matplotlib` prevents importing * from pylab and numpy\"\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "Cannot run multiple SparkContexts at once; existing SparkContext(app=pyspark-shell, master=local[*]) created by __init__ at <ipython-input-1-04367661b21a>:11 ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-36-04367661b21a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mpyspark\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0msc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpyspark\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSparkContext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/usr/local/opt/apache-spark/libexec/python/pyspark/context.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, master, appName, sparkHome, pyFiles, environment, batchSize, serializer, conf, gateway, jsc, profiler_cls)\u001b[0m\n\u001b[1;32m    113\u001b[0m         \"\"\"\n\u001b[1;32m    114\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_callsite\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfirst_spark_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mCallSite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 115\u001b[0;31m         \u001b[0mSparkContext\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_ensure_initialized\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgateway\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgateway\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconf\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mconf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    116\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    117\u001b[0m             self._do_init(master, appName, sparkHome, pyFiles, environment, batchSize, serializer,\n",
      "\u001b[0;32m/usr/local/opt/apache-spark/libexec/python/pyspark/context.pyc\u001b[0m in \u001b[0;36m_ensure_initialized\u001b[0;34m(cls, instance, gateway, conf)\u001b[0m\n\u001b[1;32m    273\u001b[0m                         \u001b[0;34m\" created by %s at %s:%s \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    274\u001b[0m                         % (currentAppName, currentMaster,\n\u001b[0;32m--> 275\u001b[0;31m                             callsite.function, callsite.file, callsite.linenum))\n\u001b[0m\u001b[1;32m    276\u001b[0m                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    277\u001b[0m                     \u001b[0mSparkContext\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_active_spark_context\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minstance\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: Cannot run multiple SparkContexts at once; existing SparkContext(app=pyspark-shell, master=local[*]) created by __init__ at <ipython-input-1-04367661b21a>:11 "
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%pylab inline\n",
    "\n",
    "import os.path as op\n",
    "import findspark\n",
    "import os\n",
    "findspark.init()\n",
    "\n",
    "import pyspark\n",
    "sc = pyspark.SparkContext()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'PrJ9mN8MomRmLASZNa65qe'"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import shortuuid\n",
    "shortuuid.uuid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "assembly = 'hg19'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data_dir = op.expanduser(\"~/data\")\n",
    "output_dir = op.join(data_dir, assembly)   # where all of the intermediate output will be stored\n",
    "base_ucsc_dir = op.join(data_dir, 'ucsc-data/{}'.format(assembly))  # where all of the files downloaded from UCSC will be stored\n",
    "\n",
    "import shutil\n",
    "\n",
    "# create a directory to store intermediate output files\n",
    "def get_outfile(table_name):\n",
    "    outfile = op.join(output_dir, 'genbank-output/{}'.format(table_name))\n",
    "    if op.exists(outfile):\n",
    "        shutil.rmtree(outfile)\n",
    "    return outfile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the chromosome lengths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 249250621 492449994\n"
     ]
    }
   ],
   "source": [
    "#wget http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/chromInfo.txt.gz\n",
    "#mv chromInfo.txt.gz ~/data/ucsc-data/hg19/\n",
    "\n",
    "def get_chrom_lengths(base_dir):\n",
    "    '''\n",
    "    Get the cumulative start positions for the chromosomes in an assembly. The chromosomes\n",
    "    will be sorted alphabetically by their names.\n",
    "    \n",
    "    :param base_dir: A directory containing meta data about a genome assembly\n",
    "    :return: A dictionary of the from { 'chr2': 234323432 }, showing at which position\n",
    "             chromosomes start.\n",
    "    '''\n",
    "    chromLengths = (sc.textFile(op.join(base_dir, 'chromInfo.txt.gz'))\n",
    "                    .map(lambda x: x.split('\\t'))\n",
    "                    .map(lambda x: {'chrom': x[0], 'length': int(x[1]) })\n",
    "                    .collect())\n",
    "    \n",
    "    cum_chrom_lengths = {}\n",
    "    curr_cum_lengths = 0\n",
    "    \n",
    "    for x in sorted(chromLengths, key=lambda x: -x['length']):\n",
    "        cum_chrom_lengths[x['chrom']] = curr_cum_lengths\n",
    "        curr_cum_lengths += x['length']\n",
    "        \n",
    "    return cum_chrom_lengths\n",
    "\n",
    "cum_chrom_lengths = get_chrom_lengths(base_ucsc_dir)\n",
    "\n",
    "print cum_chrom_lengths['chr1'], cum_chrom_lengths['chr2'], cum_chrom_lengths['chr3']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the refgene data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "base_dir=op.expanduser(\"~/data/genbank-data/hg19\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[((9606, 1), 2591067)]\n",
      "time taken 10.4018189907\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "human_gene2pubmed = (sc.textFile(op.join(base_dir, \"human/human_gene2pubmed\"))\n",
    "                        .filter(lambda x: x[0] != '#')\n",
    "                        .map(lambda x: x.split('\\t'))\n",
    "                        .map(lambda x: ((int(x[0]), int(x[1])), int(x[2])))\n",
    "                     )\n",
    "print human_gene2pubmed.take(1)\n",
    "\n",
    "t1 = time.time()\n",
    "pubmeds_set = set([x[1] for x in human_gene2pubmed.collect()])\n",
    "t2 = time.time()\n",
    "print \"time taken\", t2 - t1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pyspark import SparkContext, SparkConf\n",
    "\n",
    "# wget http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/refGene.txt.gz\n",
    "# mv refGene.txt.gz ~/data/ucsc-data/\n",
    "def parse_exon_positions(exon_positions_str):\n",
    "    return map(int, exon_positions_str.strip(\",\").split(','))\n",
    "\n",
    "def load_refgene_data(base_dir):\n",
    "    '''\n",
    "    Load the UCSC refgene data for a particular assembly.\n",
    "    \n",
    "    :param base_dir: The directory which contains the refGene.txt.gz file.\n",
    "    '''\n",
    "    refGene = (sc.textFile(op.join(base_dir, 'refGene.txt.gz'))\n",
    "               .map(lambda x: x.split('\\t'))\n",
    "               .map(lambda x: {'name': x[1],\n",
    "                               'chrom': x[2],\n",
    "                               'strand': x[3],\n",
    "                               'txStart': x[4],\n",
    "                               'txEnd': x[5],\n",
    "                               'cdsStart': x[6],\n",
    "                               'cdsEnd': x[7],\n",
    "                               'exonCount': x[8],\n",
    "                               'exonStarts': x[9].strip(','),\n",
    "                               'exonEnds': x[10].strip(','),\n",
    "                               'chromOffset': cum_chrom_lengths[x[2]],\n",
    "                               'genomeTxStart': cum_chrom_lengths[x[2]] + int(x[4]),\n",
    "                               'genomeTxEnd': cum_chrom_lengths[x[2]] + int(x[5]),\n",
    "                               'geneName': x[12],\n",
    "                               'geneLength': int(x[5]) - int(x[4]),\n",
    "                               })\n",
    "               .filter(lambda x: x['chrom'].find('_') == -1)\n",
    "               .filter(lambda x: x['name'] )\n",
    "            )\n",
    "    \n",
    "    return refGene\n",
    "\n",
    "refGene = load_refgene_data(base_ucsc_dir)\n",
    "### add the genomic position\n",
    "refgene_dict = dict([(t['name'], t) for t in refGene.collect()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "refgene length: 63539\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[u'NM_001003806',\n",
       " u'NM_001253387',\n",
       " u'NM_001003805',\n",
       " u'NM_001003802',\n",
       " u'NM_001003803',\n",
       " u'NM_001003800',\n",
       " u'NM_001003801',\n",
       " u'NM_001168531',\n",
       " u'NM_001003809',\n",
       " u'NM_001304464']"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print \"refgene length:\", len(refgene_dict)\n",
    "refgene_dict.keys()[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((9606, 1), (u'A1BG', u'alpha-1-B glycoprotein', u'protein-coding'))]"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "base_dir = op.join(data_dir, 'genbank-data/hg19')\n",
    "taxid_gene_info_human = (sc.textFile(op.join(base_dir, 'human/human_gene_info'))\n",
    "                   .filter(lambda x: x[0] != '#')\n",
    "                   .map(lambda x: x.split('\\t'))\n",
    "                   .filter(lambda x: x[0] == \"9606\")\n",
    "                   .map(lambda x: ((int(x[0]), int(x[1])),(x[2], x[8], x[9])))\n",
    "                   )\n",
    "taxid_gene_info_human.take(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((9606, 1), u'NM_130786')]"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "taxid_gene_refseq_id = (sc.textFile(op.join(base_dir, \"human/human_gene2refseq\"))\n",
    "                        .filter(lambda x: x[0] != '#')\n",
    "                        .map(lambda x: x.split('\\t'))\n",
    "                        .filter(lambda x: x[3].split('.')[0] in refgene_dict)\n",
    "                        .map(lambda x: ((int(x[0]), int(x[1])), (x[3].split('.')[0])))\n",
    "                        )\n",
    "taxid_gene_refseq_id.take(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time taken 8.58014893532\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "t1 = time.time()\n",
    "taxid_gene_info_refseq = taxid_gene_info_human.join(taxid_gene_refseq_id)\n",
    "\n",
    "t2 = time.time()\n",
    "print \"time taken\", t2 - t1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2069907  in pubmeds_set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(pubmeds_set)\n",
    "human_year_pmid = (sc.textFile(op.join(base_dir, 'recent_pmid_year.ssv'))\n",
    "                  .map(lambda x: x.split())\n",
    "                  .map(lambda x: (int(x[0]), int(x[1])))\n",
    "                  .filter(lambda x: x[1] in pubmeds_set))\n",
    "human_year_pmid_collected = human_year_pmid.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pmid_year = dict([(x[1], x[0]) for x in human_year_pmid_collected])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2007, 1975, 2007, 2010, 2006, 2002, 2007, 2007, 2011, 2015]\n"
     ]
    }
   ],
   "source": [
    "print([k for k in pmid_year.values()[:10]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((9606, 100847044, 16381832),\n",
       "  ((u'MIR5087', u'microRNA 5087', u'ncRNA'), 16381832))]"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "taxid_gene_info_pubmed = (taxid_gene_info_human.join(human_gene2pubmed)\n",
    "                    .filter(lambda x: x[1][1] in pmid_year)\n",
    "                                 .map(lambda x: ((x[0][0], x[0][1], x[1][1]), x[1])))\n",
    "taxid_gene_info_pubmed.take(1)        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((9606, 7157), ((u'TP53', u'tumor protein p53', u'protein-coding'), 8253)),\n",
       " ((9606, 7124), ((u'TNF', u'tumor necrosis factor', u'protein-coding'), 5233)),\n",
       " ((9606, 1956),\n",
       "  ((u'EGFR', u'epidermal growth factor receptor', u'protein-coding'), 4525))]"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get total citation counts over all time for each taxid, gene_id combo\n",
    "gene_counts = sorted(taxid_gene_info_pubmed.map(lambda x: ((x[0][0], x[0][1]), (x[1][0], 1)))\n",
    ".reduceByKey(lambda x1, x2: (x1[0], x1[1] + x2[1]))\n",
    ".collect(), key=lambda x: -x[1][1])\n",
    "\n",
    "gene_counts[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[((9606, 2932, 2015), ((u'GSK3B', u'glycogen synthase kinase 3 beta', u'protein-coding'), 74))]\n",
      "301224\n",
      "[((9606, 109864272, 2015), ((u'RNA28SN4', u'RNA, 28S ribosomal N4', u'rRNA'), 1))]\n"
     ]
    }
   ],
   "source": [
    "taxid_gene_info_year = (taxid_gene_info_pubmed.map(lambda x: ((x[0][0], x[0][1], pmid_year[x[1][1]]), x[1]))\n",
    "                        .map(lambda x: ((x[0][0], x[0][1], pmid_year[x[1][1]]), (x[1][0], 1)))\n",
    "                                 .reduceByKey(lambda x1, x2: (x1[0], x1[1] + x2[1])))\n",
    "print taxid_gene_info_year.take(1)\n",
    "print taxid_gene_info_year.count()\n",
    "print taxid_gene_info_year.filter(lambda x: x[1][0][2] == 'rRNA').take(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "genes_per_year = dict(taxid_gene_info_pubmed.map(lambda x: ((pmid_year[x[1][1]], x[0][1]), 1))\n",
    " .reduceByKey(lambda x1,x2: x1+x2)\n",
    " .map(lambda x: (x[0][0], 1))\n",
    " .reduceByKey(lambda x1,x2: x1+x2)\n",
    " .collect())\n",
    "\n",
    "citations_per_year = dict((taxid_gene_info_pubmed.map(lambda x: (pmid_year[x[1][1]], 1))\n",
    "                      .reduceByKey(lambda x1,x2: x1+x2)\n",
    "                      .collect()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "citation_genes = sorted(taxid_gene_info_pubmed.map(lambda x: (x[1][1], 1))\n",
    " .reduceByKey(lambda x1, x2: x1 + x2)\n",
    " .collect(), key=lambda x: -x[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(12477932, 19866),\n",
       " (15489334, 11521),\n",
       " (28514442, 10943),\n",
       " (14702039, 10306),\n",
       " (26186194, 7656),\n",
       " (16344560, 6597),\n",
       " (26496610, 5279),\n",
       " (25416956, 4387),\n",
       " (26344197, 3473),\n",
       " (18029348, 3323)]"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(citation_genes[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(u'pseudo', 10997),\n",
       " (u'biological-region', 10171),\n",
       " (u'tRNA', 818),\n",
       " (u'snoRNA', 1513),\n",
       " (u'snRNA', 211),\n",
       " (u'unknown', 1893),\n",
       " (u'protein-coding', 1142021),\n",
       " (u'scRNA', 16),\n",
       " (u'rRNA', 1483),\n",
       " (u'ncRNA', 31333),\n",
       " (u'other', 2589)]"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_types_citations = (taxid_gene_info_pubmed.map(lambda x: (x[1][0][2], 1))\n",
    " .reduceByKey(lambda x1,x2: x1+x2)\n",
    " .collect()\n",
    ")\n",
    "print(gene_types_citations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pseudo 5841\n",
      "unknown 1172\n",
      "biological-region 1597\n",
      "snoRNA 540\n",
      "snRNA 64\n",
      "tRNA 510\n",
      "protein-coding 19426\n",
      "scRNA 4\n",
      "rRNA 43\n",
      "ncRNA 5656\n",
      "other 739\n"
     ]
    }
   ],
   "source": [
    "gene_types_counts = (taxid_gene_info_pubmed.map(lambda x: ((x[1][0][2], x[0][1]),1))\n",
    " .reduceByKey(lambda x1,x2: x1+x2)\n",
    " .map(lambda x: (x[0][0], 1))\n",
    " .reduceByKey(lambda x1,x2: x1+x2)\n",
    " .collect())\n",
    "\n",
    "for gtc in gene_types_counts:\n",
    "    print gtc[0], gtc[1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "string_values = (taxid_gene_info_year.map(lambda x: \"\\t\".join(map(str, \n",
    "                 [x[0][0], x[0][1], x[0][2], x[1][0][0], x[1][0][1], x[1][0][2], x[1][1]\n",
    "                 ])))\n",
    "                 .collect())\n",
    "with open('gene_info_by_year.tsv', 'w') as f:\n",
    "    for sv in string_values:\n",
    "        f.write(sv + \"\\n\")\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAF5CAYAAAB0sJvmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X+cXXV95/HXhx+TEGwSMRKgGkWRELFaE36u8stYqUIt\nDWoNzSJQ66rIslFXrRWhsK1KV0Llh0uFLZZAdlnQlQoSAX9UEKEytEUJcbXgqEwCIyFEwiQk+ewf\n50y4ucxkZu7cO3funNfz8biPzD3f7z33e76Z3PvO93y/50RmIkmSNNnt0u4GSJIkjQdDjyRJqgRD\njyRJqgRDjyRJqgRDjyRJqgRDjyRJqgRDjyRJqgRDjyRJqgRDjyRJqgRDjyRJqoS2h56I+POIuDci\nnoqItRHx1Yg4sK7O30fEtrrHLXV1pkTEZRHRFxEbIuKGiNi7rs4LI+LaiFgfEesi4sqI2LOuzksj\n4uaIeDoi1kTEhRHR9n6SJEljMxG+zI8CLgEOB94M7A58MyL2qKv3DWA2sE/5WFxXfjFwAnAycDSw\nH3BjXZ3rgHnAwrLu0cAVA4VluLkF2A04AngPcBpw/hiOT5IkTQAx0W44GhGzgMeAozPzznLb3wMz\nMnPREK+ZDjwOvDszv1pumwusAo7IzHsjYh7wY2BBZt5f1jkeuBl4SWauiYi3AjcB+2ZmX1nnPwGf\nBV6cmVtaduCSJKmlJsJIT72ZQAJP1G0/tjz99VBEXB4Re9WULaAYnbljYENmrgZ6gCPLTUcA6wYC\nT+n28r0Or6nzwEDgKa0EZgAHj+2wJElSO02o0BMRQXGa6s7MfLCm6BvAqcCbgI8BxwC3lPWhON21\nOTOfqtvl2rJsoM5jtYWZuZUiXNXWWTvIPqipU9/maRExPyKmDX+EkiRpwHh/h+42Hm8yCpcDrwbe\nULsxM6+vefrjiHgA+BlwLPDtcWvd4H4XuAvojojf1JXdSjFSJElS1R0P/H7dthcA8ym+97/f6gZM\nmNATEZcCbwOOyszendXNzIcjog84gCL0rAG6ImJ63WjP7LKM8s/61Vy7AnvV1Tm07u1m15QN5uXl\nn/MHKTsa+OudHYskSeLlVCX0lIHnD4FjMrNnBPVfArwIGAhH9wFbKFZl1U5kngPcXda5G5gZEa+v\nmdezEAjgnpo6n4yIWTXzet4CrAdqT7fVegRg+fLlzJs3b/iDFQBLly5l2bJl7W5Gx7HfRs8+a4z9\nNnr22eitWrWKJUuWQPld2mptDz0RcTnF8vO3A09HxMDIyvrM7C+vo3MuxfLzNRSjO58DfkJ56igz\nn4qIq4CLImIdsAH4AnBXZt5b1nkoIlYCX4qIDwBdFEvlV2TmwCjONynCzTUR8XFgX+AC4NLMfHaI\nQ+gHmDdvHvPnDzbYo8HMmDHD/mqA/TZ69llj7LfRs8/GpH883qTtoQd4P8UKqu/UbT8d+AdgK/Ba\nionMM4FHKcLOp+uCyNKy7g3AFIr5NGfW7fMU4FKKVVvbyrpnDxRm5raIOBH4IsUw29PA1RShS5Ik\ndbC2h57M3OkKsszs5/kTnwartwk4q3wMVedJYMkw+/kFcOJw7ydJkjrLhFqyLkmS1CqGHrXF4sX1\ndxHRSNhvo2efNcZ+Gz37bOKbcLeh6DQRMR+477777nMCmyRJo9Dd3c2CBQuguEVUd6vfz5EeSZJU\nCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYe\nSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJU\nCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYe\nSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCYYeSZJUCbu1uwGSJFVBT08P\nfX19AMyaNYs5c+a0uUXVY+iRJKnFenp6mDt3Hv39GwGYOnUaq1evMviMM09vSZLUYn19fWXgWQ4s\np79/4/ZRH40fR3okSRo389rdgEpzpEeSJFWCoUeSJFWCoUeSJFWCoUeSJFWCoUeSJFWCoUeSJFWC\noUeSJFVC20NPRPx5RNwbEU9FxNqI+GpEHDhIvfMj4tGI2BgRt0XEAXXlUyLisojoi4gNEXFDROxd\nV+eFEXFtRKyPiHURcWVE7FlX56URcXNEPB0RayLiwohoez9JkqSxmQhf5kcBlwCHA28Gdge+GRF7\nDFSIiI8DHwLeBxwGPA2sjIiumv1cDJwAnAwcDewH3Fj3XtdRXBlqYVn3aOCKmvfZBbiF4qKNRwDv\nAU4Dzm/KkUqSpLZp+xWZM/Nttc8j4jTgMWABcGe5+Wzggsz8elnnVGAtcBJwfURMB84A3p2Z3y3r\nnA6siojDMvPeiJgHHA8syMz7yzpnATdHxEczc01ZfhBwXGb2AQ9ExDnAZyPivMzc0rqekCRJrTQR\nRnrqzQQSeAIgIvYH9gHuGKiQmU8B9wBHlpsOoQhwtXVWAz01dY4A1g0EntLt5XsdXlPngTLwDFgJ\nzAAObsKxSZKkNplQoSciguI01Z2Z+WC5eR+KYLK2rvrasgxgNrC5DEND1dmHYgRpu8zcShGuausM\n9j7U1JEkSR2o7ae36lwOvBp4Q7sbMlpLly5lxowZO2xbvHgxixcvblOLJEmaOFasWMGKFSt22LZ+\n/fpxbcOECT0RcSnwNuCozOytKVoDBMVoTu0ozGzg/po6XRExvW60Z3ZZNlCnfjXXrsBedXUOrWva\n7JqyIS1btoz58+fvrIokSZU12EBAd3c3CxYsGLc2TIjTW2Xg+UOKCcQ9tWWZ+TBF4FhYU386xTyc\n75eb7gO21NWZC8wB7i433Q3MjIjX1+x+IUWguqemzu9ExKyaOm8B1gMPIkmSOlbbR3oi4nJgMfB2\n4OmIGBhZWZ+Z/eXPFwOfioifAo8AFwC/BL4GxcTmiLgKuCgi1gEbgC8Ad2XmvWWdhyJiJfCliPgA\n0EWxVH5FuXIL4JsU4eaacpn8vuV7XZqZz7asEyRJUsu1PfQA76eYqPyduu2nA/8AkJkXRsQ0imvq\nzAS+B7w1MzfX1F8KbAVuAKYAtwJn1u3zFOBSilVb28q6Zw8UZua2iDgR+CLFKNLTwNXAuWM8RkmS\n1GZtDz2ZOaJTbJl5HnDeTso3AWeVj6HqPAksGeZ9fgGcOJI2SZLULD09PfT1FVdMmTVrFnPmzGlz\niyaftoceSZKqrqenh7lz59HfvxGAqVOnsXr1KoNPk02IicySJFVZX19fGXiWA8vp79+4fdRHzeNI\njyRJE8a8djdgUnOkR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5Ik\nVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKh\nR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVYKhR5IkVcJuI6kUEf95pDvM\nzC803hxJkqTWGFHoAZbWPX8xMA14snw+E9gIPAYYeiRJ0oQzotNbmbn/wAP4C+BfgHmZuVdm7gXM\nA7qBc1rXVEmSpMY1MqfnAuCszFw9sKH8eSnw35rVMEmSpGYa6emtWvsO8bpdgdlja44kSQLo7e2l\nu7sbgFmzZjFnzpw2t6jzNRJ67gCuiIj3ZmY3QEQsAL4I3N7MxkmSVFWLFr2TzZufAWDq1GmsXr3K\n4DNGjZzeOgNYA/wwIjZFxCbgXmAt8N5mNk6SpKoqAs9yYDn9/Rvp6+trd5M63qhHejLzceBtEXEg\ncFC5+aHM/ElTWyZJUuXNa3cDJpVGTm8NeAQI4GeZuaU5zZEkSWqNUZ/eiohpEXEVxXV5fgzMKbdf\nEhGfaHL7JEmSmqKROT2fAV4HHAv012y/HfjjJrRJkiSp6Ro5vXUS8MeZ+YOIyJrtPwZe2ZxmSZIk\nNVcjIz0vprjdRL09gRxkuyRJUts1Enp+CJxQ83wg6LwXuHvMLZIkSWqBRk5vfRL4RkS8unz92eXP\n/wE4ppmNkyRJapZRj/Rk5p3A71IEngeAt1Cc7joyM+9rbvMkSZKao6Hr9GTmz4A/a3JbJEmSWqaR\n6/S8LSKOH2T78RHx1uY0S5Ikqbkamcj82SG2x07KJEmS2qqR0PMqYPUg2x8CDhhbcyRJklqjkdCz\nHnjFINsPAJ4eW3MkSZJao5HQ8zXg4ojYfvXliDgA+DxwU7MaJkmS1EyNhJ6PUYzoPBQRD0fEw8Aq\n4NfAR5vZOEmSpGZp5Do96ykuRHgCcDnFCM/CzHxTZj7ZSCMi4qiIuCkifhUR2yLi7XXlf19ur33c\nUldnSkRcFhF9EbEhIm6IiL3r6rwwIq6NiPURsS4iroyIPevqvDQibo6IpyNiTURcGBGNhENJkiaE\nnp4euru76e7upqenp93NaZtGr9OTwDfLRzPsCfwLcBXwlSHqfAM4jWKVGMCmuvKLgbcCJwNPAZcB\nNwJH1dS5DpgNLAS6gKuBK4AlAGW4uQV4FDgC2A+4BtgMfKqxQ5MkqX16enqYO3ce/f0bAZg6dRqr\nV69izpw5bW7Z+Gso9ETEQorgsDd1o0WZecZo95eZtwK3lvuOIaptyszHh2jPdOAM4N2Z+d1y2+nA\nqog4LDPvjYh5wPHAgsy8v6xzFnBzRHw0M9eU5QcBx2VmH/BARJwDfDYizsvMLaM9NkmSxltPTw99\nfX0A9Pb2loFnOQD9/Uvo6+sbU+ip3f+sWbM6JkCNOvRExLnApyluPNrL+N1Z/diIWAusA74FfCoz\nnyjLFlAcyx0DlTNzdUT0AEcC91KM3KwbCDyl2ynafzjFBO0jgAfKwDNgJfBF4GDgX1txYJIkNUv9\nyE5X19SyZF5L9t9JI0eNjPS8HzgtM69pdmN24hsUp6oeBl4JfAa4JSKOLE+17QNszsyn6l63tiyj\n/POx2sLM3BoRT9TVWTvIPgbKDD2SpAmtr69vh5GdzZuXtHT/zRg5Gi+NhJ4u4PvNbsjOZOb1NU9/\nHBEPAD8DjgW+PZ5tGcrSpUuZMWPGDtsWL17M4sWL29QiSVK1NWdkp1n7X7FiBStWrNhh2/r165vZ\noGE1EnquBE4BLmhyW0YsMx+OiD6KCyJ+G1gDdEXE9LrRntllGeWf9au5dgX2qqtzaN3bza4pG9Ky\nZcuYP3/+aA9FkqRKGGwgoLu7mwULFoxbGxoJPVOB90XEm4F/A56tLczMDzejYTsTES8BXkQxpwjg\nPmALxeTqr5Z15gJzgLvLOncDMyPi9TXzehZSrAa7p6bOJyNiVs28nrdQXIX6wdYdkSRJarVGQs9r\nKZaXA7ymrqyhSc3ltXIO4Lnl6K+IiNcBT5SPcynm9Kwp630O+AnFJGMy86mIuAq4KCLWARuALwB3\nZea9ZZ2HImIl8KWI+ADFabpLgBXlyi0oluA/CFwTER8H9qUY0bo0M3cId5IkVUX9arBONerQk5nH\ntaAdh1Ccpsry8fly+5eBD1IErVOBmRTX0FkJfLouiCwFtgI3AFMolsCfWfc+pwCXUqza2lbWPXug\nMDO3RcSJFKu1vk9x5emrKUKXJEmVM/RqsM7T0HV6YPv9tl4J/FNmPhMRUa6kGrXy2jo7u+rx749g\nH5uAs8rHUHWepLwQ4U7q/AI4cbj3kySpClq9Gmw8NXKdnhcB1wPHUYzKvAr4d+CqiFiXmR9pbhMl\nSVL7tXo1WOs1ck+pZRSTl+cAG2u2/29GMCIjSZLUDo2c3noLcHxm/rLujhH/D3hZU1olSZLUZI2M\n9OzJjiM8A/bi+TcBlSRJmhAaCT3fo1hJNSDLu5N/jAlydWRJkqR6jZze+hhwR0QcQnGtmwspbsa5\nF/CGJrZNkiSpaRq5Ts+PIuJA4EMUFwF8AfAV4LLM7NwrFkmSVFG1Fx+cNWtWR9w8tBGjCj0RsRvw\nSeB/ZuZftaZJkiRpvNRffHDq1GmsXr1qUgafUc3pycwtFKe3Gr6ooSRJmjh2vPjgcvr7N24f9Zls\nGgkvdwDHAI80tymSJGk89Pb20t3dvf3nwrxhyjtfI6HnG8BnI+J3KO5u/nRtYWbe1IyGSZKk1li0\n6J1s3vwMMPi9tIYr71SNhJ7Lyz8/PEhZArs23hxJktRqRaAZ+l5aw5V3qkZWbzVybR9JkjShDHcv\nrc6/11a9MQWYiJg8Y16SJGlSG3XoiYhdI+KciPgV8JuIeEW5/YKI+NOmt1CSJKkJGhnp+QvgNIql\n65trtv8IeG8T2iRJktR0jYSeU4H3Zea1wNaa7f8KHNSUVkmSJDVZI6Hnt4GfDrGv3cfWHEmSpNZo\nJPQ8CBw1yPZ3APePrTmSJEmt0ch1es4HvhwRv00RmhZFxFyK014nNrNxkiRJzTLqkZ7M/BrwB8Cb\nKa7GfD7FYv4/yMzbmts8SZKk5hjRSE9E/Gfg7zKzPyLmAHdm5u+1tmmSJEnNM9LTWxcB/wvoBx4G\n9gUea1WjJElSZ+rp6dl+l/ZZs2YxZ86cNrfoOSMNPY8CJ0fELUAALxnqasyZ2dOsxkmSpM7R09PD\n3Lnz6O/fCMDUqdNYvXrVhAk+Iw09/w24BLiU4qai/zxIncAbjkqSVFl9fX1l4CluVtrfv4S+vr7O\nCj2Z+XcRsQJ4GfBvFJOYf93KhkmSpKFN5NNIE/VmpSNesp6ZG4AfRcTpwF2Zual1zZIkSUOZ6KeR\nJqpGlqx/GdgjIt4bEZ+JiL0AImJ+ee0eSZLUQjueRlpOf//G7aM+GtqoL04YEa8FbgfWAy8HvgQ8\nASwC5lBcpFCSJLXcxDyNNFE1chuKZcDVmfkqiiXsA24Bjm5KqyRJkpqskdtQHAK8b5DtvwL2GVtz\nJEmSWqORkZ5NwPRBth8IPD625kiSJLVGIyM9NwGfjoh3lc+zvDXF54Abm9YySZK0XW9vL93d3dt/\n1ug1Eno+AtxAcRuKPYDvUpzWuhv4i+Y1TZIkDVi06J1s3vwMAF1dg94UQcMYdejJzPXA70XEG4DX\nAS8AujPz9mY3TpIkFYrAs7z8eUl7G9OhRhV6ImJ34Fbg/Zl5F3BXS1olSZIG4RL1sRhV6MnMZ8vr\n9EiSpAlk1apVgPN9dqaROT3LgT8FPtHktkiSpFHrBXZhyZLilJfzfYbWSOjZDTgjIt4M3Ac8XVuY\nmR9uRsMkSdJIPAlsw/k+w2sk9LwG6C5/PrCuLMfWHEmS1Bjn+wynkdVbx7WiIZIkVYnX3Rl/jYz0\nSJKkMfK6O+PP0CNJUht43Z3xZ+iRJKltnIcznhq54agkSVLHMfRIkqRKGHXoiYj3RMQJNc8vjIgn\nI+L7EfGy5jZPkiSpORoZ6fkk8AxARBwJnAl8DOgDljWvaZIkSc3TyETmlwI/LX8+CbgxM/8uIu4C\nvtOshkmSJDVTIyM9vwFeVP78FuC28ud+YI9mNEqSJKnZGhnpuQ24MiLup7gNxS3l9oOBR5rULkmS\npKZqZKTnTOBu4MXAyZn563L7AmBFI42IiKMi4qaI+FVEbIuItw9S5/yIeDQiNkbEbRFxQF35lIi4\nLCL6ImJDRNwQEXvX1XlhRFwbEesjYl1EXBkRe9bVeWlE3BwRT0fEmnKitqvcJEnqcKP+Ms/MJzPz\nQ5n5h5l5a832czPzrxpsx57AvwAfZJCblkbEx4EPAe8DDqO4s/vKiOiqqXYxcAJwMnA0sB9wY92u\nrqO4EtTCsu7RwBU177MLxcjVbsARwHuA04DzGzwuSZI0QTQ0glGOzCwvl6n/drntP0bEGxvZX2be\nmpmfzsyvATFIlbOBCzLz65n5I+BUilBzUvne04EzgKWZ+d3MvB84HXhDRBxW1pkHHA/8aWb+MDO/\nD5wFvDsi9inf53jgIOBPMvOBzFwJnAOcGRFevVqSpA7WyHV6TgZWUixbnw9MKYtmUCxnb6qI2B/Y\nB7hjYFtmPgXcAxxZbjqEYnSmts5qoKemzhHAujIQDbidYmTp8Jo6D2RmX02dlRTHdnCTDkmSJLVB\nIyM9nwLen5l/Bjxbs/0uihDUbPtQBJO1ddvXlmUAs4HNZRgaqs4+wGO1hZm5FXiirs5g70NNHUmS\n1IEaCT1zgX8aZPt6YObYmiNJktQajcxTWQMcwPOXp78R+PexNmiI9wuK0ZzaUZjZwP01dboiYnrd\naM/ssmygTv1qrl2BverqHFr3/rNryoa0dOlSZsyYscO2xYsXs3jx4p29TJKkjrdq1SoAent7h6yz\nYsUKVqzYcZH3+vXrW9queo2Eni8BfxsRZ1CcdtqvvB3FfwcuaGbjADLz4YhYQ7Hi6t9g+8Tlw4HL\nymr3AVvKOl8t68wF5lAsr6f8c2ZEvL5mXs9CikB1T02dT0bErJp5PW+hGMV6cGftXLZsGfPnt+Ls\nniSpE/X09NDXV3yV7CwMdLZeYBeWLFkCQFfX1CFrDjYQ0N3dzYIFC1rZwB00Eno+S3Fa7A5gGsWp\nrk3Af8/MSxppRHmtnAN4buXWKyLidcATmfkLiuXon4qIn1KMMF0A/BL4GhQTmyPiKuCiiFgHbAC+\nANyVmfeWdR6KiJXAlyLiA0AXcAmwIjMHRnG+SRFurimXye9bvtelmVk7f0mSpCH19PQwd+48+vs3\nAjsPA53tSWAbsByAzZuXtLU1wxl16MnMBP4qIv6GIqi8AHgwM38zhnYcAnybYuQogc+X278MnJGZ\nF0bENIpr6swEvge8NTM31+xjKbAVuIFiRdmtFBdSrHUKcCnFqq1tZd2za45tW0ScCHwR+D7F9YCu\nBs4dw7FJkiqmr6+vDDydEQbGbl67GzAiDV97pgwcOz3lM4p9fZdhJlVn5nnAeTsp30Rx3Z2zdlLn\nSWCnv3nlyNKJO6sjSdLIdEYYqIpRh57yVNQnKObD7E1dWMnMVzSnaZIkSc3TyEjPlcAxwDUUM5ie\nd9sISZKkiaaR0PNW4ITMvKvZjZEkSWqVRi5OuI7iKsaSJEkdo5HQcw5wfrmaSpIkqSM0cnrrI8Ar\ngbUR8Qg73n+LzPQKfZIkacJpJPT836a3QpIkqcUauTjhX7aiIZIkdbJq3HaiszV0ccKImAm8g+I0\n199k5hMRMR9Ym5m/amYDJUma6Kpz24nO1sjFCV9LcRuH9cDLKW5A+gSwiOIGn6c2sX2SJE141bvt\nRGdqZPXWRcDVmfkqoL9m+y3A0U1plSRJHWke3npi4mok9BxKcePPer8C9hlbcyRJklqjkdCzCZg+\nyPYDgcfH1hxJkqTWaCT03AR8OiJ2L59nRMwBPgfc2LSWSZIkNVEjoecjwAuAx4A9gO8CPwU2AH/R\nvKZJkiQ1TyPX6VkP/F5EvBF4LUUA6s7M25vdOEmSpGZp6Do9AJl5J3BnE9siSZLUMiMOPRGxB7Aw\nM79ePv8MMKWmylbgnMzsH+z1kiRJ7TSakZ73ACcAXy+ffwj4MfBM+fwg4FFgWdNaJ0mS1CSjmcj8\nJ8Df1W07JTOPy8zjgP8KvKtpLZMkSWqi0YSeA4AHap73A9tqnt8LvLoZjZIkSWq20ZzemknNHJ7M\nfHFd+S7sOMdHkiRpwhjNSM8vgdfspPy1ZR1JkqQJZzSh5xbg/IiYWl9Qruw6F7i5WQ2TJElqptGc\n3vprionKqyPiUuAn5fa5FCu5divrSJIkTTgjDj2ZuTYi/gPwReCzQAwUAbcBH8zMtc1voiRJ0tiN\n6orMmfkw8PsRsRfFai6An2bmE01vmSRJUhM1dBuKMuTc2+S2SJIktUzD996SJGmy6enpoa+vD4BZ\ns2YxZ86cEdVftWrVeDRPY2TokSSJIsDMnTuP/v6NAEydOo3Vq1cNGXzq62viG82SdUmSJq2+vr4y\nwCwHltPfv3H7qM/w9S8Yp1ZqLBzpkSRpB/OGLKk9/dXb2ztsfU0shh5Jkkag/nRWV9fzrtWrCc7Q\nI0nSCOx4Ogs2b17S3gZp1Aw9kiSNiqezOpUTmSVJUiUYeiRJUiUYeiRJUiUYeiRJUiUYeiRJUiUY\neiRJUiUYeiRJUiUYeiRJUiUYeiRJUiUYeiRJUiUYeiRJUiUYeiRJUiUYeiRJUiUYeiRJUiUYeiRJ\nUiUYeiRJUiV0ROiJiHMjYlvd48G6OudHxKMRsTEibouIA+rKp0TEZRHRFxEbIuKGiNi7rs4LI+La\niFgfEesi4sqI2HM8jlGSJLVWR4Se0o+A2cA+5eONAwUR8XHgQ8D7gMOAp4GVEdFV8/qLgROAk4Gj\ngf2AG+ve4zpgHrCwrHs0cEULjkWS1AF6e3vp7u6mu7ub3t7edjdHY7RbuxswClsy8/Ehys4GLsjM\nrwNExKnAWuAk4PqImA6cAbw7M79b1jkdWBURh2XmvRExDzgeWJCZ95d1zgJujoiPZuaalh6dJGnC\nWbTonWze/AwAXV1T29wajVUnjfS8KiJ+FRE/i4jlEfFSgIjYn2Lk546Bipn5FHAPcGS56RCKgFdb\nZzXQU1PnCGDdQOAp3Q4kcHhrDkmSNJ56enq2j9z09PQMW78IPMuB5Wze3N/y9qm1OmWk5wfAacBq\nYF/gPOCfIuI1FIEnKUZ2aq0ty6A4Lba5DEND1dkHeKy2MDO3RsQTNXUkSR2qp6eHuXPn0d+/EYCp\nU6exevUq5syZM8wr57W+cRoXHRF6MnNlzdMfRcS9wM+BdwEPtadVkqRO0tfXVwae5QD09y+hr69v\nBKFHk0VHhJ56mbk+In4CHAB8BwiK0Zza0Z7ZwMCpqjVAV0RMrxvtmV2WDdSpX821K7BXTZ0hLV26\nlBkzZuywbfHixSxevHiERyVJGh+O3LTDihUrWLFixQ7b1q9fP65t6MjQExEvoAg8X87MhyNiDcWK\nq38ry6dTzMO5rHzJfcCWss5XyzpzgTnA3WWdu4GZEfH6mnk9CykC1T3DtWnZsmXMnz+/CUcnSdLk\nM9hAQHd3NwsWLBi3NnRE6ImIvwH+keKU1m8Dfwk8C/yvssrFwKci4qfAI8AFwC+Br0ExsTkirgIu\nioh1wAb77fH4AAAOLUlEQVTgC8BdmXlvWeehiFgJfCkiPgB0AZcAK1y5JUlS5+uI0AO8hOIaOi8C\nHgfuBI7IzF8DZOaFETGN4po6M4HvAW/NzM01+1gKbAVuAKYAtwJn1r3PKcClFKu2tpV1z27RMUmS\npHHUEaEnM4edGJOZ51Gs6hqqfBNwVvkYqs6TwJLRt1CSJE10nXSdHkmSpIYZeiRJUiUYeiRJUiUY\neiRJUiUYeiRJUiV0xOotSZJaobe3l+7u7u0/a3Iz9EiSKmvRoneWd1KHrq6pbW6NWs3QI0mqrCLw\nLC9/9jJtk52hR5JUcd6AtCqcyCxJkirB0CNJkirB0CNJkirB0CNJkirB0CNJkirB0CNJkirB0CNJ\nkirB0CNJkirB0CNJkirB0CNJkirB21BIkiatnp4e+vr6AO+iLkOPJGmS6unpYe7cefT3bwS8i7oM\nPZKkSaqvr68MPN5FXQVDjyRpkvMu6io4kVmSJFWCIz2SpEnDicvaGUOPJGlScOKyhmPokSRNCk5c\n1nAMPZKkScaJyxqcE5klSVIlGHokSVIlGHokSVIlOKdHktTRBpapr1q1qt1N0QRn6JEkdaz6ZerS\nzhh6JEkdpf4ChM8tU38YOKedTdMEZ+iRJHWMoS9A6DJ1Dc/QI0nqGF6AUGNh6JEkdSBHdjR6LlmX\nJEmV4EiPJGlC887pahZDjyRpwvLO6WomQ48kaUIZekm6E5c1NoYeSdKE4ZJ0tZKhR5I0YbgkXa1k\n6JEktdXgE5Ud2VHzGXokSW3jRGWNJ0OPJGlcOVFZ7WLokSSNGycqq50MPZKkceNEZbWToUeS1FJO\nVNZEYeiRJLWME5U1kRh6JElN5URlTVTeZX0QEXFmRDwcEc9ExA8i4tB2t2myWbFiRbub0JHst9Gz\nzxozXL/19PTQ3d1Nd3c3PT0925/ffPPNHHjgQSxYsIAFCxawaNE7ylfMw1NaajdHeupExB8Dnwfe\nB9wLLAVWRsSBmdnX1sZNIitWrGDx4sXtbkbHsd9Gzz5rTH2/1Y7ebNq0iTe96c3bT1lNmTIVCDZt\neqZmD47saOIx9DzfUuCKzPwHgIh4P3ACcAZwYTsbJkk7Ux9MpkyZAsCsWbOYM2fOqPb1zDPP0N3d\nvX1ftSGnq2sqmzf3MxBsNm0aCDbLgYeBc3BURxORoadGROwOLAD+emBbZmZE3A4c2baGSaqk2hAz\nXHCpnzAMuwJbAZg6dRrf+tbtO4QgYId91z7ftGkT3/rWt1mwYAHw/JDz3OhNfbAx6GhiM/TsaBbF\nJ8Xauu1rgbmN7nQ0H1xSu4z193Rnrx9u3/XlMPQX8mifr1u3bvuIRbP3PZLwMNKgMdhra0dXhgsu\nO04YHhhtKUJKf/8Sjj12IZs3F6ef6k9H1T/v6prKtm1bGT7kSJ3F0DN2UwFWrVo1aGFvby+LFr2j\n/F9S8WHyla/cwL777jt+LZyA1q9fv/2LSCPXqn4b6+/pzl4/3L7ry3ffvQsInn12U1OeA9tHLJq9\n7+HfexdgW0Ov3W23LrZs2Qz8KQD9/Vdx9NHHsWXL0PULDwOP1vxcKAJPsa9Nm64qtw7+fPPmq573\n+sItI3j+6CjqNvu57z0R33uo78e6snG5lkFk5ni8T0coT29tBE7OzJtqtl8NzMjMPxrkNacA145b\nIyVJmnz+JDOva/WbONJTIzOfjYj7gIXATQAREeXzLwzxspXAnwCPAP3j0ExJkiaLqcDLKb5LW86R\nnjoR8S7gauD9PLdk/R3AQZn5eBubJkmSxsCRnjqZeX1EzALOB2YD/wIcb+CRJKmzOdIjSZIqwdtQ\nSJKkSjD0SJKkSjD0ABFxVETcFBG/iohtEfH2uvK9I+LqsvzpiLglIg6oqzMlIi6LiL6I2BARN0TE\n3nV1XhgR10bE+ohYFxFXRsSe43GMzTbWPiv74gsR8VBEbIyIn0fE30bE9Lr9TJo+g+b8rtXV/8YQ\n+5k0/dasPouIIyPijoj4Tdkv34mIKTXlk6bPoGmfa/tGxHURsabst/siYlFdnUnTbxHx5xFxb0Q8\nFRFrI+KrEXHgIPXOj4hHy8+u26r8fdCMPhvP7wNDT2FPignLHwQGm+T0NYoldX8A/C7QA9weEXvU\n1LmY4h5dJwNHA/sBN9bt5zqKS5ouLOseDVzRrIMYZ2Pts/2AfYEPAwcD7wF+H7iybj+Tqc+gOb9r\nAETEUor7DAy2n8nUb2Pus4g4EvgGcCtwSPm4lIErBxYmU59Bc37XrgVeQdEfrwG+AlwfEa+rqTOZ\n+u0o4BLgcODNwO7AN+t+lz4OfIjiptSHAU9T3JS6q2Y/Vfo+aEafjd/3QWb6qHlQfAi+veb5q8pt\nB9VsC4pbU5xRPp8ObAL+qKbO3PJ1h5XP55XPX19T53hgC7BPu497vPtsiP28A3gG2KV8ftBk7bOx\n9hvPfUntPch+Jm2/NdpnwN3AeTvZ76TtszH22waKi8bV7quv5rNv0n6ulccyqzy+N9ZsexRYWvN8\nevm59a6a51X+Phh1nw2xn5Z8HzjSM7wpFP9L2n5N+yx6exPwxnLTIRTL/++oqbOa4ktp4EalRwDr\nMvP+mn3fXu778FY1vk1G0meDmQk8lZkD//s+kur0GYyw38r/QV0LfDAzHxtkP1Xqt2H7LCJeTHHc\nfRFxV3mq5jsR8Yaa/VSpz2Dk/0ZvAf64PK0QEfHu8rXfKcsn++faTIpjeQIgIvYH9mHHz/qngHt4\n7rO+6t8HjfTZUPtp+veBoWd4DwG/AD4TETMjoqscqnsJxXAcFNfz2Vz+RdZaS/GXTfnnDl9QmbmV\n4hdjHyaXkfTZDqK4NtKn2HGoskp9BiPvt2XAnZn59SH2U6V+G0mfvaL881yK36/jgW7gjoh4ZVlW\npT6Dkf+unQZMA35NEYi+SDGC8e9l+aTtt4gIitNUd2bmg+XmfSi+ZAe7KfXA8Vb2+2AMfVa/n5Z9\nHxh6hpGZW4A/Ag6k6NzfAMdQ/A9o205eWlmj7bOI+C3gZuBHwF+OX0snlpH0WzkZ9U0UVwqvvBH+\nrg18zv2PzPyHzPzXzPwwsBo4Y5ybPCGM4t/otRRzg94ELAAuAv5PRBw8rg1uj8uBVwPvbndDOsiY\n+6zV3wdekXkEyuG0+eVfRldm/joifgD8c1llDdAVEdPr0v3ssmygTv3s/V2BvWrqTBoj6DMAIuIF\nFPdceRJYVCb3AZXqMxhRvx1HMXKxvvhP1XZfiYh/ysw3UbF+G0Gf9ZZ/1t/qeRUwp/y5Un0Gw/db\nRBwEnAQcnJkDffdARBwNnEkxQXpS9ltEXAq8DTgqM3tritZQzH2azY4jF7OB+2vqVO77YIx9NrCP\nln8fONIzCpm5ofxgeBXFedv/WxbdRzGZauFA3YiYS/GBene56W5gZkS8vmaXCyl+Ge5pddvbZSd9\nNpDov0kxWe3tmbm57uWV7DPYab99Bngt8LqaB8DZwOnlz5Xst6H6LDMfoZhIObfuJQcCPy9/rmSf\nwU5/13ahOC2xte4lW3nuu2PS9Vv55f2HwHGZ2VNblpkPU3zB1n7WT6eYU/L9clPlvg+a0Gfj930w\nHrO5J/qDYvj2dRQrYrYB/6V8/tJ8bhb5McD+5V/sw8D1dfu4vNx+LMUw8F3A9+rq3AL8EDgUeAPF\n8Po17T7+dvQZ8FvADyiW1O5PkfoHHrtMxj5r1u/aIPvcYWXOZOu3Jv37PBtYR7GE+JXABRTLZvef\njH3WjH4DdgUepJi0fCjFCONHKL7Qj5+M/UbxOb6OYhl27WfS1Jo6H6OY4/QHwO9QhMT/RzFaVruf\nSnwfNKPPGMfvg7Z32ER4lP/wt1H8D6b28T/L8rMoZt73l7/I5wG71e1jCsW1Cvoolnn+H2Dvujoz\ngeXA+vKX5EvAtHYffzv6rHx9/WsH9jdnMvZZs37XBtnnVp4feiZNvzWrz8oP3p+X/z7vBI6crH3W\nrH6j+AK6nmKkbAPF6YhTJmu/DdFfW4FT6+qdV/bJRorTMQfUlVfm+6AZfcY4fh94w1FJklQJzumR\nJEmVYOiRJEmVYOiRJEmVYOiRJEmVYOiRJEmVYOiRJEmVYOiRJEmVYOiRJEmVYOiRJEmVYOiRJEmV\nYOiR1BEi4raIuHWQ7R+MiHURsV872iWpcxh6JHWK04HDIuLPBjZExP7A54AzM/PRVrxpROzaiv1K\nGn+GHkkdITN/CfwX4PMR8bJy81XArZl5HUBEHB0Rd0bExoh4JCIuiog9BvYREadGxA8jYkNE9EbE\nNRExq6Z8YURsi4jjI+K+iNgEHD6OhymphbzLuqSOEhFfAWYCXwE+Bbw6M5+IiAOB+4BPALcA+wCX\nAf+cmf+pfO0ZwC+BnwCzgWXAY5l5Ulm+ELgNuB/4KPAI8ERmrh+3A5TUMoYeSR0lIl4M/Bh4IbAo\nM/+x3P73wG8y86yausdShJg9MnPLIPs6ArgLmJaZm2pCz9sy83nzhyR1Nk9vSeoomfk4cAWwaiDw\nlF4HvLc8dbUhIjYAXwcCeBlARBwaEf8YET+PiKeA28vXvrT2LShGjCRNMru1uwGS1IAt5aPWCyhO\nZ11GEXRq9UTEbwG3AjcBpwCPAQdQBKOuuvpPN7vBktrP0CNpsugGDs7MhwcrjIh5FHOBPpGZa8tt\nbxjH9klqM09vSZosPgMcExF/GxGvjYgDIuKkiPjbsvznwLPA2RGxf0ScBPx521oradwZeiRNCpn5\nr8AxwEHAnRTzcj5NsVqLcnTnDODdFBOhPwx8pC2NldQWrt6SJEmV4EiPJEmqBEOPJEmqBEOPJEmq\nBEOPJEmqBEOPJEmqBEOPJEmqBEOPJEmqBEOPJEmqBEOPJEmqBEOPJEmqBEOPJEmqhP8PrnYVu1Mw\nC04AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x13153b610>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "years = genes_per_year.keys()\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar(years, [genes_per_year[y] for y in years])\n",
    "plt.xlabel(\"Year\")\n",
    "plt.ylabel(\"Genes referenced\")\n",
    "plt.savefig(\"genes_per_year.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAF5CAYAAAB6GeWJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XuYnWV97//3BzSkWAPaSMBd0qJoiNZTIger4oFWPKC1\noJVY6llbixTTWlGrW4TdaunWoELd/hSr5ZD9UzxWqSjaekAKlaGKEtKiyKgkwSkxICEJge/+43kG\nVpaTkFmzZtYc3q/rWlfWup97Pc933QyzvnM/9yFVhSRJkmCPQQcgSZI0XZgYSZIktUyMJEmSWiZG\nkiRJLRMjSZKklomRJElSy8RIkiSpZWIkSZLUMjGSJElqmRhJkiS1pkVilOTJST6X5KdJ7kryvDHq\nnJbkxiSbk3w5ycFdx/dKcnaSkSS3JrkwyX5ddR6Q5Pwkm5JsTPLhJPfrqnNgki8kuS3J+iRnJNmj\nq86jk3w9ye1Jbkjyl/1sD0mSNBjTIjEC7gf8B/CnwC9t3pbkFOB1wGuAw4DbgIuTzOuodibwHOA4\n4EjgwcAnu051AbAUOKqteyTwwY7r7AFcBNwHOAJ4KfAy4LSOOvcHLgauB5YBfwmcmuRVvXxwSZI0\nfWS6bSKb5C7g+VX1uY6yG4G/q6pV7esFwAbgpVX18fb1z4Djq+rTbZ0lwBrgiKq6IslS4PvA8qq6\nqq1zNPAF4Neran2SZwGfAw6oqpG2zh8D7wIeVFXbk7wWOB3Yv6q2t3XeCfxeVT1ikptHkiRNounS\nY7RTSQ4C9ge+MlpWVbcAlwNPaIseT9PL01lnLTDcUecIYONoUtS6hKaH6vCOOlePJkWti4F9gEd2\n1Pn6aFLUUWdJkn16/JiSJGkamPaJEU1SVDQ9RJ02tMcAFgHb2oRpZ3X2B27qPFhVdwI3d9UZ6zqM\ns84OkuydZFmSvcc6LkmSxjbV36H3mYqLiMcClwJDSX7RdeyLND1OkiTNdUcDz+wq+1WaMb1PBL41\n2QHMhMRoPRCaXqHOnppFwFUddeYlWdDVa7SoPTZap3uW2p7AA7vqHNp1/UUdx0b/XXQvdbr9Zvvv\nsjGOHQn8zU7eJ0mSGr+JiRFU1fVJ1tPMJPsu3D34+nDg7LbalcD2tk7n4OvFwGVtncuAfZM8rmOc\n0VE0SdflHXXekmRhxzijZwCbgGs66vyvJHu2t+JG66ytqk07+Rg/AjjvvPNYunTp+Bthjlq5ciWr\nVq0adBgzju02frZZb2y38bPNxm/NmjWccMIJ0H6XTrZpkRi1awkdTJOkADwkyWOAm6vqxzRT8d+a\n5Dqahjkd+AnwWWgGYyc5B3hPko3ArcD7gEur6oq2zrVJLgY+1M4smwe8H1hdVaM9PV+iSYDObZcI\nOKC91llVdUdb5wLgfwIfSfK3wKOAPwNO3sVH3AKwdOlSli0bq9NIY9lnn31srx7YbuNnm/XGdhs/\n22xCtkzFRaZFYkQzq+xfaAZZF/DutvxjwCuq6ox20NUHgX2BbwDPqqptHedYCdwJXAjsRTN258Su\n67wYOItmNtpdbd27E5qquivJMcAHaLrrbgM+Cry9o84tSZ5B01v1bWAEOLWqzplYE0iSpEGbFolR\nVX2Ne5khV1WnAqfu4vhW4KT2sbM6PwdOuJfr/Bg45l7qfA94yq7qSJKkmWcmTNeXJEmaEiZGmrZW\nrFgx6BBmJNtt/Gyz3thu42ebTX/TbkuQ2SjJMuDKK6+80kF3kiSNw9DQEMuXL4dmS6+hyb6ePUaS\nJEktEyNJkqSWiZEkSVLLxEiSJKllYiRJktQyMZIkSWqZGEmSJLVMjCRJklomRpIkSS0TI0mSpJaJ\nkSRJUsvESJIkqWViJEmS1DIxkiRJapkYSZIkte4z6AAkSdLgDQ8PMzIyAsDChQtZvHjxgCMaDBMj\nSZLmuOHhYZYsWcqWLZsBmD9/b9auXdNzcjSTkywTI0mS5riRkZE2KToPgC1bTmBkZKSnhKbfSdZU\nc4yRJElqLW0fvdsxyTqPLVs23917NBPYYyRJkibBxBKsQbHHSJIkqWViJEmS1DIxkiRJapkYSZIk\ntUyMJEmSWiZGkiRJLRMjSZKklomRJElSy8RIkiSpZWIkSZLUMjGSJElqmRhJkiS1TIwkSZJaJkaS\nJEktEyNJkqSWiZEkSVLLxEiSJKllYiRJktS6z6ADkCRJU294eJiRkREA1q1bN+Bopg8TI0mS5pjh\n4WGWLFnKli2bAZg3b/6AI5o+TIwkSZpjRkZG2qToPAC2bTthwucc7YFas2bNhM81SCZGkiTNWUv7\ncpbuHqiZzMHXkiRpQnbsgTp90OFMiD1GkiSpT/rTAzVI9hhJkiS1TIwkSZJaJkaSJEktEyNJkqTW\njEiMkuyZ5J1Jrk+yOcl1Sd46Rr3TktzY1vlykoO7ju+V5OwkI0luTXJhkv266jwgyflJNiXZmOTD\nSe7XVefAJF9IcluS9UnOSDIj2lKSJO3cTPky/yvglcBrgUOANwJvTPK60QpJTgFeB7wGOAy4Dbg4\nybyO85wJPAc4DjgSeDDwya5rXUAzrP6otu6RwAc7rrMHcBHNjL4jgJcCLwNO68snlSRJAzNTpusf\nCny2qr7Yvh5O8mKaBGjUycDpVfV5gCQvATYAzwc+nmQB8Arg+Kr6Wlvn5cCaJIdV1RVJlgJHA8ur\n6qq2zknAF5K8oarWt8cPAZ5WVSPA1UneBrwryalVtX1SW0KSpGlgtu61NlMSo38G3pDkYVX1X0ke\nAzwRWAmQ5CBgf+Aro2+oqluSXA48Afg48Hiaz9tZZ22S4bbOFTQ9QBtHk6LWJUABhwOfbetc3SZF\noy4GPgA8EvhOPz+4JEnTzWzea23CiVGSPYFHATdU1caJh/TLqurvkxwIrE2yneYW4F9V1f9tq+xP\nk7xs6HrrhvYYwCJgW1Xdsos6+wM3dV37ziQ3d9UZ6zqjx0yMJEmz2mTstTZdjDsxSnImTY/JOW1S\n9DXgt4HNSY6pqn/tc4wk+TOasTwvAq4BHgu8N8mNVXVuv683WVauXMk+++yzQ9mKFStYsWLFgCKS\nJGki+rvS9erVq1m9evUOZZs2berrNe5NLz1GL2A0RYTnAgfRjLn5I+CvaW5x9dtbgHdU1Sfa199P\n8pvAm4FzgfVAaHqFOntzFgGjt8XWA/OSLOjqNVrUHhut0z1LbU/ggV11Du2Kb1HHsZ1atWoVy5Yt\n21UVSZKmpakYUzRWZ8HQ0BDLly+flOuNpZfEaCH3JADPBj5RVf+Z5CM0A6Anwx7AnV1ld7XlVNX1\nSdbTzCT7LkA72Ppw4Oy2/pXA9rbOp9s6S4DFwGVtncuAfZM8rmOc0VE0SdflHXXekmRhxzijZwCb\naHqzJEmaVWbzmKJuvSRGG4BHJFkHPJNmCj3A3vxy8tIvnwHemuQnwPeBZTQDrz/cUefMts51wI9o\ntvf9Cc2A6dHB2OcA70myEbgVeB9waVVd0da5NsnFwIeSvBaYB7wfWN3OSAP4Ek0CdG67RMAB7bXO\nqqo7JunzS5I0MLN5TFG3XhKjf6CZ5bWOZsDzJW354cC1fYqr2+uBdwBn0dy2upFmFtjpoxWq6owk\ne9OsObQv8A3gWVW1reM8K2mStwuBvYAvAid2XevF7XUuoemVupCOnrCquivJMe31v0WzXtJHgbf3\n56NKkjRd9XdM0XQ07sSoqk5N8j3gQJrbaFvbQ3cC7+pncB3X3Az8ZfvYZWzAqbs4vhU4qX3srM7P\ngV2mwlX1Y+CYXdWRJEkzT0/T9avqQoAk8zvKPtavoCRJkgZh3FuCtPuWvS3JT4FfJHlIW356klf2\nPUJJkqQp0steaX9FszfYG4HO8TvfA17Vh5gkSZIGopfE6CXAa6rqfHachfYdmvWMJEmSZqReEqP/\nAVy3k3Pdd2LhSJIkDU4vidE1wJPHKH8B96wyLUmSNOP0MivtNOBjSf4HTWJ1bLuC9EtwCrskSZrB\nelnH6LNJngv8T5rFDU8DhoDnVtWX+xyfJEkagHXr1jE0NHT387mi13WMvgH8bp9jkSRJ08Sxx76Q\nbdtuB2b33mjdekqMAJLMo9mJfodxSlU1PNGgJEnSYDVJ0ezfG63buBOjJA8DPgL8dvchmr3T9uxD\nXJIkaeBm/95o3XrpMfoosJ1moPXoRrKSJEkzXi+J0WOB5VV1bb+DkSRJGqRe1zFa2O9AJEmSBq2X\nxOgU4IwkT03ya0kWdD76HaAkSdJU6eVW2iXtv1/pKnfwtSRJmtF6SYye1vcoJEmSpoFeVr7+2mQE\nIkmSNGi9jDEiyZOTnJfkW+2eaST5oyRP6m94kiRJU2fciVGS44CLgduBZcBe7aF9gLf0LzRJkqSp\n1UuP0VuBP6mqVwN3dJRfSpMoSZIkzUi9JEZLgK+PUb4J2Hdi4UiSJA1OL4nReuDgMcqfBPxwYuFI\nkiQNTi+J0YeA9yY5nGbdogcn+UPgfwMf6GdwkiRJU6mXdYzeRZNQfQXYm+a22lbgf1fV+/sYmyRJ\n0pTqZR2jAv46yd/R3FL7VeCaqvpFv4OTJEmaSr30GAFQVdtoNpSVJEmaFcadGCX5NM3Yom4FbAGu\nAy6oqrUTjE2SJHUYHh5mZGQEgIULF7J48eIBRzT79NJjtAl4PvBz4Mq2bBnNVP0vAS8CTklyVFVd\n2pcoJUma44aHh1myZClbtmwGYP78vVm7do3JUZ/1Mivtp8AFwEOq6riqOg54KHAezXT9pcDHgL/t\nW5SSJM1xIyMjbVJ0HnAeW7Zsvrv3SP3TS2L0auDMqrprtKB9/n7g1e3g7LOA3+pPiJIk6R5L24cm\nQy+J0X2BQ8YoPwTYs32+hbHHIUmSJE1bvYwxOhc4J8nfAP/elh1Ks4HsP7avnwJ8f+LhSZIkTZ1e\nEqOVwAbgjcCitmwDsIp7xhV9CfjihKOTJEmaQr0s8Hgn8Nc0izwuaMtu6aoz3J/wJElSP3RO9V+3\nbt2Ao5m+el7gEX45IZIkSdNP91T/efPmDzii6WtCiZEkSZr+dpzqD9u2nTDYgKYxEyNJkuYMp/nf\nm16m60uSJM1KJkaSJEmtXjaR/bOdHOrcRPbr7ew1SZKkGaPXdYweBOwNbGzLHgBsBn4B7Af8MMnT\nqurHfYlSkiRpCvRyK+0UmhWvH1ZVv1ZVvwY8HLgceD2wGFhPs+CjJEnSjNFLj9HfAC+oqh+MFlTV\ndUneAHyyqh6S5I3AJ/sVpCRJ0lTopcfowYydUN0H2L99fiNw/16DkiRJGoReEqN/AT6Y5HGjBe3z\nDwBfbYseBVw/8fAkSZKmTi+J0SuBm4Erk2xNshX4dlv2yrbOL4C/6E+IkiRJU6OXTWTXA7+b5BCa\nQdcAa6tqbUedf+lTfJIkSVOm5y1Bqupa4No+xiJJkjRQvSzwuCfwMuAomjWLdrgdV1VP70tkkiRJ\nU6yXHqP30iRGXwC+R7PitSRJ0ozXy+Dr44E/qKoXVdXrq2pl56PfAY5K8uAk5yYZSbI5yXeSLOuq\nc1qSG9vjX05ycNfxvZKc3Z7j1iQXJtmvq84DkpyfZFOSjUk+nOR+XXUOTPKFJLclWZ/kjCTuOydJ\n0gzXy5f5Npr90KZMkn2BS4GtwNHAUppZbxs76pwCvA54DXAYcBtwcZJ5Hac6E3gOcBxwJM2aTN0L\nUV7Qnv+otu6RwAc7rrMHcBFNb9sRwEtpetBO68dnlSRJg9NLYvRu4OQk6Xcwu/AmYLiqXlVVV1bV\nDVV1SVV1rpV0MnB6VX2+qr4HvIQm8Xk+QJIFwCuAlVX1taq6Cng58MQkh7V1ltIkXq+sqm9X1beA\nk4Djk4wuXnk0cAjwh1V1dVVdDLwNODFJz4PZJUnS4PWSGD0J+EPgB0n+KcmnOh99jm/Uc4FvJ/l4\nkg1JhpK8avRgkoNoVt3+ymhZVd1Cs3/bE9qix9P08nTWWQsMd9Q5AtjYJk2jLqEZR3V4R52rq2qk\no87FwD7AIyf6QSVJ0uD00sPxc+DT/Q7kXjwEeC1Nb9Vf09wqe1+SrVV1Lk1SVMCGrvdt4J5tShYB\n29qEaWd19gdu6jxYVXcmubmrzljXGT32nfF9NEmSGsPDw4yMNH93L1y4kMWLF/flXOvWretLfHNB\nLws8vnwyArkXewBXVNXb2tffSfJbwJ8A5w4gnp6sXLmSffbZZ4eyFStWsGLFigFFJEmaLoaHh1my\nZClbtmwGYP78vVm7dk1PyVH3uebNm9/XWCfL6tWrWb169Q5lmzZtmtIYZsqYmHXAmq6yNcCx7fP1\nQGh6hTp7cxYBV3XUmZdkQVev0aL22Gid7llqewIP7KpzaFcsizqO7dSqVatYtmzZrqpIkuaokZGR\nNpE5D4AtW05gZGSkp8So+1zbtp3Qx0gnz1idBUNDQyxfvnzKYtitxCjJEHBUVW1MchW7WLuoqibj\nm/9SYElX2RLghvaa1ydZTzOT7LttzAtoxgWd3da/Etje1vl0W2cJsBi4rK1zGbBvksd1jDM6iibp\nuryjzluSLOwYZ/QMYBNwTV8+rSRpDls6Tc81N+xuj9FnaabKA3xmkmLZlVXApUneDHycJuF5FfDq\njjpnAm9Nch3wI+B04Cc0sVNVtyQ5B3hPko3ArcD7gEur6oq2zrVJLgY+lOS1wDzg/cDqdo84gC/R\nJEDntksEHNBe66yqumOyGkCSJE2+3UqMquodYz2fKlX17SS/D7yLZmr89cDJVfV/O+qckWRvmjWH\n9gW+ATyrqrZ1nGolcCdwIbAX8EXgxK7LvRg4i2Y22l1t3ZM7rnNXkmOADwDfolkv6aPA2/v1eSVJ\n0mD0PMaoXThxrL3Shica1Fiq6iKahRV3VedU4NRdHN9Ksy7RSbuo83Nglzdjq+rHwDG7qiNJkmae\nXjaRfThwDvDb3Ydoxh7t2Ye4JEnSvVi3bh1DQ0PAxKf3q9FLj9E/0AxiPoZmtpibyEqSNADHHvtC\ntm27HZjY9H7do5fE6LHA8qq6tt/BSJKk3dckRROf3q979JIYXQMs7HcgkiSpF07J76de9ko7BTgj\nyVOT/FqSBZ2PfgcoSZI0VXrpMbqk/fcrXeUOvpYkSTNaL4nR0/oehSRJ0jTQyyayX5uMQCRJkgat\npwUek+wLvJJ7Rnx9H/hIVU3tFriSJEl9NO7B10keD/yAZnuNB7aPPwd+kMSt4yVJ0ozVS4/RKuBz\nwKurajtAkvsAH6bZyPXI/oUnSZI0dXpJjB5PR1IEUFXbk5wBfLtvkUmSJE2xXtYxugUYa1nNA4Fb\nJxaOJEnS4PSSGP3/wDlJXpTkwPZxPM2ttNX9DU+SJGnq9HIr7Q00Czn+Y8f77wA+ALypT3FJkiRN\nuV7WMdoGnJzkzcBD2+IfVNXmvkYmSZI0xXpaxwigTYSu7mMskiRJA7VbiVGSTwEvq6pb2uc7VVXH\n9iUySZKkKba7PUabaMYVQTMrrXZRV5IkaUbarcSoql7e8fxlkxaNJEnSAPWyJchX273SussXJPlq\nf8KSJEmaer2sY/RUYN4Y5fOBJ08oGkmSpAHa7VlpSR7d8fIRSfbveL0n8Ezgp/0KTJIkaaqNZ7r+\nf9AMui5grFtmtwMn9SMoSZKkQRhPYnQQEOCHwGHAzzqObQNuqqo7+xibJEnSlNrtxKiqbmif9jIu\nSZIkadrraeXrJA8DngbsR1eiVFWn9SEuSZKkKTfuxCjJq2k2jB0B1rPjYo8FmBhJkqQZqZceo7cC\nf1VVf9vvYCRJkgapl/FCDwA+0e9AJEmSBq2XxOgTwDP6HYgkSdKg9XIr7Trg9CRHAFcDd3QerKr3\n9SMwSZKkqdZLYvQa4BfAU9pHpwJMjCRJGoB169YxNDTEmjVrBh3KjDXuxKiqDpqMQCRJ0sQce+wL\n2bbt9kGHMaP1vFhjknlJliTpaS0kSZLUX01SdB5w+qBDmbHGnRgl2TvJOcBm4PvA4rb8/Une1Of4\nJEnSuCyl2cVLveilx+idwGOApwJbOsovAV7Uh5gkSZIGopfbYM8HXlRV/5akc9Xr7wMP7U9YkiRJ\nU6+XHqMHATeNUX4/dtweRJIkaUbpJTH6NvCcjtejydCrgMsmHJEkSdKA9HIr7S3APyd5RPv+k9vn\nv80vr2skSZI0Y4y7x6iqvgk8liYpuppme5CbgCdU1ZX9DU+SJGnq9LQGUVX9AHh1n2ORJEkaqJ4X\nZ0yyH7AfXb1OVfXdiQYlSZI0CONOjJIsBz5Gs4JUug4XsGcf4pIkSZpyvfQYfQT4T+CVwAacoi9J\nkmaJXhKjhwDHVdV1/Q5GkiTdY926dQwNDd39XJOvl8ToX2i2BDExkiRpEh177AvbjWFh3rz5A45m\nbuglMXoVcG6S3wK+B9zRebCqPtePwCRJmuuapOi89vkJgw1mjuglMToUOAz43TGOOfhakqS+Wjro\nAOaUXrYEOQs4HzigqvboepgUSZKkGauXxGghsKqqNvQ7mN2V5E1J7krynq7y05LcmGRzki8nObjr\n+F5Jzk4ykuTWJBe26zF11nlAkvOTbEqyMcmHk9yvq86BSb6Q5LYk65OckaSXtpQkSdNIL1/mnwKe\n1u9AdleSQ4HXAN/pKj8FeF177DDgNuDiJPM6qp1JswHuccCRwIOBT3Zd4gKafsuj2rpHAh/suM4e\nwEU0tyGPAF4KvAw4rR+fT5IkDU4vY4z+E3hnkifR7JXWPfj6ff0IbCxJfpVmFNqrgLd1HT4ZOL2q\nPt/WfQnNOkvPBz6eZAHwCuD4qvpaW+flwJokh1XVFUmWAkcDy6vqqrbOScAXkryhqta3xw8BnlZV\nI8DVSd4GvCvJqVW1fbI+vyRJmly99Bi9CvgF8BSaHpqVHY/X9y+0MZ0N/FNVfbWzMMlBwP7AV0bL\nquoW4HLgCW3R42kSwc46a4HhjjpHABtHk6LWJTSDyg/vqHN1mxSNuhjYB3jkRD6cJEkarHH3GFXV\nQZMRyL1JcjzwWJoEp9v+NMlL97inDe0xgEXAtjZh2lmd/YGbOg9W1Z1Jbu6qM9Z1Ro99B0mSdLfO\nhSoXLlzI4sWLBxzRzvW8iexUSvLrNOODfqeq7ri3+pIkafroXKhy/vy9Wbt2zbRNjmZEYgQsBx4E\nDCUZ3bh2T+DIJK+jGfMTml6hzt6cRcDobbH1wLwkC7p6jRa1x0brdM9S2xN4YFedQ7viW9RxbKdW\nrlzJPvvss0PZihUrWLFixa7eJknSjNa5UOWWLScwMjIyZmK0evVqVq9evUPZpk2bpiLEu82UxOgS\n4FFdZR8F1gDvqqofJllPM5PsuwDtYOvDacYlAVwJbG/rfLqtswRYDFzW1rkM2DfJ4zrGGR1Fk3Rd\n3lHnLUkWdowzegawCbhmVx9i1apVLFu2bBwfW5Kk2eLeF6ocq7NgaGiI5cuXT1ZQv2RGJEZVdRtd\nSUeS24D/rqo1bdGZwFuTXAf8CDgd+Anw2fYctyQ5B3hPko3ArcD7gEur6oq2zrVJLgY+lOS1wDzg\n/cDqdkYawJfaWM5tlwg4oL3WWd7mkyRpZpsRidFO1A4vqs5IsjfNmkP7At8AnlVV2zqqrQTuBC4E\n9gK+CJzYdd4X06zufQlwV1v35I7r3JXkGOADwLdo1kv6KPD2fn0wSZI0GONOjJI8E/hFVX2zfX0i\n8GqaXpQTq2pjf0McW1U9fYyyU4FTd/GercBJ7WNndX4O7HKnvqr6MXDMboYqSZJmiF7WMfo7YAFA\nkkcB76ZZCfog4D27eJ8kSdK01suttIO4Z7zPccDnq+otSZbRJEiSJEkzUi89RtuAvdvnv0MzGBng\nZtqeJEmSpJmolx6jb9LM7LqUZrPWF7XlD6eZBSZJkjQj9dJj9Dqa9YBeALy2qn7alj+LZpaXJEnS\njNTLXmnDjDEjq6pW9iUiSZKkAelpHaMkewAH02yfsUOvU1V9vQ9xSZIkTble1jE6ArgA+A2arTI6\nFc0eZpIkSTNOLz1G/wf4NvAcYB1dK1BLkiTNVL0kRg8DXlBV1/U7GEmSpEHqZVba5TTjiyRJkmaV\nXnqM3g+8O8n+wNXADjvKV9V3+xGYJEnSVOslMfpk++9HOsqKZiC2g68lSdKM1eteaZIkSbNOLws8\n3jAZgUiSJA1arws8PhR4PbC0LboGeG9V/aBfgUmSJE21cc9KS3I0TSJ0GPDd9nE48P0kv9vf8CRJ\nkqZOLz1G7wJWVdWbOguTvAv4W+DL/QhMkiRpqvWyjtFS4Jwxyj8CPGJi4UiSJA1OL4nRz4DHjlH+\nWOCmiYUjSZI0OL3cSvsQ8P8leQjwrbbsicApwHv6FZgkSdJU6yUxOh24FfgL4J1t2Y3AqcD7+hOW\nJEnS1OtlHaMCVgGrkty/Lbu134FJkiRNtZ7WMRplQiRJkmaT3UqMkgwBR1XVxiRX0eyJNqaqWtav\n4CRJkqbS7vYYfRbY2vF8p4mRJEnSTLVbiVFVvaPj+amTFo0kSdIA9bIlyA+T/NoY5fsm+WF/wpIk\nSZp6vSzw+JvAnmOU7wX8+oSikSRJGqDdnpWW5HkdL49Osqnj9Z7AUcD1/QpMkiRpqo1nuv5n2n8L\n+FjXsTuAH9Es+ihJkjQj7XZiVFV7ACS5Hji0qkYmLSpJkqQB6GXl64MmIxBJkqRB62nl6yT3A54C\nLAbmdR6rKvdLkyRJM9K4E6MkjwMuAvYG7gfcDCwENgM34UaykiRphupluv4q4J+ABwC3A0cAvwFc\nCbyhf6FJkiRNrV5upT0W+OOquivJncBeVfXDJG+kma32qb5GKEnSLDU8PMzISDOXad26dQOORtBb\nYnQHcFf7/CaacUZrgE3AgX2KS5KkWW14eJglS5ayZctmAObNmz/giAS9JUZXAYcC/wV8DTgtyULg\nj4Dv9TE2SZJmrZGRkTYpOg+AbdtOGGxAAnobY/QWYLS/76+AjcAHgAcBf9ynuCRJmiOWtg9NB72s\nY/Ttjuc3Ac/sa0SSJEkDMu4eoyRfTbLvGOULkny1P2FJkiRNvV5upT2VrkUdW/OBJ08oGkmSpAHa\n7VtpSR7d8fIRSfbveL0nzS21n/YrMEmSpKk2njFG/wFU+xjrltntwEn9CEqSJGkQxpMYHQQE+CFw\nGPCzjmPzllXyAAASfklEQVTbgJuq6s4+xiZJkjSldjsxqqob2qe9jEuSJEma9nYrMUryPOCfq+qO\n9vlOVdXn+hKZJEnSFNvdHqPPAPvTbAHymV3UK5qB2JIkSTPObiVGVbXHWM8lSZJmE5McSZKk1m4n\nRkmenuSaJAvGOLZPku8nObq/4d19/jcnuSLJLUk2JPl0koePUe+0JDcm2Zzky0kO7jq+V5Kzk4wk\nuTXJhUn266rzgCTnJ9mUZGOSDye5X1edA5N8IcltSdYnOSOJSaYkSTPceL7MXw98qKpu6T5QVZuA\nDzJ56xg9GXg/cDjwO8B9gS8l+ZXRCklOAV4HvIZmOYHbgIuTdK7SfSbwHOA44EjgwcAnu651Ac1u\nfke1dY+k+Wyj19kDuIjmNuQRwEuBlwGn9eWTSpJmreHhYYaGhhgaGmLdunX3/gZNufGsY/QY4JRd\nHP8S8IaJhTO2qnp25+skL6MZCL4c+GZbfDJwelV9vq3zEmAD8Hzg421P1yuA46vqa22dlwNrkhxW\nVVckWQocDSyvqqvaOicBX0jyhqpa3x4/BHhaVY0AVyd5G/CuJKdW1fbJaANJ0sw2PDzMkiVL2bJl\nMwDz5s0fcEQay3h6jBYBd+zi+HbgQRMLZ7ftSzMD7maAJAfRzJr7ymiFtmfrcuAJbdHjaRLBzjpr\ngeGOOkcAG0eTotYl7bUO76hzdZsUjboY2Ad4ZB8+myRpFhoZGWmTovOA89i2bcugQ9IYxpMY/RT4\nrV0cfzQw6f2CSUJzS+ybVXVNW7w/TfKyoav6hvYYNIndtjFuBXbWGV2S4G7tat43d9UZ6zp01JEk\naSeWtg9NR+NJjC4CTk/yS31/7VifdwCf71dgu/D3wCOA46fgWpIkaQ4Zzxij/wUcC/xnkrOAtW35\nIcCJNAs7/nV/w9tRe91nA0+uqs7eqfU0+7gtYsfenEXAVR115iVZ0NVrtKg9Nlqne5bansADu+oc\n2hXaoo5jO7Vy5Ur22WefHcpWrFjBihUrdvU2SZLmhNWrV7N69eodyjZt2jSlMYxnr7QNSX4b+ADw\nTppEBJpbWBcDJ1ZV9y2mvmmTot8DnlJVw12xXZ9kPc1Msu+29RfQjAs6u612Jc04qKOAT7d1lgCL\ngcvaOpcB+yZ5XMc4o6NoPuvlHXXekmRhxzijZwCbgNFbe2NatWoVy5YtG+9HlyRpThirs2BoaIjl\ny5dPWQzj6TEa3Uj22UkeABxMkzD8V1VtnIzgRiX5e2AF8DzgtiSjPTSbqmp09NqZwFuTXAf8CDgd\n+Anw2Tb2W5KcA7wnyUbgVuB9wKVVdUVb59okFwMfSvJaYB7NMgGr2xlp0My+uwY4t10i4ID2WmdV\n1a4Gp0uSpGluXInRqDYR+vc+x7Irf0LTM/WvXeUvB/6xjemMJHvTrDm0L/AN4FlVta2j/krgTuBC\nYC/gizS3ATu9GDiLZjbaXW3dk0cPVtVdSY6h6Tn7Fs16SR8F3j7BzyhJ0pywbt06hoaGAFi4cCGL\nFy8ecET36Ckxmmq7uz9bVZ0KnLqL41tpFqHc6UKUVfVz4IR7uc6PgWN2JyZJkrSjY499Idu23Q7A\n/Pl7s3btmmmTHLmNhSRJmlJNUtSs57Rly2ZGRkbu7S1TZkb0GEmSpNlmeq7lZI+RJElSy8RIkiSp\n5a00SZImyfDw8N3jZ9atm/Rds9QHJkaSJE2C4eFhlixZ2m4cC/Pm/dKOWpqGTIwkSZoEIyMjbVJ0\nHgDbtu1yJRhNEyZGkiRNquk5+0pjMzGSJKmPRscVrVmzZtChqAcmRpIk9Un3uCLNPE7XlySpT3Yc\nV3T6oMNRD+wxkiSp7xxXNFPZYyRJktQyMZIkSWp5K02SpHHoXM164cKFLF68eMARqZ9MjCRJ2k3d\ns87mz9+btWvXmBzNIt5KkyRpN+046+w8tmzZfHfvkWYHe4wkSRo3Z53NVvYYSZIktUyMJEmSWt5K\nkyRpAtatW8fQ0NDdzzWzmRhJkjQBxx77QrZtux2AefPmDzgaTZSJkSRJE9AkRee1z08YbDCaMBMj\nSZImzFlqs4WDryVJklomRpIkSS0TI0mSpJaJkSRJUsvB15IkdRgeHr57/7OFCxe6QewcY2IkSVJr\neHiYJUuWthvFwvz5e7N27RqToznExEiSpNbIyEibFDXrEm3ZcgJXX3313T1Irmw9+5kYSZL0S+5Z\nl8iVrecWEyNJknbBla3nFhMjSZLulStbzxVO15ckSWqZGEmSJLVMjCRJklomRpIkSS0TI0mSpJaz\n0iRJc1rnFiAu4CgTI0nSnNW9BYgLOMrESJI0p3T3EHVuAeICjjIxkiTNGTvvIXIBRzVMjCRJc0b3\nJrH2EKmbiZEkaVYbe3C1PUQam4mRJGnWcnC1xsvESJI064z2Eq1Zs8ZbZxoXEyNJ0ozXebts69at\nPP3pv3N3L1HDW2faPSZGkqQZbazbZdu2baHpJboeeNsgw9MMY2IkSZpxdm8tInuJNH4mRpKkaW9X\nt8pci0j95CayPUpyYpLrk9ye5N+SHDromGab1atXDzqEGcl2Gz/brDeT2W7Dw8MMDQ0xNDTEZZdd\nxpIlS1m+fDnLly/nqU99ekcP0XntbTOpP0yMepDkRcC7gbcDjwO+A1ycZOFAA5tl/LLqje02frZZ\nb/rZbr0lQkuxl0j95q203qwEPlhV/wiQ5E+A5wCvAM4YZGCSNBPc262xewZPO2ZIU8vEaJyS3BdY\nDvzNaFlVVZJLgCcMLDBJGrDOZGfhwqYDvTP52Wuvve5+biKk6crEaPwWAnsCG7rKNwBLej1p9y+U\nxYsX9xygNFkm+nO6q/ff27l39aU73tedX9ILFy7k9ttvZ2hoqC/nHu+1J/Nak33tznbrTnb22ms+\nELZuvZ3GnsCdgImQpjcTo6kxH2DNmjVjHly3bh3HHvuCu++bz5s3n0996kIOOOCAqYtwGtq0adPd\nv3S1+yar3Sb6c7qr99/bubuP3/e+84Bwxx1be3rdDK+86+5j27dvZ/ny5X0593ivPZnXmuxrd7bb\nfe4zj+3btwGvBGDr1nPaa74SGAE+e/exbdtGj13Pji6ahNc3TuK5vXa/rr2z78euY1Oyn0uqaiqu\nM2u0t9I2A8dV1ec6yj8K7FNVvz/Ge14MnD9lQUqSNPv8YVVdMNkXscdonKrqjiRXAkcBnwNIkvb1\n+3bytouBPwR+BDivVJKk3Tcf+E2a79JJZ49RD5L8AfBR4E+AK2hmqb0AOKSqfjbA0CRJ0gTYY9SD\nqvp4u2bRacAi4D+Ao02KJEma2ewxkiRJarnytSRJUsvESJIkqWVitJuSPDnJ55L8NMldSZ7XdXy/\nJB9tj9+W5KIkB3fV2SvJ2UlGktya5MIk+3XVeUCS85NsSrIxyYeT3G8qPmO/TbTN2rZ4X5Jrk2xO\nckOS9yZZ0HWeWdNm0J+fta76/7yT88yadutXmyV5QpKvJPlF2y7/mmSvjuOzps2gb7/XDkhyQZL1\nbbtdmeTYrjqzpt2SvDnJFUluSbIhyaeTPHyMeqclubH93fXlufx90I82m8rvAxOj3Xc/mkHWfwqM\nNTDrszTTCZ8LPBYYBi5J8isddc6k2VPtOOBI4MHAJ7vOcwHN8q9HtXWPBD7Yrw8xxSbaZg8GDgD+\nHHgk8FLgmcCHu84zm9oM+vOzBkCSlTTLDY91ntnUbhNusyRPAP4Z+CLw+PZxFqMrIjZmU5tBf37W\nzgceQtMevwV8Cvh4ksd01JlN7fZk4P3A4cDvAPcFvtT1s3QK8DrgNcBhwG00G43P6zjPXPo+6Eeb\nTd33QVX5GOeD5hfl8zpeP6wtO6SjLDTbhLyifb0A2Ar8fkedJe37DmtfL21fP66jztHAdmD/QX/u\nqW6znZznBcDtwB7t60Nma5tNtN2454tsvzHOM2vbrdc2Ay4DTt3FeWdtm02w3W6lWXiv81wjHb/7\nZu3vtfazLGw/35M6ym4EVna8XtD+3vqDjtdz+ftg3G22k/NMyveBPUb9sRfNX1uj6+5TzX+RrcCT\n2qLH0yyP8JWOOmtpvrhGN589AthYVVd1nPuS9tyHT1bwA7I7bTaWfYFbqmr0r/gnMHfaDHaz3dq/\nxM4H/rSqbhrjPHOp3e61zZI8iOZzjyS5tL0t9K9JnthxnrnUZrD7/49eBLyovYWRJMe37/3X9vhs\n/722L81nuRkgyUHA/uz4u/4W4HLu+V0/178PemmznZ2n798HJkb9cS3wY+CdSfZNMq/tFvx1mq4/\naNY72tb+x+60geYHgvbfHb7EqupOmh+e/ZlddqfNdpBm7ai3smO36FxqM9j9dlsFfLOqPr+T88yl\ndtudNntI++/baX6+jgaGgK8keWh7bC61Gez+z9rLgL2B/6ZJmj5A0xPyw/b4rG23JKG5JfbNqrqm\nLd6f5ot4rI3GRz/vnP0+mECbdZ9n0r4PTIz6oKq2A78PPJzmP8AvgKfQ/CV11y7eOmeNt82S3B/4\nAvA94B1TF+n0sjvt1g6gfTrNiuxz3m7+rI3+Lvw/VfWPVfWdqvpzYC3wiikOeVoYx/+j59OMVXo6\nsBx4D/CJJI+c0oAH4++BRwDHDzqQGWTCbTbZ3weufN0nbdfdsvY/2Lyq+u8k/wb8e1tlPTAvyYKu\nvxIWtcdG63TPStgTeGBHnVljN9oMgCS/SrNHzs+BY9u/AEbNqTaD3Wq3p9H0gGxq/ji726eSfL2q\nns4ca7fdaLN17b/dW3yvARa3z+dUm8G9t1uSQ4DnA4+sqtG2uzrJkcCJNIO6Z2W7JTkLeDbw5Kpa\n13FoPc1YrEXs2AOyCLiqo86c+z6YYJuNnmPSvw/sMeqzqrq1/eXxMJr7yJ9pD11JMwDsqNG6SZbQ\n/NK9rC26DNg3yeM6TnkUzQ/M5ZMd+6Dsos1G/zL4Es0Au+dV1baut8/JNoNdtts7gUcDj+l4AJwM\nvLx9PifbbWdtVlU/ohn8uaTrLQ8Hbmifz8k2g13+rO1Bcwvkzq633Mk93y+zrt3aL/jfA55WVcOd\nx6rqepov4c7f9Qtoxrh8qy2ac98HfWizqfs+mIoR6LPhQdNV/BiamT53Aa9vXx9Y94yOfwpwUPsf\n/3rg413n+Pu2/Kk0Xc6XAt/oqnMR8G3gUOCJNF355w768w+izYD7A/9GM534IJq/HkYfe8zGNuvX\nz9oY59xhxtFsa7c+/f95MrCRZvr0Q4HTaaYMHzQb26wf7QbsCVxDM9D6UJqeyr+g+dI/eja2G83v\n8Y00U9A7fyfN76jzRpoxV88FHkWTSP4XTa9b53nmxPdBP9qMKfw+GHiDzZRH+8vhLpq/hDofH2mP\nn0Qzo2BL+8N+KnCfrnPsRbOWwwjNFNdPAPt11dkXOA/Y1P4gfQjYe9CffxBt1r6/+72j51s8G9us\nXz9rY5zzTn45MZo17davNmt/Od/Q/v/5TeAJs7XN+tVuNF9SH6fpcbuV5tbHi2dru+2kve4EXtJV\n79S2TTbT3Po5uOv4nPk+6EebMYXfB24iK0mS1HKMkSRJUsvESJIkqWViJEmS1DIxkiRJapkYSZIk\ntUyMJEmSWiZGkiRJLRMjSZKklomRJElSy8RIkiSpZWIkadZI8uUkXxyj/E+TbEzy4EHEJWnmMDGS\nNJu8HDgsyatHC5IcBPwtcGJV3TgZF02y52ScV9LUMzGSNGtU1U+A1wPvTvIbbfE5wBer6gKAJEcm\n+WaSzUl+lOQ9SX5l9BxJXpLk20luTbIuyblJFnYcPyrJXUmOTnJlkq3A4VP4MSVNolTVoGOQpL5K\n8ilgX+BTwFuBR1TVzUkeDlwJvAm4CNgfOBv496r64/a9rwB+AvwnsAhYBdxUVc9vjx8FfBm4CngD\n8CPg5qraNGUfUNKkMTGSNOskeRDwfeABwLFV9U9t+T8Av6iqkzrqPpUm0fmVqto+xrmOAC4F9q6q\nrR2J0bOr6pfGM0ma2byVJmnWqaqfAR8E1owmRa3HAK9qb5PdmuRW4PNAgN8ASHJokn9KckOSW4BL\n2vce2HkJmp4nSbPMfQYdgCRNku3to9Ov0tw6O5smGeo0nOT+wBeBzwEvBm4CDqZJnuZ11b+t3wFL\nGjwTI0lzyRDwyKq6fqyDSZbSjE16U1VtaMueOIXxSRowb6VJmkveCTwlyXuTPDrJwUmen+S97fEb\ngDuAk5MclOT5wJsHFq2kKWdiJGnOqKrvAE8BDgG+STNO6H/SzEKj7SV6BXA8zeDtPwf+YiDBShoI\nZ6VJkiS17DGSJElqmRhJkiS1TIwkSZJaJkaSJEktEyNJkqSWiZEkSVLLxEiSJKllYiRJktQyMZIk\nSWqZGEmSJLVMjCRJklr/D3/61Mz1QLjPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1173de0d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "years = citations_per_year.keys()\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar(years, [citations_per_year[y] for y in years])\n",
    "plt.xlabel(\"Year\")\n",
    "plt.ylabel(\"Citations mentioning genes\")\n",
    "plt.savefig(\"citations_per_year.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finished creating gene_info_by_year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[((9606, 8604, 9722566), (((u'SLC25A12', u'solute carrier family 25 (aspartate/glutamate carrier), member 12', u'protein-coding'), u'NM_003705'), 9722566))]\n",
      "2457278\n"
     ]
    }
   ],
   "source": [
    "joined = (taxid_gene_info_refseq.join(human_gene2pubmed)\n",
    "                                 .filter(lambda x: x[1][1] in pmid_year)\n",
    "                                 .map(lambda x: ((x[0][0], x[0][1], x[1][1]), x[1])))\n",
    "print joined.take(1)\n",
    "print joined.count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((9606, 2475),\n",
       "  (((u'MTOR',\n",
       "     u'mechanistic target of rapamycin (serine/threonine kinase)',\n",
       "     u'protein-coding'),\n",
       "    u'NM_004958'),\n",
       "   20038818)),\n",
       " ((9606, 414332),\n",
       "  (((u'LCN10', u'lipocalin 10', u'protein-coding'), u'NM_001001712'),\n",
       "   14702039)),\n",
       " ((9606, 3356),\n",
       "  (((u'HTR2A',\n",
       "     u'5-hydroxytryptamine (serotonin) receptor 2A, G protein-coupled',\n",
       "     u'protein-coding'),\n",
       "    u'NM_000621'),\n",
       "   24257701)),\n",
       " ((9606, 161835),\n",
       "  (((u'FSIP1', u'fibrous sheath interacting protein 1', u'protein-coding'),\n",
       "    u'NM_152597'),\n",
       "   20379614)),\n",
       " ((9606, 4318),\n",
       "  (((u'MMP9',\n",
       "     u'matrix metallopeptidase 9 (gelatinase B, 92kDa gelatinase, 92kDa type IV collagenase)',\n",
       "     u'protein-coding'),\n",
       "    u'NM_004994'),\n",
       "   19383343))]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rjoined = (joined.reduceByKey(lambda x1, x2: x1)\n",
    "                 .map(lambda x: ((x[0][0], x[0][1]), x[1]))\n",
    "           )\n",
    "rjoined.take(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[((9606, 6597, 1996), (((u'SMARCA4', u'SWI/SNF related, matrix associated, actin dependent regulator of chromatin, subfamily a, member 4', u'protein-coding'), u'NM_001128844'), 5))]\n",
      "time: 14.0963561535\n"
     ]
    }
   ],
   "source": [
    "taxid_gene_info_refseq_pubmed = (rjoined\n",
    "                                 .map(lambda x: ((x[0][0], x[0][1], pmid_year[x[1][1]]), (x[1][0], 1)))\n",
    "                                 .reduceByKey(lambda x1, x2: (x1[0], x1[1] + x2[1]))\n",
    "                                 )\n",
    "\n",
    "t1 = time.time()\n",
    "print taxid_gene_info_refseq_pubmed.take(1)\n",
    "taxid_gene_info_collected1 = taxid_gene_info_refseq_pubmed.collect()\n",
    "t2 = time.time()\n",
    "print \"time:\", t2 - t1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('by_year.tsv', 'w') as f:\n",
    "    for x in taxid_gene_info_collected1:\n",
    "        #print x\n",
    "        refseq_id = x[1][0][1].split('.')[0]\n",
    "        if refseq_id not in refgene_dict:\n",
    "            continue\n",
    "        refseq_data = refgene_dict[refseq_id]\n",
    "        #print refseq_data\n",
    "        out_str = \"\\t\".join(map(str, [x[0][0], x[0][1], x[0][2], refseq_id, refseq_data['chrom'],\n",
    "                                  refseq_data['txStart'], refseq_data['txEnd'],\n",
    "                                  x[1][0][0][0], x[1][0][0][1], x[1][0][0][2], x[1][1]]))\n",
    "        f.write(out_str + \"\\n\")\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taxid_gene_info_all = (taxid_gene_info_refseq_pubmed.map(lambda x: ((x[0][0], x[0][1]), x[1]))\n",
    "                       .reduceByKey(lambda x1,x2: (x1[0], x1[1] + x2[1]))\n",
    "                       )\n",
    "\n",
    "print taxid_gene_info_all.take(1)\n",
    "taxid_gene_info_all_collected = taxid_gene_info_all.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with open('all_gene_counts.tsv', 'w') as f:\n",
    "    for x in taxid_gene_info_all_collected:\n",
    "        #print x\n",
    "        refseq_id = x[1][0][1].split('.')[0]\n",
    "        if refseq_id not in refgene_dict:\n",
    "            continue\n",
    "        refseq_data = refgene_dict[refseq_id]\n",
    "        #print refseq_data\n",
    "        out_str = \"\\t\".join(map(str, [x[0][0], x[0][1], refseq_id, refseq_data['chrom'],\n",
    "                                  refseq_data['txStart'], refseq_data['txEnd'],\n",
    "                                  x[1][0][0][0], x[1][0][0][1], x[1][0][0][2], x[1][1]]))\n",
    "        f.write(out_str + \"\\n\")\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# head all_gene_counts.tsv > short_all_gene_counts.tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import streamgraph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "sorted_non_protein = pd.read_csv('sorted_non_protein_gene_info_by_year.tsv', names=['taxid', 'geneid', 'year', 'symbol', 'name', 'type', 'citations'],\n",
    "                                delimiter='\\t')\n",
    "sorted_non_protein.head()\n",
    "min_year = min(sorted_non_protein['year'])\n",
    "max_year = max(sorted_non_protein['year'])\n",
    "types = dict([(x,i) for i,x in enumerate(set(sorted_non_protein['type']))])\n",
    "\n",
    "print min_year, max_year\n",
    "print \"types:\", types\n",
    "print sorted_non_protein.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import collections as col\n",
    "year_type_counts = col.defaultdict(lambda: col.defaultdict(int))\n",
    "\n",
    "for ix, row in sorted_non_protein.iterrows():\n",
    "    year_type_counts[row['year']][row['type']] += row['citations']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layers = []\n",
    "\n",
    "for type in types:\n",
    "    type_counts = []\n",
    "    for year in range(min_year, max_year+1):\n",
    "    \n",
    "        type_counts += [year_type_counts[year][type]]\n",
    "    layers += [type_counts]\n",
    "\n",
    "print layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as pl\n",
    "import numpy as np\n",
    "\n",
    "def symmetric(sorted_streams, stream_bounds):\n",
    "    \"\"\"Symmetric baseline\"\"\"\n",
    "    lb, ub = np.min(stream_bounds[:,0,:],axis=0), np.max(stream_bounds[:,1,:],axis=0)\n",
    "    return .5*(lb+ub)\n",
    "def pos_only(sorted_streams, stream_bounds):\n",
    "    \"\"\"Lumps will only be positive\"\"\"\n",
    "    lb, ub = np.min(stream_bounds[:,0,:],axis=0), np.max(stream_bounds[:,1,:],axis=0)\n",
    "    return lb\n",
    "def zero(sorted_streams, stream_bounds):\n",
    "    \"\"\"Zero baseline\"\"\"\n",
    "    return np.zeros(stream_bounds.shape[2])\n",
    "def min_weighted_wiggles(sorted_streams, stream_bounds):\n",
    "    \"\"\"Baseline recommended by Byron and Wattenberg\"\"\"\n",
    "    lb, ub = np.min(stream_bounds[:,0,:],axis=0), np.max(stream_bounds[:,1,:],axis=0)\n",
    "    weight = ub-lb\n",
    "    \n",
    "    sorted_streams = np.abs(sorted_streams)\n",
    "    for i in xrange(len(sorted_streams)):\n",
    "        sorted_streams[i,:] *= (-1)**i\n",
    "    cusum_f = np.vstack((np.zeros(sorted_streams.shape[1]),\n",
    "                        np.cumsum(sorted_streams[:-1,:], axis=0)))\n",
    "    f_prime = np.diff(sorted_streams, axis=1)\n",
    "    cusum_f_prime = np.diff(cusum_f, axis=1)\n",
    "    g_prime = np.hstack(([0],-np.sum((f_prime*.5  + cusum_f_prime)*sorted_streams[:,1:],axis=0) / weight[1:]))\n",
    "    g_prime[np.where(weight==0)] = 0\n",
    "    g = np.cumsum(g_prime)\n",
    "    \n",
    "    return g\n",
    "    \n",
    "\n",
    "def stacked_graph(streams, cmap=pl.cm.bone, color_seq='linear', \n",
    "                  baseline_fn=min_weighted_wiggles, ax=None, \n",
    "                  labels=None, xvalues=None):\n",
    "    \"\"\"\n",
    "    Produces stacked graphs using matplotlib.\n",
    "    \n",
    "    Reference: 'Stacked graphs- geometry & aesthetics' by Byron and Wattenberg\n",
    "    http://www.leebyron.com/else/streamgraph/download.php?file=stackedgraphs_byron_wattenberg.pdf\n",
    "    \n",
    "    Parameters:\n",
    "      - streams: A list of time-series of positive values. Each element must be of the same length.\n",
    "      - cmap: A matplotlib color map. Defaults to 'bone'.\n",
    "      - colo_seq: 'linear' or 'random'.\n",
    "      - baseline_fn: Current options are symmetric, pos_only, zero and min_weighted_wiggles.\n",
    "    \"\"\"\n",
    "    # Sort by onset times\n",
    "    onset_times = [np.where(np.abs(stream)>0)[0][0] for stream in streams]\n",
    "    order = np.argsort(onset_times)\n",
    "    streams = np.asarray(streams)\n",
    "    if labels is not None:\n",
    "        labels = [labels[i] for i in order]\n",
    "    sorted_streams = streams[order]\n",
    "    \n",
    "    if xvalues is None:\n",
    "        t = np.arange(streams.shape[1])\n",
    "    else:\n",
    "        t = xvalues\n",
    "    \n",
    "    # Establish bounds\n",
    "    stream_bounds = [np.vstack((np.zeros(streams.shape[1]), sorted_streams[0])),\n",
    "                    np.vstack((-sorted_streams[1], (np.zeros(streams.shape[1]))))]\n",
    "    side = -1\n",
    "    for stream in sorted_streams[2:]:\n",
    "        side *= -1\n",
    "        if side==1:\n",
    "            stream_bounds.append(np.vstack((stream_bounds[-2][1], stream_bounds[-2][1]+stream)))\n",
    "        else:\n",
    "            stream_bounds.append(np.vstack((stream_bounds[-2][0]-stream, stream_bounds[-2][0])))\n",
    "            \n",
    "    stream_bounds = np.array(stream_bounds)\n",
    "    \n",
    "    # Compute baseline\n",
    "    baseline = baseline_fn(sorted_streams, stream_bounds)\n",
    "    \n",
    "    print \"streams.shape\", streams.shape\n",
    "    if labels is None:\n",
    "        labels = ['' for i in range(streams.shape[0])]\n",
    "    \n",
    "    # Choose colors\n",
    "    t_poly = np.hstack((t,t[::-1]))\n",
    "    if color_seq=='linear':\n",
    "        colors = np.linspace(0,1,streams.shape[0])\n",
    "    elif color_seq=='random':\n",
    "        colors = np.random.random(size=streams.shape[1])\n",
    "    else:\n",
    "        raise ValueError, 'Color sequence %s unrecognized'%color_seq\n",
    "    \n",
    "    # Plot    \n",
    "    #pl.axis('off') \n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "    for i in xrange(len(stream_bounds)):\n",
    "        bound = stream_bounds[i]\n",
    "        color = cmap(colors[i])\n",
    "        ax.fill(t_poly, np.hstack((bound[0]-baseline,(bound[1]-baseline)[::-1])), \n",
    "                facecolor=color, linewidth=0.,edgecolor='none',label=labels[i])\n",
    "        \n",
    "        \n",
    "# Demo\n",
    "if __name__ == '__main__':\n",
    "    pl.clf()\n",
    "    N_dsets = 50\n",
    "    T = 100\n",
    "    amp = 1\n",
    "    fade = .15\n",
    "    dsets = []\n",
    "    for i in xrange(N_dsets):\n",
    "        this_dset = np.zeros(T)\n",
    "        t_onset = np.random.randint(.9*T)-T/3\n",
    "        if t_onset >= 0:   \n",
    "            remaining_t = np.arange(T-t_onset)     \n",
    "        else:\n",
    "            remaining_t = np.arange(T)-t_onset\n",
    "        this_dset[max(t_onset,0):]=np.exp(-.15*np.random.gamma(10,.1)*remaining_t)\\\n",
    "                            * remaining_t * np.random.gamma(6,.2) # * np.cos(-fade*remaining_t*np.random.gamma(10,.1))**2\n",
    "        dsets.append(this_dset)\n",
    "    \n",
    "    stacked_graph(dsets, baseline_fn = min_weighted_wiggles, color_seq='linear')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(10,8))\n",
    "stacked_graph(np.array(layers), baseline_fn = symmetric,\n",
    "              cmap=pl.cm.Set2 , color_seq='linear', ax=ax, labels=[t for t in types],\n",
    "             xvalues=range(min_year, max_year+1))\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(handles, labels, loc=3)\n",
    "\n",
    "labels = [item.get_text() for item in ax.get_xticklabels()]\n",
    "print \"labels:\", labels\n",
    "#a =  [x for x in range(min_year, max_year, (max_year - min_year) / len(labels))]\n",
    "print \"len a:\", len(layers[0])\n",
    "ax.get_yaxis().set_visible(False)\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['left'].set_visible(False)\n",
    "\n",
    "print max_year\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Genes by year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1902 2017\n",
      "   taxid  geneid  year    symbol                      name            type  \\\n",
      "0   9606      12  1975  SERPINA3  serpin family A member 3  protein-coding   \n",
      "1   9606      12  1991  SERPINA3  serpin family A member 3  protein-coding   \n",
      "2   9606      12  1987  SERPINA3  serpin family A member 3  protein-coding   \n",
      "3   9606      12  1983  SERPINA3  serpin family A member 3  protein-coding   \n",
      "4   9606      12  1995  SERPINA3  serpin family A member 3  protein-coding   \n",
      "\n",
      "   citations  \n",
      "0          1  \n",
      "1          1  \n",
      "2          2  \n",
      "3          3  \n",
      "4          3  \n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "sorted_gene_info = pd.read_csv('gene_info_by_year.tsv', names=['taxid', 'geneid', 'year', 'symbol', 'name', 'type', 'citations'],\n",
    "                                delimiter='\\t')\n",
    "min_year = min(sorted_gene_info['year']) \n",
    "#min_year = 2000\n",
    "max_year = max(sorted_gene_info['year'])\n",
    "\n",
    "print min_year, max_year\n",
    "print sorted_gene_info.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/pkerp/.virtualenvs/stuff/lib/python2.7/site-packages/ipykernel/__main__.py:10: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)\n"
     ]
    }
   ],
   "source": [
    "import collections as col\n",
    "\n",
    "all_genes = set()\n",
    "all_year_genes = col.defaultdict(lambda: col.defaultdict(int))\n",
    "all_decade_genes = col.defaultdict(lambda: col.defaultdict(int))\n",
    "all_genes_citations = col.defaultdict(int)\n",
    "\n",
    "for year in range(min_year, max_year+1):\n",
    "    for ix, row in  (sorted_gene_info[sorted_gene_info['year'] == year]\n",
    "                    .sort('citations', ascending=[False]).iterrows()):\n",
    "        \n",
    "        all_year_genes[year][row['symbol']] = row['citations']\n",
    "        all_genes.add(row['symbol'])\n",
    "        all_decade_genes[year / 10][row['symbol']] += row['citations']\n",
    "        all_genes_citations[row['symbol']] += row['citations']\n",
    "        \n",
    "        \n",
    "with open('all_year_genes.tsv', 'w') as f:\n",
    "    f.write(\"gene\\t{}\\n\".format(\"\\t\".join(map(str, range(min_year, max_year+1)))))\n",
    "\n",
    "    for gene in all_genes:\n",
    "        str_counts = []    \n",
    "        for year in range(min_year, max_year+1):\n",
    "            if gene in all_year_genes[year]:\n",
    "                str_counts += [str(all_year_genes[year][gene])]\n",
    "            else:\n",
    "                str_counts += ['']\n",
    "            #print year / 10, top_half_decades[year / 10]\n",
    "        out_str = \"{gene}\\t{counts}\".format(gene=gene, counts=\"\\t\".join(str_counts))\n",
    "        f.write(out_str + \"\\n\")\n",
    "        #print out_str"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Genes by decade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('all_decade_genes.tsv', 'w') as f:\n",
    "    f.write(\"gene\\t{}\\n\".format(\"\\t\".join(map(str, range(min_year, max_year+1, 10)))))\n",
    "    for gene in all_genes:\n",
    "        str_counts = []    \n",
    "        for year in range(min_year, max_year+1, 10):\n",
    "            if gene in all_decade_genes[year / 10]:\n",
    "                str_counts += [str(all_decade_genes[year / 10][gene])]\n",
    "            else:\n",
    "                str_counts += ['']\n",
    "            #print year / 10, top_half_decades[year / 10]\n",
    "        out_str = \"{gene}\\t{counts}\".format(gene=gene, counts=\"\\t\".join(str_counts))\n",
    "        f.write(out_str + \"\\n\")\n",
    "        print out_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "layers = []\n",
    "all_genes = set()\n",
    "top_n = 20\n",
    "\n",
    "# tabulate the top ten genes for each year\n",
    "top_tens = col.defaultdict(lambda: col.defaultdict(int))\n",
    "top_half_decades = col.defaultdict(lambda: col.defaultdict(int))\n",
    "        \n",
    "for year in range(min_year, max_year+1):\n",
    "    for ix, row in  (sorted_gene_info[sorted_gene_info['year'] == year]\n",
    "                    .sort('citations', ascending=[False]).head(top_n).iterrows()):\n",
    "        top_tens[year][row['symbol']] = row['citations']\n",
    "        top_half_decades[year / 10][row['symbol']] += row['citations']\n",
    "        all_genes.add(row['symbol'])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stacks = [[0 for i in range(min_year, max_year+1)] for i in range(len(all_genes))]\n",
    "stack_idxs = dict([(x,i) for (i,x) in enumerate(all_genes)])\n",
    "\n",
    "print \"len:\", len(stacks)\n",
    "\n",
    "for i,year in enumerate(range(min_year, max_year+1)):\n",
    "    for gene in top_tens[year]:\n",
    "        #print stack_idxs[gene], stacks[stack_idxs[gene]]\n",
    "        stacks[stack_idxs[gene]][i] = top_tens[year][gene]\n",
    "        all_genes.add(gene)\n",
    "\n",
    "\n",
    "with open('top_ns_year.tsv', 'w') as f:\n",
    "    f.write(\"gene\\t{}\\n\".format(\"\\t\".join(map(str, range(min_year, max_year+1, 10)))))\n",
    "    for gene in all_genes:\n",
    "        str_counts = []    \n",
    "        for year in range(min_year, max_year+1, 10):\n",
    "            if gene in top_half_decades[year / 10]:\n",
    "                str_counts += [str(top_half_decades[year / 10][gene])]\n",
    "            else:\n",
    "                str_counts += ['']\n",
    "            #print year / 10, top_half_decades[year / 10]\n",
    "        out_str = \"{gene}\\t{counts}\".format(gene=gene, counts=\"\\t\".join(str_counts))\n",
    "        f.write(out_str + \"\\n\")\n",
    "        print out_str\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(10,8))\n",
    "stacked_graph(np.array(stacks), baseline_fn = symmetric,\n",
    "              cmap=pl.cm.Set3 , color_seq='linear', ax=ax, labels=None,\n",
    "             xvalues=range(min_year, max_year+1))\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(handles, labels, loc=3)\n",
    "\n",
    "labels = [item.get_text() for item in ax.get_xticklabels()]\n",
    "print \"labels:\", labels\n",
    "#a =  [x for x in range(min_year, max_year, (max_year - min_year) / len(labels))]\n",
    "ax.get_yaxis().set_visible(False)\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['left'].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Count the references\n",
    "\n",
    "For this exercise, the result should be a tsv with the following columns:\n",
    "\n",
    "```\n",
    "<taxId> <geneId> <refSeqID> <citation count>\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_per_year(genbank_dir):\n",
    "    '''\n",
    "    gene2pubmed = (sc.textFile(op.join(genbank_dir, \"gene2pubmed\"))\n",
    "                     .filter(lambda x: x[0] != '#')\n",
    "                     .map(lambda x: x.split('\\t'))\n",
    "                     .map(lambda p: {'taxid': int(p[0]), 'geneid': int(p[1]), 'pmid': int(p[2]), 'count': 1})\n",
    "                     .map(lambda x: ((x['taxid'], x['geneid']), {'count': x['count']}))\n",
    "                     )\n",
    "    '''\n",
    "    pubmed_gene = (sc.textFile(op.join(genbank_dir, \"gene2pubmed\"))\n",
    "                   .filter(lambda x: x[0] != '#')\n",
    "                   .map(lambda x: x.split('\\t'))\n",
    "                   .map(lambda x: (int(x[1]), int(x[0]))\n",
    "                  ))\n",
    "    \n",
    "    pmid_year = (sc.textFile(op.join(genbank_dir, 'pmid_year.ssv'))\n",
    "                 .map(lambda x: x.split())\n",
    "                 .map(lambda x: (int(x[1]), int(x[0].split('-')[0])))\n",
    "                 )\n",
    "    print pmid_year.take(1)\n",
    "    print pubmed_gene.take(1)\n",
    "    \n",
    "    pmid_year_gene = pubmed_gene.join(pmid_year)\n",
    "    return pmid_year_gene\n",
    "    \n",
    "pmid_year_gene = count_per_year(op.join(data_dir, 'genbank-data/'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_gene_counts(genbank_dir):\n",
    "    gene2pubmed = (sc.textFile(op.join(genbank_dir, \"gene2pubmed\"))\n",
    "                     .filter(lambda x: x[0] != '#')\n",
    "                     .map(lambda x: x.split('\\t'))\n",
    "                     .map(lambda p: {'taxid': int(p[0]), 'geneid': int(p[1]), 'pmid': int(p[2]), 'count': 1})\n",
    "                     .map(lambda x: ((x['taxid'], x['geneid']), {'count': x['count']}))\n",
    "                     )\n",
    "    \n",
    "    def reduce_count(r1, r2):\n",
    "        '''\n",
    "        A reduce function that simply counts the number of elements in the table.\n",
    "        \n",
    "        @param r1: A Row\n",
    "        @param r2: A Row\n",
    "        @return: A new Row, equal to the first Row with a summed count.\n",
    "        '''\n",
    "        #print >>sys.stderr, \"r1:\", r1\n",
    "        r1['count'] += r2['count']\n",
    "        return r1\n",
    "\n",
    "    print gene2pubmed.take(1)\n",
    "    reduced_gene2pubmed = gene2pubmed.reduceByKey(reduce_count)\n",
    "    \n",
    "    outfile = get_outfile('taxid-geneid-count')\n",
    "\n",
    "    (reduced_gene2pubmed\n",
    "        .map(lambda x: \"{}\\t{}\\t{}\".format(x[0][0], x[0][1], x[1]['count']))\n",
    "        .saveAsTextFile(outfile)\n",
    "        )\n",
    "    \n",
    "    return reduced_gene2pubmed\n",
    "\n",
    "\n",
    "taxid_geneid_count = load_gene_counts(op.join(data_dir, 'genbank-data/'))\n",
    "print taxid_geneid_count.take(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the gene2refseq annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def take_one(r1, r2):\n",
    "    return r1\n",
    "\n",
    "def load_refseq2gene(genbank_base_dir):\n",
    "    '''\n",
    "    Get the mapping from refseq IDs to gene IDs\n",
    "    \n",
    "    :param genbank_base_dir: The directory that contains all of the genbank files.\n",
    "    :return: A set of tuples of the form (refseq_id, (taxid, geneid))\n",
    "    '''\n",
    "    gene2refseq = (sc.textFile(op.join(genbank_base_dir, 'gene2refseq'))\n",
    "                   .filter(lambda x: x[0] != '#')\n",
    "                   .map(lambda x: x.split('\\t'))\n",
    "                   .map(lambda p: {'taxid': int(p[0]), 'geneid': int(p[1]), 'refseqid': p[3] })\n",
    "                   .map(lambda x: (x['refseqid'].split('.')[0], (x['taxid'], x['geneid'])))\n",
    "                   )\n",
    "    \n",
    "    def reduce_by_refseq_id(r1, r2):\n",
    "        # because we're just looking for a mapping from geneId to refseqId, we just need to throw\n",
    "        # away single entries with identical refseq ids\n",
    "        return r1\n",
    "    \n",
    "    print gene2refseq.take(10)\n",
    "    refseq2gene = gene2refseq.reduceByKey(take_one)\n",
    "    print refseq2gene.take(1)\n",
    "    \n",
    "    outfile = get_outfile('refseqid-taxid-geneid')\n",
    "\n",
    "    (refseq2gene.map(lambda x: \"{}\\t{}\\t{}\".format(x[0], x[1][0], x[1][1]))\n",
    "         .saveAsTextFile(outfile)\n",
    "    )\n",
    "    return refseq2gene\n",
    "\n",
    "refseqid_taxid_geneid = load_refseq2gene(op.join(data_dir, 'genbank-data'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Join the reference counts and refseq to geneid translations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print refseqid_taxid_geneid.take(1)\n",
    "\n",
    "def take_one(r1, r2):\n",
    "    return r1\n",
    "\n",
    "def take_max(r1, r2):\n",
    "    if r1 > r2:\n",
    "        return r1\n",
    "    else:\n",
    "        return r2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def join_counts_and_ids(refseqid_taxid_geneid, taxid_geneid_count):\n",
    "    taxid_geneid_refseq = refseqid_taxid_geneid.map(lambda x: (x[1], x[0]))\n",
    "    print \"count1:\", taxid_geneid_refseq.count()\n",
    "    taxid_geneid_refseq = taxid_geneid_refseq.reduceByKey(take_one)\n",
    "    print \"count2:\", taxid_geneid_refseq.count()\n",
    "    \n",
    "    \n",
    "    '''    \n",
    "    taxid_geneid_refseq = (sc.textFile(op.join(output_dir, 'genbank-output/refseqid-taxid-geneid'))\n",
    "                   .map(lambda x: x.split())\n",
    "                   .map(lambda x: ((int(x[1]), int(x[2])), x[0]))\n",
    "                        )\n",
    "    '''\n",
    "    print taxid_geneid_refseq.take(1)\n",
    "    \n",
    "    '''\n",
    "    (sc.textFile(op.join(output_dir, 'genbank-output/taxid-geneid-count'))\n",
    "                          .map(lambda x: x.split())\n",
    "                          .map(lambda x: ((int(x[0]), int(x[1])), int(x[2])))\n",
    "                          )\n",
    "    '''\n",
    "    print \"taxid_geneid_count\", taxid_geneid_count.take(1)\n",
    "    print \"1. taxid_geneid_count.count():\", taxid_geneid_count.count()\n",
    "    taxid_geneid_count = taxid_geneid_count.reduceByKey(take_max)\n",
    "    print \"2. taxid_geneid_count.count():\", taxid_geneid_count.count()\n",
    "    \n",
    "    taxid_geneid_count_refseq = taxid_geneid_count.join(taxid_geneid_refseq)\n",
    "    print taxid_geneid_count_refseq.take(1)\n",
    "    print taxid_geneid_count.take(1)\n",
    "    \n",
    "    outfile = get_outfile('taxid-geneid-refseqid-count')\n",
    "\n",
    "    (taxid_geneid_count_refseq.map(lambda x: \"{}\\t{}\\t{}\\t{}\".format(x[0][0],\n",
    "                                                                  x[0][1],\n",
    "                                                                  x[1][1],\n",
    "                                                                  x[1][0]))\n",
    "     .saveAsTextFile(outfile)\n",
    "     )\n",
    "    return taxid_geneid_count_refseq\n",
    "\n",
    "taxid_geneid_count_refseq = join_counts_and_ids(refseqid_taxid_geneid, taxid_geneid_count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Join the refgene data with the count data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def join_refgene_and_counts(refGene, taxid_geneid_count_refseq):\n",
    "    '''\n",
    "    Combine the refGene information about the genes with the citation\n",
    "    count information.\n",
    "    '''\n",
    "    refseqid_refgene = refGene.map(lambda x: (x['name'], x))\n",
    "    \n",
    "    print refseqid_refgene.take(1)\n",
    "    \n",
    "    refseqid_count = taxid_geneid_count_refseq.map(lambda x: (x[1][1], x[1][0]))\n",
    "    \n",
    "    print refseqid_count.take(1)\n",
    "    \n",
    "    refseqid_refgene_count = refseqid_refgene.join(refseqid_count)\n",
    "    print refseqid_refgene_count.take(1)\n",
    "    \n",
    "\n",
    "    return refseqid_refgene_count\n",
    "\n",
    "refseqid_refgene_count = join_refgene_and_counts(refGene, taxid_geneid_count_refseq)\n",
    "\n",
    "outfile = get_outfile('refgene-count')\n",
    "(refseqid_refgene_count.map(lambda x: \"{name}\\t{chrom}\\t{strand}\\t{txStart}\\t{txEnd}\\t{genomeTxStart}\\t{genomeTxEnd}\\t{cdsStart}\\t{cdsEnd}\\t{exonCount}\\t{exonStarts}\\t{exonEnds}\\t{geneName}\\t{count}\\t{uid}\"\n",
    "                            .format(count=x[1][1]['count'],uid=shortuuid.uuid(), **x[1][0]))\n",
    " .saveAsTextFile(outfile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print refseqid_refgene_count.take(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "refseqid_refgene_count_plus = refseqid_refgene_count.filter(lambda x: x[1][0]['strand'] == '+')\n",
    "\n",
    "outfile = get_outfile('refgene-count-plus')\n",
    "(refseqid_refgene_count_plus.map(lambda x: \"{name}\\t{chrom}\\t{strand}\\t{txStart}\\t{txEnd}\\t{genomeTxStart}\\t{genomeTxEnd}\\t{cdsStart}\\t{cdsEnd}\\t{exonCount}\\t{exonStarts}\\t{exonEnds}\\t{geneName}\\t{count}\\t{uid}\"\n",
    "                            .format(count=x[1][1]['count'],uid=shortuuid.uuid(), **x[1][0]))\n",
    " .saveAsTextFile(outfile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "refseqid_refgene_count_minus = refseqid_refgene_count.filter(lambda x: x[1][0]['strand'] == '-')\n",
    "\n",
    "outfile = get_outfile('refgene-count-minus')\n",
    "(refseqid_refgene_count_minus.map(lambda x: \"{name}\\t{chrom}\\t{strand}\\t{txStart}\\t{txEnd}\\t{genomeTxStart}\\t{genomeTxEnd}\\t{cdsStart}\\t{cdsEnd}\\t{exonCount}\\t{exonStarts}\\t{exonEnds}\\t{geneName}\\t{count}\\t{uid}\"\n",
    "                            .format(count=x[1][1]['count'],uid=shortuuid.uuid(), **x[1][0]))\n",
    " .saveAsTextFile(outfile))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the entire pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assembly = 'hg19'\n",
    "output_dir = op.join(data_dir, assembly)    # where all of the intermediate output will be stored\n",
    "base_ucsc_dir = op.join(data_dir, 'ucsc-data/{}'.format(assembly))  # where all of the files downloaded from UCSC will be stored\n",
    "\n",
    "cum_chrom_lengths = get_chrom_lengths(base_ucsc_dir)\n",
    "refGene = load_refgene_data(base_ucsc_dir)\n",
    "taxid_geneid_count = load_gene_counts(op.join(data_dir, 'genbank-data/'))\n",
    "refseqid_taxid_geneid = load_refseq2gene(op.join(data_dir, 'genbank-data'))\n",
    "taxid_geneid_count_refseq = join_counts_and_ids(refseqid_taxid_geneid, taxid_geneid_count)\n",
    "refseqid_refgene_count = join_refgene_and_counts(refGene, taxid_geneid_count_refseq)\n",
    "\n",
    "outfile = get_outfile('refgene-count')\n",
    "(refseqid_refgene_count.map(lambda x: \"{name}\\t{chrom}\\t{strand}\\t{txStart}\\t{txEnd}\\t{genomeTxStart}\\t{genomeTxEnd}\\t{cdsStart}\\t{cdsEnd}\\t{exonCount}\\t{exonStarts}\\t{exonEnds}\\t{geneName}\\t{count}\\t{uid}\"\n",
    "                            .format(count=x[1][1]['count'],uid=shortuuid.uuid(), **x[1][0]))\n",
    " .saveAsTextFile(outfile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "refseqid_refgene_count_plus = refseqid_refgene_count.filter(lambda x: x[1][0]['strand'] == '+')\n",
    "\n",
    "outfile = get_outfile('refgene-count-plus')\n",
    "print outfile\n",
    "print refseqid_refgene_count_plus.take(1)\n",
    "(refseqid_refgene_count_plus.map(lambda x: \"{name}\\t{chrom}\\t{strand}\\t{txStart}\\t{txEnd}\\t{genomeTxStart}\\t{genomeTxEnd}\\t{cdsStart}\\t{cdsEnd}\\t{exonCount}\\t{exonStarts}\\t{exonEnds}\\t{geneName}\\t{count}\\t{uid}\"\n",
    "                            .format(count=x[1][1]['count'],uid=shortuuid.uuid(), **x[1][0]))\n",
    " .saveAsTextFile(outfile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/peter/data/hg19/genbank-output/refgene-count-minus\n",
      "[(u'NR_031592', ({'exonEnds': u'10514214', 'geneName': u'MIR1181', 'chromOffset': 2937113968, 'name': u'NR_031592', 'txStart': u'10514133', 'exonCount': u'1', 'strand': u'-', 'cdsEnd': u'10514214', 'genomeTxStart': 2947628101, 'exonStarts': u'10514133', 'cdsStart': u'10514214', 'chrom': u'chr19', 'genomeTxEnd': 2947628182, 'txEnd': u'10514214', 'geneLength': 81}, {'count': 3}))]\n"
     ]
    }
   ],
   "source": [
    "refseqid_refgene_count_plus = refseqid_refgene_count.filter(lambda x: x[1][0]['strand'] == '-')\n",
    "\n",
    "outfile = get_outfile('refgene-count-minus')\n",
    "print outfile\n",
    "print refseqid_refgene_count_plus.take(1)\n",
    "(refseqid_refgene_count_plus.map(lambda x: \"{name}\\t{chrom}\\t{strand}\\t{txStart}\\t{txEnd}\\t{genomeTxStart}\\t{genomeTxEnd}\\t{cdsStart}\\t{cdsEnd}\\t{exonCount}\\t{exonStarts}\\t{exonEnds}\\t{geneName}\\t{count}\\t{uid}\"\n",
    "                            .format(count=x[1][1]['count'],uid=shortuuid.uuid(), **x[1][0]))\n",
    " .saveAsTextFile(outfile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}