---
title: "Selection tests"
author: "Karl Schmid, University of Hohenheim"
date: today
subtitle: 'Module: Plant Genetic Resources (3502-470)'
---

# R packages for this computer lab

For this computer lab the following R packages are required:

* `pegas` 
* `phyclust`
* `gap`

Make sure that they are installed and loaded

```{r}
library(pegas)
library(phyclust)
library(gap)
```


# The data set

We continue to use our data set with teosinte individuals, maize landraces and maize improved varieties. 

Download the data from the website:

- [https://web350b.ipsp.uni-hohenheim.de/pgr/data/zea_snps.RData](https://web350b.ipsp.uni-hohenheim.de/pgr/data/zea_snps.RData)
- [https://web350b.ipsp.uni-hohenheim.de/pgr/data/teosinte_acc.Rdata](https://web350b.ipsp.uni-hohenheim.de/pgr/data/teosinte_acc.Rdata)
- [https://web350b.ipsp.uni-hohenheim.de/pgr/data/landraces_acc.Rdata](https://web350b.ipsp.uni-hohenheim.de/pgr/data/landraces_acc.Rdata)
- [https://web350b.ipsp.uni-hohenheim.de/pgr/data/varieties_acc.Rdata](https://web350b.ipsp.uni-hohenheim.de/pgr/data/varieties_acc.Rdata)

First, load the data set and sample lists that you saved as an `Rdata` objects in the first computer lab:



```{r}
load("data/zea_snps.RData")
load("data/teosinte_acc.Rdata")
load("data/landraces_acc.Rdata")
load("data/varieties_acc.Rdata")
```

# Background

The goal of this computer lab is to conduct a test of neutral evolution with our real data set, i.e. we want to know if this 3 Mb region on chromosome 1 evolves neutrally or whether it is influenced by either demographic events or selection. 

Next to the estimators of the mutation parameter $\theta=4N_e\mu$, the Watterson's $\theta_W$ and the nucleotide diversity $\Pi$, we will also calculate the neutrality test statistic Tajima’s $D$. 
The latter will allow us to detect departures from a coalescent model with neutral mutation.

# Tajima's $D$

We will calculate first Tajima's $D$ for our sample sets:

```{r}
pegas::tajima.test(zea_snps[teosinte_acc,])
tajima.test(zea_snps[landraces_acc,])
tajima.test(zea_snps[varieties_acc,])
```

Write the values into a table like this one:

|data set          | Tajima's $D$ | P-value |
|:-----------------|-------------:|--------:|
|teosinte          |   -0.8907504  | 0.3730631 |
|landraces         |   -0.7142114 | 0.4750964 |
|improved varieties|    0.3525113 | 0.7244548 |

:::{.callout-note title="Exercise"}

1. What conclusions can you draw for the three sample sets from the Tajima's $D$ values?
2.Do those values rather allow you to draw conclusions about the demographic past of the samples or about selection? Discuss!

:::

# Simulations

We will now simulate data sets with the parameters from our true data sample sets (sample size, $\theta_W$), calculate $\theta_W$, $\pi$ and Tajima's $D$ for all simulated data sets and then plot the distributions and compare those to the observed values. 
This is analogous to what we did in the last computer lab.

(*Note that although we do only 300 simulations, run times of some of the commands can still be a few minutes.*)

::: {.callout-warning}
If Rstudio crashes with 300 simulations, reduce the number of simulations to 200 or even 100 or 50... by setting the variable `sim_count`below.
If it keeps crashing, move on to the next part on page 4, in which we do the same for a specific region within out data set.
:::

First we simulate the data sets with `ms`:

```{r}
sim_count <- 50 # number of simulations to run

# teosinte
sim_teosinte <- ms(length(teosinte_acc), sim_count, paste("-t", theta.s(zea_snps[teosinte_acc,])))
seq_teosinte <- read.ms.output(sim_teosinte,is.file = FALSE)
gametes_teosinte <- seq_teosinte$gametes

# landraces
sim_landraces <- ms(length(landraces_acc), sim_count, paste("-t",theta.s(zea_snps[landraces_acc,])))
seq_landraces <- read.ms.output(sim_landraces,is.file = FALSE)
gametes_landraces <- seq_landraces$gametes

# varieties
sim_varieties <- ms(length(varieties_acc), sim_count, paste("-t",theta.s(zea_snps[varieties_acc,])))
seq_varieties <- read.ms.output(sim_varieties, is.file = FALSE)
gametes_varieties <- seq_varieties$gametes

```

And then we recode the output into a `DNAbin` objects:

```{r}
f.recode <- function(m){ # Function to recode sequence as mock DNA sequence
  m[m==0] <- "a"
  m[m==1] <- "c"
  as.DNAbin(t(m))
}

# teosinte
gametesseq_teosinte <- lapply(gametes_teosinte, f.recode)
# landraces
gametesseq_landraces <- lapply(gametes_landraces, f.recode)
# varieties
gametesseq_varieties <- lapply(gametes_varieties, f.recode)
```

We can now plot the distribution of $\theta_W$ values; the red lines will show the mean of the distributions, the blue lines the observed values:

```{r}
library(ggplot2)

# 1. Gather your theta estimates into one data frame
df <- data.frame(
  theta = c(
    sapply(gametesseq_teosinte, theta.s),
    sapply(gametesseq_landraces, theta.s),
    sapply(gametesseq_varieties, theta.s)
  ),
  group = factor(rep(
    c("teosinte", "landraces", "varieties"),
    times = c(
      length(gametesseq_teosinte),
      length(gametesseq_landraces),
      length(gametesseq_varieties)
    )
  ), levels = c("teosinte", "landraces", "varieties"))
)

# 2. Compute the two reference lines per group
summary_df <- data.frame(
  group = factor(c("teosinte", "landraces", "varieties"),
                 levels = c("teosinte", "landraces", "varieties")),
  mean_theta = c(
    mean(sapply(gametesseq_teosinte, theta.s)),
    mean(sapply(gametesseq_landraces, theta.s)),
    mean(sapply(gametesseq_varieties, theta.s))
  ),
  theta_pop = c(
    theta.s(zea_snps[teosinte_acc,]),
    theta.s(zea_snps[landraces_acc,]),
    theta.s(zea_snps[varieties_acc,])
  )
)

# 3. Plot with facets
ggplot(df, aes(x = theta)) +
  geom_histogram(bins = 20, fill = "grey80", color = "black") +
  geom_vline(
    data = summary_df,
    aes(xintercept = mean_theta),
    color = "red", size = 1
  ) +
  geom_vline(
    data = summary_df,
    aes(xintercept = theta_pop),
    color = "blue", linetype = "dashed", size = 1
  ) +
  facet_wrap(~ group, ncol = 1, scales = "free_y") +
  labs(
    x = expression(theta),
    y = "Count"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    strip.text = element_text(face = "bold"),
    panel.spacing = unit(1, "lines")
  )


```


These are the plots for the $\pi$ values:

```{r}
library(ggplot2)

# 1. Define your per-locus $\pi$ function
locus_nuc_div <- function(x) {
  nuc.div(x) * ncol(x)
}

# 2. Build a single data frame of all $\pi$-values
df_pi <- data.frame(
  pi = c(
    sapply(gametesseq_teosinte,    locus_nuc_div),
    sapply(gametesseq_landraces,   locus_nuc_div),
    sapply(gametesseq_varieties,   locus_nuc_div)
  ),
  group = factor(rep(
    c("teosinte", "landraces", "varieties"),
    times = c(
      length(gametesseq_teosinte),
      length(gametesseq_landraces),
      length(gametesseq_varieties)
    )
  ), levels = c("teosinte", "landraces", "varieties"))
)

# 3. Compute the two vertical‐line positions per group
summary_pi <- data.frame(
  group   = factor(c("teosinte", "landraces", "varieties"),
                   levels = c("teosinte", "landraces", "varieties")),
  mean_pi = c(
    mean(sapply(gametesseq_teosinte,  locus_nuc_div)),
    mean(sapply(gametesseq_landraces, locus_nuc_div)),
    mean(sapply(gametesseq_varieties, locus_nuc_div))
  ),
  theta_ref = c(
    theta.s(zea_snps[teosinte_acc,]),
    theta.s(zea_snps[landraces_acc,]),
    theta.s(zea_snps[varieties_acc,])
  )
)

# 4. Plot
ggplot(df_pi, aes(x = pi)) +
  geom_histogram(bins = 20, fill = "grey80", color = "black") +
  # red line = sample mean $\pi$
  geom_vline(data = summary_pi,
             aes(xintercept = mean_pi),
             color = "red", size = 1) +
  # blue dashed = reference $\hat{\theta}$
  geom_vline(data = summary_pi,
             aes(xintercept = theta_ref),
             color = "blue", linetype = "dashed", size = 1) +
  facet_wrap(~ group, ncol = 1, scales = "free_y") +
  labs(
    x = expression(pi),
    y = "Count"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    strip.text = element_text(face = "bold"),
    panel.spacing = unit(1, "lines")
  )

```


And for the Tajima's $D$ values:

```{r}
library(ggplot2)

# 1. Define your Tajima’s D extractor
tajD <- function(v) {
  tajima.test(v)$D
}

# 2. Build a single data frame of all D-values
df_D <- data.frame(
  D = c(
    sapply(gametesseq_teosinte, tajD),
    sapply(gametesseq_landraces,  tajD),
    sapply(gametesseq_varieties,  tajD)
  ),
  group = factor(rep(
    c("teosinte", "landraces", "varieties"),
    times = c(
      length(gametesseq_teosinte),
      length(gametesseq_landraces),
      length(gametesseq_varieties)
    )
  ), levels = c("teosinte", "landraces", "varieties"))
)

# 3. Compute per-group reference lines
summary_D <- data.frame(
  group    = factor(c("teosinte", "landraces", "varieties"),
                    levels = c("teosinte", "landraces", "varieties")),
  mean_D   = c(
    mean(sapply(gametesseq_teosinte, tajD)),
    mean(sapply(gametesseq_landraces,  tajD)),
    mean(sapply(gametesseq_varieties,  tajD))
  ),
  D_ref    = c(
    tajD(zea_snps[teosinte_acc,]),
    tajD(zea_snps[landraces_acc,]),
    tajD(zea_snps[varieties_acc,])
  )
)

# 4. Plot with facets
ggplot(df_D, aes(x = D)) +
  geom_histogram(bins = 20, fill = "grey80", color = "black") +
  # red: group mean
  geom_vline(data = summary_D,
             aes(xintercept = mean_D),
             color = "red", size = 1) +
  # blue dashed: reference Tajima's D
  geom_vline(data = summary_D,
             aes(xintercept = D_ref),
             color = "blue", linetype = "dashed", size = 1) +
  facet_wrap(~ group, ncol = 1, scales = "free_y") +
  labs(
    x = "Tajima’s D",
    y = "Count"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    strip.text = element_text(face = "bold"),
    panel.spacing = unit(1, "lines")
  )

```


:::{.callout-note Title="Exercise"}

1. How well do the observed values fit with the simulated values? Are there any differences between the sample sets?
2. Can you reject the null hypothesis (= neutral evolution) for any of the sample sets? Why or why not?
:::

# Fine-scale signatures of selection

So far, we have only looked at the entire 3Mb region on chromosome 1. 
Sliding windows calculated with the R package `PopGenome` for nucleotide diversity, $\pi$, and Tajima's $D$ show fluctuations along this region (see figure). We will now closer investigate a small region which shows a reduced genetic diversity and a lower Tajima's $D$ in the landraces and the varieties but not in teosinte. The dashed line in the figure marks this region.

![Changes in diversity along the chromosome.](img/div_sw_clab.png){#fig-divsw}

For this we will first extract the respective region from our `DNAbin` object:

```{r}
zea_region <- zea_snps[,10735:11426] # extracts the SNPs from the DNAbin object
```

We can calculate Tajima'S $D$ for this region:

```{r}
tajima.test(zea_region[teosinte_acc,])
tajima.test(zea_region[landraces_acc,])
tajima.test(zea_region[varieties_acc,])
```

Write the values into a table like this one:

|data set          | Tajima's $D$ | p-value
|:-----------------|-------------:|----------|
|teosinte          |    0.5925368   |   0.5534912     |
|landraces         |   -2.335788    |  0.01950229      |
|improved varieties|   -2.329903   |  0.0198113       |


:::{.callout-note title="Exercise"}

What conclusions can you draw for this region for the three sample sets?

:::


And then we will repeat the simulations from before:

```{r}
## simulations
sim_count <- 50 # number of simulations to run

# teosinte
sim_teosinte <- ms(length(teosinte_acc), sim_count, paste("-t",theta.s(zea_region[teosinte_acc,])))
seq_teosinte <- read.ms.output(sim_teosinte, is.file = FALSE)
gametes_teosinte <- seq_teosinte$gametes

# landraces
sim_landraces <- ms(length(landraces_acc), sim_count, paste("-t",theta.s(zea_region[landraces_acc,])))
seq_landraces <- read.ms.output(sim_landraces, is.file = FALSE)
gametes_landraces <- seq_landraces$gametes

# varieties
sim_varieties <- ms(length(varieties_acc),sim_count, paste("-t",theta.s(zea_region[varieties_acc,])))
seq_varieties <- read.ms.output(sim_varieties, is.file = FALSE)
gametes_varieties <- seq_varieties$gametes

## recoding in DNAbin objects

f.recode <- function(m){ # Function to recode sequence as mock DNA sequence
  m[m==0] <- "a"
  m[m==1] <- "c"
  as.DNAbin(t(m))
}

# teosinte
gametesseq_teosinte <- lapply(gametes_teosinte, f.recode)
# landraces
gametesseq_landraces <- lapply(gametes_landraces, f.recode)
# varieties
gametesseq_varieties <- lapply(gametes_varieties, f.recode)
```

We plot the distribution of Tajima's $D$ values; again the red lines will show the mean of the distributions and the blue lines the observed values:


```{r}
library(ggplot2)

# 1. Tajima’s D extractor
tajD <- function(v) {
  tajima.test(v)$D
}

# 2. Stack all D-values into one data frame
df_D <- data.frame(
  D = c(
    sapply(gametesseq_teosinte, tajD),
    sapply(gametesseq_landraces,  tajD),
    sapply(gametesseq_varieties,  tajD)
  ),
  group = factor(rep(
    c("teosinte", "landraces", "varieties"),
    times = c(
      length(gametesseq_teosinte),
      length(gametesseq_landraces),
      length(gametesseq_varieties)
    )
  ), levels = c("teosinte", "landraces", "varieties"))
)

# 3. Compute per‐group means and reference D from zea_region
summary_D <- data.frame(
  group    = factor(c("teosinte", "landraces", "varieties"),
                    levels = c("teosinte", "landraces", "varieties")),
  mean_D   = c(
    mean(sapply(gametesseq_teosinte, tajD)),
    mean(sapply(gametesseq_landraces,  tajD)),
    mean(sapply(gametesseq_varieties,  tajD))
  ),
  D_ref    = c(
    tajD(zea_region[teosinte_acc,]),
    tajD(zea_region[landraces_acc,]),
    tajD(zea_region[varieties_acc,])
  )
)

# 4. Plot with vertical facets and reference lines
ggplot(df_D, aes(x = D)) +
  geom_histogram(bins = 20, fill = "grey80", color = "black") +
  # red = mean of sample
  geom_vline(data = summary_D,
             aes(xintercept = mean_D),
             color = "red", size = 1) +
  # blue dashed = tajD from zea_region
  geom_vline(data = summary_D,
             aes(xintercept = D_ref),
             color = "blue", linetype = "dashed", size = 1) +
  facet_wrap(~ group, ncol = 1, scales = "free_y") +
  labs(
    x = "Tajima’s D",
    y = "Count"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    strip.text   = element_text(face = "bold"),
    panel.spacing = unit(1, "lines")
  )

```


:::{.callout-note title="Exercise"}

1. How well do the observed values fit with the simulated values? Are there any differences between the sample sets?
2. Can you reject the null hypothesis (= neutral evolution) for any of the sample sets? Why?
:::

# Further reading

The concept of this computer lab has been used in the following two publications:

- Schmid et al. (2005) A multilocus sequence survey in *Arabidopsis thaliana* reveals a genome-wide departure from a neutral model of DNA sequence polymorphism. [Genetics](http://www.genetics.org/content/169/3/1601.short)
- Wright et al. (2005) The effects of artificial selection on the maize genome.  [Science](http://science.sciencemag.org/content/308/5726/1310)


## References {.unnumbered}
\begingroup
\raggedright
\small
::: {#refs}
:::
\endgroup


