---
title: 'Analysis of the genetic structuration of bakery yeasts using microsatellite genotyping data'
author: '[Frédéric Bigey](frederic.bigey@inra.fr)'
date: "`r Sys.Date()`"
output:
  pdf_document: null
  df_print: kable
  html_document:
    df_print: paged
fontsize: 12pt
documentclass: article
classoption: a4paper
---
    
```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(
    include   = TRUE,
    echo      = TRUE,
    warning   = TRUE,
    message   = FALSE, 
    fig.align = "center", 
    fig.width = 3.0,
    fig.asp   = 1.0,
    fig.path  = 'Figs/',
    cache     = TRUE,
    size      = "tiny"
)
options(knitr.table.format = "latex")
```


# MAT & MET

## Materials

The microsatellite genotype file in GenAlEx format was received from Delphine SICARD. $323$ yeasts were genotyped at $15$ loci.

```{r echo=TRUE}
input.file = "../Pierre GOMES/bakery_microsat_alphanumeric_sorted_samples.genalex.tsv"
info.file  = "../Pierre GOMES/Collection_microsat_bakery_23-mai-2019.xls"
```


## Packages

Most of analysis performed in this report were obtained using the `poppr` and `adegenet` packages.

```{r}
library(dplyr)
library(poppr)
```

## Import of the dataset 

```{r eval=TRUE}
# read in a genAlEx-formated file, 
# return a genind or genclone object
data = poppr::read.genalex(
    input.file, 
    ploidy = 4, 
    sep = "\t",
    recode = FALSE
)

# read informations stored in Excel file
info.data = readxl::read_excel(path = info.file, sheet = 1, col_names = TRUE, na = "NA")

# strain name subtitution
indNames(data) %>% head()
indNames(data) = info.data[["Name"]]
indNames(data) %>% head()

# select strata informations to keep
strata.data = info.data %>% 
  select(Name, Habitat, Isolation, Country, Bakery, Sourdough, Ploidy, Include)

# set strata using the data frame
strata(data) = strata.data
strata(data) %>% head(4)

# clean
rm(strata.data, input.file, info.file)
```


Keep only strains containing the I(nclude) label

```{r}
# filter out some samples (Include == "I")
nInd(data)

setPop(data) = ~Include
data = popsub(data, sublist = c("I"))

nInd(data)
```


Here are the length of the microsatellite repeat unit:

```{r}
# set the length of the microsatellite repeats
loci.nucl.repeat = data.frame(
  length = c(2,3,3,2,2,3,3,3,3,3,3,3,3,2,3),
  #length = c(2, 3, 3, 2, 2, 3, 3, 3, 6, 3, 9, 3, 3, 2, 3),
  row.names = locNames(data))
loci.nucl.repeat

# testing the repeat length 
test_replen(data, replen = loci.nucl.repeat$length)
```

## Sampling 3 yeasts per sourdough

```{r, eval= TRUE}
library(tidyr)
library(purrr)

# sampling function used in map2()
sampling = function(tbl.df, nb) {
  if(nb > 2) {
    sample(tbl.df$Name, 3)    
  } else if (nb > 1) {
    sample(tbl.df$Name, 2)
  } else {
    tbl.df$Name
  }
}

# Sample max 3 yeasts per sourdough 
sampled.sourdough.yeasts = info.data %>%
  filter(Isolation == "sourdough" & Include == "I") %>%
  select(Name, Sourdough) %>%
  group_by(Sourdough) %>%
  nest() %>%
  mutate(n = map_dbl(data, nrow)) %>%
  mutate(samples = map2(data, n, sampling))

sampled.sourdough.yeasts.name = sampled.sourdough.yeasts$samples %>% unlist()

# get commercial yeasts name vector
sampled.commercial.yeast.name = info.data %>%
  filter(Isolation == "commercial" & Include == "I") %>%
  select(Name) %>%
  flatten_chr()

# all sampled yeast
sampled.yeast.name = c(sampled.sourdough.yeasts.name, sampled.commercial.yeast.name)

# clean
rm(info.data, sampled.sourdough.yeasts, sampled.sourdough.yeasts.name, sampled.commercial.yeast.name)
```


## Reformat of diploid genotypes

For diploids, it is not sufficiant to change ploidy like `ploidy(diploid) = 2`, it is necessery to convert the fictive tetraploid genotype into a new diploid genotype. In this dataset, diploid genotypes are encoded as autotetraploids. An issue comes from some individuals that have `3` different alleles at some loci (the missing allele is coded `0`). This individuals will be removed from the analysis.

+ Previous encoding:

```{r}
# # filter 2N yeasts
# setPop(data) = ~Ploidy
# diploid = popsub(data, sublist = c("2N"))
# 
# # filter out ind with genotype containing a "0"
# df = genind2df(diploid, sep = "/")
# df %>% head(4)
# 
# tmp = apply(df, MARGIN = 1, FUN = function(x) {sum(unlist(strsplit(x, "/")) == "0", na.rm = TRUE)}) == 0
# tmp = tmp %>% as.data.frame()
# diploids = diploid[tmp[,1],]
# diploids.strata = strata(diploids)
# 
# # divide internal genotype count per 2
# # 4 is tranformed in 2
# # 2 in 1
# diploids@tab = floor(diploids@tab / 2)
# 
# # set ploidy to 2
# df = genind2df(diploids, sep = "/")
# diploids = df2genind(df[,-1], sep = "/")
# ploidy(diploids) = 2
# strata(diploids) = diploids.strata
# diploids
# 
# # clean
# rm(df, tmp, diploid, diploids.strata)

# select 2N and store strata
setPop(data) = ~Ploidy
diploids = popsub(data, sublist = c("2N"))
diplo.strata = strata(diploids)

# recode function
recode.2n = function(x) {
  
    #x = "218/218/210/210"
    x = strsplit(x, split = "/") %>% unlist() %>% unique()
    
    if ( length(x) == 1 && !is.na(x) ) {
      x = c(x,x)
    } else if (length(x) == 2) {
      x= c(x)
    } else if (length(x) >= 3) {
      x = c(NA)
    }
    
    paste(x, collapse = "/")
}

df = genind2df(diploids, sep = "/")[,-1] %>% 
  as_tibble(rownames = "Strain")

df.recoded = df %>% 
  mutate(C11 = purrr::map_chr(C11, recode.2n)) %>%
  mutate(C3 = purrr::map_chr(C3, recode.2n)) %>%
  mutate(C4 = purrr::map_chr(C4, recode.2n)) %>%
  mutate(C5 = purrr::map_chr(C5, recode.2n)) %>%
  mutate(C6 = purrr::map_chr(C6, recode.2n)) %>%
  mutate(C8 = purrr::map_chr(C8, recode.2n)) %>%
  mutate(C9 = purrr::map_chr(C9, recode.2n)) %>%
  mutate(SCAAT1 = purrr::map_chr(SCAAT1, recode.2n)) %>%
  mutate(SCAAT3 = purrr::map_chr(SCAAT3, recode.2n)) %>%
  mutate(SCYOR267C = purrr::map_chr(SCYOR267C, recode.2n)) %>%
  mutate(YKL172W = purrr::map_chr(YKL172W, recode.2n)) %>%
  mutate(YKR072C = purrr::map_chr(YKR072C, recode.2n)) %>%
  mutate(YLL049 = purrr::map_chr(YLL049, recode.2n)) %>%
  mutate(YLR = purrr::map_chr(YLR, recode.2n)) %>%
  mutate(YPL009C = purrr::map_chr(YPL009C, recode.2n))

# convert to data frame and in Genind object
df.recoded = data.frame(df.recoded[,-1], row.names = df.recoded$Strain)
diploids = df2genind(df.recoded, sep = "/", ploidy = 2)
strata(diploids) = diplo.strata
diploids

rm(df, df.recoded, diplo.strata)
```

+ New encoding:

```{r, warning=FALSE}
setPop(diploids) = ~Isolation
genind2df(diploids, sep = "/") %>% head(4)

# export to genAlEx formatted *.csv file
# genind2genalex(
#   diploids, 
#   filename = "data.diploid.genalex.csv", 
#   pop = ~Isolation/Country/Bakery, 
#   allstrata = FALSE, 
#   overwrite = TRUE)
```

This results from a new dataset containing `r adegenet::nInd(diploids)` diploid genotypes.


## Reformat tetraploid genotypes


+ Reformat microsat genotypes

The probleme is that genotypes with biallelic genotype are encoded without missing data. Exemple "123/123/160/160" should be encoded in "123/160/0/0" since we have no info for the last two alleles. Examples:



```{r}
# select 4N and store strata
setPop(data) = ~Ploidy
tetraploids = popsub(data, sublist = c("4N", "aneuploid"))
tetra.strata = strata(tetraploids)
#rm(data)

# exemple
genind2df(tetraploids, sep = "/") %>% head

# recode function
recode = function(x) {
  
    #x = "218/218/210/210"
    x = strsplit(x, split = "/") %>% unlist() %>% unique()
    
    if (length(x) == 2) {
      x= c(000,000,x)
    } else if (length(x) == 1 && !is.na(x)) {
        x = c(x,x,x,x)
    } else if (length(x) == 3) {
      x = c(000,x)
    } 
    
    paste(x, collapse = "/")
}

df = genind2df(tetraploids, sep = "/")[,-1] %>% 
  as_tibble(rownames = "Strain")

df.recoded = df %>% 
  mutate(C11 = purrr::map_chr(C11, recode)) %>%
  mutate(C3 = purrr::map_chr(C3, recode)) %>%
  mutate(C4 = purrr::map_chr(C4, recode)) %>%
  mutate(C5 = purrr::map_chr(C5, recode)) %>%
  mutate(C6 = purrr::map_chr(C6, recode)) %>%
  mutate(C8 = purrr::map_chr(C8, recode)) %>%
  mutate(C9 = purrr::map_chr(C9, recode)) %>%
  mutate(SCAAT1 = purrr::map_chr(SCAAT1, recode)) %>%
  mutate(SCAAT3 = purrr::map_chr(SCAAT3, recode)) %>%
  mutate(SCYOR267C = purrr::map_chr(SCYOR267C, recode)) %>%
  mutate(YKL172W = purrr::map_chr(YKL172W, recode)) %>%
  mutate(YKR072C = purrr::map_chr(YKR072C, recode)) %>%
  mutate(YLL049 = purrr::map_chr(YLL049, recode)) %>%
  mutate(YLR = purrr::map_chr(YLR, recode)) %>%
  mutate(YPL009C = purrr::map_chr(YPL009C, recode))

df.recoded = data.frame(df.recoded[,-1], row.names = df.recoded$Strain)
tetraploids = df2genind(df.recoded, sep = "/", ploidy = 4)

# set strata 
strata(tetraploids) = tetra.strata
strata(tetraploids) %>% head(4)
setPop(tetraploids) = ~Ploidy

tetraploids

# clean
rm(df, df.recoded, tetra.strata)
```

This results from a new dataset containing `r adegenet::nInd(tetraploids)` tetraploids genotypes.

+ New encoding

```{r}
genind2df(tetraploids, sep = "/") %>% head

# see summary
# Count ploidy per Origin
strata(tetraploids) %>% 
  group_by(Ploidy, Isolation) %>% summarise(Count = n()) %>%
  tidyr::spread(key = Ploidy, value = Count, fill = 0)

# Filter only commercial and sourdough
setPop(tetraploids) = ~Isolation
tetraploids = popsub(tetraploids, sublist = c("commercial", "sourdough"))

# see summary
strata(tetraploids) %>% 
  group_by(Ploidy, Isolation) %>% summarise(Count = n()) %>%
  tidyr::spread(key = Ploidy, value = Count, fill = 0)

# export to genAlEx formatted *.csv file
# genind2genalex(
#   tetraploids, 
#   filename = "data.tetraploid.genalex.csv", 
#   pop = ~Isolation/Country/Bakery, 
#   allstrata = FALSE, 
#   overwrite = TRUE)

#poppr(tetraploids, strata = ~Isolation)
```


## Merging 2N+4N genotypes

```{r}
# merge 2n + 4n
df.2n = genind2df(diploids, sep = "/")[,-1]
df.4n = genind2df(tetraploids, sep = "/")[,-1]
df.total = rbind(df.2n, df.4n)

# generate total genotype file
total = df2genind(df.total, sep = "/")
ploidy(total) = c(ploidy(diploids), ploidy(tetraploids))
strata(total) = rbind(strata(diploids), strata(tetraploids))

# clean
rm(df.2n, df.4n, df.total)
```

