{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import pandas as pd\n",
"import numpy as np\n",
"import re\n",
"import os\n",
"import math\n",
"from multiprocessing import Pool\n",
"from tqdm import tqdm\n",
"from scipy import stats\n",
"## init\n",
"mySpecie='Homo_sapiens'\n",
"#prealigned_dir='/cellar/users/btsui/all_seq_snp/Homo_sapiens_all_merged_snp.TCGA.prealigned.pickle'\n",
"targetted_align_dir='/cellar/users/btsui/all_seq_snp/Homo_sapiens_all_merged_snp.TCGA.pickle'\n",
"manifest_dir='/cellar/users/btsui/Project/METAMAP/notebook/RapMapTest/XGS_WGS/./tcga_lgg_wgs_bams.df.wxs_rnaseq.pickle'\n",
"\n",
"targetted_df=pd.read_pickle(targetted_align_dir).loc[\"TCGA\"]"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"all_UUIDs=targetted_df.index.get_level_values('Run_digits').unique()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 112,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/cellar/users/btsui/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:20: SettingWithCopyWarning: \n",
"A value is trying to be set on a copy of a slice from a DataFrame.\n",
"Try using .loc[row_indexer,col_indexer] = value instead\n",
"\n",
"See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n"
]
}
],
"source": [
"manifest_df=pd.read_pickle(manifest_dir)\n",
"\n",
"manifest_df['processed']=manifest_df.file_id.isin(all_UUIDs)\n",
"\n",
"uuid_barcode_mapDf=pd.read_csv('/cellar/users/andreabc/GDC_barcodes/uuid_barcode_map.txt',sep='\\t').set_index('file_id')\n",
"\n",
"manifest_df['sample_barcode']=uuid_barcode_mapDf.loc[manifest_df.file_id]['sample_barcode'].values\n",
"\n",
"m_data_category=manifest_df.data_category=='Raw Sequencing Data'\n",
"m_experimental_strategy=manifest_df['experimental_strategy'].isin(['RNA-Seq','WXS'])\n",
"\n",
"manifest_df_sub=manifest_df[manifest_df['processed']&m_data_category&m_experimental_strategy]\n",
"\n",
"tmpVC=manifest_df_sub['sample_barcode'].value_counts()\n",
"\n",
"#identify data with both RNAseq and WXS\n",
"with_both=tmpVC.index[tmpVC==2]\n",
"\n",
"\n",
"manifest_df_sub['is_tumor']=manifest_df_sub['sample_barcode'].str.contains('TCGA-\\w+-\\w+-01')\n",
"manifest_df_w_RNA_WXS=manifest_df_sub[manifest_df_sub.sample_barcode.isin(with_both)&manifest_df_sub['is_tumor']]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"inVcfDir='/data/cellardata/users/btsui/dbsnp/Homo_sapiens/All_20170710.f1_byte2_not_00.vcf.gz' \n",
"vcfDf=pd.read_csv(inVcfDir,sep='\\t',header=None)\n",
"vcfDf.columns=['Chr','Pos','RsId','RefBase','AltBase','','','Annot']\n",
"vcfDf['Chr']=vcfDf['Chr'].astype(np.str)"
]
},
{
"cell_type": "code",
"execution_count": 270,
"metadata": {},
"outputs": [],
"source": [
"#queryChr,querySite,refBase='17',7673803,'G','A' #TP53\n",
"#8\t142877758\n",
"#6, 31356729\n",
"#G\tA\n",
"#1\t237591774\n",
"#\n",
"queryChr,querySite,refBase='1',237591774,'C',\n",
"m_chrom=(targetted_df.index.get_level_values('Chr')==queryChr)"
]
},
{
"cell_type": "code",
"execution_count": 271,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Chr | \n",
" Pos | \n",
" RsId | \n",
" RefBase | \n",
" AltBase | \n",
" | \n",
" | \n",
" Annot | \n",
"
\n",
" \n",
" \n",
" \n",
" 28989 | \n",
" 1 | \n",
" 237591774 | \n",
" rs75901196 | \n",
" C | \n",
" A | \n",
" . | \n",
" . | \n",
" RS=75901196;RSPOS=237591774;dbSNPBuildID=131;S... | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Chr Pos RsId RefBase AltBase \\\n",
"28989 1 237591774 rs75901196 C A . . \n",
"\n",
" Annot \n",
"28989 RS=75901196;RSPOS=237591774;dbSNPBuildID=131;S... "
]
},
"execution_count": 271,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vcfDf[vcfDf.Pos==querySite]"
]
},
{
"cell_type": "code",
"execution_count": 272,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" | \n",
" features | \n",
" ReadDepth | \n",
" AverageBaseQuality | \n",
"
\n",
" \n",
" Run_digits | \n",
" Chr | \n",
" Pos | \n",
" base | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" 08ce1dd9-3167-4fe4-8619-724727a32a36 | \n",
" 1 | \n",
" 14727 | \n",
" A | \n",
" 18 | \n",
" 28 | \n",
"
\n",
" \n",
" G | \n",
" 280 | \n",
" 32 | \n",
"
\n",
" \n",
" 630825 | \n",
" T | \n",
" 11 | \n",
" 32 | \n",
"
\n",
" \n",
" 630833 | \n",
" C | \n",
" 9 | \n",
" 34 | \n",
"
\n",
" \n",
" 833068 | \n",
" G | \n",
" 2 | \n",
" 34 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"features ReadDepth \\\n",
"Run_digits Chr Pos base \n",
"08ce1dd9-3167-4fe4-8619-724727a32a36 1 14727 A 18 \n",
" G 280 \n",
" 630825 T 11 \n",
" 630833 C 9 \n",
" 833068 G 2 \n",
"\n",
"features AverageBaseQuality \n",
"Run_digits Chr Pos base \n",
"08ce1dd9-3167-4fe4-8619-724727a32a36 1 14727 A 28 \n",
" G 32 \n",
" 630825 T 32 \n",
" 630833 C 34 \n",
" 833068 G 34 "
]
},
"execution_count": 272,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"targetted_df.head()"
]
},
{
"cell_type": "code",
"execution_count": 273,
"metadata": {},
"outputs": [],
"source": [
"m_site=(targetted_df.index.get_level_values('Pos')==querySite)"
]
},
{
"cell_type": "code",
"execution_count": 274,
"metadata": {},
"outputs": [],
"source": [
"hitDf=targetted_df[m_site&m_chrom]"
]
},
{
"cell_type": "code",
"execution_count": 275,
"metadata": {},
"outputs": [],
"source": [
"hitDfResetDf=hitDf.reset_index()"
]
},
{
"cell_type": "code",
"execution_count": 276,
"metadata": {},
"outputs": [],
"source": [
"hitDfResetDf['is_ref']=hitDfResetDf.base==refBase"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 277,
"metadata": {},
"outputs": [],
"source": [
"uuidToExperimentS=manifest_df_w_RNA_WXS.set_index('file_id')['experimental_strategy']"
]
},
{
"cell_type": "code",
"execution_count": 278,
"metadata": {},
"outputs": [],
"source": [
"uuidToBarcodeS=manifest_df_w_RNA_WXS.set_index('file_id')['sample_barcode']"
]
},
{
"cell_type": "code",
"execution_count": 279,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/cellar/users/btsui/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:1: FutureWarning: \n",
"Passing list-likes to .loc or [] with any missing label will raise\n",
"KeyError in the future, you can use .reindex() as an alternative.\n",
"\n",
"See the documentation here:\n",
"https://pandas.pydata.org/pandas-docs/stable/indexing.html#deprecate-loc-reindex-listlike\n",
" \"\"\"Entry point for launching an IPython kernel.\n",
"/cellar/users/btsui/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:2: FutureWarning: \n",
"Passing list-likes to .loc or [] with any missing label will raise\n",
"KeyError in the future, you can use .reindex() as an alternative.\n",
"\n",
"See the documentation here:\n",
"https://pandas.pydata.org/pandas-docs/stable/indexing.html#deprecate-loc-reindex-listlike\n",
" \n"
]
}
],
"source": [
"hitDfResetDf['Strategy']=uuidToExperimentS.loc[hitDfResetDf['Run_digits']].values\n",
"hitDfResetDf['Barcode']=uuidToBarcodeS.loc[hitDfResetDf['Run_digits']].values"
]
},
{
"cell_type": "code",
"execution_count": 280,
"metadata": {},
"outputs": [],
"source": [
"validHitResetDf=hitDfResetDf.dropna()"
]
},
{
"cell_type": "code",
"execution_count": 281,
"metadata": {},
"outputs": [],
"source": [
"unstackDf=validHitResetDf.set_index(['Barcode','base','Strategy'])['ReadDepth'].unstack().unstack().fillna(0)#.set_index(['Strategy'])"
]
},
{
"cell_type": "code",
"execution_count": 282,
"metadata": {},
"outputs": [],
"source": [
"import seaborn as sns"
]
},
{
"cell_type": "code",
"execution_count": 283,
"metadata": {},
"outputs": [],
"source": [
"#unstackDf"
]
},
{
"cell_type": "code",
"execution_count": 284,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'ref at 1:237591774'"
]
},
"execution_count": 284,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"'ref at {}:{}'.format(queryChr,querySite)"
]
},
{
"cell_type": "code",
"execution_count": 292,
"metadata": {},
"outputs": [],
"source": [
"altBase='A'"
]
},
{
"cell_type": "code",
"execution_count": 293,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/cellar/users/btsui/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
" warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n",
"/cellar/users/btsui/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
" warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 293,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X94VOWd9/HP1wCSIiyiQCXhhyiNKD81Cj5BRXcptGpL0Vop3UJbRZ9V649tWtza1lp70V167bZ7bd0+qK37A9GqFG1tRWoFVqpgEBQsRIpSCPiUUM0jYJAQvs8fc5KGZCbMhDn3ZCbv13XNxcw997nP956Q85lz5uSMubsAAAjphFwXAADoeggfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4LrluoBWuNwCgHxnuS4gH7DnAwAIrrPt+XTIw2t2tGn77IQhOagEAJAO9nwAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQXLdcFxCXh9fsSNr+2QlDAlcCAGiNPR8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRXs3/mkkuzvf/jbHwAIiz0fAEBwXW7PJ5lUV0NIhr0kADh+7PkAAIJjzydDXDMOAI4f4ZMlmRy6y0eZhGsmJ3V05kOenbm2TGTyhqmrvrniRKTwzN1zXUMzM3tG0qkdWPRUSXuzXE5nxVwLU1eaq1TY893r7tNyXURn16nCp6PMrMrdy3NdRwjMtTB1pblKXW++aIsTDgAAwRE+AIDgCiV8Fua6gICYa2HqSnOVut580UpBfOYDAMgvhbLnAwDII4QPACA4wgcAEBzhAwAIrlOFz7Rp01wSN27cuOXzLW0Fus1LS6cKn717C/VqGwDQVlfe5nWq8AEAdA2EDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAguG5xDWxmZZIebdE0XNI33f0Hca0TCGnp+l1asKxau+vqNahvsSqnlmn6+JJclwXkhdjCx92rJY2TJDMrkrRL0s/jWh8Q0tL1u3Tnko2qb2iUJO2qq9edSzZKEgEEpCHUYbe/lrTN3f8YaH1ArBYsq24Onib1DY1asKw6RxUB+SVU+FwraXGyJ8xsrplVmVlVbW1toHKA47O7rj6jdqAJ27yE2MPHzHpI+oSkx5I97+4L3b3c3cv79+8fdzlAVgzqW5xRO9CEbV5CiD2fj0l6xd3/FGBdQBCVU8tU3L3oqLbi7kWqnFqWo4qA/BLbCQctzFSKQ25Avmo6qYCz3YCOiTV8zOxDkqZIuiHO9QC5MH18CWEDdFCs4ePu70s6Jc51AADyD1c4AAAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABBcrOFjZn3N7HEz22Jmm83swjjXBwDID91iHv+Hkp5x96vNrIekD8W8PgBAHogtfMysj6SLJc2RJHc/JOlQXOsDAOSPOA+7DZdUK+mnZrbezB4ws16tO5nZXDOrMrOq2traGMsBgNxjm5cQZ/h0k3SupH939/GSDkia17qTuy9093J3L+/fv3+M5QBA7rHNS4gzfGok1bj7mujx40qEEQCgi4stfNz9/0raaWZlUdNfS/p9XOsDAOSPuM92u0XSouhMtzclfSHm9QEA8kCs4ePuGySVx7kOAED+4QoHAIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAILrFufgZrZd0j5JjZIOu3t5nOsDAOSHWMMncqm77w2wHuSZpet3acGyau2uq9egvsWqnFqm6eNLUrYDKBwhwgdoY+n6XbpzyUbVNzRKknbV1evOJRtV9cd39MS6XW3aJRFAQAGJ+zMfl/Ssma0zs7kxrwt5ZMGy6uaAaVLf0KjFa3YmbV+wrDpkeQBiFveeT4W77zazAZKWm9kWd1/VskMUSnMlaciQITGXg85id1190vZG94z6A/mGbV5CrHs+7r47+nePpJ9LuiBJn4XuXu7u5f3794+zHHQig/oWJ20vMsuoP5Bv2OYlxBY+ZtbLzHo33Zf0UUmb4lof8kvl1DIVdy86qq24e5FmThictL1yalnI8gDELM7DbgMl/dwS72S7SXrY3Z+JcX3II00nDyQ7q618aD/OdgMKnHmKY+y5UF5e7lVVVbkuAwCOR/Jjx0kU6DYvrflzhQMAQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgusW9AjMrklQlaZe7XxH3+rLlrqUbtXjNTjW6q8hMMycM1r3TR8e+3qXrd2nBsmrtrqvXoL7FqpxapunjS3IyTrZqQf7gZ45QYg8fSbdK2iypT4B1ZcVdSzfqv1/a0fy40b35cZwBtHT9Lt25ZKPqGxolSbvq6nXnko2SlNEGIBvjZKsW5A9+5ggp1sNuZlYq6XJJD8S5nmxbvGZnRu3ZsmBZdfMvfpP6hkYtWFYdfJxs1YL8wc8cIcX9mc8PJH1V0pFUHcxsrplVmVlVbW1tzOWkp9E9o/Zs2V1Xn1F7nONkqxbkD37mYXTGbV4uxBY+ZnaFpD3uvq69fu6+0N3L3b28f//+cZWTkSKzjNqzZVDf4oza4xwnW7Ugf/AzD6MzbvNyIc49nwpJnzCz7ZIekXSZmf13jOvLmpkTBmfUni2VU8tU3L3oqLbi7kWqnFoWfJxs1YL8wc8cIaV1woGZDZU0wt1/Y2bFkrq5+772lnH3OyXdGS0/WdJX3P1zx1lvEE0nFYQ+263pQ93jPdsoG+NkqxbkD37mCMn8GJ9jmNn1kuZK6ufuZ5jZCEk/dve/Tnslfwmfdk+1Li8v96qqqnSHBYDOKO3j8wW6zUtr/ukcdrtJiUNo70mSu2+VNCCTStx9RT79jQ8AIF7phM8H7n6o6YGZdZMU72lfAICClk74rDSzf5BUbGZTJD0m6RfxlgUAKGTphM88SbWSNkq6QdKvJN0VZ1EAgMJ2zLPd3P2IpPujGwAAxy1l+JjZRrXz2Y67j4mlIgBAwWtvz4ez0wAAsUgZPu7+x5CFAAC6jvYOu+1T8sNuJsndPW++IgEA0Lm0t+fTO2QhAICuo709n37tLeju72S/HABAV9DeCQfrlDjsluw6PS5peCwVAQAKXnuH3U4PWQgAoOtI9ysVTpY0QlLPpjZ3XxVXUQCAwnbM8DGz6yTdKqlU0gZJEyW9KOmyeEsDABSqdK7tdquk8yX90d0vlTReiWu9AQDQIemEz0F3PyhJZnaiu2+RxPfqAgA6LJ3PfGrMrK+kpZKWm9m7knbHWxYAoJClc1XrT0V37zaz5yX9laRnYq0KAFDQ0jrbrYm7r4yrEABA15HOZz4AAGQV4QMACO6Y4WNm/5hOGwAA6Upnz2dKkraPZbsQAEDX0d5Vrf+3pL+TNNzMXmvxVG9Jq+MuDABQuNo72+1hSb+WNF/SvBbt+/g6BQDA8WgvfIokvSfpptZPmFk/AggA0FHpfJ+P1PY7ffg+HwBAh/F9PgCQI+8cOJTrEnKG7/MBAATH9/kAAILj+3wAAMGlc9jtoLsfNLPm7/Mxs2N+n4+Z9ZS0StKJ0Xoed/dvHWe9HbZ0/S4tWFat3XX1GtS3WJVTyzR9fEnW+sddT6bjZNo+6/4XtXrbX05grDijnxZdf2FO5tTZFOq8gFwyd2+/g9nPJX1B0m1KHGp7V1J3d//4MZYzSb3cfb+ZdZf0gqRb3f2lVMuUl5d7VVVVhlM4tqXrd+nOJRtV39DY3FbcvUjzZ4xOuhHJtH/c9WQ6zlXnleiJdbvSbi89uae27jnQZvxMAiju1yxXCnVeiFXrs4NTGj5yjL+5+bVjd8wvac3/mIfd3P1T7l7n7ndL+oakByVNT2M5d/f90cPu0a39pIvJgmXVR208JKm+oVELllVnpX/c9WQ6zuI1OzNqTxY8ko7aE+poLdl6zXKlUOcF5Fqs3+djZkVK/L3QmZJ+5O5rkvSZK2muJA0ZMiST4dO2u64+1va468l0nMYUe7Op2rMh7tcsVwp1Xsidltu8Uz/cdfeeY/1KBXdvdPdxSpwpd4GZjUrSZ6G7l7t7ef/+/WOpY1Df4ljb464n03GKLPleb6r2bIj7NcuVQp0XcqflNq933365Lidngnyfj7vXSVohaVqI9bVWObVMxd2Ljmor7l6kyqnJz5vItH/c9WQ6zswJgzNqHzGgV9LxK85I/xcj7tcsVwp1XkCuZXTYLRNm1l9Sg7vXmVmxpL+RlJPvAWr6YDjdM5Yy7R93PR0Zp3xov4zaj/dst7hfs1wp1HkBuXbMs906PLDZGEn/ocQFSk+Q9DN3v6e9ZeI62w0AAuJstzTEtufj7q8p8QepAAAcJchnPgCAtvr16pHrEnKG8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AGAHHnnwKFcl5AzhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCiy18zGywmT1vZpvN7HUzuzWudQEA8ku3GMc+LOnv3f0VM+staZ2ZLXf338e4zpSWrt+lBcuqtbuuXoP6Fqtyapmmjy9J2T7ln1do654DzcuPGNBLy++YnPE4cdefrfnGXT+QL/hdCMPcPcyKzJ6U9G/uvjxVn/Lycq+qqsr6upeu36U7l2xUfUNjc1tx9yJddV6Jnli3q017n55F+tO+tuffD+zdQ+8dbEx7nPkzRmflP22q+lONn+l8464fyBeZ/q6lYOl2HD5yjL+5+bUMq+z00pp/kM98zGyYpPGS1oRYX2sLllUf9Z9JkuobGrV4zc6k7cmCR5L+tO9QRuMsWFadhepT159q/EznG3f9QL7I9HcNHRd7+JjZSZKekHSbu7+X5Pm5ZlZlZlW1tbWx1LC7rj5pe2OW9vpSjZNqvZlKNU6m7anqjLt+IF9k+jvVES23efvq3snauPkm1vAxs+5KBM8id1+SrI+7L3T3cncv79+/fyx1DOpbnLS9yNLeO25XqnFSrTdTqcbJtD1VnXHXD+SLTH+nOqLlNq93335ZGzffxHm2m0l6UNJmd//nuNaTjsqpZSruXnRUW3H3Is2cMDhp+8DePZKOM7B3j4zGqZxaloXqU9efavxM5xt3/UC+yPR3DR0X59luFZL+VtJGM9sQtf2Du/8qxnUm1fRBYbIzWMqH9svK2W6pxom7/mzNN876gXyR6e/a8erXK/kb3a4g2Nlu6YjrbDcACCjt4/kFus3rPGe7AQDQEuEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABEf4AACCI3wAAMERPgCA4AgfAEBwhA8AILhuuS4AALqqdw4c0sNrdrRp/+yEITmoJiz2fAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAEF9vf+ZjZTyRdIWmPu4+KYx3D5j3dpq2bSYf96Md/mH+5zrzz6aTts+5/Uau3vdPcXnFGPy26/sKU7WO+9Yze+6Cxub3PiUV67dvTUvZfun6XFiyr1u66eg3qW6zKqWWaPr4k4/6ZtqeSaf9cyIca8wWvJTorc/dj9+rIwGYXS9ov6T/TDZ/y8nKvqqpKa/xkwZMtfU4sOipgmpikZK9WqvYRA3qp5t2Dqm/4y1jF3YtUenJPbd1zIO3+V51XoifW7Uq7ff6M0Uk3MEvX79KdSzam3T8X8qHGfMFrmTOWbsfhI8f4vQ/9sk17nv+RaVrzj+2wm7uvkvTOMTt2QsmCR0oeMO21b91z4KhffEmqb2hMGjzt9V+8ZmdG7QuWVScdf8Gy6oz650I+1JgveC3RmeX8Mx8zm2tmVWZWVVtbm+tyOqXGFHunqdp319VnpT0X8qHGfMFr2Tm13Obtq8vL9+dZkfPwcfeF7l7u7uX9+/fPdTmdUpEl34tN1T6ob3FW2nMhH2rMF7yWnVPLbV7vvv1yXU7O5Dx8OqM+JxYlbU91IDNV+4gBvVTc/eixirsXacSAXhn1nzlhcEbtlVPLko5fObUso/65kA815gtey86vX68e+uyEIW1uXUHehs/2712etL2btX28/XuXp2yvOOPodx4VZ/TTa9+elrT9re9d3iaY+pxYpLdSjLP8jsmaP2O0SvoWyySV9C3W/BmjtfyOyRn1v3f66IzaU32YPH18SUb9cyEfaswXvJbozOI8222xpMmSTpX0J0nfcvcH21smk7PdAKCTSvtstwLd5qU1/9j+zsfdZ8Y1NgAgv+XtYTcAQP4ifAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAgOMIHABAc4QMACI7wAQAER/gAAIIjfAAAwRE+AIDgCB8AQHCEDwAguG65LgAAuqp3DhzSw2t25LqM4/bZCUMyXoY9HwBAcIQPACA4wgcAEBzhAwAIjhMOELuGhgbV1NTo4MGDuS4FyJqePXuqtLRU3bt3z3UpeYnwQexqamrUu3dvDRs2TGaW63KA4+bu+vOf/6yamhqdfvrpuS4nL3HYDbE7ePCgTjnlFIIHBcPMdMopp7A3fxwIHwRB8KDQ8H/6+BA+AIDgYv3Mx8ymSfqhpCJJD7j797I5/rB5T6fVb/v3Lk/atyPtE767XH/ad6i5bWDvHlrz9Slaun6XFiyr1u66eg3qW6zKqWWaPr5Es+5/Uau3vdPcv+KMflp0/YW6a+lGLV6zU43uKjLTzAmDde/00SnnkGr8VO3oup555hndeuutamxs1HXXXad58+Yl7fezn/1Md999t8xMY8eO1cMPPyxJ+upXv6qnn35aR44c0ZQpU/TDH/4wZ+/yv//976uyslK1tbU69dRT2zxfVFSk0aMTvzdDhgzRU089lfbY1dXV+sxnPtP8+M0339Q999yj2267TY899pjuvvtubd68WWvXrlV5efnxTwZHMXePZ2CzIklvSJoiqUbSy5JmuvvvUy1TXl7uVVVVaY2fbvCE0OfEIjUckeobGpvbirsXqfTkntq650Cb/gN79zgqwJp8buKQpAG0dP0u3blkY5vxrzqvRE+s29Wmff6M0Z0qgDZv3qyRI0fmuoysO3z4sLp1y/77t8bGRhUVFXV42Y985CNavny5SktLdf7552vx4sU6++yzj+q3detWXXPNNfrtb3+rk08+WXv27NGAAQP0u9/9TpWVlVq1apUkadKkSZo/f74mT558vNPK2M6dO3Xddddpy5YtWrduXdLwOemkk7R///7jXldjY6NKSkq0Zs0aDR06VJs3b9YJJ5ygG264Qd///vdThk+K/9tpJ/XwkWP83od+eRyVdw6tLq+T1vzjPOx2gaQ/uPub7n5I0iOSPhnj+nLmvQ8ajwoAKRFEyYJHUtLgkaTFa3YmbV+wrDrp+IvX7EzavmBZdbqldwnbt2/XWWedpdmzZ2vMmDG6+uqr9f7770uS1q1bp0suuUTnnXeepk6dqrfffluSdP/99+v888/X2LFjddVVVzX3nzNnju644w5deuml+trXvqaVK1dq3LhxGjdunMaPH699+/bJ3VVZWalRo0Zp9OjRevTRRyVJK1as0OTJk3X11VfrrLPO0qxZs9T05m/YsGG65557NGnSJD322GMdnuvatWt15plnavjw4erRo4euvfZaPfnkk2363X///brpppt08sknS5IGDBggKfE5xsGDB3Xo0CF98MEHamho0MCBAyVJ1113nZK9OZwzZ45uvPFGXXTRRfrIRz6iX/4yOxvT22+/Xf/0T//Uob2uVD/XVJ577jmdccYZGjp0qCRp5MiRKisr61DdSE+ch91KJLXcmtZImtC6k5nNlTRXSuw2d2WNKfZCd9fVZ6V/V1ZdXa0HH3xQFRUV+uIXv6j77rtPt956q2655RY9+eST6t+/vx599FF9/etf109+8hPNmDFD119/vSTprrvu0oMPPqhbbrlFkvTGG2/oN7/5jYqKinTllVfqRz/6kSoqKrR//3717NlTS5Ys0YYNG/Tqq69q7969Ov/883XxxRdLktavX6/XX39dgwYNUkVFhVavXq1JkyZJSvzdyAsvvNCm9kWLFmnBggVt2s8880w9/vjjR7Xt2rVLgwcPbn5cWlqqNWvWtFn2jTfekCRVVFSosbFRd999t6ZNm6YLL7xQl156qU477TS5u26++ebmd/YPPPBAytd3+/btWrlypbZt26ZLL71Uf/jDH9SzZ8/m5/ft26eLLroo6bIPP/xwmz2zp556SiUlJRo7dmzKdUqJMynLy8vVrVs3zZs3T9OnT1dDQ0PKn2sqjzzyiGbOnNnuurKl9TavIxflLARxhk+ytytttpbuvlDSQilx2C3Gejq9ohTv8Ab1LdauJIFSZJY0gAb1Lc56bflu8ODBqqiokCR97nOf07/+679q2rRp2rRpk6ZMmSIpcejltNNOkyRt2rRJd911l+rq6rR//35NnTq1eaxPf/rTzYfFKioqdMcdd2jWrFmaMWOGSktL9cILL2jmzJkqKirSwIEDdckll+jll19Wnz59dMEFF6i0tFSSNG7cOG3fvr05fFp+/tDSrFmzNGvWrLTmmewwerI9h8OHD2vr1q1asWKFampqdNFFF2nTpk3au3evNm/erJqaGknSlClTtGrVqubwTOWaa67RCSecoBEjRmj48OHasmWLxo0b1/x87969tWHDhrTm8P777+u73/2unn322WP23bFjhwYNGqQ333xTl112mUaPHq36+vqUP9dkDh06pKeeekrz589Pq77jxTYvIc7wqZE0uMXjUkm7Y1xfzmTrM5+ZEwa3aZOkyqllGX3mUzmVwwWttd4Am5ncXeecc45efPHFNv3nzJmjpUuXauzYsXrooYe0YsWK5ud69erVfH/evHm6/PLL9atf/UoTJ07Ub37zm6QB0OTEE09svl9UVKTDhw8nHbelTPZ8SktLtXPnXw441NTUaNCgQW2WLS0t1cSJE9W9e3edfvrpKisraw6jiRMn6qSTTpIkfexjH9NLL710zPBJ9vq2lMmez7Zt2/TWW2817/XU1NTo3HPP1dq1a/XhD3/4qGWb5jZ8+HBNnjxZ69evV1lZWdKf686dO3XllVdKkm688UbdeOONkqRf//rXOvfcc5sPLyKMOD/zeVnSCDM73cx6SLpWUvqnohzD9u9dftx9O9I+sHePo9oG9u6h1749TfNnjFZJ32KZpJK+xZo/Y7SW3zFZFWf0O6p/xRn9tObrU/S5iUOa93SKzFKebCBJ08eXJB3/3umjk7Z3ppMNOosdO3Y0b4wWL16sSZMmqaysTLW1tc3tDQ0Nev311yUlNpannXaaGhoatGjRopTjbtu2TaNHj9bXvvY1lZeXa8uWLbr44ov16KOPqrGxUbW1tVq1apUuuOCCDtc+a9Ysbdiwoc2tdfBI0vnnn6+tW7fqrbfe0qFDh/TII4/oE5/4RJt+06dP1/PPPy9J2rt3r9544w0NHz5cQ4YM0cqVK3X48GE1NDRo5cqVzYfdPv/5z2vt2rVJa3zsscd05MgRbdu2TW+++Wabz0ua9nyS3Vofchs9erT27Nmj7du3a/v27SotLdUrr7zSJnjeffddffDBB81zWL16tc4+++yUP9fBgwc3r7MpeKTE/4dQh9zQgrvHdpP0cSXOeNsm6evH6n/eeec5Cs/vf//7nK7/rbfe8pEjR/oNN9zgo0eP9hkzZviBAwfc3X39+vV+0UUX+ZgxY/zss8/2hQsXurv7fffd58OGDfNLLrnEb775Zp89e7a7u8+ePdsfe+yx5rFvvvlmP+ecc3zMmDF+7bXX+sGDB/3IkSP+la98xc855xwfNWqUP/LII+7u/vzzz/vll1/evOxNN93kP/3pT93dfejQoV5bW5uV+T799NM+YsQIHz58uN97773N7d/4xjf8ySefdHf3I0eO+O233+4jR470UaNG+eLFi93d/fDhwz537lw/66yzfOTIkX777bc3Lz927FjfsWNHm/XNnj3bb7vtNp80aZKPGDHCf/GLX2RlHk1avjYvv/yyf+lLX3J399WrV/uoUaN8zJgxPmrUKH/ggQeal0n1c23twIED3q9fP6+rqzuqfcmSJV5SUuI9evTwAQMG+Ec/+tGky6f4v532NrJAt3lpzT22U607IpNTrZE/cn2q9fbt23XFFVdo06ZNOash37333nv60pe+lPRMvDlz5uiKK67Q1VdfnYPKcut4T7Uu0G1ezk+1BlAg+vTpc1yngAOtcVVrFLwXtPm/AAAKa0lEQVRhw4ax1xOjhx56KNclIA+x54MgOtPhXSAb+D99fAgfxK5nz57685//zC8rCoZH3+fT8g9pkRkOuyF2paWlqqmpUW1tba5LAbKm6ZtM0TGED2LX9IeMANCEw24AgOAIHwBAcIQPACC4TnWFAzOrlfTHDix6qqS9WS6ns2KuhakrzVUq7Pnudfdp6XQ0s2fS7VtoOlX4dJSZVbl7l/ieW+ZamLrSXKWuN1+0xWE3AEBwhA8AILhCCZ+FuS4gIOZamLrSXKWuN1+0UhCf+QAA8kuh7PkAAPII4QMACC6vw8fMpplZtZn9wczm5bqebDOzn5jZHjPb1KKtn5ktN7Ot0b8n57LGbDGzwWb2vJltNrPXzezWqL3g5mtmPc1srZm9Gs3121H76Wa2Jprro2bWI9e1ZouZFZnZejP7ZfS4YOeK9ORt+JhZkaQfSfqYpLMlzTSzs3NbVdY9JKn1H6DNk/Scu4+Q9Fz0uBAclvT37j5S0kRJN0U/z0Kc7weSLnP3sZLGSZpmZhMl/aOkf4nm+q6kL+Wwxmy7VdLmFo8Lea5IQ96Gj6QLJP3B3d9090OSHpH0yRzXlFXuvkrSO62aPynpP6L7/yFpetCiYuLub7v7K9H9fUpsqEpUgPP1hP3Rw+7RzSVdJunxqL0g5ipJZlYq6XJJD0SPTQU6V6Qvn8OnRNLOFo9rorZCN9Dd35YSG2xJA3JcT9aZ2TBJ4yWtUYHONzoMtUHSHknLJW2TVOfuh6MuhfT/+QeSvirpSPT4FBXuXJGmfA4fS9LGeeN5zsxOkvSEpNvc/b1c1xMXd29093GSSpXYix+ZrFvYqrLPzK6QtMfd17VsTtI17+eKzOTzl8nVSBrc4nGppN05qiWkP5nZae7+tpmdpsQ754JgZt2VCJ5F7r4kai7Y+UqSu9eZ2QolPufqa2bdoj2CQvn/XCHpE2b2cUk9JfVRYk+oEOeKDOTzns/LkkZEZ830kHStpKdyXFMIT0maHd2fLenJHNaSNdHnAA9K2uzu/9ziqYKbr5n1N7O+0f1iSX+jxGdcz0u6OupWEHN19zvdvdTdhynxO/pbd5+lApwrMpPXVziI3k39QFKRpJ+4+3dzXFJWmdliSZOVuPz8nyR9S9JSST+TNETSDkmfdvfWJyXkHTObJOl/JG3UXz4b+AclPvcpqPma2RglPmQvUuIN4M/c/R4zG67EiTP9JK2X9Dl3/yB3lWaXmU2W9BV3v6LQ54pjy+vwAQDkp3w+7AYAyFOEDwAgOMIHABAc4QMACI7wAQAER/jkETO7KLoK8obo70PSWWa7mZ0a3d+fRv9j9olTrtZvZreZ2YcyXGaOmQ1K8dyno5/VETMrT9En6ZW8o+e+Y2avRT/rZ5vWY2aVUdsGM9tkZo1m1i96bruZbYyeq2ox1lgzezF67hdm1idqPyVa/34z+7cW/Xu3WMcGM9trZj+InvuXFu1vmFldqzn1MbNdLccDknJ3bp3kpsRlR05o5/kfS/pChmNul3RqdH9/Gv2P2aedZYuy8Bp0eP3Hud7m1ymDZVZIKk/x3EhJZcfoc5qkc6P7vSW9Iens6HGfFv2+LOnHSZa/Uok/2mx3Dkr8QfYl0f0vSvpOdL+XpEmSbpT0b+3Mc52ki5O036LE39e1bPuhpIfbG48bN3dnzyfXzGxY9M73PkmvSBpsZh+N3qm+YmaPmdlJZnadpGskfdPMFiUZZ6mZrYveQc9NY72VZvZy9O7628fRZ7+Z3WNmayRdaGbnmdnKqJZl0SVxZGbXR2O9amZPNO1lRFeoeDF67jvt1Pv5qI5Xzey/orahZvZc1P6cmQ2J2h8ys6tbLLs/+neyma0ws8fNbIuZLbKEL0saJOl5M3s+ybq/GdW3ycwWRstcLalc0iJLsifq7pvdvbqdH4E89ZW85Udf166Xkl/7bKakxe2tI1ImaVV0f7mkq6J1HHD3FyQdTLWgmY1Q4mKu/3Os9ZvZeZIGSno2jZrQ1eU6/br6TdIwJf6if2L0+FQlNhS9osdfk/TN6P5Dkq5OMU6/6N9iSZsknRI93q5Wez6SPippoaI9LUm/VPTONp0+rdbrkq6J7neX9DtJ/aPHn1H0zripnuj+vZJuie4/Jenz0f2blGTPR9I5kqpbzKNprr+QNDu6/0VJS5O9Ti3mNFnS/1PiWmInSHpR0qTWr1Oq1za6/1+Srozur1CKvZoW/Y/qo0TI/SrF/4MdOnqP57tKXLl9U9Nr2uK5DynxdRsta3tLiTcw6yTNbdH+O0mfjO7fIWlfq7HmKMWeiqRvSvp+kvahkt5WtLcbvZ4rlLjeYsrxuHFrurHn0zn80d1fiu5PVOLL8VZb4pL7s5X4RT+WL5vZq5JeUmIDMKKdvh+NbuuV2FidlaR/On0kqVGJi4FKiXfYoyQtj2q/S4kNvSSNMrP/MbONkmYpEShS4sKTTe+e/ytFvZdJetzd90qS/+XyOhcqcYinadlJqafcbK2717j7EUkblNjoH8ullvjWzY1RLecca4FU3H23u3+8ZZuluJK3u3/d3QdLWiTp5lZDXSlptR99qaEKdz9XiS9YvMnMLo7avxg9XqfE4b1DGZR8rZLvXV2rxM+kMXr8d0qE6s4kfYE28vmq1oXkQIv7Jmm5u89Md2FLXDPrbyRd6O7vW+IqyT3bW0TSfHf/P8fZR5IOttgAmaTX3f3CJP0ekjTd3V81szlK7IU0OdY1niyNPi3HOazoZBozM0ktv6K55fXDGnWM3wEz6ynpPiX2Xnaa2d1q/7XNiCW/kndrD0t6Wolr+zVpEwruvjv6d4+Z/VyJr2pY5e5blHgjITP7iBJf7JZObWMldfOjvw6h5fpvavH4QkkXmdnfSTpJUg8z2+/uhfDNs4gBez6dz0uSKszsTEkysw9FG4z2/JWkd6PgOUuJvaf2LJP0xegdt8ysxMxaf0lbOn1aq5bU38wujJbpbmZNewm9Jb0dbWxntVhmtRIbMrVqb+k5SdeY2SnRuP2i9t+1WvaF6P52SedF9z+pxOHAY9kX1dhaU9DsjV6Lq1s8l2qZtETBmOxK3k2ftTT5hKQtLZ77K0mXqMWVoM2sl5n1brqvRNhsih4PiP49QYm90R+nWWLSz5TMrEzSyUoctpQkufssdx/iiatXf0XSfxI8aA/h08m4e60Sx8wXm9lrSoTRWcdY7BlJ3aL+34mWaW8dzyrxbvrF6FDS42q1EU2nT5JxDymxcf7H6BDgBkn/K3r6G0pcoXq5WmxIJd2qxCGhl5UI0WTjvq7E5x8ro3GbNtRflvSFaN5/G40lSfdLusTM1kqaoKP3LFNZKOnXrU84cPe6aLyNSlxR/OUWTz8k6cfJTjgws0+ZWY0SewRPm9myqH2Qmf0q6lYR1X2Z/eX05aZDct+LTnB4TYkgubXF8J+S9Ky7t5zXQEkvRK/PWklPu/sz0XMzzewNJV733ZJ+2qLO7Uq8nnPMrMbMzm4x5jVKfshtpqRH3J2rEqPDuKo1ACA49nwAAMERPgCA4AgfAEBwhA8AIDjCBwAQHOEDAAiO8AEABPf/AWxk6+S0Sie/AAAAAElFTkSuQmCC\n",
"text/plain": [
"