*Analysing tree sequences*#

Note

This tutorial is a work in progress.

## Computing statistics#

Tskit provides an extensive and flexible interface for computing population genetic statistics, which is documented in detail in the general statistics section of the offical documentation. This tutorial aims to give a quick overview of how the APIs work and how to use them effectively.

First, let’s load a tree sequence to work with which has roughly human parameters for 10 thousand samples and 10Mb chromosomes:

```
ts = tskit.load("data/computing_statistics.trees")
ts
```

Tree Sequence | |
---|---|

Trees | 36599 |

Sequence Length | 10000000.0 |

Time Units | unknown |

Sample Nodes | 10000 |

Total Size | 9.5 MiB |

Metadata | No Metadata |

Table | Rows | Size | Has Metadata |
---|---|---|---|

Edges | 157894 | 4.8 MiB | |

Individuals | 0 | 24 Bytes | |

Migrations | 0 | 8 Bytes | |

Mutations | 38778 | 1.4 MiB | |

Nodes | 44864 | 1.2 MiB | |

Populations | 1 | 16 Bytes | |

Provenances | 1 | 1.0 KiB | |

Sites | 38778 | 946.7 KiB |

This tree sequence has ~36.6 thousand trees & ~39 thousand segregating sites. We’d now like to compute some statistics on this dataset.

### One-way statistics#

We refer to statistics that are defined with respect to a single set of
samples as “one-way”. An example of such a statistic is diversity, which
is computed using the `TreeSequence.diversity()`

method:

```
d = ts.diversity()
print("Average diversity per unit sequence length = {d:.3G}")
```

```
Average diversity per unit sequence length = {d:.3G}
```

This tells the average diversity across the whole sequence and returns a single
number. We’ll usually want to compute statistics in
windows along the genome and we
use the `windows`

argument to do this:

```
windows = np.linspace(0, ts.sequence_length, num=5)
d = ts.diversity(windows=windows)
print(windows, d, sep="\n")
```

```
[ 0. 2500000. 5000000. 7500000. 10000000.]
[0.00040717 0.00040891 0.00041644 0.00039899]
```

The `windows`

argument takes a numpy array specifying the breakpoints
along the genome. Here, we use numpy to create four equally spaced windows
of size 2.5 megabases (the windows array contains k + 1 elements to define
k windows). Because we have asked for values in windows, tskit now returns
a numpy array rather than a single value. (See
Output dimensions for a full description of how the output
dimensions of statistics are determined by the `windows`

argument.)

Suppose we wanted to compute diversity within a specific subset of samples.
We can do this using the `sample_sets`

argument:

```
A = ts.samples()[:100]
d = ts.diversity(sample_sets=A)
print(d)
```

```
0.00041007806060604034
```

Here, we’ve computed the average diversity within the first hundred samples across the whole genome. As we’ve not specified any windows, this is again a single value.

We can also compute diversity in *multiple* sample sets at the same time by providing
a list of sample sets as an argument:

```
A = ts.samples()[:100]
B = ts.samples()[100:200]
C = ts.samples()[200:300]
d = ts.diversity(sample_sets=[A, B, C])
print(d)
```

```
[0.00041008 0.00040893 0.00040853]
```

Because we’ve computed multiple statistics concurrently, tskit returns a numpy array of these statistics. We have asked for diversity within three different sample sets, and tskit therefore returns an array with three values. (In general, the dimensions of the input determines the dimensions of the output: see Output dimensions for a detailed description of the rules.)

We can also compute multiple statistics in multiple windows:

```
d = ts.diversity(sample_sets=[A, B, C], windows=windows)
print("shape = ", d.shape, "\n", d)
```

```
shape = (4, 3)
[[0.00041199 0.00040542 0.00040588]
[0.00040889 0.00041216 0.00041419]
[0.00041636 0.00041463 0.0004128 ]
[0.00040307 0.0004035 0.00040126]]
```

We have computed diversity within three different sample sets across four genomic windows, and our output is therefore a 2D numpy array with four rows and three columns: each row contains the diversity values within A, B and C for a particular window.

### Multi-way statistics#

Many population genetic statistics compare multiple sets of samples to
each other. For example, the `TreeSequence.divergence()`

method computes
the divergence between two subsets of samples:

```
A = ts.samples()[:100]
B = ts.samples()[:100]
d = ts.divergence([A, B])
print(d)
```

```
0.0004059772800000012
```

The divergence between two sets of samples A and B is a single number, and we we again return a single floating point value as the result. We can also compute this in windows along the genome, as before:

```
d = ts.divergence([A, B], windows=windows)
print(d)
```

```
[0.00040787 0.00040481 0.00041219 0.00039904]
```

Again, as we have defined four genomic windows along the sequence, the result is numpy array with four values.

A powerful feature of tskit’s stats API is that we can compute the divergences
between multiple sets of samples simultaneously using the `indexes`

argument:

```
d = ts.divergence([A, B, C], indexes=[(0, 1), (0, 2)])
print(d)
```

```
[0.00040598 0.00040943]
```

Here, we’ve specified three sample sets A, B and C and we’ve computed the
divergences between A and B, and between A and C. The `indexes`

argument is used
to specify which pairs of sets we are interested in. In this example
we’ve computed two different divergence values and the output is therefore
a numpy array of length 2.

As before, we can combine computing multiple statistics in multiple windows to return a 2D numpy array:

```
windows = np.linspace(0, ts.sequence_length, num=5)
d = ts.divergence([A, B, C], indexes=[(0, 1), (0, 2)], windows=windows)
print(d)
```

```
[[0.00040787 0.00040834]
[0.00040481 0.00041207]
[0.00041219 0.00041464]
[0.00039904 0.00040267]]
```