This results from a new dataset containing `r adegenet::nInd(total)` genotypes.


## Summary of the new dataset

```{r}
# Count ploidy per Origin
counts = strata(total) %>% 
  group_by(Ploidy, Isolation) %>% summarise(Count = n()) %>%
  tidyr::spread(key = Ploidy, value = Count, fill = 0)

counts
counts = counts[,c(2,3,4)] %>% as.matrix()

# Fisher's exact test for contingency table
fisher.test(counts)
rm(counts)

# poppr
setPop(total) = ~Isolation
poppr(total, sample = 999, plot = TRUE)

# poppr
#setPop(total) = ~Ploidy
#poppr(total)
```


```{r}
# Filter only commercial and sourdough
setPop(total) = ~Isolation/Country
sourdough.europe = popsub(total, sublist = c("sourdough_France", "sourdough_Italy", "sourdough_Belgium"))

# Count ploidy per Origin
counts = strata(sourdough.europe) %>% 
  group_by(Ploidy, Isolation, Country) %>% summarise(Count = n()) %>%
  tidyr::spread(key = Ploidy, value = Count, fill = 0)

counts
counts = counts[,c(3,4,5)] %>% as.matrix()

# Fisher's exact test for contingency table
fisher.test(counts)
  rm(counts)

# poppr
poppr(sourdough.europe, sample = 999, plot = TRUE)

# clean
rm(sourdough.europe)
```


# RESULTS: analysis of the yeasts dataset

## Analysis of Molecular Variance (AMOVA) 

### Analysis of *sourdough* yeasts from *France*

Analysis of Molecular Variance (AMOVA) is a method of estimating population differentiation directly from molecular data and testing hypotheses about such differentiation. In panmictic populations, we would expect to see most of the variance arise from within samples. On the contrary, if we see that the most of the variance occurs among samples within populations or among populations, then there is evidence that we have some sort of population structure. In the case of clonal organisms, this would help support a hypothesis of clonal reproduction. 

Let's analyse dataset structuration using AMOVA and the hierarchical model `~Sourdough`. We restrict this analysis on yeasts isolated from *France*.

```{r}
# select yeasts isolated from France
setPop(total) = ~Isolation
sourdough = popsub(total, sublist = c("sourdough"))
setPop(sourdough) = ~Country
sourdough.france = popsub(sourdough, sublist = c("France"))

# compute AMOVA
amova.sourdough.france = poppr.amova(sourdough.france, hier = ~Sourdough, cutoff = .20, freq = TRUE)
amova.sourdough.france
```

The `$componentsofcovariance` table shows how much variance is detected at each stratification. We expect variations within samples to give the greatest amount of variation for populations that are not significantly differentiated. 

To test if populations are significantly different, we perform a randomization test using the function `randtest()` from the *ade4* package. This will randomly permute the sample matrices as described in (Excoffier et al., 1992).

```{r fig.asp=.8}
library(ade4)
set.seed(5893)
amova.france.sourdough.test = randtest(amova.sourdough.france, nrepet = 1000)
plot(amova.france.sourdough.test)
amova.france.sourdough.test

#clean
rm(sourdough, sourdough.france, amova.sourdough.france, amova.france.sourdough.test)
```

### Analysis of total *sourdough* yeasts from Europe by Country

```{r}
# select sourdough yeasts
setPop(total) = ~Isolation
sourdough.yeast = popsub(total, sublist = c("sourdough"))

# select European yeasts
setPop(sourdough.yeast) = ~Country
sourdough.europe = popsub(sourdough.yeast, sublist = c("France", "Italy", "Belgium"))

# compute AMOVA
amova.sourdough.europe = poppr.amova(sourdough.europe, hier = ~Country, cutoff = .20, freq = TRUE)
amova.sourdough.europe

library(ade4)
set.seed(5893)
amova.sourdough.europe.test = randtest(amova.sourdough.europe, nrepet = 1000)
plot(amova.sourdough.europe.test)
amova.sourdough.europe.test

# clean
rm(amova.sourdough.europe, amova.sourdough.europe.test)
```

### Analysis of sampled *sourdough* yeasts from Europe by Country

```{r}
# get sampled yeasts
setPop(total) = ~Name
sampled.yeasts = popsub(total, sublist = sampled.yeast.name)

# select sourdough yeasts
setPop(sampled.yeasts) = ~Isolation
sampled.sourdough = popsub(sampled.yeasts, sublist = c("sourdough"))

# select European yeasts
setPop(sampled.sourdough) = ~Country
sampled.sourdough.europe = popsub(sampled.sourdough, sublist = c("France", "Italy", "Belgium"))

# compute AMOVA
amova.sourdough.europe = poppr.amova(sampled.sourdough.europe, hier = ~Country, cutoff = .20, freq = TRUE)
amova.sourdough.europe

library(ade4)
set.seed(5893)
amova.sourdough.europe.test = randtest(amova.sourdough.europe, nrepet = 1000)
plot(amova.sourdough.europe.test)
amova.sourdough.europe.test

# clean
rm(sampled.yeasts, amova.sourdough.europe, amova.sourdough.europe.test)
```

### Analysis of total *sourdough* yeasts vs *commercial* by Isolation

```{r}
# AMOVA
amova.total.yeasts = poppr.amova(total, hier = ~Isolation, cutoff = .20, freq = TRUE)
amova.total.yeasts

# test
library(ade4)
#set.seed(5893)
amova.total.yeasts.test = randtest(amova.total.yeasts, nrepet = 1000)
plot(amova.total.yeasts.test)
amova.total.yeasts.test

# clean
rm(amova.total.yeasts, amova.total.yeasts.test)
```

### Analysis of *sourdough* yeasts vs *commercial* by Isolation/Ploidy

```{r}
# AMOVA
amova.total.yeasts = poppr.amova(total, hier = ~Isolation/Ploidy, cutoff = .20, freq = TRUE)
amova.total.yeasts

# test
library(ade4)
#set.seed(5893)
amova.total.yeasts.test = randtest(amova.total.yeasts, nrepet = 1000)
plot(amova.total.yeasts.test)
amova.total.yeasts.test

# clean
rm(amova.total.yeasts, amova.total.yeasts.test)
```


### Analysis of sampled *sourdough* yeasts vs *commercial* by Isolation

```{r}
# samples yeasts
setPop(total) = ~Name
sampled.yeasts = popsub(total, sublist = sampled.yeast.name)

# AMOVA
amova.sampled.yeasts = poppr.amova(sampled.yeasts, hier = ~Isolation, cutoff = .20, freq = TRUE)
amova.sampled.yeasts

# test
library(ade4)
#set.seed(5893)
amova.sampled.yeasts.test = randtest(amova.sampled.yeasts, nrepet = 1000)
plot(amova.sampled.yeasts.test)
amova.sampled.yeasts.test

# clean
rm(amova.sampled.yeasts, amova.sampled.yeasts.test)
```


## NJ tree using Bruvo's genetic distance

### Total yeasts set

```{r total.tree, fig.width=10, fig.asp=1.5}
# compute distance matrix using Bruvo's distance
bruvo.dist.matrix = bruvo.dist(total, replen = loci.nucl.repeat, add = TRUE, loss = TRUE)

# compute NJ tree (use package ape)
library(ape)
nj.data = nj(bruvo.dist.matrix) %>% ladderize()
plot(nj.data, cex = 0.5)

# save NJ tree
write.tree(nj.data, file = "total.yeasts.bruvo.nj-tree.nwk")
```

### Diploid yeasts set

```{r diplo.tree, fig.width=10, fig.asp=1.5}
# compute distance matrix using Bruvo's distance
diplo.bruvo.dist.matrix = bruvo.dist(diploids, replen = loci.nucl.repeat, add = TRUE, loss = TRUE)

# compute NJ tree (use package ape)
library(ape)
nj.data = nj(diplo.bruvo.dist.matrix) %>% ladderize()
plot(nj.data, cex = 0.5)

# save NJ tree
write.tree(nj.data, file = "diploid.yeasts.bruvo.nj-tree.nwk")
```

### Tetraploid yeasts set

```{r tetra.tree, fig.width=10, fig.asp=1.5}
# compute distance matrix using Bruvo's distance
tetra.bruvo.dist.matrix = bruvo.dist(tetraploids, replen = loci.nucl.repeat, add = TRUE, loss = TRUE)

# compute NJ tree (use package ape)
library(ape)
nj.data = nj(tetra.bruvo.dist.matrix) %>% ladderize()
plot(nj.data, cex = 0.5)

# save NJ tree
write.tree(nj.data, file = "tetraploid.yeasts.bruvo.nj-tree.nwk")
rm(diplo.bruvo.dist.matrix, tetra.bruvo.dist.matrix, nj.data)
```


## Analysis of the index of association $I_a$

Now we test if populations are in linkage disequilibrium or not. This test is useful to determine if populations are clonal (where significant disequilibrium is expected due to linkage among loci) or sexual (where linkage among loci is not expected). The null hypothesis tested is that alleles observed at different loci are not linked if populations are sexual while alleles recombine freely into new genotypes during the process of sexual reproduction. In molecular ecology we typically use the index of association or related indices to test this phenomenon.

+ Here is the test for the *sourdough* yeasts:

```{r}
# compute index of association
ia(sampled.sourdough.europe, sample = 999)
```

Let’s look at the clone-corrected data and make sure this is not an artifact of clonality:

```{r}
# correct for clonality
sampled.sourdough.europe %>% clonecorrect(strata= ~Sourdough) %>% ia(sample=999)

# calculate IA and rd over all pairs of loci
sourdough.pair <- sampled.sourdough.europe %>% clonecorrect(strata = ~Sourdough) %>% pair.ia
rm(sourdough.pair)
```

Here we can reject the hypothesis that alleles are not linked across loci with $p=0.001$. Thus *sourdough* yeasts seem to be clonal.

+ Here for commercial yeasts

```{r}
setPop(total) = ~Isolation
commercial.yeasts = popsub(total, sublist = "commercial")

# compute index of association
ia(commercial.yeasts, sample = 999)

# with clone correction
commercial.yeasts %>% clonecorrect() %>% ia(sample=999)
```

The same for commercial yeasts that seemed to be clonal.


## Analysis of diploid yeasts genetic structuration using DAPC

The Discriminant Analysis of Principal Components (DAPC) is designed to investigate the genetic structure of biological populations. This multivariate method consists in a two-steps procedure. First, genetic data are transformed (centred, possibly scaled) and submitted to a Principal Component Analysis (PCA). Second, principal components of PCA are submitted to a Linear Discriminant Analysis (LDA). DAPC in itself requires prior groups to be defined. 

We need to remove clones and filter the dataset:


### Clone correction and dataset filtration

+ Clone correction

Since we want to only observe unique genotypes per population, clone correction requires specification of the stratifications at which clones should be censored. 

We go about correcting for the clones:

```{r}
# censor clones with respect to Origin
clone.corrected = clonecorrect(pop = total, strata = ~Sourdough, combine = TRUE)
poppr(clone.corrected)
```

It is important to analyze each locus independently to make sure that one locus is not introducing bias or spurious errors into the analysis.


+ Filtering uninformative loci

```{r}
# removing uninformative loci
cutoff.diff = 2/nInd(clone.corrected)
MAF = 0.01
```

We try to remove loci that contain less than a given percentage of divergent individuals (min differentiating samples: `r cutoff.diff`; min MAF: `r MAF`):

```{r}
inform = informloci(clone.corrected, quiet = FALSE, cutoff = cutoff.diff, MAF = MAF)
locNames(inform)
cat("Number of loci: ", nLoc(inform), sep = "")
```

All loci seem informative and none were removed using default setings.


+ Filtering loci with high level of missing data

```{r}
cutoff.loci = .2
```

Here we try to remove loci specifying a cutoff (max `r cutoff.loci`) representing 
the percent of missing genotype allowed:

```{r}
# filtering missing data for loci, max 20% of missing genotypes
filtered.loci = missingno(pop = inform,  type = "loci", cutoff = cutoff.loci, quiet = FALSE)
locNames(filtered.loci)
n.loci = nLoc(filtered.loci)
cat("Number of loci:", n.loci)
```

After filtering, `r n.loci` loci were retained.


+ Filtering yeasts with high level of missing genotypes

