---
title: "Tools for population genetics"
author: "Matheus Januario and Jennifer Auler"
date: "Oct/2023"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{popgen_old}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r include=FALSE}
options(rmarkdown.html_vignette.check_title = FALSE)
```


```{r}
library(evolved)
```

## Hardy-Weinberg equilibrium

The Hardy-Weinberg equilibrium (HWE) is a deterministic model for the maintenance of allele-level biodiversity. However, stochastic thinking is already well-established as an important aspect of studying evolutionary biology - and so there are clear advantages in introducing this kind of thinking to students.

However, typical HWE equations or study material related to it do not address stochastic variation coming from finite numbers of individuals. The function `OneGenHWSim()` addresses this issue. With it, students can graph and explore stochasticity under HWE:

```{r fig.height=5, fig.width=6, fig.align='center'}
nsims = 1000
sims <- OneGenHWSim(n.ind = 57, p = .2, n.sim = nsims)

boxplot(sims, frame.plot=F,
        ylab="Number of individuals", xlab="Genotype")
for(i in 1:nrow(sims)){
  lines(x=1:3, y=sims[i, ], col="#A020F064")  
}
```

Of course, the object generated by the simulations can also be manipulated as if it was "data"

```{r}
head(sims)
```

With this type of object, students can also simulate closed populations of a specific size, thus for instance testing (using e.g. chi-square tests) actual hypotheses about HWE adequacy to explain a certain genotypic frequency.

## Genetic drift

Students can also simulate genetic drift in Wright-Fisher populations, by using the function `WFDriftSim()`

```{r fig.height=5, fig.width=6, fig.align='center'}
drifting = WFDriftSim(Ne = 10, n.gen = 100, p0 = .5, n.sim = 1000, 
                        print.data = T, plot.type = "static", knitr = TRUE)

#View the first 10 generation of a few simulations.
drifting[1:6, 1:11]
```

The argument `printData = T` creates an object that can also be handled as it was "data", so students can use it to calculate drift processes, for instance, the amount of time needed until fixation:

```{r,fig.height=5, fig.width=6, fig.align='center'}
fixchecker <- function(x){
  res = suppressWarnings(
      min(which(x %in% 0:1), na.rm = T)
    )
    return(res)  
}

hist(x = apply(X = drifting, MARGIN = 1, fixchecker),
     main="Time until fixation", xlab = "Fixation time")
```

Or other types of manipulations or measurements.

## Natural selection

Just like the previously mentioned functions, `NatSelSim()` can generate objects that can be handled by students as if they were "data". 

```{r,fig.height=5, fig.width=6, fig.align='center'}
selected = NatSelSim(w11 = .2, w12 = .3, w22 = .1,
            n.gen = 20, print.data = T, plot.type = "static", knitr= TRUE)
```

Note that by default the function provides 4 animated plots, but this can be changed using the argument `plot_type` (see function help for details).

Now, we can ask: Assuming we simulated the dynamics for enough time, which genotype in the simulation above reached its equilibrium frequency faster?

```{r}
# a numeric tolerance to accept some frequency as close enough to equilibrium
tol = 0.0001 

aux = abs(apply(selected, 2, diff))

# answer (in number of generations)
apply(aux>tol, 2, which.min)
```

In the above result, we are tabulating the number of dominant homozygotes ("AA"), hetrozygotes ("Aa"), and recessive homozygotes ("aa"). Students can manipulate the objects freely, and thus educators can allow them to use simulations to explore relevant scenarios in an investigative manner (see types of inquiries in Banchi & Bell, 2008)

## References:

Banchi, H., & Bell, R. (2008). The many levels of inquiry. Science and children, 46(2), 26.