{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Collapse trees to shared clades for comparison" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The script prunes and collapses two or more trees to a shared set of clades. Note that the goals are \"**clades**\" rather than \"**taxa**\". A shared clade is defined as a clade present in all trees, with the same set of descendants. The resulting trees can then be used for back-to-back comparison (e.g., using [Dendroscope](http://dendroscope.org/)'s tanglegram function).\n", "\n", "Finding a shared set of clades is not trivial. In computer science it is a special form of the [covering problem](https://en.wikipedia.org/wiki/Covering_problems). While exploring the optimal solution is beyond the scope of this study, I designed and implemented a working heuristic to achieve the goal." ] }, { "cell_type": "code", "execution_count": 112, "metadata": {}, "outputs": [], "source": [ "import scipy as sp\n", "import matplotlib.pyplot as plt\n", "from skbio import TreeNode" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Number of clades to retain." ] }, { "cell_type": "code", "execution_count": 113, "metadata": {}, "outputs": [], "source": [ "clades_to_retain = 50" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Minimum number of descendants for a retained clade." ] }, { "cell_type": "code", "execution_count": 114, "metadata": {}, "outputs": [], "source": [ "clade_min_size = 50" ] }, { "cell_type": "code", "execution_count": 115, "metadata": {}, "outputs": [], "source": [ "trees = [\n", " TreeNode.read('astral.cons.nid.e5p68057.nwk'),\n", " TreeNode.read('concat.cons.nid.b50.nwk')\n", "]" ] }, { "cell_type": "code", "execution_count": 116, "metadata": {}, "outputs": [], "source": [ "node2taxa = []\n", "for tree in trees:\n", " n2t = {}\n", " for node in tree.postorder():\n", " if node.is_tip():\n", " n2t[node.name] = [node.name]\n", " else:\n", " n2t[node.name] = []\n", " for child in node.children:\n", " n2t[node.name].extend(n2t[child.name])\n", " node2taxa.append(n2t)" ] }, { "cell_type": "code", "execution_count": 117, "metadata": {}, "outputs": [], "source": [ "for i in range(len(node2taxa)):\n", " node2taxa[i] = {k: v for k, v in node2taxa[i].items() if len(v) >= clade_min_size}" ] }, { "cell_type": "code", "execution_count": 118, "metadata": {}, "outputs": [], "source": [ "for i in range(len(node2taxa)):\n", " for node in node2taxa[i]:\n", " node2taxa[i][node] = set(node2taxa[i][node])" ] }, { "cell_type": "code", "execution_count": 119, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Matching pairs: 273.\n" ] } ], "source": [ "matches = {}\n", "for node1, taxa1 in node2taxa[0].items():\n", " for node2, taxa2 in node2taxa[1].items():\n", " if taxa1 == taxa2:\n", " matches[','.join((node1, node2))] = taxa1\n", "print('Matching pairs: %d.' % len(matches))" ] }, { "cell_type": "code", "execution_count": 120, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.764080433286628e+55" ] }, "execution_count": 120, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sp.special.comb(len(matches), clades_to_retain)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is my heuristic for the covering problem. It starts at a given size of clade: *k*, and progressively recruits clades from left and right, until the desired number of clades: *n*, is reached." ] }, { "cell_type": "code", "execution_count": 121, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Taxa under matching pairs: 10575.\n" ] } ], "source": [ "universe = set().union(*[v for k, v in matches.items()])\n", "total_size = len(universe)\n", "print('Taxa under matching pairs: %d.' % total_size)" ] }, { "cell_type": "code", "execution_count": 122, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "212" ] }, "execution_count": 122, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mean_clade_size = round(len(universe) / clades_to_retain)\n", "mean_clade_size" ] }, { "cell_type": "code", "execution_count": 123, "metadata": { "scrolled": true }, "outputs": [], "source": [ "def select_clades(start_size):\n", " res = []\n", " for match, taxa in sorted(matches.items(), key=lambda x: abs(start_size - len(x[1]))):\n", " is_unique = True\n", " for m in res:\n", " if not taxa.isdisjoint(matches[m]):\n", " is_unique = False\n", " break\n", " if is_unique is True:\n", " res.append(match)\n", " if len(res) == clades_to_retain:\n", " break\n", " return res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test multiple starting sizes." ] }, { "cell_type": "code", "execution_count": 124, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting at 212, covering 6754 taxa, CV = 0.481\n", "Starting at 213, covering 6754 taxa, CV = 0.481\n", "Starting at 214, covering 6754 taxa, CV = 0.481\n", "Starting at 215, covering 6754 taxa, CV = 0.481\n", "Starting at 216, covering 6754 taxa, CV = 0.481\n", "Starting at 217, covering 6754 taxa, CV = 0.481\n", "Starting at 218, covering 6754 taxa, CV = 0.481\n", "Starting at 219, covering 6754 taxa, CV = 0.481\n", "Starting at 220, covering 6754 taxa, CV = 0.481\n", "Starting at 221, covering 6754 taxa, CV = 0.481\n", "Starting at 222, covering 6754 taxa, CV = 0.481\n", "Starting at 223, covering 6808 taxa, CV = 0.487\n", "Starting at 224, covering 6808 taxa, CV = 0.487\n", "Starting at 225, covering 6808 taxa, CV = 0.487\n", "Starting at 226, covering 6808 taxa, CV = 0.487\n", "Starting at 227, covering 6808 taxa, CV = 0.487\n", "Starting at 228, covering 6904 taxa, CV = 0.516\n", "Starting at 229, covering 6904 taxa, CV = 0.516\n", "Starting at 230, covering 6904 taxa, CV = 0.516\n", "Starting at 231, covering 6904 taxa, CV = 0.516\n", "Starting at 232, covering 6904 taxa, CV = 0.516\n", "Starting at 233, covering 6904 taxa, CV = 0.516\n", "Starting at 234, covering 6904 taxa, CV = 0.516\n", "Starting at 235, covering 6904 taxa, CV = 0.516\n", "Starting at 236, covering 6904 taxa, CV = 0.516\n", "Starting at 237, covering 6904 taxa, CV = 0.516\n", "Starting at 238, covering 7003 taxa, CV = 0.539\n", "Starting at 239, covering 7003 taxa, CV = 0.539\n", "Starting at 240, covering 7003 taxa, CV = 0.539\n", "Starting at 241, covering 7003 taxa, CV = 0.539\n", "Starting at 242, covering 7003 taxa, CV = 0.539\n", "Starting at 243, covering 7095 taxa, CV = 0.550\n", "Starting at 244, covering 7168 taxa, CV = 0.557\n", "Starting at 245, covering 7168 taxa, CV = 0.557\n", "Starting at 246, covering 7168 taxa, CV = 0.557\n", "Starting at 247, covering 7168 taxa, CV = 0.557\n", "Starting at 248, covering 7168 taxa, CV = 0.557\n", "Starting at 249, covering 7168 taxa, CV = 0.557\n", "Starting at 250, covering 7168 taxa, CV = 0.557\n", "Starting at 251, covering 7168 taxa, CV = 0.557\n", "Starting at 252, covering 7168 taxa, CV = 0.557\n", "Starting at 253, covering 7178 taxa, CV = 0.558\n", "Starting at 254, covering 7178 taxa, CV = 0.558\n", "Starting at 255, covering 7178 taxa, CV = 0.558\n", "Starting at 256, covering 7178 taxa, CV = 0.558\n", "Starting at 257, covering 7178 taxa, CV = 0.558\n", "Starting at 258, covering 7178 taxa, CV = 0.558\n", "Starting at 259, covering 7178 taxa, CV = 0.558\n", "Starting at 260, covering 7178 taxa, CV = 0.558\n", "Starting at 261, covering 7178 taxa, CV = 0.558\n", "Starting at 262, covering 7178 taxa, CV = 0.558\n", "Starting at 263, covering 7178 taxa, CV = 0.558\n", "Starting at 264, covering 7178 taxa, CV = 0.558\n", "Starting at 265, covering 7178 taxa, CV = 0.558\n", "Starting at 266, covering 7178 taxa, CV = 0.558\n", "Starting at 267, covering 7178 taxa, CV = 0.558\n", "Starting at 268, covering 7178 taxa, CV = 0.558\n", "Starting at 269, covering 7178 taxa, CV = 0.558\n", "Starting at 270, covering 7178 taxa, CV = 0.558\n", "Starting at 271, covering 7178 taxa, CV = 0.558\n", "Starting at 272, covering 7178 taxa, CV = 0.558\n", "Starting at 273, covering 7178 taxa, CV = 0.558\n", "Starting at 274, covering 7178 taxa, CV = 0.558\n", "Starting at 275, covering 7183 taxa, CV = 0.559\n", "Starting at 276, covering 7183 taxa, CV = 0.559\n", "Starting at 277, covering 7183 taxa, CV = 0.559\n", "Starting at 278, covering 7292 taxa, CV = 0.575\n", "Starting at 279, covering 7292 taxa, CV = 0.575\n", "Starting at 280, covering 7292 taxa, CV = 0.575\n", "Starting at 281, covering 7292 taxa, CV = 0.575\n", "Starting at 282, covering 7292 taxa, CV = 0.575\n", "Starting at 283, covering 7348 taxa, CV = 0.586\n", "Starting at 284, covering 7403 taxa, CV = 0.614\n", "Starting at 285, covering 7417 taxa, CV = 0.616\n", "Starting at 286, covering 7417 taxa, CV = 0.616\n", "Starting at 287, covering 7417 taxa, CV = 0.616\n", "Starting at 288, covering 7417 taxa, CV = 0.616\n", "Starting at 289, covering 7417 taxa, CV = 0.616\n", "Starting at 290, covering 7472 taxa, CV = 0.645\n", "Starting at 291, covering 7472 taxa, CV = 0.645\n", "Starting at 292, covering 7541 taxa, CV = 0.652\n", "Starting at 293, covering 7541 taxa, CV = 0.652\n", "Starting at 294, covering 7541 taxa, CV = 0.652\n", "Starting at 295, covering 7555 taxa, CV = 0.654\n", "Starting at 296, covering 7555 taxa, CV = 0.654\n", "Starting at 297, covering 7555 taxa, CV = 0.654\n", "Starting at 298, covering 7555 taxa, CV = 0.654\n", "Starting at 299, covering 7555 taxa, CV = 0.654\n", "Starting at 300, covering 7611 taxa, CV = 0.660\n", "Starting at 301, covering 7611 taxa, CV = 0.660\n", "Starting at 302, covering 7611 taxa, CV = 0.660\n", "Starting at 303, covering 7611 taxa, CV = 0.660\n", "Starting at 304, covering 7611 taxa, CV = 0.660\n", "Starting at 305, covering 7611 taxa, CV = 0.660\n", "Starting at 306, covering 7611 taxa, CV = 0.660\n", "Starting at 307, covering 7611 taxa, CV = 0.660\n", "Starting at 308, covering 7611 taxa, CV = 0.660\n", "Starting at 309, covering 7611 taxa, CV = 0.660\n", "Starting at 310, covering 7639 taxa, CV = 0.664\n", "Starting at 311, covering 7649 taxa, CV = 0.665\n", "Starting at 312, covering 7649 taxa, CV = 0.665\n", "Starting at 313, covering 7649 taxa, CV = 0.665\n", "Starting at 314, covering 7649 taxa, CV = 0.665\n", "Starting at 315, covering 7649 taxa, CV = 0.665\n", "Starting at 316, covering 7649 taxa, CV = 0.665\n", "Starting at 317, covering 7649 taxa, CV = 0.665\n", "Starting at 318, covering 7649 taxa, CV = 0.665\n", "Starting at 319, covering 7649 taxa, CV = 0.665\n", "Starting at 320, covering 7649 taxa, CV = 0.665\n", "Starting at 321, covering 7649 taxa, CV = 0.665\n", "Starting at 322, covering 7649 taxa, CV = 0.665\n", "Starting at 323, covering 7649 taxa, CV = 0.665\n", "Starting at 324, covering 7649 taxa, CV = 0.665\n", "Starting at 325, covering 7649 taxa, CV = 0.665\n", "Starting at 326, covering 7649 taxa, CV = 0.665\n", "Starting at 327, covering 7649 taxa, CV = 0.665\n", "Starting at 328, covering 7649 taxa, CV = 0.665\n", "Starting at 329, covering 7649 taxa, CV = 0.665\n", "Starting at 330, covering 7697 taxa, CV = 0.671\n", "Starting at 331, covering 7697 taxa, CV = 0.671\n", "Starting at 332, covering 7697 taxa, CV = 0.671\n", "Starting at 333, covering 7697 taxa, CV = 0.671\n", "Starting at 334, covering 7699 taxa, CV = 0.672\n", "Starting at 335, covering 7699 taxa, CV = 0.672\n", "Starting at 336, covering 7701 taxa, CV = 0.672\n", "Starting at 337, covering 7701 taxa, CV = 0.672\n", "Starting at 338, covering 7701 taxa, CV = 0.672\n", "Starting at 339, covering 7725 taxa, CV = 0.675\n", "Starting at 340, covering 7725 taxa, CV = 0.675\n", "Starting at 341, covering 7725 taxa, CV = 0.675\n", "Starting at 342, covering 7725 taxa, CV = 0.675\n", "Starting at 343, covering 7725 taxa, CV = 0.675\n", "Starting at 344, covering 7725 taxa, CV = 0.675\n", "Starting at 345, covering 7725 taxa, CV = 0.675\n", "Starting at 346, covering 7725 taxa, CV = 0.675\n", "Starting at 347, covering 7725 taxa, CV = 0.675\n", "Starting at 348, covering 7764 taxa, CV = 0.681\n", "Starting at 349, covering 7764 taxa, CV = 0.681\n", "Starting at 350, covering 7764 taxa, CV = 0.681\n", "Starting at 351, covering 7764 taxa, CV = 0.681\n", "Starting at 352, covering 7764 taxa, CV = 0.681\n", "Starting at 353, covering 7764 taxa, CV = 0.681\n", "Starting at 354, covering 7764 taxa, CV = 0.681\n", "Starting at 355, covering 7764 taxa, CV = 0.681\n", "Starting at 356, covering 7764 taxa, CV = 0.681\n", "Starting at 357, covering 7764 taxa, CV = 0.681\n", "Starting at 358, covering 7764 taxa, CV = 0.681\n", "Starting at 359, covering 7764 taxa, CV = 0.681\n", "Starting at 360, covering 7764 taxa, CV = 0.681\n", "Starting at 361, covering 7764 taxa, CV = 0.681\n", "Starting at 362, covering 7764 taxa, CV = 0.681\n", "Starting at 363, covering 7764 taxa, CV = 0.681\n", "Starting at 364, covering 7764 taxa, CV = 0.681\n", "Starting at 365, covering 7764 taxa, CV = 0.681\n", "Starting at 366, covering 7764 taxa, CV = 0.681\n", "Starting at 367, covering 7764 taxa, CV = 0.681\n", "Starting at 368, covering 7764 taxa, CV = 0.681\n", "Starting at 369, covering 7764 taxa, CV = 0.681\n", "Starting at 370, covering 7764 taxa, CV = 0.681\n", "Starting at 371, covering 7764 taxa, CV = 0.681\n", "Starting at 372, covering 7764 taxa, CV = 0.681\n", "Starting at 373, covering 7764 taxa, CV = 0.681\n", "Starting at 374, covering 7764 taxa, CV = 0.681\n", "Starting at 375, covering 7764 taxa, CV = 0.681\n", "Starting at 376, covering 7764 taxa, CV = 0.681\n", "Starting at 377, covering 7764 taxa, CV = 0.681\n", "Starting at 378, covering 7764 taxa, CV = 0.681\n", "Starting at 379, covering 7764 taxa, CV = 0.681\n", "Starting at 380, covering 7764 taxa, CV = 0.681\n", "Starting at 381, covering 7764 taxa, CV = 0.681\n", "Starting at 382, covering 7764 taxa, CV = 0.681\n", "Starting at 383, covering 7764 taxa, CV = 0.681\n", "Starting at 384, covering 7764 taxa, CV = 0.681\n", "Starting at 385, covering 7764 taxa, CV = 0.681\n", "Starting at 386, covering 7817 taxa, CV = 0.692\n", "Starting at 387, covering 7817 taxa, CV = 0.692\n", "Starting at 388, covering 7817 taxa, CV = 0.692\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Starting at 389, covering 7817 taxa, CV = 0.692\n", "Starting at 390, covering 7817 taxa, CV = 0.692\n", "Starting at 391, covering 7817 taxa, CV = 0.692\n", "Starting at 392, covering 7817 taxa, CV = 0.692\n", "Starting at 393, covering 7817 taxa, CV = 0.692\n", "Starting at 394, covering 7817 taxa, CV = 0.692\n", "Starting at 395, covering 7869 taxa, CV = 0.718\n", "Starting at 396, covering 7869 taxa, CV = 0.718\n", "Starting at 397, covering 7869 taxa, CV = 0.718\n", "Starting at 398, covering 7869 taxa, CV = 0.718\n", "Starting at 399, covering 7869 taxa, CV = 0.718\n", "Starting at 400, covering 7869 taxa, CV = 0.718\n", "Starting at 401, covering 7983 taxa, CV = 0.757\n", "Starting at 402, covering 7983 taxa, CV = 0.757\n", "Starting at 403, covering 7983 taxa, CV = 0.757\n", "Starting at 404, covering 7983 taxa, CV = 0.757\n", "Starting at 405, covering 7983 taxa, CV = 0.757\n", "Starting at 406, covering 7983 taxa, CV = 0.757\n", "Starting at 407, covering 7983 taxa, CV = 0.757\n", "Starting at 408, covering 7983 taxa, CV = 0.757\n", "Starting at 409, covering 7983 taxa, CV = 0.757\n", "Starting at 410, covering 7983 taxa, CV = 0.757\n", "Starting at 411, covering 7983 taxa, CV = 0.757\n", "Starting at 412, covering 7983 taxa, CV = 0.757\n" ] } ], "source": [ "x, y, z = [], [], []\n", "selections = []\n", "start_size = mean_clade_size\n", "while True:\n", " selected_matches = select_clades(start_size)\n", " if len(selected_matches) < clades_to_retain:\n", " break\n", " selections.append(selected_matches)\n", " covered = set().union(*[matches[x] for x in selected_matches])\n", " cv = sp.stats.variation([len(matches[x]) for x in selected_matches])\n", " x.append(start_size)\n", " y.append(len(covered) / total_size)\n", " z.append(cv)\n", " print('Starting at %d, covering %d taxa, CV = %.3f' % (start_size, len(covered), cv))\n", " start_size += 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot them out." ] }, { "cell_type": "code", "execution_count": 134, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEKCAYAAAD+XoUoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYFNXZ9/Hv7bC5gQiYGBZBg7uAMm6ocUdiFBSjwqMGNYZoombTxMTEBfNGjTEuCWrUEHcw7iRqfFzAjSgMiAooCogwIjIzbCrIer9/nOpnimGG6Rlmqnr5fa6rr+46darr7mK4+/SpqnPM3RERkeKwRdoBiIhIcpT0RUSKiJK+iEgRUdIXESkiSvoiIkVESV9EpIgo6YuIFBElfRGRIqKkLyJSRFqkHUBNHTt29O7du6cdhohIXpk8eXKlu3eqr17OJf3u3btTVlaWdhgiInnFzD7Opp66d0REioiSvohIEVHSFxEpIkr6IiJFRElfRKSIKOmLiBQRJX0RkSKSc9fpi4jkmwUL4O67Ye3azXufLl1g+PCmiakuSvoiIpvpvrtXc+WVrTDbvDnHD9xjOcOHt2uiqGqXVfeOmQ0ws5lmNsvMLqtl/U1mNjV6fGBmS2Pr1sXWjW3K4EVEcsHSFyfTmq9Y71ts1uO/2x7X7LHW29I3sxJgJHAsUA5MMrOx7j4jU8fdfxarfxGwb+wtVrp7n6YLWUQktyyft5S2W3wJU6Zu3htttVXTBLQJ2XTvHADMcvc5AGY2BhgEzKij/lDgyqYJT0Qk9y1ftJK2W62B3r3TDqVe2XTvdAbmx5bLo7KNmNlOQA/gpVhxGzMrM7M3zOykRkcqIpKLKipYvqIFbdtZ2pFkJZuWfm2fpK6zFUOAR919Xaysm7svMLOdgZfM7F13n73BDsyGA8MBunXrlkVIIiI5YsoUltGOth1bpR1JVrJp6ZcDXWPLXYAFddQdAoyOF7j7guh5DjCeDfv7M3XudPdSdy/t1Kne4aBFRHLHlCkspy1td9wm7Uiykk3SnwT0NLMeZtaKkNg3ugrHzHYD2gP/jZW1N7PW0euOwCHUfS5ARCT/TJ7M8hbb07ZDy7QjyUq93TvuvtbMLgSeA0qAUe4+3cxGAGXunvkCGAqMcfd4188ewN/MbD3hC+a6+FU/IiI578034d57614/fjzLt2hP27bJhbQ5sro5y92fAZ6pUXZFjeWratluArDPZsQnIpIed7jgApgxgzqzeosWLF+/Ne2a956qJqM7ckVE6vLaa/DWW/C3v9U5PsKqVbC6Td3fCblGSV9EiseXX8KaNdnXv+kmaN8ezjyzzirLl4dnJX0RkVwyfjwcdVTosmmISy/d5J2yy5aFZyV9EZFccuON0LEj/OY32W/TsuUmW/mglr6ISO6ZMweefhouvxx++tMmfWslfRHJL+vXw+jRMH9+/XXz1auvwhZbwPnnN/lbK+mLSH556CE466y0o2h+55wDnWsdNmyzKOmLSP744gv41a+gb1945ZXQGi5UrVs3y9sq6YtIOtauhYsvDl0Z2fr88zDX3yOPJDKWeyHKJH3dnCUiTW/9+pBl3MO1gkuWVF+CeOutYbiAAQMalsAvuQT69WueeIvAsmXhIp9m+iHR5JT0RXLRvfeGZLx+fXVZJtHHy2q6+mq44oq610uTW748dO1Yfgynr6QvknNWr4bf/S5cU37MMRuua9cOtt8+9L1vuy106FDdD9+hg1rsKcgk/XyhpC+SazKXTz7zDHz722lHI/VQ0hcpJKtWhST87LOwbl399ZvChAnQq1fom5ecp6QvUigqK0N3yYcfQrduoTslCR07wvXX508ncZFbvhx23DHtKLKnpC+Fb/lyWLGiYdu4hzFX5s2DsWPhhBOUhKVWy5bBbrulHUX2lPSlsL3/fugqachwunH33AMnntikIUlhUfeOSC65447w/Ne/QklJw7bt3l396sIdd8Bll9U9IvPy5flzYxYo6UshW7kS7rsPBg+GH/847WgkT738crgqdtiw2tdvsQV8//vJxrQ5lPQl/330EUyfvnH5m2+GO1Z/+MPkY5KCUVkJu+8eJtEqBEr6kl9Wr4Z//Ss8A7z9Nvz5z3X32e+xBxxxRGLhSeGpqAgXbxUKJX3JL6NGwQUXbFh25pmh+6ZFLX/OPXroqhvZLBUVYRDSQqGkL/nlkUdg113DZZQAW28NXbqkG5MULPfQvdOpU9qRNB0lfckflZXhrNqvfpVfF0ZL3vr889CTWEhJv4BnTJCC89RTYSiEU05JOxIpEhUV4bljx3TjaEpq6UvumzUrjDr5xhuhj37ffdOOSIpEJumrpS+SpGuvhccfD/33l12mE7OSmMrK8Fx0Sd/MBpjZTDObZWaX1bL+JjObGj0+MLOlsXXDzOzD6FHH7Q0idViyJIxyOWwYTJsGw4enHZEUkUJs6dfbvWNmJcBI4FigHJhkZmPdfUamjrv/LFb/ImDf6PX2wJVAKeDA5GjbJU36KaRw3X9/uLO25mWaIgko1j79A4BZ7j4HwMzGAIOAGXXUH0pI9ADHAc+7++Jo2+eBAcDozQlaCszTT4epARcs2HjdihVw4IHqx5dUVFRAmzahZ7FQZJP0OwPzY8vlwIG1VTSznYAewEub2LZzw8OUgnXbbeHGqj32gHPPrb3OmWcmG5NIJHONfiGdRsom6df2cesYb44hwKPunpliKKttzWw4MBygWyHd7yybtn49/OlPYaKSceOgVau0IxLZQEVFYfXnQ3YncsuBrrHlLkAtv8OBkPTjXTdZbevud7p7qbuXdiq0Iyx1e/75MFjaRRcp4UtOqqgorP58yC7pTwJ6mlkPM2tFSOxja1Yys92A9sB/Y8XPAf3NrL2ZtQf6R2UiYaDyTp3g5JPTjkSkVkXZ0nf3tcCFhGT9HvBPd59uZiPMbGCs6lBgjHv1VAPRCdxrCF8ck4ARmZO6UuQ++SSMlnnOOdC6ddrRiNSq0MbdgSzvyHX3Z4BnapRdUWP5qjq2HQWMamR8Uqj+/vcwpIKuu5eEuMPs2eEK4GysWRPG3inKpC/SpNauhbvugv79YZdd0o5GisTLL8ORRzZ8u84Fdr2hkr40n/feC4Ok1VReHh633JJ8TFK0XnstPI8eDS1bZrdNq1ahbVJIlPSl+fzyl/Dvf9e+rmdPOPHEZOORojZ5cpiKYciQtCNJlwZck+bhDhMnwllnhU7Umo/338++uSXSBKZMgf32SzuK9CnpS/OYPx8WLYKDDgr3sdd8bKE/PUlOZSXMm1dY0x42lv7nSfMoKwvPpaXpxiFC6NoBJX1Q0pfmMmlS6L7p3TvtSESYMiU8a9w+nciV5jJpEvTqpRuv8oA73Hln9TDChejxx8PVwdttl3Yk6VPSl6azbFm4o8U9dO8U+2USeWL8eDj//LSjaH4XX5x2BLlBST/XfPlluHkpW4sWwe23w9tvN19M9XGHuXPD4GlxBxyQSjjSMPffD9tuG0bG2HLLtKNpPi2U7QAl/dzyxBMweHDDt2vRAvbfH0pKmj6mbO23H/zwh9WzTbRuDUOHphePZGXFCnj0Ufjud0Pil8KnpJ9LHnwQdtghTP6drVat4KSTCu9ecWmU8vLqC6eyMXlyGF/mrLOaLybJLUr6ueKrr+A//wmzRP3sZ/XXF6mhqipcnVJZ2bDtdt4ZDj+8eWKS3KOknyvGjQv9+QMH1l9XpBa//CUsWRKmHP7GN7LfrmtX3StXTJT0c8VTT4X+8KOOSjsSyQHz58NDD4UZJTPWr4eFC8Nj/Xr47LMwl3xmBou5c0PiP/74VEKWPKGknwvcw8BkAwaEIQqk6F1xBdxzz8blbduGVnxJSZjGr1+/6vP3p50GV16ZaJiSh5T0c8Hs2eF6uWOOSTsSyQFr1oQffmecEeaaiWvVCszSiUsKg5J+LnjllfCss2lCOL2zZAmceqpuaJamp9M3ueDll8OcbLvvnnYkkgMeeyyc3im0yTskN6ilnwteeQW+9S39bs9BX30FI0fCI4+EKX2TMGNGmF+mkO+OlfQo6adt3rxw2YWuzc85X34JBx4I06eH5x12SGa/O+4Iv/hFMvuS4qOkn5THH4dLLtnwGjwIs0iB+vNz0O9+FxL+Y481bnQMkVykpJ+UBx8Mo1DWNi9s586wzz7Jx1QEVq0Ko1osWdKw7datC/9kF1yghC+FRUk/Ce4wYUK4a6a2i6+l2Tz9NNx8c/hebegoi0cdBddd1zxxiaRFST8JH30UbqPs1y/tSIrOE09Ahw7htImG1hXRJZvJmDAhPCvpJ2rNmnCj84knKuGLZCjpJ2HChDBY+d57px1JURk/HpYuhZNPTjsSkdyRVfvHzAYAtwAlwN3uvlFPp5mdBlwFOPC2u/9PVL4OeDeqNs/d83cYyZUrQz9BQ73yChx0ULqTnBSoJUvg0kvDlLw1LVoUbnI69tjk4xLJVfUmfTMrAUYCxwLlwCQzG+vuM2J1egK/Bg5x9yVmFr+ieaW792niuNNx2mmhv6Cx28pmmzULfvSjcFUOwAcfhPHjjzsOWrbcsO7OO4fhjHSTk0i1bFr6BwCz3H0OgJmNAQYBM2J1fgCMdPclAO6+qKkDTd2aNWFQlBNPDCNhNURJSchKstnuuQdefDHcwAxQWhpGliwtTTUskbyRTdLvDMyPLZcDB9aosyuAmb1O6AK6yt3/E61rY2ZlwFrgOnd/cvNCTsnbb4dbNM84A04/Pe1oitazz8Ihh4TvXxFpuGySfm0Dwngt79MTOALoArxqZnu7+1Kgm7svMLOdgZfM7F13n73BDsyGA8MBunXr1sCPkJDXXgvPhxySbhxFbOFCmDIF/vCHtCMRyV/ZXL1TDnSNLXcBFtRS5yl3X+PuHwEzCV8CuPuC6HkOMB7Yt+YO3P1Ody9199JOnTo1+EMk4vXXoXt36NIl7UiK1n+i346aGUqk8bJJ+pOAnmbWw8xaAUOAsTXqPAkcCWBmHQndPXPMrL2ZtY6VH8KG5wLyg3to6auVn4ply8JdtbffHmaN6tUr7YhE8le93TvuvtbMLgSeI/TXj3L36WY2Aihz97HRuv5mNgNYB1zq7lVm1g/4m5mtJ3zBXBe/6ifnzJ8Pv/0trF69Yfnq1aFv4dBD04mryD3ySPUgpL/4hUagFtkcWV2n7+7PAM/UKLsi9tqBn0ePeJ0JQP6MJPbHP4ZRtnbZZeN1++4LJ5yQfExCRUV4rqqC7bdPNxaRfKeb0zO++ALuvReGDoX77087GompqgrX2ivhi2w+DcOQ8cAD8Pnn4c4fySlVVWHQNBHZfIXT0q+ogK99rfHbu4cunIMOarqYpEksXqykL9JUCifpb7VVOAm7OU4+WWcJc5D68kWaTuEk/a23hhEj0o5CmkFVlQYoFWkq6tOXnKfuHZGmo6QvOc09JH1174g0DSV9yWmffw5r16qlL9JUlPQlp1VVhWe19EWahpK+5LRM0ldLX6RpKOlLTlu8ODwr6Ys0DSV9yWnq3hFpWkr6ktPUvSPStJT0JadlunfU0hdpGkr6ktOqqqBtW2hROPeOi6RKSV9ymkbYFGlaSvqS0zQEg0jTUtKXnKYRNkWalnpKJXWrVsFNN8G8eRuv++AD+Pa3k49JpFAp6UuqPvsMBg6EiROhY8eNpzNo1QqOPDKd2EQKkZK+pGrkSCgrg0cfhVNOSTsakcKnPn1J1cKFsMMOSvgiSVHSl1RVVoZuHRFJhpK+pErX4YskS0lfUqWWvkiylPQlVWrpiyRLSV9S4x6Svlr6IsnJKumb2QAzm2lms8zssjrqnGZmM8xsupk9FCsfZmYfRo9hTRW45L/lyzX/rUjS6r1O38xKgJHAsUA5MMnMxrr7jFidnsCvgUPcfYmZ7RCVbw9cCZQCDkyOtl3S9B9F8k1lZXhWS18kOdm09A8AZrn7HHdfDYwBBtWo8wNgZCaZu/uiqPw44Hl3Xxytex4Y0DShS77TBCkiycsm6XcG5seWy6OyuF2BXc3sdTN7w8wGNGBbzGy4mZWZWVlFRUX20UteU0tfJHnZJH2rpcxrLLcAegJHAEOBu81suyy3xd3vdPdSdy/t1KlTFiFJIVBLXyR52ST9cqBrbLkLsKCWOk+5+xp3/wiYSfgSyGZbKVJq6YskL5ukPwnoaWY9zKwVMAQYW6POk8CRAGbWkdDdMwd4DuhvZu3NrD3QPyoToaoKSkqgXbu0IxEpHvVevePua83sQkKyLgFGuft0MxsBlLn7WKqT+wxgHXCpu1cBmNk1hC8OgBHuvrg5Pojkn8rKMEFKzeGURaT5mPtGXeypKi0t9bKysrTDkASceipMnw4zZtRfV0Q2zcwmu3tpffV0R66kprJSJ3FFkqakL6nREAwiyVPSl9SopS+SPCV9SYUGWxNJh5K+pKKyElavhq9/Pe1IRIqLkr6k4u23w3OvXunGIVJslPQlFVOnhufevdONQ6TYKOlLKqZOha5ddSJXJGlK+pKKqVOhT5+0oxApPkr6kriVK+H995X0RdKgpC+Jmz4d1q1T0hdJQ70Drok0VkUFTJiwcfn48eFZSV8keUr60mx+/nN44IHa1+2wA3Tvnmg4IoKSvjSj996Dfv3gr3/deN03vgFbqHNRJHFK+tJs5syB00+HffdNOxIRyVBbS5rFkiXhscsuaUciInFK+tIs5swJzzvvnG4cIrIhJX1pFpmkr5a+SG5R0pdmoZa+SG5S0pdmMXs2dOoE226bdiQiEqekL81izhy18kVykZK+NAslfZHcpKQvTW7NGpg3TydxRXKRbs4qAK+/DpMmpbd/d5g4EZ55Br76KiyvW6eWfl6qqICTT4alS+H3v4eTTgrlgwbB7beHW6lrevJJ2HVX2HPPZGNN29y5cMIJMG1a2pE0iJJ+Dpk7F0aMCC3lbM2fDy+/3GwhZa1DBzj11HDyFqBNm5A7JM+MHg3DhsGQITBgQEj6//oX7Ldf7QkfQtI/4YTCSfrr1kFJSdpRNBsl/Rxy111wzz3Qo0f227RuDX/8I5x9NrRs2VyR1W+bbaCF/pryX8uWYcKDVavC4Ehr18LNN4fEX5sJE2Ds2NDy+P3v4bHH4KWX4M47w8z33/wm3H8/bLVV+LVwyinwve/B3/4Gr7wCDz4Y/vBrqx/3xRdw0UVQVgZmcOWV4b1Gj4Y//CH8vPzOd+D668Mvko8+Cv8xIPynmjwZ/vKXMALgrbeGfR14INx2W0jw22wTRgh87jm48UbYcsuw/MUX0LFjeI8ddwzvc+65Ib5DD23Of4nm4+459ejbt68Xq9JS90MPTTsKKWpLl7off7x7377uL7zgfsst7vfcs+lthg1zf+SR6uXKyurXl1/ufuut4fXChe677OL+yivuPXu6V1Vtun7cL3/p/pOfVC8vXuz+ySfuXbu6L1rkvmaN+5FHuj/xRFjeZZfqugMGuL/6qvuMGe4nnOC+enUov+AC93vvDa/B/eGHw+vVq90PPji8j7v7mDHu55wTXu+zj/v48eH1JZe477XXpo9NgoAyzyLHZtU2M7MBwC1ACXC3u19XY/3ZwA3AJ1HRX9397mjdOuDdqHyeuw/c/K+qwlNZGRoRV1+ddiRS1Nq1g6efDq+XLAkt58cfhx/8ICz/4hdw8MGbfo9p0+C3vw3nBb74Ao47LpR/7Wuh//LII+GJJ2D77TddP+6FF2DMmOrl9u3DL4UjjqjuUzzjjFB20knhhNIbb0DPnjBzJhxyCIwcGf6T7b9/qL9yZRjjG0Jr/5RTwuuZM0NMxx4bltetC638ZctCjIcfHsrPOguefTbrQ5sr6k36ZlYCjASOBcqBSWY21t1n1Kj6sLtfWMtbrHR3TZdRjxdfDL9Q+/dPOxKRyIgRcPnloQulb1/4n/8JXTTjxm16u7PPDv38vXuHbpHMrDkA774bTgAtWJBd/Qz30K1Ts6wup58O//wn7L57OLlkFuoPGwbXXrtx/TZtqvvx3WGvveC//92wztKlG8eQh7K5ZPMAYJa7z3H31cAYYFDzhlV8/vd/YbvtoLQ07UhEgA8/DIn58MNhxYrQv28WLs+qadtt4fPPq5c//zy0jNesCX32GRMnhpbxW2/Bn/4U+t03VT+uf/8NJ2ZYsiT0yb/8cviZvG5d+HLKtMIHDw5fJKNHhy8AgKOPhkcfhUWLwvLixfDxxxvva7fdwlVMmaS/Zk2Y43O77cIvoddeC+V1xZrjskn6nYH5seXyqKymU8zsHTN71My6xsrbmFmZmb1hZifVtgMzGx7VKauoqMg++gLywgvhb7KALxqQfHL55eHELMDQoaEFftBBcMklG9cdMgRuuCFMnDB7NlxzTUjIxx4bWtoQTgz/4AcwalS4CujGG8MJUffa69f029+GRL/33uEXwbhx4Yvi2mtDd1Hv3uEKo0FRe7R9+3A10ccfwwEHhLI99wyfqX9/6NUr7O/TTzfeV6tW4cvhV78K79unT/W8n//4B/z4x6GLa8stG31402S+qZ9IgJmdChzn7udFy2cBB7j7RbE6HYAv3H2VmZ0PnObuR0XrvuHuC8xsZ+Al4Gh3n13X/kpLS72srGyzP1g+WbAAOncOF0n85CdpRyMi+cjMJrt7vX0F2bT0y4F4y70LsCBewd2r3H1VtHgX0De2bkH0PAcYD2gepRoyN1ZlGiQiIs0lm6Q/CehpZj3MrBUwBBgbr2BmO8YWBwLvReXtzax19LojcAhQ8wRw0Zs4MVzj3kenu0WkmdV79Y67rzWzC4HnCJdsjnL36WY2gnBd6FjgYjMbCKwFFgNnR5vvAfzNzNYTvmCuq+Wqn6I3cWLoYszTLkIRySP19uknrdj69NevD5crDx0abiQUEWmMpuzTl2b04Yfhng/154tIEjRaSkL+9a8wlMf69RuWr1gRnpX0RSQJSvoJee45KC+H735343WdO8MeeyQfk4gUHyX9hHz2WRg98/77045ERIqZ+vQTsnBhGG9KRCRNSvoJWbgQvv71tKMQkWKnpJ+Qzz5TS19E0qekn4AVK8JAgmrpi0jalPQT8Nln4VktfWlOq1bBMceE4TwefhhefTUMC9+nD3zySe1XjsWddx7MaOT98uPHVw9E2dh4m9Idd8B99226zpNPbvh5r7gijHZb6HT1TgIWLgzPaulLc3rrrTD0+9SpYfn888NIyOecE5YffXTT2999d+P3PX58mGa2X7/st6kZb1NZuzZ89vrUnM99xIimjSNXqaWfgExLX0lf6nLffWH8pd69wyx8EIaCP/roUH700TBvXiivqAgz++2/f3i8/nqYF+TMM0MC7dMnzDv+z3+GRHbGGTB3bhiKHsJ8I5dcAvvsE977L38J5UccEeYdhzCpz8EHhyHqTz01zGII0L17mJN8v/3C9u+/H977jjvgppvCvl99dcPPtnhxmMGwV68wJP8772wc7+zYYOvvvbfhzYpz54ZtIXye/fcPn2X48OrJs444An7zmzCHyi23wFVXhXlaIMy7vv/+4dieckrobs3M537ppdX7P/vs6i/GF18M0wPss08Y9n/Vqro/f97JZiLdJB+FODH67beHeZc/+STtSCQXTZvmvuuu7hUVYTkzX/gJJ1TPSf73v7sPGhReDx0a5vl2d//4Y/fddw+vx41z/853qt83Pl/5Rx9Vz+F9223ugweHucTj+zv8cPdJk0Ichx3m/sUXofy669yvvjq83mmn6nnLR450//73w+srr3S/4YbaP9+FF7pfdVV4/eKL7r171x5vXO/e7rNnV+//mms2jNXd/cwz3ceOrY79gguq18XjqWve9ZrzuWeWV65079LFfebMUH7WWe433bTpz58LyHJidLX0E7BwYZhpLjN/s0jcSy+F/vaOHcNyZr7w//43TEsLofWfmaXvhRfgwgtDC3XgQFi+fMPZCuvzwguh+6NFiw33l/HGG6Gv+5BDwj7uvXfDWQUHDw7PffuGVnh9Xnut+tfLUUdBVVUYb2pTTjst/FKB0N+fmfFw3LgwydY++4TjNn169TaZOjVNmwaHHRa2efDBDbepzcyZ4UbKXXcNy8OGhfnWMxr6+XON+vQT8NlnYS7oli3TjkRyUW1zftcmU2f9+vCF0NihuOvbn3uYSXD06NrXt24dnktKQv95Nvurqb7Pe/rpoVtp8OBQt2fPMD3vj34UuqC6dg1dOPEpe7feuvb3ymbe9frijWvo5881auknQDdmyaYcfXRo1VZVheXFi8Nzv34wZkx4/eCDcOih4XXNOcIbeiK0f//QB59JWJn9ZRx0UDhPMGtWWF6xAj74YNPvWXNu9Lhvfat6DvHx48MvmrZtN/1+u+wSkuo111S34DMJvmPHcI6hvhPTGXXNu15XzLvvHlrwmc9///3V860XAiX9BOjGLNmUvfYK85Affnhojf7856H81lvDPNy9eoXEc8st1eVlZaF8zz1DAm+I886Dbt2qTxw/9NCG6zt1Ci3ioUOrT77Wd8LyxBPhiSdqP5F71VXV8V52Weguysbpp8MDD4SuHoDttgtzq++zTzgxvP/+2b1PXfOu15zPPaNNm3DcTz017GuLLbK7GihfaBKVBOy8c2i1PfBA2pGISKHKdhIV9ek3wJw5Db8BBeDTT9W9IyK5QUm/AS64IFy/3Bi77da0sYiINIaSfgN8+mnoF7zttoZt16IF7LRT88QkItIQSvoNUFUVTh5985tpRyIi0ji6eqcBFi/e+EYWEZF8oqSfpRUrwnXCHTqkHYmISOMp6Wcpc+OMkr6I5DMl/Swp6YtIIVDSz1LmVnX16YtIPlPSz5Ja+iJSCLJK+mY2wMxmmtksM7uslvVnm1mFmU2NHufF1g0zsw+jx7CmDD5JSvoiUgjqvU7fzEqAkcCxQDkwyczGunvN2TQfdvcLa2y7PXAlUAo4MDnadkmTRJ8gJX0RKQTZtPQPAGa5+xx3Xw2MAQZl+f7HAc+7++Io0T8PDGhcqOlavDiM150ZS1tEJB9lk/Q7A/Njy+VRWU2nmNk7ZvaomXVtyLZmNtzMysysrKKiIsvQk1VVpZO4IpL/skn6tc1xU3M85n8B3d29F/ACkBkxO5ttcfc73b3U3Us75eicglVV6toRkfySRqixAAAJfklEQVSXTdIvB7rGlrsAC+IV3L3K3aP54rkL6JvttvlCSV9ECkE2SX8S0NPMephZK2AIMDZewcx2jC0OBN6LXj8H9Dez9mbWHugfleUdJX0RKQT1Xr3j7mvN7EJCsi4BRrn7dDMbAZS5+1jgYjMbCKwFFgNnR9suNrNrCF8cACPcffFGO8kDGmxNRApBVkMru/szwDM1yq6Ivf418Os6th0FjNqMGFO3fn1I+mrpi0i+0x25WVi2LCR+JX0RyXdK+lnQjVkiUig0c1bMokVw7bWwevWG5ZlbB5T0RSTfKenH3Hkn3HwzdOy48bru3WHvvRMPSUSkSSnpxzz2GPTrB6+/nnYkIiLNQ336kTlzYOpUGDw47UhERJqPkn7kiSfCs5K+iBSyguneWbwYDjus8duXl8O++0KPHk0Xk4hIrimYpF9SAnvu2fjt99wTzj236eIREclFBZP027WDRx5JOwoRkdymPn0RkSKipC8iUkSU9EVEioiSvohIEVHSFxEpIkr6IiJFRElfRKSIKOmLiBQRc/e0Y9iAmVUAH6cdR6QjUJl2ELXI1bggd2NTXA2juBomF+Layd071Vcp55J+LjGzMncvTTuOmnI1Lsjd2BRXwyiuhsnVuGqj7h0RkSKipC8iUkSU9DftzrQDqEOuxgW5G5viahjF1TC5GtdG1KcvIlJE1NIXESkiRZ30zayrmY0zs/fMbLqZ/SQqv8HM3jezd8zsCTPbLirvbmYrzWxq9Lgj4biuMrNPYvs/PrbNr81slpnNNLPjEo7r4VhMc81salSe1PFqY2YTzeztKK6ro/IeZvammX0YxdgqKm8dLc+K1ndPOK4Ho3+naWY2ysxaRuVHmNmy2PG6IuG47jGzj2L77xOVm5ndGh2vd8xsv4TjejUW0wIzezIqT+R4xeIrMbO3zOzf0XKqf1+N5u5F+wB2BPaLXm8LfADsCfQHWkTl1wPXR6+7A9NSjOsq4JJa6u8JvA20BnoAs4GSpOKqUedG4IqEj5cB20SvWwJvAgcB/wSGROV3ABdEr38E3BG9HgI8nHBcx0frDBgdi+sI4N8pHq97gO/WUv944Nlou4OAN5OMq0adx4DvJXm8Yvv+OfBQZp9p/3019lHULX13/9Tdp0SvPwfeAzq7+/+6+9qo2htAl1yIaxObDALGuPsqd/8ImAUckHRcZmbAaYRElhgPvogWW0YPB44CHo3K7wVOil4PipaJ1h8dxZ5IXO7+TLTOgYkk//dV1/GqyyDgvmi7N4DtzGzHpOMys20J/6ZPNvW+62NmXYDvAHdHy0bKf1+NVdRJPy76CbYvoXURdy6hlZPRI/qJ97KZbcZU7I2O68LoJ/YoM2sflXUG5sc2K2fTXxLNERfAYcBn7v5hrCyR4xX99J4KLAKeJ/zaWRr78o4fk/87XtH6ZUCHJOJy9zdj61oCZwH/iW1ycNS98ayZ7dUcMdUT1/+L/r5uMrPWUVlif1+bOl7AycCL7r48VpbI8QJuBn4JrI+WO5ADf1+NoaQPmNk2hJ+NP43/QZnZ5cBa4MGo6FOgm7vvS/RTz8zaJhjX7cAuQJ8olhszVWvZvNkuy6rreAFD2bCVn9jxcvd17t6H0Go+ANijtmrRc2LHq2ZcZrZ3bPVtwCvu/mq0PIVwK31v4C80Y4u2jrh+DewO7A9sD/wqqp4rx6vm31cix8vMTgAWufvkeHEtVRP/+2qMok/6UWvrMeBBd388Vj4MOAE4I/oZTtR9UhW9nkxoTe6aVFzu/ln0n2I9cBfVXTjlQNfY5l2ABUnFFZW3AAYDD2fKkjxesX0uBcYT+p63i+KCDY/J/x2vaH07YHFCcQ2I9nsl0InwZZipszzTveHuzwAtzaxjUnFF3Xfu7quAf5DC31dtcQGYWYconqdjdZI6XocAA81sLjCG0K1zMzn099UQRZ30o362vwPvufufY+UDCK2cge6+IlbeycxKotc7Az2BOQnGFe9HPRmYFr0eCwyJrhroEcU1Mam4IscA77t7eax+Userk1VfYbVlFMt7wDjgu1G1YcBT0eux0TLR+pcyX+wJxPW+mZ0HHAcMjb7AM/W/nun7NbMDCP8/qxKMa8eozAj90/G/r+9ZcBCwzN0/TSquaPWphBOoX8XqJ3K83P3X7t7F3bsTTsy+5O5nkPLfV6PVd6a3kB/AoYSfXe8AU6PH8YQTofNjZZkz8acA0wlXykwBTkw4rvuBd6PyscCOsW0uJ7SkZwLfTjKuaN09wPk16id1vHoBb0VxTaP66qGdCV9+s4BHgNZReZtoeVa0fueE41ob/VtljmGm/MLY8XoD6JdwXC9Ff1/TgAeovpLGgJFRzO8CpUnGFa0bT/g1Eq+fyPGqsc8jqL56J9W/r8Y+dEeuiEgRKeruHRGRYqOkLyJSRJT0RUSKiJK+iEgRUdIXESkiSvpSNMzsp2a2VSO2O9vMvtGA+gPN7LKG7kckCbpkU4pGdEdlqbtXNmCbEuBFwuimZc0Vm0hS1NKXgmNmW5vZ09FAXNPM7HQzuxj4BjDOzMZF9W43szKLjd0elc81syvM7DXCeC+lwIMWxmzfssa+LjazGdEgZWOisrPN7K/R66mxx0ozOzyKb5SZTbIwGN2ghA6NCC3qryKSdwYAC9z9OwBm1s7dl5nZz4EjYy39y919caY1b2a93P2daN1X7n5otP151N3Svwzo4e6rMkMIxHkYPAwzO5EwSuME4GrCrfnnRttMNLMX3P3LJjsCInVQS18K0bvAMWZ2vZkd5u7L6qh3mplNIdz6vxdhMpqMh2vfZCPvEH4FnEkYXmEjZtYTuAE43d3XECbpuSwaQng84bb9blnuT2SzKOlLwXH3D4C+hOR/rdUyjV40MN0lwNHu3oswemObWJVsW93fIYxL0xeYHBt1MbOfrQkzLP3A3TOjMBpwirv3iR7d3P297D+hSOMp6UvBia60WeHuDwB/AjJzun5OmOYRoC0hsS8zs68B397EW8a3i+9nC6Cru48jdN1sB2xTo9o/gH949Zj5AM8BF8VGiNy3AR9PZLOoT18K0T7ADWa2HlgDXBCV3wk8a2afuvuRZvYWYZTGOcDrm3i/e4A7zGwlcLC7r4zKS4AHzKwdofV+k7svjXI5ZrYTYWjdXc3s3Gib84BrCOOxvxMl/rmEuRtEmp0u2RQRKSLq3hERKSJK+iIiRURJX0SkiCjpi4gUESV9EZEioqQvIlJElPRFRIqIkr6ISBH5/0qvNVb62T1jAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(x, y, 'r')\n", "plt.plot(x, z, 'b')\n", "plt.xlabel('start size')\n", "plt.text(400, 0.55, '% taxa covered', color='r', ha='right')\n", "plt.text(400, 0.52, 'coefficient of variation', color='b', ha='right');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final starting size *k* is manually determined based on the plot." ] }, { "cell_type": "code", "execution_count": 135, "metadata": {}, "outputs": [], "source": [ "start_size = 348" ] }, { "cell_type": "code", "execution_count": 136, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "N1081,N646: 350 taxa.\n", "N3758,N8774: 346 taxa.\n", "N2668,N2313: 353 taxa.\n", "N191,N458: 358 taxa.\n", "N1640,N2157: 336 taxa.\n", "N2666,N1610: 367 taxa.\n", "N146,N162: 323 taxa.\n", "N5329,N3630: 315 taxa.\n", "N6294,N6394: 315 taxa.\n", "N4527,N1716: 300 taxa.\n", "N4526,N2077: 291 taxa.\n", "N1838,N2733: 252 taxa.\n", "N4162,N2076: 244 taxa.\n", "N4999,N4055: 199 taxa.\n", "N207,N90: 197 taxa.\n", "N3529,N4053: 197 taxa.\n", "N4218,N3629: 184 taxa.\n", "N250,N428: 156 taxa.\n", "N604,N531: 149 taxa.\n", "N4008,N1609: 131 taxa.\n", "N3029,N1363: 129 taxa.\n", "N3855,N2903: 116 taxa.\n", "N4256,N4934: 109 taxa.\n", "N2699,N3445: 108 taxa.\n", "N4628,N4938: 108 taxa.\n", "N233,N1442: 107 taxa.\n", "N1604,N1435: 106 taxa.\n", "N1636,N2153: 106 taxa.\n", "N117,N21: 105 taxa.\n", "N3757,N8773: 97 taxa.\n", "N2422,N3819: 86 taxa.\n", "N375,N117: 85 taxa.\n", "N1832,N2000: 80 taxa.\n", "N752,N681: 73 taxa.\n", "N4244,N4051: 70 taxa.\n", "N1306,N63: 69 taxa.\n", "N930,N680: 68 taxa.\n", "N1287,N825: 66 taxa.\n", "N1603,N1159: 65 taxa.\n", "N54,N315: 64 taxa.\n", "N923,N354: 64 taxa.\n", "N1602,N522: 63 taxa.\n", "N33,N23: 61 taxa.\n", "N2809,N2483: 60 taxa.\n", "N891,N642: 59 taxa.\n", "N3656,N2312: 57 taxa.\n", "N2116,N2727: 56 taxa.\n", "N190,N457: 55 taxa.\n", "N4616,N4464: 55 taxa.\n", "N1304,N152: 54 taxa.\n" ] } ], "source": [ "selected_matches = select_clades(start_size)\n", "for m in selected_matches:\n", " print('%s: %d taxa.' % (m, len(matches[m])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write selected clades." ] }, { "cell_type": "code", "execution_count": 137, "metadata": {}, "outputs": [], "source": [ "with open('selected_matches.txt', 'w') as f:\n", " for m in selected_matches:\n", " f.write('%s\\n' % m.replace(',', '\\t'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Export collapsed trees. Each selected clade becomes a tip, with its height equaling the median of its descendants." ] }, { "cell_type": "code", "execution_count": 138, "metadata": {}, "outputs": [], "source": [ "def get_clade_height(node):\n", " return np.median([x.accumulate_to_ancestor(node) for x in node.tips()])" ] }, { "cell_type": "code", "execution_count": 139, "metadata": {}, "outputs": [], "source": [ "def collapse_clades(tree, clades):\n", " tcopy = tree.copy()\n", " clades = set(clades)\n", " nodes_to_remove = []\n", " for node in tcopy.non_tips():\n", " if node.name in clades:\n", " node.length += get_clade_height(node)\n", " nodes_to_remove.extend(node.children)\n", " tcopy.remove_deleted(lambda x: x in nodes_to_remove)\n", " tcopy = tcopy.shear(clades)\n", " tcopy.prune()\n", " return tcopy" ] }, { "cell_type": "code", "execution_count": 140, "metadata": {}, "outputs": [], "source": [ "for i, tree in enumerate(trees):\n", " clades = [x.split(',')[i] for x in selected_matches]\n", " tcopy = collapse_clades(tree, clades)\n", " tcopy.write('tree%d.nwk' % (i + 1))\n", " for tip in tcopy.tips():\n", " tip.name = 'X%d' % (clades.index(tip.name) + 1)\n", " tcopy.write('tree%d.trans.nwk' % (i + 1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "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.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }