---
title: 'Analysis of the genetic structuration of bakery yeasts using Adegenet'
author: '[Frédéric Bigey](frederic.bigey@inrae.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   = FALSE,
    message   = FALSE, 
    fig.align = "center", 
    fig.width = 4.0,
    fig.asp   = 1.0,
    fig.path  = 'Figs/',
    cache     = TRUE,
    size      = "tiny"
)
options(knitr.table.format = "latex")
```


# Materials & Methods #

## Materials

```{r echo=TRUE}
input.file = "Bakery.SNP.BIALLELIC.mind-0.5.geno-0.1.maf-0.03.ld-0.5.raw"
name_file="liste_Bakery_Fred_10Octobre2017.csv"
```


## Methods

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

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

# Analysis of the dataset

## Import the dataset 

```{r import, eval=TRUE}
# read in a Plink file, return a genlight object
data = read.PLINK(input.file)

# convert to snpclone object
data = as.snpclone(data)

# import strain informations
info = readr::read_csv(name_file)
sample.id = tibble("Sample ID" = indNames(data))
info = left_join(sample.id, info, by = "Sample ID")

# change strain name
indNames(data) = info[["Strain name"]]
indNames(data)
```

Import strata of the dataset:

```{r set.strata}
# set strata using the data frame
strata(data) = info[c("Isolation", "Type","Cluster")]
strata(data) %>% head(6)
```


# Principal Componant Analysis (PCA)

```{r pca}
# compute PCA
pca.res = glPca(data, nf = 50)
```

In the following plots, colors are based on the first 3 PCs giving the 3 RGB channels. Genetic diversity is now represented in two complementary ways: by distances and by colors.

```{r plot.pca}
# axis 1/2
myColor=colorplot(pca.res$scores[,c(1,2)], pca.res$scores, transp=TRUE, cex=2, xlab="PC 1", ylab="PC 2")
title("PCA of SNP dataset\naxes 1-2")
abline(v=0, h=0, col="grey", lty=2)

# axis 2/3
colorplot(pca.res$scores[,c(2,3)], pca.res$scores, transp=TRUE, cex=2, xlab="PC 2", ylab="PC 3")
title("PCA of SNP dataset\naxes 2-3")
abline(v=0, h=0, col="grey", lty=2)
```

## Neighbor-joining tree on pairwise distances between individuals

```{r plot.tree, fig.width=7, eval=TRUE}
library(ape)
dist.mat = dist(as.matrix(data))
tree = nj(dist.mat)
plot(tree, type = "fan", show.tip = TRUE, cex = .75)
tiplabels(pch = 20, col = myColor, cex = 3)
title("NJ tree of pairwise distances between individuals")
```


# 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

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 specifying the `Type` (Sourdough vs Commercial) strata:

```{r clone.correct}
# censor clones with respect to Origin
data.clone.cor = clonecorrect(pop = data, strata = ~Type, combine = TRUE)

# output filtered individuals
clone.rm = setdiff(x = data@ind.names, y = data.clone.cor@ind.names)
print(paste("Individuals removed:", length(clone.rm)))
clone.rm

# DO NOT clone correct!!!!
#data.clone.cor = data
```


## Discrimination using the `Type` group: *Sourdough*, *Dough*, *Commercial*

```{r dapc.1}
# population structurationn
setPop(data.clone.cor) = ~Type
pop.names = popNames(data.clone.cor)
print(paste("Number of population:", length(pop.names)))
pop.names

# vector of population colors
pop.colors = nPop(data.clone.cor) %>% funky()

# Determine number of PCA to keep using cross validate
set.seed(999)
dapc.xval = xvalDapc(
  tab(data.clone.cor), 
  pop(data.clone.cor), 
  n.pca = 3:20, 
  n.rep = 50, 
  parallel = "multicore", 
  ncpus = 4L)

dapc.xval$`Number of PCs Achieving Lowest MSE`
dapc.res = dapc.xval$DAPC

# The slot assign.per.pop indicates the proportions of successful reassignment
#  (based on the discriminant functions) of individuals to their original clusters. 
summary(dapc.res)

# how well they have been retrieved by the procedure
assign.tab = table(pop(data.clone.cor), dapc.res$assign)
assign.tab
#table.value(assign.tab, col.lab=paste("infer", 1:6), row.lab=paste0("origin", 1:6))

#This figure is the simple graphical translation of the posterior table above. Heat colors
#represent membership probabilities (red=1, white=0); blue crosses represent the prior
#cluster provided to DAPC
assignplot(dapc.res)

# Compute DAPC
# dapc.res = dapc(
#     data.clone.cor, 
#     pop = pop(data.clone.cor), 
#     n.pca = 50, 
#     n.da = nPop(data.clone.cor), 
#     scale = FALSE, 
#     var.contrib = TRUE, 
#     var.loadings = FALSE, 
#     pca.info = TRUE)

