Skip to content

Instantly share code, notes, and snippets.

@SilasK

SilasK/Atlas_gene_graph.ipynb

Created October 21, 2021 13:02
Show Gist options
  • Save SilasK/5932a6887ee4d2520b5a59cec06d09b7 to your computer and use it in GitHub Desktop.
Save SilasK/5932a6887ee4d2520b5a59cec06d09b7 to your computer and use it in GitHub Desktop.
Download ZIP
Linking Atlas genes to contigs and bins
Display the source blob
Display the rendered blob
Raw
Atlas_gene_graph.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 5,
"id": "opposite-model",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import pandas as pd\n",
"import numpy as np\n",
"import matplotlib.pylab as plt\n"
]
},
{
"cell_type": "code",
"execution_count": 303,
"id": "dominant-monitoring",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"check if file contigs2bins exists: True\n",
"check if file contigs2mags exists: True\n",
"check if file allbins2genome exists: True\n",
"check if file old2newID exists: True\n",
"check if file orf2gene exists: True\n",
"check if file taxonomy exists: True\n",
"check if file egg_Nog exists: True\n"
]
}
],
"source": [
"\n",
"atlas_dir= \"/Users/silas/Documents/scripts_baobab/scratch/Debug_atlas/WD/\"\n",
"\n",
"\n",
"\n",
"\n",
"files=dict(\n",
"contigs2bins=\"genomes/clustering/all_contigs2bins.tsv.gz\",\n",
"contigs2mags=\"genomes/clustering/contig2genome.tsv\",\n",
"allbins2genome=\"genomes/clustering/allbins2genome.tsv\",\n",
"old2newID=\"genomes/clustering/old2newID.tsv\",\n",
"orf2gene=\"Genecatalog/clustering/orf2gene.tsv.gz\",\n",
"taxonomy=\"genomes/taxonomy/gtdb_taxonomy.tsv\",\n",
"egg_Nog= \"Genecatalog/annotations/eggNog.tsv.gz\"\n",
")\n",
"\n",
"for key in files:\n",
" files[key]= os.path.join(atlas_dir,files[key])\n",
"\n",
"# check if files exists\n",
"for key in files:\n",
" print(f\"check if file {key} exists: {os.path.exists(files[key])}\")"
]
},
{
"cell_type": "code",
"execution_count": 112,
"id": "frank-sword",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>ORF</th>\n",
" <th>Gene</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>sample1_1_30</td>\n",
" <td>Gene0001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>sample2_172_4</td>\n",
" <td>Gene0001</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>sample1_2_33</td>\n",
" <td>Gene0002</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>sample2_72_9</td>\n",
" <td>Gene0002</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>sample1_3_8</td>\n",
" <td>Gene0003</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" ORF Gene\n",
"0 sample1_1_30 Gene0001\n",
"1 sample2_172_4 Gene0001\n",
"2 sample1_2_33 Gene0002\n",
"3 sample2_72_9 Gene0002\n",
"4 sample1_3_8 Gene0003"
]
},
"execution_count": 112,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"orfs = pd.read_table(files['orf2gene'])\n",
"orfs.head()"
]
},
{
"cell_type": "code",
"execution_count": 113,
"id": "regulated-description",
"metadata": {},
"outputs": [],
"source": [
"orfs['GeneId']= orfs.Gene.str[4:].astype(int)\n",
"orfs[['Contig','GenePosition']]= orfs.ORF.str.rsplit('_',n=1,expand=True)\n",
"orfs['Contig']= orfs.Contig.astype('category')\n",
"orfs['GenePosition']= orfs.GenePosition.astype(np.uint16)\n",
" \n",
" \n",
"orfs.drop(['ORF','Gene'],axis=1,inplace=True)\n"
]
},
{
"cell_type": "code",
"execution_count": 114,
"id": "alone-employer",
"metadata": {},
"outputs": [],
"source": [
"orfs.index= pd.MultiIndex.from_frame(orfs[['Contig','GenePosition']])\n",
"#orfs.sort_index(inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "overhead-modern",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 181,
"id": "proprietary-pendant",
"metadata": {},
"outputs": [],
"source": [
"right_neighbour_index= pd.MultiIndex.from_arrays([orfs.Contig, orfs.GenePosition + 1])\n",
"valid_neighbour= right_neighbour_index.isin(orfs.index)\n",
"right_neighbour_index= right_neighbour_index[valid_neighbour]\n",
"left_neighbour_index= orfs.index[valid_neighbour]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "presidential-hollywood",
"metadata": {
"tags": []
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 184,
"id": "skilled-principal",
"metadata": {},
"outputs": [],
"source": [
"edgelist= pd.DataFrame({'LeftGene': orfs.loc[left_neighbour_index,'GeneId'].values,\n",
" 'RightGene' : orfs.loc[right_neighbour_index,'GeneId'].values,\n",
" 'Contig': orfs.loc[left_neighbour_index,'Contig'].values\n",
" })\n",
"\n",
"support= edgelist.iloc[:,[0,1]].value_counts()\n",
"support.name='N_connections'"
]
},
{
"cell_type": "code",
"execution_count": 185,
"id": "protective-danger",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>LeftGene</th>\n",
" <th>RightGene</th>\n",
" <th>Contig</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>3940</td>\n",
" <td>sample1_1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>1</td>\n",
" <td>689</td>\n",
" <td>sample2_172</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2</td>\n",
" <td>3366</td>\n",
" <td>sample1_2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>3367</td>\n",
" <td>sample1_3</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>3</td>\n",
" <td>3367</td>\n",
" <td>sample2_36</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" LeftGene RightGene Contig\n",
"0 1 3940 sample1_1\n",
"1 1 689 sample2_172\n",
"2 2 3366 sample1_2\n",
"3 3 3367 sample1_3\n",
"4 3 3367 sample2_36"
]
},
"execution_count": 185,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"edgelist.head()"
]
},
{
"cell_type": "code",
"execution_count": 188,
"id": "resistant-punch",
"metadata": {},
"outputs": [],
"source": [
"Gene_of_interest='Gene3367'\n",
"GeneID_of_interest = int(Gene_of_interest[4:])\n",
"\n",
"Gene_Neighbourhood = edgelist.query(\"LeftGene == @GeneID_of_interest | RightGene == @GeneID_of_interest\")"
]
},
{
"cell_type": "code",
"execution_count": 189,
"id": "deadly-locking",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>LeftGene</th>\n",
" <th>RightGene</th>\n",
" <th>Contig</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>3</td>\n",
" <td>3367</td>\n",
" <td>sample1_3</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>3</td>\n",
" <td>3367</td>\n",
" <td>sample2_36</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3643</th>\n",
" <td>3367</td>\n",
" <td>2025</td>\n",
" <td>sample1_3</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3644</th>\n",
" <td>3367</td>\n",
" <td>2025</td>\n",
" <td>sample2_36</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" LeftGene RightGene Contig\n",
"3 3 3367 sample1_3\n",
"4 3 3367 sample2_36\n",
"3643 3367 2025 sample1_3\n",
"3644 3367 2025 sample2_36"
]
},
"execution_count": 189,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Gene_Neighbourhood"
]
},
{
"cell_type": "markdown",
"id": "million-enlargement",
"metadata": {},
"source": [
"## Graph Plot \n",
"Plot gene neighbourhood"
]
},
{
"cell_type": "code",
"execution_count": 267,
"id": "otherwise-warren",
"metadata": {},
"outputs": [],
"source": [
"import networkx as nx\n",
"\n",
"G = nx.from_pandas_edgelist(edgelist,'LeftGene','RightGene',edge_attr='Contig', create_using=nx.MultiGraph)"
]
},
{
"cell_type": "code",
"execution_count": 293,
"id": "atomic-attraction",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 299,
"id": "brutal-announcement",
"metadata": {},
"outputs": [],
"source": [
"neighborhood_size=5\n",
"\n",
"Gneighborhood= G.subgraph( nx.single_source_shortest_path_length(G, GeneID_of_interest, cutoff=neighborhood_size).keys() )"
]
},
{
"cell_type": "code",
"execution_count": 300,
"id": "likely-norfolk",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Legend contigs: {'red': 'sample1_3', 'blue': 'sample2_36'}\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAFMCAYAAAAUZHNQAAAAAXNSR0IArs4c6QAAQABJREFUeAHtXQncTdXaf6RooswiQ4ZkShQiRNH4ZWgSLqlug5sGuimXut24zQPVpa4UFU2KBtJkSlHJUJGoSJShSdHM+p7/Ou1zz3vec/ZZa5+999nnvM/z+73vGfYa/2vv56z1jKUUEwkJAoKAIBAcAuP3CK5taVkQEAQEgRgCwmjkThAEBIHAEdgz8B6kA0FAEAgFgV9//ZU+/fRT+uqrr/Tfd999R7/88ov+K126NO299976r3r16nTQQQdRrVq16OCDDw5lbKVERhMKztKJIOArArt27aJly5bRggUL6M0336QPP/yQPv/8c2rQoAFVrVpVM5ADDjhAM5b99tuPfvvtN81w8Lpp0ybNiPAKZtS8eXNq2bIldezYkTp16kQ1a9b0dazc2HhhNH5DKu0JAgEhAOYya9YsevbZZ+mFF17QzOSYY46hzp07U9OmTenQQw+lPfe0O6Rs376dVq5cSW+//Ta98cYbtHDhQt1ur1696JxzzqGGDRv6MRthNH6gKG0IAkEisGXLFho7diw9/PDD1KhRIzrjjDOoZ8+e+ugTRL/YJT3//PP06KOPUuPGjWnQoEF01lln0R57eBbpCqMJYqGkTUHADwS2bt1Ko0aNoilTptDAgQPpwgsv1A++H22btPH7779rhnPvvffqo9aIESOof//+VKpUKZPqiWVEvZ2IhrwXBKKAAI5Id911FzVp0oT23Xdf+uijj/Rn7C7CpL322kvvnubNm0f3338/TZgwgdq1a0fvvvuu9TDsDnTWzUsFQUAQsEHgk08+oXPPPZf2339/LeTFUSkK1KVLFy3DmTRpEnXv3p0GDx5M1157LUGbZUKeD10mjUsZQUAQMEdgzpw51KFDBzr99NPp5Zdf1vIY89rhlMQRbsWKFfTaa6/R8OHD6YcffjDqWBiNEUxSSBAIFgGoqI8//niaPHkyXXXVVcF2lmXrUJ+DKe7YsUMfpSDLyUTCaDIhJNcFgYARgIYHx5H33nuPTjzxxIB786d5CITHjRund2Dly5en77//3rVhsaNxhUcuCgLBIrBkyRI65ZRTCAJXCH/zkSCvOfDAA+mGG25IZ8cjWqd8XFgZc2Eg8PPPP2t18ZgxY/KWyWAl7rnnHi0ovvPOO9MujOxo0kIjFwSBYBEYPXo0ffDBB/Tkk08G21EIrW/cuJFatGihBcUp/KfEYC+ENZAuBIFiCHzzzTfaZQA+SnBw9EI7d+4k+DElE3ZKoH322Sf5kvZtgkoacpVUlK7NVGWTvwPjhFMnLJiTSBhNEiDyURAIBQEcN+BfBKtfW4Lg9fbbb9d1169fH68OQ78BAwZoy114clerVo3uu+8+fR11+vbtS+XKldO2Lz/99BNNnTpVGwSiQLo2440bvPn666+pfv369MUXXyQzMpHRGOAnRQQB3xF47LHH6OKLL/bU7iuvvEJgMFAvJxKEsXCQRNs4js2cOVNb9KIMrsGBEt+DuYEZwA7GoXRtOtdNXitXrkwnn3yydvpMLi/q7WRE5LMgEDACiBEDhtC6dWtPPZ199tnUpk2bYnXh1d2nTx/9PRwg4Qj53HPP6c8zZsyI14FqGq4EEydOpD/++ENfT9dmsU4yfAFv8kWLFhUrJYymGCTyhSAQLAIff/yxls+kkqF47RlyGbSbGEsGsh/HLwkyIUd2gz6OOOIIgjwGMWz8pKOOOkrHxkluUxhNMiLyWRAIGAE89BUqVPC1l/fff58go8HxxaEqVaoQ+vrss8+oW7duWibkXCtTpox+axu/xqmf7hXzQp/JJE6VyYjIZ0EgDxGA8BcEj2uHypYtq9/iqHTeeefpGDYINVGvXj269dZbtWC4du3aTvFAX4XRBAqvNC4IFEegUqVKWs1c/Ir3bxzbFcSwcQjCYsQJrlOnDh1yyCHaG3z+/PlaM4SQE9jxeIgt4zSf8hWhQTG/ZBJGk4yIfBYEAkYAoR/WrFmjZSZ+yWmwM4E6+8svv4yPHpopMBQnMt7RRx9N+HvppZe07AYxh/0muFQ0a9asWLMioykGiXwhCASLAHYZiPHrCGq99AbNFWQyDkHWAnV5ol3Oq6++SoiKl0jr1q3TkfpuuukmHZQ88Vpym4nXTN/DCx0arWKELAhCgoAgEC4CbLCnWBXtqdMnnnhC8a4IGWbVsGHDFKdX0e1whgPVr18/1aNHD8VxfhWHm4i3zwJhdeONN6pWrVqpxYsXx7933qRr07lu8rpt2zbFFseK7XWSi48TX6dirFe+EASCRwCaGWQtyMYFId0of/zxR4Ig2NEsodzy5cv15yA9xMUFId2KyPeCQA4RKElOlSKjyeGNJl2XbAQQSQ/2L/A5ymfavXs38ZGN+BiXNvOlaJ3yeYVl7HmNADROyJ106qmnapVz27Zt83I+I0eO1AGv3EKQyo4mL5dWBl0oCMBkH2lMoHZeunRp3k1ryJAh2nnzmWeeSRddT89JhMF5t7Qy4EJEAGphxA2GxzWYTtSJ1UpUt25d2rBhA0H4jPQwLiRhIlzAkUuCQGgIwOv56aef1m4CbiExQxuQS0ewPj7uuON0OhjEtcnAZHRLcnRyAVQuCQJhIoCHd+HChTqeC7IhwBs7aoQEcgjZ2bVrV21hbGrZLIwmaisp4ynRCDRo0IAWLFig065gl3PNNdfQli1bco7J3LlzqWPHjjqQFtLDwOLYNEslBi+MJudLKAMQBIoigAd46NChtGrVKh0zBv5K0OggB3eYhMRw06ZNo86dO9Mll1yiXRcQ1MpLwC4RBoe5ctKXIOABAexoxo4dq4N+wyETKXN79epFtWrV8tBa5irYUSEyH0KCwpIYTAbR+hznzMwtFCshwcmLQSJfCAIRRQBOlLNmzaLp06cTji8IDYHjFXYc8Jhu2LChq4o51bQQlBw7JwRKB4OB9gvtgpn17t1bt5mqnuV3wmgsAZPigkAkEADTQZgHhzkgP9RmDstZj7MQVKxeXYf0POCAA3TKFaRkYYdLHZYCrwglgb9Nm7Zz9oNN7MXdjFq2bEmdOnXSfzVq1PB7jsJo/EZU2hMEcoXALk6nsrl5c1rNgcs3b96sg2shTjDCP0DuAw0R/qpWrUpgJoMGteT82WU5rEPgIx4vLgiBYywdCALhIFCa8yrVbNWKah5/vFGHiOIZlkJLtE5GSyKFBIE8QABcg6PsmRKKCqMxRUvKCQKCQAwBYTRyJwgCgkDgCHz7LbEAxrgblhkTn7ZCITk6hQKzdCIIBIwAJ4PjXCvEUl/jjjgNN2udjItnVVAYTVbwSWVBICIIgGMceKDVYFBcGI0VZFJYECjhCHhgNGxmQ9u3h4Ob7GjCwVl6EQSCRQAco2JFqz5QXBiNFWRSWBAo4QhARsP5omwIxVEtDJIdTRgoSx+CQNAIsAUwm/1a9YLiqBYGCaMJA2XpQxAIGgFhNEEjLO0LAoIAOzRZH51kRyP3jSAgCNgh8OuvRPvua1UHjIaduUMhOTqFArN0IggEjMAffxCVKmXVCYpztIlQSBhNKDBLJ4JAwAiAY1hYBTujQRXwqKBJGE3QCEv7gkAYCGTBaMLY1QijCeMmkD4EgaAREEYTNMLSviAgCGj5zB72+wYPpy1PYNuPzFM3UkkQEAQCRQBMhtOj2BK0TmEwG2E0tisj5QWBKCIAbuFB2OLxxGWNgDAaa8ikgiAQQQSyYDR7hhA5XBhNBO8ZGZIgYI0AuIVSVtVQPIxjEwYljMZqaaSwIBBRBMqWJfrpJ6vBwT2qTBmrKp4LC6PxDJ1UFAQihIAHxyUPfpieJyyMxjN0UlEQiBACwmgitBgyFEGgUBEQRlOoKyvzEgQihAA8t+HBbUGILMFpuUMhOTqFArN0IggEjAAijSOvkwV99x1R+fIWFbIoKowmC/CkqiAQGQQ85E7xkDjB83SF0XiGTioKAhFCQBhNhBZDhiIIFCoCELbAccnCDQGpVixzznlGT3Y0nqGTioJAxBCoVIlo61bjQW3bRoQqYZAwmjBQlj4EgTAQqFaNaMsW4542byaqXt24eFYFhdFkBZ9UFgQihAC4BriHIYEngTeFQSH4bYYxDelDEBAEdlWuTJuXLqXVe+1FX331FX3H+utf2FgGf6XZe3JvTk2Jv+rMkA466CDasKEVMxr2kQqBSimmEPqRLgQBQcBHBHax0HfZsmW0YMECevPNN+nDDz+kXevXU826dal0zZp08MEH0wFsWwPGsh8Lin9jQTEYDl43bdqkGdG6deVpx44ldPjhjahly5bUsWNH6tSpE9Xk+j7TeGE0PiMqzQkCQSEA5jJr1ix69tln6YUXXtDM5JhjjqHOnTtT06ZN6dBDD6U9LYPLbGfV08qVK+ntt9+mN954gxYuXKjb7dWrF51zzjnUsGFDP6YjjMYPFKUNQSBIBLawMGXs2LH08MMPU6NGjeiMM86gnj17Uq1atQLpFruk559/nh599FFq3LgxDRo0iM466yzaw0NM4j8HKIwmkJWSRgUBHxDYyqrqUaNG0ZQpU2jgwIF04YUX6gffh6aNmvidYxCD4dx77736qDVixAjq378/x0G3S1THnY0XrZMR5FJIEAgPARyR7rrrLmrSpAlnud2XPvroI/0Zu4swaS8WKmP3NG/ePLr//vtpwoQJ1K5dO3r33XethyFaJ2vIpIIgEBwCn3zyCZ177rm0//77ayEvjkpRoC5dumgZzqRJk6h79+40ePBguvbaa7U2y2R8sqMxQUnKCAIhIDBnzhzq0KEDnX766fTyyy9reUwI3Vp1gSPcihUr6LXXXqPhw4fTDz/8YFRfGI0RTFJIEAgWAaiojz/+eJo8eTJdddVVwXaWZetVq1YlMMUdO3booxRkOZlIGE0mhOS6IBAwAhC44jjy3nvv0Yknnhhwb/40D4HwuHHj9A6sPAe1+R4xJ1xI7GhcwJFLgkDQCCxZsoROOeUULXCF8DcfCfKaA9kN/IYbbkhnxyNap3xcWBlzYSDwM6chgLp4zJgxWsOUr7O65557tKD4zjvvTDsF2dGkhUYuCALBIjB69Gj64IMP6Mknnwy2oxBa37hxI7Vo0UILiuH+kERisJcEiHwUBEJB4JtvvtEuA/BRgoOjF9q5c6f2Y0quC2dKOFFCdpJMcDn4448/OA5N8UA02GGB9kFGBQ8Exvnpp59qC+ak6sJokgCRj4JAKAjguAH/Ilj92hIEr7fffruuu54dKR3C93379qVy5cppRvMTZ66cOnWqNvrbvXu3Pqa9+uqr2qu7d+/e9Nhjj+mqMBAcMGCAtvj9lTMpVOPYEffdd5/TrPHr119/TfXr16cvvvgimcmN53S9QoKAIBA2Aq1bt1bz58/31C0ftRQzFMW7kiL1r7jiCtW+fXv9HTMW1apVK3X55ZfrzxMnTlSox0xFPfjgg4jYoJYvX66vjRw5UvGxR7/H9bp166rx48frz7b/mIEp9slKrjZO1NvG/FoKCgL+IIBwDfCYZmbjqcGzzz6b2rRpU6zujBkz4t9D/Qx3AWYw+qh05JFHEurBMRJaLlyvzPFrQPAG79Onj36P63CgfO655/Rn23/wJl+0aFGxasJoikEiXwgCwSLw8ccfa/mMV1lIutFB7uPIWVDmiCOOIMhxPv/8cy2oderdfffd2igQcWdQHuNJjEEDmZEXfya0f9RRR+nYOE5fzqswGgcJeRUEQkIADKFChQq+99atWzct93EaLlOmjH7rxKjZsGEDYccB+c7MmTMJ43j//fc5ccKu+O4GFapUqaKvffbZZ05Txq+YF9pNJmE0yYjIZ0EgTxE477zzNONAOImbb76ZWD6jBcO1a9fWM4Laedq0aQRB9Nq1a3W8GQh/QfDUdqhs2bL6rYdwEE4TxV6F0RSDRL4QBIJFAKplqKD9ptNOO017fDdo0EBHxkNYic6dO8fjx0D+gmPRZZddRtA6TZ8+XUfTwzgQ+8Yh+DAhBGidOnWcr4xfMa9UqnMJE2EMoRQUBPxBAKEf1qxZo+Ujfstpjj76aMLfSy+9pOUsiCucig455BAtGMZuB+rsL7/8Ml4MKnMwKS8R9eBS0axZs3hbzhvZ0ThIyKsgEBIC2C0gxq9XgSuGCc0VZCupaN26dToa30033UTNmzfXReBtDeYGgoB49uzZNGTIEO2bdPHFFxex54GtDaLpeSF4oUPbVYySFd7yWRAQBIJHgOUkilXKnjp64oknFO+KtC3MsGHDFKdW0e2w8FbdeOON2n5m8eLFRdqG3Q3LXBQfpRQLjdXcuXPj1zkzgurXr5/q0aOH4vjAisNUxK/ZvNm2bZtia2TF1sfJ1caJr1Mx1itfCALBIwDNDLIWZOOCkDxKNsAjaJrSeYHDchcpWBIFv4lt/PjjjwRBsKOtSrxm8l5cEExQkjKCQMgIlCSnSpHRhHxzSXeCgIMAIunBjgX+SPlM8KPioxfxMS6uxUqej2idkhGRz4JASAhA44TcSaeeeqp2Rmzbtm1IPfvbDftKaaGyWwhS2dH4i7m0JghYIQCTfaQxgUp6KefNzjeC5gpWxs8880y66Hp6SiIMzreVlfEWJAJQCyNuMB5aMJ2oE6uViL28CW4NECIjPYwLSShPF3DkkiAQGgLwQXr66ad1qlu3kJihDcilI1gRH3fccTodDGLeZGAyuiU5OrkAKpcEgTARwMO7cOFCHbYB2RDgVR01QgI5hOzs2rWrtj42tWwWRhO1lZTxlGgE4Ke0YMECnXYFu5xrrrmGtmzZknNM2MCPOnbsqFPjIj0MLIcRLtSUhNGYIiXlBIGQEMADPHToUFq1apV2F4DfETQ6yMEdJiExHLy94Zh5ySWXaLcGBLXyErBLhMFhrpz0JQh4QAA7mrFjx+qg33DIRMrcXr16Ua1atTy0lrkKdlSIsIeYwrAyBpNB1D0vTpZ/9ibByTPDLiUEgWggACfKWbNm6fAOOL4gvgyOV9hxwGO6YcOGrirmVLNAQHPsnBAoHQwG2i+0C2aGUBJo0wcSRuMDiNKEIBA6AmA6CAHhMAfkh9rMITvrcRaCitWr69Cc8GuCsHa//fYjdpzUYSnwipAQ+Nu0aTunst3EHt7NqGXLltSpUyf9V6NGDb/nI4zGb0SlPUEgVwjs4lQrmzksxGoOXL5582YdXAsxgRFSAnIfMB38Va1alcBMBg1qyfmzy3JYh8BHPF5cEALHWDoQBMJBoDR7Z9ds1YpqHn+8UYeI8BmWQku0TkZLIoUEgTxAAFyDo+WZEooKozFFS8oJAoJADAFhNHInCAKCQOAIfPstsQDGuBuWGROftkIhOTqFArN0IggEjADHAebQecRSX+OOOEU3a52Mi2dVUBhNVvBJZUEgIgiAYxx4oNVgUFwYjRVkUlgQKOEIeGA0bGZD27eHg5vsaMLBWXoRBIJFAByjYkWrPlBcGI0VZFJYECjhCEBGw/mibAjFUS0Mkh1NGChLH4JA0AiwBTCb/Vr1guKoFgYJowkDZelDEAgaAWE0QSMs7QsCggA7NFkfnWRHI/eNICAI2CHw669E++5rVQeMhp25QyE5OoUCs3QiCASMwB9/EJUqZdUJinO0iVBIGE0oMEsngkDACIBjWFgFO6NBFfCooEkYTdAIS/uCQBgIZMFowtjVCKMJ4yaQPgSBoBEQRhM0wtK+ICAIaPnMHvb7Bg+nLU9g24/MUzdSSRAQBAJFAEyG06PYErROYTAbYTS2KyPlBYEoIgBu4UHY4vHEZY2AMBpryKSCIBBBBLJgNHuGEDlcGE0E7xkZkiBgjQC4hVJW1VA8jGMTBiWMxmpppLAgEFEEypYl+uknq8HBPapMGasqngsLo/EMnVQUBCKEgAfHJQ9+mJ4nLIzGM3RSURCIEALCaCK0GDIUQaBQERBGU6grK/MSBCKEADy34cFtQYgswWm5QyE5OoUCs3QiCASMACKNI6+TBX33HVH58hYVsigqjCYL8KSqIBAZBDzkTvGQOMHzdIXReIZOKgoCEUJAGE2EFkOGIggUKgIQtsBxycINAalWLHPOeUZPdjSeoZOKgkDEEKhUiWjrVuNBbdtGhCphkDCaMFCWPgSBMBCoVo1oyxbjnjZvJqpe3bh4VgWF0WQFn1QWBCKEALgGuIchgSeBN4VBIfhthjENuz5+ZXuDTz/9lL766iv99x3r+X5howL8lWYvs705hR/+qvPCHXTQQVSrVi06+OCD7TqR0r4i8AcHtsWaffnll/rve1aZ/MS+PYlrti/bkmDNatSoQbVr19bvfR1ExBvbVbkybV66lFbvtZfRfb1hQytmNOwjFQKVUkwh9JOzLnaxcGzZsmW0YMECevPNN+nDDz+kzz//nBo0aEBVq1bVDOQAtkEAY9mPBWq/sUANNy9eN23apBcMr2BGzZs3p5YtW1LHjh2pU6dOVLNmzZzNq9A7/vjjj/WaYd2WL19Oa9eupbp168YZSWV+qMqyI2HimmHdsFZgRl988QXHgfqdmjVrRm3bttXrhXXDWhcCpbqvd61fTzUZo9J8X+KHMdN9vW5dedqxYwkdfnijoO/r8QXJaLAIs2bNomeffZZeeOEFDfoxxxxDnTt3pqZNm9Khhx5Ke1oG4djOIvqVK1fS22+/TW+88QYtXLhQt9urVy8655xzqGHDhoVw/+Z0DvhBeOKJJ/S6YQ07dOhAxx13HB1xxBHUuHFjzVhsBriVBaOrVq3S6wWG9e6771KbNm0Ia9anTx/WuBxo01zOy+bxfT2eQ1gUDm3evFkNHz5c8fZZHXvsseqee+5RGzZsCGyC8+fPV1dddZXinZHujx8SxTdDYP0VYsO861APPfSQat26tTrkkEPUyJEj1dKlSwOZKu941DPPPKP+8pe/qPLly6tzzz1X8W4pkL78bLQA7utxBcFotmzZogYPHqwqVKighgwZovhXzM91ztgWH7PUtGnTNLPh3ZKaPHmy2r17d8Z6JbkAGPKDDz6o6tSpo0466SQ1c+bMUDH7+uuv1Z133ql/lHiHo/hIHbnlKKD7Or8ZDQsI9c1SqVIlNWzYMAXOn2uaM2eO4i2/YrmAeuedd3I9nEj2D1ywg+GjrGK5WU7H+PPPP6uxY8cqlvmooUOHqh07duR0POi8AO/r/GU0LBxU7du3VyeccIJavXp1zm+O5AE8/PDD+tdy9OjR+sZJvl4SP2MXc8stt2hcJk2aFCkIWJ6jzj//fMVKAsVyuJyNrUDv6/xkNK+//rqqVq2auuOOO3J2Q5h0jK0vfrWvvvpqxcJkkyoFWwbHy+7du+vdXhR2numAhpwN99bjjz+erkhg3xfwfZ1/jIa1PVDHq9mzZwe24H42DFnNoEGDVJMmTRQetpJIOJ7geMvmAHmxu8OuAgqFp556KrTlKvD7Or8YzXPPPacqVqyo3nvvvdBuAL86uuiiixTb6ii2x/GrybxoB5oeyKz69euXF+N1BsnGgfoHDUfgoKkE3Nf5w2jYBkJVqVJFsS1L0OseWPuXXnqpGjFihIJKtyQQdnNgMDfffHNeTheyP7YMV/PmzQts/CXkvs4PRsOm5uqwww5TU6ZMCWzBw2gYwlC2KNYC0TD6y3UfTz75pF63nTt35noonvtng09Vr169QGRsJei+HpcXlsGsuaEPPviA+MbNuXVmtgPYuHEjtWjRglasWKEti7NtL6r1WS5D9evXJz4WEKuyozpMo3GxNkq7Ptx0001G5U0LlaD7OvqWwTCsglyG/Vc8/Sqx851CG26USgsCe4pMNhXQJC1atMit6ZTXRo0apQYOHJjyWqF8yQ+l6t+/f6DTgb3JunXrAu0DjUP1jXuQnXB96yvI+zrTPZ/uvsYOC39eyeW+jv7RCcZUffv2tZ47jimoB7kO+zWlFUZia4ybyCEsAvsuqX322Uftv//+2rYiWVsEy+OuXbsqyFygLbClbdu2aRP4QlZ5824mUHsUWBWjD/a21xq9oBkO3BVgSewXBXFfZ7rnoYj4xz/+oa2xE+cBho1nBfK0M888U9/XiddN37vc19FnNLAghU+RLU2cOFFBRgDwcVPydraYX8u3336r1ZhwXXAIvktwZ4BmC+4MqDd9+nTnsnr11VdVmTJldNvxLz286d27twpDo+FhaFlXgcHb4YcfnnU76Rrg46e68sortdyEvbRVuXLl9AOUrrwf37NTpjrqqKP8aEq3EcR9nemex/MAhgJTg0SCfxkf5/VXeF7YS16NHz8+sYjx+zT3dbQZDewvOMaIp+1corMcjl2lSpVSuEETCccXaIESdzRQQzuEnQyY0AUXXKC/whELxlx+HHvg8JnYl9NnIbwGPbePPvpIcUyhOFTsma/dCOJfBPAmm3sxeTjZtOV2X7tdc8YwZsyYYowGNl6w2HYIBqbwP/NCadZ+XKQj7CEmCUI68DHGVL4WLweBq0N333038U6lSPyYGTNmEDv0acGsUw6vDzzwQPzjXhxACHFrEGLCucbWvlrIee2119J///tfHfMkXsHiDf866tg4FlXypihi/mB+QRFrIDk5fSw7PcJ2sFaL+GgTVHe6XcQrwr24Zs2arPsJ6r7OdM+nGjiE9hhPYmwlBHtDSA0vlO6+jjSj+eabb4h3FF7mq+twiAjNJG6//XZi72BCeyAWxBFvDYl3M/pzun8ojxurZ8+eushrr72mX3lnQwi8xNt3uvjii9NVd/0e83LG41owDy/ykTSUWC+XX365Xl+25NVBzYKGCvFrMLdsKaj7GuNKd8+nG/P777/PiRN26fvZKcNyTX1vfvbZZ85Xxq/p7utIMxrj2aUpiChjHL6BeDunI7Q9+uijuiQYBG8VCTsWN2IpOl133XVxZocgShzLhO677z76+9//TixYI45voqPxubUj14JBAOphlpkRHox///vfwXQSwVbT3dcYqtu1VFNBWFtQ4rOAyIUgFjfoVz/+RZrRsNBKh9D0OtE99thDx/y97LLLiIVUxEJdvSXEzgTMgmUvNG7cOPrxxx/1ew5fEO8Kkfnw64Ujl0MI/YkttEMI5/nDDz/oqHvOd6avCA2K+RUiscyLENM3aOLgVdSlSxdiDQ699dZbnADAPAOAl7FhTphbthTEfe2MKdU971xL9QrGBEI0QodY86rvc4gWbCndfR3p4OSNGjXSRxecI73IaRJB4uhthEXAkee8886LXwIwCEgOJuKc+1nYSKwR0mElUZDVf8SqOx0zmLVR8boO5/dy8y1ZskTHs403VkBvEKcX87vwwgtDmRViCSN2MNY2KEI8YhyjIafJloK4r1ONybnnU11zvkMQd1Zw6DjLznfrOfYwQqfiebGltPe1F8lymHW8qgHhcs9CLj1U2MZANZkqRCTHFS6idYJRFt8IaurUqTrqGwuNteEZx7NV8EuBatsJD8pCZsXxbD1FhkujBgwT2sD6Clq9zbID9corr8RxR/jW66+/PrD5oOGoqLfd7mu3aw440C6xrMn5qF+BHQfdj3/Xrl07HTEy/oXFmzT3dbTV25gf1GUcSNpiqrGisBeAShvxYLp166bmzp2bso1kRgPGxly8yJ9jY4AGWM6jo+chSBLUgl7ChroYNqUcYz5+GaTBHuL+8q+t9kGCGva2224LPPyE3wZ7QdzXme55xNrBjyjub0SkdCydYcYBY70ePXrokCawJfNCLvd19BlNNqbamHiyVa8XAJPrwPua03okf2382cVU27iNqBcM2gWB5WqBODqmwjVqLghu97XbtVRzS/yO5Y1F7JMSr5m8d7mvxanS9gyabXlxqswWwfDri1NlZswz3NfRd6oEJy1B7vQmPxx5U0bCRLgvVQm6r6N/dHKWCoJY5E9avHix81XevUJoyQnRSlzgKwgI85FY+xJK4KsScF/nD6PBjYqQh7yJy8tQnnAChKMhPGhLEiGUJ4KWITB5PhFCecLzPwzH1xJwX+cXo8GNirAMcIL0EgcmFzc6wlmyrYJmkBBglkQCBkceeaT2HEZIgqgTB1nT64WjX1hU4Pd1tJ0qU4mg4OD49NNPa/8jjg+SqkhkvoO1JXJHw0CLz+PE8W0iM7YwBwJTduQrhxU1mxsEbsGbzdyYueg143QrdPbZZ2fTlFXdgr+vw+LYfveDlBiSQM5vVINtD7FOJIGcO8YFel/n39EpcZmwDUfUMwTykZS4ichE+72TEpf9lCKREhcxWpASt0GDFWrNmtwHUi/A+zq/GY3zOCEjJMJqIkgV8id7sdZ12vLyCqNAPs6pY489VrEvjJo8eXLcPN5LeyWhDnY3EyZM0GElYd3LYTxCxQyGoMh0ikBmp59+uuIYOuqhh5SqUUOxsiEaK1BA93VhMBrntkAEPKiQkWWwd9u2ahybpjt+SU4ZP18RYnTQoOt4R9VMuzrAxBsPkJA5ArCyfoifcLh+sBOgQljJVD5p5i2mLwkNGIcNURzqQ8dsZudaxdkoilRg1zY2o1AcvrXI1zn9kHhf48cM8YaDvq/xgw21O1x4fLiv88My2EqqxoV3c2ySbRxBb3jz5vT888/rGB0QtkEQCc/ihg0bEgcst2oWIQIQjwYR3djBTgdagov9aaedTQ89dBVNmLAXnXSSVZNSOAkBdlzVHvPsf6aDMXGGSx0Ggh1XtTdxYoiOpKopPyJsBNaMNTrEPwq0gcOA1Dv6aOp1xhk6bAgCWaWi5cuJTj2V6OqriYObpSqRm+8QoGrWrFk63EnQ9zXv8jRGeFZ8oPF5kdfJaqLs4k5HHkmEcA4cPgCLgxvYYQ7ID4UoZAjRidAQBx1Um8qVq8B/e+pQA3wMIoSlKMXhIT9nrdHGr74idj7TcXHApNjLlRCHBn81atTQQ0Pgvb/+lWjlSuI2rEYrhdMggPCSWDP8IQdWBf58YPXq9D2vKXBHsCuE6UB4CGfNEMSJfdB0yIM1a1rzj8ksDu1Rl9q0aUO8E6DjmdnszYGymOsQcSwbN/riC6KTTyYaOJBoyBDiUCJupcO/luq+3vz551SPc2lVZJwQmvOAAw7Q4VUSMQJWHENb/23atJ3jBm1ijFLf1z7OqgAZDe4ODoZEw4alxQk3JBtkaQYycOCRrCqfR7VrryXEHEFsGsS+ab14MdXlu23jXXfpRatVq1ba9nABIWsRx4qLCwWAwG6ObLiVd6ir27bVDwniCMFkIHHNsG58bNaM6JJLjqD//KcMtWuXNJjBg4k4fCW9/DLxQiddLPqRtfF01lnEQaCIWNtN++5b9HrUPu3q25c2AyNmrAg3C4zwo5mMEX5gwawHDWrJgd/KFsfI/4nlh6+T8QH5qaeUatpUsY2/URWkVGaj1dQE4zqkYdmwIfX1pG+Row5n+yVLki7IR38QOOEEpV56ybgtjnigONRQahowQHEuHcWu/amvJ3wL+0JO6aTatFGKgwFEmzgcipo923iMMNYOSRaVfwZ7aZkt/7px3E2iBx8k3jOnLZZ44d57iTjKZ2qCcR2nQuWAw6mvJ33r7GYuuYR9JGBXKuQvAvwLzdsV4zY5aBwbBqYp/tBDiOJNbJGH8IlpCsW+xpFp0iTS8jfsjj75xLV4bi9iwpi4IbliZNiGabE9TAtGvhyCU/M5nFjYZ0KcAptefz125ElbHlwINyXLa0yoX7/YNpuzsAj5jQBi2nIaEFOC+CwtowH3wFno999jAhiDX4Z//Yto+HAilk9z3GnTUYRcThhNwIAjLcT99xPddptxR5xthSDocxXeIjjz8ccTBxC2anfkSOJ0FcZVpKAJApwihyODm5TUZZClx3UNkAGDM1gQx4fWN4IBs8EGd+JEolNOIXrpJeOhhFcQqWBY/mJK2CAC1jCoMHY0UAtA+Gv4i8dB3jVfuuIKA4jRLhhYhi220xIrpjglCxHnlxPyCwHsKJEwLkN6nMTuoLnOmIgBaUU42wVBU4n8XAbMBmrvF18kzppB9Gf2nsRuc/ceGAEf7NYMqVw5A4wM28pULP8ZzSuvEKfaIxo6NNNc49cnTCDq1k1rv+PfpX2DjIuIfD91atoiyRduvDH2i8eJAIT8QAAcI43NS7rmjRgNKkPzxLYp2jbh1lvTNVfke1Z80bx5RNi5RkbLGCRGRWbv7UN+Mxq2kdEWVXfcYfxrh40Jbg4X7XdxJHE454RzJr94qIxfiptuIoIm1eBHsnh/8k1RBIJ+iHB+nj2biI089aIV7T3lJ/z2vPlmTPcQid2rB4zYzIa2b085Pd+/zG9GAw0TDpr/93/GwHAUALYMJmrVyrhKTE6Dm5GtjE2pf39iY0GrjZBp0yWvHAxaLFMjoziqGRN+HWDMB/sanI8MCLnX2OhY724gv8F654zAMSxkWBgnigujybRinF2S/vlP670rp+G228044+D0t3Tzzc6njK/IJgr1OXZO0LwLZYEA5A8JGUJNWkJxQ2Xh/5oDs+EkcWzt5mL38L/ieIeHdc4cYiNCIrbaZ+O4otdD+xQWRh4nlL87Ghxl4FzEfjCmhA3JHjxjTz5JPXvGOAbO84YETTs07oZHf8NWS2Axtm7NZMWbjApEL6hmTfiFmDs3psP+29+Mzr6wGIaAGKZXI0ZY7qSsB5imQpgYpRmC29f5yWjYn0WrjSx2GACBc8JrAZ4bIGmv4Qa8/noiGFRYEMRH//kPEex2hDwigC2hpf0/inveScIPCscoeFcOGmTEbGAj+thjMdMceMCEpTaOIyqMJg6Ff29wZII60lCdjY5x30Ct3atXFsNgr1+9H4fg0JBgOAZr4euuM6wgxYojgPOIq8FT8SpgNOzS5p1wjIJG88MPiS66yKgd/BbBkBx2Nh07Ejt4GlXzpxAwsjxeet71eRhx/u1oVq8meu45omuusZoudjPY1uJm8EwedzUYKk5c7Dgu5AUBcAxDtxKneZiUZC0vwVkIwmGYT1x+udHOBv2PGhXz5ufIJOy864wo4FdgZLnrA6NhZ+5QKP8YDbgFJKzQzRkSx8XWwro+fQwruBWDOy8Ebxa7Gvw44tSF+CZCHhCATYIlo4Hd2u7dHvpKruKovvErMXCgsWoJbnfYxWJng01R4ASMLH9FUTwsTVl+MRo4mXD4hvSekKmXEwojnLYgCM6asDqwyINtjYWRDE56cMiDf5WQJQJ4GiwsXp3WUcXQoNupkvoVOwX4HEC1xKEYtI9U6pJFvoX18NixMeNQjrkVLOUaowyz8+PRy9CFj5dhGYUzkMVZFEGp4I+H+8M3ggYKv7DTphk3ieKeVevGvRRoQQ87GiABzH1hNGgM99zMmTEJ85lnGguAsAGeNCkmt8HOOjDKgtGEsavJH0aDqGjr1xMNGGC1Vti+QlHky24msWeY/qJxi1WCIBrj4LRUQjYIYOdoeXRC85DT+HJ8csYKf6vp04lD+xF1726sPz/xxNiaY/3xwxcICaPxCVY4loBjWGyhIUaBdWggecDgLAWt1yOPWE3QA3+yar9gC1swdAcDD2ILp2r6VzA8hJhAMBdEczTUoUPlDTsuyAlhc+M74Ujv4dfU4nHKasj5saPxeP7BhgPiFA/4m4EKOx4wPws9qkf+ZDaeQi2Fp8EDo8FuJpC1x3gmT44520KlCD8jA2rfPibqgewGESp8JUwU8XUsCVqnMJhNfjAaD+efp56KyWph+hIYwfQXgdAh8bMg7GpuuCE81aLF0KJZFA+RhzOQx9OEGQbYQSDCGXY4nTvHBIEGNREMADZdl15KNGWKQQXTIh6ZcaAYJYw9+owGRlOwtOvdO2HY7m+xZYamySIOlnuDblfhCoGOEHTIkBASskUL4hQthhVKejGPjCawHU3ietx9d8zJCaH3kDrBgA4/POaICd0GAjj6QlkwGg/iL+shR5/R4GgC+Qx+QQwJi1e3LnGydsMK2RSDKziYILRhFoQdDXY2Fqcui9YLrKgHqS7kx+BP+AucYCSFmCCw0INTpgEddhhxKpmYfRWCEGRN4BYW5hboD8XDODahrzCWAf14o3nziLZti+W8MGwBLh94iEN1ZISRDs7s69cbjjIWpgKnLtnVGEDmwRUb9wGURKERLIdhEgwvWvhIGdAhh8SYDX5LH3jAoIJbEWjCDAXTTjNhYhRtRoOA4zCMs/hZgriEc7tp0YkDaOCviNMKU1DLCEjgT5Anh2UGHjgOQXUARmPpT+DBxzD70SO517hxRCecEDMsNWixXr0Ys8GtjqqeyYPjUpgYRZfRcOpZvQ1FBClDgmEejOJwJAmdEEr0rbeILKyysKNBAC7Z1WRYLQ+MxoOPYYZBGF6GsQzCvsLOxjCCOXY2uG0g6kMMI08kjMYTbDFuAYmuhaQKyqnzziPCr0TohIcBXM7C+Q5jhPod8mTfLFhDn3gIHXrw/vPgY+jfRLp2jRnLDBwYix1h0DISbsC+BrcPfIatSRiNNWSxFBjY0WArakjIew2jTciNc0YQCmPBJ00yHgKn8qamTWMiHuNKJa0gvKgtY07CUNMysoS/qHJaWp3jGz+W0EwZEDJowB8Oedzvu8+gQmIRD3ExsOsLC6NoHp0QLQpSfOwSDAme0XDstgyWb9i6RTHsfSFXgkrekCDawbbZg6mIYQ95XsxDFG3wJQsH/2AAgmoJx2nY2+CeMKD69YmQPQOPgJXMBpO1MLHAUDg1NyHGVxgUPUaDeK3PPkuEMIqGhJAhyCEHI6icEwQvCJaOM5Ehde4cY5AIwi+UAoF8ZTSYCiKYI10CAgvDJNhABY1jFFTfOIkba6OM88v8D18PiRP+V9nyXfQYDUKUIQObYUR3aGyw+YG2yUKcYwmTZXGokpDdctUq44rYYedEiG08whwWxM+uVUqD2Ekr5zsaBzLcy4hDjJB7UDVC3ZOBateO8SaYZxklShVGkwHRxMuwA4AK5u9/T/zW9f2YMUSNGxPBQzYyVKVKzJjHYovVo0fs/ps3LzKziM5AEDkM0l0LfyfwpZwfoxMRhAwFWTHXrYul7zE45kAbhY0QRAIZ8xdC2IJfXQuMcLwMC6No7WgQ3RnWldg7GtBXX8UM88BsIkcIao29KTx9DenKK62zxxi2XADFKleOpUExnApO4OD3kSJYOCOPLgy94Idi4LIAw3NshiCDhETBlSpVMva5QjuwhUWVMChajAbHJuj3DAmOswj8nRN1dqYxwshw/PjYHWIoGIbJ0KJFMc1DpuZL3HUYRcJQypBQ1CLfvWGrPhWDPcNllxHBnXvZsoyNNmoUM8nBve6a7QehK7ZsydieUwDMGPkXw6DoMBqcGaB2MXRQwpYSsbCwrYwswbsbWQ8N8+9CyYabKZI7tFyDXEiMBlhCsAgNJayI4TicgeCICSYDO7G0x2twDXAPQwJPAm8Kg/YMo5Nf+Xz9KYeD/4rPOvj7jvVqv7ASH3+l2atrb37CTuMjxg+dO9Nv/JNeq1YtFtaztD4N4SiKGLzIl2QZ+D1NiwF+Db011JzYrmC7nIb+YIs9YNSy5Ta2o2jFTqFTWEHxfRGM9uXJVuebqQbncKnN0kK8jyQhtq6xusRwBrClQehUqIsN6Phf2lDLeWuIVvDx1S+qWdM49YpRlwgJizXEKzQByKvrQggxgeiMCH0CptO6ddHCu/h4uXnpUlrNR7R0zxmeNdw3B3HQtg0bWjGjKVu0kYA+lVJMfra9i4VRy3g7uID1c2+yWu9DDgH/+eefU4MGDXgrW1UzkANYHYAJ78cCrN+Ya4Dh4HUTS+UBEF7BjJo3b84PXkuOJN+Rj7WdqCYWmglOaCtWGJxZ/ZxYNm3hAYGHLwaNczrTx5zCAxjhbzk74a1du5aZS904I6nMN01ZdpRLxAg4AZsv+UH+gs/3v3Ogo2Zs5dW2bVuND3ACtjknqPhxdLTIiZ7zMZsM4IYbYgJdv+eFnCzY2ZxzDhGcnjIQmMzAgYo3RCv5fngl/pztWr+eavI9VJqfE/xQZ3rO1q0rz+ZeS+jwwxulfM4yDMPm8nhfGA2Yyyye/bMsrXqBJeuY5DEs1O3MO5SmbPZ66KGHsurZbvO0nUXiK9nc9222EH6DHUEWcjZ1tNup0wU0adIlfK00Mx6buea47Gmn0Ze8UxvLGhTgBMw6cAyT4/ioeASn9W3MqjMwFhvayoKIVaxCBz5gWO9ylog2bJHai/1t+nDMyAPDUikkDxoPIs6Afj+Qyf2E/TnIeX3zDRHfI1p0AF+aFPdC4nP2/PNz+Ye7KRevG8hzhnvoHGZ8DSGNzp7G8/bcO23evFkNHz5c8VZMHXvsseqee+5RGzZs8N5ghprz589XgwZdpypXrqv7e+KJJxSDn6FWbi/zrkM99NBDqu1RR6k2tWqpkSNHqqVLlwYyKN7xqGeeeUb95S9/UeXLl1fnnnuu4t1SIH25NnrqqUq98IJrkby8GPS8eP1U795KtWun1NatcYhy8ZxdddVVik8gfj1n4zwxmi1btqjBgwerChUqqCFDhij+VY2DEsYbPmapadOmaRB4t6QmT56sdu/eHUbXxn2AAT744IOqTp066qSTTlIzZ84MdYxff/21uvPOO/WPAP86KT7CGo8964JBP5BZD9BjA2HN67rrlKpbV21du7ZQnjM7RsMCS33zVqpUSQ0bNkyB0+aa5syZo/gIolhOod55551cD0f3j3G0bt1a8dFRsZwqp2P6+eef1dixY3kXWFkNHTpU7dixI/jxhPVABj+Toj2ENC88ZzP79lUX7L9/oTxn5oyGhZWqffv26oQTTlCrV68uugAR+PTwww/rX+/Ro0crLFQuCLuYW265RY9j0qRJuRhC2j5ZnqPOP/98xUJ5xXKvtOV8uRDSA+nLWG0aCWFeBfqcmTGa119/XVWrVk3dcccdNssSelkc6bCLuPrqqxULk0PtH8e57t27691VFHZ66SYPuRbO3o8//ni6Itl/H8IDmf0gPbQQ8LwK+DnLzGhY2wP1t5o9e7aHlQm/CmQ1gwYNUk2aNFF4+MMgHE9wnGT1e852UzbzxK8mBPhPPfWUTTXzsgE/kOYD8blkgPMq8OdsHBs7pKfnObUe/0rTe++9x06LJ6YvGKErpUqV4jge47TqmDUv7G7ko8FWinnCGLEbZ4VjgS9t3LhRGyCmKBapr2DTBBunszmFJx/xIjW2kjiYEvGcpeP5bJOhqlSpotiWJV2RyH9/6aWXqhEjRiiomIMg7J769eunbr755iCaD7xNyNrYQlTNmzfP374C/OX3d6CWrQUwrxLynKU+Ov3000/qsMMOU1OmTLFciWgVh3CWLYq1gDaIkT355JMap507dwbRfChtsoGlqlevnr8yrQAeyFDAyNSJz/MqQc9Z6qMT21+wWfLh1Ldv37zeye6xxx6cdnQKh8m8TR9r/JwMy2XoSo7r8Mgjj7C/FccayVP6P7Z2ZWNLDpB+S57OIH+HXaKes2QmDkOvihUrKvanSb5k9Bl2GulsNcDB8ZeKILj99ttvU13S32ejRRo1apQayM4hftJNN92k+vfv72eTGdtieZPC+vhNUH1jzdnPzJ+ms/zld5un2zVn8LhXFi1a5Hz07zXLeSUOJNvnzGkrlYbT7TnDM5buWXKr5/Tn9urynBU/OsG4i3cybu2lvAbmwr4Rap999lH7s6ERbDYcrQ/sWtAm5BlnnnmmguzEIchPwASgkr7wwgu1Fe22bdv0ZXasVCeffLJioaVi3x3Vo0cP5eWYgvZgkp8OYGcsNq/169f31R7l+uuvVyzILvIHHHH8wx/wg8yMfcY0jjZjNSl7LrsrwJLYF/L4QLrN0+2aM2ZYqHft2lXfX9DigNxwdeoZv3qcV6r2vT5niW3h2IsfCIfcnjO3Z8mtntO2yavLc1ac0cCiFT5FtgTfCLglsIZKuyXwhlZNnz5dNwP/nhYtWuj3uGHYS1mNHz9ef77//vsVO2DGuzv99NMVOCPoiiuu0EaCeA/Ba6tWrdTll1+Oj9bUm31IYNTnB8HgjY+WfjSl2wAmMIaEbQs7SOo/PpZpxooCEydOVJAHoRzcGoCt3z5M7JSpjmJ/LF/I4wPpNk+3axjzq6++qsqUKaNxcuaQCVennPGrx3mlat/rc+a0hZ0JTBTgBuSQ23Pm9iy51XPaNn1N85wVZTSwB2F5Q9rjjVtnF110UfwydjIA4IILLtDfwaYFFrMOYfcC/x8QfnFgQObsfk477TR144036mvwE8ID5xB2Qhw2wZMWCQ6fiWN02vTy6mdb6B+7wWTGAf8kPFygxGs40mLnw6p0fc2vf9msfbExeHwg3ebpdg3HBxiUJh+PM+FabNyZvvA4r+Rm/cAac4VGNXFH4/acuT1LbvWSx57pc5pno6gwGDFSENKBjz/WErYHEgId7bXXXjr+DEJFQGiKdp1YMmgYQXcQ0gB0Bkfx4S0X8bGI4xpNo48++ojY4E5f+4Zd51HfIYRT4KOTjm/jfGf6yr/WOjaOaXm3coixg/b8IsSc4R1fvDnEnXnttdc4OB9H52NKvHY3JyPj3WMRPOMVs3iD+EBY+zVrOFhUjshtnm7XcO+xVTjxcZbTn1/LaZT+q2P1ZMI1R9PUz4PX5wxjnsF5eZhxFLkvMj1n6Z4lPJtuz6ctRumesyIGexgM70Rs2y5WHu3ghu3JkcPef/99HXsFgZwcYlkDocxnnIwJ2i0+r3Ka4pforLPOIg5zQE5ZGMIhHo1DvDXWb21j26AS5oU+/SDetgYa64WPATqOD/9Kx4fL4Td0jJ/bOdkPe4L7Npd4B/wG8Wswt1yS2zzTXQNTBvHORt870AZejBCMSZQK16QioXzM5jljITKHoh7PIWyLxrDN9Jyle5bYStz1+bQFJN1zVoTR2DaarjzLWOg6Dt6DTmE5C8IuxyEnwBOseFn2oqPLsdBXR5MDc1rPkcJA53GAVADIQmJioziOW345lePAUQhjWcj0HCdfhkV2IiHoF3Z8vDXVeD2KaPoFSG7zTHcNwb84Bg+nkb2PM/X8nf7BSbLwg8XH8SIIpcK1SIE8+AAmClOExOcJw870nKV7lpyNRWJ7ic+nX5AUYTTsr6NDaGbTOCLs4VcR23sQbg4QosE5xGdnHcoT2z/sZDjUg47Mhyh6P/74I0czjIUzZHmNNpWHyTwifSEKHTtNEhiULSE0KObnB/G5ODDXBjBeYJjMaGAThCPnZRw9nwVunGd8uh9TKdIG3DUwt1yS2zzTXUOIWBz9HELY1x84sVPibjgdrk6dMF+9PmcQN2D3BobK8k/taoPnBe+d3Xq65yzds4T43KB09WxxSfecFYmv2YjzOuDIg/OeFzkN5Cus2SH2ENbjY7WZ5rw4AiDOrUPYsYBp4MbBQ4XjE4KUQwZzCYeAhBGcQ0dzJgH8gSEBaMQj9kJLOKEx4uv6QWgH7WGn5Tex/QexWlsfndK1fQhnFgN2fhLkQlh7yA6iQm7zTLyG2NLwx3PI+UVOZJomuDr1g371+pxBpICdiUN4qPHcgNGCYbg9Z6iT6lnC6SBTPac/k9e0z1myFNmr2g3GXgygmjp1qo4mxwIrbdDGjEFrljjIeLyrdhyqEBHyQAhbwCDFr0HjhFCUicSyHO0ZfeuttyZ+bfU+jdrNqg2nsN/qbaddvPLWv4imDd8hfAAL7PBWa6ighvY7HGgU1Ntu83S7Bn8hqLadMLIsMFf8o1UkomEqXDWgNv980jqhS6/PWeJwOfZ0Ea0TNLjpnjOnXqpnyaSeUz/Ta5rnrKh6G41APQXjOFsCcMzxivyxpkA3A9U1jPVgcIcQDrC5cQjGQgMGDNCGeQgLinL8y6ovAxQwHtjPLF682Kli/epiSGTdllPBb4M9p10+IipEDUwkGOtBpc3HRsVCPTV37tzEy768j4LBnts83a4BAJZZ6SiLMBSFujY5vGwqXK2B85HReH3OEseczGjcnjO3Z8mtXmJ/md67PGfFGY1fptGpBsXnZsVCq1SXFAJrw3oxkbAb8sN73MU0OrE7q/dBuSDwsVIb5nvEougAABfWSURBVCUPBouIGyIIipILgts83a4BF1iZczqalBClwzVl4XRf+shown7OTJ4lt+czHSSJ37s8Z8UZDSqiAsz+C4E4/5HeXuLVT4JfCEIsRCVOcTZz43O/zmaRTRtF6vr4QBZpN9cffJ5XCXrOihrsOcIeaIygVmZ5i/NVXr5C08BHMc5IOyyu/fJrIhCWj+HctXzsI2Y6fjUbejsvvvgipxaerw3dQu+8hHdYkp6zlKoLPESw02CZSREVYb7dF+zDoRPXOap2v8cPA8MjOSsjywX8bjqU9pBBlN0yiPNOEaIRCoWLQIl6ztx2o2zgpIW7cJTMN4KPFBwfk+U+fs8DsiUECUNg8nwizvOtvcH9cjQtMnefjxhF2s7lh4DmVQKes9QymsS1hLs9HLcCie+R2JFP7/m4pNg2QDNINmbyqVX3ZtAn72x0KAdo0aJOH3zwgcYHHuGBUEAPZCBjtWk0wHkV+HOWWkaTuIGEY+TTTz+t/ZYQESzKBOtG5LKGQRTkJjB8C4NgqYz817BGZRW0dvALo18vfTBz0RhxSAodnNxLG1LHfwQK/TlLKaNJhhEPL9wDkJwe2RDg7Rk1QjR/ePhy4CNtRezFsjmbOaE/+NIgNCYsnDlNbzbN+V6XVcPaVB1yKwiAkcBdKFoIFPRzZrNzxLEAUdiQw0hS4qZHzkmJ26VLl0ikxGXtmKTETb9c5lcCPDolDqIAn7PMMppEAJz3yAiJIFTs+anzOSdbYTrlgnqF4Rof5xQH1Vbsm6N491DE3Dyofm3aRXS3CRMmKHYc1UG+OLRDqGOEQRgyi7Ifi0LUQo6hYzP87MqG9EBmN0gPtUOeVwE9Z94YjbNEiGw2fPhwHVIQDz3ioDr+Jk4ZP18RYhSJ6hGRD+b48JPCAx1lgrUqq4+1bws7AiqETfTbT8mZPzRg8CHr1esyVa5cbQVDvBUrVjiXw3sN+YEMbWI5mlcBPGfjSmGRsj2p8sNOs2bN0qELkHUPoSEg3IJgFJ7OCPFgG6wKIQsQZwSu/uzwp8NFoF3+ddZhEtBmvhE8z+HZDlkXMOvQoQPx8UrLdODNnhjqwGRuiCoHjCA/g9EdPGfbtm3L2A+ku+/uzREF9+BIfCYt+VyG5VTshk8ssPK54Rw3l+N55fFzNt4XRpO4/AADD5TDHFiVSpvZMKweh1msWL26DkF5wAEH6DAUCLXIxyAdlgKvCCWBv02btnO8l03UvHkzYm9UQnwR/NWoUSOxq7x+D4E6MMIf7zp0MCu47HPAaT1PRCFEuINEjDhYD33GQl1gxC4VxGd5zcjbtGmjczOBcQFb0L/+RdwuMVPLAUytWhG1b08cACgHnQfY5Y03Esc1iQQDTfWcIQIhYjchbARC52Z6zjjigo4/hc1AwM+Z/4wm1TLv4kR0mzlmyGp+IBBuEXE0EPMGMVAQTwMaG/wBIDCTQYNaclCfstSuXarWCvM7MA32sI0zW2AEFX0iRj04Ts8PjOFvHOjIYUrp0EBwuaZNiXc2OXgu+AbmwL3phlb8+7lzEUia6LDDil9L8c3ls06iC1otoxbVt6S4GuBX2B7+9a8BdpBd04iyx4aYBAZi8pyBGTmBr7LrOWPt8RTK+ZZDG6jZs427gpHtn5lajOuUiIIID9G4sfFUOYwNGy8qzoVlXCX8gs88o1TTpnC9NuobacLZEFsovxDIbLCXkVeZFGBZAofxMimpy6AoqgglIcAyLw6tRxz7NOlC6o9s/sRyMqKbbkp9Peffwhn1iiuIo20TC/GMhnPvvcThTI2KSqEIIcB3bQgkjMY/kDlAO0coN24PKbXvv594S21cJbyCo0cTC5eIOnY06pPFUprHnnuuUXEpFCEEwmE0SOHB8hdTYnkocVYJoVQIcLR/VsERJ7dKdbXYdxB9cJQM9sQvdim3X3zyCdGDDxLddZfxOMAwBw4kFpAbV5GCEUEgeEbDCd84Qjmx1Nd4ypxRhbVOxsVLVsF99yWOik50xx3G8+YMHZyYj+jll42rBF9w6FCiq682/gHixBmEHIVybAp+aYLoIXhGA47BiclsCMWF0bgghuPTY48R57VxKfS/S8i7xxk6OC8WsUr8f9/n7B0Svq1cGZPPGA5iwgRiPzYiTgAhlIcIRJLRwBRk+/Y8RDOsIeNsyYaLbANg3CP7wnLKWKL//Me4SjAF2c5KC4BxZPoz82imjsAcURxHQKH8RCB4RgOOYZmUDMWF0WS4oTgjI0EF82cm0Ayl9WWctiB/zeluEfY1MLzkXOumxJEt2LqcCHaAQvmJQPCMBjKahCyCJjChOKoJuSDALgucEYxo0iSXQkUvNWlCdMYZRJyxODfEls10ww1ElnGNwCCvuSY3Q5Ze/UEgeEbDFsBs9ms1WhRHNaEMCHCOaW0kYyF4gRU9eBMbIYdPt95KxDnWOTWpcd/sOscpkInjIBlXkYIRREAYTQQXxXhI7EDJOWyJI8kbV4GVAU5dw4cbV/GnIFIiQ6aEs5sFoTjH6hLKcwSCZzTsz2R7dJIdjcVdNWIE0c03E3FqGVOCupsjj3K+atMaPpS77jqiiy8mdu03buzVV4mg1u7Vy7iKFIwoAsEzGggrYfthQWA0cAoUMkCgc2ci9vTmwM4GhWNFgC9EJaHJPWDEM3MmEZiiBWE3gyo4OgnlNwLBMxrIDyzvFBSHFlTIEIF//jN2JLEILYRUVDAufv11wz6yKQZuAeM8i2DxHGKHw4UQ9emTTcdSNyoIBM9owDEsrIIdYFDFQsbpVCuZryecEHuIoQc2JPgwwtkycNuUd98leucdoksvNRxZrBjk3OCf8CEVyn8Egl/GLBiN7GosbjDorPFkWshqONGmfpCfecaiH9ui4BjXX28lp8MuC364HMZIqEAQEEZTIAup7fOhUrLMlw7+BD5gwZ/MEZs3L6ZHt0wZDLkxZEgeNsLmY5OSoSIQPKOBwMXD/lduMg/3AaSn2NVYbAVPOok45CORxanLfGDgGDif4ZxmSLNnx6zCe/c2rCDF8gKB4BkNmMzvv1uDAa2TMBtL2BDbBV6HlsnrYMSHGMMW/CnzwKCbRqwPnM8sCLwJ4/Hw22TRixQNG4HgGQ24hYc7GFWE0Xi4HWBTg6cV9kuGBK9oD6cu99axs7KU5kJDjyMc3CSECguBSDMaix13Ya1KNrNp3TrmAwWHSwvCCQe7Gl9kNa+8EvPcPPts4xFAwwi58W23GVeRgnmEQPCMBtzCwr4D2KG47GayuIvANeBXZOECz5ladA6oxx/Pol+nKjgW/izOPw89RFSnDtHxxzuNyGshIRA8o+HcRJw3xAozOFQahiqxarfEFG7UiKhnzxizsZg0fIr+/W/r34WiPSB1yjffWJ1/sN433GA93KL9yqdIIxA8o/HguOTB4TvSIOdkcNhRIPYLciwZUrduMbu/6dMNK6QqBk4Fj02L3cyYMbH45EcemapB+a4QEBBGUwirmGoOSHYGJ0ZLN23saqAl90SwAF67lqhfP+PqW7cS3X47EbI1CBUuAsJoCndtia69lgiGKcuXG8+ye/fY0QnVrAm7GXhqWkjxYSx43nkxrbx1f1IhbxAIntHAc9si3CSQg2ZWUmr4cA8hnQSOUMg4YEGIV2Ot/Vm1igh+TRYpY1EF7g/QxgsVNgLBMxqYnRpG63eg5rTTVL6880les0IADz7nO6cXXzRuBla5OAEtXWpcJRae829/s5Liw6EbTMYySYbFoKRoVBAIntF4yJ3iIUNLVPCM3jhgJzB2bGxXY2ihjZMPNkHGuxowsmefJQKjMaRZs4iQQ27QIMMKUiyvERBGk9fLZzh4qJMQzNwiKyRy1MGLYP16gz6QohcZNA2zXcC9BCm3UQ25BYUKH4HgGQ2ELbizLNwQYGcm22mfb767745tUQzV3YhRdcklsc2Q60hgiwA1OgQ7hgR1NvieBBw3BKwAigXPaABSpUpE0GMa0rZtsSqGxaWYCQL16sU4h0WkKxxrJk3KkPrmkUdiRjAw6zUg8DkcycD3hEoOAuEwmmrVYpGMDHHFkR/JGIV8RgDORPPnEy1aZNQw4ojDJcDVGRznH+TaNSRovy+6KJY107CKFCsABMJhNOAa4B6GhOhq4E1CPiOAYyzkNNBEGcZJveyyWELMlCOBuwGoS5fYa4b/c+YQLVgg6VMywFSQl1m/EDztqlyZNrOudDVL/r7ivfN3rL/+hY1l8FeatSJ7c2pK/FVnhnTQQQfRhg2tmNGUDX5gJbGHM88kmjgxYzLrX9n26dNPP2V+9BWL2FrwpmUB1ar1aZE1O4mNYHbylucP3iHVqlWLM6nwFigNQUwHQ2Xk/rZMipGmRfk6nxAopZj8HPAuFvouW7aMf7kW0Jtvvkkffvgh7Vq/nmrWrUul2SweN+MBbFsDxrIf/8L+xncgGA5eN3HYezCidevKcz6fJZzQsBG1bNmSOnbsSJ06dWLv4pp+DrXktsXrQXAsQmInXpdUa/Y5p0ho0KABx6mparRmWDv8gDRv3jzlmiGYFQyUoQUXKnEIjPeF0eBGncWGEc/yXfTCCy/oG/OYY46hzp07U9OmTTmZ4qFslW63edrOqqeVK1fS22+/zcnO3qCFnH8DTKoXZxM755xzOOl7wxK3Wn5OeDe7C2ybMYOubdYs8DU77rgB9NZbF9GMGfvqAFt+zkPaygsExnPsF++0efNmNXz4cMVHHnXssceqe+65R23YsMF7gxlqzp8/X1111VWKf2V1f0888YRiJpehllxORCBxzXq2aydrlgiOvA8KgXGeGM2WLVvU4MGDVYUKFdSQIUPUqlWrghpgynb5mKWmTZummQ3vltTkyZPV7t27U5aVL2MIyJrJnZBDBOwYzR9//KHuvPNOValSJTVs2DCFX8dc05w5c1SHDh1U27Zt1TvvvJPr4USuf1mzyC1JSRyQOaNZu3atat++vTrhhBPU6tWrIwfWww8/rI9wo0ePVni4hJSSNZO7ICIImDGa119/XVWrVk3dcccdERl36mHgeMACaHX11VcrFianLlRCvpU1KyELnR/TzMxoWNsD9beaPXt2XkwJsppBgwapJk2aKMhySiLJmpXEVY/0nN0ZzXPPPacqVqyo3nvvvUjPItXgLrroIsW2OoptO1JdLtjvZM0KdmnzeWLj0trRLFmyhE455RSaN28e8e4gL5T1yYNkzRh7gR/IEfZvsLbjSW4rHz7LmuXDKpXIMaa2o/npp5/UYYcdpqZMmZLPXFTb2LBFsbrlllvyeh4mg5c1M0FJyuQIgdQ7Gtbc0AcffMCJ35/Me/a7ceNGatGiBa1YsUJbFuf9hNJMQNYsDTDydRQQKO6C8A0n/4LLAHyU4ODohXbu3Kn9mJLrwq2AVc8cnobj07jQxx9/TI2QBO1P+p1DUO7YsYPYQND5yuoVDyEcBFkFblUvXwpnu2ZYLxB8z5LpZwS2YtoH+bmSyK1eYtHk9Uy8lu59oa9ZunkX6PfFj05jx45Vffv29bTDguD1H//4h6pTp06R+nATQJtVqlRR7POk+vXrV+Q6yxZUqVKl4n8DBw7U15nBKLyHuvrCCy9UJ510ktq2bVuRuiYfUKd8+fIFq/L2umbMvBX7jSlmImr//fdX559/flxTB1skrBnW6swzz1SXXnppHGq3eiiUbj3jDRi8KfQ1M4CgkIoU1zq1bt1awafIC/FRS9+csBxOpIkTJypcA8N58MEHtbp8+fLl8SL9+/fX38+YMUPhjz2B9bX7779fsXNmvNzpp5+uRo0aFf9s86Z3796KdzQ2VfKmrNc1g98YXEmgVYQrCf+aqunTp+t5jxw5UvGRU7/HutWtW1eNHz9ef3arhwLp1lNXtvhXyGtmAUMhFB1XJPAVwjXAY5pvXE87uLPPPpvatGlTrO6RHJIA1/bgNKnQZPHuhSpzjBoQuzFwjrPZ+kjVqlUr6tGjB9WoUUNf+/LLLzntx1rC0QmEV0Zdv7f9B2/yRYaR5WzbzmX5bNbsxx9/pHvvvZeA+6233qqPpi/+mZYFnvh9+vTRU8O6nXXWWcSqc/3ZrZ7betriVKhrZotDIZQvwmhwloZ8JtV5PJvJQhjr0N0cLJZ/EeOxZRBWAjKGSzgSNuQyuO7QGWecQbyF1syHnSjpo48+4vQcHMjWAx111FFa7uShaqSrZLNmDzzwQHxue3FQMsSfwcMNuQzaTYz/A3ndu0gQx5SuHq65rSeu21ChrpkNBgVTNnFfBrP1Ll26JH5l/X7MmDHa6TK5IgdS0r5SDJxq3Lix+vrrr+NFvv/+e71lr127tt6+v/LKK/FrCD2BOvhjzVH8e9s3zKQUMzLbapEv78eaYZJYDw5Ipr799lu1ePFijfdLL70Un/+jjz6qv2Ohevy75HrOBbf1dMqYvBbqmpnMvcDKFD06Bck9EbQKuxJmHPo4xDduvDtE3OvZsydnRlxKzGxo6tSp+hq7E+iyJ598staIoMx6RIcT8h0Bln1x1sjr9PEJYTxB2OU4VLZsWf0Wx95ESqznfJ9uPZ3r8lryEChydILaGeEYgyCc87H9voyjXbOQj1joWKwb9N+ds8zD9gXEv6jEYSD0dhwR9iAb+DcSyXsgzCuTWt1Dszmv4sea4bjDOxl9pMWEnNi/WxNS5MC8AOFXWaMYn3NyvfiFP98kr2fy9UyfC3XNMs27EK8XYTSQkaxZs0af0YOc7CGHHFLkhk3sC8JNxJ0F4UY+/PDDdQDzI444Qstx+FiVWNz4Pczzm3HYykKjbNcMci/YF7E2UEMDOyfsZNhbnyCMdwg7ST7yaoE+vktVD/GekylxPZOvZfpcqGuWad6FeL0Io8EvFmL8OkI/LxPGjcXq0CJVsSsBAwPByAtaJlan6s+PcAIyjtCn30Mo/NprrxHbbOjPLC/SsYL1B/5XpkwZHaTc+WzzikDp7dq1s6mSF2WzWTNoiBCDGRol4A6tEtvSaOH8xZyygF1Q4hi8yvlxR4wYoT+nq8dhOshtPeONGb4p1DUznH5hFUsWOkH4ymrN5K+NPiOGLwSujJCOwMe/cLoeDL9gkIdYMd26dVNz586Nt3f00Ufra127dlXwS3rrrbfi12A0NmDAAMUyGm3nAeMxZljx66ZvCt34y+uawf4Ga5X459jOIMQG8GZzAx12A7YzDrnVc1tPp77Ja6GvmQkGBVSmuK9Ttubs6dgwazV0mpVEASPKYvfDAc05P3xFfT1VfQgnoXKFJ7YXKnRz9qDWDFhDLgZBMHaTJmSynibtFPqamWBQQGWK+zphclhkcarMr2WWNcuv9Sphox1fREbjTB4Gde+//35czex8n2+vUI/z9p84kHpck5JvczAdr6yZKVJSLhcIpMzqBstg2LmceuqpVL9+feIMA7kYW9Z9sr+ODniFh7DQSdas0Fc4z+fnJnBCWEieXl6G8rzyyisVq8ZLZChPWTO3u1qu5QCB4sLgZL4JFSOM6GbOnEmsUUi+HLnPDCKnk66rBcwQZHL4g8iNMegByZoFjbC0b4lAahlNYiNwsnv66ae1iwAnj0u8FLn3sGQ97rjjtHMmh7YskUwGiyJrFrlbs8QPKKUwOBkVPLxwAUDogBNPPFF79iaXyfXnSZMm6ZCdbI+jXRf89kDP9fxs+5c1s0VMygeKgM15DQZ0khLXBrHcl5U1y/0ayAhU8Qh7JqAgIyRCO3IMXzV06FDFLgQm1XwrA6tVPs6pY489VnH8HDV58mSFxHFC6RGQNUuPjVwJHAFvjMYZFvu8qOHDh+uc13joEbuWrXydy76/IsQoGFvVqlW1OwNcHhBmUsgcAVkzc6ykpG8IZNY6mZzbYHY+a9YsHfrh+eef18ZxEEiyb5P2mG7YsKF1AjcOnqSdLd9++21asGABQZOC8AUcN1iHmUCbQt4RkDXzjp3UtEYgtQuCdTMJFXADL1u2LM4c4MoAXyaEieSdiA4PicBIENYivQcfg7QfE14RlgB/CDeAWCQI69CyZUvtsc0Ol/FYwgndyVsfEJA18wFEacINAf8ZTare4BSJvEpgIAgxACYCJ0mElChdurRmOmA8YEQITI5YtbVq1UrVlHwXEgKyZiEBXTK6CYfRlAwsZZaCgCCQBoHMBntpKsrXgoAgIAgYI2BksGfcmhQUBAQBQSAFAv8Pu2fiME7me2IAAAAASUVORK5CYII=\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"execution_count": 300,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import networkx as nx\n",
"import matplotlib.pyplot as plt\n",
"from IPython.display import Image\n",
"\n",
"def pydot_plot_graph(Graph,higlight_node=None):\n",
"\n",
" G = Graph.copy()\n",
" \n",
" Colors=['red','blue','black','orange','green','yellow']\n",
" Contigs=[]\n",
" \n",
" for edge in G.edges(data=True): \n",
" contig= edge[2]['Contig']\n",
" if contig not in Contigs:\n",
" Contigs.append(contig)\n",
" #edge[2]['label'] = contig\n",
" edge[2]['color'] = Colors[Contigs.index(contig)]\n",
" \n",
" \n",
" print(\"Legend contigs: \", dict(zip( Colors,Contigs )))\n",
" \n",
"\n",
" \n",
" p=nx.drawing.nx_pydot.to_pydot(G)\n",
" if higlight_node is not None:\n",
" node_of_interest= p.get_node(str(higlight_node))[0]\n",
" node_of_interest.set_color('red')\n",
" node_of_interest.set_shape('box')\n",
" \n",
" p.write_png('multi.png')\n",
" \n",
" \n",
"pydot_plot_graph(Gneighborhood, GeneID_of_interest)\n",
"Image(filename='multi.png')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "honey-provision",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "shaped-choice",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 286,
"id": "prescription-canvas",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "incorporated-technician",
"metadata": {},
"source": [
"## Contigs and Bins\n",
"\n",
"Needs the file \"genomes/clustering/all_contigs2bins.tsv.gz\" \n",
"To re-create it run\n",
"`atlas run None \"genomes/clustering/all_contigs2bins.tsv.gz\"`"
]
},
{
"cell_type": "code",
"execution_count": 334,
"id": "sitting-opposition",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Bin</th>\n",
" </tr>\n",
" <tr>\n",
" <th>0</th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>sample2_0</th>\n",
" <td>sample2_maxbin_2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sample2_1</th>\n",
" <td>sample2_maxbin_2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sample2_2</th>\n",
" <td>sample2_maxbin_2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sample2_3</th>\n",
" <td>sample2_maxbin_2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sample2_4</th>\n",
" <td>sample2_maxbin_2</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Bin\n",
"0 \n",
"sample2_0 sample2_maxbin_2\n",
"sample2_1 sample2_maxbin_2\n",
"sample2_2 sample2_maxbin_2\n",
"sample2_3 sample2_maxbin_2\n",
"sample2_4 sample2_maxbin_2"
]
},
"execution_count": 334,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"contigs= pd.read_table(files['contigs2bins'],index_col=0,header=None)\n",
"contigs.columns= ['Bin']\n",
"contigs.head()"
]
},
{
"cell_type": "code",
"execution_count": 316,
"id": "accessible-keyboard",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>MAG</th>\n",
" <th>Domain</th>\n",
" <th>phylum</th>\n",
" <th>class</th>\n",
" <th>order</th>\n",
" <th>family</th>\n",
" <th>genus</th>\n",
" <th>species</th>\n",
" <th>Label</th>\n",
" </tr>\n",
" <tr>\n",
" <th>genome</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>sample1_maxbin_2</th>\n",
" <td>MAG1</td>\n",
" <td>Bacteria</td>\n",
" <td>Firmicutes</td>\n",
" <td>Bacilli</td>\n",
" <td>Mycoplasmatales</td>\n",
" <td>Metamycoplasmataceae</td>\n",
" <td>Mesomycoplasma</td>\n",
" <td>Mesomycoplasma hyorhinis</td>\n",
" <td>Mesomycoplasma hyorhinis</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sample2_maxbin_3</th>\n",
" <td>MAG1</td>\n",
" <td>Bacteria</td>\n",
" <td>Firmicutes</td>\n",
" <td>Bacilli</td>\n",
" <td>Mycoplasmatales</td>\n",
" <td>Metamycoplasmataceae</td>\n",
" <td>Mesomycoplasma</td>\n",
" <td>Mesomycoplasma hyorhinis</td>\n",
" <td>Mesomycoplasma hyorhinis</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sample2_maxbin_2</th>\n",
" <td>MAG2</td>\n",
" <td>Bacteria</td>\n",
" <td>Firmicutes</td>\n",
" <td>Bacilli</td>\n",
" <td>Lactobacillales</td>\n",
" <td>Streptococcaceae</td>\n",
" <td>Streptococcus</td>\n",
" <td>Streptococcus thermophilus</td>\n",
" <td>Streptococcus thermophilus</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sample2_maxbin_1</th>\n",
" <td>MAG3</td>\n",
" <td>Bacteria</td>\n",
" <td>Firmicutes</td>\n",
" <td>Bacilli</td>\n",
" <td>Mycoplasmatales</td>\n",
" <td>Mycoplasmoidaceae</td>\n",
" <td>Ureaplasma</td>\n",
" <td>Ureaplasma urealyticum</td>\n",
" <td>Ureaplasma urealyticum</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" MAG Domain phylum class order \\\n",
"genome \n",
"sample1_maxbin_2 MAG1 Bacteria Firmicutes Bacilli Mycoplasmatales \n",
"sample2_maxbin_3 MAG1 Bacteria Firmicutes Bacilli Mycoplasmatales \n",
"sample2_maxbin_2 MAG2 Bacteria Firmicutes Bacilli Lactobacillales \n",
"sample2_maxbin_1 MAG3 Bacteria Firmicutes Bacilli Mycoplasmatales \n",
"\n",
" family genus \\\n",
"genome \n",
"sample1_maxbin_2 Metamycoplasmataceae Mesomycoplasma \n",
"sample2_maxbin_3 Metamycoplasmataceae Mesomycoplasma \n",
"sample2_maxbin_2 Streptococcaceae Streptococcus \n",
"sample2_maxbin_1 Mycoplasmoidaceae Ureaplasma \n",
"\n",
" species Label \n",
"genome \n",
"sample1_maxbin_2 Mesomycoplasma hyorhinis Mesomycoplasma hyorhinis \n",
"sample2_maxbin_3 Mesomycoplasma hyorhinis Mesomycoplasma hyorhinis \n",
"sample2_maxbin_2 Streptococcus thermophilus Streptococcus thermophilus \n",
"sample2_maxbin_1 Ureaplasma urealyticum Ureaplasma urealyticum "
]
},
"execution_count": 316,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bins = pd.read_table(files['allbins2genome'],index_col=0)\n",
"Tax= pd.read_table(files['taxonomy'],index_col=0)\n",
"\n",
"# create label also for unnabed species\n",
"Tax['Label'] = Tax.ffill(axis=1).species\n",
"Tax.loc[Tax.species.isnull(),'Label']+= ' '+ Tax.index[Tax.species.isnull()]\n",
"\n",
"bins= bins.join(Tax,on='MAG')\n",
"\n",
"\n",
"\n",
"bins.head()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ahead-defendant",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 336,
"id": "japanese-reputation",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Bin</th>\n",
" <th>MAG</th>\n",
" <th>Domain</th>\n",
" <th>phylum</th>\n",
" <th>class</th>\n",
" <th>order</th>\n",
" <th>family</th>\n",
" <th>genus</th>\n",
" <th>species</th>\n",
" <th>Label</th>\n",
" </tr>\n",
" <tr>\n",
" <th>0</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>sample1_3</th>\n",
" <td>sample1_maxbin_2</td>\n",
" <td>MAG1</td>\n",
" <td>Bacteria</td>\n",
" <td>Firmicutes</td>\n",
" <td>Bacilli</td>\n",
" <td>Mycoplasmatales</td>\n",
" <td>Metamycoplasmataceae</td>\n",
" <td>Mesomycoplasma</td>\n",
" <td>Mesomycoplasma hyorhinis</td>\n",
" <td>Mesomycoplasma hyorhinis</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sample2_36</th>\n",
" <td>sample2_maxbin_3</td>\n",
" <td>MAG1</td>\n",
" <td>Bacteria</td>\n",
" <td>Firmicutes</td>\n",
" <td>Bacilli</td>\n",
" <td>Mycoplasmatales</td>\n",
" <td>Metamycoplasmataceae</td>\n",
" <td>Mesomycoplasma</td>\n",
" <td>Mesomycoplasma hyorhinis</td>\n",
" <td>Mesomycoplasma hyorhinis</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Bin MAG Domain phylum class \\\n",
"0 \n",
"sample1_3 sample1_maxbin_2 MAG1 Bacteria Firmicutes Bacilli \n",
"sample2_36 sample2_maxbin_3 MAG1 Bacteria Firmicutes Bacilli \n",
"\n",
" order family genus \\\n",
"0 \n",
"sample1_3 Mycoplasmatales Metamycoplasmataceae Mesomycoplasma \n",
"sample2_36 Mycoplasmatales Metamycoplasmataceae Mesomycoplasma \n",
"\n",
" species Label \n",
"0 \n",
"sample1_3 Mesomycoplasma hyorhinis Mesomycoplasma hyorhinis \n",
"sample2_36 Mesomycoplasma hyorhinis Mesomycoplasma hyorhinis "
]
},
"execution_count": 336,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Contig and Bin information for gene of interest\n",
"\n",
"contigs_from_gene_of_interest= orfs.query(\"GeneId == @GeneID_of_interest\").Contig.unique()\n",
"\n",
"\n",
"contigs.loc[contigs_from_gene_of_interest].join(bins,on='Bin')\n",
"\n",
"#bins.loc[contigs.loc[contigs_from_gene_of_interest].values]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "economic-andorra",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
},
"widgets": {
"application/vnd.jupyter.widget-state+json": {
"state": {},
"version_major": 2,
"version_minor": 0
}
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment