---
title: "Analyse effet ploidie sur la fermentation"
author: "F. Bigey"
date: "27/02/2020"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Importation des donnees et calcul paramètres cinetique fermentaire

```{r}
library(tidyverse)
library(readxl)
library(broom)
library(patchwork)
library(knitr)

#setwd("/home/bigey/Documents/Genomics/Bakery/Fermentation")

# Import cytometry data
all.data <- read_excel("../Data/2020-02-05_Results_4n_vs_2n.xlsx", na="na")
all.data$CO2max <- 0
all.data$seuil_1g <- 0
all.data$Vmax <- 0
all.data$tVmax <- 0
all.data$Date <- gsub("-","",all.data$Date)

# Import fermentation kinetic data (Bloc 1)
mydata1 <- read_csv2("../Data/2018-10-17 Levain 4n vs 2n Bloc1.csv") %>%
    select(
        code_fermentation = `code fermentation`, 
        time = `temps manip (heures)`, 
        cumulCO2, 
        debitCO2 = `debitCO2 (g/l/h)`
    )

mydata1$time <- as.numeric(mydata1$time)
mydata1$cumulCO2 <- as.numeric(mydata1$cumulCO2)
mydata1$debitCO2 <- as.numeric(mydata1$debitCO2)
mydata1$code_fermentation <- as.factor(mydata1$code_fermentation)

for(i in unique(mydata1$code_fermentation)) {
    tmp.data <- mydata1 %>% filter(code_fermentation == i)
    
    # Compute kinetic ferm. parameters
    seuil_1g <- tmp.data$time[tmp.data$cumulCO2 >= 1][1]
    CO2max <- max(tmp.data$cumulCO2)
    Vmax <- max(tmp.data$debitCO2)
    tVmax <- tmp.data$time[tmp.data$debitCO2 == Vmax]
    date = str_split(i, "-", simplify = T)[1,2]
    ID = as.numeric(str_split(i, "-", simplify = T)[1,3])
    
    if(nrow(all.data[all.data$Bloc == 1 & all.data$Date == date & all.data$Position == ID,]) != 0){
        all.data$seuil_1g[all.data$Bloc == 1 & all.data$Date == date & all.data$Position == ID] <- seuil_1g
        all.data$CO2max[all.data$Bloc == 1 & all.data$Date == date & all.data$Position == ID] <- CO2max
        all.data$Vmax[all.data$Bloc == 1 & all.data$Date == date & all.data$Position == ID] <- Vmax
        all.data$tVmax[all.data$Bloc == 1 & all.data$Date == date & all.data$Position == ID] <- tVmax
    }
}


# Import fermentation kinetic data (Bloc 2)
mydata2 <- read_csv2("../Data/2018-10-23 Levain 4n vs 2n Bloc2.csv") %>%
    select(
        code_fermentation = `code fermentation`, 
        time = `temps manip (heures)`,
        cumulCO2, 
        debitCO2 =  `debitCO2 (g/l/h)`
    )

mydata2$time <- as.numeric(mydata2$time)
mydata2$cumulCO2 <- as.numeric(mydata2$cumulCO2)
mydata2$debitCO2 <- as.numeric(mydata2$debitCO2)
mydata2$code_fermentation <- as.factor(mydata2$code_fermentation)

for(i in unique(mydata2$code_fermentation)){
    tmp.data <- mydata2 %>% filter(code_fermentation == i)
    
    seuil_1g <- tmp.data$time[tmp.data$cumulCO2 >= 1][1]
    CO2max <- max(tmp.data$cumulCO2)
    Vmax <- max(tmp.data$debitCO2)
    tVmax <- tmp.data$time[tmp.data$debitCO2 == Vmax]
    date = str_split(i, "-", simplify = T)[1,2]
    ID = as.numeric(str_split(i, "-", simplify = T)[1,3])
    
    if(nrow(all.data[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID,]) != 0){
        #print(i)
        all.data$seuil_1g[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID] <- seuil_1g
        #print(all.data$seuil_1g[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID])
        all.data$CO2max[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID] <- CO2max
        #print(all.data$CO2max[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID])
        all.data$Vmax[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID] <- Vmax
        #print(all.data$Vmax[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID] )
        all.data$tVmax[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID] <- tVmax
        #print(all.data$tVmax[all.data$Bloc == 2 & all.data$Date == date & all.data$Position == ID])
    }
}
```

## Analyse des fermentations en culture isolee (2n et 4n)

