{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Notebook 16: Simulating RADseq data\n",
"#### Eaton et al. 2015"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"## Requirements\n",
"## - Python 2.7\n",
"## - pyrad v.3.1.0 (http://github.com/dereneaton/pyrad)\n",
"## - simrrls v.0.0.7 (http://github.com/dereneaton/simrrls)"
]
},
{
"cell_type": "code",
"execution_count": 150,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import itertools\n",
"import ete2\n",
"import numpy as np\n",
"import toyplot\n",
"from collections import OrderedDict, Counter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generate trees for simulations\n",
"Make a balanced tree with 64 tips and a tree length = 10"
]
},
{
"cell_type": "code",
"execution_count": 151,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for n, node in enumerate(Tbal.get_leaves()):\n",
" node.name = \"t\"+str(n) "
]
},
{
"cell_type": "code",
"execution_count": 152,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"((((((t0:1.66667,t1:1.66667)1:1.66667,(t2:1.66667,t3:1.66667)1:1.66667)1:1.66667,((t4:1.66667,t5:1.66667)1:1.66667,(t6:1.66667,t7:1.66667)1:1.66667)1:1.66667)1:1.66667,(((t8:1.66667,t9:1.66667)1:1.66667,(t10:1.66667,t11:1.66667)1:1.66667)1:1.66667,((t12:1.66667,t13:1.66667)1:1.66667,(t14:1.66667,t15:1.66667)1:1.66667)1:1.66667)1:1.66667)1:1.66667,((((t16:1.66667,t17:1.66667)1:1.66667,(t18:1.66667,t19:1.66667)1:1.66667)1:1.66667,((t20:1.66667,t21:1.66667)1:1.66667,(t22:1.66667,t23:1.66667)1:1.66667)1:1.66667)1:1.66667,(((t24:1.66667,t25:1.66667)1:1.66667,(t26:1.66667,t27:1.66667)1:1.66667)1:1.66667,((t28:1.66667,t29:1.66667)1:1.66667,(t30:1.66667,t31:1.66667)1:1.66667)1:1.66667)1:1.66667)1:1.66667)1:1.66667,(((((t32:1.66667,t33:1.66667)1:1.66667,(t34:1.66667,t35:1.66667)1:1.66667)1:1.66667,((t36:1.66667,t37:1.66667)1:1.66667,(t38:1.66667,t39:1.66667)1:1.66667)1:1.66667)1:1.66667,(((t40:1.66667,t41:1.66667)1:1.66667,(t42:1.66667,t43:1.66667)1:1.66667)1:1.66667,((t44:1.66667,t45:1.66667)1:1.66667,(t46:1.66667,t47:1.66667)1:1.66667)1:1.66667)1:1.66667)1:1.66667,((((t48:1.66667,t49:1.66667)1:1.66667,(t50:1.66667,t51:1.66667)1:1.66667)1:1.66667,((t52:1.66667,t53:1.66667)1:1.66667,(t54:1.66667,t55:1.66667)1:1.66667)1:1.66667)1:1.66667,(((t56:1.66667,t57:1.66667)1:1.66667,(t58:1.66667,t59:1.66667)1:1.66667)1:1.66667,((t60:1.66667,t61:1.66667)1:1.66667,(t62:1.66667,t63:1.66667)1:1.66667)1:1.66667)1:1.66667)1:1.66667)1:1.66667);\n"
]
}
],
"source": [
"print Tbal.write()"
]
},
{
"cell_type": "code",
"execution_count": 153,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"## base tree\n",
"Tbal = ete2.Tree()\n",
"\n",
"## branch lengths\n",
"bls = 10/6.\n",
"\n",
"## first nodes\n",
"n1 = Tbal.add_child(dist=bls)\n",
"n2 = Tbal.add_child(dist=bls)\n",
"\n",
"## make balanced tree\n",
"while len(Tbal.get_leaves()) < 64:\n",
" thisrep = Tbal.get_descendants()\n",
" for node in thisrep:\n",
" if len(node.get_children()) < 1:\n",
" node.add_child(dist=bls)\n",
" node.add_child(dist=bls)\n",
"\n",
"## set leaf names\n",
"for n, node in enumerate(Tbal.get_leaves()):\n",
" node.name = \"t\"+str(n) \n",
" \n",
"## Save newick string to file\n",
"Tbal.write(outfile=\"Tbal.tre\", format=3)"
]
},
{
"cell_type": "code",
"execution_count": 154,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"((((((t0:1.66667,t1:1.66667)NoName:1.66667,(t2:1.66667,t3:1.66667)NoName:1.66667)NoName:1.66667,((t4:1.66667,t5:1.66667)NoName:1.66667,(t6:1.66667,t7:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667,(((t8:1.66667,t9:1.66667)NoName:1.66667,(t10:1.66667,t11:1.66667)NoName:1.66667)NoName:1.66667,((t12:1.66667,t13:1.66667)NoName:1.66667,(t14:1.66667,t15:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667,((((t16:1.66667,t17:1.66667)NoName:1.66667,(t18:1.66667,t19:1.66667)NoName:1.66667)NoName:1.66667,((t20:1.66667,t21:1.66667)NoName:1.66667,(t22:1.66667,t23:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667,(((t24:1.66667,t25:1.66667)NoName:1.66667,(t26:1.66667,t27:1.66667)NoName:1.66667)NoName:1.66667,((t28:1.66667,t29:1.66667)NoName:1.66667,(t30:1.66667,t31:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667,(((((t32:1.66667,t33:1.66667)NoName:1.66667,(t34:1.66667,t35:1.66667)NoName:1.66667)NoName:1.66667,((t36:1.66667,t37:1.66667)NoName:1.66667,(t38:1.66667,t39:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667,(((t40:1.66667,t41:1.66667)NoName:1.66667,(t42:1.66667,t43:1.66667)NoName:1.66667)NoName:1.66667,((t44:1.66667,t45:1.66667)NoName:1.66667,(t46:1.66667,t47:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667,((((t48:1.66667,t49:1.66667)NoName:1.66667,(t50:1.66667,t51:1.66667)NoName:1.66667)NoName:1.66667,((t52:1.66667,t53:1.66667)NoName:1.66667,(t54:1.66667,t55:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667,(((t56:1.66667,t57:1.66667)NoName:1.66667,(t58:1.66667,t59:1.66667)NoName:1.66667)NoName:1.66667,((t60:1.66667,t61:1.66667)NoName:1.66667,(t62:1.66667,t63:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667)NoName:1.66667);"
]
}
],
"source": [
"cat Tbal.tre"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Make an imbalanced tree of the same total tree length with 64 tips"
]
},
{
"cell_type": "code",
"execution_count": 155,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"## base tree\n",
"Timb = ete2.Tree()\n",
"\n",
"## scale branches to match balanced treelength\n",
"brlen = (bls*6.)/63\n",
"\n",
"## first nodes\n",
"n1 = Timb.add_child(dist=brlen)\n",
"n2 = Timb.add_child(dist=brlen)\n",
"\n",
"while len(Timb.get_leaves()) < 64: \n",
" ## extend others\n",
" for tip in Timb.get_leaves()[:-1]:\n",
" tip.dist += brlen\n",
" ## extend the last node\n",
" Timb.get_leaves()[-1].add_child(dist=brlen)\n",
" Timb.get_leaves()[-1].add_sister(dist=brlen)\n",
"\n",
"## set leaf names\n",
"for n, node in enumerate(Timb.get_leaves()):\n",
" node.name = \"t\"+str(n) \n",
" \n",
"## write to file\n",
"Timb.write(outfile=\"Timb.tre\", format=3)"
]
},
{
"cell_type": "code",
"execution_count": 156,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(t0:10,(t1:9.84127,(t2:9.68254,(t3:9.52381,(t4:9.36508,(t5:9.20635,(t6:9.04762,(t7:8.88889,(t8:8.73016,(t9:8.57143,(t10:8.4127,(t11:8.25397,(t12:8.09524,(t13:7.93651,(t14:7.77778,(t15:7.61905,(t16:7.46032,(t17:7.30159,(t18:7.14286,(t19:6.98413,(t20:6.8254,(t21:6.66667,(t22:6.50794,(t23:6.34921,(t24:6.19048,(t25:6.03175,(t26:5.87302,(t27:5.71429,(t28:5.55556,(t29:5.39683,(t30:5.2381,(t31:5.07937,(t32:4.92063,(t33:4.7619,(t34:4.60317,(t35:4.44444,(t36:4.28571,(t37:4.12698,(t38:3.96825,(t39:3.80952,(t40:3.65079,(t41:3.49206,(t42:3.33333,(t43:3.1746,(t44:3.01587,(t45:2.85714,(t46:2.69841,(t47:2.53968,(t48:2.38095,(t49:2.22222,(t50:2.06349,(t51:1.90476,(t52:1.74603,(t53:1.5873,(t54:1.42857,(t55:1.26984,(t56:1.11111,(t57:0.952381,(t58:0.793651,(t59:0.634921,(t60:0.47619,(t61:0.31746,(t62:0.15873,t63:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873)NoName:0.15873);"
]
}
],
"source": [
"cat Timb.tre"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Check tree lengths"
]
},
{
"cell_type": "code",
"execution_count": 157,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"set([10.0]) treelength\n",
"64 tips\n",
"set([9.999999999999998]) treelength\n",
"64 tips\n"
]
}
],
"source": [
"print set([i.get_distance(Tbal) for i in Tbal]), 'treelength'\n",
"print len(Tbal), 'tips'\n",
"\n",
"print set([i.get_distance(Timb) for i in Timb]), 'treelength'\n",
"print len(Timb), 'tips'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get topology from empirical study of Viburnum with 64 tips\n",
"Make tree ultrametric using the penalized likelihood"
]
},
{
"cell_type": "code",
"execution_count": 158,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The rpy2.ipython extension is already loaded. To reload it, use:\n",
" %reload_ext rpy2.ipython\n"
]
}
],
"source": [
"%load_ext rpy2.ipython"
]
},
{
"cell_type": "code",
"execution_count": 159,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAJYCAMAAAB7MkC6AAADAFBMVEUAAAABAQECAgIDAwMEBAQF\nBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcY\nGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKior\nKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+\nPj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBR\nUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2Nk\nZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29wcHBxcXFycnJzc3N0dHR1dXV2dnZ3\nd3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGCgoKDg4OEhISFhYWGhoaHh4eIiIiJiYmK\nioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJyd\nnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+w\nsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLD\nw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW\n1tbX19fY2NjZ2dna2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp\n6enq6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8\n/Pz9/f3+/v7////isF19AAAgAElEQVR4nO19CXwUVfZ1K4IsyqKiIzBuuIBBEjSyJRAgISgkYMJu\nUBYBBcEgqFERwRFZdNSIooDBMCOMAgPKKAyLMOCCOiDhD6OogPkUcBCUZVBiCN33q3ffra27qnqr\n7rzurvMbi+rqqk53zqTr3PvuPdcFDoSCq7rfgAM9HEIEg0OIYHAIEQwOIYLBIUQwOIQIBocQweAQ\nIhgcQgSDQ4hgcAgRDA4hgsEhRDA4hAgGhxDB4BAiGBxCBINDiGBwCBEMDiGCwSFEMDiECAaHEMHg\nECIYHEIEg0OIYHAIEQwOIYLBIUQwOIQIBocQweAQIhgcQgSDQ4hgcAgRDA4hgsEhRDA4hAgGhxDB\n4BAiGBxCBINDiGBwCBEMDiGCwSFEMDiECAaHEMHgECIYHEIEg0OIYHAIEQwOIYLBIUQwOIQIBocQ\nweAQIhgcQgSDQ4hgcAgRDA4hgsEhRDA4hAgGhxDB4BAiGBxCBINDiGBwCBEMDiGCwSFEMDiECIYE\nJ6TgGP6z+pj+8P2DAQ7VKgOYUzozqUPSRjw4p9T4NXZ0SfnMvneU4IT04v/0OaM/3C0XYGLT3wGG\n7xz4LexpgweH7zR+jS7/2zXWvneUKISc6gsfPV2ZnvrNY1nZVTPTO/zkmZ7d8deKfJh2R/sf1jR5\nhT19JnVI8lPZzwO4u9wBR4ZkSJe1O5PqBuCEtDvDzt2b2a/lIbyq8tb+qXuk423cxSX2vdF4IGTd\n4/P9Ytqf4OX3ypJPLHkail8b7nnttfefgJeWbZ/y2zyYtv7w0I+flp7efbtn2JoTdwDsvSff/fjK\nuwAqO/yWLvHTif0UaZ+du7zAs+g13CvreHbBmwCe7KtanLXvw8QDIde4/OOcRTD0QGkJDPgFYGyP\ngT3fHfZfeG3VwhXf3ld0/c/ril9eLT395pvQ07PtSYDlL9+3Z/Cq51fDjns+HwOw8WHpLiPt47mT\nP4GVC3Bv0SIokr7FZizzjP3Yvg8TD4TkJR3zi2Ff/3SFp3AP9D1yYsKob91z3bnfV2SfGL8v+9jR\n9jB7yyj29KQvz+bA6ysBJn80PfeTpza0gkcWz58Hx7ocBA/bx3N7bXbnfY97D34JvaV7T9bxqoGH\n7Psw8UBIfmv/57yedlcv6OWG9a07bN5xc59FsKZ19geQ4y7qOuLq/cWj8encs3smwv3lADknF3WF\n/KOjU4ZUjW2Z1m2bdD3bx3OzumctBtzLOQs9pWfe6tp+uY0fJiKEHK0RwJeIjWgYiQ9hjCPDI/wD\nIkLIQVePoijihqsDeE8UamQd9H7C9wjh9/RbG6WNaFAB6e/BwqfhxxHsYL9x8EWPo88lpaU95pnZ\n7dq32TH2zIYerb+E0qTkobApNXlliL83hogQcsg1PxIva4ZAvrJ8Qo1AsE2KLxpU7Gi2GDr/+Mol\nT7FDWcc/yTkBBd9Iu9Of9Zxo7gFgz5zKPbN2AgyokCRZysmfuwX/oxTEAyEZjf3/Fd19+SKMHfJe\nn7rfdXzKZiWQqMxTg4pTumgEoGQBI2Rk8dyv+sPZl96Rjnhu3dTnN4C2ks79JVN6vEuKUtgz65+H\nsieOtRgwrOrkBvh0QhgfJh4IaVerkV807s9jh8nLpxa1+annr0ogUTZZDSr00QjAuM8lQg4M2vzM\npI0AI/dLRw42yb0ToKJRRsYj782Qf770zLaxnhGlXy3z3LMO4L3McERXPBASCNYVY+xQuvyDRwcU\nzH1TDSRKl6v7+mgEoIv0x9Bg6sadEzpJ30yZ0h8DrB5wNu0wfDFa2i1+FeDrN9hp0jPuu3LTd0m7\nf95SVTjmZDhvNFEImb0FY4fCvV/csLSge6UaSBTuVff10Qh42kmbBlme8iYvAFSxbyh4ZiXMnw4L\nJS7gHzmeM91ZVpE9M/PH4709mUcq+/5WFKYGThRCikfz2MH93WWVBTM0gQSFIrivj0bgu4HSpsE8\nOF5bCu+/uo+9zIByONmi6oFPpV33iLYZS9kx9sz6Tnn74a0bc7b+fnlycjj3dHsJebjBNYgrXZPt\nfNmoQhazYzdLm+du7NhqFcjKt/y29m12S0dR2OIGxa7A6fc+tfojclxT7XzZqILE7KqaP0kPCr6G\n3X0V5XvHPlg1iYQtbrjYFTj9/sBF/N+IfWWtvr5f/9Aw7BZJ1npY3r0yPXXmHHh1CUu/nx6YOUGv\ndUnMHuqXId3HoaMbnn9aUb4p71Wd+RUAhS1uUOyKnH6POCEPupJvCQ2tb5Vk7UaWdy9LPrGtqKLn\nP1j6/d2HoESvdbmYded/lyX9U3FxRvro3xTlu7P/1dMYTVzYShsUuyKn3yNOyCuuIyFeibIW8+6l\nJXB0yJxVmH4/PqbjKr3W5WL22cX7RkkXfXE/XkvK9/Nf4WS7fZKqYsIWNyh2RU6/C0wIl7Us7164\nBzxZeR5Mvxe5v+2n17pczOZkp7SQfskL5+G1pHxH/RN+6VYBgMIWNyh2RU6/C0xIzp3HegLm3Xu5\nYXWLD3j6fWZGl+1M62Z1l7Ru7thyVcymN++QtPGBT1jesLRGs6GS8p3Uom69W/4NgMKWq9vpHZjY\n1affw8wtJgwhcjkDwjDR6MnVPuK1DZhRZClDYGlFWPy0/hKD1wk3txgZQn5w3RKiGPKDNq7DIb41\nKmfAHOKaJouUfax1qMyrKux69/0owpjwks7ntQ0sb4gpQ4mvNuC5/0M8hb9SeuqcJqXsap6uXM9L\nHsLNLUaGkOMNW4QohvzgCpeyeHHzVcFd+rCaQzw8VN3HWoeyySUvwfjXUYQx4QUg1zaM3A+UMjx4\nUcaNt7nxFH518onDQ/Fqnq7cwUsews0tRoaQiEHzlXXRlUH9bbVfrOYQKdGI+yi2SpffcRoe3YYi\njAkvALm2ATOK8OctUgj0JHhyduMpPE1ZAuuK8WqerqSSh3Bzi7FLSPO7grqSyhlQbFGiEfdRbBXu\n7fXzieanUYQx4QVAtQ0sb4gpQ4Bn3oHTKcfxFJ6m3AOzt+DVPF1JJQ/h5hYThhBtDlFONLJ9FFuX\nVd10/ahMSYQ1m4vCC6D9+RfWkyKSLS0kbVuv/uvSC9zQPO365KMTal903WNF11xec44k1oqHZLfO\nfrZHi4srM5uldDq7qUHyyrBziwlDiBVUfaVIMY3IolSV1dqtu+3vWx8JW19xxBYhU1295ZvCBVlB\nXcnE08apMPtdnrTCglBUWKivaB23Ip+rKK3IolSV5drtr1vhhcVh6yuO2CLklbpXX0Oo2S6oK5l4\n2jv6ZO4unrTCglDUSKivaB13+xRUUTqRRakqy7VbgJ4tq8LWVxyxRYgG+q+sitoXWK+pNxvWcVVF\n7z9toqQVFoSiRkJ9Reu4C1egitKJLEpVWa7d/gKeEV+Era844oSQk65u1lUnHX+TxFOnvkBJKywI\nRY2E+orWccfvQxWlE1mUqrJcu+106mzvY2HrK464IeR569NRPKVulRdosSAUFdayJElfyRrMjdku\nACogZYuzlKqyXLst6XTbqvD1FUfMEnJ5kvYPYIJrmv9LDt8p70W8IDR0xCwhjbzqe4f5Of/jPluz\nyuUHG4rZ1ruRzRu9WJjuVWma6w7lzQaBmCVED79fWVBw1OeQn+pST7bBsVzfY/YibgiZbt0fsq7p\n85qmtrOGjWys242FJnSkPDU3+XBlHk8Ks/gEQ5ZIf5I4IeS4v46FWr1w7Zya2sqMGtnYEy+w0ISO\nvDPYs3hu2WRMCmN8giFLpD9JnBACDz5n3WM45UlcO6emNsNGNvYEhiZ0ZOoWePfV0uXqajyGLJH+\nIPFCiD8seHeUpqnNsJGNPYGhCR3ps+nsHT8U7lVX4zFkifQbTRRCxpZrm9pY1NFnRtPnbmqW9IpS\nOjo9qcNmDE3oSKsLGjTbmLO+YfKOnrCj1aXN2zyfNKpTvbUA948DKGIRCa64syadw307rYV9t7Xe\nyTcb2qbsCPWNJgohvmCZXZ308hZdmPHFRK70gKoWi/rPB+h6UyXc/itoVtxz9/84A/r8/MGDuJGv\nCQkJQwjTVWc0qd6K/Mq8LU3n9/+fmejiGV9M5EoPeNXikcy1j4O788PvQgaAuuK+n4Xo34yD7Y/g\nRr4mJESJkKp/bq9e7GS6arcm1bt9Stnk3/pBL1PRJWd8MZFLVYtFa78sgH3Dducf7s8+lbzi/mrv\n3qOqSpbByrm4ka8JCVEiZEyE+m8DxjkFkq7SpnoXrihd/vm0k4NMRRfP+GIiV3qAVYtHmmZ3S4cV\nxdD+rWfYp5JX3B9+tKpo45TPoXAXbuia0BAlQl5w/X1DteJWpqu0qd7x+wr3Lnj3o5mmootnfDGR\nC1S1WLQGPDfClM1QfN37AOqK+7itZwvKH113vCfghq4JDVEiJIwKN3uA2VxtqjfHneMeW/7yP2XR\n5eMewDO+LJErXb4rpV3nfx/p7AFodabPMThy3gGgZPCNOVth5y0dl8B/bsr+mm/omtCQKIRoIHCq\nFxKIENWqjKd6jcHtyjCOwADDsokKY47bk5Ob/Ygbe95nwhBiYlXmBU3sQSW95k1UPPAA2H2vvLED\nMUbI4WeXhYb387FTak2TUrmA1y1bmWE6V1vSi3EElfRaNFFhzAFQmXeSNrYgxgi5LWTZ+wBGGYeH\nKgW8ipUZpnP1Jb1SHEElvRZNVBRzPMm6Q3FjC2KMkMdrhBgYzl6BUca6YqWAV7Eyw3SutqRXjiNY\nSa9FExXGHFCeLz2NG3sQY4Q8cV6IF47fh1GGFGzIBbyKlRmmc7UlvRhHUEmvRRMVxhww7kOgjT1I\nFEJy3BhlFI/GuhJWg6JYmWGNLy/p1cQeGGBIH+kTvJzsAzAcoSYqHnMc7SJ9g+HGJlQ/IZ8G881z\nT6iEwNCfgjtf9uv1ZyDg5ZYVNuwn5NR0A5Fzjykh9wV5dw71vXUP8nzy6/VnIODllhU+7CdkmvEv\n8oDJNbNcpUGo1741lAurhj4zK3DMbTn4lj1qV5uieU3Nsni5tV8DAS+3rPB/h/YT8pZr3X4fTDP9\nC3nRdTyIn6C5h6wM6g+rZurZeW+qnWyK5jUzy+ICOBADAS+3rHBhPyFvu77yfcb8HhIyIVtdf/fv\nDqvgjUXw6E61k03RvGZmWVwA+zcQ8HLLCuKjmCCWCVkbxIWobdVONkXzmpllcQHs30DAyy0riHdk\nAnsI+aVVH9ZEc11diDAhudqvoQVBXIjatuiygjpbcY9p3jnJkubVmmW1ufjaW87uuXcES763qp/W\nbU2tMpj4SOnMpMuvZnprTp2+N7Y9ZwnTvu17A6yrWeeqbl5uWUG8ITPYQ8g61x9YE02DWhBhQtYN\nUG7TY13vBP0+te4BPnlG3qlG7qMG8xGG7+Rii1c32NTB5gN7CNnmeo/9E/mvLA2C+8pSChvcZkUN\n1KnG3UeN5iO0O8PFFq9usKmDzQexS8hmV8/RgeOOwVTYYFrUIHeqofuowXyEyg4ktuRzbOlg80Hs\nEvL5efX9u8OqeI8KG0yLGuRONaxc8JmPwPa+UCqtpXNs6mDzQewSEhyUwgbTogbqVOPuo97zEQ6w\nPRJbXFnZ1MHmg0QhRClsMC1qoE417j7qMx+B7VGiEc+xq4PNBwlCiNiFDVpEiZBZrikmSaacqBCi\nFjZ49aStPqb2rJnOSdAjsp1uUSKkyCLN5NtrFkF496QFPzQhwp1uUSLEHP6+sobelGUH7nxSttHg\nlmXpqXt5MLIgj9U3sP8q8ygmOVWdnW7CE9KsgS3WZ20ekm00uGVZ8gmaz1Y2mdU3sP/KJlNMUq2d\nbsIT0u720N+WBqqNBrcsK5GDkdLlrL6B/Ve6nGKSau10E5GQn0o1VdIt7BGXqo0GtyzbAxSMFO5l\n9Q3sv8K9FJNUa6ebiIRk6e75TWx5h9irhjZla1o33sCUEm9rG53jbpnRaM217btszzlcY40Uk4yo\nO4HFJ9n5t+TD/Z+PgPU3Xtz6aenqw5envz5xUuOURmlZS1lEwzLEWCSB5hy2QURCJp6vKWtoZeOH\n5dDIJJb9VRUS9qvRlnepYeoX9+gfA69YmyEiIUV1NA9suodwGYWqCGUSjcIl+7LKPOkU7FejLWCX\nGqZ1+R7/h3vF2vKGzCA8Ic3qXmMHbhjHZBSqIpRJNAqX7MvK2LwT7FejLWCXGqZ1aQ//4V6x4fym\n/EJ4Qoa0s8WCefBIJqNQFaFMolG4ZF9Wulz608B+Nb4FwC41TOvSHv7DvWLD+lX5g/CE2AQuo1AV\noUyiUbhkX1a4VzoD+9X4FgC71DCty/f4P9wr1v53p0GiEDLzso7byZ+sbw9JJtEoXKa7Zl7Y4cIP\nYOWFySvh6uR6KwBaNaiANtd2HP2058Gm176985aOk1OTV+5sUO+KZr2bX1DnIfvfnBaJQohF5SJv\n0lHmGIE63ZP8YOXDu++N1EK6BtVOSKErxTvN8Yca/i8LFpWsctGj9OdUpqfOnwqz32HqCtfM1TlG\noEz3lP1g6XBl3slILaRrUO2ELGmW6Z0IvCrkimpz7Ox4dh6NPML5nsknmGnsjslqk448xwiU6Z6q\nHywexqacyCyka1DthBhA85U15pygCkbNUX8OPLpT7c8pLQFmGsvUFV8zV+cYgVKiKK+y88OsKSdS\nC+kaCE7IqBrhzPXWoMMOOfXE53vuAWYay9QVr1FU5xiBUqIo+8Hyw6wpJ1IL6RoITshj54f5WjLU\n1BPusVW/1K2SzKImHc0cI1Dac2iVnR9mTTkRW0jXIFEI8YVqGisUEpMQ6Q/kM9U01ggvt0xrM4zm\nsMqdU+gngF5lEURCEmK0LO6Fe3ZLAWIVD0QoLDmIfgKU+40YRCRkYi2lYyqvFh07GkQviBF+pamr\nZWy9nOd7LR0E0qug/HYKRCgsQT+B/Yl4D+mhkavn8EMPhyt7zxvAp67iejnmey0dBM5cktGxzn4K\nROTOKeYnQLnfyEFEQk59oCxPDae/kJmu561tYP2hWzafuorr5ZjvtXQQ2D0SYPJSCkQoLEE/Acr9\nRg4iEqKBfA95OdzqLXnqKq6XY77X0kHgzTkAkz6kQITCEvQToNxv5JAohMhTV7GeF/O9lg4Ck0Z3\nu7bdh280qnuZFIjUuhWOn1+rrKRTywdmNq9Xly0oYne6uuH2TOqxMJAohAQJLq+o89ygmwoX29UN\n2jNpH4aORCFEGWXPdvhIMJ2wMhpz/zPvPDfqpsI+K3WD9kyah2FAcEKGnysPwTWqnj+4KWBPjrny\nKHumotAndqNOWBl1U1HnuUE3FfVQqRvMBKsPw4DghDxUj9coXOIyctC7InDZe84QGmXPVBT6xOqF\nlXE3FXae+3RTsR+NfVbKhtszqcfCgOCEyDD+yurSImCb2DvlUfZMRWE7lV5YGXVTUee5dzcVa0TA\nPit1g/ZM6sNwENOE3B74LEN5lD3uYDuVXlgZdVNR57lPNxVQn5Wy4fZM6rFwkCiEaCF0O1XCEPJx\nH+XLXW2nYn1QRnld1Qkr6ogRQua4PjCQTh1vCfwVDIaCcRjldRUnrOgjRggZaaydGionuGe8bWm0\n9VzT+eiMlZ76ZXrqt5T6/ZIVnRjkdVUnrOgjRgg5+rKRdGqbopzwrj/Z25M7YyWfkP7HQxBpd7JZ\nXld2woo+YoQQY2juIZ+53vC1TdPgvWnYMVVaAtL/MARhu8uN87qqE1b0YSshd9edNWvWna7/2PKa\nAUBHyBrLUxe8ix1ThXvY/+TULys6Mcrrqk5Y0YethJDv9Me2vGYACIKQseVo39DLzVK+cuo3J8ed\nde8rfVa1+ZAKTHi5CRWXfNUhKS3tblJhcISNAyOrUk1aV9Zluh8WTte6rYQgoviV1b2VIrgW+SHE\nENg8VXD06U+Nny74hm25CuONVTzZq03rymXAPi8bKmKakGba2/aLlqd6JmfeUyR3p2OOl5qn2GCw\nU9S6XqY3Jm17ll2JKowaq3iyV5vWJV2mpo7xZTdOhdmrQ/tMMU3Ip39VBNfLrvf5sZODDN1jR8yA\n+96Wu9Mxx0vNU2wwmNy6rk/5VjTKyMg4zFUYb6ySk72atC7pMjV1jC/L6oZDNLuOaUI0UO4hfzWW\nved+DWO/lrvTMcdLzVNsMJjcuq5P+X4xmr0eV2G8sUpO9mrSuqTL1NQxviyrGw7xg8QdIZtchpMj\nep08ft1ZuTsdc7zUPMUGg8mt6/qU78LX2IVchfHGKkr2atO6slOpkjrmXeud+ob6QeKHkGW8AGuV\nMSEr2o3IlLvTuQsptpfzwWDYTiUvuCsp3wdapqUNh0nSy/VdwceBUbJXm9alMmA1dcy71lNDNiiN\nF0LeU7+d/mbw9PHVsOzZ6L2bMOqG44WQ3+56mt+9R7k2GDztHph5xa/Ko9U+Uwafswo5cPhXMPBT\nN2yJeCFEgck9RJft9TXJsgw5aPhXdBB3hLzvyjWyiR3YZB4LNjBeoMlg2lJeq5CDhn9FCXFHyL9q\nXWHk5HB9DgYbGC/QZDBNKa9lyEHDv6KEuCPEBJ9Pw2AD4wWaDKYp5bUMOfjwr2ghUQhZ8C4GGxgv\n0GQwTSmvZcjBh39FC4lCyNhyDDYwXqDJYJpSXiXkGJojhxwN2skhBw7/ItXFF9uVHC9fiadpVDu6\npHwGNLQqdCQKIQFC9Xvo5fWM6vcAJlYPXf63a6w8tCp0JAohcg6XJXtZy9Tp9NRvmexS2qew3rdS\nmVRVkW80k8rS6qGNu7hEHloVOmKUkGeuDNImNrM75nAx2ctappYmn0DZpbRPYb0v+j3gOdunGM2k\nsrJ68GRf1eIsDa0KAzFKSO55QRKSPA9zuJjsZS1TpSXcalRpn8J6378ok6oWrjCaSWVl9TBjmWfs\nxzS0KgzEKCETGvo/RwfK4WKyl7VMFe7hRqVK+xTW+6qTqmhUq9dMKiurh6zjVQMP0dCqMJAohFAO\nF5O9rGWql5tbZintU3x8Ak6qwoSwW5/6Vf0e2KsZWT281bU942XijvA+WaIQIlvo87yi7xxWy3rf\nyNrv6xFrhMzJwtTUTXWDvE72CuB5RV83M7Xe1/xa3TH77Pf1iDVCWtXAAUbn1wzyOrLQX9PkDaZ5\nmbpFWavva8NHUbff1yPWCOnXCv8J6CurRwM1uXhVPlroHx7KdO4ypm5R1uqzjPgo6vb7esQzIamX\nqOn3tuvQQn9dMdO5qG5R1uqzjPgo6vb7esQzIbmabgWy0J+9helcddSqV18bexR1+309YpSQkXUC\nGNndsqV6HXoF5LiLRzOdS6NWJVmrzzLiI0ApjJqYRC9dGyn7fT1ilJDsgFpvG5i8SIG6pp51kMlZ\nbc0uTGrTLnU3PHdjx1argBbU+ap6VDqqYpQQ938DsGS67WaTF/HO5GoX0IEmIMhzb3FBHTfR6aiK\nUUICgvYeoozArUxP3Z2veGXRYARNza48AYHm3uKCOl9Vj05HVTwT0lzz7VX7ERqBW5Z8YvsUxSuL\nD0bQ1ezyCQjy3FtcUKdV9ah0VFUbIT8uCrjpX4vONwT+ruaOUG/wKf+lEbilJbBwheKVxQcj6Gp2\n+QQEee4tLqjTqnpUOqqqjZD0gG7LvqgX2jtURuAW7mGZXNkriw9G0NbsAp+AIM+9xQV1vqoenY6q\naiNkTP2AjWO0yGwR2jtURuCyxim34pXFpuNSKlcu1x2QnJa8DOS5t7igzlfVwzZpCAjVRkhho5B+\nUjD3EB00UteglFSLaOZ2fZEwhGilruW81ajmdn0Rw4Qc/sayEVqP7/NBaWZb02SRsl/NuV1fxDAh\nFwSjBc59QG1mOzxU3a/m3K4vYpiQP7YLwiZ2zN/VZrZ1xep+Ned2fRHDhFxXEMR14/epzWyzt6j7\n1Zzb9UWiEJLjLup63RW8ma34tsJGg67en5PjZjW7PLfb5O+U270j9/KM6OV2fZEohDAoK+lcZemk\nk/JcL9/UY1QRa4S00tyo2wRxHTNjYivprKV/XJMFeeqse7WGlD23It+niDS6iDVCXslTphg16mF6\n1k+LvXNgS5JPYJ0oa+n/79Cyyeqse7WGlD23bYp3EWmUEWuEaGDxldXOV/fO5ivprKWfTbVXZ92r\nNaTsuYUrvItIo4z4JOSeRt5/IYP3kKrq1Ben2quz7tUaUvbc+H3eRaRRRnwSMuYy7yMsp8h0E2vp\nLx6d40bphAvnag0pHxPmXUQaZcQwIY0u8h4RqqCx6fBWn5Z+TQ2pEGPCYpiQ7A6mw7qbm62afJZV\nrqvvzTp481Oa58J9SzYghgmxgO9XlhZm9b1CwH5C5rkGGHXue6NVsOXSwWBofe9c1sQ58OoSLOTV\n1vdSjYNIsJ+Qv9a6yqhz3xsNgi2XDgadfLO9Yyp6bsBCXm19L69xEAr2ExIgIvqVVfmV93rIt0Pm\n8NJdXX0vr3EQCvFJiC88WXm8kFdX38trHIRCohACSeQaoKnvvfUY1jh4r7Br6ko39Gj9pXSEl5XK\nbgHlt7VvsxvILQCLTPE0PhgsbCQMIQagpXXvFXa1rvRU7pm1rC2dl5XKbgF37INVk4DcAliRqXza\n7ntteFOJQgjmdDGNi4vlbCH9Pb60vqbJKwYWAVg4uv55KHtCOsLLSsktAFLeqzrzK3C3ACwypdPY\nYLDwIRgh+/IshtSHg6F3eha9hmlcXCxnOutDvrR+eKixRcDI/bBtrGdEKQCVlcpO/Tv7Xz0N3Rre\nyzyERaZ02pPhmZwQBCMk7Jm3ZqjzD1i5ANO4uFjOdBYtra8rNrIIwMJR91256bsAqKyU3AI+/xVO\ntttHbgFYZMpP44PBwoZghJS4fgj/3RiBzwFjaVyUV0xn0dL67C1GFgFYODrzx+O9pb8FKiuVHWH/\nCb90qyC3ACwy5afhYLDwIRAh/zt27NicSBGCOV1M42Jil+ms4tFYW1o82sAigBeOru+UxwreqayU\n3AJ2pbTr/G/ZLQCLTPE0PhgsfIhDiOdc/G6JTC+D0HPAdBCHELcrd/78u13f2/KGvGHgC6Ct3sV9\nXf0u940t+koGlagAACAASURBVJe0y+MOsi6LOEQi5KkI3kP8Q1+/y31je0j3fsC4Q7YuizgEI+RZ\nVy/TRY6wMPZhamljjWzUzIbVuzipgopQ5EY3kH1j27INxh2ydVnEIRghJedfHUiqOHg0e5Ra2lgj\nGzWzsYCET6rgRShKo5vsG/tf/KOhuINbl0UcghESMbDAA8MO1sjGm9kwIOGTKngRitLoJvvG/nOK\ntOFxh2xdFnEkCiFKSxtrZOPNbBiQ8EkVvAhFaXSTfWNn/l3a8LijKFpp+kQhZH3DW3lLG2tkG5aG\nzWwsIMHCkxz3sssaNUzLufm6R2FmRvuO7dsUtExpcNFnA/9Yt15/jDsU67KIo9oI6XFOIy+4gq3V\nDQp+Oqi8VRU3fR1QEcm3ZIhqI2RRJ69F9lGuSLYpVeSr1ljUQYVzwlBlga+qQtPXYy0GDKuK4Jsy\nQrUR4oMQv7IebmJanaVFx0dUayzqoHprFtz3NqosA1WFpq9fLfPcsy7sDxYcYp6QrrUDsolNeVO1\nxqI0b/4pGPs1qiwDVYWmr9KjP28J+4MFh5gnZPiVAZ1GZgHaDiqcE4Yqy0BVoelr5pHKvr+F8qbC\nQKIQQmYBvIOKp3lxThiqLNmNf2DrtLQDPJuLpq9v3ZgT8nCvUJEohCj4uPecV4BlEnFOWIFXgYNx\nTrE0KXloKO8tFCQcIXwWlScX3AMHjjjt3b5mnFOMpvoViZBrgvRzRzQxs43TQ550y0bfnmCZRCx4\n8GlfM8wpRlX9ikMIXHVzSITU177Gfed7h5uEq5L4pFs2+hYzifjIu33NOKcYVfUrECGhQf+V1a+2\nSW139qt80i0bfYuZRHzk3b5mmlOMnvqNM0LGX2xymjzplo2+5ZlE9si7fc04pxhV9ZsohMiTbtno\nW2xnez3t8syP/1DFDhbclrlZO/5WUr9t686T1O+6ehexnOITderP1MgtnDV1uG+ntTQT1F4kCiGG\noOmfBk4BOOyTtqi0NHIL0465+3+cIZ9gK2KekG661MkfzZoLuT0ASS1MMVbmbWmK0z8r8ivzvMQW\nH/bJt7xmVCO3WNpxP+bi6QRbEfOEjNFNo/pjbZPT0B6gLPmEWk1aNvm3fpha3D6lbLKXVwAO+6Qt\n0OqtLLcw7fhq796jqtQTbETME6KH9isrTVtKWmsVrFxQWgJqNWnpcpr+uXBF6XK92OLDPvmWakZV\nuYVpx4cfrSraKJ9gK+KYkO6NNT7wN25w530vSy1q1qHpn+P3Fe7Viy0+7JNvqWZUlVuYdhy39WxB\nuXyCrYhjQgZpDUyxlFSWWphizHHT9E8fr4AjOOzzEB/5yWtGacwUSzZi2nHnLR2XHKET7EWcETJI\nExg2b6IeNyolzTr4cR8sViw45uVB6ptjRL1rU4+UNeKMkCE11WxJLU2Wy3jElKns9Vlhp/SiLT1S\n1ogzQrTQfWUFJ3u9c4yUXrSnR8oacUzIgOaatuiSYGSvT46R0ov29EhZI44JuUEre89ZGYTs9ckx\nAqYXbeqRskYcE7Jp9jIVqf8KQvb65Bh5etGmHilrxDEhOgQjew1W2JnetatHyhoJQojTQSUYNhSv\n9gk2ZM+GDW1TdvANKG4Nnpndrn0byK0B44+oBCGQMIQYT0RAzwZ329+3PoIbAMWtYfqznhPNPSD3\nTWH8EYUgBIQm5PgbIQ1FMsaLTRbwlimMOKhNCj0bft0KLyzGDYDs1vALs7Le5Qbqm8L4IxpBCAhN\nyCBbnQP6UssURhy8TUr2bOjZsoo2slvDezPkN4GZd4w/ohGEgNCETHd9EtJUJEPMK6aWKYw4eJsU\n92z4BTwjvsANKG4Nxa8CfP0GZd55/BGVIASEJmSWy8aZBMyrVx2By9ukuGdDp1Nnex/DDShuDf/I\n8Zzp/hkl3Xn8EZUgBBKHEObVq47AZW1UsmdDSafbVvENKG4N7hFtM5bKbg0Yf0QnCIHEIYRB4wxA\nTgHaSbgofMkewFD9YqEJWZpFDlEkJNhRkiNtJsR7spfyWKN+SeYaql9WaCJ7lUUOUSTE1yrUH4Kp\nTvvsIcuRR2/iZC/WybZ8Dry6iDsFaCfhovCV7QGM1C8WmsiWZpFDFAkZe2Fwwmh8UH8ht1pze840\nGP86drJtK6ro+QV3CtBPwkXhy+0BjNQvFpqQV1kEEUVCJjQM7kWDu4fc19hylHcOq+fFTrajQ+bI\nTgHaSbgofOXaEiP1i4UmZGkWQcQNIdbu4oD1vNjJxgxjySlAOwkXha9sD2CkfrHQhCzNIohEIQTr\nebGTjRnG4riKiw5qJ+Ey4avYAzD1e16LDkm3t722aVqbuy6vfWO3ipRWHZeM/XOnvP3mCUhe8Bse\nBCbkGdf7QWiy3EtCe+um4G6x9+wGaFU1+0V44G8Aq2qyiQrmCUhe8BseBCYkyFzWeZYvJhc5fKvO\nAePJRvfDWePGYMbxlJFbbHoVlN8OhzqtHQFwqF8G+7oyTUDut8N8Q2BCjv81mKgF/0LeuMPMJnbw\nYF7koM4Bo2Rj6Yvw4DzMOBq5xZ65JKNjnf0At6efAnf+d1nsfZkmIHnBb5gQmJDggPeQ1jVMOtoa\n1aPaXnUOGCUbe1fAI59jxtHILXb3SIDJS6VftxS/P7t43yj2k0wTkKjDwv0c8UVIwXVmT6NNbOEe\nUOeAUbIx5/jxa37DjKORW+ybcwAmfQjQ9wBATnZKC2buYJqARB0W7ueIH0Iu2b9/f59rzZ6Wixw0\nc8B4svEfrQra84yjkVvspM0SGdLvuh1q3Ylsldc8AckKfsP+HHFDSHu8s5v16+iKHLwqHn4WYRaY\ngrgh5KNJ8+fPb9fM5Fle20tlDl6Fvpue0z6yKsE2mFrB08NstAUOrsCi7HuXwLedQ+0SjRtCEOb3\nkGBexLQE23dqBQ9C2GgLHFzBi7J3djvd7btQf3p8EdLlAjOb2OEPZWWf5pEHbjDxS0EHmpjRvlUJ\ntsHUCh6EsNEWOLiCirLT71gf8keIL0IGNDYzMLvuMSheyiMP3GDil4IOjE1o36oE23BqBQYhI/cD\nDq6gouwH+ob+EeKLEHMwm1iKPHCDiV8KOjA2oX2rEmyDqRWUHs5088EVgEXZ64e2+l/IbzRRCGE2\nsRR54AYTvxR0NP78na4naJ9KsLvU7bOvcO9ddd7GyKT7CBaZdJ2IKcT6WYcGsKkVr3S5KbVi4LUp\nO2BSSv3U3Ti4Aouyyzv/+uTckN9oohCCc8B45IEbTPxS0LGm9aVr5QCESrAH9W/tznHfdkl7FplM\nqfOUFJmsOvd9YCnEBoMueYxNrbgg9ZZ+p2u2vvZKyPr8vsVP4+AKVpR9Ov1LKE8KOUefKISwvCKz\nJWXDp6rYTRs9SamVymcMlc44AO/ZlFlk5Yv4ENA/ANd9PW3Ac79tRUKJQkhZ8gm0JWXDp15gN230\nJKVWKu8xVN7GAdI9mzKLWL44kk15Qf8AXPc9eFHGjbe5LX52UIh9QiovaxGATWyLOYC2pOzmjjdt\n9CSlVirvMVR64wC8Z/PMIi9fzGS/ffQPwHXf1U+CJ2e3XR8n9gn5xZUSwLiKFtsAbUnZzZ0nGJkn\nKbVSeY+h0hsH4DgqnlnE8kX2UIromX8Arvs+8w6cTjlu18eJfUKOu14M4Kxebj7ok93cMcGInqTU\nSuU1hkpvHMDHUWFmkZcv8ofoH4AFjwOS05KX2fZxYoyQdTN8Kq5edM0O4YXQkzT892M/YoyQi4wW\nbwN15NEOmXo/b+CI0+AzB9fHwgFgbHTaEGTEGCG3dvaey72/zPV8YNfqSklN5uD6WDhQYUP0EGOE\ntL/N51Bg9xB1yJR2QIJPAOIzJIEKG6KH2CfkZ1dds4V0DS6aLg+Z0g5I8A5AfCwc5MKG6CH2CYEe\n9wYwnPiGY/KQKe2ABO8AxMfCQS5siB7sIuRp/4U6+fX9v5A/GBESENQhU9oBCd4BiI+Fg1zYED3Y\nQ8i6gErZzg3/B4VMyLJ6d2aCz4AEbQDyx0+zDj7Q8sq6lyZ1aZ1268BdKS0vvE7iaOIEg4Vb4Ku2\nvIWHNfrYUUTKYQ8hsGS94H8hgYIXkLI1WblE1GDhVjqKq7asdJR3uNtRRMphEyGBIFL3kIBQmadb\ntJVnevJBnzznqzbvYEJXbt4xWLgFngHG0lFM+NpSRMoRbULcRyzbOPwhtUcgP+mU74WfPqFdtFVm\nevJBnzznq23ewYQurs4aL9ziGVQ62rNllS1FpBzRJqRNcBXUPjCr89HhfN/rzpuvXbRVZnpi4pdy\nvmrzDiZ0aXXWaOGWHZDOwNJRTPjaUkTKEW1Csi+dFQ6u7BjIT6qX7nNh5jfaRVtlpicmfinnqzbv\nYEKXmncWzsN/qH6Uay/gKV8sHcWEry1FpBzRJmRAeEa3gd1DGvo2yua6tYu2KLdYszo3zOI5X7V5\nhyV05eadBz7B66l+lNtnAU/5YukoJnxtKSLliANCnuviPSyk1mDvc6pyMI1ItYhD5fRUwTHvIVS8\nZx1FLTk2GdUrAj8D26ZwWAKdawPigJDWtbxndZ/b3fuc4Yu1aUTl6V4+9Yka2UuOTRayF9um2LAE\n+VwbEAeEDGzpfcTgK2tNkwXpqV/m8d6ploNv2bNxKsxeke8zhEqVvbJjk7ns5TUPSrFDWB9MQVwS\n0uBOn6b3T4eWJZ8om4y9Uzs7np335t7RJ3O3TfEeQqWTvdi4biF7seZBLXYI64MpiEtCavjK3nOm\nlZZA6XLsnfrLInh0Z0XvP21auMJ7CJVG9pJjk4XsxZoHtdghrA+mIC4JeXvhMm+M2lK4Bwr3Yu8U\nylzo1BfG7/MeQqWRveTYZCF7seZBLXYI64MpiEtCDFA8updbkra8d4rJXEjdamAQq8pe2bHJQvZi\nzYNa7GAPYoyQqy72nWR4WeOAL9f0Th0Wqm9KRYwRkn2DASGBWwaw3ikeeXyWVe5b4kDGWdwS634p\nmDlUq4yXPchNUdgmhRv7Ig89YowQIwT0laVCE3l4lzhgECJbYnXLBZjY9Hde9iA3RWGbFNvYGHno\nEUeEHHzF51auYt3AzAk0/pbVW581LHHAIIQssdxd7oAjQzKo7IGaorBNCjc2Rh56xBEhHa2yxOfc\nCyU0/pbVW5cZlTjwIIQssfbek+9+fOVdvOxBborCNinc2Bh56BFHhIxpaOGGtnNMx1U0/pbVWxuW\nOPAghCyxlr98357Bq57nZQ9yUxS2SeHGxshDjzgiZJzVzZ3l2mn8Lau3Nixx4EEIWWJN/mh67idP\nfcDLHuSmKGyTwo2NkYcecUBIm5o8pVj/fIuTWK6dj6pws3prnxIHNQhZ3ymPZUVyTi7qCvlHcXKF\n0hSFbVK4sTHy0CMOCHkinSfdk+oFfg2blWBqEMDV78st09oM87Q6ABVp+0nj+lHCtpSexAEhMiy/\nsrzABK+pQYCPdZmsca2VsD2lJ3FEyLB6Fmu/b0yVdC4TvCeY0l3TZIGFQYCPdRlpXD9K2J7Skzgi\nJNtK9taaJelcJnjfYEr38FALgwBf6zLSuH6UsD2lJ3FECFjVD33FdC4TvKh01xVbGAT4WpeRxvWj\nhO0pPYknQqyAOpcJXlS6bFaC6Yw2X+sy0rh+lLA9pSeJQgjqXCZ4UekWj276s2ZGW37dxhduzD2U\n3XStpH4b3NuyfuttC9P+NbPbBdMA2m1MTR7b6aYrkodCt6w/3mqhhO0pPUkUQryh01ZcVpFIMnSB\nRdsllFKRRhQJGc+837vfEL0fqIVHsYelFKPPLCpZJBm5wKLtEkqpiL/RKBKSjnKnboReffgA7+os\nHQY9JtvDUorRZxaVLJKMXGDRdomSihFGFAn5dsGGDRsyTG0qw0SN2pYtbXV3yfawlGL0mUUliyRD\nF1hmuyQbMEUWUSQEEbF7SL2HLZ9W7WEpxegzi0oWSUYusGi7hFIq4qgWQr7f69PcHDbqPmT5c1V7\nWEox+syikkWSkQsszsJFKRVxVAshBnVT4SPD7882G4jLMo0swYjpROBDD+RaXl7lG+nBUxpUCyG1\n7rBYbA0R54/1+7O97GEVyEvrPJ3I5a4sflkNb+QHT2lQLYTUftT+F/ZzD0GByxbSmZMZ8y7Tmggs\n4O1smE4kuSvX8rIa3sgPntIgbgipeaVvgZAGnQuhBJvaypJPoHeZ1kSAt7PxdCLJXbmWl9XwRn7w\nlAZxQ8iVzS0NzDoM77gKF9JLSwC9y7QmArydjacTSe6S+MUa3sgPntIgbgjxAyZwMcFYuAfQu0xr\nIsDb2Xg6keQuiV+s4Y384CkNEoUQJnAxwdjLzb3LtCYCvJ2NT0IguUviF2t4aZU9OkgUQjQ4vjpr\nXjDeZb3Y0rvXunuubaaX3khAQtwD0bssUHiyDY7l+h6zCYlCiOLXW8k621AA4yEmfvVL6sqZdKQ8\nNTf5ME0QQ62M1lsRe6NxR8g/LrzKuweU4drHya8XO9tQALMV9teY+NUvqSvOvnTkncGexXNpghhq\nZZzaGrH3H3eEvOi6zcgm9sr/R3692Nn2rbzCjuJXv6SuOPvSkalb4N1XyccfT8eprRF7/3FHyN9c\nhhlAxa8XO9vUFXYmfvVL6sqZdKTPprN3/EA+/ng6Wm9F7P2LTsin1lO5fTHSmBDFr7fpz1lZblph\nb7RJEr9JB/UFpezMYbUbXrUx96XmTdq0+K5+RtbitG7u7sMlhbyjRu6i9KKkztfoxq5ixSL6B9gA\n0QlJDT7tu9XyBTWL6eZayXfsKhkFGC6442I88w8I+FNZQHRC7r00SAOn143/Qph2Wo9eASSY+IR7\nJrmUekarsas0EsFowZ0vxjP/gPB/PeITct8fgvwBJvcQpp1Wo1cACSY+4Z5JLqWe0bJ2EY0CDBfc\ncTGe/APCR9wRMt1V33BNXaKJewWQYOIT7pnkUuoZLWoXaSSC4YI7Lsajf0CQb9UQcUdIWaeHjWxi\nWx44MYF7BZBg4hPumeRS6hktahdpJILhgjsuxqN/QJBv1RBxR4gJmHaSvQIwpYgT79FMQKlntBq7\nykciGC6442I8+gfYgUQhhEFUrwAdRCLkw3W+3rI5gds0+AN6BSjopcvX6pO5vIUKYwxsjTJzMfOa\nu0rnhodoE5J7uRS81TQsGagwjCrOi8z7MMrhKsAwhMcYVM9r6GLmPXfVltrfaBOShL/lTKOnfnWN\nMvgLUfvUjgYZkejws1zkoM3hYpqXhSE0lBVPAApDMMaQ63kNXcy85q7aU/sbbUJO79m/f//5RUZP\n/eoyWDZS7yHTwqjZkv7Q7qUiB20OF9O8LAyhoax4gupiJsUYVM9r5mKmm7tqT+1vtAlBGN9D/BBS\n7JoVZFpLi7n3UZGDNoeL0YdmKCueIIchGGNQPa+xi5nX3FV7an9jh5DXXQfC+JFKkYM2h4vRh2Yo\nK54ghyEYY1A9r7GLGXf2ZRGK2kYVNhKFEKXIgYKNvj2ylvIyX81QVjyBwhAeY1A9r6GLGTn7sghF\nbaMKG4lCiJfM1YIX9hqe71PaILsFRA7CE6Jke18KixBLmRvg+Z5cxS0gchCdkLYapeT1heBea+H+\n442NTOZ65PHEZ2TLLMy/s8LeQEsbuFtABCE6IesKZZ10l+sH/VOTg5G95/SWZK4ynniqbJmF+XdW\n2BtoaQN3C4ggRCdEhc89ZKHrr/7H+sgYxmSuMp5YsczC/Dsr7A2wtIHcAiKIGCbkL64gCjy5zJXH\nEyuWWZh/Z4W9gZU2yG4BEUSiEIIyVxlPrFhmYf6dFfYqRRBaWUxqGC9j6XrFLSCCiCVC/qNPT70a\nDCG+MNK6CjT5Xszm+s/34p4N3rGxQ8iDvnfqssi8PdDle6nexE++F/fs8I6NHUL2jZynT08Nd+3T\nPP3CEssWxMlz4NUl2MUmCV0+lU1RuqdlC1mvshM1mwt+8724Z4d3bOwQ4gPdPWSPP9k7pqLnBuxi\nk4Qun8qmKF3m52A8tULO5vrP99Je+N6x8UJIuetZyzb2/xsyZxXKV0no0lQ2Rekely1kvcpO1Gyu\n/3wv7tnhHRs/hJRanuzJyvOgfJWELk1lU5SuYiHrVXaiZnP953txzw7v2Bgm5E+uBmrdVQPXVOuz\nkz6Q5GuzuVho0tPdk6pNUOnObNRhu0HZye8X1q5XbyMsvJQppwfyN0sv8ty5zVv9rXb9il31buj8\nxNPl2TXr3NANoLRx8tLLr6l7wY6SmxpcHa7UimFCtrcbqbr9DHa9EMAlvLA30EQjl1pcOa2qySa7\nmUstPM0OqRXDhOjg9yuLmQVU5OOqebmSaGSSi6UXSV4ZSS2unA71y2D3cHOphafZIbXihZDvXP0t\nR4SOmAH3vb19Cq6aY+Jwoyy5WHqR5JWx1JKUkzv/uyzwI7VQYIUvtaqFkPNSjDzGRris3TEsUe5H\n9p77NYz9euEKzChOVRKNLLfI0oskr4ykFiqnZxfvGwWWUgtPs0NqVQshjS41agO82mXi1mMH0Cxg\n/D7MKKqJRpZbZOlFeUSugdRC5ZSTndLiY0uphafZIbWqhRBjhPWV5Q9oFpDDW6fURKMkuTC9SPLK\nqMKXXPcnMuFksbTOTlP3wnijCULI8dWwLIJ024gEIUQ1C9AmeWkomDKg2HcsAs/mymMRAOaU0vRi\ngNd6JN19FvikBM+zWaljz5IDGj1T+jIcqRf8tIQEIcQY5FymtB36jkXAbK4yFgFg+E5+DsAzz3g8\no1cDn5RQ9LLHfddSupKeKfwIivrPD/o9xRgh21NHWbrBmmLSNNnFQetctoiFJRX5eAgMWgt5NlcZ\niwDQ7gw/B368XRK8pyv4pIT9kkiGd6bwK+Vnup48krn28aB/CzFGyBOuK4wEmn80e1B2cdA6l2FY\nsn0KHjJuLcTAgsYiAFR2kIOTJS/B+xlJX/NJCSULpTMXLOBX0jOeVCha+2VB0L+FGCNkkSvEqijV\nxUHrXIZhycIVeMiotZCyuTQWYTXsuEeeXjyLOZh2/p1PSnhKuhhy9/Ir6Znv+x5pmt0tPeg3KiAh\n/2dRXjUtVEJUFwetcxmGJeP34SGj1kKezaWxCBtawSOLaXoxrB7k9ixIAj4pYeXdbs+c8eSARs+s\nerpoDXiCdykWj5DN1iF3iD1KmNbFmEPrXIZhCR0yaC2k6l15QNjolCFVNL0YPONvvn3hCMBJCeC+\n/9Zuf3LTlfTMn97oLN1LWp3x8758IB4h/3YVma/E3h/qX4gvTGeCcSVb9C9pl9cyoGuGXdYZ/iAe\nIdtc75mfE/I9xBemM8G4ku3xi7ThCXZ0zbDLOsMfYoyQWa40SzdYU+Q9prUsk/TvQtOZYHzmV1u2\n4Ql2dM2wyzrDH2KMkOWNu4VGSI8XtJZlkv41nQnGlex/0ZgGE+zommGbdYY/xBghIUNvWSbpX9OZ\nYFzJ/nOKtOEJdnTNsM06wx8ShRC9ZZmkf01ngnElO/Pv0oYn2NE1wzbrDH9IFEKKbr6CDwVjEnd1\ndubBseVr6nfnM8HkrHt60h+u3wj5teu1GA4DL26VdmDXpS07/7XHguZ1WqT16XpFK7Qs045LGLuZ\nakltRKIQIosqBE8q+viXaepHqagh6/gnOSe47CLLMtCMS5DOkbOP9iFRCKFBq5hGZOMQsCHKvH6U\nFzV4bt3U5zcuu2TLMs24BHYOZR9tROwS8nmdYPKMN+WiqMI0IhuHgA1RpvWjVNRwsEnunVx2PaJa\nlsnjEugczD7aiNglZLUrw8gP1gR5j6OoQo3FxiFgQ5Rp/SgVNawecDbtMJddsmWZOi4Bz5FHhtmH\n2CVkiyuYwbMkqlBjsXEI2BBlWj9KRQ3PrIT502GhxIViWaaOS8BziuwxydIgUQihQauYRmTjENC/\nzKh+lPKJWNQwoBxOtqh64FNplyzLdOMSJu6Qz7YRAhFySk7ovhzQ6cERosPQn+S9gmNeC+ooYymR\nKLdL6Ux6b09ObvYjKt7IQCBCYNAY5o441LVUOZKRap4MucX1bqg/qLuy18trQZ3LWEokGpr0Auy+\nV1bFkYBIhHBov7JqXmY+xOgGDXH+oa3erWw5+BbujVWRT+W+oBO9lEg0MukFqMw7KZf6RgJiE3Lh\nJPPz5K+sdD9FpBwXdlOrd3d2PDvvTVpOp3JfneiVE4lGJr0AT24GUrwRQewT0rWxkS2sNwaVqtW7\nf1kEj+6k5XQq99WJXjmRaDgVtzwfZFUcEcQ+If1aBfKq2upd1Rtr/D4q99WJXjmRaGTSC+M+VFRx\nRBC7hKx3DcL/799wRSCvqq3eVb2xWGrR1yRL9uA1Muk92gVvHhPDbUc3RewS8sG53FO8VoOAXlb1\ny1IKSLMOmq6rVx9ilxAZgX1ladrYNAWkpuvq1YfYJyT/+kCme29T2tjWNHmFaV+2mL7FdF29+hD7\nhFwfkOytcYfcxnZ4KGpftphuuq5ejYh9QtY/EYhNbP/NchvbumLUvmwx3XRdvRoR+4QEBrWNbfYW\n1L5sMd10Xb0akSiE9L24E7axTTxWPJprX0nxji1v07JtRzbAwt2u3eWZauKXQZNsLE1KHiodmfQv\n+KLHUVbY+Jicb/RaYcfxYOFBPEI+cF2r5Ktq2P593se72NZUaGmSjQMq8IjfFXY+Hiw8iEfID1fn\nKOt85/Wz61XJYHRNk1KmtNBGANt3zIWWmmw81mLAsKpAVtj327A2Ih4hWph8ZX20NnD3S8KyydSj\ng0rrrVlw39t8CKuZ0NIkG79a5rlnXSAr7DgeLMyPHIuEnDknIKWrwznzqEcHlVb+KRj7NbbvmAot\nnUMT/HlLICvsOB4szI8ci4T87hoR9F/I4L3Uo8MtRpmNALbvmAotTbIx80hl398CWWHH8WBhfuTY\nJCT4OycZjBaPxlV1tBHA9h1aa9cWMPqssL91Yw6bHup3hZ2PBwsTiUKIrj9dtREwSipq9C7JWPPF\ndRtsSL2QMIRoodoIGEGjd0nGmi6u22FD6oVEIaQyD9fO92b2a3nIw4RupWpNajQIVzPl1mJx3Q4b\nUi/EASrzSQAAC01JREFUMCFvvBjEyKMXnsC18+UFnkWvodDVWJMaDcJVp9yCxeK6HTakXohhQs4N\nTvbi2vnkT2DlAhS6GmtSU70ry1jTxXU7bEi9EMOE1L43kIUQwr187bzXZnfe9yh0NdakpnpXlrGm\ni+t22JB6IYYJqWM4fM8EuXztPKt71mLes44q2HQQrmbKLVgsrofpjWUEsQk5r7bRKO5GrrvYk8EQ\nUpWD/xwZrihdLCddfYwPoOJKl2tY7r1vYMg0p5Qb8oN8TkQgNiGDhhnVWD3kwgmnwRAynGuhDcXK\nESwnlXO/GkNY8t43MGQavpO7xIJ8TkQgNiHGCOUr6/Wp+13Hp6xHfyZM9bJyUtZKxQZQaQ1hufe+\ngSGTtMddYkH2548I4oKQH/Ie8Ve5OPTJojY/9dzBEryY6sVyUtZKxQZQaQxhyVnfwJCpsgO5xAKd\nExnEBSFP+Ze9548bUDD3TUzwYqoXy0lZKxUbQKUxhCXvfR9DJrZHLrFA50QGcUHIEtc3/q754oal\nBd0rMcGLqV7K/fIBVBpDWPLe9zZkOsD2yCVW9uePDBKFkO8uqyyYwRO83DGW5375ACqtISzL3xoY\nMrE9comVz4kMEoUQFYI7xsY4Id/NYi5aDwRGSK4762DBsczenQsw1bvaK8b28cnCQhKqOIkaYpyQ\nTLpnfxrAZejcwKpLKPrwLkDx8cnCQhKqOIkaYpOQm1hJSg0p7H6oNktUvSj/hfxfcl+TPvUx47re\nfX9lHqsuIXdYqvG18MlihSRUcRI9xCIhZ+v/gZkznNMboKgOO6DcQ2a7/mji5NB0Oox/vWwyqy4h\nd1iq8TX3ycJCEqo4iR5ikRAC+8ryIuRt11cmJ6NzQ+lyVl1C7rBU42vukyU3rrGKk+ghLghZ7FrJ\nXWRnmBKCzg2Fe1l1CbnDUo2vuU8WFpJQxUn0EBeEaAZ5f2RyMjo35GAZL7nDXl9QdNEQlnbPHZs1\nrBzI8XVg67S0A4W1L6j7KmxKbfw8PH5BfZaBtyzrRSFmV7lDXBBy6iWqvpps+hdihADKeqlu17qs\nlwkx28od4oIQBeb3EKxuwDqHyvTUmXPg1UUBlfXyul3rsl4UYraVO8QXIS+5kkx8H9oM8Cx6DdVV\nWfKJbUUVPb8IqKwX63b9lPWSELOp3CG+CFnduKuJM8q1a2HlAlRXpSVwdMicVYGV9WLdrp+yXmBC\nzLZyh/gixBxY3YDqqnAPzmENrKwX63b9lPWiELOt3CFRCMHqBlRXrGE96YMAy3qpbteyrBdLf20r\nd0gQQo74m8inKeglAWtd5hAxw6wEIURT3WAMVenKAtayzCFyhlkJQgjIo2+/ZUr39MDMCaZKlwSs\ndZlD5AyzYpiQGg2vuaZhDXpwTWNrm9jmg7kHPyrddx+CElOlSwLWsswhgoZZMUxIVp/Ro1Nq0oOG\nLayntN2UxT34UekeH9NxlanSJQFrWeYQQcOsGCaEQfnKajbS+kTZgx+VbpH7236mSpcErGWZQwQN\nsxKFEHn0LSrdmRldtpsqXRKw1mUOkTPMShRC1PmqWvjpaCNLWEzpKioYn9l3W+udfGP3mns8ELKl\nSdes87v7OTNgTyztjAS0hOUpXVkF82f6/PzBg3xj95p7PBAyx9U+6/yu1idW5KPOPda3b8NvA+1o\n45awPKWrqGB85ptxsP0R3Ni+5h4PhCx37aavrAtMS0nrTkCdC/Dn1wPuaKOuNUzpKioYnylZBivn\n4sb2Nff4IuTCjrNMMHAp6lwonQEBd7Rh1xrwlC6pYOpnm/I5FO7CDdi95h5fhFwyzuzE8ftQ566a\n6IGAO9qwa41SulwFf0D9bI+uO96Tb2xfc48XQhphYFh7kNmJPKN7/MK8vksD7mgjS1hM6ZIKpn62\n/9yU/TXfkMuDfYgXQi6vy3rdzjW/tWsKR7GbTVW8Xj1tKGbh9/RbG6WNwOpSvgueVgegIm3/kXpr\npahxDcC92yLyieKFEA7zryxt4aiBQtb0tKGYBXYvl6tLcRdg9ovwwN+gqP98SUNLN/x0cy+IcBBf\nhFx8t1lb9N94/SjL82I3W2Ue93NgPW6g7WlDMcteq2SBXF2Ku1JY2GntCDiSufZxgBQPuNtF5hPF\nFyHnmcrec/thhQPL82I3W9lkrn5Zj5uupw3FLHutcZ9TdSnuMtyefgqK1n5ZAD83yMhoOzgynyjG\nCRlXU5K0Ba5d9HDBa2Zzv6e9iBUOLM+L3Wyly1H9Yo+brqdNFrPQ5TeqLsVdhidXwpGm2d3SYeMk\ngL/PjswninFCevH//2/yeyLVj7I8LyleVL/Y46braUMxK8HTTq4uxV2GvgegaA14boTnl0j39XDt\nR00Q44RU/Xjs2LFFyleWOah+dGZGl5d6ZB5kihfVL/a46XraUMxK+G6gXF2KuwztPEc6S0K41Zm7\nv5L+r/DfyHyiGCcEsTwAQlT4GMWKhUQhhMkrZe4qWgfIq+xlotjwc8QLIaV+TDCfleSVMncVrQM2\n0qR7YWz4OeKEkLf8OgcM6ajOXUXrAHnSvTA2/BxxQsgW12vbLTGCLaMrc1eZdYC8yi6MDT9H3BDi\nx8CYLaMrc1fROkBeZdfnFqsdiUKIFsbWAQGNocJJVl5n5tpqfBIfhPzzrWAIsXaJtYJmkpV6LDek\nlzJDPBAyg920F1qfo8heErjcL5aUL++oYoWmeIAtt7PzJeoey8quYklIuqqcTbKiQQt4cVVh17vv\nt/7JQSIeCPllwryHXOvlRy/ONrKJnSTLXhK46BcrK1/eUSXt8ANsuZ1lIQHYkvsLLAlJV70z2LN4\nbhkftIDnlrwE41+39cPEAyGgm+VdwzjbexfJXhK46BcrK1/eUSXt4AFcbmdZSAC25I5JSLpq6hZ4\n99XS5XgGnovt77Z+kvgjpOF9xwzwoCx7SeBiR5WsfHlHFZPA7AAmH1kWEoAtuWMSkq7CSVaFfNAC\nH7PA2t9t/STxR0ijQqMTFNlLAhc7qmTlyzuqpB08gMvt7HwAtuSOSUi6qm+PrKXyoAU8F9vfbUWc\nELJM8+0USKOA346qakOcEPJL14my3WXtggDON+6okgMMbtZAsyp4F5s8sUJ20lIrfNUZbbYgTgjR\nwPgrKyhwswaaVcG72OSJFeSkpVb4ama02YL4I6ReqpFN7PQpzMSBxQ9o50BxhUcXd1Tm8eiDmzXs\npVkVvIuNJlaQk5amwled0WbP248/Qi6sa+ROfulsZuLA4gcsdqC4Qh93lE3m0Qc3a6BZFdTFRhMr\nuJMWe1qu8FVntNmD+CPEGFjigPED7lFcoY87Spfz6IObNdCsCupio4kV3ElLW+GrzmizB4lCCJY4\nYPyAexRX6OOOwr08+uBmDTSrgrrYaGIFd9LSVviqM9rseaOJQgiWOGD8gHsUV+jjjhw3jz64WQPN\nqqAuNppYsSulXed/g7bCV53RZg8ShZCYgUOIYHAIEQwOIYLBIUQwOIQIBocQweAQIhgcQgSDQ4hg\ncAgRDA4hgsEhRDA4hAgGhxDB4BAiGBxCBINDiGBwCBEMDiGCwSFEMDiECAaHEMHgECIYHEIEg0OI\nYHAIEQwOIYLBIUQwOIQIBocQweAQIhgcQgSDQ4hgcAgRDA4hgsEhRDA4hAgGhxDB4BAiGBxCBIND\niGBwCBEMDiGCwSFEMDiECAaHEMHgECIYHEIEg0OIYHAIEQwOIYLBIUQwOIQIBocQweAQIhgcQgSD\nQ4hgcAgRDA4hgsEhRDA4hAgGhxDB4BAiGBxCBINDiGBwCBEMDiGCwSFEMDiECAaHEMHgECIYHEIE\ng0OIYHAIEQwOIYLBIUQwOIQIBocQweAQIhgcQgSDQ4hgcAgRDA4hgsEhRDA4hAgGhxDB8P8B+VhZ\nP9a9/ukAAAAASUVORK5CYII=\n"
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%R -w 400 -h 600\n",
"library(ape)\n",
"\n",
"## make tree ultrametric using penalized likelihood\n",
"Vtree <- read.tree(\"~/Dropbox/RAxML_bestTree.VIB_small_c85d6m4p99\")\n",
"Utree <- ladderize(chronopl(Vtree, 0.5))\n",
"Utree <- drop.tip(Utree, \"clemensiae_DRY6_PWS_2135\")\n",
"\n",
"## multiply bls so tree length=6\n",
"Utree$edge.length <- Utree$edge.length*3\n",
"\n",
"## save the new tree\n",
"write.tree(Utree, \"Tvib.tre\")\n",
"plot(Utree, cex=0.7, edge.width=2)\n",
"#edgelabels(round(Utree$edge.length,3))"
]
},
{
"cell_type": "code",
"execution_count": 160,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(((((((((((((anamensis_C6_PWS_2094:0.2976564833,sempervirens_combined:0.2976564833):0.141895001,setigerum:0.4395514843):0.04243491827,((tashiori_D30_TET_YAH:0.3525423068,luzonicum_D27_9M_2005:0.3525423068):0.08217954946,erosum_D23_MJDJP_4:0.4347218562):0.04726454628):0.07166703939,((parvifolium_D29_KFC_1953:0.3445468166,dilitatum_ELS45:0.3445468166):0.1240260218,wrightii_D31_MJDJP_1:0.4685728384):0.08508060353):0.06099278658,japonicum_D26_WC_273:0.6146462285):0.08345145638,integrifolium_D25_KFC_1946:0.6980976849):0.3154380055,((betulifolium:0.5142650528,muhalla_D28_WC_274:0.5142650528):0.1041896534,foetidum_D24_KFC_1942:0.6184547062):0.3950809842):0.1548603709,(acerfolium_ELS88:0.7716194155,orientale_DRY2_MJD_GEORGIA:0.7716194155):0.3967766459):0.09417177866,((((glaberrimum_D34_PWS_2323:0.4241929153,vernicosum_D21_PWS_2123:0.4241929153):0.1822087602,beccarii_combinedX:0.6064016754):0.2342803007,sambucina_D20_PWS_2100:0.8406819761):0.140414823,(cylindricum_DRY1_WC_268:0.3033948796,coriaceum_combined:0.3033948796):0.6777019196):0.2814710408):0.3155053291,(sargentii_RCW19:0.4218466302,opulus_D6_WC_250:0.4218466302):1.156226539):0.06460952668,(((acutifolium_DRY3_MEX_006:0.3891310523,sulcatum_D9_MEX_003:0.3891310523):0.1111911591,(triphyllum_D13_PWS_1783:0.2564483885,jamesonii_D12_PWS_1636:0.2564483885):0.2438738229):0.2769195884,(dentatum_ELS4:0.5058594461,recognitum_AA_1471_83B:0.5058594461):0.2713823538):0.8654408959):0.06743502437,(((cinnamomifolium_PWS2105X:0.3487440909,propinquum_DRY4_WC_276:0.3487440909):0.06803468874,davidii_D32_WC_269:0.4167787796):0.6313883324,tinus_D33_WC_277:1.048167112):0.6619506082):0.5042243598,(((((((prunifolium_ELS57:0.1699967845,prunifolium_AA_22586A:0.1699967845):0.1640201856,rufidulum_ELS25:0.3340169701):0.193303849,lentago_ELS85:0.5273208191):0.7484887363,(cassinoides_ELS2:1.212204852,punctatum_D19_PWS_2097:1.212204852):0.06360470335):0.2416862958,(((veitchii:0.2378830264,rhytidophyllum:0.2378830264):0.3752483916,lantan_combined:0.6131314181):0.1488822671,((carlesii_D1_BP_001:0.1568053155,bitchuense_combined:0.1568053155):0.507361206,macrocephalum_D2_WC_284:0.6641665215):0.0978471637):0.755482166):0.533597105,((((((((foetens_ERAD10:0.3020842316,grandiflorum_ERAD11_Wendy:0.3020842316):0.1232923404,farreri_RCW21:0.425376572):0.1582832473,suspensum_C5_MJD_111711:0.5836598192):0.2987513558,erubescens_RCW36:0.882411175):0.1115548786,henryi_D22_WC_272:0.9939660536):0.1515051274,(sieboldii_AA_616_6B:0.9700046214,odoratissimum_combined:0.9700046214):0.1754665596):0.2601275914,((plicatum_C1_MJDJP_12:0.6181807158,hanceanum_D4_PWS_2195:0.6181807158):0.1554133605,lutescens_D35_PWS_2077:0.7735940763):0.632004696):0.3575494042,amplificatum_D3_SAN_156003:1.763148177):0.2879447796):0.07196593945,((((furcatum_combined:0.4369797903,sympodiale_D18_KFC_1932:0.4369797903):0.2037853101,lantanoides_D15_Beartown_2:0.6407651003):0.1452765729,nervosum_C4_PWS_2298:0.7860416732):1.220475926,(taiwanianum_TW1_KFC_1952:0.4868234324,urceolatum_MJD_Japan_8:0.4868234324):1.519694167):0.1165412966):0.09128318429);\r\n"
]
}
],
"source": [
"#### load TVib tree into Python and print newick string\n",
"Tvib = ete2.Tree(\"Tvib.tre\")\n",
"! cat Tvib.tre"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simulate sequence data on each tree\n",
"Here I use the _simrrls_ program to simulate RADseq data on each input topology with locus dropout occurring with respect to phylogenetic distances. Find simrrls in my github profile."
]
},
{
"cell_type": "code",
"execution_count": 161,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%%bash\n",
"## balanced tree \n",
"mkdir -p simulations/Tbal_full/\n",
"mkdir -p simulations/Tbal_mut/\n",
"mkdir -p simulations/Tbal_cov/\n",
"\n",
"## imbalanced tree\n",
"mkdir -p simulations/Timb_full/\n",
"mkdir -p simulations/Timb_mut/\n",
"mkdir -p simulations/Timb_cov/\n",
"\n",
"## empirical Viburnum tree\n",
"mkdir -p simulations/Tvib_full/\n",
"mkdir -p simulations/Tvib_mut/\n",
"mkdir -p simulations/Tvib_cov/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### show simrrls options"
]
},
{
"cell_type": "code",
"execution_count": 163,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"usage: simrrls [-h] [--version] -o outname [-mc dropout] [-ms dropout] [-e error] [-f datatype] [-I indels]\n",
" [-l length] [-L nLoci] [-n Ninds] [-N Ne] [-t tree] [-u mu] [-df depthfunc] [-dm depthmean]\n",
" [-ds depthstd] [-c1 cut_1] [-c2 cut_2] [-i1 min_insert] [-i2 max_insert] [-r1 seed_1] [-r2 seed_2]\n",
"\n",
"optional arguments:\n",
" -h, --help show this help message and exit\n",
" --version show program's version number and exit\n",
" -o outname [str] output file name prefix (default 'out')\n",
" -mc dropout [0/1] allelic dropout from mutation to cut sites (default 0)\n",
" -ms dropout [0/1] allelic dropout from new cut sites in seq (default 0)\n",
" -e error [float] sequencing error rate (default 0.0005)\n",
" -f datatype [str] datatype (default rad) (options: rad, gbs, ddrad, pairddrad, pairgbs)\n",
" -I indels [float] rate of indel mutations (default 0) ex: 0.001\n",
" -l length [int] length of simulated sequences (default 100)\n",
" -L nLoci [int] number of loci to simulate (default 100)\n",
" -n Ninds [int] N individuals from each taxon (default 1)\n",
" -N Ne [int] pop size (Ne for all lineages; default 5e5)\n",
" -t tree [str] file name or newick string of ultrametric tree (default 12 taxon balanced tree w/ bls=1)\n",
" -u mu [float] per site mutation rate (default 1e-9)\n",
" -df depthfunc [str] model for sampling copies (default norm, other=exp)\n",
" -dm depthmean [int] mean sampled copies in norm, 1/m for exp (default 10)\n",
" -ds depthstd [int] stdev sampled copies, used with norm model (default 0)\n",
" -c1 cut_1 [str] restriction site 1 (default CTGCAG)\n",
" -c2 cut_2 [str] restriction site 1 (default CCGG)\n",
" -i1 min_insert [int] total frag len = (2*l)+insert (default 100)\n",
" -i2 max_insert [int] total frag len = (2*l)+insert (default 400)\n",
" -r1 seed_1 [int] random seed 1 (default 1234567)\n",
" -r2 seed_2 [int] random seed 2 (default 7654321)\n",
"\n",
" Example usage: \n",
"\n",
" simrrls -o test1 -N 1e6 -n 4 -L 1000 \n",
" simrrls -o test2 -f ddrad -c1 CCTGCAGG -c2 CCGG \n",
" simrrls -o test3 -f pairgbs -c1 CTGCAG -D 1 -i1 100 -i2 400 \n",
" echo \"((a:1,b:1):1,c:2);\" > treefile \n",
" simrrls -o test4 -t treefile -df norm -dm 5 -ds 1\n",
" \n"
]
}
],
"source": [
"%%bash\n",
"simrrls -h"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simulate RAD data on the balanced tree without missing data"
]
},
{
"cell_type": "code",
"execution_count": 98,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\tsimulating rad data\n",
"\tTHETA= 0.002\n",
"\tcreating new barcode map\n"
]
}
],
"source": [
"%%bash\n",
"simrrls -mc 0 -ms 0 -t Tbal.tre \\\n",
" -L 1000 -l 100 \\\n",
" -u 1e-9 -N 5e5 \\\n",
" -f rad -c1 CTGCAG \\\n",
" -o simulations/Tbal_full/Tbal"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### And with missing data from mutation-disruption"
]
},
{
"cell_type": "code",
"execution_count": 138,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\tsimulating rad data\n",
"\tTHETA= 0.002\n",
"\tcreating new barcode map\n"
]
}
],
"source": [
"%%bash\n",
"simrrls -mc 1 -ms 1 -t Tbal.tre \\\n",
" -L 1000 -l 100 \\\n",
" -u 1e-9 -N 5e5 \\\n",
" -f rad -c CTGCAG \\\n",
" -s 300,600 \\\n",
" -o simulations/Tbal_mut/Tbal"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### And with missing data from low sequencing coverage"
]
},
{
"cell_type": "code",
"execution_count": 147,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\tsimulating rad data\n",
"\tTHETA= 0.002\n",
"\tcreating new barcode map\n"
]
}
],
"source": [
"%%bash\n",
"simrrls -D 0 -t Tbal.tre \\\n",
" -L 1000 -l 100 \\\n",
" -u 1e-9 -N 5e5 \\\n",
" -f rad -c CTGCAG \\\n",
" -d 5,5 \\\n",
" -o simulations/Tbal_cov/Tbal"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Assemble data sets in _pyRAD_ (v. 3.1.0)"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\tnew params.txt file created\n"
]
}
],
"source": [
"%%bash\n",
"## new params file\n",
"pyrad -n"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%%bash\n",
"## add phy and nex outputs\n",
"sed -i '/## 1. /c\\simulations/Tbal_full/ ## 1. working dir ' params.txt\n",
"sed -i '/## 2. /c\\simulations/Tbal_full/*.gz ## 2. data loc ' params.txt\n",
"sed -i '/## 3. /c\\simulations/Tbal_full/*barcodes.txt ## 3. Bcode ' params.txt\n",
"sed -i '/## 6. /c\\TGCAG ## 6. Cutter ' params.txt\n",
"sed -i '/## 7. /c\\4 ## 7. Nproc ' params.txt\n",
"sed -i '/## 10. /c\\.82 ## 10. clust thresh' params.txt\n",
"sed -i '/## 11. /c\\rad ## 11. datatype ' params.txt\n",
"sed -i '/## 13. /c\\6 ## 13. maxSH' params.txt\n",
"sed -i '/## 14. /c\\Tbal ## 14. outname' params.txt\n",
"sed -i '/## 16. /c\\ ## 16. addon taxa' params.txt\n",
"sed -i '/## 24./c\\99 ## 24. maxH' params.txt\n",
"sed -i '/## 30./c\\n,p,s ## 30. out format' params.txt"
]
},
{
"cell_type": "code",
"execution_count": 99,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%%bash\n",
"pyrad -p params.txt -q"
]
},
{
"cell_type": "code",
"execution_count": 139,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%%bash\n",
"sed -i '/## 1. /c\\simulations/Tbal_mut ## 1. working dir ' params.txt\n",
"sed -i '/## 2. /c\\simulations/Tbal_mut/*.gz ## 2. data loc ' params.txt\n",
"sed -i '/## 3. /c\\simulations/Tbal_mut/*barcodes.txt ## 3. Bcode ' params.txt"
]
},
{
"cell_type": "code",
"execution_count": 140,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%%bash\n",
"pyrad -p params.txt -q"
]
},
{
"cell_type": "code",
"execution_count": 148,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%%bash\n",
"sed -i '/## 1. /c\\simulations/Tbal_cov ## 1. working dir ' params.txt\n",
"sed -i '/## 2. /c\\simulations/Tbal_cov/*.gz ## 2. data loc ' params.txt\n",
"sed -i '/## 3. /c\\simulations/Tbal_cov/*barcodes.txt ## 3. Bcode ' params.txt"
]
},
{
"cell_type": "code",
"execution_count": 149,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%%bash\n",
"pyrad -p params.txt -q"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Functions for measuring shared data"
]
},
{
"cell_type": "code",
"execution_count": 141,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def getarray(loci, tree):\n",
" \"\"\" parse the loci list and return a \n",
" presence/absence matrix ordered by \n",
" the tips on the tree\"\"\"\n",
" ## order (ladderize) the tree\n",
" tree.ladderize()\n",
" ## get tip names\n",
" names = tree.get_leaf_names()\n",
" ## make empty matrix\n",
" lxs = np.zeros((len(names), len(loci)))\n",
" ## fill the matrix\n",
" for loc in xrange(len(loci)):\n",
" for seq in loci[loc].split(\"\\n\"):\n",
" if \">\" in seq:\n",
" tname = seq.split()[0][1:-1]\n",
" lxs[names.index(tname),loc] += 1\n",
" return lxs"
]
},
{
"cell_type": "code",
"execution_count": 142,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def countmatrix(lxsabove, lxsbelow, max=0):\n",
" \"\"\" fill a matrix with pairwise data sharing\n",
" between each pair of samples. You could put\n",
" in two different 'share' matrices to have\n",
" different results above and below the diagonal.\n",
" Can enter a max value to limit fill along diagonal.\n",
" \"\"\"\n",
" share = np.zeros((lxsabove.shape[0], \n",
" lxsbelow.shape[0]))\n",
" ## fill above\n",
" names = range(lxsabove.shape[0])\n",
" for row in lxsabove:\n",
" for samp1,samp2 in itertools.combinations(names,2):\n",
" shared = lxsabove[samp1, lxsabove[samp2,]>0].sum()\n",
" share[samp1,samp2] = shared\n",
" ## fill below\n",
" for row in lxsbelow:\n",
" for samp2,samp1 in itertools.combinations(names,2):\n",
" shared = lxsabove[samp1, lxsabove[samp2,]>0].sum()\n",
" share[samp1,samp2] = shared\n",
" ## fill diagonal\n",
" if not max:\n",
" for row in range(len(names)):\n",
" share[row,row] = lxsabove[row,].sum()\n",
" else:\n",
" for row in range(len(names)):\n",
" share[row,row] = max\n",
" return share"
]
},
{
"cell_type": "code",
"execution_count": 143,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(64, 883)\n",
"[[ 1. 1. 1. ..., 0. 1. 1.]\n",
" [ 1. 1. 1. ..., 0. 1. 1.]\n",
" [ 1. 1. 1. ..., 0. 1. 1.]\n",
" ..., \n",
" [ 1. 1. 1. ..., 1. 1. 1.]\n",
" [ 1. 1. 1. ..., 1. 1. 1.]\n",
" [ 1. 1. 1. ..., 1. 1. 1.]]\n"
]
}
],
"source": [
"locidata = open(\"simulations/Tbal_mut/outfiles/Tbal.loci\")\n",
"loci = locidata.read().split(\"|\\n\")[:-1]\n",
"lxs = getarray(loci, Tbal)\n",
"print lxs.shape\n",
"print lxs"
]
},
{
"cell_type": "code",
"execution_count": 144,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(64, 64)\n",
"[[ 773. 756. 733. ..., 689. 696. 691.]\n",
" [ 756. 770. 729. ..., 686. 693. 688.]\n",
" [ 733. 729. 759. ..., 678. 688. 682.]\n",
" ..., \n",
" [ 689. 686. 678. ..., 786. 757. 750.]\n",
" [ 696. 693. 688. ..., 757. 797. 776.]\n",
" [ 691. 688. 682. ..., 750. 776. 788.]]\n"
]
}
],
"source": [
"share = countmatrix(lxs, lxs)\n",
"print share.shape\n",
"print share"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plotting function"
]
},
{
"cell_type": "code",
"execution_count": 145,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"## Get ordered names of tips on the tree\n",
"tree = Tbal\n",
"tree.ladderize()\n",
"names = tree.get_leaf_names()\n",
"floater = [\"Taxon: %s\" % i for i in names]"
]
},
{
"cell_type": "code",
"execution_count": 146,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"