Each row again corresponds to a window, which contains the average divergence values between the chosen sets.

If the `indexes`

parameter is 1D array, we interpret this as specifying
a single statistic and remove the empty outer dimension:

```
d = ts.divergence([A, B, C], indexes=(0, 1))
print(d)
```

```
0.0004059772800000012
```

It’s important to note that we don’t **have** to remove empty dimensions: tskit
will only do this if you explicitly ask it to. Here, for example, we can keep the
output as an array with one value if we wish:

```
d = ts.divergence([A, B, C], indexes=[(0, 1)])
print(d)
```

Please see Sample sets and indexes for a
full description of the `sample_sets`

and `indexes`

arguments.

## Allele frequency spectra#

The allele frequency spectrum is a fundamental tool in population genetics, and tskit provides a flexible and powerful approach to computing such spectra. Suppose we have simulated the following tree sequence:

```
from IPython.display import display
ts = tskit.load("data/afs.trees")
tree = ts.first()
display(tree.draw_svg())
ts.tables.sites
```

id | position | ancestral_state | metadata |
---|---|---|---|

0 | 0.30043643 | 0 | |

1 | 0.32220794 | 0 | |

2 | 0.36507027 | 0 | |

3 | 0.50940255 | 0 | |

4 | 0.51327137 | 0 | |

5 | 0.51400861 | 0 | |

6 | 0.54796110 | 0 | |

7 | 0.75929404 | 0 | |

8 | 0.80591800 | 0 | |

9 | 0.92324208 | 0 |

Computing the allele frequency spectrum is then easy:

```
afs = ts.allele_frequency_spectrum(polarised=True, span_normalise=False)
print(afs)
```

```
[0. 2. 6. 1. 1. 0. 0.]
```

This tells us that we have two singletons, six doubletons and one 3-ton and
one 4-ton. Note that the first element of the returned AFS array does *not* correspond
to the singletons (see below for why). Because we have simulated these mutations,
we know the ancestral and derived states we have set `polarised`

to True. We
can get the “folded” AFS by setting polarised to False. Because we want simple
counts here and not averaged values, we set `span_normalise=False`

: by
default, windowed statistics are divided by the sequence length, so they are
comparable between windows.

The returned value here is actually a 2D array, and this is because we can also perform these computations in windows along the genome:

```
afs = ts.allele_frequency_spectrum(windows=[0, 0.5, 1], span_normalise=False, polarised=True)
print(afs)
```

```
[[0. 1. 1. 1. 0. 0. 0.]
[0. 1. 5. 0. 1. 0. 0.]]
```

This time, we’ve asked for the number of sites at each frequency in two equal windows. Now we can see that in the first half of the sequence we have three sites (compare with the site table above): one singleton, one doubleton and one tripleton.

### Joint spectra#

We can also compute allele frequencies within multiple sets of samples,
the *joint allele frequency spectra*.

```
node_colours = {0: "blue", 2: "blue", 3: "blue", 1: "green", 4: "green", 5: "green"}
styles = [f".n{k} > .sym {{fill: {v}}}" for k, v in node_colours.items()]
tree.draw_svg(style = "".join(styles))
```

Here we’ve marked the samples as either blue or green (we can imagine these belonging to different populations, for example). We can then compute the joint AFS based on these two sets:

```
afs = ts.allele_frequency_spectrum([[0, 2, 3], [1, 4, 5]], polarised=True)
print(afs)
```

```
[[0. 2. 0. 0.]
[0. 6. 0. 0.]
[0. 1. 1. 0.]
[0. 0. 0. 0.]]
```

Now, each window in our AFS is a 2D numpy array, where each dimension corresponds to frequencies within the different sets. So, we see for example that there are six sites that are singletons in both sets, 1 site that is a doubleton in both sets, and 2 sites that singletons in \([1, 4, 5]\) and not present in the other sample set.

### Branch length spectra#

Up to now we’ve used the `allele_frequency_spectrum()`

method
to summarise the number of sites that occur at different frequencies. We can also
use this approach to compute the total branch lengths subtending a given
number of samples by setting `mode="branch"`

:

```
afs = ts.allele_frequency_spectrum(mode="branch", polarised=True, span_normalise=False)
print(afs)
```

```
[0. 4.86089166 5.39638988 2.55239269 2.07444286 0.
0. ]
```

Thus, the total branch length over example one sample is 4.86, over two is 5.39, and so on.

### Zeroth and final entries in the AFS#

The zeroth element of the AFS is significant when we are working with sample sets that are a subset of all samples in the tree sequence. For example, in the following we compute the AFS within the sample set [0, 1, 2]:

```
afs = ts.allele_frequency_spectrum([[0, 1, 2]], mode="branch", polarised=True)
print(afs)
```

```
[4.33184862 5.30022646 5.252042 0. ]
```

Thus, the total branch length over 0, 1 and 2 is 5.3, and over pairs from this set
is 5.25. What does the zeroth value of 4.33 signify? This is the total branch length
over all samples that are **not** in this sample set. By including this value, we
maintain the property that for each tree, the sum of the AFS for any sample set
is always equal to the total branch length. For example, here we compute:

```
print("sum afs = ", np.sum(afs))
print("total branch len = ", tree.total_branch_length)
```

```
sum afs = 14.884117086717392
total branch len = 14.884117086717396
```

The final entry of the AFS is similar: it counts alleles (for mode=”site”) or branches (for mode=”branch”) that are ancestral to all of the given sample set, but are still polymorphic in the entire set of samples of the tree sequence. Note, however, that alleles fixed among all the samples, e.g., ones above the root of the tree, will not be included.