```{r}
all.data.isolee <- all.data %>%
    select(Souche_nonGFP,
           Souche_GFP,
           Ploidie,
           Type,
           celltotalt0,
           pur_t0nonGFP,
           pur_t0GFP,
           celltotalt24,
           pur_t24nonGFP,
           pur_t24GFP,
           CO2max,
           seuil_1g,
           Vmax,
           tVmax) %>%
    filter(Type == "isolee") %>%
    as.data.frame() %>%
    mutate(Souche = ifelse(is.na(Souche_nonGFP), Souche_GFP, Souche_nonGFP))# Generate plots

my.theme =  function() {
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))
}

g.CO2max <- all.data.isolee %>%
    ggplot(aes(x = reorder(Souche,CO2max), y = CO2max, fill = Ploidie)) +
    geom_boxplot() +
    labs(
        x = "Strains",
        y = expression(paste("max ", CO["2"], " (g)"))
    ) +
    my.theme()

g.seuil_1g <- all.data.isolee %>%
    ggplot(aes(x = reorder(Souche,seuil_1g), y = seuil_1g, fill = Ploidie)) +
    geom_boxplot() +
    labs(
        x = "Strains",
        y = expression(paste(t["1g"], " (h)"))
    ) +
    my.theme()

g.Vmax <- all.data.isolee %>%
    ggplot(aes(x = reorder(Souche,Vmax), y = Vmax, fill = Ploidie)) +
    geom_boxplot() +
    labs(
        x = "Strains",
        y = "Vmax (g/L/h)"
    ) +
    my.theme()

g.tVmax <- all.data.isolee %>%
    ggplot(aes(x = reorder(Souche,tVmax), y = tVmax, fill = Ploidie)) +
    geom_boxplot() +
    labs(
        x = "Strains",
        y = expression(paste(t[Vmax], " (h)"))
    ) +
    my.theme()

# Assemble plot (patchwork package)
cinet.param.plots = g.CO2max + g.seuil_1g + g.Vmax + g.tVmax + plot_annotation(tag_levels = 'A')
cinet.param.plots
ggsave(filename = "effet_ploidy_cinet_param.pdf")

# Compute ANOVA
param.kinetic.anova = rbind.data.frame(
    tidy(anova(lm(CO2max~Ploidie, data = all.data.isolee))),
    tidy(anova(lm(seuil_1g~Ploidie, data = all.data.isolee))),
    tidy(anova(lm(Vmax~Ploidie, data = all.data.isolee))),
    tidy(anova(lm(tVmax~Ploidie, data = all.data.isolee)))
) 

param.kinetic.anova %>% knitr::kable()
write_csv(param.kinetic.anova, path = "effet_ploidy_param_cinet.csv")
```

## Analyse des fermentations en competition : effet 2n *vs* 4n

```{r}
# Competition experiment => Type == "competition"
all.data.comp <- all.data %>%
    select(
        Souche_nonGFP,
        Souche_GFP,
        Ploidie,
        Type,
        celltotalt0,
        pur_t0nonGFP,
        pur_t0GFP,
        celltotalt24,
        pur_t24nonGFP,
        pur_t24GFP,
        CO2max,
        seuil_1g, 
        Vmax,
        tVmax
    ) %>%
    filter(Type == "competition", Ploidie %in% c("2n_4n", "4n_2n")) %>%  # remove "4n_4n" data
    mutate(couple = paste(Souche_nonGFP, "vs", Souche_GFP)) %>% # create strain couple
    mutate(lnGFP = log((pur_t24GFP*celltotalt24)/(pur_t0GFP*celltotalt0))) %>%
    mutate(lnnonGFP = log((pur_t24nonGFP*celltotalt24)/(pur_t0nonGFP*celltotalt0))) %>%
    mutate(W = if_else(Ploidie == "2n_4n", lnnonGFP/lnGFP, lnGFP/lnnonGFP))

# ANOVA 
#fitness.anova = anova(lm(W~Ploidie, data = all.data.comp)) %>% tidy()
#fitness.anova %>% kable()
#write_csv(fitness.anova, path = "effet_ploidy_sur_fitness.csv")

# t-test sur la totalité des valeurs 
# One Sample t-test two.sided

# Check t-test assumptions
ggplot(data = all.data.comp, aes(sample = W)) +
    geom_qq() + stat_qq_line()
shapiro.test(all.data.comp$W)

# t-test 
fitness.t.test = t.test(all.data.comp$W, mu = 1) %>% tidy()
fitness.t.test %>% kable()
write_csv(fitness.t.test, path = "effet_ploidy_fitness_global.csv")

# t-test par couple de souche
t.test.data = all.data.comp %>% 
    select(Souche_nonGFP, Souche_GFP, Ploidie, W) %>%
    group_by(Souche_nonGFP, Souche_GFP, Ploidie) %>% 
    nest()

t.test.data = t.test.data %>% mutate(results = map(data, function(df) t.test(df$W, mu = 1))) %>% 
    mutate(p.val = map_dbl(results, function(x) x[["p.value"]])) %>%
    mutate(df = map_dbl(results, function(x) x[["parameter"]])) %>%
    mutate(t = map_dbl(results, function(x) x[["statistic"]])) %>%
    mutate(mean = map_dbl(data, function(x) mean(x$W))) %>%
    select(Souche_nonGFP, Souche_GFP, Ploidie, df, t, p.val, mean) 
t.test.data %>% kable()
write_csv(t.test.data, path = "effet_ploidy_fitness_per_couple.csv")

# plot Strain vs w
couple.plot = all.data.comp %>%
    ggplot(aes(x = reorder(couple, W), y = W)) +
    geom_boxplot() +
    labs(
        x = "Couple of strains",
        y = "w, relative fitness"
    ) +
    my.theme()

couple.plot
ggsave(filename = "effet_ploidy_fitness.pdf", plot = couple.plot)
```