```{r}
# for genotypes, max 4 loci not genotyped
min.genotype = 4
cutoff.genotype = min.genotype/nLoc(filtered.loci)
```

Here we try to remove individuals specifying a cutoff (max `r round(cutoff.genotype, 2)`) representing the percent of missing genotypes allowed for an individual.

```{r}
# filter
total.filtered = missingno(pop = filtered.loci, type = "geno", cutoff = cutoff.genotype)

# summary
n.individuals = nInd(total.filtered)
cat("Number of individuals:", n.individuals)
```

Filtering removes `r adegenet::nInd(filtered.loci) - n.individuals` individuals, resulting in a dataset 
containing `r n.individuals` individuals. The shape of the rarefaction curve above doen't changes significantly.

```{r}
# accumulation curve
genotype_curve(total.filtered, sample = 1000, quiet = TRUE)

# set repeat length for the remaining loci
loci.nucl.repeat = loci.nucl.repeat[locNames(total.filtered),]

# clean
rm(cutoff.diff, MAF, cutoff.loci, n.loci, min.genotype, cutoff.genotype, n.individuals)
rm(clone.corrected, inform, filtered.loci)
```

Summary statistics

```{r}
# population structuration
setPop(total.filtered) = ~Isolation

# 
poppr(total.filtered)


# Filter only sourdough
setPop(total.filtered) = ~Isolation/Country
sourdough.europe = popsub(total.filtered, sublist = c("sourdough_France", "sourdough_Italy", "sourdough_Belgium"))

# poppr
poppr(sourdough.europe)
```

## Discrimination of the origin group: *Sourdough* vs *Commercial*

We first use population strata as prior groups:

```{r}
# population structuration
setPop(total.filtered) = ~Isolation

# vector of population colors
pop.colors = nPop(total.filtered) %>% funky()

# Compute DAPC
dapc.res = dapc(
    total.filtered, 
    pop = pop(total.filtered), 
    n.pca = 10, 
    n.da = nPop(total.filtered), 
    scale = FALSE, 
    var.contrib = TRUE, 
    var.loadings = FALSE, 
    pca.info = TRUE)

# barplot of eigenvalues
barplot(dapc.res$eig, main = "Eigenvalues", col = RColorBrewer::brewer.pal(10, name = "Set3")[6])

# discriminant 1/2
scatter.dapc(
    x = dapc.res,
    xax = 1,
    yax = 2,
    scree.da = FALSE,
    col = pop.colors,
    bg = "white",
    pch = 20,
    cell = 1.5,
    cstar = 1,
    solid = .6,
    cex = 2,
    clab = .6,
    legend = F,
    xlab = "Discriminant Function 1",
    ylab = "Discriminant Function 2")

rm(dapc.res, pop.colors)
```


## Discrimination of the geographical isolation

```{r}
# population structuration
setPop(total.filtered) = ~Country
#diplo.na.rm = popsub(diploids.cc.filter, blacklist = c("ND"))

# vector of population colors
pop.colors = nPop(total.filtered) %>% funky()

# Compute DAPC
dapc.res = dapc(
    total.filtered, 
    pop = pop(total.filtered), 
    n.pca = 10, 
    n.da = nPop(total.filtered), 
    scale = FALSE, 
    var.contrib = TRUE, 
    var.loadings = FALSE, 
    pca.info = TRUE)

# barplot of eigenvalues
barplot(dapc.res$eig, main = "Eigenvalues", col = RColorBrewer::brewer.pal(10, name = "Set3")[6])

# discriminant LD1/LD2
scatter.dapc(
    x = dapc.res,
    xax = 1,
    yax = 2,
    scree.da = FALSE,
    col = pop.colors,
    bg = "white",
    pch = 20,
    cell = 1.5,
    cstar = 1,
    solid = .6,
    cex = 2,
    clab = .6,
    legend = F,
    xlab = "Discriminant Function 1",
    ylab = "Discriminant Function 2")

# discriminant LD2/LD3
scatter.dapc(
    x = dapc.res,
    xax = 2,
    yax = 3,
    scree.da = FALSE,
    col = pop.colors,
    bg = "white",
    pch = 20,
    cell = 1.5,
    cstar = 1,
    solid = .6,
    cex = 2,
    clab = .6,
    legend = FALSE, 
    xlab = "Discriminant Function 2",
    ylab = "Discriminant Function 3")

rm(dapc.res, pop.colors)
```


