testdocs/analysis/cookbook-tetrad-ipcoal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><span style=\"color:gray\">ipyrad-analysis toolkit:</span> tetrad</h1>\n",
    "\n",
    "The `tetrad` tool is a framework for inferring a species tree topology using quartet-joining on large unlinked SNP data sets. It is particularly optimized for RAD-seq type datasets that are likely to involve a lot of missing data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# conda install ipyrad -c conda-forge -c bioconda\n",
    "# conda install tetrad -c conda-forge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipyrad.analysis as ipa\n",
    "import toytree\n",
    "import ipcoal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate a random tree with 20 tips and crown age of 10M generations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div class=\"toyplot\" id=\"t2e81144ec985480f8f7d032e4c5694a1\" style=\"text-align:center\"><svg class=\"toyplot-canvas-Canvas\" height=\"360.0px\" id=\"ta278a41b6e5d46d88a8df63527b63a4d\" preserveAspectRatio=\"xMidYMid meet\" style=\"background-color:transparent;border-color:#292724;border-style:none;border-width:1.0;fill:rgb(16.1%,15.3%,14.1%);fill-opacity:1.0;font-family:Helvetica;font-size:12px;opacity:1.0;stroke:rgb(16.1%,15.3%,14.1%);stroke-opacity:1.0;stroke-width:1.0\" viewBox=\"0 0 265.0 360.0\" width=\"265.0px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:toyplot=\"http://www.sandia.gov/toyplot\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><g class=\"toyplot-coordinates-Cartesian\" id=\"t21d9ee58a9b7444bb94e8001de5b12cc\"><clipPath id=\"ta0f2915eb36049fd9cb26102375c951d\"><rect height=\"300.0\" width=\"205.0\" x=\"30.0\" y=\"30.0\"></rect></clipPath><g clip-path=\"url(#ta0f2915eb36049fd9cb26102375c951d)\"><g class=\"toytree-mark-Toytree\" id=\"tfe8d84f7034b4f2698e5cce84a80798c\"><g class=\"toytree-Edges\" style=\"fill:none;stroke:rgb(14.9%,14.9%,14.9%);stroke-linecap:round;stroke-opacity:1;stroke-width:2\"><path d=\"M 51.7 124.6 L 51.7 70.8 L 88.6 70.8\" id=\"38,37\"></path><path d=\"M 51.7 124.6 L 51.7 178.4 L 130.7 178.4\" id=\"38,36\"></path><path d=\"M 88.6 70.8 L 88.6 80.2 L 155.4 80.2\" id=\"37,35\"></path><path d=\"M 130.7 178.4 L 130.7 127.8 L 152.3 127.8\" id=\"36,34\"></path><path d=\"M 130.7 178.4 L 130.7 229.1 L 164.8 229.1\" id=\"36,33\"></path><path d=\"M 152.3 127.8 L 152.3 108.3 L 157.9 108.3\" id=\"34,32\"></path><path d=\"M 152.3 127.8 L 152.3 147.3 L 158.7 147.3\" id=\"34,31\"></path><path d=\"M 164.8 229.1 L 164.8 197.2 L 167.1 197.2\" id=\"33,30\"></path><path d=\"M 164.8 229.1 L 164.8 261.1 L 169.4 261.1\" id=\"33,29\"></path><path d=\"M 157.9 108.3 L 157.9 117.6 L 181.8 117.6\" id=\"32,28\"></path><path d=\"M 158.7 147.3 L 158.7 158.2 L 184.8 158.2\" id=\"31,27\"></path><path d=\"M 167.1 197.2 L 167.1 208.1 L 178.3 208.1\" id=\"30,26\"></path><path d=\"M 169.4 261.1 L 169.4 242.4 L 186.5 242.4\" id=\"29,25\"></path><path d=\"M 169.4 261.1 L 169.4 279.8 L 172.4 279.8\" id=\"29,24\"></path><path d=\"M 184.8 158.2 L 184.8 167.5 L 188.1 167.5\" id=\"27,23\"></path><path d=\"M 178.3 208.1 L 178.3 217.4 L 179.9 217.4\" id=\"26,22\"></path><path d=\"M 172.4 279.8 L 172.4 267.3 L 179.4 267.3\" id=\"24,21\"></path><path d=\"M 172.4 279.8 L 172.4 292.3 L 176.9 292.3\" id=\"24,20\"></path><path d=\"M 88.6 70.8 L 88.6 61.5 L 189.6 61.5\" id=\"37,19\"></path><path d=\"M 155.4 80.2 L 155.4 74.0 L 189.6 74.0\" id=\"35,18\"></path><path d=\"M 155.4 80.2 L 155.4 86.4 L 189.6 86.4\" id=\"35,17\"></path><path d=\"M 157.9 108.3 L 157.9 98.9 L 189.6 98.9\" id=\"32,16\"></path><path d=\"M 181.8 117.6 L 181.8 111.4 L 189.6 111.4\" id=\"28,15\"></path><path d=\"M 181.8 117.6 L 181.8 123.9 L 189.6 123.9\" id=\"28,14\"></path><path d=\"M 158.7 147.3 L 158.7 136.3 L 189.6 136.3\" id=\"31,13\"></path><path d=\"M 184.8 158.2 L 184.8 148.8 L 189.6 148.8\" id=\"27,12\"></path><path d=\"M 188.1 167.5 L 188.1 161.3 L 189.6 161.3\" id=\"23,11\"></path><path d=\"M 188.1 167.5 L 188.1 173.8 L 189.6 173.8\" id=\"23,10\"></path><path d=\"M 167.1 197.2 L 167.1 186.2 L 189.6 186.2\" id=\"30,9\"></path><path d=\"M 178.3 208.1 L 178.3 198.7 L 189.6 198.7\" id=\"26,8\"></path><path d=\"M 179.9 217.4 L 179.9 211.2 L 189.6 211.2\" id=\"22,7\"></path><path d=\"M 179.9 217.4 L 179.9 223.7 L 189.6 223.7\" id=\"22,6\"></path><path d=\"M 186.5 242.4 L 186.5 236.1 L 189.6 236.1\" id=\"25,5\"></path><path d=\"M 186.5 242.4 L 186.5 248.6 L 189.6 248.6\" id=\"25,4\"></path><path d=\"M 179.4 267.3 L 179.4 261.1 L 189.6 261.1\" id=\"21,3\"></path><path d=\"M 179.4 267.3 L 179.4 273.6 L 189.6 273.6\" id=\"21,2\"></path><path d=\"M 176.9 292.3 L 176.9 286.0 L 189.6 286.0\" id=\"20,1\"></path><path d=\"M 176.9 292.3 L 176.9 298.5 L 189.6 298.5\" id=\"20,0\"></path></g><g class=\"toytree-TipLabels\" style=\"fill:rgb(14.9%,14.9%,14.9%);fill-opacity:1.0;font-family:helvetica;font-size:11px;font-weight:normal;stroke:none;white-space:pre\"><g transform=\"translate(189.63,298.51)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r0</text></g><g transform=\"translate(189.63,286.04)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r1</text></g><g transform=\"translate(189.63,273.56)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r2</text></g><g transform=\"translate(189.63,261.09)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r3</text></g><g transform=\"translate(189.63,248.61)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r4</text></g><g transform=\"translate(189.63,236.14)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r5</text></g><g transform=\"translate(189.63,223.66)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r6</text></g><g transform=\"translate(189.63,211.19)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r7</text></g><g transform=\"translate(189.63,198.71)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r8</text></g><g transform=\"translate(189.63,186.24)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r9</text></g><g transform=\"translate(189.63,173.76)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r10</text></g><g transform=\"translate(189.63,161.29)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r11</text></g><g transform=\"translate(189.63,148.81)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r12</text></g><g transform=\"translate(189.63,136.34)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r13</text></g><g transform=\"translate(189.63,123.86)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r14</text></g><g transform=\"translate(189.63,111.39)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r15</text></g><g transform=\"translate(189.63,98.91)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r16</text></g><g transform=\"translate(189.63,86.44)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r17</text></g><g transform=\"translate(189.63,73.96)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r18</text></g><g transform=\"translate(189.63,61.49)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r19</text></g></g></g></g><g class=\"toyplot-coordinates-Axis\" id=\"tf397413142f940a7916991f56370bc63\" transform=\"translate(50.0,310.0)translate(0,20.0)\"><line style=\"\" x1=\"1.672198434214192\" x2=\"139.62856925688516\" y1=\"0\" y2=\"0\"></line><g><line style=\"\" x1=\"139.62856925688516\" x2=\"139.62856925688516\" y1=\"0\" y2=\"-5\"></line><line style=\"\" x1=\"70.65038384554968\" x2=\"70.65038384554968\" y1=\"0\" y2=\"-5\"></line><line style=\"\" x1=\"1.672198434214192\" x2=\"1.672198434214192\" y1=\"0\" y2=\"-5\"></line></g><g><g transform=\"translate(139.62856925688516,6)\"><text style=\"fill:rgb(16.1%,15.3%,14.1%);fill-opacity:1.0;font-family:helvetica;font-size:10.0px;font-weight:normal;stroke:none;vertical-align:baseline;white-space:pre\" x=\"-2.78\" y=\"8.555\">0</text></g><g transform=\"translate(70.65038384554968,6)\"><text style=\"fill:rgb(16.1%,15.3%,14.1%);fill-opacity:1.0;font-family:helvetica;font-size:10.0px;font-weight:normal;stroke:none;vertical-align:baseline;white-space:pre\" x=\"-19.459999999999997\" y=\"8.555\">5000000</text></g><g transform=\"translate(1.672198434214192,6)\"><text style=\"fill:rgb(16.1%,15.3%,14.1%);fill-opacity:1.0;font-family:helvetica;font-size:10.0px;font-weight:normal;stroke:none;vertical-align:baseline;white-space:pre\" x=\"-22.24\" y=\"8.555\">10000000</text></g></g><g class=\"toyplot-coordinates-Axis-coordinates\" style=\"visibility:hidden\" transform=\"\"><line style=\"stroke:rgb(43.9%,50.2%,56.5%);stroke-opacity:1.0;stroke-width:1.0\" x1=\"0\" x2=\"0\" y1=\"-3.0\" y2=\"4.5\"></line><text style=\"alignment-baseline:alphabetic;fill:rgb(43.9%,50.2%,56.5%);fill-opacity:1.0;font-size:10px;font-weight:normal;stroke:none;text-anchor:middle\" x=\"0\" y=\"-6\"></text></g></g></g></svg><div class=\"toyplot-behavior\"><script>(function()\n",
       "{\n",
       "var modules={};\n",
       "modules[\"toyplot/canvas/id\"] = \"ta278a41b6e5d46d88a8df63527b63a4d\";\n",
       "modules[\"toyplot/canvas\"] = (function(canvas_id)\n",
       "    {\n",
       "        return document.querySelector(\"#\" + canvas_id);\n",
       "    })(modules[\"toyplot/canvas/id\"]);\n",
       "modules[\"toyplot.coordinates.Axis\"] = (\n",
       "        function(canvas)\n",
       "        {\n",
       "            function sign(x)\n",
       "            {\n",
       "                return x < 0 ? -1 : x > 0 ? 1 : 0;\n",
       "            }\n",
       "\n",
       "            function mix(a, b, amount)\n",
       "            {\n",
       "                return ((1.0 - amount) * a) + (amount * b);\n",
       "            }\n",
       "\n",
       "            function log(x, base)\n",
       "            {\n",
       "                return Math.log(Math.abs(x)) / Math.log(base);\n",
       "            }\n",
       "\n",
       "            function in_range(a, x, b)\n",
       "            {\n",
       "                var left = Math.min(a, b);\n",
       "                var right = Math.max(a, b);\n",
       "                return left <= x && x <= right;\n",
       "            }\n",
       "\n",
       "            function inside(range, projection)\n",
       "            {\n",
       "                for(var i = 0; i != projection.length; ++i)\n",
       "                {\n",
       "                    var segment = projection[i];\n",
       "                    if(in_range(segment.range.min, range, segment.range.max))\n",
       "                        return true;\n",
       "                }\n",
       "                return false;\n",
       "            }\n",
       "\n",
       "            function to_domain(range, projection)\n",
       "            {\n",
       "                for(var i = 0; i != projection.length; ++i)\n",
       "                {\n",
       "                    var segment = projection[i];\n",
       "                    if(in_range(segment.range.bounds.min, range, segment.range.bounds.max))\n",
       "                    {\n",
       "                        if(segment.scale == \"linear\")\n",
       "                        {\n",
       "                            var amount = (range - segment.range.min) / (segment.range.max - segment.range.min);\n",
       "                            return mix(segment.domain.min, segment.domain.max, amount)\n",
       "                        }\n",
       "                        else if(segment.scale[0] == \"log\")\n",
       "                        {\n",
       "                            var amount = (range - segment.range.min) / (segment.range.max - segment.range.min);\n",
       "                            var base = segment.scale[1];\n",
       "                            return sign(segment.domain.min) * Math.pow(base, mix(log(segment.domain.min, base), log(segment.domain.max, base), amount));\n",
       "                        }\n",
       "                    }\n",
       "                }\n",
       "            }\n",
       "\n",
       "            var axes = {};\n",
       "\n",
       "            function display_coordinates(e)\n",
       "            {\n",
       "                var current = canvas.createSVGPoint();\n",
       "                current.x = e.clientX;\n",
       "                current.y = e.clientY;\n",
       "\n",
       "                for(var axis_id in axes)\n",
       "                {\n",
       "                    var axis = document.querySelector(\"#\" + axis_id);\n",
       "                    var coordinates = axis.querySelector(\".toyplot-coordinates-Axis-coordinates\");\n",
       "                    if(coordinates)\n",
       "                    {\n",
       "                        var projection = axes[axis_id];\n",
       "                        var local = current.matrixTransform(axis.getScreenCTM().inverse());\n",
       "                        if(inside(local.x, projection))\n",
       "                        {\n",
       "                            var domain = to_domain(local.x, projection);\n",
       "                            coordinates.style.visibility = \"visible\";\n",
       "                            coordinates.setAttribute(\"transform\", \"translate(\" + local.x + \")\");\n",
       "                            var text = coordinates.querySelector(\"text\");\n",
       "                            text.textContent = domain.toFixed(2);\n",
       "                        }\n",
       "                        else\n",
       "                        {\n",
       "                            coordinates.style.visibility= \"hidden\";\n",
       "                        }\n",
       "                    }\n",
       "                }\n",
       "            }\n",
       "\n",
       "            canvas.addEventListener(\"click\", display_coordinates);\n",
       "\n",
       "            var module = {};\n",
       "            module.show_coordinates = function(axis_id, projection)\n",
       "            {\n",
       "                axes[axis_id] = projection;\n",
       "            }\n",
       "\n",
       "            return module;\n",
       "        })(modules[\"toyplot/canvas\"]);\n",
       "(function(axis, axis_id, projection)\n",
       "        {\n",
       "            axis.show_coordinates(axis_id, projection);\n",
       "        })(modules[\"toyplot.coordinates.Axis\"],\"tf397413142f940a7916991f56370bc63\",[{\"domain\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 1839090.909090908, \"min\": -10121212.121212121}, \"range\": {\"bounds\": {\"max\": Infinity, \"min\": -Infinity}, \"max\": 165.0, \"min\": 0.0}, \"scale\": \"linear\"}]);\n",
       "})();</script></div></div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tree = toytree.rtree.bdtree(ntips=20, seed=555)\n",
    "tree = tree.mod.node_scale_root_height(10e6) \n",
    "tree.draw(scalebar=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate SNPs with missing data and write to database (.seqs.hdf5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wrote 259208 SNPs to /tmp/test-tet-miss50.snps.hdf5\n"
     ]
    }
   ],
   "source": [
    "# init simulator with one diploid sample from each tip\n",
    "model = ipcoal.Model(tree, Ne=1e6, nsamples=2, recomb=0)\n",
    "\n",
    "# simulate sequence data on 10K loci\n",
    "model.sim_loci(10000, 50)\n",
    "\n",
    "# add missing data (50%)\n",
    "model.apply_missing_mask(0.5)\n",
    "\n",
    "# write results to database file\n",
    "model.write_snps_to_hdf5(name=\"test-tet-miss50\", outdir='/tmp', diploid=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Infer tetrad tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNPS = \"/tmp/test-tet-miss50.snps.hdf5\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading snps array [20 taxa x 259208 snps]\n",
      "max unlinked SNPs per quartet [nloci]: 10000\n",
      "quartet sampler [full]: 4845 / 4845\n"
     ]
    }
   ],
   "source": [
    "tet = ipa.tetrad(\n",
    "    data=SNPS,\n",
    "    name=\"test-tet-miss50\",\n",
    "    workdir=\"/tmp\",\n",
    "    nboots=10, \n",
    "    nquartets=1e6,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parallel connection | latituba: 8 cores\n",
      "initializing quartet sets database\n",
      "[####################] 100% 0:00:14 | full tree * | mean SNPs/qrt: 3167 \n",
      "[####################] 100% 0:00:12 | boot rep. 1 | mean SNPs/qrt: 3163 \n",
      "[####################] 100% 0:00:12 | boot rep. 2 | mean SNPs/qrt: 3161 \n",
      "[####################] 100% 0:00:12 | boot rep. 3 | mean SNPs/qrt: 3170 \n",
      "[####################] 100% 0:00:12 | boot rep. 4 | mean SNPs/qrt: 3144 \n",
      "[####################] 100% 0:00:11 | boot rep. 5 | mean SNPs/qrt: 3172 \n",
      "[####################] 100% 0:00:11 | boot rep. 6 | mean SNPs/qrt: 3191 \n",
      "[####################] 100% 0:00:11 | boot rep. 7 | mean SNPs/qrt: 3144 \n",
      "[####################] 100% 0:00:11 | boot rep. 8 | mean SNPs/qrt: 3153 \n",
      "[####################] 100% 0:00:11 | boot rep. 9 | mean SNPs/qrt: 3180 \n",
      "[####################] 100% 0:00:11 | boot rep. 10 | mean SNPs/qrt: 3128 \n"
     ]
    }
   ],
   "source": [
    "tet.run(auto=True, force=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Draw the inferred tetrad tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div class=\"toyplot\" id=\"t8e1e426b34f24ebbbf53a24c1eec84ee\" style=\"text-align:center\"><svg class=\"toyplot-canvas-Canvas\" height=\"360.0px\" id=\"t2117e0dac77847e6b57f17bdfd096e54\" preserveAspectRatio=\"xMidYMid meet\" style=\"background-color:transparent;border-color:#292724;border-style:none;border-width:1.0;fill:rgb(16.1%,15.3%,14.1%);fill-opacity:1.0;font-family:Helvetica;font-size:12px;opacity:1.0;stroke:rgb(16.1%,15.3%,14.1%);stroke-opacity:1.0;stroke-width:1.0\" viewBox=\"0 0 265.0 360.0\" width=\"265.0px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:toyplot=\"http://www.sandia.gov/toyplot\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><g class=\"toyplot-coordinates-Cartesian\" id=\"tcb3e257b1bf546dba1bdeb797a3de19a\"><clipPath id=\"t3235cd5270e04079b8040a6e1f9ad82a\"><rect height=\"300.0\" width=\"205.0\" x=\"30.0\" y=\"30.0\"></rect></clipPath><g clip-path=\"url(#t3235cd5270e04079b8040a6e1f9ad82a)\"><g class=\"toytree-mark-Toytree\" id=\"ted1fa960fe414d4dbd32641bea39e1e0\"><g class=\"toytree-Edges\" style=\"fill:none;stroke:rgb(40%,76.1%,64.7%);stroke-linecap:round;stroke-opacity:1;stroke-width:2\"><path d=\"M 51.7 124.6 L 51.7 70.8 L 143.6 70.8\" id=\"38,37\"></path><path d=\"M 51.7 124.6 L 51.7 178.4 L 74.7 178.4\" id=\"38,36\"></path><path d=\"M 143.6 70.8 L 143.6 80.2 L 166.6 80.2\" id=\"37,35\"></path><path d=\"M 74.7 178.4 L 74.7 127.8 L 97.7 127.8\" id=\"36,34\"></path><path d=\"M 74.7 178.4 L 74.7 229.1 L 97.7 229.1\" id=\"36,33\"></path><path d=\"M 97.7 127.8 L 97.7 108.3 L 143.6 108.3\" id=\"34,32\"></path><path d=\"M 97.7 127.8 L 97.7 147.3 L 120.7 147.3\" id=\"34,31\"></path><path d=\"M 97.7 229.1 L 97.7 197.2 L 120.7 197.2\" id=\"33,30\"></path><path d=\"M 97.7 229.1 L 97.7 261.1 L 120.7 261.1\" id=\"33,29\"></path><path d=\"M 143.6 108.3 L 143.6 117.6 L 166.6 117.6\" id=\"32,28\"></path><path d=\"M 120.7 147.3 L 120.7 158.2 L 143.6 158.2\" id=\"31,27\"></path><path d=\"M 120.7 197.2 L 120.7 208.1 L 143.6 208.1\" id=\"30,26\"></path><path d=\"M 120.7 261.1 L 120.7 242.4 L 166.6 242.4\" id=\"29,25\"></path><path d=\"M 120.7 261.1 L 120.7 279.8 L 143.6 279.8\" id=\"29,24\"></path><path d=\"M 143.6 158.2 L 143.6 167.5 L 166.6 167.5\" id=\"27,23\"></path><path d=\"M 143.6 208.1 L 143.6 217.4 L 166.6 217.4\" id=\"26,22\"></path><path d=\"M 143.6 279.8 L 143.6 267.3 L 166.6 267.3\" id=\"24,21\"></path><path d=\"M 143.6 279.8 L 143.6 292.3 L 166.6 292.3\" id=\"24,20\"></path><path d=\"M 143.6 70.8 L 143.6 61.5 L 189.6 61.5\" id=\"37,19\"></path><path d=\"M 166.6 80.2 L 166.6 74.0 L 189.6 74.0\" id=\"35,18\"></path><path d=\"M 166.6 80.2 L 166.6 86.4 L 189.6 86.4\" id=\"35,17\"></path><path d=\"M 143.6 108.3 L 143.6 98.9 L 189.6 98.9\" id=\"32,16\"></path><path d=\"M 166.6 117.6 L 166.6 111.4 L 189.6 111.4\" id=\"28,15\"></path><path d=\"M 166.6 117.6 L 166.6 123.9 L 189.6 123.9\" id=\"28,14\"></path><path d=\"M 120.7 147.3 L 120.7 136.3 L 189.6 136.3\" id=\"31,13\"></path><path d=\"M 143.6 158.2 L 143.6 148.8 L 189.6 148.8\" id=\"27,12\"></path><path d=\"M 166.6 167.5 L 166.6 161.3 L 189.6 161.3\" id=\"23,11\"></path><path d=\"M 166.6 167.5 L 166.6 173.8 L 189.6 173.8\" id=\"23,10\"></path><path d=\"M 120.7 197.2 L 120.7 186.2 L 189.6 186.2\" id=\"30,9\"></path><path d=\"M 143.6 208.1 L 143.6 198.7 L 189.6 198.7\" id=\"26,8\"></path><path d=\"M 166.6 217.4 L 166.6 211.2 L 189.6 211.2\" id=\"22,7\"></path><path d=\"M 166.6 217.4 L 166.6 223.7 L 189.6 223.7\" id=\"22,6\"></path><path d=\"M 166.6 242.4 L 166.6 236.1 L 189.6 236.1\" id=\"25,5\"></path><path d=\"M 166.6 242.4 L 166.6 248.6 L 189.6 248.6\" id=\"25,4\"></path><path d=\"M 166.6 267.3 L 166.6 261.1 L 189.6 261.1\" id=\"21,3\"></path><path d=\"M 166.6 267.3 L 166.6 273.6 L 189.6 273.6\" id=\"21,2\"></path><path d=\"M 166.6 292.3 L 166.6 286.0 L 189.6 286.0\" id=\"20,1\"></path><path d=\"M 166.6 292.3 L 166.6 298.5 L 189.6 298.5\" id=\"20,0\"></path></g><g class=\"toytree-AlignEdges\" style=\"stroke:rgb(66.3%,66.3%,66.3%);stroke-dasharray:2, 4;stroke-linecap:round;stroke-opacity:1.0;stroke-width:2\"><path d=\"M 189.6 298.5 L 189.6 298.5\"></path><path d=\"M 189.6 286.0 L 189.6 286.0\"></path><path d=\"M 189.6 273.6 L 189.6 273.6\"></path><path d=\"M 189.6 261.1 L 189.6 261.1\"></path><path d=\"M 189.6 248.6 L 189.6 248.6\"></path><path d=\"M 189.6 236.1 L 189.6 236.1\"></path><path d=\"M 189.6 223.7 L 189.6 223.7\"></path><path d=\"M 189.6 211.2 L 189.6 211.2\"></path><path d=\"M 189.6 198.7 L 189.6 198.7\"></path><path d=\"M 189.6 186.2 L 189.6 186.2\"></path><path d=\"M 189.6 173.8 L 189.6 173.8\"></path><path d=\"M 189.6 161.3 L 189.6 161.3\"></path><path d=\"M 189.6 148.8 L 189.6 148.8\"></path><path d=\"M 189.6 136.3 L 189.6 136.3\"></path><path d=\"M 189.6 123.9 L 189.6 123.9\"></path><path d=\"M 189.6 111.4 L 189.6 111.4\"></path><path d=\"M 189.6 98.9 L 189.6 98.9\"></path><path d=\"M 189.6 86.4 L 189.6 86.4\"></path><path d=\"M 189.6 74.0 L 189.6 74.0\"></path><path d=\"M 189.6 61.5 L 189.6 61.5\"></path></g><g class=\"toytree-NodeLabels\" style=\"fill:rgb(14.9%,14.9%,14.9%);fill-opacity:1.0;font-size:9px;stroke:none\"><g transform=\"translate(159.13,294.58)\"><text>100</text></g><g transform=\"translate(159.13,269.63)\"><text>100</text></g><g transform=\"translate(161.63,219.72)\"><text>70</text></g><g transform=\"translate(159.13,169.82)\"><text>100</text></g><g transform=\"translate(136.14,282.10)\"><text>100</text></g><g transform=\"translate(159.13,244.67)\"><text>100</text></g><g transform=\"translate(136.14,210.37)\"><text>100</text></g><g transform=\"translate(136.14,160.47)\"><text>100</text></g><g transform=\"translate(159.13,119.92)\"><text>100</text></g><g transform=\"translate(113.14,263.39)\"><text>100</text></g><g transform=\"translate(113.14,199.45)\"><text>100</text></g><g transform=\"translate(113.14,149.55)\"><text>100</text></g><g transform=\"translate(136.14,110.57)\"><text>100</text></g><g transform=\"translate(90.15,231.42)\"><text>100</text></g><g transform=\"translate(90.15,130.06)\"><text>100</text></g><g transform=\"translate(159.13,82.50)\"><text>100</text></g><g transform=\"translate(67.16,180.74)\"><text>100</text></g><g transform=\"translate(136.14,73.14)\"><text>100</text></g></g><g class=\"toytree-TipLabels\" style=\"fill:rgb(90.6%,54.1%,76.5%);fill-opacity:1.0;font-family:helvetica;font-size:11px;font-weight:normal;stroke:none;white-space:pre\"><g transform=\"translate(189.63,298.51)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r3</text></g><g transform=\"translate(189.63,286.04)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r2</text></g><g transform=\"translate(189.63,273.56)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r0</text></g><g transform=\"translate(189.63,261.09)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r1</text></g><g transform=\"translate(189.63,248.61)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r4</text></g><g transform=\"translate(189.63,236.14)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r5</text></g><g transform=\"translate(189.63,223.66)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r7</text></g><g transform=\"translate(189.63,211.19)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r6</text></g><g transform=\"translate(189.63,198.71)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r8</text></g><g transform=\"translate(189.63,186.24)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r9</text></g><g transform=\"translate(189.63,173.76)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r11</text></g><g transform=\"translate(189.63,161.29)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r10</text></g><g transform=\"translate(189.63,148.81)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r12</text></g><g transform=\"translate(189.63,136.34)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r13</text></g><g transform=\"translate(189.63,123.86)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r15</text></g><g transform=\"translate(189.63,111.39)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r14</text></g><g transform=\"translate(189.63,98.91)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r16</text></g><g transform=\"translate(189.63,86.44)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r18</text></g><g transform=\"translate(189.63,73.96)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r17</text></g><g transform=\"translate(189.63,61.49)rotate(0)\"><text style=\"\" x=\"15.00\" y=\"2.81\">r19</text></g></g></g></g></g></svg><div class=\"toyplot-behavior\"><script>(function()\n",
       "{\n",
       "var modules={};\n",
       "})();</script></div></div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tre = toytree.tree(tet.trees.cons)\n",
    "rtre = tre.root([\"r19\", \"r18\", \"r17\"])\n",
    "rtre.draw(ts='d', use_edge_lengths=False, node_labels=\"support\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Does this tree match the *true* species tree?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rfdist = rtre.treenode.robinson_foulds(tree.treenode)[0]\n",
    "rfdist == 0"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}