# barplot of eigenvalues
barplot(dapc.res$eig, main = "Eigenvalues")

# 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")
```


## Discrimination using the `Cluster` grouping (P1..P6)

```{r dapc.2}
# set population structurationn
setPop(data.clone.cor) = ~Cluster
pop.names = popNames(data.clone.cor)
print(paste("Number of population:", length(pop.names)))
pop.names

# vector of population colors
pop.colors = nPop(data.clone.cor) %>% funky()

# Determine number of PCA to keep using cross validate
set.seed(999)
dapc.xval = xvalDapc(
  tab(data.clone.cor), 
  pop(data.clone.cor), 
  n.pca = 3:20, 
  n.rep = 50, 
  parallel = "multicore", 
  ncpus = 4L)

dapc.xval$`Number of PCs Achieving Lowest MSE`
dapc.res = dapc.xval$DAPC

# The slot assign.per.pop indicates the proportions of successful reassignment
#  (based on the discriminant functions) of individuals to their original clusters. 
summary(dapc.res)

# how well they have been retrieved by the procedure
assign.tab = table(pop(data.clone.cor), dapc.res$assign)
assign.tab
#table.value(assign.tab, col.lab=paste("infer", 1:6), row.lab=paste0("origin", 1:6))

#This figure is the simple graphical translation of the posterior table above. Heat colors
#represent membership probabilities (red=1, white=0); blue crosses represent the prior
#cluster provided to DAPC
assignplot(dapc.res)

# barplot of eigenvalues
barplot(dapc.xval$DAPC$eig, main = "Eigenvalues")

# 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,
    clabel = .6,
    legend = F,
    xlab = "Discriminant Function 1",
    ylab = "Discriminant Function 2")

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

## 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, and computes associated summary statistics (by default, BIC).

```{r find.k}
# search the best k
#set.seed(999)
#clust.res = find.clusters(x = data.clone.cor, n.pca = 60, max.n.clust = 40)
```

The lowest BIC values were obtained with $k$ between $5$ and $7$. Test for:


### Discriminate analysis

```{r dapc.3}
# set population structuration
setPop(data.clone.cor) = ~Cluster
pop.names = popNames(data.clone.cor)
print(paste("Number of population:", length(pop.names)))
pop.names

# k-means clustering
k = 7
set.seed(999)
clust.res = find.clusters(x = data.clone.cor, n.pca = 60, n.clust = k)

# how well they have been retrieved by the procedure
assign.tab = table(pop(data.clone.cor), clust.res$grp)
assign.tab
table.value(assign.tab, col.lab=paste("infer", 1:k), row.lab=paste0("origin", 1:6))
    
# Determine number of PCA to keep using cross validate
set.seed(999)
dapc.xval = xvalDapc(
  tab(data.clone.cor), 
  clust.res$grp, 
  n.pca = 3:15, 
  n.rep = 50, 
  parallel = "multicore", 
  ncpus = 4L)

dapc.xval$`Number of PCs Achieving Lowest MSE`
dapc.res = dapc.xval$DAPC

# The slot assign.per.pop indicates the proportions of successful reassignment
#  (based on the discriminant functions) of individuals to their original clusters. 
summary(dapc.res)

# how well they have been retrieved by the procedure
assign.tab = table(pop(data.clone.cor), dapc.res$assign)
assign.tab
#table.value(assign.tab, col.lab=paste("infer", 1:k), row.lab=paste("origin", 1:6))

#This figure is the simple graphical translation of the posterior table above. Heat colors
#represent membership probabilities (red=1, white=0); blue crosses represent the prior
#cluster provided to DAPC
assignplot(dapc.res)

# barplot of eigenvalues
barplot(dapc.res$eig, main = "Eigenvalues")

# vector of population colors
pop.colors = funky(k)

# 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,
    clabel = F,
    legend = T,
    posi.pca = "bottomleft",
    posi.da = "bottomright",
    txt.leg=paste("Cluster",1:k),
    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,
    clabel = F,
    legend = T, 
    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("diploid.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 in the clustering algorithm match up on a phylo tree:

```{r tree.dapc, fig.width=10, fig.asp=1.5}
# compute NJ tree
dist.mat = dist(as.matrix(data.clone.cor))
tree = ape::nj(dist.mat)

# plot
ape::plot.phylo(
    tree, 
    type = "fan",
    cex = 0.8, 
    font = 2, 
    adj = 0, 
    tip.color = pop.colors[clust.res$grp],
    label.offset = 0.0125)
ape::nodelabels(
    tree$node.label, 
    adj = c(1.3, -0.5), 
    frame = "n", 
    cex = 0.8,
    font = 3, 
    xpd = TRUE)
ape::axisPhylo(3)
```

Structure like plot:

```{r fig.asp=.4, fig.width=8.0}
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)
```