## Cluster identification using the K-means

Here the identification of the clusters is achieved by `find.clusters`. This function first transforms the
data using PCA, asking the user to specify the number of retained PCs interactively unless
the argument n.pca is provided. Then, it runs k-means algorithm with increasing values of k, 
unless the argument `n.clust` is provided, and computes associated summary statistics (by default, BIC).

```{r}
clust.res = find.clusters(x = total.filtered, 
                          n.pca = 5,
                          max.n.clust = 20)
```

The lowest BIC values were obtained with $k$ between $6$ and $8$.


```{r dapc, fig.width=5}
k.values = 6:8

for (k in k.values) {

  print(paste("Discriminate analysis using k=", k, sep = ""))
  
  clust.res = find.clusters(x = total.filtered, n.pca = 5, n.clust = k)
  
  # compute DAPC
  dapc.res = dapc(
    total.filtered, 
    pop = clust.res$grp, 
    n.pca = 5, 
    n.da = 10, 
    scale = FALSE, 
    var.contrib = TRUE, 
    var.loadings = FALSE, 
    pca.info = TRUE
  )
  
  # barplot of eigenvalues
  barplot(dapc.res$eig, main = "Eigenvalues", col = RColorBrewer::brewer.pal(10, name = "Set3")[6])
  
  # vector of population colors
  pop.colors = funky(k)
  
  # discriminant LD1/LD2
  scatter(
    x = dapc.res,
    xax = 1,
    yax = 2,
    scree.da = FALSE,
    col = pop.colors,
    bg = "white",
    pch = 20,
    cell = 1.5,
    cstar = 1,
    solid = .6,
    cex = 2,
    clab = .6,
    legend = F,
    posi.pca = "bottomleft",
    posi.da = "bottomright",
    xlab = "Discriminant Function 1",
    ylab = "Discriminant Function 2"
  )
  
  # discriminant LD2/LD3
  scatter(
    x = dapc.res,
    xax = 2,
    yax = 3,
    scree.da = FALSE,
    col = pop.colors,
    bg = "white",
    pch = 20,
    cell = 1.5,
    cstar = 1,
    solid = .6,
    cex = 2,
    clab = .6,
    legend = FALSE, 
    xlab = "Discriminant Function 2",
    ylab = "Discriminant Function 3"
  )
  
  
  # exporting membership probability
  struct.export = as.data.frame(dapc.res$posterior)
  colnames(struct.export) = paste("P", 1:length(struct.export), sep = "")
  init.columns = colnames(struct.export)
  
  # create Group with missing values
  struct.export$Group = as.character(NA)
  
  # create Ind column
  struct.export$Ind = row.names(struct.export)
  rownames(struct.export) = NULL
  
  # order columns to export
  struct.export = struct.export[,c("Group","Ind", init.columns)]
  
  # output file
  out_fic = paste("microsat.data.adegenet.K=", k, ".csv", sep = "")
  write.csv(struct.export, file = out_fic, row.names = FALSE, quote = FALSE)
  
  # # Let’s see how the above groups found ith the clustering algorithm match up on a phylo tree:
  # tree = bruvo.boot(total.filtered, replen = loci.nucl.repeat, sample = 100, cutoff = 50, quiet = TRUE, showtree = FALSE)
  # plot.phylo(tree,
  #            cex = 0.6, font = 2, adj = 0, tip.color = pop.colors[clust.res$grp],
  #            label.offset = 0.0125)
  # nodelabels(tree$node.label, adj = c(1.3, -0.5), frame = "n", cex = 0.8,
  #            font = 3, xpd = TRUE)
  # axisPhylo(3)
  
  # Structure like plot:
  compoplot(
    dapc.res, 
    legend = TRUE, 
    posi = list(x=-0.25, y=1.3), 
    ncol = k, 
    col = pop.colors, 
    txt.leg = paste("Cluster", 1:k), 
    cleg = .7)
}

# clean
rm(init.columns, struct.export, out_fic, clust.res, dapc.res, tree)
```


```{r}
sessionInfo()
```
