ts_workshop_part1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "vO9IP0P3aeir"
   },
   "source": [
    "# Tree sequence workshop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_JrVmHtqaeit"
   },
   "source": [
    "### Friday 14 December"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5FHoBkSsaeit"
   },
   "source": [
    "Welcome to the tree sequence workshop. Today we will give you an overview of the tree sequence data structure, the $\\texttt{msprime}$ coalescent simulator, and our scalable inference framework: $\\texttt{tsinfer}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ixMoi6pcaeiu"
   },
   "source": [
    "# Outline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xMWuwAnDaeiv"
   },
   "source": [
    "\n",
    "1. [The tree sequence format](#ts)\n",
    "2. [<span style=\"color: gray\">Simulations (msprime & slim)</span>](ts_workshop_part2.ipynb#msprime)\n",
    "3. [<span style=\"color: gray\">Calculating statistics with tree sequences</span>](ts_workshop_part2.ipynb#stats)\n",
    "4. [<span style=\"color: gray\">Inferring tree sequences (tsinfer)</span>](ts_workshop_part3.ipynb#tsinfer)\n",
    "5. [<span style=\"color: gray\">Metadata in tree sequences</span>](ts_workshop_part3.ipynb#metadata)\n",
    "6. [<span style=\"color: gray\">References and Documentation</span>](ts_workshop_part3.ipynb#ref)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "IIt2MpE3aeix",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from IPython.display import SVG, display, HTML\n",
    "import msprime as tskit # NB allows use of future terminology e.g. tskit.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>aside {color:grey;}</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<style>aside {color:grey;}</style>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NOtPQ3W3aei1"
   },
   "source": [
    "<a id='ts'></a>\n",
    "# 1. The tree sequence format\n",
    "\n",
    "A succint tree sequence is an efficient way of storing the ancestry of different parts of the genome. It can be used to store the output of evolutionary simulations, or to store an ancestry inferred from real data.\n",
    "\n",
    "The primary benefits to this are:\n",
    "1. Tree sequences provide an extremely efficient way to store *and* to access genetic data (i.e. they are [succinct data structures](https://en.wikipedia.org/wiki/Succinct_data_structure))\n",
    "\n",
    "2. Tree sequences store ancestry rather than variation data (although the second can be trivially generated from the first). In other words, they reflect the evolutionary process rather than the genetic sequence which results from this process. Not only is this more meaningful, but for simulated tree sequences it also results in precision and efficiency gains, and for inferred tree sequences we hope it will lead to novel analytical approaches.\n",
    "\n",
    "3.  There now exist efficient tools for evolutionary simulation and ancestral inference which produce tree sequences. In particular, these include the following\n",
    "    -  Software for *simulating* tree sequences:\n",
    "        -  Reverse time (coalescent) simulators: **`msprime`**\n",
    "        -  Forward time simulators: **`SLiM`** (easiest), **`fwdpp`** (experimental, faster)\n",
    "    -  Software for *inferring* tree sequences from real data\n",
    "        -  **`tsinfer`**\n",
    "\n",
    "4. There are built-in routines that simplify tree sequences to the minimum knowable ancestral history. This turns out to be a powerful operation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "L7_4gmLcaei2"
   },
   "source": [
    "## Terminology\n",
    "\n",
    "A tree sequence is, as you might expect, a set of evolutionary trees, together with some associated data that specifies, for example, which branches in the trees contain individual mutations. Tree sequences can be accessed using the `tskit` python library. <aside>Currently the tskit functionality is part of the msprime library, it can be accessed by specifying `import msprime as tskit` at the top of your python script</aside>\n",
    "\n",
    "### Nodes\n",
    "Trees in a tree sequence are represented in a slightly unusual way. In the majority of libraries dealing with trees, each node is represented as an object in memory and the relationship between nodes as pointers between these objects. In tree sequences, however, all nodes, both internal and terminal, are represented by an integer ID, unique over the entire tree sequence.\n",
    "\n",
    "The tips (leaves) of the tree, are usually marked as *sample nodes*: In most cases tips will correspond to samples with a full genetic sequence\n",
    "\n",
    "Internal nodes represent ancestors in which a single DNA sequence was duplicated (in forwards-time terminology) or in which multiple sequences coalesced (in backwards-time terminology). \n",
    "\n",
    "### Edges\n",
    "Edges specify a directional connection between node IDs in the tree. Since trees in a tree sequence are usually taken to represent marginal trees along a genome, as well as branches on the tree having a length (a time dimension), each branch also has a *span*, giving the length of a genome it covers.\n",
    "\n",
    "### Genomes\n",
    "A tree sequence is usually taken to cover a contiguous section of a single chromosome. For clarity of explanation, we use the word \"genome\" to refer to the section spanned by the whole tree sequence (if you find it easier, imagine that the examples in this notebook refer to to an organism with a single, short chromosome)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cXUuOh1hEC_f"
   },
   "source": [
    "The tree below shows a pre-loaded example of a tree from a tree sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "cellView": "both",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 221
    },
    "colab_type": "code",
    "id": "BPjJtzqZ3-MP",
    "outputId": "6f79b67c-7799-4d91-b463-8fd96771a295"
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"200\" version=\"1.1\" width=\"200\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs/><g id=\"lines\" stroke=\"black\"><line x1=\"37.5\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"20.0\"/><line x1=\"37.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\"/><line x1=\"25.0\" x2=\"25.0\" y1=\"180.0\" y2=\"156.91628089290032\"/><line x1=\"25.0\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\"/><line x1=\"50.0\" x2=\"50.0\" y1=\"180.0\" y2=\"156.91628089290032\"/><line x1=\"50.0\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\"/><line x1=\"112.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"20.0\"/><line x1=\"112.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\"/><line x1=\"87.5\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"95.28485719410656\"/><line x1=\"87.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"95.28485719410656\"/><line x1=\"75.0\" x2=\"75.0\" y1=\"180.0\" y2=\"164.32147375967384\"/><line x1=\"75.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\"/><line x1=\"100.0\" x2=\"100.0\" y1=\"180.0\" y2=\"164.32147375967384\"/><line x1=\"100.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\"/><line x1=\"137.5\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"95.28485719410656\"/><line x1=\"137.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"95.28485719410656\"/><line x1=\"125.0\" x2=\"125.0\" y1=\"180.0\" y2=\"128.866241702499\"/><line x1=\"125.0\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"128.866241702499\"/><line x1=\"150.0\" x2=\"150.0\" y1=\"180.0\" y2=\"128.866241702499\"/><line x1=\"150.0\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"128.866241702499\"/></g><g font-size=\"14\" text-anchor=\"start\"><text x=\"117.5\" y=\"90.28485719410656\">11</text><text x=\"142.5\" y=\"123.866241702499\">9</text></g><g font-size=\"14\" text-anchor=\"end\"><text x=\"32.5\" y=\"151.91628089290032\">7</text><text x=\"82.5\" y=\"159.32147375967384\">6</text></g><g font-size=\"14\" text-anchor=\"middle\"><text x=\"75.0\" y=\"15.0\">12</text><text x=\"25.0\" y=\"200.0\">0</text><text x=\"50.0\" y=\"200.0\">1</text><text x=\"75.0\" y=\"200.0\">4</text><text x=\"100.0\" y=\"200.0\">5</text><text x=\"125.0\" y=\"200.0\">2</text><text x=\"150.0\" y=\"200.0\">3</text></g><circle cx=\"75.0\" cy=\"20.0\" fill=\"black\" r=\"3\"/><circle cx=\"37.5\" cy=\"156.91628089290032\" fill=\"black\" r=\"3\"/><circle cx=\"25.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"50.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"112.5\" cy=\"95.28485719410656\" fill=\"black\" r=\"3\"/><circle cx=\"87.5\" cy=\"164.32147375967384\" fill=\"black\" r=\"3\"/><circle cx=\"75.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"100.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"137.5\" cy=\"128.866241702499\" fill=\"black\" r=\"3\"/><circle cx=\"125.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"150.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"start\"/><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"end\"/></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Ignore this code for the moment: it loads a tree sequence from an existing file and plots the 1st tree\n",
    "my_tree_seq_object = tskit.load(\"data/tree_seq_file1.trees\")\n",
    "SVG(my_tree_seq_object.first().draw())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "hev6coAfUlhr"
   },
   "source": [
    "A genealogical tree like this one records the lines of descent along which genomes have been inherited. For example, ancestral genome 7 gave rise to two lineages which ended in sample genome 0 and sample genome 1 in the current population. \n",
    "\n",
    "#### A note on ploidy\n",
    "Since diploids have two copies of each autosomal chromosome, diploid individuals contain two such lines of descent. The tree above might therefore representing three diploid (or 6 haploid) individuals. <aside>Note that tree sequences do have a concept of an individual that can contain multiple genomes (see [here](https://msprime.readthedocs.io/en/stable/interchange.html#nodes-genomes-or-individuals)), but currently this is only used when attaching metadata to a tree sequence (something detailed in part 3: [Metadata in tree sequences](ts_workshop_part3.ipynb#metadata)).</aside>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "g96G8L2Haei3"
   },
   "source": [
    "### Trees\n",
    "The example above is a single tree. We might wish to inspect all the trees in a tree sequence. The simpest way to do this is to print them one after another. Our pre-loaded example was actually from a 1kilobase tree sequence which happens to have 3 trees along the genome, one to the left, one in the middle of the DNA region, and one to the right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "cellView": "both",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 268
    },
    "colab_type": "code",
    "id": "rYsNAhmSaei3",
    "outputId": "8cb702b7-a8f8-43bb-af6a-3427af4057a3",
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border:0; text-align: center\"><tr><td><svg baseProfile=\"full\" height=\"200\" version=\"1.1\" width=\"200\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs /><g id=\"lines\" stroke=\"black\"><line x1=\"37.5\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"20.0\" /><line x1=\"37.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\" /><line x1=\"25.0\" x2=\"25.0\" y1=\"180.0\" y2=\"156.91628089290032\" /><line x1=\"25.0\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\" /><line x1=\"50.0\" x2=\"50.0\" y1=\"180.0\" y2=\"156.91628089290032\" /><line x1=\"50.0\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\" /><line x1=\"112.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"20.0\" /><line x1=\"112.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\" /><line x1=\"87.5\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"95.28485719410656\" /><line x1=\"87.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"95.28485719410656\" /><line x1=\"75.0\" x2=\"75.0\" y1=\"180.0\" y2=\"164.32147375967384\" /><line x1=\"75.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\" /><line x1=\"100.0\" x2=\"100.0\" y1=\"180.0\" y2=\"164.32147375967384\" /><line x1=\"100.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\" /><line x1=\"137.5\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"95.28485719410656\" /><line x1=\"137.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"95.28485719410656\" /><line x1=\"125.0\" x2=\"125.0\" y1=\"180.0\" y2=\"128.866241702499\" /><line x1=\"125.0\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"128.866241702499\" /><line x1=\"150.0\" x2=\"150.0\" y1=\"180.0\" y2=\"128.866241702499\" /><line x1=\"150.0\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"128.866241702499\" /></g><g font-size=\"14\" text-anchor=\"start\"><text x=\"117.5\" y=\"90.28485719410656\">11</text><text x=\"142.5\" y=\"123.866241702499\">9</text></g><g font-size=\"14\" text-anchor=\"end\"><text x=\"32.5\" y=\"151.91628089290032\">7</text><text x=\"82.5\" y=\"159.32147375967384\">6</text></g><g font-size=\"14\" text-anchor=\"middle\"><text x=\"75.0\" y=\"15.0\">12</text><text x=\"25.0\" y=\"200.0\">0</text><text x=\"50.0\" y=\"200.0\">1</text><text x=\"75.0\" y=\"200.0\">4</text><text x=\"100.0\" y=\"200.0\">5</text><text x=\"125.0\" y=\"200.0\">2</text><text x=\"150.0\" y=\"200.0\">3</text></g><circle cx=\"75.0\" cy=\"20.0\" fill=\"black\" r=\"3\" /><circle cx=\"37.5\" cy=\"156.91628089290032\" fill=\"black\" r=\"3\" /><circle cx=\"25.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"50.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"112.5\" cy=\"95.28485719410656\" fill=\"black\" r=\"3\" /><circle cx=\"87.5\" cy=\"164.32147375967384\" fill=\"black\" r=\"3\" /><circle cx=\"75.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"100.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"137.5\" cy=\"128.866241702499\" fill=\"black\" r=\"3\" /><circle cx=\"125.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"150.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"start\" /><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"end\" /></svg></td><td><svg baseProfile=\"full\" height=\"200\" version=\"1.1\" width=\"200\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs /><g id=\"lines\" stroke=\"black\"><line x1=\"37.5\" x2=\"37.5\" y1=\"128.866241702499\" y2=\"20.0\" /><line x1=\"37.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\" /><line x1=\"25.0\" x2=\"25.0\" y1=\"180.0\" y2=\"128.866241702499\" /><line x1=\"25.0\" x2=\"37.5\" y1=\"128.866241702499\" y2=\"128.866241702499\" /><line x1=\"50.0\" x2=\"50.0\" y1=\"180.0\" y2=\"128.866241702499\" /><line x1=\"50.0\" x2=\"37.5\" y1=\"128.866241702499\" y2=\"128.866241702499\" /><line x1=\"112.5\" x2=\"112.5\" y1=\"127.25197685970193\" y2=\"20.0\" /><line x1=\"112.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\" /><line x1=\"87.5\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"127.25197685970193\" /><line x1=\"87.5\" x2=\"112.5\" y1=\"127.25197685970193\" y2=\"127.25197685970193\" /><line x1=\"75.0\" x2=\"75.0\" y1=\"180.0\" y2=\"164.32147375967384\" /><line x1=\"75.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\" /><line x1=\"100.0\" x2=\"100.0\" y1=\"180.0\" y2=\"164.32147375967384\" /><line x1=\"100.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\" /><line x1=\"137.5\" x2=\"137.5\" y1=\"156.91628089290032\" y2=\"127.25197685970193\" /><line x1=\"137.5\" x2=\"112.5\" y1=\"127.25197685970193\" y2=\"127.25197685970193\" /><line x1=\"125.0\" x2=\"125.0\" y1=\"180.0\" y2=\"156.91628089290032\" /><line x1=\"125.0\" x2=\"137.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\" /><line x1=\"150.0\" x2=\"150.0\" y1=\"180.0\" y2=\"156.91628089290032\" /><line x1=\"150.0\" x2=\"137.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\" /></g><g font-size=\"14\" text-anchor=\"start\"><text x=\"117.5\" y=\"122.25197685970193\">10</text><text x=\"142.5\" y=\"151.91628089290032\">7</text></g><g font-size=\"14\" text-anchor=\"end\"><text x=\"32.5\" y=\"123.866241702499\">9</text><text x=\"82.5\" y=\"159.32147375967384\">6</text></g><g font-size=\"14\" text-anchor=\"middle\"><text x=\"75.0\" y=\"15.0\">12</text><text x=\"25.0\" y=\"200.0\">2</text><text x=\"50.0\" y=\"200.0\">3</text><text x=\"75.0\" y=\"200.0\">4</text><text x=\"100.0\" y=\"200.0\">5</text><text x=\"125.0\" y=\"200.0\">0</text><text x=\"150.0\" y=\"200.0\">1</text></g><circle cx=\"75.0\" cy=\"20.0\" fill=\"black\" r=\"3\" /><circle cx=\"37.5\" cy=\"128.866241702499\" fill=\"black\" r=\"3\" /><circle cx=\"25.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"50.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"112.5\" cy=\"127.25197685970193\" fill=\"black\" r=\"3\" /><circle cx=\"87.5\" cy=\"164.32147375967384\" fill=\"black\" r=\"3\" /><circle cx=\"75.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"100.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"137.5\" cy=\"156.91628089290032\" fill=\"black\" r=\"3\" /><circle cx=\"125.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"150.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"start\" /><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"end\" /></svg></td><td><svg baseProfile=\"full\" height=\"200\" version=\"1.1\" width=\"200\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs /><g id=\"lines\" stroke=\"black\"><line x1=\"62.5\" x2=\"62.5\" y1=\"127.25197685970193\" y2=\"20.0\" /><line x1=\"62.5\" x2=\"100.0\" y1=\"20.0\" y2=\"20.0\" /><line x1=\"37.5\" x2=\"37.5\" y1=\"164.32147375967384\" y2=\"127.25197685970193\" /><line x1=\"37.5\" x2=\"62.5\" y1=\"127.25197685970193\" y2=\"127.25197685970193\" /><line x1=\"25.0\" x2=\"25.0\" y1=\"180.0\" y2=\"164.32147375967384\" /><line x1=\"25.0\" x2=\"37.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\" /><line x1=\"50.0\" x2=\"50.0\" y1=\"180.0\" y2=\"164.32147375967384\" /><line x1=\"50.0\" x2=\"37.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\" /><line x1=\"87.5\" x2=\"87.5\" y1=\"156.91628089290032\" y2=\"127.25197685970193\" /><line x1=\"87.5\" x2=\"62.5\" y1=\"127.25197685970193\" y2=\"127.25197685970193\" /><line x1=\"75.0\" x2=\"75.0\" y1=\"180.0\" y2=\"156.91628089290032\" /><line x1=\"75.0\" x2=\"87.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\" /><line x1=\"100.0\" x2=\"100.0\" y1=\"180.0\" y2=\"156.91628089290032\" /><line x1=\"100.0\" x2=\"87.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\" /><line x1=\"137.5\" x2=\"137.5\" y1=\"138.6226520567097\" y2=\"20.0\" /><line x1=\"137.5\" x2=\"100.0\" y1=\"20.0\" y2=\"20.0\" /><line x1=\"125.0\" x2=\"125.0\" y1=\"180.0\" y2=\"138.6226520567097\" /><line x1=\"125.0\" x2=\"137.5\" y1=\"138.6226520567097\" y2=\"138.6226520567097\" /><line x1=\"150.0\" x2=\"150.0\" y1=\"180.0\" y2=\"138.6226520567097\" /><line x1=\"150.0\" x2=\"137.5\" y1=\"138.6226520567097\" y2=\"138.6226520567097\" /></g><g font-size=\"14\" text-anchor=\"start\"><text x=\"92.5\" y=\"151.91628089290032\">7</text><text x=\"142.5\" y=\"133.6226520567097\">8</text></g><g font-size=\"14\" text-anchor=\"end\"><text x=\"57.5\" y=\"122.25197685970193\">10</text><text x=\"32.5\" y=\"159.32147375967384\">6</text></g><g font-size=\"14\" text-anchor=\"middle\"><text x=\"100.0\" y=\"15.0\">12</text><text x=\"25.0\" y=\"200.0\">4</text><text x=\"50.0\" y=\"200.0\">5</text><text x=\"75.0\" y=\"200.0\">0</text><text x=\"100.0\" y=\"200.0\">1</text><text x=\"125.0\" y=\"200.0\">2</text><text x=\"150.0\" y=\"200.0\">3</text></g><circle cx=\"100.0\" cy=\"20.0\" fill=\"black\" r=\"3\" /><circle cx=\"62.5\" cy=\"127.25197685970193\" fill=\"black\" r=\"3\" /><circle cx=\"37.5\" cy=\"164.32147375967384\" fill=\"black\" r=\"3\" /><circle cx=\"25.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"50.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"87.5\" cy=\"156.91628089290032\" fill=\"black\" r=\"3\" /><circle cx=\"75.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"100.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"137.5\" cy=\"138.6226520567097\" fill=\"black\" r=\"3\" /><circle cx=\"125.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><circle cx=\"150.0\" cy=\"180.0\" fill=\"black\" r=\"3\" /><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"start\" /><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"end\" /></svg></td></tr><tr><td colspan=\"3\"><hr></td></tr><tr style=\"font-weight: bold; font-family: serif\"><td style=\"text-align: left\">0</td><td style=\"text-align: center\">genome position</td><td style=\"text-align: right\">1000</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Ignore this code for the moment\n",
    "svg_trees=[tree.draw() for tree in my_tree_seq_object.trees()]\n",
    "display(HTML('<table style=\"border:0; text-align: center\"><tr><td>{}</td></tr><tr><td colspan=\"{}\">{}</td></tr><tr style=\"font-weight: bold; font-family: serif\"><td style=\"text-align: left\">{}</td><td style=\"text-align: center\">{}</td><td style=\"text-align: right\">{}</td></tr></table>'.format(\n",
    "    \"</td><td>\".join(svg_trees),len(svg_trees),\"<hr>\",0, \"genome position\", 1000)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ST__17zM7lem"
   },
   "source": [
    "Of course, the *sample nodes* 0..5 are present in all trees (i.e. at all genomic positions) since we have their full genomes. What you should note is that most of the *internal nodes* are also present in multiple trees: the exceptions are node 11 which is only present in the first tree, node 9 which is only present in the first two trees, and node 8 which is only present in the last tree.\n",
    "\n",
    "However, in tree sequence terminology, we don't explictly keep track of where nodes start and end. Only edges (not nodes) possess a genomic span. So for example, this tree sequence is defined using edges like (0->7) which spans the entire genome, edges like (9->11) which only span the leftmost section of the genome, and edges like (10->12) which span the remain part of the genome. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "LQ5iEQKFDGiN"
   },
   "source": [
    "###Number of edges: a key tree sequence measure\n",
    "Because edges can span multiple trees, they are an efficient way of storing ancestral information. In fact, the size of a tree sequence file is largely determined by the number of edges. This number can be used as a rough measure of many things such as the compressibility of a simulation, the speed with which analyses are likely to run, and even the effectiveness of ancestral inference. To give some sense of scale our inferred ancestry of human chromosome 20 in the UK Biobank dataset of a million genomes consists of 62 million edges which define 16 thousand trees; the uncompressed tree sequence takes 2.1 gigabytes of disk space, and on a reasonably specced computer it takes ~11 seconds to iterate over all these trees.\n",
    "\n",
    "The `tskit` library allow us to use tree sequence objects in Python; these have specific properties and can also be queried using various methods. Currently these are documented in the [msprime API documentation](https://msprime.readthedocs.io/en/stable/api.html). For example, we can easily show the number of edges in a tree sequence as follows "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "XNT1c8m00OpS",
    "outputId": "02df15be-9f53-437d-aa1a-8f1fcc831e7a"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_tree_seq_object.num_edges # a lot less than the UK Biobank!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "18CLa9sTKwAM"
   },
   "source": [
    "and we can iterate over the edges using the `.edges()` iterator, which returns an `Edge` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 306
    },
    "colab_type": "code",
    "id": "rfDkbp_4KvO3",
    "outputId": "a693b02a-4d7b-44a3-e7cd-9409e9cce3e5"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{left=0.000, right=1000.000, parent=6, child=4}\n",
      "{left=0.000, right=1000.000, parent=6, child=5}\n",
      "{left=0.000, right=1000.000, parent=7, child=0}\n",
      "{left=0.000, right=1000.000, parent=7, child=1}\n",
      "{left=480.525, right=1000.000, parent=8, child=2}\n",
      "{left=480.525, right=1000.000, parent=8, child=3}\n",
      "{left=0.000, right=480.525, parent=9, child=2}\n",
      "{left=0.000, right=480.525, parent=9, child=3}\n",
      "{left=249.237, right=1000.000, parent=10, child=6}\n",
      "{left=249.237, right=1000.000, parent=10, child=7}\n",
      "{left=0.000, right=249.237, parent=11, child=6}\n",
      "{left=0.000, right=249.237, parent=11, child=9}\n",
      "{left=0.000, right=249.237, parent=12, child=7}\n",
      "{left=480.525, right=1000.000, parent=12, child=8}\n",
      "{left=249.237, right=480.525, parent=12, child=9}\n",
      "{left=249.237, right=1000.000, parent=12, child=10}\n",
      "{left=0.000, right=249.237, parent=12, child=11}\n",
      "[    0.             0.             0.             0.\n",
      "     0.             0.          2311.62912472  3403.44472285\n",
      "  6100.64244197  7539.11962928  7777.12552145 12490.33157918\n",
      " 23590.26953715]\n"
     ]
    }
   ],
   "source": [
    "for edge in my_tree_seq_object.edges():\n",
    "  print(edge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly we can iterate over other components of a tree sequence, like the trees themselves. The `.trees()` iterator  returns `SparseTree` objects, which we can `.draw()` in SVG or in unicode text format; this is how we produced the plots above. <aside>To encourage efficiency, there no in-built provision to jump directly into the middle of a tree sequence. To obtain a tree at a particular point in the genome, you need to iterate through the trees until you identify the one that you want.</aside>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"200\" version=\"1.1\" width=\"200\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs/><g id=\"lines\" stroke=\"black\"><line x1=\"37.5\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"20.0\"/><line x1=\"37.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\"/><line x1=\"25.0\" x2=\"25.0\" y1=\"180.0\" y2=\"156.91628089290032\"/><line x1=\"25.0\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\"/><line x1=\"50.0\" x2=\"50.0\" y1=\"180.0\" y2=\"156.91628089290032\"/><line x1=\"50.0\" x2=\"37.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\"/><line x1=\"112.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"20.0\"/><line x1=\"112.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\"/><line x1=\"87.5\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"95.28485719410656\"/><line x1=\"87.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"95.28485719410656\"/><line x1=\"75.0\" x2=\"75.0\" y1=\"180.0\" y2=\"164.32147375967384\"/><line x1=\"75.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\"/><line x1=\"100.0\" x2=\"100.0\" y1=\"180.0\" y2=\"164.32147375967384\"/><line x1=\"100.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\"/><line x1=\"137.5\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"95.28485719410656\"/><line x1=\"137.5\" x2=\"112.5\" y1=\"95.28485719410656\" y2=\"95.28485719410656\"/><line x1=\"125.0\" x2=\"125.0\" y1=\"180.0\" y2=\"128.866241702499\"/><line x1=\"125.0\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"128.866241702499\"/><line x1=\"150.0\" x2=\"150.0\" y1=\"180.0\" y2=\"128.866241702499\"/><line x1=\"150.0\" x2=\"137.5\" y1=\"128.866241702499\" y2=\"128.866241702499\"/></g><g font-size=\"14\" text-anchor=\"start\"><text x=\"117.5\" y=\"90.28485719410656\">11</text><text x=\"142.5\" y=\"123.866241702499\">9</text></g><g font-size=\"14\" text-anchor=\"end\"><text x=\"32.5\" y=\"151.91628089290032\">7</text><text x=\"82.5\" y=\"159.32147375967384\">6</text></g><g font-size=\"14\" text-anchor=\"middle\"><text x=\"75.0\" y=\"15.0\">12</text><text x=\"25.0\" y=\"200.0\">0</text><text x=\"50.0\" y=\"200.0\">1</text><text x=\"75.0\" y=\"200.0\">4</text><text x=\"100.0\" y=\"200.0\">5</text><text x=\"125.0\" y=\"200.0\">2</text><text x=\"150.0\" y=\"200.0\">3</text></g><circle cx=\"75.0\" cy=\"20.0\" fill=\"black\" r=\"3\"/><circle cx=\"37.5\" cy=\"156.91628089290032\" fill=\"black\" r=\"3\"/><circle cx=\"25.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"50.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"112.5\" cy=\"95.28485719410656\" fill=\"black\" r=\"3\"/><circle cx=\"87.5\" cy=\"164.32147375967384\" fill=\"black\" r=\"3\"/><circle cx=\"75.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"100.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"137.5\" cy=\"128.866241702499\" fill=\"black\" r=\"3\"/><circle cx=\"125.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"150.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"start\"/><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"end\"/></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"200\" version=\"1.1\" width=\"200\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs/><g id=\"lines\" stroke=\"black\"><line x1=\"37.5\" x2=\"37.5\" y1=\"128.866241702499\" y2=\"20.0\"/><line x1=\"37.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\"/><line x1=\"25.0\" x2=\"25.0\" y1=\"180.0\" y2=\"128.866241702499\"/><line x1=\"25.0\" x2=\"37.5\" y1=\"128.866241702499\" y2=\"128.866241702499\"/><line x1=\"50.0\" x2=\"50.0\" y1=\"180.0\" y2=\"128.866241702499\"/><line x1=\"50.0\" x2=\"37.5\" y1=\"128.866241702499\" y2=\"128.866241702499\"/><line x1=\"112.5\" x2=\"112.5\" y1=\"127.25197685970193\" y2=\"20.0\"/><line x1=\"112.5\" x2=\"75.0\" y1=\"20.0\" y2=\"20.0\"/><line x1=\"87.5\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"127.25197685970193\"/><line x1=\"87.5\" x2=\"112.5\" y1=\"127.25197685970193\" y2=\"127.25197685970193\"/><line x1=\"75.0\" x2=\"75.0\" y1=\"180.0\" y2=\"164.32147375967384\"/><line x1=\"75.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\"/><line x1=\"100.0\" x2=\"100.0\" y1=\"180.0\" y2=\"164.32147375967384\"/><line x1=\"100.0\" x2=\"87.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\"/><line x1=\"137.5\" x2=\"137.5\" y1=\"156.91628089290032\" y2=\"127.25197685970193\"/><line x1=\"137.5\" x2=\"112.5\" y1=\"127.25197685970193\" y2=\"127.25197685970193\"/><line x1=\"125.0\" x2=\"125.0\" y1=\"180.0\" y2=\"156.91628089290032\"/><line x1=\"125.0\" x2=\"137.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\"/><line x1=\"150.0\" x2=\"150.0\" y1=\"180.0\" y2=\"156.91628089290032\"/><line x1=\"150.0\" x2=\"137.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\"/></g><g font-size=\"14\" text-anchor=\"start\"><text x=\"117.5\" y=\"122.25197685970193\">10</text><text x=\"142.5\" y=\"151.91628089290032\">7</text></g><g font-size=\"14\" text-anchor=\"end\"><text x=\"32.5\" y=\"123.866241702499\">9</text><text x=\"82.5\" y=\"159.32147375967384\">6</text></g><g font-size=\"14\" text-anchor=\"middle\"><text x=\"75.0\" y=\"15.0\">12</text><text x=\"25.0\" y=\"200.0\">2</text><text x=\"50.0\" y=\"200.0\">3</text><text x=\"75.0\" y=\"200.0\">4</text><text x=\"100.0\" y=\"200.0\">5</text><text x=\"125.0\" y=\"200.0\">0</text><text x=\"150.0\" y=\"200.0\">1</text></g><circle cx=\"75.0\" cy=\"20.0\" fill=\"black\" r=\"3\"/><circle cx=\"37.5\" cy=\"128.866241702499\" fill=\"black\" r=\"3\"/><circle cx=\"25.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"50.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"112.5\" cy=\"127.25197685970193\" fill=\"black\" r=\"3\"/><circle cx=\"87.5\" cy=\"164.32147375967384\" fill=\"black\" r=\"3\"/><circle cx=\"75.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"100.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"137.5\" cy=\"156.91628089290032\" fill=\"black\" r=\"3\"/><circle cx=\"125.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"150.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"start\"/><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"end\"/></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"200\" version=\"1.1\" width=\"200\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs/><g id=\"lines\" stroke=\"black\"><line x1=\"62.5\" x2=\"62.5\" y1=\"127.25197685970193\" y2=\"20.0\"/><line x1=\"62.5\" x2=\"100.0\" y1=\"20.0\" y2=\"20.0\"/><line x1=\"37.5\" x2=\"37.5\" y1=\"164.32147375967384\" y2=\"127.25197685970193\"/><line x1=\"37.5\" x2=\"62.5\" y1=\"127.25197685970193\" y2=\"127.25197685970193\"/><line x1=\"25.0\" x2=\"25.0\" y1=\"180.0\" y2=\"164.32147375967384\"/><line x1=\"25.0\" x2=\"37.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\"/><line x1=\"50.0\" x2=\"50.0\" y1=\"180.0\" y2=\"164.32147375967384\"/><line x1=\"50.0\" x2=\"37.5\" y1=\"164.32147375967384\" y2=\"164.32147375967384\"/><line x1=\"87.5\" x2=\"87.5\" y1=\"156.91628089290032\" y2=\"127.25197685970193\"/><line x1=\"87.5\" x2=\"62.5\" y1=\"127.25197685970193\" y2=\"127.25197685970193\"/><line x1=\"75.0\" x2=\"75.0\" y1=\"180.0\" y2=\"156.91628089290032\"/><line x1=\"75.0\" x2=\"87.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\"/><line x1=\"100.0\" x2=\"100.0\" y1=\"180.0\" y2=\"156.91628089290032\"/><line x1=\"100.0\" x2=\"87.5\" y1=\"156.91628089290032\" y2=\"156.91628089290032\"/><line x1=\"137.5\" x2=\"137.5\" y1=\"138.6226520567097\" y2=\"20.0\"/><line x1=\"137.5\" x2=\"100.0\" y1=\"20.0\" y2=\"20.0\"/><line x1=\"125.0\" x2=\"125.0\" y1=\"180.0\" y2=\"138.6226520567097\"/><line x1=\"125.0\" x2=\"137.5\" y1=\"138.6226520567097\" y2=\"138.6226520567097\"/><line x1=\"150.0\" x2=\"150.0\" y1=\"180.0\" y2=\"138.6226520567097\"/><line x1=\"150.0\" x2=\"137.5\" y1=\"138.6226520567097\" y2=\"138.6226520567097\"/></g><g font-size=\"14\" text-anchor=\"start\"><text x=\"92.5\" y=\"151.91628089290032\">7</text><text x=\"142.5\" y=\"133.6226520567097\">8</text></g><g font-size=\"14\" text-anchor=\"end\"><text x=\"57.5\" y=\"122.25197685970193\">10</text><text x=\"32.5\" y=\"159.32147375967384\">6</text></g><g font-size=\"14\" text-anchor=\"middle\"><text x=\"100.0\" y=\"15.0\">12</text><text x=\"25.0\" y=\"200.0\">4</text><text x=\"50.0\" y=\"200.0\">5</text><text x=\"75.0\" y=\"200.0\">0</text><text x=\"100.0\" y=\"200.0\">1</text><text x=\"125.0\" y=\"200.0\">2</text><text x=\"150.0\" y=\"200.0\">3</text></g><circle cx=\"100.0\" cy=\"20.0\" fill=\"black\" r=\"3\"/><circle cx=\"62.5\" cy=\"127.25197685970193\" fill=\"black\" r=\"3\"/><circle cx=\"37.5\" cy=\"164.32147375967384\" fill=\"black\" r=\"3\"/><circle cx=\"25.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"50.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"87.5\" cy=\"156.91628089290032\" fill=\"black\" r=\"3\"/><circle cx=\"75.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"100.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"137.5\" cy=\"138.6226520567097\" fill=\"black\" r=\"3\"/><circle cx=\"125.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><circle cx=\"150.0\" cy=\"180.0\" fill=\"black\" r=\"3\"/><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"start\"/><g alignment-baseline=\"middle\" font-size=\"14\" font-style=\"italic\" text-anchor=\"end\"/></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for tree in my_tree_seq_object.trees():\n",
    "    display(SVG(tree.draw()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zixkspB9DQ5y"
   },
   "source": [
    "### Tree changes, ancestral recombinations, and SPRs\n",
    "\n",
    "Note that adjacent trees in this tree sequence differ by a single \"tree edit\" or SPR (subtree-prune-and-regraft) operation: the subtree defined by node 6 in the first tree has been pruned and regrafted onto the edge from (12->7) to create the second tree. See if you can work out the SPR that converts the second tree to the third tree... \n",
    "\n",
    "In general, each detectable recombination occurring in ancestral history results in a single SPR in the tree sequence. If, as above, recombination breakpoints occurs at unique positions (an \"infinite sites\" model of breakpoints), then the number of trees in a tree sequence equals the number of ancestral recombination events minus one. If recombinations can occur at the same physical position (e.g. if a simulation is [run using a recombination map with finite loci](https://msprime.readthedocs.io/en/stable/api.html#msprime.RecombinationMap.uniform_map), then moving from one tree to the next in a tree sequence might require multiple SPRs, because there are multiple, overlaid ancestral recombination events.\n",
    "\n",
    "If you are moving from tree to tree, but only need to work out what changed each time, you can use the `.edge_diffs()` method on the tree sequence. We will encounter this later, when calculating incremental statistics over a tree sequence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SbYTYkpiMShU"
   },
   "source": [
    "## Tables: the underlying data store\n",
    "\n",
    "The underlying data in tree sequences is stored as a set of tables. For example, there is a Nodes table and an Edges table, as well as others like Sites, Mutations, Populations, etc. The API for accessing the underlying tables of a tree sequence is documented [here](https://msprime.readthedocs.io/en/stable/api.html#tables-api).\n",
    "\n",
    "Accessing tree sequence data via tables is an alternative to using the built-in methods and properties of `TreeSequence` objects: it requires a deeper working knowledge of the tree sequence implementation, but it is more powerful. For example, it can be used to edit tree sequences, which are otherwise immutable objects. As the documentation says:  \n",
    "<blockquote>To change the data held in a particular tree sequence, first get the table information as a <code>TableCollection</code> instance (using <code>dump_tables()</code>), edit those tables using the tables api, and create a new tree sequence using <code>TableCollection.tree_sequence()</code>.</blockquote>\n",
    "\n",
    "In summary, since a tree sequence object is entirely specified by tables, these provide access to the guts of a stored tree sequence. Below is a simple example of inspecting all tables; later in this workshop we will meet examples where accessing the raw table data allows simple and fast extraction of statistics for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "syNxQgNLaei8"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "############################################################\n",
      "#   Individuals                                            #\n",
      "############################################################\n",
      "id\tflags\tlocation\tmetadata\n",
      "############################################################\n",
      "#   Nodes                                                  #\n",
      "############################################################\n",
      "id\tflags\tpopulation\tindividual\ttime\tmetadata\n",
      "0\t1\t0\t-1\t0.00000000000000\t\n",
      "1\t1\t0\t-1\t0.00000000000000\t\n",
      "2\t1\t0\t-1\t0.00000000000000\t\n",
      "3\t1\t0\t-1\t0.00000000000000\t\n",
      "4\t1\t0\t-1\t0.00000000000000\t\n",
      "5\t1\t0\t-1\t0.00000000000000\t\n",
      "6\t0\t0\t-1\t2311.62912471615846\t\n",
      "7\t0\t0\t-1\t3403.44472285223856\t\n",
      "8\t0\t0\t-1\t6100.64244196680283\t\n",
      "9\t0\t0\t-1\t7539.11962928479534\t\n",
      "10\t0\t0\t-1\t7777.12552144685378\t\n",
      "11\t0\t0\t-1\t12490.33157918277175\t\n",
      "12\t0\t0\t-1\t23590.26953715075069\t\n",
      "############################################################\n",
      "#   Edges                                                  #\n",
      "############################################################\n",
      "id\tleft\t\tright\t\tparent\tchild\n",
      "0\t0.00000000\t1000.00000000\t6\t4\n",
      "1\t0.00000000\t1000.00000000\t6\t5\n",
      "2\t0.00000000\t1000.00000000\t7\t0\n",
      "3\t0.00000000\t1000.00000000\t7\t1\n",
      "4\t480.52523483\t1000.00000000\t8\t2\n",
      "5\t480.52523483\t1000.00000000\t8\t3\n",
      "6\t0.00000000\t480.52523483\t9\t2\n",
      "7\t0.00000000\t480.52523483\t9\t3\n",
      "8\t249.23677096\t1000.00000000\t10\t6\n",
      "9\t249.23677096\t1000.00000000\t10\t7\n",
      "10\t0.00000000\t249.23677096\t11\t6\n",
      "11\t0.00000000\t249.23677096\t11\t9\n",
      "12\t0.00000000\t249.23677096\t12\t7\n",
      "13\t480.52523483\t1000.00000000\t12\t8\n",
      "14\t249.23677096\t480.52523483\t12\t9\n",
      "15\t249.23677096\t1000.00000000\t12\t10\n",
      "16\t0.00000000\t249.23677096\t12\t11\n",
      "############################################################\n",
      "#   Sites                                                  #\n",
      "############################################################\n",
      "id\tposition\tancestral_state\tmetadata\n",
      "############################################################\n",
      "#   Mutations                                              #\n",
      "############################################################\n",
      "id\tsite\tnode\tderived_state\tparent\tmetadata\n",
      "############################################################\n",
      "#   Migrations                                             #\n",
      "############################################################\n",
      "id\tleft\tright\tnode\tsource\tdest\ttime\n",
      "############################################################\n",
      "#   Populations                                            #\n",
      "############################################################\n",
      "id\tmetadata\n",
      "0\t\n",
      "############################################################\n",
      "#   Provenances                                            #\n",
      "############################################################\n",
      "id\ttimestamp\trecord\n",
      "0\t2018-12-13T11:34:14.290625\t{\"schema_version\": \"1.0.0\", \"software\": {\"name\": \"msprime\", \"version\": \"0.6.2\"}, \"parameters\": {\"command\": \"simulate\", \"random_seed\": 52, \"TODO\": \"add other simulation parameters\"}, \"environment\": {\"libraries\": {\"gsl\": {\"version\": \"2.5\"}, \"kastore\": {\"version\": \"0.1.0\"}}, \"os\": {\"system\": \"Darwin\", \"node\": \"Yans-MacBook-Air.local\", \"release\": \"18.2.0\", \"version\": \"Darwin Kernel Version 18.2.0: Mon Nov 12 20:24:31 PST 2018; root:xnu-4903.231.4~1/RELEASE_X86_64\", \"machine\": \"x86_64\"}, \"python\": {\"implementation\": \"CPython\", \"version\": \"3.7.1\"}}}\n"
     ]
    }
   ],
   "source": [
    "print(my_tree_seq_object.tables)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "w9_aIYpXMITy"
   },
   "source": [
    "## Tree sequences and ARGs\n",
    "Much of the literature on full ancestral inference concentrates on the Ancestral Recombination Graph, or ARG. A tree sequence is essentially a compact representation of an ARG. The main difference is that ARGs store details of the position and potentially the timing of recombination events, whereas tree sequences simply store the trees created by these events. Many (but not all) details of the actual recombination events can be inferred from the changes in edges along a tree sequence. As we state in the [tsinfer paper](https://www.biorxiv.org/content/early/2018/11/01/458067)\n",
    "<blockquote>an ARG encodes the events that occurred in the history of a sample, whereas a tree sequence encodes the outcome of those events</blockquote>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "ts_workshop.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}