Tskit and R

Tskit and R#

To interface with tskit in R, we can use the reticulate R package, which lets you call Python functions within an R session. In this tutorial, we’ll go through a couple of examples to show you how to get started. If you haven’t done so already, you’ll need to install reticulate in your R session via install.packages("reticulate").

We’ll begin by simulating a small tree sequence using msprime.

msprime <- reticulate::import("msprime")

ts <- msprime$sim_ancestry(80, sequence_length=1e4, recombination_rate=1e-4, random_seed=42)
ts  # See "Jupyter notebook tips", below for how to render this nicely

<tskit.trees.TreeSequence object at 0x7f2dfc9dcf50>

Attributes and methods#

reticulate allows us to access a Python object’s attributes via the $ operator. For example, we can access (and assign to a variable) the number of samples in the tree sequence:

n <- ts$num_samples
n
160

The $ operator can also be used to call methods, for example, the simplify() method associated with the tree sequence. The method parameters are given as native R objects (but note that object IDs still use tskit’s 0-based indexing system).

reduced_ts <- ts$simplify(0:7)  # only keep samples with ids 0, 1, 2, 3, 4, 5, 6, 7
reduced_ts <- reduced_ts$delete_intervals(list(c(6000, 10000)))  # delete data after 6kb
reduced_ts <- reduced_ts$trim()  # remove the deleted region
paste(
    "Reduced from", ts$num_trees, "trees over", ts$sequence_length/1e3, "kb to",
    reduced_ts$num_trees, "trees over", reduced_ts$sequence_length/1e3, "kb.")
'Reduced from 26 trees over 10 kb to 6 trees over 6 kb.'

IDs and indexes#

Note that if a bare digit is provided to one of these methods, it will be treated as a floating point number. This is useful to know when calling tskit methods that require integers (e.g. object IDs). For example, the following will not work:

ts$node(0)  # Will raise an error

In this case, to force the 0 to be passed as an integer, you can either coerce it using as.integer or simply prepend the letter L:

ts$node(as.integer(0))
# or
ts$node(0L)

Node(id=0, flags=1, time=0.0, population=0, individual=0, metadata=b'')

Node(id=0, flags=1, time=0.0, population=0, individual=0, metadata=b'')

Coercing in this way is only necessary when passing parameters to those underlying tskit methods that expect integers. It is not needed e.g. to index into numeric arrays. However, when using arrays, very careful attention must be paid to the fact that tskit IDs start at zero, whereas R indexes start at one:

root_id <- ts$first()$root
paste("Root time via tskit method:", ts$node(root_id)$time)
# When indexing into tskit arrays in R, add 1 to the ID
paste("Root time via array access:", ts$nodes_time[root_id + 1])
'Root time via tskit method: 3.4634393266485'
'Root time via array access: 3.4634393266485'

Analysis#

From within R we can use tskit’s powerful Statistics framework to efficiently compute many different summary statistics from a tree sequence. To illustrate this, we’ll first add some mutations to our tree sequence with the msprime.sim_mutations() function, and then compute the genetic diversity for each of the tree sequence’s sample nodes:

ts_mut = msprime$sim_mutations(reduced_ts, rate=1e-4, random_seed=321)

paste(ts_mut$num_mutations, "mutations, genetic diversity is", ts_mut$diversity())
'13 mutations, genetic diversity is 0.000720238095238095'

Numerical arrays and matrices work as expected. For instance, we can use the tree sequence genotype_matrix() method to return the genotypes of the tree sequence as a matrix object in R.

G = ts_mut$genotype_matrix()
G
A matrix: 13 × 8 of type int
00010000
00000010
10000000
00100000
00110110
01110110
01111111
01111111
10000000
00001001
10000000
01110110
00001000

We can then use R functions directly on the genotype matrix:

allele_frequency = rowMeans(G)
allele_frequency
  1. 0.125
  2. 0.125
  3. 0.125
  4. 0.125
  5. 0.5
  6. 0.625
  7. 0.875
  8. 0.875
  9. 0.125
  10. 0.25
  11. 0.125
  12. 0.625
  13. 0.125

Jupyter notebook tips#

When running R within a Jupyter notebook, a few magic functions can be defined that allow tskit objects to be rendered within the notebook:

# Define some magic functions to allow objects to be displayed in R Jupyter notebooks
repr_html.tskit.trees.TreeSequence <- function(obj, ...){obj$`_repr_html_`()}
repr_html.tskit.trees.Tree <- function(obj, ...){obj$`_repr_html_`()}
repr_svg.tskit.drawing.SVGString <- function(obj, ...){obj$`__str__`()}

This leads to much nicer tabular summaries:

ts_mut
Tree Sequence
Trees6
Sequence Length6000.0
Time Unitsgenerations
Sample Nodes8
Total Size6.1 KiB
MetadataNo Metadata
Table Rows Size Has Metadata
Edges 32 1.0 KiB
Individuals 4 136 Bytes
Migrations 0 8 Bytes
Mutations 13 497 Bytes
Nodes 20 568 Bytes
Populations 1 224 Bytes ✅
Provenances 5 3.1 KiB
Sites 13 341 Bytes

It also allows trees and tree sequences to be plotted inline:

ts_mut$draw_svg(y_axis=TRUE, y_ticks=0:10)
'<svg baseProfile="full" height="200" version="1.1" width="1200" 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><style type="text/css"><![CDATA[.background path {fill: #808080; fill-opacity: 0}.background path:nth-child(odd) {fill-opacity: .1}.axes {font-size: 14px}.x-axis .tick .lab {font-weight: bold; dominant-baseline: hanging}.axes, .tree {font-size: 14px; text-anchor: middle}.axes line, .edge {stroke: black; fill: none}.axes .ax-skip {stroke-dasharray: 4}.y-axis .grid {stroke: #FAFAFA}.node > .sym {fill: black; stroke: none}.site > .sym {stroke: black}.mut text {fill: red; font-style: italic}.mut.extra text {fill: hotpink}.mut line {fill: none; stroke: none}.mut .sym {fill: none; stroke: red}.mut.extra .sym {stroke: hotpink}.node .mut .sym {stroke-width: 1.5px}.tree text, .tree-sequence text {dominant-baseline: central}.plotbox .lab.lft {text-anchor: end}.plotbox .lab.rgt {text-anchor: start}]]></style></defs><g class="tree-sequence"><g class="background"><path d="M56.8,0 l187.2,0 l0,138.2 l-144.144,25 l0,5 l-43.056,0 l0,-5 l0,-25 l0,-138.2z" /><path d="M244,0 l187.2,0 l0,138.2 l123.739,25 l0,5 l-455.083,0 l0,-5 l144.144,-25 l0,-138.2z" /><path d="M431.2,0 l187.2,0 l0,138.2 l217.714,25 l0,5 l-281.174,0 l0,-5 l-123.739,-25 l0,-138.2z" /><path d="M618.4,0 l187.2,0 l0,138.2 l184.954,25 l0,5 l-154.44,0 l0,-5 l-217.714,-25 l0,-138.2z" /><path d="M805.6,0 l187.2,0 l0,138.2 l129.168,25 l0,5 l-131.414,0 l0,-5 l-184.954,-25 l0,-138.2z" /><path d="M992.8,0 l187.2,0 l0,138.2 l0,25 l0,5 l-58.032,0 l0,-5 l-129.17,-25 l0,-138.2z" /></g><g class="axes"><g class="x-axis"><g class="title" transform="translate(618.4 200)"><text class="lab" text-anchor="middle" transform="translate(0 -11)">Genome position</text></g><line class="ax-line" x1="56.8" x2="1180" y1="163.2" y2="163.2" /><g class="ticks"><g class="tick" transform="translate(56.8 163.2)"><line x1="0" x2="0" y1="0" y2="5" /><g transform="translate(0 6)"><text class="lab">0</text></g></g><g class="tick" transform="translate(99.856 163.2)"><line x1="0" x2="0" y1="0" y2="5" /><g transform="translate(0 6)"><text class="lab">230</text></g></g><g class="tick" transform="translate(554.939 163.2)"><line x1="0" x2="0" y1="0" y2="5" /><g transform="translate(0 6)"><text class="lab">2661</text></g></g><g class="tick" transform="translate(836.114 163.2)"><line x1="0" x2="0" y1="0" y2="5" /><g transform="translate(0 6)"><text class="lab">4163</text></g></g><g class="tick" transform="translate(990.554 163.2)"><line x1="0" x2="0" y1="0" y2="5" /><g transform="translate(0 6)"><text class="lab">4988</text></g></g><g class="tick" transform="translate(1121.97 163.2)"><line x1="0" x2="0" y1="0" y2="5" /><g transform="translate(0 6)"><text class="lab">5690</text></g></g><g class="tick" transform="translate(1180 163.2)"><line x1="0" x2="0" y1="0" y2="5" /><g transform="translate(0 6)"><text class="lab">6000</text></g></g></g><g class="site s0" transform="translate(144.222 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m0"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s1" transform="translate(334.979 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m1"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s2" transform="translate(413.603 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m2"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s3" transform="translate(433.259 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m3"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s4" transform="translate(663.702 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m4"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s5" transform="translate(749.627 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m5"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s6" transform="translate(842.478 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m6"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s7" transform="translate(878.795 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m7"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s8" transform="translate(916.984 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m8"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s9" transform="translate(934.955 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m9"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s10" transform="translate(1059.63 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m10"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s11" transform="translate(1090.71 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m11"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g><g class="site s12" transform="translate(1107.37 163.2)"><line class="sym" x1="0" x2="0" y1="0" y2="-10" /><g class="mut m12"><polyline class="sym" points="2.5,-6.5 0,-1.5 -2.5,-6.5" /></g></g></g><g class="y-axis"><g class="title" transform="translate(0 65.7)"><text class="lab" text-anchor="middle" transform="translate(11) rotate(-90)">Time (generations)</text></g><line class="ax-line" x1="56.8" x2="56.8" y1="121.4" y2="10" /><g class="ticks"><g class="tick" transform="translate(56.8 121.4)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">0</text></g></g><g class="tick" transform="translate(56.8 108.736)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">1</text></g></g><g class="tick" transform="translate(56.8 96.073)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">2</text></g></g><g class="tick" transform="translate(56.8 83.4095)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">3</text></g></g><g class="tick" transform="translate(56.8 70.746)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">4</text></g></g><g class="tick" transform="translate(56.8 58.0824)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">5</text></g></g><g class="tick" transform="translate(56.8 45.4189)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">6</text></g></g><g class="tick" transform="translate(56.8 32.7554)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">7</text></g></g><g class="tick" transform="translate(56.8 20.0919)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">8</text></g></g><g class="tick" transform="translate(56.8 7.42839)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">9</text></g></g><g class="tick" transform="translate(56.8 -5.23512)"><line x1="0" x2="-5" y1="0" y2="0" /><g transform="translate(-6 0)"><text class="lab" text-anchor="end">10</text></g></g></g></g></g><g class="plotbox trees"><g class="tree t0" transform="translate(56.8 0)"><g class="plotbox"><g class="c2 node n18 p0 root" transform="translate(63.9875 77.5407)"><g class="a18 c2 node n13 p0" transform="translate(25.5875 30.7712)"><g class="a13 i0 leaf node n1 p0 sample" transform="translate(-23.575 13.0881)"><path class="edge" d="M 0 0 V -13.0881 H 23.575" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">1</text></g><g class="a13 c2 node n12 p0" transform="translate(23.575 2.03424)"><g class="a12 i1 leaf node n2 p0 sample" transform="translate(-28.75 11.0538)"><path class="edge" d="M 0 0 V -11.0538 H 28.75" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">2</text></g><g class="a12 c2 node n10 p0" transform="translate(28.75 5.24241)"><g class="a10 i2 leaf node n4 p0 sample" transform="translate(16.1 5.81144)"><path class="edge" d="M 0 0 V -5.81144 H -16.1" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">4</text></g><g class="a10 c2 node n9 p0" transform="translate(-16.1 1.39353)"><g class="a9 i2 leaf node n5 p0 sample" transform="translate(13.8 4.41791)"><path class="edge" d="M 0 0 V -4.41791 H -13.8" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">5</text></g><g class="a9 c2 node n8 p0" transform="translate(-13.8 0.8617)"><g class="a8 i1 leaf node n3 p0 sample" transform="translate(-9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">3</text></g><g class="a8 i3 leaf node n6 p0 sample" transform="translate(9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">6</text></g><path class="edge" d="M 0 0 V -0.8617 H 13.8" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">8</text></g><path class="edge" d="M 0 0 V -1.39353 H 16.1" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">9</text></g><path class="edge" d="M 0 0 V -5.24241 H -28.75" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">10</text></g><path class="edge" d="M 0 0 V -2.03424 H -23.575" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">12</text></g><path class="edge" d="M 0 0 V -30.7712 H -25.5875" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">13</text></g><g class="a18 c2 node n14 p0" transform="translate(-25.5875 21.7136)"><g class="a14 i0 leaf node n0 p0 sample" transform="translate(-9.2 22.1457)"><path class="edge" d="M 0 0 V -22.1457 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">0</text></g><g class="a14 i3 leaf node n7 p0 sample" transform="translate(9.2 22.1457)"><path class="edge" d="M 0 0 V -22.1457 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">7</text></g><path class="edge" d="M 0 0 V -21.7136 H 25.5875" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">14</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab" transform="translate(0 -11)">18</text></g></g></g><g class="tree t1" transform="translate(244 0)"><g class="plotbox"><g class="c2 node n18 p0 root" transform="translate(63.9875 77.5407)"><g class="a18 c2 node n13 p0" transform="translate(25.5875 30.7712)"><g class="a13 i0 leaf node n1 p0 sample" transform="translate(-23.575 13.0881)"><path class="edge" d="M 0 0 V -13.0881 H 23.575" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">1</text></g><g class="a13 c2 node n12 p0" transform="translate(23.575 2.03424)"><g class="a12 i1 leaf m3 node n2 p0 s3 sample" transform="translate(-28.75 11.0538)"><path class="edge" d="M 0 0 V -11.0538 H 28.75" /><g class="mut m3 s3" transform="translate(0 -3.79834)"><line x1="0" x2="0" y1="0" y2="3.79834" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab lft" transform="translate(-5 0)">3</text></g><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">2</text></g><g class="a12 c2 node n10 p0" transform="translate(28.75 5.24241)"><g class="a10 i2 leaf node n4 p0 sample" transform="translate(16.1 5.81144)"><path class="edge" d="M 0 0 V -5.81144 H -16.1" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">4</text></g><g class="a10 c2 node n9 p0" transform="translate(-16.1 1.39353)"><g class="a9 i2 leaf node n5 p0 sample" transform="translate(13.8 4.41791)"><path class="edge" d="M 0 0 V -4.41791 H -13.8" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">5</text></g><g class="a9 c2 node n8 p0" transform="translate(-13.8 0.8617)"><g class="a8 i1 leaf m0 node n3 p0 s0 sample" transform="translate(-9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H 9.2" /><g class="mut m0 s0" transform="translate(0 -2.73828)"><line x1="0" x2="0" y1="0" y2="2.73828" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab lft" transform="translate(-5 0)">0</text></g><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">3</text></g><g class="a8 i3 leaf m1 node n6 p0 s1 sample" transform="translate(9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H -9.2" /><g class="mut m1 s1" transform="translate(0 -1.38406)"><line x1="0" x2="0" y1="0" y2="1.38406" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab rgt" transform="translate(5 0)">1</text></g><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">6</text></g><path class="edge" d="M 0 0 V -0.8617 H 13.8" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">8</text></g><path class="edge" d="M 0 0 V -1.39353 H 16.1" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">9</text></g><path class="edge" d="M 0 0 V -5.24241 H -28.75" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">10</text></g><path class="edge" d="M 0 0 V -2.03424 H -23.575" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">12</text></g><path class="edge" d="M 0 0 V -30.7712 H -25.5875" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">13</text></g><g class="a18 c2 node n16 p0" transform="translate(-25.5875 6.71959)"><g class="a16 i0 leaf m2 node n0 p0 s2 sample" transform="translate(-9.2 37.1397)"><path class="edge" d="M 0 0 V -37.1397 H 9.2" /><g class="mut m2 s2" transform="translate(0 -29.5048)"><line x1="0" x2="0" y1="0" y2="29.5048" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab lft" transform="translate(-5 0)">2</text></g><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">0</text></g><g class="a16 i3 leaf node n7 p0 sample" transform="translate(9.2 37.1397)"><path class="edge" d="M 0 0 V -37.1397 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">7</text></g><path class="edge" d="M 0 0 V -6.71959 H 25.5875" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">16</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab" transform="translate(0 -11)">18</text></g></g></g><g class="tree t2" transform="translate(431.2 0)"><g class="plotbox"><g class="c2 node n18 p0 root" transform="translate(73.475 77.5407)"><g class="a18 c2 m5 node n13 p0 s5" transform="translate(30.475 30.7712)"><g class="a13 i0 leaf node n1 p0 sample" transform="translate(-19.55 13.0881)"><path class="edge" d="M 0 0 V -13.0881 H 19.55" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">1</text></g><g class="a13 c2 m4 node n12 p0 s4" transform="translate(19.55 2.03424)"><g class="a12 i1 leaf node n2 p0 sample" transform="translate(-20.7 11.0538)"><path class="edge" d="M 0 0 V -11.0538 H 20.7" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">2</text></g><g class="a12 c2 node n9 p0" transform="translate(20.7 6.63594)"><g class="a9 i2 leaf node n5 p0 sample" transform="translate(13.8 4.41791)"><path class="edge" d="M 0 0 V -4.41791 H -13.8" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">5</text></g><g class="a9 c2 node n8 p0" transform="translate(-13.8 0.8617)"><g class="a8 i1 leaf node n3 p0 sample" transform="translate(-9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">3</text></g><g class="a8 i3 leaf node n6 p0 sample" transform="translate(9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">6</text></g><path class="edge" d="M 0 0 V -0.8617 H 13.8" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">8</text></g><path class="edge" d="M 0 0 V -6.63594 H -20.7" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">9</text></g><path class="edge" d="M 0 0 V -2.03424 H -19.55" /><g class="mut m4 s4" transform="translate(0 -1.54933)"><line x1="0" x2="0" y1="0" y2="1.54933" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab rgt" transform="translate(5 0)">4</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">12</text></g><path class="edge" d="M 0 0 V -30.7712 H -30.475" /><g class="mut m5 s5" transform="translate(0 -28.0667)"><line x1="0" x2="0" y1="0" y2="28.0667" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab rgt" transform="translate(5 0)">5</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">13</text></g><g class="a18 c2 node n16 p0" transform="translate(-30.475 6.71959)"><g class="a16 i0 leaf node n0 p0 sample" transform="translate(-13.8 37.1397)"><path class="edge" d="M 0 0 V -37.1397 H 13.8" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">0</text></g><g class="a16 c2 node n11 p0" transform="translate(13.8 29.6414)"><g class="a11 i2 leaf node n4 p0 sample" transform="translate(-9.2 7.49833)"><path class="edge" d="M 0 0 V -7.49833 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">4</text></g><g class="a11 i3 leaf node n7 p0 sample" transform="translate(9.2 7.49833)"><path class="edge" d="M 0 0 V -7.49833 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">7</text></g><path class="edge" d="M 0 0 V -29.6414 H -13.8" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">11</text></g><path class="edge" d="M 0 0 V -6.71959 H 30.475" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">16</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab" transform="translate(0 -11)">18</text></g></g></g><g class="tree t3" transform="translate(618.4 0)"><g class="plotbox"><g class="c2 node n19 p0 root" transform="translate(68.5875 26.8)"><g class="a19 i0 leaf m8 node n0 p0 s8 sample" transform="translate(-39.3875 94.6)"><path class="edge" d="M 0 0 V -94.6 H 39.3875" /><g class="mut m8 s8" transform="translate(0 -30.0795)"><line x1="0" x2="0" y1="0" y2="30.0795" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab lft" transform="translate(-5 0)">8</text></g><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">0</text></g><g class="a19 c2 m6 m7 node n18 p0 s6 s7" transform="translate(39.3875 50.7407)"><g class="a18 c2 node n13 p0" transform="translate(-40.825 30.7712)"><g class="a13 i0 leaf node n1 p0 sample" transform="translate(-19.55 13.0881)"><path class="edge" d="M 0 0 V -13.0881 H 19.55" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">1</text></g><g class="a13 c2 node n12 p0" transform="translate(19.55 2.03424)"><g class="a12 i1 leaf node n2 p0 sample" transform="translate(-20.7 11.0538)"><path class="edge" d="M 0 0 V -11.0538 H 20.7" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">2</text></g><g class="a12 c2 node n9 p0" transform="translate(20.7 6.63594)"><g class="a9 i2 leaf node n5 p0 sample" transform="translate(13.8 4.41791)"><path class="edge" d="M 0 0 V -4.41791 H -13.8" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">5</text></g><g class="a9 c2 node n8 p0" transform="translate(-13.8 0.8617)"><g class="a8 i1 leaf node n3 p0 sample" transform="translate(-9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">3</text></g><g class="a8 i3 leaf node n6 p0 sample" transform="translate(9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">6</text></g><path class="edge" d="M 0 0 V -0.8617 H 13.8" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">8</text></g><path class="edge" d="M 0 0 V -6.63594 H -20.7" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">9</text></g><path class="edge" d="M 0 0 V -2.03424 H -19.55" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">12</text></g><path class="edge" d="M 0 0 V -30.7712 H 40.825" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">13</text></g><g class="a18 c2 m9 node n11 p0 s9" transform="translate(40.825 36.361)"><g class="a11 i2 leaf node n4 p0 sample" transform="translate(-9.2 7.49833)"><path class="edge" d="M 0 0 V -7.49833 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">4</text></g><g class="a11 i3 leaf node n7 p0 sample" transform="translate(9.2 7.49833)"><path class="edge" d="M 0 0 V -7.49833 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">7</text></g><path class="edge" d="M 0 0 V -36.361 H -40.825" /><g class="mut m9 s9" transform="translate(0 -5.53803)"><line x1="0" x2="0" y1="0" y2="5.53803" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab rgt" transform="translate(5 0)">9</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">11</text></g><path class="edge" d="M 0 0 V -50.7407 H -39.3875" /><g class="mut m7 s7" transform="translate(0 -45.4606)"><line x1="0" x2="0" y1="0" y2="45.4606" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab rgt" transform="translate(5 0)">7</text></g><g class="mut m6 s6" transform="translate(0 -7.38293)"><line x1="0" x2="0" y1="0" y2="7.38293" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab rgt" transform="translate(5 0)">6</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">18</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab" transform="translate(0 -11)">19</text></g></g></g><g class="tree t4" transform="translate(805.6 0)"><g class="plotbox"><g class="c2 node n19 p0 root" transform="translate(98.4875 26.8)"><g class="a19 c2 node n11 p0" transform="translate(50.3125 87.1017)"><g class="a11 i2 leaf m12 node n4 p0 s12 sample" transform="translate(-9.2 7.49833)"><path class="edge" d="M 0 0 V -7.49833 H 9.2" /><g class="mut m12 s12" transform="translate(0 -5.04465)"><line x1="0" x2="0" y1="0" y2="5.04465" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab lft" transform="translate(-5 0)">12</text></g><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">4</text></g><g class="a11 i3 leaf node n7 p0 sample" transform="translate(9.2 7.49833)"><path class="edge" d="M 0 0 V -7.49833 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">7</text></g><path class="edge" d="M 0 0 V -87.1017 H -50.3125" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">11</text></g><g class="a19 c2 node n17 p0" transform="translate(-50.3125 51.63)"><g class="a17 i0 leaf m10 node n0 p0 s10 sample" transform="translate(-18.975 42.97)"><path class="edge" d="M 0 0 V -42.97 H 18.975" /><g class="mut m10 s10" transform="translate(0 -10.9638)"><line x1="0" x2="0" y1="0" y2="10.9638" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab lft" transform="translate(-5 0)">10</text></g><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">0</text></g><g class="a17 c2 m11 node n13 p0 s11" transform="translate(18.975 29.8819)"><g class="a13 i0 leaf node n1 p0 sample" transform="translate(-19.55 13.0881)"><path class="edge" d="M 0 0 V -13.0881 H 19.55" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">1</text></g><g class="a13 c2 node n12 p0" transform="translate(19.55 2.03424)"><g class="a12 i1 leaf node n2 p0 sample" transform="translate(-20.7 11.0538)"><path class="edge" d="M 0 0 V -11.0538 H 20.7" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">2</text></g><g class="a12 c2 node n9 p0" transform="translate(20.7 6.63594)"><g class="a9 i2 leaf node n5 p0 sample" transform="translate(13.8 4.41791)"><path class="edge" d="M 0 0 V -4.41791 H -13.8" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">5</text></g><g class="a9 c2 node n8 p0" transform="translate(-13.8 0.8617)"><g class="a8 i1 leaf node n3 p0 sample" transform="translate(-9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">3</text></g><g class="a8 i3 leaf node n6 p0 sample" transform="translate(9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">6</text></g><path class="edge" d="M 0 0 V -0.8617 H 13.8" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">8</text></g><path class="edge" d="M 0 0 V -6.63594 H -20.7" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">9</text></g><path class="edge" d="M 0 0 V -2.03424 H -19.55" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">12</text></g><path class="edge" d="M 0 0 V -29.8819 H -18.975" /><g class="mut m11 s11" transform="translate(0 -24.632)"><line x1="0" x2="0" y1="0" y2="24.632" /><path class="sym" d="M -3,-3 l 6,6 M -3,3 l 6,-6" /><text class="lab rgt" transform="translate(5 0)">11</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">13</text></g><path class="edge" d="M 0 0 V -51.63 H 50.3125" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">17</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab" transform="translate(0 -11)">19</text></g></g></g><g class="tree t5" transform="translate(992.8 0)"><g class="plotbox"><g class="c2 node n19 p0 root" transform="translate(95.325 26.8)"><g class="a19 c2 node n17 p0" transform="translate(-48.875 51.63)"><g class="a17 i0 leaf node n0 p0 sample" transform="translate(-17.25 42.97)"><path class="edge" d="M 0 0 V -42.97 H 17.25" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">0</text></g><g class="a17 c2 node n13 p0" transform="translate(17.25 29.8819)"><g class="a13 i0 leaf node n1 p0 sample" transform="translate(-16.1 13.0881)"><path class="edge" d="M 0 0 V -13.0881 H 16.1" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">1</text></g><g class="a13 c2 node n12 p0" transform="translate(16.1 2.03424)"><g class="a12 i1 leaf node n2 p0 sample" transform="translate(-13.8 11.0538)"><path class="edge" d="M 0 0 V -11.0538 H 13.8" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">2</text></g><g class="a12 c2 node n8 p0" transform="translate(13.8 7.49764)"><g class="a8 i1 leaf node n3 p0 sample" transform="translate(-9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">3</text></g><g class="a8 i3 leaf node n6 p0 sample" transform="translate(9.2 3.55621)"><path class="edge" d="M 0 0 V -3.55621 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">6</text></g><path class="edge" d="M 0 0 V -7.49764 H -13.8" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">8</text></g><path class="edge" d="M 0 0 V -2.03424 H -16.1" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">12</text></g><path class="edge" d="M 0 0 V -29.8819 H -17.25" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">13</text></g><path class="edge" d="M 0 0 V -51.63 H 48.875" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">17</text></g><g class="a19 c2 node n15 p0" transform="translate(48.875 57.9197)"><g class="a15 i2 leaf node n5 p0 sample" transform="translate(13.8 36.6803)"><path class="edge" d="M 0 0 V -36.6803 H -13.8" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">5</text></g><g class="a15 c2 node n11 p0" transform="translate(-13.8 29.182)"><g class="a11 i2 leaf node n4 p0 sample" transform="translate(-9.2 7.49833)"><path class="edge" d="M 0 0 V -7.49833 H 9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">4</text></g><g class="a11 i3 leaf node n7 p0 sample" transform="translate(9.2 7.49833)"><path class="edge" d="M 0 0 V -7.49833 H -9.2" /><rect class="sym" height="6" width="6" x="-3" y="-3" /><text class="lab" transform="translate(0 11)">7</text></g><path class="edge" d="M 0 0 V -29.182 H 13.8" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab lft" transform="translate(-3 -7.0)">11</text></g><path class="edge" d="M 0 0 V -57.9197 H -48.875" /><circle class="sym" cx="0" cy="0" r="3" /><text class="lab rgt" transform="translate(3 -7.0)">15</text></g><circle class="sym" cx="0" cy="0" r="3" /><text class="lab" transform="translate(0 -11)">19</text></g></g></g></g></g></svg>'

Interaction with R libraries#

R has a number of libraries to deal with genomic data and trees. Below we focus on the phylogenetic tree representation defined in the the popular ape package, taking all the trees exported in Nexus format, or individual trees exported in Newick format:

file = tempfile()
ts_mut$write_nexus(file)
# Warning - ape trees are stored independently, so this will use much more memory than tskit
trees <- ape::read.nexus(file, force.multi = TRUE)  # return a set of trees

# Or simply read in a single tree
tree <- ape::read.tree(text=ts_mut$first()$as_newick())

# Now we can plot the tree in tskit style, but using the ape library
plot(tree, direction="downward", srt=90, adj=0.5)  # or equivalently use trees[[1]]
_images/ba4e21555374a0a3002f3cbfaba7386fcc049d206857d1b97cda5bae4dc2c46e.png

Note that nodes are labelled with the prefix n, so that nodes 0, 1, 2, … become n0, n1, n2 … etc. This helps to avoid confusion between the the zero-based counting system used natively by tskit, and the one-based counting system used in R.

Further information#

Be sure to check out the reticulate documentation, in particular on Calling Python from R, which includes important information on how R data types are converted to their equivalent Python types.

Contents