pN_pS.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import vcf\n",
      "from collections import defaultdict\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "import pandas as pd\n",
      "import seaborn as sns\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "codon_site_dict =  {'CTT': [3, 0, 3, 0, 0, 3], 'ACC': [3, 0, 3, 0, 0, 3], 'ACA': [3, 0, 3, 0, 0, 3],\n",
      "                    'ACG': [3, 0, 3, 0, 0, 3], 'ATC': [3, 0, 3, 0, 1, 2], 'AAC': [3, 0, 3, 0, 2, 1],\n",
      "                    'ATA': [3, 0, 3, 0, 1, 2], 'AGG': [2, 1, 3, 0, 2, 1], 'CCT': [3, 0, 3, 0, 0, 3],\n",
      "                    'CTC': [3, 0, 3, 0, 0, 3], 'AGC': [3, 0, 3, 0, 2, 1], 'AAG': [2, 0, 3, 0, 2, 1],\n",
      "                    'AGA': [1, 1, 3, 0, 2, 1], 'CAT': [3, 0, 3, 0, 2, 1], 'AAT': [3, 0, 3, 0, 2, 1],\n",
      "                    'ATT': [3, 0, 3, 0, 1, 2], 'CTG': [2, 1, 3, 0, 0, 3], 'CTA': [2, 1, 3, 0, 0, 3],\n",
      "                    'ACT': [3, 0, 3, 0, 0, 3], 'CAC': [3, 0, 3, 0, 2, 1], 'AAA': [2, 0, 3, 0, 2, 1],\n",
      "                    'CCG': [3, 0, 3, 0, 0, 3], 'AGT': [3, 0, 3, 0, 2, 1], 'CCA': [3, 0, 3, 0, 0, 3],\n",
      "                    'CAA': [2, 0, 3, 0, 2, 1], 'CCC': [3, 0, 3, 0, 0, 3], 'TAT': [3, 0, 3, 0, 0, 1],\n",
      "                    'GGT': [3, 0, 3, 0, 0, 3], 'TGT': [3, 0, 3, 0, 1, 1], 'CGA': [1, 1, 3, 0, 0, 3],\n",
      "                    'CAG': [2, 0, 3, 0, 2, 1], 'TCT': [3, 0, 3, 0, 0, 3], 'GAT': [3, 0, 3, 0, 2, 1],\n",
      "                    'CGG': [2, 1, 3, 0, 0, 3], 'TTT': [3, 0, 3, 0, 2, 1], 'TGC': [3, 0, 3, 0, 1, 1],\n",
      "                    'GGG': [3, 0, 3, 0, 0, 3], 'GGA': [2, 0, 3, 0, 0, 3], 'TGG': [3, 0, 2, 0, 2, 0],\n",
      "                    'GGC': [3, 0, 3, 0, 0, 3], 'TAC': [3, 0, 3, 0, 0, 1], 'TTC': [3, 0, 3, 0, 2, 1],\n",
      "                    'TCG': [3, 0, 2, 0, 0, 3], 'TTA': [2, 1, 1, 0, 2, 1], 'TTG': [2, 1, 2, 0, 2, 1],\n",
      "                    'TCC': [3, 0, 3, 0, 0, 3], 'GAA': [2, 0, 3, 0, 2, 1], 'TCA': [3, 0, 1, 0, 0, 3],\n",
      "                    'GCA': [3, 0, 3, 0, 0, 3], 'GTA': [3, 0, 3, 0, 0, 3], 'GCC': [3, 0, 3, 0, 0, 3],\n",
      "                    'GTC': [3, 0, 3, 0, 0, 3], 'GCG': [3, 0, 3, 0, 0, 3], 'GTG': [3, 0, 3, 0, 0, 3],\n",
      "                    'GAG': [2, 0, 3, 0, 2, 1], 'GTT': [3, 0, 3, 0, 0, 3], 'GCT': [3, 0, 3, 0, 0, 3],\n",
      "                    'GAC': [3, 0, 3, 0, 2, 1], 'CGT': [3, 0, 3, 0, 0, 3], 'ATG': [3, 0, 3, 0, 3, 0],\n",
      "                    'CGC': [3, 0, 3, 0, 0, 3]}\n",
      "\n",
      "comp_nuc_dict = {'A':'T',\"T\":\"A\",\"G\":\"C\",\"C\":\"G\"}\n",
      "codon_dict = {'TTT':'F', 'TTC':'F', 'TTA':'L', 'TTG':'L',\n",
      "              'TCT':'S', 'TCC':'S', 'TCA':'S', 'TCG':'S',\n",
      "              'TAT':'Y', 'TAC':'Y', 'TAA':'*', 'TAG':'*',\n",
      "              'TGT':'C', 'TGC':'C', 'TGA':'*', 'TGG':'W',\n",
      "              'CTT':'L', 'CTC':'L', 'CTA':'L', 'CTG':'L',\n",
      "              'CCT':'P', 'CCC':'P', 'CCA':'P', 'CCG':'P',\n",
      "              'CAT':'H', 'CAC':'H', 'CAA':'Q', 'CAG':'Q',\n",
      "              'CGT':'R', 'CGC':'R', 'CGA':'R', 'CGG':'R',\n",
      "              'ATT':'I', 'ATC':'I', 'ATA':'I', 'ATG':'M',\n",
      "              'ACT':'T', 'ACC':'T', 'ACA':'T', 'ACG':'T',\n",
      "              'AAT':'N', 'AAC':'N', 'AAA':'K', 'AAG':'K',\n",
      "              'AGT':'S', 'AGC':'S', 'AGA':'R', 'AGG':'R',\n",
      "              'GTT':'V', 'GTC':'V', 'GTA':'V', 'GTG':'V',\n",
      "              'GCT':'A', 'GCC':'A', 'GCA':'A', 'GCG':'A',\n",
      "              'GAT':'D', 'GAC':'D', 'GAA':'E', 'GAG':'E',\n",
      "              'GGT':'G', 'GGC':'G', 'GGA':'G', 'GGG':'G',\n",
      "              '---':'-'}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#All the genome scaffold information\n",
      "from collections import defaultdict\n",
      "genome_scaf_info = defaultdict(list)\n",
      "\n",
      "genomes_order = [\"J07HWQ1\", \"J07HWQ2\", \"J07HQX50\", \"J07AB56\", \"J07AB43\", \"J07HN4\", \"J07HN6\", \"J07HX64\",\n",
      "                 \"J07HX5\", \"J07HB67\", \"J07HR59\", \"A07HB70\", \"A07HR67\", \"A07HN63\", \"A07HR60\", \"J07SB\", \"G22\"]\n",
      "\n",
      "sample_ids = [\"Jan23\", \"Jan25\", \"Aug07\", \"Aug09\"]\n",
      "\n",
      "for line in open(\"genome_scaffold_list.txt\"):\n",
      "    if line.strip():\n",
      "        line=line.rstrip()\n",
      "        genome,scaffold = line.split(\"\\t\")\n",
      "        genome_scaf_info[genome].append(scaffold)\n",
      "\n",
      "scaffold_length = defaultdict(int)\n",
      "\n",
      "for line in open(\"scaffold_length.txt\"):\n",
      "    if line.strip():\n",
      "        line=line.rstrip()\n",
      "        scaffold,length = line.split(\"\\t\")\n",
      "        scaffold_length[scaffold] = length\n",
      "        \n",
      "genome_length = defaultdict(int)\n",
      "\n",
      "for line in open(\"genome_lengths.txt\"):\n",
      "    if line.strip():\n",
      "        line=line.rstrip()\n",
      "        genome,length=line.split(\"\\t\")\n",
      "        genome_length[genome]=length"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Genome name to JGI ID number\n",
      "name_jgi_id = {\"G22\": \"2512047006\",\n",
      "\"J07SB\": \"2513020022\",\n",
      "\"A07HR60\": \"2518645513\",\n",
      "\"A07HN63\": \"2528768014\",\n",
      "\"A07HR67\": \"2518645514\",\n",
      "\"A07HB70\": \"2518645512\",\n",
      "\"J07HR59\": \"2512875011\",\n",
      "\"J07HB67\": \"2506783034\",\n",
      "\"J07HX5\": \"2512875010\",\n",
      "\"J07HX64\": \"2502082092\",\n",
      "\"J07HN6\": \"2512875008\",\n",
      "\"J07HN4\": \"2512875007\",\n",
      "\"J07AB43\": \"2502422326\",\n",
      "\"J07AB56\": \"2502422327\",\n",
      "\"J07HQX50\": \"2512875009\",\n",
      "\"J07HWQ2\": \"2512875006\",\n",
      "\"J07HWQ1\": \"2512875005\"}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jgi_annotation_folder = \"/home/ubuntu/LT_data/New_Reads/jgi_annotation/\"\n",
      "vcf_folder = \"/home/ubuntu/LT_data/New_Reads/Mapping/variation/separated_results/by_month/\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def calculate_dn_ds(input_gff,input_vcf,input_fasta):\n",
      "    import re\n",
      "    from Bio import SeqIO\n",
      "    from Bio.Seq import Seq\n",
      "    from cogent import LoadSeqs,DNA\n",
      "    \n",
      "    #Read the GFF\n",
      "    gene_info = defaultdict(lambda:defaultdict())\n",
      "    for line in open(input_gff):\n",
      "        if line.strip():\n",
      "            if line.startswith(\"#\"):\n",
      "                continue\n",
      "            \n",
      "            line = line.rstrip()\n",
      "            elements = line.split(\"\\t\")\n",
      "            #save type, start, stop, orientation\n",
      "            gene_info[elements[0]][elements[8].split(\";\")[0].split(\"=\")[1]]=[elements[2],elements[3],elements[4],elements[6]]\n",
      "    \n",
      "    #Read the VCF\n",
      "    snp_gene_info = defaultdict(lambda:defaultdict())\n",
      "    \n",
      "    vcf_reader = vcf.Reader(open(input_vcf, 'r'))\n",
      "    \n",
      "    for record in vcf_reader:\n",
      "        if record.INFO[\"TYPE\"][0] == \"snp\":\n",
      "            snp_gene_info[record.CHROM][record.POS] = (str(record.REF), str(record.ALT[0]))\n",
      "    \n",
      "    #Read the fasta file\n",
      "    record_dict = SeqIO.index(input_fasta, \"fasta\")\n",
      "    \n",
      "    #Generate snps by gene\n",
      "    snps_by_gene = defaultdict(lambda:defaultdict(list))\n",
      "    \n",
      "    for scaffold in gene_info:\n",
      "        for gene in gene_info[scaffold]:\n",
      "            \n",
      "            start = int(gene_info[scaffold][gene][1])\n",
      "            stop = int(gene_info[scaffold][gene][2])\n",
      "            \n",
      "            for position in range(start,stop+1):\n",
      "                if position in snp_gene_info[scaffold].keys():\n",
      "                    \n",
      "                    #Only store coding sequences\n",
      "                    if not gene_info[scaffold][gene][0] == \"CDS\":\n",
      "                        continue\n",
      "\n",
      "                    position_in_gene = int(position) - start + 1\n",
      "                    snps_by_gene[gene][position_in_gene] = snp_gene_info[scaffold][position], position\n",
      "            \n",
      "    #########\n",
      "    #Calculate pN/pS\n",
      "    mixture_dict={'N':'ATGC'}\n",
      "    mixture_re = re.compile(\"[N]\")\n",
      "    stop_codons=[\"TAA\",\"TGA\",\"TAG\"]\n",
      "    gene_dn_ds_list=[]\n",
      "    \n",
      "    #Get the fasta sequence per gene\n",
      "    fasta_sequence_gene = defaultdict(str)\n",
      "    for scaffold in gene_info:\n",
      "        scaf_sequence = record_dict[scaffold].seq\n",
      "        \n",
      "        for gene in gene_info[scaffold]:\n",
      "            feat_type, start, stop, strand = gene_info[scaffold][gene]\n",
      "            \n",
      "            if feat_type == \"CDS\":\n",
      "                \n",
      "                dna_sequence = str(scaf_sequence[int(start)-1:int(stop)])\n",
      "                \n",
      "                if strand == \"-\":\n",
      "                    fasta_sequence_gene[gene] = str(Seq(dna_sequence).reverse_complement()) # Use Biopython to reverse complement\n",
      "                else:\n",
      "                    fasta_sequence_gene[gene] = str(scaf_sequence[int(start)-1:int(stop)])\n",
      "                \n",
      "    \n",
      "    #Calculate pn_ps\n",
      "    gene_pn_ps = defaultdict(list)\n",
      "    \n",
      "    for gene in snps_by_gene:\n",
      "        if len(snps_by_gene[gene]) == 0:\n",
      "            continue\n",
      "        \n",
      "        #Get sequence and codon list\n",
      "        seq_object_pycogent = DNA.makeSequence(fasta_sequence_gene[gene], gene)\n",
      "        codon_list = seq_object_pycogent.getInMotifSize(3)\n",
      "        \n",
      "        #Calculate expected NS and S values\n",
      "        gene_expected_NS = 0\n",
      "        gene_expected_S = 0\n",
      "        gene_observed_NS = 0\n",
      "        gene_observed_S = 0\n",
      "        \n",
      "        for codon in codon_list:\n",
      "            #Check for stop codons\n",
      "            if codon in stop_codons:\n",
      "                continue\n",
      "            \n",
      "            #Deal with Ns\n",
      "            mixtures = mixture_re.findall(codon)\n",
      "            \n",
      "            if len(mixtures) > 1:\n",
      "                #codon contains more than one mixture, no way to resolve\n",
      "                continue\n",
      "            \n",
      "            elif len(mixtures)==1:\n",
      "                resolutions = mixture_dict[mixtures[0]]\n",
      "                means = [0,0]\n",
      "                \n",
      "                for r in resolutions:\n",
      "                    rcodon = codon.replace(mixtures[0],r)\n",
      "                    \n",
      "                    if rcodon in stop_codons:\n",
      "                        break\n",
      "                    else:\n",
      "                        position = 0\n",
      "                        \n",
      "                        while position < 3:\n",
      "                            index = position*2\n",
      "                            means[0] += float(codon_site_dict[rcodon][index] / float(len(resolutions)))\n",
      "                            means[1] += float(codon_site_dict[rcodon][index+1]) / float(len(resolutions))\n",
      "                            position += 1\n",
      "                \n",
      "                \n",
      "                gene_expected_NS =+ means[0]\n",
      "                gene_expected_S =+ means[1]\n",
      "                continue                \n",
      "        \n",
      "            #Iterate over each position of the codon\n",
      "            position = 0\n",
      "            while position < 3:\n",
      "                index = position*2\n",
      "                expected_NS,expected_S = tuple(codon_site_dict[codon][index:index+2])\n",
      "                position += 1\n",
      "                gene_expected_NS += expected_NS\n",
      "                gene_expected_S += expected_S  \n",
      "\n",
      "        #Calculate the observed frequencies  \n",
      "    \n",
      "        codon_errors = 0\n",
      "        errors = []\n",
      "        \n",
      "        for snp in snps_by_gene[gene]:\n",
      "            codon_division = 0\n",
      "            codon_number = 0\n",
      "            snp_in_codon_position = 0\n",
      "            codon_in_list = 0\n",
      "            take_complementary_base = 0\n",
      "            codon_division = snp%3\n",
      "            codon_number = snp/3\n",
      "            \n",
      "            if codon_division==1:\n",
      "                snp_in_codon_position=1\n",
      "                codon_number==codon_number+1\n",
      "            if codon_division==2:\n",
      "                snp_in_codon_position=2\n",
      "                codon_number=codon_number+1\n",
      "            if codon_division==0:\n",
      "                snp_in_codon_position=3\n",
      "                \n",
      "            #Get the codon information\n",
      "            try:\n",
      "                consensus_codon=codon_list[codon_number-1]\n",
      "            except IndexError:\n",
      "                codon_errors += 1\n",
      "    \n",
      "            bases_codon=list(consensus_codon)\n",
      "                                    \n",
      "            if consensus_codon in codon_dict:\n",
      "                consensus_aminoacid=codon_dict[consensus_codon]\n",
      "            else:\n",
      "                #Skips N and failed results\n",
      "                errors.append([gene,consensus_codon])\n",
      "                continue\n",
      "                        \n",
      "            new_codon = []\n",
      "            try:\n",
      "                ref_base,new_base = snps_by_gene[gene][snp][0]\n",
      "            except ValueError:\n",
      "                continue\n",
      "                            \n",
      "            \n",
      "            if snp_in_codon_position==1:\n",
      "                new_codon=str(new_base),bases_codon[1],bases_codon[2]\n",
      "            elif snp_in_codon_position==2:\n",
      "                new_codon=bases_codon[0],str(new_base),bases_codon[2]\n",
      "            elif snp_in_codon_position==3:\n",
      "                new_codon=bases_codon[0],bases_codon[1],str(new_base)\n",
      "                              \n",
      "            try:\n",
      "                new_aminoacid = codon_dict[''.join(new_codon)]\n",
      "            except KeyError:\n",
      "                continue\n",
      "            \n",
      "            #Add into observed NS or S\n",
      "            \n",
      "            if new_aminoacid==consensus_aminoacid:\n",
      "                gene_observed_S += 1\n",
      "            else:\n",
      "                gene_observed_NS += 1\n",
      "    \n",
      "        #Calculate dn/ds ratio\n",
      "    \n",
      "        if gene_expected_NS > 0:\n",
      "            dn=float(gene_observed_NS) / gene_expected_NS\n",
      "        else:\n",
      "            dn = 0\n",
      "            \n",
      "        if gene_expected_S > 0:\n",
      "            ds=float(gene_observed_S) / gene_expected_S\n",
      "        else:\n",
      "            ds = 0\n",
      "        \n",
      "        if ds>0:\n",
      "            dn_ds=round(dn/ds,3)\n",
      "        else:\n",
      "            if dn == 0:\n",
      "                dn_ds = \"NaN\"\n",
      "            else:\n",
      "                dn_ds = \"Inf\"\n",
      "        \n",
      "        gene_pn_ps[gene] = [dn,ds,dn_ds]\n",
      "        \n",
      "    return gene_pn_ps\n",
      "    \n",
      "        \n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# import cPickle as pickle\n",
      "# pN_pS_Jan = defaultdict()\n",
      "# pN_pS_Aug = defaultdict()\n",
      "\n",
      "# #test_genome = [\"A07HR60\"]\n",
      "\n",
      "# for genome in genomes_order:\n",
      "#     print \"Analyzing genome: \" + genome\n",
      "    \n",
      "#     gff_file = \"/home/ubuntu/LT_data/New_Reads/jgi_annotation/\" + genome + \"/\" + name_jgi_id[genome] + \".gff\"\n",
      "#     vcf_file_jan = \"/home/ubuntu/LT_data/New_Reads/Mapping/variation/separated_results/by_month/\" + genome + \"_Jan.vcf\"\n",
      "#     vcf_file_aug = \"/home/ubuntu/LT_data/New_Reads/Mapping/variation/separated_results/by_month/\" + genome + \"_Aug.vcf\"\n",
      "#     fasta_file = \"/home/ubuntu/LT_data/New_Reads/jgi_annotation/\" + genome + \"/\" + name_jgi_id[genome] + \".fna\"\n",
      "    \n",
      "#     print \"January\"\n",
      "#     pN_pS_Jan[genome] = calculate_dn_ds(gff_file,vcf_file_jan,fasta_file)\n",
      "#     print \"August\"\n",
      "#     pN_pS_Aug[genome] = calculate_dn_ds(gff_file,vcf_file_aug,fasta_file)\n",
      "\n",
      "    \n",
      "# pickle.dump(pN_pS_Jan, open(\"pN_pS_Jan.dict\", 'wb'))\n",
      "# pickle.dump(pN_pS_Aug, open(\"pN_pS_Aug.dict\", 'wb'))\n",
      "    \n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cPickle as pickle\n",
      "pN_pS_Jan = pickle.load(open(\"pN_pS_Jan.dict\", 'r'))\n",
      "pN_pS_Aug = pickle.load(open(\"pN_pS_Aug.dict\", 'r'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Print a list of the genomes, adding the product\n",
      "annotation_files = \"/home/ubuntu/LT_data/CompGenomics/LT_Genomes/v2/processed_genomes/annotation/\"\n",
      "\n",
      "#Count the number of CDS, based on the annotation\n",
      "cds_per_genome = defaultdict(int)\n",
      "for genome in genomes_order:\n",
      "    for line in open(annotation_files + genome + \".txt\", 'r'):\n",
      "        if line.strip():\n",
      "            line = line.rstrip()\n",
      "            elements = line.split(\"\\t\")\n",
      "            \n",
      "            if elements[2] == \"CDS\":\n",
      "                cds_per_genome[genome] += 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print cds_per_genome"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "defaultdict(<type 'int'>, {'J07SB': 1641, 'J07HWQ1': 3584, 'G22': 3525, 'A07HB70': 2514, 'J07HB67': 2847, 'J07AB43': 1678, 'J07HN6': 2914, 'J07HN4': 3230, 'J07AB56': 1411, 'J07HX64': 3049, 'A07HR67': 2891, 'J07HR59': 1841, 'A07HR60': 2861, 'J07HWQ2': 3856, 'J07HQX50': 2872, 'J07HX5': 2139, 'A07HN63': 2507})\n"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for genome in genomes_order:\n",
      "    count = 0\n",
      "    for gene in pN_pS_Jan[genome]:\n",
      "        if float(pN_pS_Jan[genome][gene][2]) > 1 and not np.isinf(float(pN_pS_Jan[genome][gene][2])):\n",
      "                count += 1\n",
      "    \n",
      "    print genome, count, cds_per_genome[genome]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "J07HWQ1 61 3584\n",
        "J07HWQ2 46 3856\n",
        "J07HQX50 19 2872\n",
        "J07AB56 96 1411\n",
        "J07AB43 96 1678\n",
        "J07HN4 198 3230\n",
        "J07HN6 59 2914\n",
        "J07HX64 191 3049\n",
        "J07HX5 137 2139\n",
        "J07HB67 196 2847\n",
        "J07HR59 27 1841\n",
        "A07HB70 154 2514\n",
        "A07HR67 167 2891\n",
        "A07HN63 42 2507\n",
        "A07HR60 83 2861\n",
        "J07SB 91 1641\n",
        "G22 0 3525\n"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib import rcParams\n",
      "rcParams['axes.labelsize'] = 18 #18\n",
      "rcParams['xtick.labelsize'] = 8 # 15\n",
      "rcParams['ytick.labelsize'] = 8 #15\n",
      "rcParams['legend.fontsize'] = 14 # 18\n",
      "rcParams['font.family'] = 'serif'\n",
      "rcParams['font.serif'] = ['Computer Modern Roman']\n",
      "rcParams['text.usetex'] = False\n",
      "rcParams['figure.figsize'] = 5.8,3.5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot_dict = defaultdict(list)\n",
      "jan_genes = defaultdict(list)\n",
      "for genome in genomes_order:\n",
      "    count = 0\n",
      "    \n",
      "    for gene in pN_pS_Jan[genome]:\n",
      "        if float(pN_pS_Jan[genome][gene][2]) > 1 and not np.isinf(float(pN_pS_Jan[genome][gene][2])):\n",
      "            count += 1\n",
      "            jan_genes[genome].append(gene)\n",
      "            \n",
      "    plot_dict[\"CDS\"].append(cds_per_genome[genome])\n",
      "    plot_dict[\"PS\"].append(count)\n",
      "    plot_dict[\"pct\"].append(float(count)/cds_per_genome[genome]*100)\n",
      "    \n",
      "df = pd.DataFrame(plot_dict, index=genomes_order)\n",
      "df.plot(kind=\"bar\",stacked=True,legend=True)\n",
      "    \n",
      "print df"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "           CDS   PS       pct\n",
        "J07HWQ1   3584   61  1.702009\n",
        "J07HWQ2   3856   46  1.192946\n",
        "J07HQX50  2872   19  0.661560\n",
        "J07AB56   1411   96  6.803685\n",
        "J07AB43   1678   96  5.721097\n",
        "J07HN4    3230  198  6.130031\n",
        "J07HN6    2914   59  2.024708\n",
        "J07HX64   3049  191  6.264349\n",
        "J07HX5    2139  137  6.404862\n",
        "J07HB67   2847  196  6.884440\n",
        "J07HR59   1841   27  1.466594\n",
        "A07HB70   2514  154  6.125696\n",
        "A07HR67   2891  167  5.776548\n",
        "A07HN63   2507   42  1.675309\n",
        "A07HR60   2861   83  2.901084\n",
        "J07SB     1641   91  5.545399\n",
        "G22       3525    0  0.000000\n",
        "\n",
        "[17 rows x 3 columns]\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/home/ubuntu/anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:1236: UserWarning: findfont: Font family ['serif'] not found. Falling back to Bitstream Vera Sans\n",
        "  (prop.get_family(), self.defaultFamily[fontext]))\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAEPCAYAAABm//5NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX++PHXDMjNYbiJmIq3BE0N8ZK7oGVecrNQtLQI\n0/K6Wm1bSaWZ/rL1mrabbStLmVkmyWbmmpoa2qaG18xrqaAJaYoX0GG4KDDn9wdfJlCYCxxgBt7P\nx6PHAw5n3vM5x9N7znzO+/P5aBRFURBCCOHwtHXdACGEELaRhC2EEE5CErYQQjgJSdhCCOEkJGEL\nIYSTkIQthBBOwqaEvWLFCmJjYwFYtmwZsbGxxMXFUVRUBMD69euJiYlh8uTJGI1GAHbv3k1MTAxj\nxowhMzOzhpovhBANh9WEffPmTU6cOIFGoyErK4t9+/aRmJhIhw4dSE5OprCwkKSkJBITExk6dChJ\nSUkAxMfHs3z5cuLi4khISKjxAxFCiPrOasL+/PPPGTZsGIqicPToUXr16gVAZGQkhw4dIiMjg9DQ\nULRarXlbQUEBHh4eeHl5ERYWRmpqao0fiBBC1HcWE3ZhYSH79+/nj3/8IwA5OTnodDoAdDodBoMB\ng8FgcRuAyWSqqfYLIUSDYTFh//e//yUqKsr8u7e3t7mP2mg0otfrrW4D0Gqtd5UXFRVX6QCEEKKh\ncLX0x7Nnz/Lzzz+zevVq0tLSOHr0KEePHmXChAmkpKQQHh5OmzZtSE1NxWQymbd5enpSUFBAXl4e\naWlphISEWG1Idnae1X0CA725fDnH9qNzsliO2CZHjeWIbXLUWI7YpoYQq6pxAgO9K/2bxYQdFxdn\n/nnUqFE899xzfPDBB8TGxtK8eXPGjh2Lq6srI0eOJDY2Fl9fXxYvXgzAlClTGDduHO7u7ixcuNDu\nRgshhCjPYsIua9WqVQBMnDiRiRMnlvtbdHQ00dHR5bZFREQQERGhQhOFEEKAHQnbWRUXF7Njx7e3\nbffx8eL69fLdMPfd1w8XF5faapoQQtil3ifss2fP8Na2JXj5N7a4X15WLq1atebOO633twshRF2o\n9wm7uNiEV4COxhY68gHQaCgulvJDIYTjqvcJGxSuHWjLDW9/i3vl52TBw7L4jhDCcdX7hO3i4kJA\ny7vQ+bWwuJ8x+7z0XwshHJrM1ieEEE5CErYQQjgJSdhCCOEkJGELIYSTkIQthBBOQhK2EEI4CUnY\nQgjhJByyDru4uJizZ8/ctj07W0dWlrHctjZt2kn9tBCiQXDIhH36dBrrnplIoIenxf0uF+QzbOkH\nhIZ2qKWWCSFE3XHIhA0K3/7Rx4YJm9wYhgwnF0I0DA6ZsF1cXAjseAfezX0t7pfz2zXpDhFCNBgO\nmbCLi03k2rC0Tu7lHJlhTwjRYFhM2KmpqcycORMXFxfat2/P7Nmz6dGjB506dUKj0fDee++h1+tZ\nv349iYmJ5iXCdDodu3fvZsmSJbi5ubFo0SKCgoLsaJbMsCeEELeymLDbtm3L6tWrAZg+fTo//fQT\nHTp0YOXKleZ9CgsLSUpKIjExkc2bN5OUlMT48eOJj49n+fLlpKWlkZCQwKxZs2xulMywJ4RwJhVV\nttVEVZvFhO3q+vufb9y4gbe3N6dPn2bUqFF0796dqVOnkp6eTmhoKFqtlsjISGbOnElBQQEeHh54\neXkRFhbGokWLqtxAIYRwdGfPnuGvi9bj5dO00n3yrl9iyctDq7WqldU+7G3btvHOO+/QpUsXgoOD\n+eabb9Dr9cyaNYvt27fj5+eHTqcDQKfTYTAYMBgM5m0AJpP0Mwsh6jcvn6ZWewWqy2rCHjBgAAMG\nDGDOnDl8//339O7dG4CBAwfy888/M2DAAIzGktt+o9GIXq/H29vbvA1Aq7U+oNLPzwtX15KvCtnZ\nOit7/87fX0egheW/1IxlSVVfV1NxGkIsR2yTo8ZyxDbVp1i25pnq5BiwkrBv3ryJm5sbAI0bN6ag\noACTyYRWq+WHH36gY8eOtGnThtTUVEwmEykpKYSHh+Pp6UlBQQF5eXmkpaUREmL9K0B29u8rmN/a\n72NJVpaRyxYqStSMVZnAQO8qva6m4jSEWI7YJkeN5Yhtqm+xbM0ztuQYSwndYsLeuXMnK1asQFEU\nWrZsSVBQECNGjMDLy4vg4GBeeOEFNBoNI0eOJDY21lwlAjBlyhTGjRuHu7s7CxcutOlgRN0qLi5m\nx45vb9vu4+PF9et55bbdd18/eeArRC2zmLBLu0PKWrt27W37RUdHEx0dXW5bREQEERERKjRR1Jaz\nZ8/w1rYlNowwzaVVq9bVengihLCfQw6cEXXH1hGmQojaJ9OrCiGEk5CELYQQTkISthBCOAlJ2EII\n4SQkYQshhJOQhC2EEE5CyvqEw6toQI8M5hENkSRsYeaoC0fYMqBHBvOIhkAStijDcReOsDagRwbz\niIZAErYwk4UjhHBs8tBRCCGchCRsIYRwEpKwhRDCSUjCFkIIJyEJWwghnITFKpHU1FRmzpyJi4sL\n7du3Z/bs2Sxbtozt27fTvHlzFixYgKurK+vXrycxMdG84oxOp2P37t0sWbIENzc3Fi1aRFBQUG0d\nkxBCBTJgyfFYTNht27Zl9erVAEyfPp0jR46wb98+EhMT+eCDD0hOTmbAgAEkJSWRmJjI5s2bSUpK\nYvz48cTHx7N8+XLS0tJISEhg1qxZtXJAov6xZUBPbQ/maQjLqcmAJcdjMWG7uv7+5xs3bnD06FF6\n9eoFQGRkJF999RUhISGEhoai1WqJjIxk5syZFBQU4OHhgZeXF2FhYSxatKhmj0LUc9YH9NT2YJ6G\nspyaDFhyLFYHzmzbto133nmHzp074+PjY75T0Ol0GAwGDAYDOp2u0m0AJlPt3fmI+seWAT11MZhH\nllMTtc1qwi5diHfOnDl4enpy8eJFAIxGI3q9Hm9vb4xGY6XbALRa6882/fy8cHUt+R8uO1tnZe/f\n+fvrLC4Lr2YsS6r6upqKU5VYap6r4uJikpOTrcYZOHCg1URra7tq89+vIVxXjnjeHTVWbZwrsJKw\nb968iZubGwCNGzcmJyeH/fv3M2HCBFJSUggPD6dNmzakpqZiMpnM2zw9PSkoKCAvL4+0tDRCQqx/\nHczO/r3fLyvLaGHP8rKyjFy20L+pZqzKBAZ6V+l1NRWnqrHUPFenT6cyY+18q/2fPj6BVrsLbG1X\nbf77NYTryhHPu6PGUvNcWUroFhP2zp07WbFiBYqi0LJlS55//nkuX75MbGwszZs3Z+zYsbi6ujJy\n5EhiY2PNVSIAU6ZMYdy4cbi7u7Nw4UKbDkbUL9L/KdRm68NeZ33Qa43FhF3aHVLWxIkTmThxYrlt\n0dHRREdHl9sWERFBRESESs0UzsYRKzuE82volSsyW58dGvqnu30cr7JD2MdRP3Qb8jc3Sdh2OH06\njdlr5uGh96p0nwJDHu+3CCY0tEMttszxOGplh7CHfOg6GknYdlFwufZHXIsrv4BdcrIAuYCF83PE\nD11HveuvLZKw7eCIF7CoG466nFr917Dv+iVhC1EljrucWn3W0G+aJGELUQWynJqoCzK9qhBCOAlJ\n2EII4SQkYQshhJOQhC2EEE5CErYQQjgJqRIRoh6RZb3qN0nYQtQjDX1ypPpOErYQ9UxDnhypvpM+\nbCGEcBKSsIUQwklYTNiHDx8mJiaG2NhY5s+fD0CPHj0YPXo0Y8aMwWAwALB+/XpiYmKYPHmyeS3H\n3bt3ExMTw5gxY8jMzKzhwxBCiPrPYh92ixYt+OSTT3BzcyMuLo5Tp07RoUMHVq5cad6nsLCQpKQk\nEhMT2bx5M0lJSYwfP574+HiWL19OWloaCQkJzJo1q8YPRgjheGxd+AOkesUaiwm7SZMm5p8bNWqE\ni4sLp0+fZtSoUXTv3p2pU6eSnp5OaGgoWq2WyMhIZs6cSUFBAR4eHnh5eREWFsaiRYtq/ECEEI7J\nlsoVkOoVW9hUJXLixAmysrK48847+eabb9Dr9cyaNYvt27fj5+eHTleyxLtOp8NgMGAwGMzbAEwm\nmQ9YiNrgiBP8Fxeb8ArQ0djCauAAaDQyd7gVVhP2tWvXmDNnDkuWLAFAr9cDMHDgQH7++WcGDBhg\n7rc2Go3o9Xq8vb3N2wC0WuvPNv38vHB1LfkqlJ2ts7L37/z9dRaXha+LWNbiWFLV16kRy1HPlSOe\ndzXPlSX2vu7KlcY2TfDvP7pxrZ13W9pka7tq41qAmrseqtsuiwm7qKiIl19+mVdeeYWAgADy8/Nx\nc3PDxcWFH374gY4dO9KmTRtSU1MxmUykpKQQHh6Op6cnBQUF5OXlkZaWRkiI9a842dm/92VlZRkt\n7FleVpaRyxbuKOoilrU4lQkM9K7S69SK5ajnyhHPu5rnqjJVadf16/k2TfB//Xp+rZ13W9pka7tq\n+lqAmr0ebGmXpYRuMWFv3ryZY8eOmfugp06dyptvvomXlxfBwcG88MILaDQaRo4cSWxsLL6+vixe\nvBiAKVOmMG7cONzd3Vm4cKFNByOEEKJyFhN2VFQUUVFR5batXbv2tv2io6OJjo4uty0iIoKIiAgV\nmiiEEAJk4IwQQjgNSdhCCOEkJGELIYSTkIQthBBOQhK2EEI4CUnYQgjhJGQBA9FgyCREwtlJwhYN\nhkxCJJydJGzRoFhbPgtkCS3huKQPWwghnIQkbCGEcBKSsIUQwklIwhZCCCchDx1Fg2HLaixQ+yuy\nCGErSdiiAVFsXvmEh5VaapMQtpOELRoMFxcXm1c+kUEzwhFJH7YQQjgJiwn78OHDxMTEEBsby/z5\n8wFYtmwZsbGxxMXFUVRUBMD69euJiYlh8uTJ5sV3d+/eTUxMDGPGjCEzM7OGD0MIIeo/iwm7RYsW\nfPLJJyQmJnL16lX279/Pvn37SExMpEOHDiQnJ1NYWEhSUhKJiYkMHTqUpKQkAOLj41m+fDlxcXEk\nJCTUysEIIUR9ZjFhN2nSBDc3NwAaNWpEamoqvXr1AiAyMpJDhw6RkZFBaGgoWq3WvK2goAAPDw+8\nvLwICwsjNTW15o9ECCHqOZseOp44cYKsrCz0ej1abUmO1+l0GAwGDAYDOp2u0m0AJpOUSAkhRHVZ\nTdjXrl1jzpw5LFmyhGPHjnHx4kUAjEYjer0eb29vc791RdsAc5K3xM/PC1fXkifz2dk6K3v/zt9f\nR2Cgd6V/r4tY1uJYUtXXqRHLUc+VWrEc9bqypKb+Dev7ea/OOYe6Pe+WWEzYRUVFvPzyy7zyyisE\nBATQpUsXEhMTmTBhAikpKYSHh9OmTRtSU1MxmUzmbZ6enhQUFJCXl0daWhohIdanqczO/n0+4qws\no4U9y8vKMnLZwmCIuohlLU5lAgO9q/Q6tWI56rlSK5ajXleVqcl/w/p+3qt6zqHuz7ulhG4xYW/e\nvJljx46xaNEiAKZOnUrPnj2JjY2lefPmjB07FldXV0aOHElsbCy+vr4sXrwYgClTpjBu3Djc3d1Z\nuHChTQcjhBCichYTdlRUFFFRUeW2hYeHM3HixHLboqOjiY6OLrctIiKCiIgIlZophBBCBs4IIYST\nkKHpQtQxW9ealHUmhSRsIeqYLWtNyjqTAiRhC1HniotNeAXoaGyp3EujkSlfhSRsIeqe9WlfZcpX\nAZKwhahztkz7KlO+CpCE7fRsfWAF8tBKCGcnCdvJ2fLACuShlRD1gSTseiCw4x14N/e1uE/Ob9dq\nqTVCOAdnLKeUhC2EaJCcsZxSErYQosGy9u3U0b6ZSsJ2csXFJnJtmFks93KO1PEK4eQkYTs96zW8\nIHW8QtzKlpsdR7vRkYTt5Gyp4QWp4xXids43YEkSthCiQXLGAUsyvaoQQjgJqwn70qVLDB8+nLCw\nMPNiuj169GD06NGMGTMGg8EAwPr164mJiWHy5Mnm9Rx3795NTEwMY8aMITMzswYPQwgh6j+rCdvX\n15ePP/6Yrl27mrd16NCBlStX8sknn6DX6yksLCQpKYnExESGDh1KUlISAPHx8Sxfvpy4uDgSEhJq\n7iiEEKIBsJqw3dzc0Ov15badPn2aUaNG8fbbbwOQnp5OaGgoWq2WyMhIDh06REFBAR4eHnh5eREW\nFkZqamrNHIEQQjQQVXro+M0336DX65k1axbbt2/Hz88Pna5kmXedTofBYMBgMJi3AebuFFGiomGx\nMmGTEMKSKiXs0jvugQMH8vPPPzNgwABzv7XRaESv1+Pt7W3eBqDVWr6Z9/PzwtW1JDFlZ+ss7luW\nv7/O4rLwdRHLWhyAU6dO2TQsNjy8M6GhodVuky3tctRzpVYsOVd1E0vOlX2xLLErYSuKQn5+Pm5u\nbri4uPDDDz/QsWNH2rRpQ2pqKiaTiZSUFMLDw/H09KSgoIC8vDzS0tIICbE8Fj87+/c7y6wso4U9\ny8vKMnLZQvF7XcSyFqd0H1uGxTri8akZy9ZzVZttctRYtXmu1Iwl58q+WJYSutWEXVRUxIQJEzh5\n8iQTJkzgxRdf5I033sDLy4vg4GBeeOEFNBoNI0eOJDY2Fl9fXxYvXgzAlClTGDduHO7u7ixcuNCm\nAxJCCFExqwnb1dWVFStWlNu2du3a2/aLjo4mOjq63LaIiAgiIiKq10IhhBCAjHSsM844j4EQom5J\nwq4zzjePgRC2KC4u5uzZM+bfMzLSybtufeBc3vVMMjLSLe5TWazGPs3QaOt/NZUk7DrijPMYCGGL\ns2fP8NdF6/HyaWre5uUTZPV1Xj5BxG++AFywul9ZedcvAVidAO1W169f4/qZ77h8/VeKCq7j4tYY\nN+9m+LXtTeOmHTmzbT5F+dk8nfIujRq54evrS8eOnXj44aH07n1vuVgnTvzElWNfkJmfjan4Jq7u\nejz8WhF093C0rh52tcsSSdhCCNV5+TS1O4HWpgsXfmPWrNe4cUND4F0P4a6/A1DIu5JG5tEvaTdg\nOgDewX9g3itjadEimMzMi/zvf9t47bU4hg17lBdffAWAX345w+LFC3AL7MId4Y+hdXHjZu5VjBeP\no5iKVW23JGwhRIPz9tsL0Go1NOkag3eTNubtbrqmeLfobv5d4+KGXu9D06ZBNG0axN13d6V167Ys\nWjSPvn370717T/bv34O3tzcebe/D3bsZAI28/GkcqP6yYjJbnxCiQTEYrrNv3x4GDBiE1qXRbX93\naWS5CyMqKhpvbz3ffbcdgICAJuTk5HDj+rkaaW9ZcocthGhQzp37FUVRaN68BfxWZPfrtVotwcGt\nuHDhNwD69RtIcvIWdu78gmupW/DwDcYr4E70LXvg4lb5SOaqkIQthGhQFBUKr0rmRtIAJQl8woTJ\nnMpvj6Yol/zsDLJOf8fV1O0ER07B3dv6A1dbSZeIEKJBCQ4ORqPRcP581bowiouL+fXXjJI79DJc\n3HXoW3Yn6O5htLk/Do1GQ/bp79RospkkbCFEg6LX+9Cr1x9JTt6Kqbjwtr8XF+ZbfP1XX60jN9dI\nv34DKt3HpZEnLu7emIpvVru9ZUmXiBCiwXnppVeZNOkprhz+DM1dD+OuL6nuyLtymqzT39JuwGsA\nmIpvcu3aNTIzL3LpUibffpvMF1/8h0cffZyuXbsBsG7dFxw8uJ+Cazrc3DxQTIUYzh3kZk4m/u37\nqdpuSdhCCNWVDmaprfcqO0jHFs2bt2D27Pm8tuDfXDmxiaICA9pGXrjrmxHYaah5P+Ove3nhhX00\natQIX18/OnbsxLx5i8sNnOncuQt79nzP9dPfkv3zRjQujXDTBdIs/HH0LbqpdpwgCVsIobI2bdqx\n5OXfk15GRjpvJx2yOtox73omUx8Pp1Wr1pXuU1EsL5+mNPZpZnc7fX198Wl3f6UDfNoNmI4x+zzz\nJ/2RO++svKY6JKQDEydO4Yyyp8YHC0nCFkKoysXF5bYE5+VzwaZk1qpVa4vJ0Z5Y9ZE8dBRCCCch\nCVsIIZyEJGwhhHASVhP2pUuXGD58OGFhYeaVz5ctW0ZsbCxxcXEUFZUM7Vy/fj0xMTFMnjzZvPju\n7t27iYmJYcyYMWRmWp8PVwghROWsJmxfX18+/vhjunbtCsDVq1fZt28fiYmJdOjQgeTkZAoLC0lK\nSiIxMZGhQ4eSlJQEQHx8PMuXLycuLo6EhISaPRIhhKjnrCZsNzc39Ho9ULJq+rFjx+jVqxcAkZGR\nHDp0iIyMDEJDQ9FqteZtBQUFeHh44OXlRVhYGKmpqTV7JEIIUc/Z3Yedk5ODTqcDQKfTYTAYMBgM\nFrcB5u4UIYQQVWNXHbZGo8Hb25uLFy8CYDQa0ev1eHt7m/utK9oGJTNaWeLn54Wra8lyWNnZOov7\nluXvryMw0LvSv9dFLGtx1Iwl50r9NjlqrNo8V2rGknNlXyxL7ErYiqLQpUsXEhMTmTBhAikpKYSH\nh9OmTRtSU1MxmUzmbZ6enhQUFJCXl0daWhohIZaL4bOz88w/Z2UZLexZXlaWkcsWVh+vi1jW4qgZ\nS86V+m1y1Fi1ea6qE6s6i/AeOnTc4vvaugivs5yrilhK6FYTdlFRERMmTODkyZNMmDCBF198kZ49\nexIbG0vz5s0ZO3Ysrq6ujBw5ktjYWHx9fVm8eDEAU6ZMYdy4cbi7u7Nw4UKbDkgI4dzOnj3DK+tn\n0bhM4vHvA3DK4us8gI/On4LzluPfGiv3cg65Fx+yefTj3LlvsHnzRgAUNDTy9EXXrAsBHQahdXHD\nePEYWWn/42buZZTiIqad+ZyePXvx6quv2xS/JllN2K6urqxYsaLctrCwMCZOnFhuW3R0NNHR0eW2\nRUREEBERUf1WCiGcSuNAb7yb+9ba+xVctH1fjUbDPff8gdGjn2bOx/vRFOeReWQNpuJCvO/owm8/\nrKJJh0HomnUhz3CJmL5NSUtzjKIJmUtECNGgKIqCq2sj9HofXNx16Pw6kH/1DMaLx9FoXfD0CzZP\ni3qz8Cbh4d159NHH67jVJWSkoxCiwdFobvld64piKsLV3ZubxsvcMFyom4ZZIXfYQogGp+y6jvnZ\nGRjO/0jjwBB82/YmP+sX0ne8g6uHD66Nm/Dtt0aaN2+Jp6dn3TX4/0jCFkI0OHv3ppSsEnOzEBQF\nXVBnmnYZhtbFjRa9xnEz9yr5V0+Tk/kzn3++ms2bN/HBBx/j5+dfp+2WhC2EaHDCw3sQExPL4qQj\n6INC0WjK9w67NQ7ArXEALt4teGn4nUyb9hJffrmGceMm1VGLS0gfthCiwXF3d6dp0yBc3L1vS9a3\natKkCe7u7hQUWF6ctzbIHbYQQvyfKye3opiKaNy0A408/Sg0XmLZsgQKCgro3btvXTdPErYQQn25\nVkbzqf1eLtZ3M9NoNLdViZTyCriTa+m7uXjoPxTfyEHj0oiskHYsWPB3unYNV6W91SEJWwihqjZt\n2vHW0DfNv9f0IrwuYNcivK+99v8AOH369sEwXk3uxKvJnebfjdnnedXKIry1SRK2EEJVsghvzZGH\njkII4SQkYQshhJOQhC2EEE5CErYQQjgJSdhCCOEkJGELIYSTqFLCPnfuHJGRkYwePZrx48cDsGzZ\nMmJjY4mLi6OoqAiA9evXExMTw+TJk8ut7yiEEMJ+Vb7D7t27NytXruTDDz/k6tWr7Nu3j8TERDp0\n6EBycjKFhYUkJSWRmJjI0KFDSUpKUrPdQgjR4FR54MzevXsZNWoUDzzwAG3btqVXr14AREZG8tVX\nXxESEkJoaCharZbIyEhmzpypWqOFEI6rOovwZmSkW9zH1kV466sqJeymTZuydetWGjVqxDPPPENu\nbi4BAQEA6HQ6DAYDBoMBnU5XbpsQov47e/YMKS8+zx1eXuZtLwBcs/7a4oTd/GJln1tjXcjL46u7\nYmtl9OPBgwf461+nsHFjMnq9T42/362qlLDd3NzMP99///3odDoyM0s+9YxGI3q9Hm9vb3O/dek2\nS/z8vHB1LfmEzM7W2dwWf3+dxWXh6yKWtThqxpJzpX6bHDVWbZ6r6sTKztZxh5cXrXSW49ekmjo+\nX1+v/9vWGD+/8sesVrssqVLCzs3NpXHjxgAcPHiQ0aNHs2HDBiZMmEBKSgrh4eG0adOG1NRUTCaT\neZsl2dl55p+zsmx/QJmVZeSyhZnB6iKWtThqxpJzpX6bHDVWbZ6r6sSy59zUFEvH99xzk2jTpi25\nublc2PstGo0Wn+BeNLnrITQaDYqpiKunvsFw/hBFBQbGjfuMUaPG0KdPX5566ikAIiIiABg8OIrX\nXvt/qp53Swm9Sgn7wIEDLFmyBDc3N+655x7CwsLo2bMnsbGxNG/enLFjx+Lq6srIkSOJjY3F19eX\nxYsXV+WthBBCdVu3bqZ373sJDHscF0xkHlmDq4c3fu3u4+KhJPKzzhLYOZpiXJj4YCtcXFxp2jSI\nOXPe4vXXX+HTTz9Hr9fj7u5eq+2uUsLu27cvffuWn8x74sSJTJw4sdy26OhooqOjq946IYSoAU2a\nNOHJJ5/i+Pt70Pm14GbuZbLP7KRxUCdyfjtCiz+Mp3FgKMbs84SGdjTPIOjtXXL36+fnVyd92DJw\nRgjRoGg0Gjp3vrvcNk+/VhQVGCjITgeNBq+AOyt5dd2ShC2EaHAURanrJlSJJGwhRIOiKAo//XSs\n3Lb87AxcPfR4+LYCRSHvSlqFr23UqBEAxcWmGm9nRSRhCyEanCtXrrBq1ccU5WWT89sRss/swK/d\nvbjpAvFuHkbmkTXkXDhKUcF1Tp48wZYtmwBo1uwONBoNKSk7yc7OJj+/dldSlyXChBCqu5CXZ32n\nOnovjUbDoEGDMRoNXD6ShEarxSf4Hnzb9gGgWfjjXDm5lcvH11N0w8iyX5vw5JNPAxAY2JRx4ybx\n/vtLWbhwDg8++LB5jcjaIAlbCKGqNm3awT/eNf9e04vwgn2L8AK4urowevRYfsq/67YRkhqtK4F3\nPUTgXQ9hzD7P/FsW4X366Qk8/fQEu95PLZKwhRCqcvRFeBVFkYeOQgjhDDQaDRqNpq6bUSVyhy2E\naFD++c8EAE6fTq3jlthP7rCFEMJJSMIWQggnIQlbCCGchCRsIYRwEpKwhRDCSUjCFkIIJyEJWwgh\nnIQkbCG4HTBDAAAdoElEQVSEcBI1PnBm3rx5HD9+nE6dOjFjxoyafjshhKi3avQO+/jx4+Tn57Nq\n1SoKCws5evRoTb6dEELUazWasA8fPkzv3r0BiIyM5NChQzX5dkIIUa/VaJdITk4OwcHBQMnilamp\nlY/d79Gji/nnwsJCsgx5aLQuRIz8W4X77/58JoqpmOFfe5lXgQD44Ydjt+2bd/0Suz+fWWGc0vh5\n1y9V2p6y7WoVOanS9gC3tami9gAMHx5lPsaK2nNrmypqD8B//vPlbW0v255Spe06cuRkhXF69OhS\n7rxX1J6y7Ro+PKrceS9V9njLtqui86+YimHS9krbU8qW6yEl6bXbroVb21O2TZauh4rO563nv7RN\nkY/PqzCOvddDStJrt10Lpe0p2+7K2lPKluuhbJsqa4+t10Ppe1XWntL4t7apqtfDrW2qzvVZtl2V\nXQ9dB/3FYntKFRYWciXbwF33PV3h/id2foKimBj+td5qvrJEo9TgPIOrVq3C39+fwYMHs3XrVjIz\nMxk9enRNvZ0QQtRrNdol0q1bN3bv3g3A7t27CQ8Pr8m3E0KIeq1GE3anTp1wd3dn1KhRuLi4cPfd\nd1t/kRBCiArVaJeIEEII9cjAGSGEcBKSsIUQwklIwhZCCCchCVsIIZyEyxtvvPFGXTfCHqdPn8bf\n39/u1xUWFuLiUn5gQlZWFp6ennbF2bFjBx9++CFJSUls2LCBffv24eLiQuvWrW2OcfPmzXIrN+/e\nvZutW7eSk5NDmzZt7GoPQHFxMVptyWdvTk4OqampeHl54ebmZnesW23bto127dpVO87Zs2fZs2cP\nQJX+/cq6evUqXl5e1W7TgQMH2Lx5M7m5ubRq1cqu1w4bNoysrCyCgoLw8fGpdlvquwMHDpCXl4e/\nvz979+5ly5Yt3Lhxwzywztnt2LGDffv24ePjg16vN29fs2YNnTp1Uu19nK5K5P777+d///ufzfvv\n2bOHV155hRs3btC5c2dmz55tvkiGDRvGunXrbI41Z84c0tPTGTZsGEFBQQBcvHiR//73v7Ru3ZrX\nX3/dpjhDhgzh008/xcfHh2XLlpGcnMx9993H/v376dy5M3FxcTa3adOmTcyePRtvb2+mTZvGvHnz\nCA4OJj09nTfffJP77rvP5lhbt26l9HLQaDQoisLs2bMp/UwfNGiQzbFGjx7NkiVL8Pf3Z926dcTH\nx9OzZ0+OHDnCyJEjGTNmjE1xrl27Vu53RVF49NFHWbt2LQC+vr42t2nEiBGsWbMGgP/85z+sWrWK\nBx54gF27dtGvXz/+/Oc/2xyrf//+/OlPf+Lrr78mICCAIUOGMHjwYPN1YQ+j0ciHH37I1q1buXDh\nAo0aNaJVq1Y88cQTPPLIIzbHMRgMvP/++yQnJ3P16lU0Gg3+/v4MHDiQSZMmlUsk1gwbNoxBgwYR\nFRVl94fZrd5++2327t1LcXExf/jDHzhw4AB9+/YlJSWFfv36MWHCBJtjqXmMv/zyCwkJCfj4+PD0\n008zc+ZMDhw4QKtWrZgzZw5hYWE2H9/Bgwfp1KkT3377LWPGjDFf3/bmGKsUB/Tmm29W+l94eLhd\nsYYPH66cOnVKMZlMytdff60MHDhQOXjwoKIoihIdHW1XrAceeKDC7SaTSRk4cKDNcR5++OFy7cvP\nz1cURVEKCwuVqKgou9oUFRWlXLp0ScnIyFDCwsKUtLQ0RVEU5dy5c8qjjz5qV6y77rpLmTRpkjJt\n2jRl2rRpyquvvqqEh4ebf7fHrceYlZWlKIqi5OXl2XWMHTp0UPr161fuv06dOin9+vVT+vfvb1eb\nyv57Dx8+XLl69aqiKIqSm5tbrr32xDKZTMq+ffuUWbNmKZGRkcqTTz6prF692q5Yf/7zn5U1a9Yo\nv/32m7J8+XLln//8p3LmzBnl5ZdfVt5++22b44wdO1ZJSEhQLl26pJhMJkVRFCUzM1NJSEhQxo4d\na1eb+vXrpyxYsEDp27ev8sgjjygfffSRcvHiRbtilBo8eLBSWFio5OXlKeHh4YrBYFAURVHy8/Pt\nvt7VPMbHH39cWb16tfLBBx8offr0UTZu3Kjk5+cru3btUkaMGGFznIcffli5efOmoiiKcv36dWX8\n+PHKnDlzFJPJZHeOscYh+7DXrl1LSEgInTt3pkuXLub/OnfubPfX/MLCQkJCQtBoNDz44IMsXbqU\nadOmkZycbHe73N3dOXLkyG3bjxw5goeHh81xGjduzMmTJXN6+Pn5UVBQAEBRUZH5DtdWLi4uBAYG\nEhwcTMuWLbnzzjsBaNGiBUVFRXbFWr16Nfn5+XTp0oV58+axYMEC/Pz8mD9/PvPnz7crlqurKxcv\nXgRKjre068nNzQ2TyWRznFdeeYW2bduydOlStm/fzvbt2wkKCmL79u1s27bNrjYVFxdz7do1srOz\nMZlM5q4ZLy+v27rLbKXRaLjnnnuYPXs23333HRMnTrR7krPz58/z6KOPcscddzB27Fi2b99O27Zt\nWbBgAVu2bLErzqRJkwgMDDR3tzVt2pRJkyZx/vx5u9qk1+t59dVX+fbbb5k2bRq//PILjzzyCKNH\njyYpKcmuWI0aNcLV1RVPT09atWqFt7c3AB4eHuauPFupeYz5+fk8/vjjTJgwAQ8PDx566CE8PDzo\n3bs3N2/etDlOcXGxeX4QvV7Pv//9b4xGI3/9618pLCy0q03W1Ph82FXRpUsXQkJC6NGjx21/e++9\n9+yK1ahRIy5fvkxgYCAAISEhfPzxx0yaNIlff/3Vrljz58/njTfeIDc3l2bNmgElXSI6nc6uhDZ7\n9mxefvllOnToQEBAAI8++ij33HMPp06dYtKkiieXssRkMqHVasu1oaioyO6EHRYWxkcffcTKlSsZ\nM2aMXV0zt5o+fTrjx49n0KBBhISE8PTTT9O7d29++OEHu77mjxs3jsGDBzN//nyaNWvG888/X+U2\n5ebmmt9bo9Fw6dIlmjZtitFotDtWRc8aXF1due++++zqhoKSD4wDBw7Qs2dPkpOT8fPzA7A7mTVv\n3pwPPviA4cOH06RJEwAuX77Ml19+yR133GFXrFKlH0j33HMPM2fOJCUlha+//prHH3/c5hhubm7k\n5+fj6enJl19+ad5uMBjq9BhLEz6U3FRU9jdrgoOD2bdvH7169QJKroP58+fzj3/8g61bt9rVJmsc\nsg/72rVruLu72/1AsCLff/89/v7+3HXXXeW2GwwGPv30U5555hm7Y166dInMzEwAgoKCaNq0qd0x\nioqK+P777zl79ixFRUXccccd9O7d2+4HWEeOHCE0NPS2O/xz587xww8/EB0dbXfbADIzM5k7dy7H\njx+3+062lMFg4KuvviI9PZ2ioiKaNWvGgAEDzN8C7JWcnExCQgLnz58nJSWlSjEqkp+fz5UrV+rs\nAdiJEyeYMWMG6enptG/fnnnz5tGuXTuysrLYsGGDXf3977//Ptu3b+fq1asABAQE0L9/fyZNmmRX\nf/8LL7zAO++8U6XjudWNGzdwd3e/bXtWVhaXLl2iY8eONsdS8xjDwsLMxQIZGRnl+uozMjI4fPiw\nTXFKvyFX9C374sWL5ps7NThkwi5148YN0tPTKS4upnXr1tWqDFAzVimj0Uh6ejrBwcF2PeyoaVev\nXiUgIKCum1EjCgoKyMjIIDQ01O7XHjlypNIHSevWrWPYsGE2xzp8+DDt2rXD29ub/Px83n//fX76\n6Sfat2/P5MmTzV/7xe3Onj3LiRMnaN++Pe3bt6+zdpw7d87i31u2bGlXvMLCwtumdM3Kyqp2VVRZ\nDpmwCwsL+cc//sEXX3xB8+bNAbhw4QJRUVG8+uqrZGRk2HyXpmasN954w1wxceDAAeLi4spVZNx/\n//02xdmxY4f5a7PBYGDBggUcPXqU0NBQpk+fbv6qZws1qyimT59u8e/2dPuoFevQoUPceeedqiTG\nIUOG0L17d6ZOnWr+gD158iRvvvkmer2e+Ph4m2M99NBDrF+/HldXV15//XU8PT3505/+xO7duzl5\n8qTdXXdlHThwwPzNqU+fPja/7pNPPuGBBx6ocvdHWWp+IKlVMQTqXg+3ysrK4sCBAzRv3pwuXSqe\n37sialaiWeOQCXvu3Lnk5eUxffp0dDodUHI3u2DBAoxGI2lpaWzYsKHWY5U9+aNHj2batGl07tyZ\nX3/9lb/+9a/mJGlPnNdee43AwEBGjhxJcnIy+/btY+nSpTbFAejYsaP5g6hUZmamubxs+/aKJ4Sv\nyObNm83lfKV9eBcuXGDFihUUFxezc+fOKsWCkj7BqsRSMzEWFhayfPlykpKSmDJlCqdOnWLHjh1M\nmzaNfv362RwHYPDgwXz99dcADB8+vFzf7NChQ1m/fr3NsdQqN+zRoweenp4EBwcTFRXF4MGDq3x3\np+Z5j4qKMv8/9sgjj/Dhhx/i5+dHfn4+jz32GF999VWdtGvSpEnExcURGhrKpUuXGD58OHfffTcZ\nGRk89thjPP300zbFeeSRR1i4cCHt27dny5YtvP3227z11lt069atYZT1DRw4UCkuLr5te1FRkdKz\nZ09zWV5txypbonNruY495Ttl9x0yZIi5PKn0d3t8+OGHyrhx45Sff/7ZvK1fv352xahIenq68tpr\nrymDBg1SVq1apdy4caNOYj344IPmn4cNG1bub/aeq1Lvv/++0qFDB+Xee++tcqnaX/7yF+Xzzz9X\nFEVRpk2bphw+fFhRFEU5c+aM8sgjj9gVS61yw+joaKW4uFjZuXOnMm3aNOUPf/iDMm7cOGXt2rVK\nTk6OXW1S87xHR0crFy5cUBRFUZ588klzGWtRUZHy0EMP1Vm7yr53fHy88vLLLyuKoig5OTl2lRve\nuu+pU6eUQYMGKd98843qZX0OWSWi1WorfHrs4uKCv78/3bp1q5NYZ86cYciQIUBJ/9f169fx8fGh\nuLjYroqMrKwsPvroIxRFua1CQbHzC4+aVRRQMpI0Pj6en376ifHjxzN79mxcXat2magRKyQkhDVr\n1jBixAg6duxo7of+5ZdfKlwCypLSrisoGXC0Y8cORo0axeTJkxkxYoRdsebMmcPcuXOJj4/H39+f\nJ554gqCgIO644w7mzp1rV6zSckNFUapdbqjVaunTpw99+vTh5s2b7Nixg40bN7JgwQL27t1rcxw1\nz7taFUNqt6vstZiSksJjjz0GgE6ns6tKRM1KNGscsktkypQpDBo0iOHDh5fbvm7dOrZs2WJXX6Oa\nsW59SNG0aVPc3NzMfV+2jgT85z//We6CeOKJJwgICODSpUssXryYt956y+Y2lVXdKornn3+e48eP\nM27cOB588EG0Wm25dtrTH65WLIPBwNy5czlw4AD+/v789NNP5sQ4c+ZMuyoMHnjgAV566SUGDx5s\n3paZmcn8+fO5ePEiq1evtjlWqZycHM6dO2eugin9n9Ye/fv3N/+s0Wj47LPPzOWGo0aN4r///a9N\ncSx9/c7Ly7PrQbua5700nhoVQ2q2689//jN9+vQhKCiIGTNmkJycjI+PD/n5+YwYMYKNGzfaFKey\nSrScnBw+/fRTpkyZYtcxWuKQCfvixYs899xzeHh40LlzZwCOHz9Ofn4+//rXv+wqk1EzVm5u7m31\nmqVuLQuqK/n5+fz6669VqqIomzgqYk9/uJqxQJ3EaDQazc8xbvX999/Tu3dvu2OWys/PJy0tjRYt\nWqhWFWBvueGZM2dUmfelLDXOe01Qo11Xrlzh3Xff5fLly4waNcr8gHfPnj0cP36c8ePH2xQnOTmZ\nixcv8uSTTwIlzySysrLQaDTExcWVu0GoLodM2FDSNbBnzx5SU1PRaDS0b9+eiIiIOo01cOBAXnrp\nJR566CHztoKCAuLj49m4caPNoycVRWH//v3o9Xo6duzIpk2b2L9/P61btyY2Ntau0ZyKopgnnalu\nLEekZine8ePHLf699APdFtu2bWPu3Ln4+PjwwgsvMHv2bJo0acK5c+eIi4uz+6t+RXJzc1m1alWV\nBlNByV316dOnVf0QqUqb1Kw+qkx1z1V1xMTE8Pe//9388D86OpoVK1aQn5/PtGnT+OSTT1R7L4fs\nw54zZw7du3enR48eVU7SNRFr+fLlzJ49mzVr1jBr1izS0tJ46623GDBggM1fW6FkpOOpU6e4ceMG\nbdu2JS8vj3vvvZcffviB6dOn8/bbb9dJLDUTmlqxZsyYYbEUz56EvWDBAvPPx44du610a+XKlTbH\neuedd/jwww/JyclhzJgxfPXVVwQHB3P16lWeeuopuxJ2ZmYm77//PhkZGYSEhPDss8/y+eefs3z5\ncrsm3FLzQ0StNgH07dvXYsWQPdRs17PPPsugQYMYOHBgpd+cbXHz5s1ylVo9evTAz8/PXAmjJoe8\nw165ciU//vgjhw4dQlEUunfvTrdu3ejevTsdO3a0azirmrFKffDBB/zjH/+gSZMmLFu2zO7uh8GD\nB7Np0yZu3LjBfffdR0pKCq6uriiKwpAhQ2wuM1Q71ujRo80/VzehqRVLzVK8sqpbblX29UOGDClX\nmmZv7KeffpoePXrQtWtXdu7cybZt2+jatau55NNWQ4YM4d133630Q8Sea0GtNt0qIyODhIQEDhw4\nwFNPPcWIESPs+haoZrvuvfdeunXrxp49e4iIiCAqKoq+ffva/a104MCBlX67HjBgQJVHCldI1ZqT\nGnDx4kXl66+/VubNm6cMGDBA6datW53FunnzpvLvf/9b6d+/v/LZZ58pU6ZMUcaMGWOeIc9WapUH\nqh1LrdfWRCw1SvHUbFNUVJSSnZ2tZGVlmX8u+7s9hg4dWu73e++9VykqKrK7TWWP6dY22Hu8arWp\nVFpamjJ16lRl8ODBypo1a5TCwsIqxVGzXaXnJCcnR/nyyy+V8ePHK7169VKmTZum7Ny50+Y4L730\nUoUzNCYmJiovvvhildpWGYfsEoGSvtmTJ09y8OBBfvzxR06fPk3r1q3t+gqsdqzhw4dzzz33sG7d\nOry9vYmJieHbb7/lmWeeYdCgQUydOtWmOGXL+sr+XPo3e6gZyxGpWYqnprITSQHV6rM2mUzmEauK\nouDj40NOTo7577ZW1JQtDwTKxbS360GtNkH5iqHp06ej1WrLlbPaE0vNdpXS6XQMGzbMvCjFli1b\nSEhIsHmU6fTp03n22WfZsGGDebGCn376iRs3bvCvf/3L7vZY4pBdImPHjsVoNHLXXXcRFhZGt27d\naNeunV21kTUR6+jRo9x99923bS998Pjiiy/aFKdsWZ9SZlRh6c/PPfeczW1SM1ZZao7Qqk4sNUvx\n/va3v5l/3rRpEw8//HC5flVbF6BQm1oVNY5a5eOosUaNGsWqVats3t8SRcUiCUscMmHPmjWLEydO\n4OHhQdeuXenWrRvh4eFVetKtZqxSv/76K2lpaRQXFxMSEmLX8mDWWKqKqOlYaiY0tWKpWYq3du1a\nix9ut9bqW6LmA1ohbOWQCbtUTk4Ohw4dMv+XnZ1N+/btqzSwRI1YRqORGTNmcOzYMXOR/M8//0xo\naCiLFi3ixIkT9OzZ0+62paamsmHDBjZt2oS3t7fNc5KoHUvNhKZWLEdNjGo+oFXrGB2xyseRY0HJ\nQJydO3eap0tu1qwZffr0cajZN8ty6IR98+ZNjhw5wsGDBzl48CCHDx8mICDArqfdasZ69dVXadmy\nJc8++6y5usRkMrF06VL27t3LtWvXbJ7I5tdff2XTpk1s2LCBRo0acf78eb744gu7p3RUO5ajUTMx\nqlXGdavqdh+pdYyOWOXjyLHWrVvHe++9R2RkZLkFSb7//nuee+45u25Qao2qjzBVMnfuXGXEiBFK\nz549ldGjRyt///vflf/973/K9evX6zSWpXUb//jHPypnzpyxKc5jjz2mDB8+XElISFAyMjIURan6\nhE1qxnrmmWeUdevWKUajsUqvr6lYpapb2dGnTx/lL3/5i3LPPfcozz//vLJ169ZqTWqlVrtqIpYj\ntsnRYg0aNKjCPHDt2rVK12+taw5ZJdKiRQuGDh1Kx44dqzzxUE3EsvSgUqfT0bZtW5viBAQEkJqa\nypUrV7h69Wq1VjpRM9aRI0dwcXFh7ty51apLVTuWWgICAnj33XcxGo0kJyeTlJTE66+/Tv/+/Xn4\n4Yftmnta1F9VKUioLQ6ZsHv06IFGozEvVHsre/qp1IwVHh7Oe++9x7PPPluuT3bp0qV2zfq3dOlS\nDAYDW7du5d133yU9PR2DwcDhw4fp2rWrzXHUjqVmQnPk5FjdMi4o/1A1MzOTOXPmOETFibDd5MmT\neeSRR4iMjDQv+nDhwgW+//77Ki0dWBscsg/bUfu8cnJymDFjBsePHy/30LFTp07MmzevyqtdXLly\nha+//pqNGzdy4cIFvvvuuyrFqW6sivpiSxPapk2b7DpXasVSs3JFzTIuNR/QqnWMjljl48ixAK5f\nv86uXbvKrdHap08fu9dWrS0OmbDLcpR64LLS09NJS0sz11uqOUvfuXPnVHtYaG8sNROaWrHUTIyO\nSq1jdMQqH0eOZelbsZubG61bt+aFF14gMjLS5pg1TRK2HY4dO1buAoHy/V22dq+oWa1QU5UP9ZVa\nZVxy3uu3oqIi0tLSmDp1qs3zYtcGlzdKV5V1UKtXryYmJsYhYk2dOpVdu3axa9cuEhISSE1NNf++\na9cum4cnL1iwgPz8fObOncuxY8dwdXUlODjY7tVF1I4FJQktOTmZXbt28eOPP3LhwgWCgoJwd3ev\nk1jPPvssRUVFBAcHV/uB5bp165g6dSparRYPDw+Kioo4ceIEixYtQq/X3zYBvSVqnne1jlHNc9UQ\nYlmi1Wpp0qQJWq22wtHNdcUh77Aduc+rVHXu1ktfW/pAbsOGDRw9erRKD+TUjKVmXapasdSaUQ3g\nT3/6E59//vltd9PXr19n5MiRbN261eZYap53tY5RzXPVEGI5I4esEuncubO5q6Hsz2X7rOoiltrU\nqFZQM1Z8fDxr166tNKHZk7DVilUb1SbVuQ7UOO9qHaOjVvk4aiynVMN13vVWdYr2Y2NjVWuHmrEq\nG0hw/fp1uwcSqBWrovN89epVJTExUXnyySftatPatWuVAQMGKDNnzlSWLl2qLF26VJk5c6bSv39/\nZc2aNXbFUvO8q3WMap6rhhDLGTlkl4ijPpRz1Nne1PLll1/yr3/9q9K61EcffbTWY6lZuQKOWcal\n1jE6YpWPI8dyRg6ZsB21z0vNkiI1J51RM5aaCc3RkqPaZVzONnGQcH4OmbAd9aGcWhzx4R6om9DU\njFUbidHeMi61Jw5S6xgd9UbAUWM5G4dO2GXV9Yg7UK97Rc1qBTVjWaJmXao9sWp7RrXPPvuMJ554\nwup+ap53tY7RUW8EHDWWU6qjvnOLHPWhnFqzvTniwz1bJSYm1mosR51RrTb+De09RjXPVUOI5Ywc\nsqxPzYcKasZSq6RIzUlnansCG1vuPmsjVl2XZNbGeVfrGNU8Vw0hliNzyC4RcMw+LzW7V+rzwz01\nqVm5oja1zrtax1gap3fv3ubugupW+dTnWM7IIRO2o/Z5qVVS5KgP9xyVI34gqX3e1TrGa9eusWvX\nLi5dumSOExQUxIYNG7B3FoqGEMvZOGTCdsaHcmqpq4d7jsoZP5DsPe81cYzHjx9nw4YNbN68mZYt\nWzJo0KByUw3boyHEchYO2YddGUfo86rpkiJXV1c6duzIk08+6VCx6sqPP/5Y6d8c9QPJ3vOu1jGe\nOXOGjRs3smnTJvz9/XnwwQdRFMWurrqGFMsZOWTCdtSHchV1r+zZs4e3335b9ZIiR3m458gc/QNJ\njfNuzzE+9NBD3H///Xz44Yc0b94cgI8++qhK79sQYjkjh+wSAcd8KOds3SuiYSkdGHbkyBH69OnD\n4MGDmTFjBtu3b5dY9YRDJmxHfShXWcI2GAyMGDFCErZwCLm5uWzbto2NGzeyd+9eoqOjeeCBB6o0\nqrchxHImDpmwLanLh3KOXGImREWuXbvGli1b2LhxI5988onEcnJOl7BL2TqEWO1YjlhiJoRoGJw2\nYdcFZywxE0LUH5KwVeKoJWZCiPpDW9cNqC8cvcRMCOH85A5bCCGchNxhCyGEk5CELYQQTkISthBC\nOAlJ2EII4SQkYQshhJP4//M19PynLW6RAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3d21250>"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot_dict = defaultdict(list)\n",
      "aug_genes = defaultdict(list)\n",
      "for genome in genomes_order:\n",
      "    count = 0\n",
      "    \n",
      "    for gene in pN_pS_Aug[genome]:\n",
      "        if float(pN_pS_Aug[genome][gene][2]) > 1 and not np.isinf(float(pN_pS_Aug[genome][gene][2])):\n",
      "            count += 1\n",
      "            aug_genes[genome].append(gene)\n",
      "            \n",
      "    plot_dict[\"CDS\"].append(cds_per_genome[genome])\n",
      "    plot_dict[\"PS\"].append(count)\n",
      "    plot_dict[\"pct\"].append(float(count)/cds_per_genome[genome]*100)\n",
      "    \n",
      "df = pd.DataFrame(plot_dict, index=genomes_order)\n",
      "df.plot(kind=\"bar\",stacked=True,legend=True)\n",
      "    \n",
      "print df"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "           CDS   PS       pct\n",
        "J07HWQ1   3584   62  1.729911\n",
        "J07HWQ2   3856   46  1.192946\n",
        "J07HQX50  2872   19  0.661560\n",
        "J07AB56   1411   79  5.598866\n",
        "J07AB43   1678   87  5.184744\n",
        "J07HN4    3230  196  6.068111\n",
        "J07HN6    2914   64  2.196294\n",
        "J07HX64   3049  189  6.198754\n",
        "J07HX5    2139  137  6.404862\n",
        "J07HB67   2847  168  5.900948\n",
        "J07HR59   1841   33  1.792504\n",
        "A07HB70   2514  156  6.205251\n",
        "A07HR67   2891  171  5.914908\n",
        "A07HN63   2507   50  1.994416\n",
        "A07HR60   2861   87  3.040895\n",
        "J07SB     1641   91  5.545399\n",
        "G22       3525    0  0.000000\n",
        "\n",
        "[17 rows x 3 columns]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAEPCAYAAABm//5NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclOX++P/XDMjmMGwipuKWoKkhLnkCLXPJk4WipUWY\nluvR6nQqqTTTX/ZxTTsnOx05lNkqycnMY2pqaCc1XDPXUkET0hQXkGFYFJj79wdfJlCYBW5gBt7P\nx6PHA27uec913969557rfl/XpVEURUEIIYTD09Z3A4QQQthGErYQQjgJSdhCCOEkJGELIYSTkIQt\nhBBOQhK2EEI4CZsS9kcffURsbCwAK1asIDY2lri4OIqLiwFYv349MTExTJ06FaPRCMDu3buJiYlh\n3LhxZGZm1lLzhRCi8bCasG/cuMGJEyfQaDRkZWWxb98+EhMT6dSpE8nJyRQVFZGUlERiYiLDhw8n\nKSkJgPj4eFauXElcXBwJCQm1fiBCCNHQWU3YX3zxBSNGjEBRFI4ePUqfPn0AiIyM5NChQ2RkZBAa\nGopWqzVvKywsxMPDAy8vL8LCwkhNTa31AxFCiIbOYsIuKipi//793H333QDk5uai0+kA0Ol0GAwG\nDAaDxW0AJpOpttovhBCNhsWE/d///peoqCjz797e3uY+aqPRiF6vt7oNQKu13lVeXFxSrQMQQojG\nwtXSH8+ePcsvv/zC6tWrSUtL4+jRoxw9epRJkyaRkpJCeHg47dq1IzU1FZPJZN7m6elJYWEh+fn5\npKWlERISYrUh2dn5VvcJDPTm8uVc24/OyWI5YpscNZYjtslRYzlimxpDrOrGCQz0rvJvFhN2XFyc\n+ecxY8bw7LPP8v777xMbG0vLli0ZP348rq6ujB49mtjYWHx9fVm6dCkA06ZNY8KECbi7u7N48WK7\nGy2EEKIiiwm7vFWrVgEwefJkJk+eXOFv0dHRREdHV9gWERFBRESECk0UQggBdiRsZ1VSUsKOHd/d\nst3Hx4ucnIrdMPfeOwAXF5e6apoQQtilwSfss2fP8Oa2ZXj5N7W4X35WHm3atOX22633twshRH1o\n8Am7pMSEV4COphY68gHQaCgpkfJDIYTjavAJGxSuHWjPdW9/i3sV5GbBQ7L4jhDCcTX4hO3i4kJA\n6zvQ+bWyuJ8x+7z0XwshHJrM1ieEEE5CErYQQjgJSdhCCOEkJGELIYSTkIQthBBOQhK2EEI4CUnY\nQgjhJByyDrukpISzZ8/csj07W0dWlrHCtnbtOkj9tBCiUXDIhH36dBrrnp5MoIenxf0uFxYwYvn7\nhIZ2qqOWCSFE/XHIhA0K393tY8OETW6MQIaTCyEaB4dM2C4uLgR2vg3vlr4W98v9/Zp0hwghGg2H\nTNglJSbybFhaJ+9yrsywJ4RoNCwm7NTUVGbPno2LiwsdO3Zk7ty59OrViy5duqDRaHj33XfR6/Ws\nX7+exMRE8xJhOp2O3bt3s2zZMtzc3FiyZAlBQUF2NEtm2BNCiJtZTNjt27dn9erVAMycOZOff/6Z\nTp068emnn5r3KSoqIikpicTERDZv3kxSUhITJ04kPj6elStXkpaWRkJCAnPmzLG5UTLDnhDCmVRW\n2VYbVW0WE7ar6x9/vn79Ot7e3pw+fZoxY8bQs2dPpk+fTnp6OqGhoWi1WiIjI5k9ezaFhYV4eHjg\n5eVFWFgYS5YsqXYDhRDC0Z09e4a/LVmPl0/zKvfJz7nEspeG12hVK6t92Nu2bePtt9+mW7duBAcH\n8+2336LX65kzZw7bt2/Hz88PnU4HgE6nw2AwYDAYzNsATCbpZxZCNGxePs2t9grUlNWEPWjQIAYN\nGsS8efP44Ycf6Nu3LwCDBw/ml19+YdCgQRiNpbf9RqMRvV6Pt7e3eRuAVmt9QKWfnxeurqVfFbKz\ndVb2/oO/v45AC8t/qRnLkuq+rrbiNIZYjtgmR43liG1qSLFszTM1yTFgJWHfuHEDNzc3AJo2bUph\nYSEmkwmtVsuPP/5I586dadeuHampqZhMJlJSUggPD8fT05PCwkLy8/NJS0sjJMT6V4Ds7D9WML+5\n38eSrCwjly1UlKgZqyqBgd7Vel1txWkMsRyxTY4ayxHb1NBi2ZpnbMkxlhK6xYS9c+dOPvroIxRF\noXXr1gQFBTFq1Ci8vLwIDg7m+eefR6PRMHr0aGJjY81VIgDTpk1jwoQJuLu7s3jxYpsORtSvkpIS\nduz47pbtPj5e5OTkV9h2770D5IGvEHXMYsIu6w4pb+3atbfsFx0dTXR0dIVtERERREREqNBEUVfO\nnj3Dm9uW2TDCNI82bdrW6OGJEMJ+DjlwRtQfW0eYCiHqnkyvKoQQTkISthBCOAlJ2EII4SQkYQsh\nhJOQhC2EEE5CErYQQjgJKesTDq+yAT0ymEc0RpKwhZmjLhxhy4AeGcwjGgNJ2KIcx104wtqAHhnM\nIxoDSdjCTBaOEMKxyUNHIYRwEpKwhRDCSUjCFkIIJyEJWwghnIQkbCGEcBIWq0RSU1OZPXs2Li4u\ndOzYkblz57JixQq2b99Oy5YtWbRoEa6urqxfv57ExETzijM6nY7du3ezbNky3NzcWLJkCUFBQXV1\nTEIIFciAJcdjMWG3b9+e1atXAzBz5kyOHDnCvn37SExM5P333yc5OZlBgwaRlJREYmIimzdvJikp\niYkTJxIfH8/KlStJS0sjISGBOXPm1MkBiYbHlgE9dT2YpzGQAUuOx2LCdnX948/Xr1/n6NGj9OnT\nB4DIyEi+/vprQkJCCA0NRavVEhkZyezZsyksLMTDwwMvLy/CwsJYsmRJ7R6FaOCsD+ip68E8jWX9\nSxmw5FisDpzZtm0bb7/9Nl27dsXHx8d84el0OgwGAwaDAZ1OV+U2AJNJ7nxE9dkyoKeuB/PI+pei\nPlhN2GUL8c6bNw9PT08uXrwIgNFoRK/X4+3tjdForHIbgFZr/dmmn58Xrq6l/8NlZ+us7P0Hf3+d\nxWXh1YxlSXVfV1txqhNLzXNVUlJCcnKy1TiDBw+2mmhtbVdd/vtlZ+tsXv/SWa8rRzzvjhqrLs4V\nWEnYN27cwM3NDYCmTZuSm5vL/v37mTRpEikpKYSHh9OuXTtSU1MxmUzmbZ6enhQWFpKfn09aWhoh\nIdbvLrKz//gamZVltLBnRVlZRi5b6N9UM1ZVAgO9q/W62opT3VhqnqvTp1OZtXah1f5PH59Aq3ef\ntrarLv/9GsN15YjnvaSkhLNnz9yy3d9fV6G97dp1qPY3rtq8Hmw5V5YSusWEvXPnTj766CMURaF1\n69Y899xzXL58mdjYWFq2bMn48eNxdXVl9OjRxMbGmqtEAKZNm8aECRNwd3dn8eLFNh2MaFik/1Oo\n7fTpNNY9PZlAD88q97lcWMCI5e8TGtqpDltWNywm7LLukPImT57M5MmTK2yLjo4mOjq6wraIiAgi\nIiJUaqZwNlLZIWqHwnd3+1j55ubGCOp2Nsm6IrP12cHWygBnrgpQj+NVdgj7OOKHrouLi03f3Brq\n/3+SsO1w+nQac9cswEPvVeU+hYZ83msV3CC/jtnDESs7hL3kQ9fRSMK2i4LLtbtxLan6AnbJzYIG\n+nVMNC6O+KHriHf9dUkSth0c8QIW9cNRl1Nr+Br3Xb8kbCGqxXGXU2vIGvtNkyRsIapBllMT9UGm\nVxVCCCchCVsIIZyEJGwhhHASkrCFEMJJSMIWQggnIVUiQjQgsqxXwyYJW4gGRJb1atgkYQvRwMi0\ntg2X9GELIYSTkIQthBBOwmLCPnz4MDExMcTGxrJw4UIAevXqxdixYxk3bhwGgwGA9evXExMTw9Sp\nU81rOe7evZuYmBjGjRtHZmZmLR+GEEI0fBb7sFu1asUnn3yCm5sbcXFxnDp1ik6dOvHpp5+a9ykq\nKiIpKYnExEQ2b95MUlISEydOJD4+npUrV5KWlkZCQgJz5syp9YMRQjgeWxf+AKlescZiwm7WrJn5\n5yZNmuDi4sLp06cZM2YMPXv2ZPr06aSnpxMaGopWqyUyMpLZs2dTWFiIh4cHXl5ehIWFsWTJklo/\nECGEY7KlcgWkesUWNlWJnDhxgqysLG6//Xa+/fZb9Ho9c+bMYfv27fj5+aHTlS7xrtPpMBgMGAwG\n8zYAk0nmAxaiLjjiBP8lJSa8AnQ0tbAaOAAajcwdboXVhH3t2jXmzZvHsmXLANDr9QAMHjyYX375\nhUGDBpn7rY1GI3q9Hm9vb/M2AK3W+rNNPz8vXF1LvwplZ+us7P0Hf3+dxWXh6yOWtTiWVPd1asRy\n1HPliOddzXNlib2vu3KlqU0T/PuPbVpn592WNtnarrq4FqD2roeatstiwi4uLuall17i5ZdfJiAg\ngIKCAtzc3HBxceHHH3+kc+fOtGvXjtTUVEwmEykpKYSHh+Pp6UlhYSH5+fmkpaUREmL9K0529h99\nWVlZRgt7VpSVZeSyhTuK+ohlLU5VAgO9q/U6tWI56rlyxPOu5rmqSnXalZNTYNME/zk5BXV23m1p\nk63tqu1rAWr3erClXZYSusWEvXnzZo4dO2bug54+fTpvvPEGXl5eBAcH8/zzz6PRaBg9ejSxsbH4\n+vqydOlSAKZNm8aECRNwd3dn8eLFNh2MEEKIqllM2FFRUURFRVXYtnbt2lv2i46OJjo6usK2iIgI\nIiIiVGiiEEIIkIEzQgjhNCRhCyGEk5CELYQQTkISthBCOAlJ2EII4SQkYQshhJOQBQxEoyGTEAln\nJwlbNBoyCZFwdpKwRaNibfkskCW0hOOSPmwhhHASkrCFEMJJSMIWQggnIQlbCCGchDx0FI2GLaux\nQN2vyCKErSRhi0ZEsXnlEx5S6qhNQthOErZoNFxcXGxe+UQGzQhHJH3YQgjhJCwm7MOHDxMTE0Ns\nbCwLFy4EYMWKFcTGxhIXF0dxcTEA69evJyYmhqlTp5oX3929ezcxMTGMGzeOzMzMWj4MIYRo+Cwm\n7FatWvHJJ5+QmJjI1atX2b9/P/v27SMxMZFOnTqRnJxMUVERSUlJJCYmMnz4cJKSkgCIj49n5cqV\nxMXFkZCQUCcHI4QQDZnFhN2sWTPc3NwAaNKkCampqfTp0weAyMhIDh06REZGBqGhoWi1WvO2wsJC\nPDw88PLyIiwsjNTU1No/EiGEaOBseuh44sQJsrKy0Ov1aLWlOV6n02EwGDAYDOh0uiq3AZhMUiIl\nhBA1ZTVhX7t2jXnz5rFs2TKOHTvGxYsXATAajej1ery9vc391pVtA8xJ3hI/Py9cXUufzGdn66zs\n/Qd/fx2Bgd5V/r0+YlmLY0l1X6dGLEc9V2rFctTrypLa+jds6Oe9Jucc6ve8W2IxYRcXF/PSSy/x\n8ssvExAQQLdu3UhMTGTSpEmkpKQQHh5Ou3btSE1NxWQymbd5enpSWFhIfn4+aWlphIRYn6YyO/uP\n+YizsowW9qwoK8vIZQuDIeojlrU4VQkM9K7W69SK5ajnSq1YjnpdVaU2/w0b+nmv7jmH+j/vlhK6\nxYS9efNmjh07xpIlSwCYPn06vXv3JjY2lpYtWzJ+/HhcXV0ZPXo0sbGx+Pr6snTpUgCmTZvGhAkT\ncHd3Z/HixTYdjBBCiKpZTNhRUVFERUVV2BYeHs7kyZMrbIuOjiY6OrrCtoiICCIiIlRqphBCCBk4\nI4QQTkKGpgtRz2xda1LWmRSSsIWoZ7asNSnrTAqQhC1EvSspMeEVoKOppXIvjUamfBWSsIWof9an\nfZUpXwVIwhai3tky7atM+SpAErbTs/WBFchDKyGcnSRsJ2fLAyuQh1ZCNASSsBuAwM634d3S1+I+\nub9fq6PWCOEcnLGcUhK2EKJRcsZySknYQohGy9q3U0f7ZioJ28mVlJjIs2FmsbzLuVLHK4STk4Tt\n9KzX8ILU8QpxM1tudhztRkcStpOzpYYXpI5XiFs534AlSdhCiEbJGQcsyfSqQgjhJKwm7EuXLjFy\n5EjCwsLMi+n26tWLsWPHMm7cOAwGAwDr168nJiaGqVOnmtdz3L17NzExMYwbN47MzMxaPAwhhGj4\nrCZsX19fPv74Y7p3727e1qlTJz799FM++eQT9Ho9RUVFJCUlkZiYyPDhw0lKSgIgPj6elStXEhcX\nR0JCQu0dhRBCNAJWE7abmxt6vb7CttOnTzNmzBjeeustANLT0wkNDUWr1RIZGcmhQ4coLCzEw8MD\nLy8vwsLCSE1NrZ0jEEKIRqJafdjffvstq1atIicnh+3bt5Obm4tOV7rMu06nw2AwYDAYzNsAc3eK\nEEKI6qlWlUjZHffgwYP55ZdfGDRokLnf2mg0otfr8fb2Nm8D0Gotfzb4+Xnh6lr6NDY7W2dx3/L8\n/XUWl4Wvj1jW4kDpPAbJyclWYw0ePNjiU+rGcK7UiiXnqn5iybmyL5YldiVsRVEoKCjAzc0NFxcX\nfvzxRzp37ky7du1ITU3FZDKRkpJCeHg4np6eFBYWkp+fT1paGiEhlsfiZ2f/MdlKVpbRwp4VZWUZ\nuWyh+L0+YlmLA3D6dCqz1i60Oo+Bj0+gxXkMGsO5UiuWnKv6iSXnyr5YlhK61YRdXFzMpEmTOHny\nJJMmTeKFF17g9ddfx8vLi+DgYJ5//nk0Gg2jR48mNjYWX19fli5dCsC0adOYMGEC7u7uLF682KYD\nakycbR4DIUT9spqwXV1d+eijjypsW7t27S37RUdHEx0dXWFbREQEERERNWuhEEIIQEY61htnnMdA\nCFG/JGHXG+ebx0AIW5SUlHD27Bnz7xkZ6eTnWB84l5+TSUZGusV9qorV1KcFGq3jDCGvLZKw64kz\nzmMghC3Onj3D35asx8unuXmbl0+Q1dd5+QQRv/kCcMHqfuXl51wCsDoB2s1ycq6Rc+Z7Luf8RnFh\nDi5uTXHzboFf+740bd6ZM9sWUlyQzVMp79CkiRu+vr507tyFhx4aTt++91SIdeLEz1w59iWZBdmY\nSm7g6q7Hw68NQXeOROvqYVe7LJGELYRQnZdPc7sTaF26cOF35sx5levXNQTe8SDu+tsAhfwraWQe\n/YoOg2YC4B38Jxa8PJ5WrYLJzLzI//63jVdfjWPEiEd44YWXAfj11zMsXboIt8Bu3Bb+KFoXN27k\nXcV48TiKqUTVdkvCFkI0Om+9tQitVkOz7jF4N2tn3u6ma453q57m3zUubuj1PjRvHkTz5kHceWd3\n2rZtz5IlC+jffyA9e/Zm//49eHt749H+Xty9WwDQxMufpoHqLysms/UJIRoVgyGHffv2MGjQELQu\nTW75u0sTy10YUVHReHvr+f777QAEBDQjNzeX6znnaqW95ckdthCiUTl37jcURaFly1bwe7Hdr9dq\ntQQHt+HChd8BGDBgMMnJW9i580uupW7BwzcYr4Db0bfuhYtb1QPjqkMSthCiUVFUKLwqnRtJA5Qm\n8EmTpnKqoCOa4jwKsjPIOv09V1O3Exw5DXdv6w9cbSVdIkKIRiU4OBiNRsP589XrwigpKeG33zJK\n79DLcXHXoW/dk6A7R9Duvjg0Gg3Zp79Xo8lmkrCFEI2KXu9Dnz53k5y8FVNJ0S1/LykqsPj6r79e\nR16ekQEDBlW5j0sTT1zcvTGV3Khxe8uTLhEhRKPz4ouvMGXKk1w5/DmaOx7CXV9a3ZF/5TRZp7+j\nw6BXATCV3ODatWtkZl7k0qVMvvsumS+//A+PPPIY3bv3AGDdui85eHA/hdd0uLl5oJiKMJw7yI3c\nTPw7DlC13ZKwhRCqKxvMUlfvVX6Qji1atmzF3LkLeXXRv7lyYhPFhQa0Tbxw17cgsMtw837G3/by\n/PP7aNKkCb6+fnTu3IUFC5ZWGDjTtWs39uz5gZzT35H9y0Y0Lk1w0wXSIvwx9K16qHacIAlbCKGy\ndu06sOylP5JeRkY6byUdsjraMT8nk+mPhdOmTdsq96kslpdPc5r6tLC7nb6+vvh0uK/KAT4dBs3E\nmH2ehVPutjjFcUhIJyZPnsYZZU+tDxaShC2EUJWLi8stCc7L54JNyaxNm7YWk6M9sRoieegohBBO\nQhK2EEI4CUnYQgjhJKwm7EuXLjFy5EjCwsLMK5+vWLGC2NhY4uLiKC4uHdq5fv16YmJimDp1qnnx\n3d27dxMTE8O4cePIzLQ+H64QQoiqWU3Yvr6+fPzxx3Tv3h2Aq1evsm/fPhITE+nUqRPJyckUFRWR\nlJREYmIiw4cPJykpCYD4+HhWrlxJXFwcCQkJtXskQgjRwFlN2G5ubuj1eqB01fRjx47Rp08fACIj\nIzl06BAZGRmEhoai1WrN2woLC/Hw8MDLy4uwsDBSU1Nr90iEEKKBs7sPOzc3F51OB4BOp8NgMGAw\nGCxuA8zdKUIIIarHrjpsjUaDt7c3Fy9eBMBoNKLX6/H29jb3W1e2DUpntLLEz88LV9fS5bCys3UW\n9y3P319HYKB3lX+vj1jW4qgZS86V+m1y1Fh1ea7UjCXnyr5YltiVsBVFoVu3biQmJjJp0iRSUlII\nDw+nXbt2pKamYjKZzNs8PT0pLCwkPz+ftLQ0QkIsF8NnZ+ebf87KMlrYs6KsLCOXLaw+Xh+xrMVR\nM5acK/Xb5Kix6vJc1SRWTRbhPXTouMX3tXURXmc5V5WxlNCtJuzi4mImTZrEyZMnmTRpEi+88AK9\ne/cmNjaWli1bMn78eFxdXRk9ejSxsbH4+vqydOlSAKZNm8aECRNwd3dn8eLFNh2QEMK5nT17hpfX\nz6FpucTj3w/glMXXeQAfnj8F5y3HvzlW3uVc8i4+aPPox/nzX2fz5o0AKGho4umLrkU3AjoNQevi\nhvHiMbLS/seNvMsoJcXMOPMFvXv34ZVXXrMpfm2ymrBdXV356KOPKmwLCwtj8uTJFbZFR0cTHR1d\nYVtERAQRERE1b6UQwqk0DfTGu6Vvnb1f4UXb99VoNNx1158YO/Yp5n28H01JPplH1mAqKcL7tm78\n/uMqmnUagq5FN/INl4jp35y0NMcompC5RIQQjYqiKLi6NkGv98HFXYfOrxMFV89gvHgcjdYFT79g\n87SoN4puEB7ek0ceeayeW11KRjoKIRodjeam37WuKKZiXN29uWG8zHXDhfppmBVyhy2EaHTKr+tY\nkJ2B4fxPNA0Mwbd9XwqyfiV9x9u4evjg2rQZ331npGXL1nh6etZfg/8fSdhCiEZn796U0lVibhSB\noqAL6krzbiPQurjRqs8EbuRdpeDqaXIzf+GLL1azefMm3n//Y/z8/Ou13ZKwhRCNTnh4L2JiYlma\ndAR9UCgaTcXeYbemAbg1DcDFuxUvjrydGTNe5Kuv1jBhwpR6anEp6cMWQjQ67u7uNG8ehIu79y3J\n+mbNmjXD3d2dwkLLi/PWBbnDFkKI/+fKya0opmKaNu9EE08/ioyXWLEigcLCQvr27V/fzZOELYRQ\nX56V0Xxqv5eL9d3MNBrNLVUiZbwCbuda+m4uHvoPJddz0bg0ISukA4sW/Z3u3cNVaW9NSMIWQqiq\nXbsOvDn8DfPvtb0IrwvYtQjvq6/+fwCcPn3rYBivZrfj1ex28+/G7PO8YmUR3rokCVsIoSpZhLf2\nyENHIYRwEpKwhRDCSUjCFkIIJyEJWwghnIQkbCGEcBKSsIUQwklUK2GfO3eOyMhIxo4dy8SJEwFY\nsWIFsbGxxMXFUVxcDMD69euJiYlh6tSpFdZ3FEIIYb9q32H37duXTz/9lA8++ICrV6+yb98+EhMT\n6dSpE8nJyRQVFZGUlERiYiLDhw8nKSlJzXYLIUSjU+2BM3v37mXMmDHcf//9tG/fnj59+gAQGRnJ\n119/TUhICKGhoWi1WiIjI5k9e7ZqjRZCOK6aLMKbkZFucR9bF+FtqKqVsJs3b87WrVtp0qQJTz/9\nNHl5eQQEBACg0+kwGAwYDAZ0Ol2FbUKIhu/s2TOkvPAct3l5mbc9D3DN+mtLEnbzq5V9bo51IT+f\nr++IrZPRjwcPHuBvf5vGxo3J6PU+tf5+N6tWwnZzczP/fN9996HT6cjMLP3UMxqN6PV6vL29zf3W\nZdss8fPzwtW19BMyO1tnc1v8/XUWl4Wvj1jW4qgZS86V+m1y1Fh1ea5qEis7W8dtXl600VmOX5tq\n6/h8fb3+37am+PlVPGa12mVJtRJ2Xl4eTZs2BeDgwYOMHTuWDRs2MGnSJFJSUggPD6ddu3akpqZi\nMpnM2yzJzs43/5yVZfsDyqwsI5ctzAxWH7GsxVEzlpwr9dvkqLHq8lzVJJY956a2WDq+Z5+dQrt2\n7cnLy+PC3u/QaLT4BPeh2R0PotFoUEzFXD31LYbzhyguNDBhwueMGTOOfv368+STTwIQEREBwNCh\nUbz66v+n6nm3lNCrlbAPHDjAsmXLcHNz46677iIsLIzevXsTGxtLy5YtGT9+PK6urowePZrY2Fh8\nfX1ZunRpdd5KCCFUt3XrZvr2vYfAsMdwwUTmkTW4enjj1+FeLh5KoiDrLIFdoynBhckPtMHFxZXm\nzYOYN+9NXnvtZT777Av0ej3u7u512u5qJez+/fvTv3/FybwnT57M5MmTK2yLjo4mOjq6+q0TQoha\n0KxZM5544kmOv7cHnV8rbuRdJvvMTpoGdSH39yO0+tNEmgaGYsw+T2hoZ/MMgt7epXe/fn5+9dKH\nLQNnhBCNikajoWvXOyts8/RrQ3GhgcLsdNBo8Aq4vYpX1y9J2EKIRkdRlPpuQrVIwhZCNCqKovDz\nz8cqbCvIzsDVQ4+HbxtQFPKvpFX62iZNmgBQUmKq9XZWRhK2EKLRuXLlCqtWfUxxfja5vx8h+8wO\n/Drcg5suEO+WYWQeWUPuhaMUF+Zw8uQJtmzZBECLFreh0WhISdlJdnY2BQV1u5K6LBEmhFDdhfx8\n6zvV03tpNBqGDBmK0Wjg8pEkNFotPsF34du+HwAtwh/jysmtXD6+nuLrRlb81ownnngKgMDA5kyY\nMIX33lvO4sXzeOCBh8xrRNYFSdhCCFW1a9cB/vGO+ffaXoQX7FuEF8DV1YWxY8fzc8Edt4yQ1Ghd\nCbzjQQLveBBj9nkW3rQI71NPTeKppybZ9X5qkYQthFCVoy/CqyiKPHQUQghnoNFo0Gg09d2MapE7\nbCFEo/KQDY0dAAAd20lEQVTPfyYAcPp0aj23xH5yhy2EEE5CErYQQjgJSdhCCOEkJGELIYSTkIQt\nhBBOQhK2EEI4CUnYQgjhJCRhCyGEk6j1gTMLFizg+PHjdOnShVmzZtX22wkhRINVq3fYx48fp6Cg\ngFWrVlFUVMTRo0dr8+2EEKJBq9WEffjwYfr27QtAZGQkhw4dqs23E0KIBq1Wu0Ryc3MJDg4GShev\nTE2teux+r17dzD8XFRWRZchHo3UhYvT/Vbr/7i9mo5hKGPmNl3kVCIAffzx2y775OZfY/cXsSuOU\nxc/PuVRle8q3q03klCrbA9zSpsraAzByZJT5GCtrz81tqqw9AP/5z1e3tL18e8qUtevIkZOVxunV\nq1uF815Ze8q3a+TIqArnvUz54y3frsrOv2IqgSnbq2xPGVuuh5SkV2+5Fm5uT/k2WboeKjufN5//\nsjZFPrag0jj2Xg8pSa/eci2Utad8u6tqTxlbrofybaqqPbZeD2XvVVV7yuLf3KbqXg83t6km12f5\ndlV1PXQf8leL7SlTVFTElWwDd9z7VKX7n9j5CYpiYuQ3eqv5yhKNUovzDK5atQp/f3+GDh3K1q1b\nyczMZOzYsbX1dkII0aDVapdIjx492L17NwC7d+8mPDy8Nt9OCCEatFpN2F26dMHd3Z0xY8bg4uLC\nnXfeaf1FQgghKlWrXSJCCCHUIwNnhBDCSUjCFkIIJyEJWwghnIQkbCGEcBIur7/++uv13Qh7nD59\nGn9/f7tfV1RUhItLxYEJWVlZeHp62hVnx44dfPDBByQlJbFhwwb27duHi4sLbdu2tTnGjRs3Kqzc\nvHv3brZu3Upubi7t2rWzqz0AJSUlaLWln725ubmkpqbi5eWFm5ub3bFutm3bNjp06FDjOGfPnmXP\nnj0A1fr3K+/q1at4eXnVuE0HDhxg8+bN5OXl0aZNG7teO2LECLKysggKCsLHx6fGbWnoDhw4QH5+\nPv7+/uzdu5ctW7Zw/fp188A6Z7djxw727duHj48Per3evH3NmjV06dJFtfdxuiqR++67j//97382\n779nzx5efvllrl+/TteuXZk7d675IhkxYgTr1q2zOda8efNIT09nxIgRBAUFAXDx4kX++9//0rZt\nW1577TWb4gwbNozPPvsMHx8fVqxYQXJyMvfeey/79++na9euxMXF2dymTZs2MXfuXLy9vZkxYwYL\nFiwgODiY9PR03njjDe69916bY23dupWyy0Gj0aAoCnPnzqXsM33IkCE2xxo7dizLli3D39+fdevW\nER8fT+/evTly5AijR49m3LhxNsW5du1ahd8VReGRRx5h7dq1APj6+trcplGjRrFmzRoA/vOf/7Bq\n1Sruv/9+du3axYABA/jLX/5ic6yBAwfy5z//mW+++YaAgACGDRvG0KFDzdeFPYxGIx988AFbt27l\nwoULNGnShDZt2vD444/z8MMP2xzHYDDw3nvvkZyczNWrV9FoNPj7+zN48GCmTJlSIZFYM2LECIYM\nGUJUVJTdH2Y3e+utt9i7dy8lJSX86U9/4sCBA/Tv35+UlBQGDBjApEmTbI6l5jH++uuvJCQk4OPj\nw1NPPcXs2bM5cOAAbdq0Yd68eYSFhdl8fAcPHqRLly589913jBs3znx925tjrFIc0BtvvFHlf+Hh\n4XbFGjlypHLq1CnFZDIp33zzjTJ48GDl4MGDiqIoSnR0tF2x7r///kq3m0wmZfDgwTbHeeihhyq0\nr6CgQFEURSkqKlKioqLsalNUVJRy6dIlJSMjQwkLC1PS0tIURVGUc+fOKY888ohdse644w5lypQp\nyowZM5QZM2Yor7zyihIeHm7+3R43H2NWVpaiKIqSn59v1zF26tRJGTBgQIX/unTpogwYMEAZOHCg\nXW0q/+89cuRI5erVq4qiKEpeXl6F9toTy2QyKfv27VPmzJmjREZGKk888YSyevVqu2L95S9/Udas\nWaP8/vvvysqVK5V//vOfypkzZ5SXXnpJeeutt2yOM378eCUhIUG5dOmSYjKZFEVRlMzMTCUhIUEZ\nP368XW0aMGCAsmjRIqV///7Kww8/rHz44YfKxYsX7YpRZujQoUpRUZGSn5+vhIeHKwaDQVEURSko\nKLD7elfzGB977DFl9erVyvvvv6/069dP2bhxo1JQUKDs2rVLGTVqlM1xHnroIeXGjRuKoihKTk6O\nMnHiRGXevHmKyWSyO8dY45B92GvXriUkJISuXbvSrVs3839du3a1+2t+UVERISEhaDQaHnjgAZYv\nX86MGTNITk62u13u7u4cOXLklu1HjhzBw8PD5jhNmzbl5MnSOT38/PwoLCwEoLi42HyHaysXFxcC\nAwMJDg6mdevW3H777QC0atWK4uJiu2KtXr2agoICunXrxoIFC1i0aBF+fn4sXLiQhQsX2hXL1dWV\nixcvAqXHW9b15ObmhslksjnOyy+/TPv27Vm+fDnbt29n+/btBAUFsX37drZt22ZXm0pKSrh27RrZ\n2dmYTCZz14yXl9ct3WW20mg03HXXXcydO5fvv/+eyZMn2z3J2fnz53nkkUe47bbbGD9+PNu3b6d9\n+/YsWrSILVu22BVnypQpBAYGmrvbmjdvzpQpUzh//rxdbdLr9bzyyit89913zJgxg19//ZWHH36Y\nsWPHkpSUZFesJk2a4OrqiqenJ23atMHb2xsADw8Pc1eerdQ8xoKCAh577DEmTZqEh4cHDz74IB4e\nHvTt25cbN27YHKekpMQ8P4her+ff//43RqORv/3tbxQVFdnVJmtqfT7s6ujWrRshISH06tXrlr+9\n++67dsVq0qQJly9fJjAwEICQkBA+/vhjpkyZwm+//WZXrIULF/L666+Tl5dHixYtgNIuEZ1OZ1dC\nmzt3Li+99BKdOnUiICCARx55hLvuuotTp04xZUrlk0tZYjKZ0Gq1FdpQXFxsd8IOCwvjww8/5NNP\nP2XcuHF2dc3cbObMmUycOJEhQ4YQEhLCU089Rd++ffnxxx/t+po/YcIEhg4dysKFC2nRogXPPfdc\ntduUl5dnfm+NRsOlS5do3rw5RqPR7liVPWtwdXXl3nvvtasbCko/MA4cOEDv3r1JTk7Gz88PwO5k\n1rJlS95//31GjhxJs2bNALh8+TJfffUVt912m12xypR9IN11113Mnj2blJQUvvnmGx577DGbY7i5\nuVFQUICnpydfffWVebvBYKjXYyxL+FB6U1HV36wJDg5m37599OnTByi9DhYuXMg//vEPtm7daleb\nrHHIPuxr167h7u5u9wPByvzwww/4+/tzxx13VNhuMBj47LPPePrpp+2OeenSJTIzMwEICgqiefPm\ndscoLi7mhx9+4OzZsxQXF3PbbbfRt29fux9gHTlyhNDQ0Fvu8M+dO8ePP/5IdHS03W0DyMzMZP78\n+Rw/ftzuO9kyBoOBr7/+mvT0dIqLi2nRogWDBg0yfwuwV3JyMgkJCZw/f56UlJRqxahMQUEBV65c\nqbcHYCdOnGDWrFmkp6fTsWNHFixYQIcOHcjKymLDhg129fe/9957bN++natXrwIQEBDAwIEDmTJl\nil39/c8//zxvv/12tY7nZtevX8fd3f2W7VlZWVy6dInOnTvbHEvNYwwLCzMXC2RkZFToq8/IyODw\n4cM2xSn7hlzZt+yLFy+ab+7U4JAJu8z169dJT0+npKSEtm3b1qgyQM1YZYxGI+np6QQHB9v1sKO2\nXb16lYCAgPpuRq0oLCwkIyOD0NBQu1975MiRKh8krVu3jhEjRtgc6/Dhw3To0AFvb28KCgp47733\n+Pnnn+nYsSNTp041f+0Xtzp79iwnTpygY8eOdOzYsd7ace7cOYt/b926tV3xioqKbpnSNSsrq8ZV\nUeU5ZMIuKiriH//4B19++SUtW7YE4MKFC0RFRfHKK6+QkZFh812amrFef/11c8XEgQMHiIuLq1CR\ncd9999kUZ8eOHeavzQaDgUWLFnH06FFCQ0OZOXOm+aueLdSsopg5c6bFv9vT7aNWrEOHDnH77ber\nkhiHDRtGz549mT59uvkD9uTJk7zxxhvo9Xri4+NtjvXggw+yfv16XF1dee211/D09OTPf/4zu3fv\n5uTJk3Z33ZV34MAB8zenfv362fy6Tz75hPvvv7/a3R/lqfmBpFbFEKh7PdwsKyuLAwcO0LJlS7p1\nq3x+78qoWYlmjUMm7Pnz55Ofn8/MmTPR6XRA6d3sokWLMBqNpKWlsWHDhjqPVf7kjx07lhkzZtC1\na1d+++03/va3v5mTpD1xXn31VQIDAxk9ejTJycns27eP5cuX2xQHoHPnzuYPojKZmZnm8rLt2yuf\nEL4ymzdvNpfzlfXhXbhwgY8++oiSkhJ27txZrVhQ2idYnVhqJsaioiJWrlxJUlIS06ZN49SpU+zY\nsYMZM2YwYMAAm+MADB06lG+++QaAkSNHVuibHT58OOvXr7c5llrlhr169cLT05Pg4GCioqIYOnRo\nte/u1DzvUVFR5v/HHn74YT744AP8/PwoKCjg0Ucf5euvv66Xdk2ZMoW4uDhCQ0O5dOkSI0eO5M47\n7yQjI4NHH32Up556yqY4Dz/8MIsXL6Zjx45s2bKFt956izfffJMePXo0jrK+wYMHKyUlJbdsLy4u\nVnr37m0uy6vrWOVLdG4u17GnfKf8vsOGDTOXJ5X9bo8PPvhAmTBhgvLLL7+Ytw0YMMCuGJVJT09X\nXn31VWXIkCHKqlWrlOvXr9dLrAceeMD884gRIyr8zd5zVea9995TOnXqpNxzzz3VLlX761//qnzx\nxReKoijKjBkzlMOHDyuKoihnzpxRHn74YbtiqVVuGB0drZSUlCg7d+5UZsyYofzpT39SJkyYoKxd\nu1bJzc21q01qnvfo6GjlwoULiqIoyhNPPGEuYy0uLlYefPDBemtX+feOj49XXnrpJUVRFCU3N9eu\ncsOb9z116pQyZMgQ5dtvv1W9rM8hq0S0Wm2lT49dXFzw9/enR48e9RLrzJkzDBs2DCjt/8rJycHH\nx4eSkhK7KjKysrL48MMPURTllgoFxc4vPGpWUUDpSNL4+Hh+/vlnJk6cyNy5c3F1rd5lokaskJAQ\n1qxZw6hRo+jcubO5H/rXX3+tdAkoS8q6rqB0wNGOHTsYM2YMU6dOZdSoUXbFmjdvHvPnzyc+Ph5/\nf38ef/xxgoKCuO2225g/f75dscrKDRVFqXG5oVarpV+/fvTr148bN26wY8cONm7cyKJFi9i7d6/N\ncdQ872pVDKndrvLXYkpKCo8++igAOp3OrioRNSvRrHHILpFp06YxZMgQRo4cWWH7unXr2LJli119\njWrGuvkhRfPmzXFzczP3fdk6EvCf//xnhQvi8ccfJyAggEuXLrF06VLefPNNm9tUXk2rKJ577jmO\nHz/OhAkTeOCBB9BqtRXaaU9/uFqxDAYD8+fP58CBA/j7+/Pzzz+bE+Ps2bPtqjC4//77efHFFxk6\ndKh5W2ZmJgsXLuTixYusXr3a5lhlcnNzOXfunLkKpux/WnsMHDjQ/LNGo+Hzzz83lxuOGTOG//73\nvzbFsfT1Oz8/364H7Wqe97J4alQMqdmuv/zlL/Tr14+goCBmzZpFcnIyPj4+FBQUMGrUKDZu3GhT\nnKoq0XJzc/nss8+YNm2aXcdoiUMm7IsXL/Lss8/i4eFB165dATh+/DgFBQX861//sqtMRs1YeXl5\nt9Rrlrm5LKi+FBQU8Ntvv1WriqJ84qiMPf3hasYCdRKj0Wg0P8e42Q8//EDfvn3tjlmmoKCAtLQ0\nWrVqpVpVgL3lhmfOnFFl3pfy1DjvtUGNdl25coV33nmHy5cvM2bMGPMD3j179nD8+HEmTpxoU5zk\n5GQuXrzIE088AZQ+k8jKykKj0RAXF1fhBqGmHDJhQ2nXwJ49e0hNTUWj0dCxY0ciIiLqNdbgwYN5\n8cUXefDBB83bCgsLiY+PZ+PGjTaPnlQUhf3796PX6+ncuTObNm1i//79tG3bltjYWLtGcyqKYp50\npqaxHJGapXjHjx+3+PeyD3RbbNu2jfnz5+Pj48Pzzz/P3LlzadasGefOnSMuLs7ur/qVycvLY9Wq\nVdUaTAWld9WnT59W9UOkOm1Ss/qoKjU9VzURExPD3//+d/PD/+joaD766CMKCgqYMWMGn3zyiWrv\n5ZB92PPmzaNnz5706tWr2km6NmKtXLmSuXPnsmbNGubMmUNaWhpvvvkmgwYNsvlrK5SOdDx16hTX\nr1+nffv25Ofnc8899/Djjz8yc+ZM3nrrrXqJpWZCUyvWrFmzLJbi2ZOwFy1aZP752LFjt5Ruffrp\npzbHevvtt/nggw/Izc1l3LhxfP311wQHB3P16lWefPJJuxJ2ZmYm7733HhkZGYSEhPDMM8/wxRdf\nsHLlSrsm3FLzQ0StNgH079/fYsWQPdRs1zPPPMOQIUMYPHhwld+cbXHjxo0KlVq9evXCz8/PXAmj\nJoe8w/7000/56aefOHToEIqi0LNnT3r06EHPnj3p3LmzXcNZ1YxV5v333+cf//gHzZo1Y8WKFXZ3\nPwwdOpRNmzZx/fp17r33XlJSUnB1dUVRFIYNG2ZzmaHascaOHWv+uaYJTa1YapbilVfTcqvyrx82\nbFiF0jR7Yz/11FP06tWL7t27s3PnTrZt20b37t3NJZ+2GjZsGO+8806VHyL2XAtqtelmGRkZJCQk\ncODAAZ588klGjRpl17dANdt1zz330KNHD/bs2UNERARRUVH079/f7m+lgwcPrvLb9aBBg6o9UrhS\nqtac1IKLFy8q33zzjbJgwQJl0KBBSo8ePeot1o0bN5R///vfysCBA5XPP/9cmTZtmjJu3DjzDHm2\nUqs8UO1Yar22NmKpUYqnZpuioqKU7OxsJSsry/xz+d/tMXz48Aq/33PPPUpxcbHdbSp/TDe3wd7j\nVatNZdLS0pTp06crQ4cOVdasWaMUFRVVK46a7So7J7m5ucpXX32lTJw4UenTp48yY8YMZefOnTbH\nefHFFyudoTExMVF54YUXqtW2qjhklwiU9s2ePHmSgwcP8tNPP3H69Gnatm1r11dgtWONHDmSu+66\ni3Xr1uHt7U1MTAzfffcdTz/9NEOGDGH69Ok2xSlf1lf+57K/2UPNWI5IzVI8NZWfSAqoUZ+1yWQy\nj1hVFAUfHx9yc3PNf7e1oqZ8eSBQIaa9XQ9qtQkqVgzNnDkTrVZboZzVnlhqtquMTqdjxIgR5kUp\ntmzZQkJCgs2jTGfOnMkzzzzDhg0bzIsV/Pzzz1y/fp1//etfdrfHEofsEhk/fjxGo5E77riDsLAw\nevToQYcOHeyqjayNWEePHuXOO++8ZXvZg8cXXnjBpjjly/qUcqMKy35+9tlnbW6TmrHKU3OEVk1i\nqVmK93//93/mnzdt2sRDDz1UoV/V1gUo1KZWRY2jVvk4aqwxY8awatUqm/e3RFGxSMISh0zYc+bM\n4cSJE3h4eNC9e3d69OhBeHh4tZ50qxmrzG+//UZaWholJSWEhITYtTyYNZaqImo7lpoJTa1Yapbi\nrV271uKH2821+pao+YBWCFs5ZMIuk5uby6FDh8z/ZWdn07Fjx2oNLFEjltFoZNasWRw7dsxcJP/L\nL78QGhrKkiVLOHHiBL1797a7bampqWzYsIFNmzbh7e1t85wkasdSM6GpFctRE6OaD2jVOkZHrPJx\n5FhQOhBn586d5umSW7RoQb9+/Rxq9s3yHDph37hxgyNHjnDw4EEOHjzI4cOHCQgIsOtpt5qxXnnl\nFVq3bs0zzzxjri4xmUwsX76cvXv3cu3aNZsnsvntt9/YtGkTGzZsoEmTJpw/f54vv/zS7ikd1Y7l\naNRMjGqVcd2spt1Hah2jI1b5OHKsdevW8e677xIZGVlhQZIffviBZ5991q4blDqj6iNMlcyfP18Z\nNWqU0rt3b2Xs2LHK3//+d+V///ufkpOTU6+xLK3bePfddytnzpyxKc6jjz6qjBw5UklISFAyMjIU\nRan+hE1qxnr66aeVdevWKUajsVqvr61YZWpa2dGvXz/lr3/9q3LXXXcpzz33nLJ169YaTWqlVrtq\nI5YjtsnRYg0ZMqTSPHDt2rUq12+tbw5ZJdKqVSuGDx9O586dqz3xUG3EsvSgUqfT0b59e5viBAQE\nkJqaypUrV7h69WqNVjpRM9aRI0dwcXFh/vz5NapLVTuWWgICAnjnnXcwGo0kJyeTlJTEa6+9xsCB\nA3nooYfsmntaNFzVKUioKw6ZsHv16oVGozEvVHsze/qp1IwVHh7Ou+++yzPPPFOhT3b58uV2zfq3\nfPlyDAYDW7du5Z133iE9PR2DwcDhw4fp3r27zXHUjqVmQnPk5FjTMi6o+FA1MzOTefPmOUTFibDd\n1KlTefjhh4mMjDQv+nDhwgV++OGHai0dWBccsg/bUfu8cnNzmTVrFsePH6/w0LFLly4sWLCg2qtd\nXLlyhW+++YaNGzdy4cIFvv/++2rFqWmsyvpiyxLapk2b7DpXasVSs3JFzTIuNR/QqnWMjljl48ix\nAHJycti1a1eFNVr79etn99qqdcUhE3Z5jlIPXF56ejppaWnmeks1Z+k7d+6cag8L7Y2lZkJTK5aa\nidFRqXWMjljl48ixLH0rdnNzo23btjz//PNERkbaHLO2ScK2w7FjxypcIFCxv8vW7hU1qxVqq/Kh\noVKrjEvOe8NWXFxMWloa06dPt3le7Lrg8nrZqrIOavXq1cTExDhErOnTp7Nr1y527dpFQkICqamp\n5t937dpl8/DkRYsWUVBQwPz58zl27Biurq4EBwfbvbqI2rGgNKElJyeza9cufvrpJy5cuEBQUBDu\n7u71EuuZZ56huLiY4ODgGj+wXLduHdOnT0er1eLh4UFxcTEnTpxgyZIl6PX6Wyagt0TN867WMap5\nrhpDLEu0Wi3NmjVDq9VWOrq5vjjkHbYj93mVqcndetlryx7IbdiwgaNHj1brgZyasdSsS1Urlloz\nqgH8+c9/5osvvrjlbjonJ4fRo0ezdetWm2Oped7VOkY1z1VjiOWMHLJKpGvXruauhvI/l++zqo9Y\nalOjWkHNWPHx8axdu7bKhGZPwlYrVl1Um9TkOlDjvKt1jI5a5eOosZxSLdd5N1g1KdqPjY1VrR1q\nxqpqIEFOTo7dAwnUilXZeb569aqSmJioPPHEE3a1ae3atcqgQYOU2bNnK8uXL1eWL1+uzJ49Wxk4\ncKCyZs0au2Kped7VOkY1z1VjiOWMHLJLxFEfyjnqbG9q+eqrr/jXv/5VZV3qI488Uuex1KxcAccs\n41LrGB2xyseRYzkjh0zYjtrnpWZJkZqTzqgZS82E5mjJUe0yLmebOEg4P4dM2I76UE4tjvhwD9RN\naGrGqovEaG8Zl9oTB6l1jI56I+CosZyNQyfs8up7xB2o172iZrWCmrEsUbMu1Z5YdT2j2ueff87j\njz9udT81z7tax+ioNwKOGssp1VPfuUWO+lBOrdneHPHhnq0SExPrNJajzqhWF/+G9h6jmueqMcRy\nRg5Z1qfmQwU1Y6lVUqTmpDN1PYGNLXefdRGrvksy6+K8q3WMap6rxhDLkTlklwg4Zp+Xmt0rDfnh\nnprUrFxRm1rnXa1jLIvTt29fc3dBTat8GnIsZ+SQCdtR+7zUKily1Id7jsoRP5DUPu9qHeO1a9fY\ntWsXly5dMscJCgpiw4YN2DsLRWOI5WwcMmE740M5tdTXwz1H5YwfSPae99o4xuPHj7NhwwY2b95M\n69atGTJkSIWphu3RGGI5C4fsw66KI/R51XZJkaurK507d+aJJ55wqFj15aeffqryb476gWTveVfr\nGM+cOcPGjRvZtGkT/v7+PPDAAyiKYldXXWOK5YwcMmE76kO5yrpX9uzZw1tvvaV6SZGjPNxzZI7+\ngaTGebfnGB988EHuu+8+PvjgA1q2bAnAhx9+WK33bQyxnJFDdomAYz6Uc7buFdG4lA0MO3LkCP36\n9WPo0KHMmjWL7du3S6wGwiETtqM+lKsqYRsMBkaNGiUJWziEvLw8tm3bxsaNG9m7dy/R0dHcf//9\n1RrV2xhiOROHTNiW1OdDOUcuMROiMteuXWPLli1s3LiRTz75RGI5OadL2GVsHUKsdixHLDETQjQO\nTpuw64MzlpgJIRoOSdgqcdQSMyFEw6Gt7wY0FI5eYiaEcH5yhy2EEE5C7rCFEMJJSMIWQggnIQlb\nCCGchCRsIYRwEpKwhRDCSfz/xS3yGT5WJuAAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x5252e90>"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Print a list of the genomes, adding the product\n",
      "annotation_files = \"/home/ubuntu/LT_data/CompGenomics/LT_Genomes/v2/processed_genomes/annotation/\"\n",
      "\n",
      "product_annotation = defaultdict(str)\n",
      "#Save product annotation\n",
      "for genome in genomes_order:\n",
      "    for line in open(annotation_files + genome + \".txt\", 'r'):\n",
      "        if line.strip():\n",
      "            line = line.rstrip()\n",
      "            elements = line.split(\"\\t\")\n",
      "            \n",
      "            if elements[1] == \"Product\":\n",
      "                product_annotation[elements[0]] += elements[2]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "outfile_jan = open(\"jan_genes_selection_unique.txt\",'w')\n",
      "outfile_aug = open(\"aug_genes_selection_unique.txt\",'w')\n",
      "for genome in genomes_order:\n",
      "    jan = set(jan_genes[genome])\n",
      "    aug = set(aug_genes[genome])\n",
      "    print genome, len(jan.difference(aug)), len(aug.difference(jan))\n",
      "    \n",
      "    for entry in jan.difference(aug):\n",
      "        dn_ds = pN_pS_Jan[genome][entry][2]\n",
      "        outfile_jan.write(genome + \"\\t\" + entry + \"\\t\" + product_annotation[entry] + \"\\t\" + str(dn_ds) + \"\\n\")\n",
      "    \n",
      "    for entry in aug.difference(jan):\n",
      "        dn_ds = pN_pS_Aug[genome][entry][2]\n",
      "        outfile_aug.write(genome + \"\\t\" + entry + \"\\t\" + product_annotation[entry] + \"\\t\" + str(dn_ds) + \"\\n\")\n",
      "\n",
      "outfile_jan.close()\n",
      "outfile_aug.close()\n",
      "    \n",
      "\n",
      "    \n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "J07HWQ1 0 1\n",
        "J07HWQ2 1 1\n",
        "J07HQX50 0 0\n",
        "J07AB56 37 20\n",
        "J07AB43 25 16\n",
        "J07HN4 3 1\n",
        "J07HN6 3 8\n",
        "J07HX64 3 1\n",
        "J07HX5 1 1\n",
        "J07HB67 54 26\n",
        "J07HR59 0 6\n",
        "A07HB70 1 3\n",
        "A07HR67 0 4\n",
        "A07HN63 4 12\n",
        "A07HR60 2 6\n",
        "J07SB 2 2\n",
        "G22 0 0\n"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print jan_genes[\"G22\"]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[]\n"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print aug_genes[\"G22\"]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[]\n"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}