slides_slendr.qmd

---
format:                         ### slides
  revealjs:                     ### slides
    echo: true                  ### slides
    code-line-numbers: false    ### slides
    fig-align: center           ### slides
    slide-number: true          ### slides
    self-contained: true        ### slides
---

```{r}
#| echo: false
library(slendr)
init_env(quiet = TRUE)
```


# Simulations with _slendr_


# Why use simulations?

# 

::: {.columns} 
::: {.column width="70%"}
> Many problems in population genetics cannot be solved by a mathematician,
no matter how gifted. \[It\] is already clear that computer methods are very
powerful. This is good. It \[...\] permits people with limited mathematical
knowledge to work on important problems \[...\].
:::

::: {.column width="30%"}
![](files/slendr/images/crow.jpeg)

[James F. Crow](https://en.wikipedia.org/wiki/James_F._Crow) -- [interview](http://www.gnxp.com/blog/2006/06/10-questions-for-jim-crow.php)
:::
:::










# Why use simulations?

1. Developing intuition into statistics
2. Estimating model parameters
3. Ground truth for method development










## Developing intuition into statistics

<center>![](files/slendr/images/fstats_sims.png)</center>

::: {.aside}
<small>Image from [Peter (2016)](https://academic.oup.com/genetics/article/202/4/1485/5930214)</small>
:::










## Developing intuition into statistics

<center>![](files/slendr/images/admixture.png)</center>

::: {.aside}
<small>Image from [Lawson _et al._ (2018)](https://www.nature.com/articles/s41467-018-05257-7)</small>
:::










## Estimating model parameters (i.e. [ABC](https://en.wikipedia.org/wiki/Approximate_Bayesian_computation))

<center>![](files/slendr/images/abc_scheme.png){width="50%"}</center>

::: {.aside}
<small>Image from [Wikipedia on ABC](https://en.wikipedia.org/wiki/Approximate_Bayesian_computation)</small>
:::










## Ground truth for method development

<center>![](files/slendr/images/mcmc.png)</center>

::: {.aside}
<small>Image from [Schiffels and Durbin (2014)](https://www.nature.com/articles/ng.3015)</small>
:::








## Simulation software

The most famous and widely used are [SLiM](https://messerlab.org/slim/) and [_msprime_](https://tskit.dev/msprime/docs/stable/intro.html).

::: {.fragment}

They are <u>very</u> powerful and (nearly) infinitely flexible.

:::

::: {.fragment}

However, they both require:

- quite a bit of code for complex simulations ("complex" is relative, of course)
- relatively high confidence in programming

<!-- (🐛🪲🐜). -->
:::

::: {.fragment}
<center>

<h4>Our exercises will focus on the [_slendr_](http://www.slendr.net)<br> simulation toolkit for population genetics in R.</h4>

</center>
:::

::: {.fragment}
<center><br>But, as a recap, let's look at _msprime_ and SLiM a little bit...</center>
:::










## 

::: {.columns} 

::: {.column width="60%"}
<h2>What is _msprime_?</h2>

- A Python module for writing **coalescent simulations**

-   Extremely fast (genome-scale, population-scale data!)

-   You should know Python fairly well to build complex models
:::

::: {.column width="40%"}
<center>

![<small>Image modified from [Alexei Drummond</small>](http://alexeidrummond.org/bayesian_phylo_lectures/lecture10/)](files/slendr/images/sim_sketches.001.png){width="100%"}


</center>
:::

:::










## Simple simulation using _msprime_

```{bash}
#| echo: false
cat files/slendr/script.py
```






## 

::: {.columns} 
::: {.column width="60%"}
<h2>What is SLiM?</h2>

- **A forward-time simulator**

- Has its own programming language

- Massive library of functions for:
  - demographic events
  - various mating systems
  - natural selection

-   More than 700 pages long [manual](https://github.com/MesserLab/SLiM/releases/download/v4.3/SLiM_Manual.pdf)!
:::

::: {.column width="40%"}
<center>

![<small>Image modified from [Alexei Drummond</small>](http://alexeidrummond.org/bayesian_phylo_lectures/lecture10/)](files/slendr/images/sim_sketches.001.png){width="100%"}

</center>
:::
:::










## Simple neutral simulation in SLiM

```{bash}
#| echo: false
cat files/slendr/script.slim
```















#

<h3>SLiM and _msprime_ are both incredible pieces of software...</h3>





## ... so why _slendr_?

<center>

![](files/slendr/images/slendr_logo.png){width="30%"}

<h4 style="font-family: 'Courier New'">[www.slendr.net](https://www.slendr.net)</h4>

</center>

<small><br></small>










## First motivation: spatial simulations!

<center>![](files/slendr/images/animation.gif){width="70%"}</center>










## A broader motivation for _slendr_

- Most researchers are not expert programmers

- All but the most trivial simulations require lots of code

::: {.fragment}
- Yet, 90% <sup><font color="blue">\[citation needed\]</font></sup> of simulations are basically the same!

  - create populations (splits and $N_e$ changes)

  - specify admixture rates and admixture times

- ... all this means duplication of code across many projects
:::

::: {.fragment}
- Computing statistics presents even more hurdles
:::

::: {.fragment}
<center><h4>_slendr_ makes this <u>very easy</u>, even for "complex models"</h4></center>
:::

# _slendr_ crash course

<br> <!-- ### remove for slides  -->
You can use this document as a cheat sheet as you work on the exercises. <!-- ### remove for slides  -->
It contains everything you need to know about _slendr_, in a compressed form. <!-- ### remove for slides  -->

## Typical _slendr_ workflow

We will always start our R scripts with this:

```{r}
#| output: false
library(slendr) # You can safely ignore any potential warnings!
init_env()      # This activates the internal Python environmet
```

<br>Followed by some combination of the following:

1. creating populations
2. programming $N_e$ size changes
3. encoding gene-flow events
4. simulating genomic data
5. computing popgen statistics










## Creating populations

At minimum, we need its name, size and "time of appearance":

```{r}
pop1 <- population("pop1", N = 1000, time = 1)
```

::: {.fragment}

This creates a normal R object! Typing it out gives a summary:

```{r}
pop1
```

:::

::: {.aside}
**Note:** Because _slendr_ uses either _msprime_ or SLiM internally for simulation
of genomic data, all individuals are assumed to be diploid.
:::









## Programming population splits

Splits are defined by providing a `parent = <pop>` argument:

```{r}
pop2 <- population("pop2", N = 100, time = 50, parent = pop1)
```

::: {.aside}
**Note:** Here `pop1` is an R object created above, not a string `"pop1"`!
:::

::: {.fragment}

The split is again reported in the "historical summary":

```{r}
pop2
```

:::









## Scheduling resize events

- Step size decrease:

```{r}
#| code-line-numbers: "1|2"
pop1 <- population("pop1", N = 1000, time = 1)
pop1_step <- resize(pop1, N = 100, time = 500, how = "step")
```

::: {.fragment}
- Exponential increase:

```{r}
#| code-line-numbers: true
pop2 <- population("pop2", N = 100, time = 50, parent = pop1)
pop2_exp <- resize(pop2, N = 10000, time = 500, end = 2000, how = "exponential")
```
:::










## Tidyverse-style [pipe](https://magrittr.tidyverse.org) `%>%` interface

The following leads to a more concise (and "elegant") code.

- Step size decrease:

```{r}
pop1 <-
  population("pop1", N = 1000, time = 1) %>%
  resize(N = 100, time = 500, how = "step")
```

- Exponential increase:

```{r}
pop2 <-
  population("pop2", N = 1000, time = 1) %>%
  resize(N = 10000, time = 500, end = 2000, how = "exponential")
```


::: {.aside}
**Note:** You can read (and understand) `a() %>% b() %>% c()` as "take the
result of the function `a`, pipe it into function `b`, and then pipe _that_ to
function `c`".
:::







## A more complex model

Using just the two functions introduced so far:

```{r}
#| code-line-numbers: "|15-18"
pop1 <- population("pop1", N = 1000, time = 1)

pop2 <-
  population("pop2", N = 1000, time = 300, parent = pop1) %>%
  resize(N = 100, how = "step", time = 1000)

pop3 <-
  population("pop3", N = 1000, time = 400, parent = pop2) %>%
  resize(N = 2500, how = "step", time = 800)

pop4 <-
  population("pop4", N = 1500, time = 500, parent = pop3) %>%
  resize(N = 700, how = "exponential", time = 1200, end = 2000)

pop5 <-
  population("pop5", N = 100, time = 600, parent = pop4) %>%
  resize(N = 50, how = "step", time = 900) %>%
  resize(N = 1000, how = "exponential", time = 1600, end = 2200)
```











## Again, each object carries its history!

For instance, this is the summary you will get from the last population from
the previous code chunk:

```{r}
pop5
```

::: {.fragment}

**This way, you can build up complex models step by step, checking things
as you go by interacting with the R console.**

:::










## Gene flow / admixture

We can schedule gene flow from `pop1` into `pop2` with:

```{r}
gf <- gene_flow(from = pop1, to = pop2, start = 2000, end = 2200, rate = 0.13)
```

::: {.aside}
**Note:** Here `rate = 0.13` means 13% migrants over the given time window will come from "pop1" into "pop2".
:::

::: {.fragment}

<br>

Multiple gene-flow events can be gathered in a list:

```{r}
#| eval: false
gf <- list(
  gene_flow(from = pop1, to = pop2, start = 500, end = 600, rate = 0.13),
  gene_flow(from = ..., to = ..., start = ..., end = ..., rate = ...),
  ...
)
```

:::








## Model compilation

<br>

This is the final step before we can simulate data.

<br>

```{r}
#| code-line-numbers: true
model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  generation_time = 1,       # (converts the all times into generations)
  simulation_length = 3000,  # (number of generations to run the simulation for)
  direction = "forward"      # (not totally necessary but good practice)
)
```

<br>

<center>**`compile_model()` takes a list of components, performs some consistency checks, and returns a single R object**</center>










## Model compilation

<br>

This is the final step before we can simulate data.

<br>

```{r}
#| code-line-numbers: "3"
model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  gene_flow = gf,      # <----- in case our model includes gene flow(s)
  generation_time = 1,
  simulation_length = 3000,
  direction = "forward"
)
```

<br>

<center>**Gene flow(s) can be included via the `gene_flow` argument.**</center>










## Model summary

Typing the compiled `model` into R prints a brief summary:

```{r}
model
```

This can be useful as a quick overview of the model we are working with. 
However, **a better way to check a model is...**












## Model visualization

```{r}
#| eval: false
plot_model(model)
```

```{r}
#| fig-width: 7
#| fig-align: center
#| echo: false
plot_model(model, proportions = TRUE)
```










# A note on units of time (and its direction)

<br>

Sometimes you want to work with time units
such as "years ago" (aDNA radio-carbon dated samples, etc.),
but you have to convert those to "generations forward" for some software.

_slendr_ helps by making it possible to use whatever time units
or time directions are more natural for a particular project.






## "Forward time units"

```{r}
#| code-line-numbers: "|1|22-23"
pop1 <- population("pop1", N = 1000, time = 1)

pop2 <-
  population("pop2", N = 1000, time = 300, parent = pop1) %>%
  resize(N = 100, how = "step", time = 1000)

pop3 <-
  population("pop3", N = 1000, time = 400, parent = pop2) %>%
  resize(N = 2500, how = "step", time = 800)

pop4 <-
  population("pop4", N = 1500, time = 500, parent = pop3) %>%
  resize(N = 700, how = "exponential", time = 1200, end = 2000)

pop5 <-
  population("pop5", N = 100, time = 600, parent = pop4) %>%
  resize(N = 50, how = "step", time = 900) %>%
  resize(N = 1000, how = "exponential", time = 1600, end = 2200)

model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  generation_time = 1,
  simulation_length = 3000, # forward-time sims need an explicit end
  direction = "forward"
)
```

::: {.aside}
**We started with `pop1` in generation 1, with later events at an
increasing time value.**
:::







## "Forward time units"

```{r}
#| fig-width: 7
#| fig-align: center
plot_model(model) # see time progressing from generation 1 forwards
```

::: {.aside}
**We started with `pop1` in generation 1, with later events at an
increasing time value.**
:::
















## "Backward time units"


```{r}
#| code-line-numbers: "|1|22-24"
pop1 <- population("pop1", N = 1000, time = 30000)

pop2 <-
  population("pop2", N = 1000, time = 27000, parent = pop1) %>%
  resize(N = 100, how = "step", time = 20000)

pop3 <-
  population("pop3", N = 1000, time = 26000, parent = pop2) %>%
  resize(N = 2500, how = "step", time = 22000)

pop4 <-
  population("pop4", N = 1500, time = 25000, parent = pop3) %>%
  resize(N = 700, how = "exponential", time = 18000, end = 10000)

pop5 <-
  population("pop5", N = 100, time = 24000, parent = pop4) %>%
  resize(N = 50, how = "step", time = 21000) %>%
  resize(N = 1000, how = "exponential", time = 14000, end = 8000)

model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  generation_time = 10 # (10 time units for each generation)
  # (we don't need to provide `simulation_length =` because
  # "backwards" models end at time 0 by default, i.e. "present-day")
)
```

::: {.aside}
**Same model as before, except now expressed in units of "years before present".**
:::








## "Backward time units"

```{r}
#| fig-width: 7
#| fig-align: center
plot_model(model) # see time progressing from "year" 30000 backwards
```

::: {.aside}
**Same model as before, except now expressed in units of "years before present".**
:::







# So we built a model...

<center>

<h3>... but how do we simulate data from it?</h3>

</center>










## Built-in simulation "engines"

_slendr_ has two simulation "engine scripts" built-in:

-   _msprime_ engine ([_slendr_ source](https://github.com/bodkan/slendr/blob/main/inst/scripts/script.py)) -- R function `msprime()`
-   SLiM engine ([_slendr_ source](https://github.com/bodkan/slendr/blob/main/inst/scripts/script.slim)) -- R function `slim()`

::: {.fragment}

They are designed to "understand" _slendr_ `model`s, meaning that you
can simulate data just with this command:

```{r}
#| eval: false
ts <- msprime(model, sequence_length = 10e6, recombination_rate = 1e-8)
```

:::

::: {.fragment}

<br>

<center><h3>**No need to write any _msprime_ or SLiM code!**</h3></center>

:::








# The result of a simulation is a tree sequence (`ts`)













## What is tree sequence?

![](files/slendr/images/tree_sequence_diagram_no_muts.png){width="80%" fig-align="center"}

-   a record of full genetic ancestry of a set of samples
-   an encoding of DNA sequence carried by those samples
-   an efficient analysis framework










# Why tree sequence?

<br>

<h3>Why not VCF or a normal genotype table?</h3>










## What we usually have

<center>![](files/slendr/images/vcf_screenshot.png){width="90%"}</center>










## What we usually _want_

An understanding of our samples' evolutionary history:

<center>![](files/slendr/images/tree_sequence_diagram_no_muts.png)</center>

::: {.fragment}
<center>

<h3>**This is exactly what a tree sequence *is*!**</h3>

</center>
:::


::: {.aside}
<small>Image from the [_tskit_ documentation](https://tskit.dev/tutorials/what_is.html)</small>
:::













## The magic of tree sequences

They allow us to compute statistics _without genotypes_!

<center>![](files/slendr/images/tree_sequence_diagram_no_muts.png)</center>

There is a "duality" between mutations and branch lengths.

::: aside
**Note:** See an amazing paper by [Ralph _et al._ (2020)](https://academic.oup.com/genetics/article/215/3/779/5930459) for more detail.
:::












## What if we need mutations though?

&nbsp;

<center>

![](files/slendr/images/tree_sequence_diagram_no_muts.png)












## What if we need mutations though?

Coalescent and mutation processes can be decoupled!

<center>

![](files/slendr/images/tree_sequence_diagram.png)

::: {.fragment}
<h3>This means we can add mutations to `ts` _after_ the simulation using `ts_mutate()`.</h3>
:::

</center>












## Let's go back to our example `model`...

```{r}
#| eval: false
plot_model(model)
```

```{r}
#| fig-align: center
#| echo: false
pop1 <- population("pop1", N = 1000, time = 1)

pop2 <-
  population("pop2", N = 1000, time = 300, parent = pop1) %>%
  resize(N = 100, how = "step", time = 1000)

pop3 <-
  population("pop3", N = 1000, time = 400, parent = pop2) %>%
  resize(N = 2500, how = "step", time = 800)

pop4 <-
  population("pop4", N = 1500, time = 500, parent = pop3) %>%
  resize(N = 700, how = "exponential", time = 1200, end = 2000)

pop5 <-
  population("pop5", N = 100, time = 600, parent = pop4) %>%
  resize(N = 50, how = "step", time = 900) %>%
  resize(N = 1000, how = "exponential", time = 1600, end = 2200)

gf <- gene_flow(from = pop1, to = pop2, start = 2000, end = 2200, rate = 0.13)

model <- compile_model(
  populations = list(pop1, pop2, pop3, pop4, pop5),
  gene_flow = gf,
  generation_time = 1,
  simulation_length = 3000,
  direction = "forward"
)

plot_model(model, proportions = TRUE)
```










## ... simulate a tree sequence...

<br>

In our script we'll have something like this:

```{r}
#| eval: false
#| code-line-numbers: "1-8|10-11"
library(slendr)
init_env()

# <... population() definitions ...>

# <... gene_flow() definition ...>

# <... compile_model(...) ...>
  
ts <-
  msprime(model, sequence_length = 50e6, recombination_rate = 1e-8)
```










## ... and overlay mutations on it

<br>

In our script we'll have something like this:

```{r}
#| eval: false
#| code-line-numbers: "10-12"
library(slendr)
init_env()

# <... population() definitions ...>

# <... gene_flow() definition ...>

# <... compile_model(...) ...>
  
ts <-
  msprime(model, sequence_length = 50e6, recombination_rate = 1e-8) %>%
  ts_mutate(mutation_rate = 1e-8)
```

```{r}
#| echo: false
ts <-
  msprime(model, sequence_length = 1e6, recombination_rate = 1e-8) %>%
  ts_mutate(mutation_rate = 1e-8)
```

::: {.aside}

**Note:** In some exercises, mutations won't be necessary. Where we will need them,
you can use `ts_mutate()` using the pattern shown here.</h4>
:::








# So we can simulate data

<center>

<h3>How do we work with this `ts` thing?</h3>

</center>

















## _slendr_'s R interface to [_tskit_](https://tskit.dev/tskit) statistics

<center>![](files/slendr/images/slendr_tskit.png)</center>

Allele-frequecy spectrum, $\pi$, $f$-statistics, $F_{ST}$, Tajima's D, etc.

**Find help at [slendr.net/reference](https://slendr.net/reference) or in R under `?ts_fst` etc.**









## Extracting names of recorded samples

1. We can get individuals recorded in `ts` with `ts_samples()`:

```{r}
#| eval: false
ts_samples(ts) %>% head(1) # returns a data frame (one row here, for brevity)
```

```{r}
#| echo: false
ts_samples(ts) %>% as.data.frame %>% head(1)
```

::: {.fragment}
2. A shortcut `ts_names()` can also be useful:

```{r}
ts_names(ts) %>% head(5) # returns a vector of individuals' names
```

:::

::: {.fragment}

3. We can get a per-population list of individuals like this:

```{r}
#| eval: false
ts_names(ts, split = "pop") # returns a named list of such vectors
```

```{r}
#| echo: false
ts_names(ts, split = "pop") %>% lapply(sample, 5) %>% .[1]
```

:::
















# 

<center>

<h2>All _slendr_ statistics take individuals' names as their function arguments.</h2>

</center>

<br>

<center>This is modelled after the `sample_sets=` argument of the respective [_tskit_ Python methods](https://tskit.dev/tskit/docs/stable/python-api.html#statistics) (except you use names of individuals directly, not tree-sequence node numbers).</center>

<br>
<center>Let's take a look at how it works in general...</center>











## _tskit_ computation -- option #1

**For a function which operates on one set of individuals**, we can first get
a vector of names to compute on like this:

```{r}
#| eval: false
# a random selection of names of three individuals in a tree sequence
samples <- c("popX_1", "popX_2", "popY_42")
```

::: {.fragment}

<br>

Then we can calculate the statistic of interest like this:

```{r}
#| eval: false
# this computes nucleotide diversity in our set of individuals
df_result <- ts_diversity(ts, sample_sets = list(samples))
```

:::

::: {.aside}
**Note:** Wherever you see `list(<vector of names>)`, you can think of it as
_"compute a statistic for the entire group of individuals"_ (you get a single
number). Without the `list()`, it would mean _"compute the statistic for each
individual separately"_ (and get a value for each of them individually).
:::








## _tskit_ computation -- option #2

For a function operating on multiple sets of individuals,
we want a list of vectors of names (one such vector per group):

```{r}
#| eval: false
# when we compute on multiple groups, it's a good idea to name them
samples <- list(
  popX = c("popX_1", "popX_2", "popX_3"),
  popY = c("popY_1", "popY_2", "popY_3"),
  popZ = c("popZ_1", "popZ_2")
)
```

<br>

Then we use this list of vectors in the same way as before:

```{r}
#| eval: false
# this computes a pairwise divergence between all three groups
df_result <- ts_divergence(ts, sample_sets = samples)
```










## _tskit_ computation -- option #3

For something like $f$ statistics, the function arguments must be more
precisely specified (here `A`, `B`, `C`, not `sample_sets`):

```{r}
#| eval: false
df_result <- ts_f3(
  ts,
  A = c("popX_1", "popX_2", "popX_3"),
  B = c("popY_1", "popY_2", "popY_3"),
  C = c("popZ_1", "popZ_2")
)
```

::: {.fragment}

Doing this manually can be annoying --- `ts_names()` helps by preparing 
the list of names in the correct format:

```{r}
#| eval: false
# get names of individuals in each population as a named list of vectors
samples <- ts_names(ts, split = "pop")

# use this list directly by specifying which vectors to take out
ts_f3(ts, A = samples$popX, B = samples$popY, C = samples$popZ)
```

:::












#

<center>

<h1>Some examples on simulated data</h1>

<br><h3>(A tree sequence `ts` we got earlier.)</h3>

</center>











## Example: [nucleotide diversity](https://en.wikipedia.org/wiki/Nucleotide_diversity)

```{r}
#| echo: false
pop1 <- population("pop1", N = 10000, time = 1)
pop2 <- population("pop2", N = 1000, time = 5000, parent = pop1)

model <- compile_model(
  populations = list(pop1, pop2),
  generation_time = 1,
  simulation_length = 10000
)

ts <-
  msprime(model, sequence_length = 10e6, recombination_rate = 1e-8) %>%
  ts_mutate(mutation_rate = 1e-8)
```

```{r}
#| echo: false
set.seed(42)
```

::: {.columns} 
::: {.column width="45%"}
Get a list of individuals in each population:

```{r}
samples <- ts_names(ts, split = "pop")

names(samples)
```

::: {.fragment}
<br>

We can index into the list via population name:

```{r}
#| eval: false
samples$pop1
```

```{r}
#| echo: false
samples$pop1 %>% head(3)
```

```{r}
#| eval: false
samples$pop2
```

```{r}
#| echo: false
samples$pop2 %>% head(3)
```


:::
:::

::: {.column width="2%"}
 
:::

::: {.column width="53%"}
::: {.fragment}
Compute nucleotide diversity (note the list `samples`):

```{r}
ts_diversity(ts, sample_sets = samples)
```

<br>

Our tree sequence had two populations, `pop1` and `pop2`, which is why we
get a data frame with diversity in each of them.

:::
:::
:::










## Example: [allele frequency spectrum](https://en.wikipedia.org/wiki/Allele_frequency_spectrum)

```{r}
#| echo: false
pop <- population("pop", N = 10000, time = 1)

model <- compile_model(pop, generation_time = 1, simulation_length = 10000)

ts <-
  msprime(model, sequence_length = 10e6, recombination_rate = 1e-8) %>%
  ts_mutate(mutation_rate = 1e-8)
```

```{r}
#| echo: false
set.seed(42)
```

::: {.columns} 
::: {.column width="45%"}
Get names of individuals:

```{r}
samples <- ts_names(ts)[1:5]
samples
```

::: {.fragment}
Compute the AFS:

```{r}
afs <- ts_afs(ts, sample_sets = list(samples))

# we skip the 1st item because it has a special meaning in tskit
afs[-1]
```
:::
:::

::: {.column width="2%"}
 
:::

::: {.column width="53%"}
::: {.fragment}
```{r}
#| eval: false
plot(afs[-1], type = "b",
     xlab = "allele count bin",
     ylab = "frequency")
```

```{r}
#| echo: false
#| fig-height: 7
plot(afs[-1], type = "b",
     xlab = "allele count bin",
     ylab = "frequency", lwd = 3,
     cex = 1.5, cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
```
:::
:::
:::


::: {.aside}
**Note:** One of the rare examples when a _slendr_ / _tskit_ statistical 
function does not return a data frame (`ts_afs()` returns a numerical vector,
not a data frame).
:::





#

<h2>More details on tree-sequences</h2>



## Tree sequence tables

<br>

This simulates 2 $\times$ 10000 chromosomes of 100 Mb:

```{r}
#| eval: false
library(slendr)
init_env(quiet = FALSE)

pop <- population("pop", time = 100e6, N = 10000)
model <- compile_model(pop, generation_time = 30, direction = "backward")
ts <- msprime(model, sequence_length = 1e6, recombination_rate = 1e-8)
```

```{r}
#| echo: false
library(slendr)
init_env(quiet = FALSE)

pop <- population("pop", time = 1e6, N = 10000)
model <- compile_model(pop, generation_time = 30, direction = "backward")
ts <- msprime(model, sequence_length = 1e6, recombination_rate = 1e-8)
```


<br>

**Runs in less than 30 seconds on my laptop!**

**Takes only about 66 Mb of memory!**

## How is this even possible?!


<center>

![](files/slendr/images/tables.jpeg)

<center>

## Tree-sequence tables

::: row

::: {.columns}

::: {.column width="60%"}

A tree can be represented by

::: incremental

-   a table of <font color="orange">n</font><font color="green">o</font><font color="darkblue">d</font><font color="green">e</font><font color="darkblue">s</font>,

-   a table of [edges]{.underline} between nodes,

-   a table of <font color="red">mutations</font> on edges

:::

:::

::: {.column width="40%"}

<br>

<center>![](files/slendr/images/tree_diagram.png)</center>

:::

:::

:::

<h3>**A collection of such tables is a tree sequence.**</h3>

::: {.aside}

**Note:** This is a huge oversimplification.
Find more information in  [tskit docs](https://tskit.dev/tutorials/tables_and_editing.html).

:::

## Tree-sequence tables in practice

```{r}
#| echo: false
#| fig-height: 5
#| fig-align: center
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(ggtree)
})

set.seed(123)

ts <- msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, random_seed = 42) %>% ts_mutate(1e-8, random_seed = 42)

# make a tiny example simplified tree sequence

ts_tiny <- ts_samples(ts) %>% sample_n(4) %>% pull(name) %>% ts_simplify(ts, simplify_to = .)

# extract tree #1 as an ape tree and also a tskit tree

t_phylo <- ts_phylo(ts_tiny, 1, quiet = TRUE)
t_tskit <- ts_tree(ts_tiny, 1)

# plot the phylo tree with ape

nodes <- ts_nodes(t_phylo) %>% as_tibble %>% dplyr::select(node = phylo_id, pop, node_id)

ggtree(t_phylo, branch.length = "none") %<+% nodes +
  geom_label(aes(label = node_id)) +
  guides(color = "none")
```

::: {.columns}

::: {.column width="30%"}

::: {.fragment}

**Nodes:**

```{r}
#| echo: false
ts_nodes(t_phylo) %>% head(3) %>% .[, c("node_id", "time")] %>% as.data.frame() %>%
  setNames(c("node", "time"))
```

:::

:::

::: {.column width="30%"}

::: {.fragment}

**Edges:**

```{r}
#| echo: false
ts_edges(t_phylo) %>% head(3) %>% .[, c("child_node_id", "parent_node_id")] %>% setNames(c("child", "parent")) %>% as.data.frame()
```

:::

:::

::: {.column width="30%"}

::: {.fragment}

**Mutations:**

```{r}
#| echo: false
ts_table(ts_tiny, "mutations") %>% filter(node %in% c(53, 22, 20, 74, 9)) %>% head(3) %>% .[, c("id", "node", "time")] %>% as.data.frame()
```

:::

:::

:::


# Other data formats

<br>

Tree sequence is a useful, cutting-edge data structure, but there are many
well-established bioinformatics tools out there.

<br>

_tskit_ and _slendr_ offer a couple of functions to help integrate their simulation
results into third-party software.



```{r}
#| echo: false
library(slendr)
init_env(quiet = TRUE)

pop <- population("pop", time = 1e6, N = 100)
model <- compile_model(pop, generation_time = 30, direction = "backward")
ts <- msprime(model, sequence_length = 10e6, recombination_rate = 1e-8) %>%
  ts_mutate(1e-8)
```



## Standard genotype formats

If a tree sequence doesn't cut it, you can always...

::: {.fragment}

-   export genotypes to a VCF file:

```{r}
#| eval: false
ts_vcf(ts, path = "path/to/a/file.vcf.gz")
```

:::

::: {.fragment}

-   export genotypes in the EIGENSTRAT format:

```{r}
#| eval: false
ts_eigenstrat(ts, prefix = "path/to/eigenstrat/prefix")
```

:::

::: {.fragment}

-   access genotypes as a data frame:

```{r}
#| eval: false
ts_genotypes(ts)
```

```{r}
#| echo: false
library(dplyr)
ts_genotypes(ts) %>% as.data.frame() %>% .[1:2, 1:7]
```

:::