# Replicate = actually corresponds to the batch of queens (2 batches collected at different times), each queen has been considered as a biological replicate. 
# Treatment = experimental condition
# Sample = unique sample number (matches tube number at -80°C, numbering is not continuous)
# Columns with gene names contains the CT averaged from qPCR duplicates

library("data.table")
library("tidyverse")
library("ggplot2")
library("outliers")
library("dplyr")
library("lme4")
library("car")
library("multcomp") # LOAD AFTER GRUBS TEST

############################################
####   Import and variable definition   ####
############################################

data_raw <- fread('Fig2B_data.csv') 

# INDICATE HOUSEKEEPING GENE NAME
HKG <- c("EF1a")
# INDICATE TARGET GENE NAME AND EFFICIENCIES - USE THE EXACT SAME NAME AS IN THE CSV FILE!
TargetList <- c("Def1", "Def2", "Hym", "Cru")
Efficiencies <- c(2, 2, 1.99, 1.68)

GeneList <- c(HKG, TargetList)

#####################################
#### Creation of the dCT columns ####
#####################################

for (i in TargetList) {
  new_col_name <- paste0("dCT_", i)  
  data_raw[[new_col_name]] <- data_raw[[i]] - data_raw$EF1a  
}

data_raw -> data_dCT


#####################################
####  Grubb's test for outliers  ####
#####################################

# The goal here will be to remove outliers from the dCT dataset. 
# Outliers will be identified with a Grubbs test, which can only detect 1 outlier. 
# The test is then performed twice, once for each side of the distribution (high or low value outliers).
# IF THERE IS AN ERROR, CLOSE:REOPEN R AND RUN THIS PART OF THE SCRIPT BEFORE LOADING THE MULTCOMP LIBRARY

#Remove raw CT columns
data_dCT_clean <- select(data_dCT, -starts_with(GeneList))

#Create the dCT name vector
dCT_nameList <- colnames(select(data_dCT_clean, starts_with("dCT")))

#Use the dCT name vector to pivot the table longer
data_dCT_clean_long <- data_dCT_clean %>%
  pivot_longer(cols = all_of(dCT_nameList), names_to = "Gene", values_to = "dCT")

# extract dCT values for each Treatment*Gene combination (Grubb's works on vectors)

treatments = unique(data_dCT_clean_long$Treatment)
genes = unique(data_dCT_clean_long$Gene)

# creates a named list where "Treatment_Gene" contains the dT value for each combination
res = list()
for(my_treatment in treatments){
  for(my_gene in genes){
    res_tmp = data_dCT_clean_long %>%
      filter(Treatment == my_treatment, Gene == my_gene) %>% pull (dCT)
    res[[paste0(my_treatment, "_", my_gene)]] = res_tmp
  }
}


# Then run the Grubb's test on each sublist of res. 
res_Grubbs = data.frame(condition = character(0),
                        p_min = numeric(0),
                        p_max = numeric(0))
for(i in names(res)){
  data_i = res[[i]]
  grubbs_i_max = outliers::grubbs.test(data_i, type = 10, opposite = FALSE)
  grubbs_i_min = outliers::grubbs.test(data_i, type = 10, opposite = TRUE)
  
  
  res_Grubbs_i =  data.frame(condition = i,
                             p_min = grubbs_i_min$p.value,
                             p_max = grubbs_i_max$p.value)
  
  res_Grubbs = rbind(res_Grubbs, res_Grubbs_i)
  
}

# To know how many outliers will be removed: 
# Low-value outliers: 
sum(res_Grubbs$p_min<0.05)
# High-value outliers : 
sum(res_Grubbs$p_max<0.05)


# now to remove outliers: 

data_dCT_clean_long_noOutliers <- data_dCT_clean_long 

for(j in 1:nrow(res_Grubbs)){
  condition_j = res_Grubbs[j, "condition"]
  treatment_j = stringr::str_split(condition_j, "_", simplify=T)[1]
  dCT_j = stringr::str_split(condition_j, "_", simplify=T)[2]
  genetemp_j = stringr::str_split(condition_j, "_", simplify=T)[3]
  gene_j = paste0(dCT_j, "_", genetemp_j)
  
  # if maximum value is an outlier, remove it
  if(res_Grubbs[j, "p_max"] < 0.05){
    vmax_j <- max(res[[condition_j]])
    data_dCT_clean_long_noOutliers %>%
      slice(-which(Treatment == treatment_j & Gene == gene_j & dCT == vmax_j)) -> data_dCT_clean_long_noOutliers
  }
  
  # if minimum value is an outlier, remove it
  if(res_Grubbs[j, "p_min"] < 0.05){
    vmin_j <- min(res[[condition_j]])
    data_dCT_clean_long_noOutliers %>%
      slice(-which(Treatment == treatment_j & Gene == gene_j & dCT == vmin_j)) -> data_dCT_clean_long_noOutliers
  }
}


#####################################
####   Calculation of the ddCT   ####  
##################################### 

# NOTE: If you want to keep the outliers, do not use data_dCT_clean_long_noOutliers as input but instead use data_dCT_clean_long after running the following:  
# data_dCT_clean_long <- data_dCT_clean %>%
#   pivot_longer(cols = all_of(dCT_nameList), names_to = "Gene", values_to = "dCT")

# Identification of the replicate #: 
# Replicates are not accounted for in the queen dataset as each queen is technically a standalone biological replicate. 
# Therefore the ddCT is calculated using all unchallenged samples as a reference. 

# Summarizing the average refence dCT of control treatment for each replicate): 

# DETERMINATION OF tHE CONTROL CONDITION TO WHICH CALCULATE DE ddCT: 
Control <- "Unchallenged"

data_dCT_clean_long_noOutliers %>%
  filter(Treatment == Control)  %>%
  group_by(Gene) %>%
  summarize(mean(dCT, na.rm = T)) -> df_mean_control_dCT

colnames(df_mean_control_dCT)[2] <- c("mean_control_dCT")

# Then create the ddCT columns: 

# First, add an empty column to the dataframe (NA_real_ populates the column with NAs)
data_dCT_clean_long_noOutliers %>%
  mutate(mean_control_dCT = NA_real_) -> data_dCT_clean_long_noOutliers

# Now use this column to paste the correct average control dCT in front of each dCT: 
for (my_gene in dCT_nameList) {
    mean_control_dCT_i <- df_mean_control_dCT %>%
      filter(Gene == my_gene) %>% pull(mean_control_dCT)
    # Note: case_when has a weird syntax: first part is the criterion for the 'case'
    # (here Replicate == my_rep & Gene == my_gene) and after the ~ what the loop must do when criterion is met. 
    # If the criterion is not met, the loop goes on to the second part (the alternative) and the TRUE here means 
    # "in all cases when the criterion is not met" = all the rest of the column, 
    # then the ~ indicates what to do when criterion not met, here, "paste the already existing value in the colum".
    data_dCT_clean_long_noOutliers %>%
      mutate(mean_control_dCT = case_when(
        Gene == my_gene ~ mean_control_dCT_i,
        TRUE ~ mean_control_dCT)) -> data_dCT_clean_long_noOutliers
}

data_dCT_clean_long_noOutliers -> data_noOutliers_dCT_FC_ddCT_long

# Now there are two columns: the dCT of the target, and the reference control dCT to be substracted to that of the target. 
# Simple substraction: 

data_dCT_clean_long_noOutliers$ddCT <- data_dCT_clean_long_noOutliers$dCT - data_dCT_clean_long_noOutliers$mean_control_dCT

data_dCT_clean_long_noOutliers %>%
  select(-mean_control_dCT) -> data_dCT_ddCT_clean_long_noOutliers


############################################
#### Calculation of the fold-changes FC ####
############################################

#Creation of a df containing the gene list and respective efficiency of the primers
d_FC = data.frame(TargetList, dCT_nameList, Efficiencies) 

data_dCT_ddCT_FC <- data_dCT_ddCT_clean_long_noOutliers

data_dCT_ddCT_FC %>% 
  pivot_wider(names_from = "Gene", values_from = c("ddCT", "dCT")) -> data_dCT_ddCT_FC_wide

# tidy the colnames: 
new_col_names <- colnames(data_dCT_ddCT_FC_wide) %>%
  str_replace("_dCT_", "_")
colnames(data_dCT_ddCT_FC_wide) <- new_col_names

#Now calculate the fold-change from the dCT, corrected by the primer efficiency: 
for (i in TargetList) {
  #FC column is created (empty)
  new_col_name <- paste0("FC_", i)
  #efficiency is collected from d_FC
  eff_i = d_FC %>%
    filter(TargetList==i) %>% pull(Efficiencies)
  
  #Column to use as input is identified: 
  dCT_name <- paste0("dCT_", i)  

  #Selection of the values matching the gene: 
    gene_i = data_dCT_ddCT_FC_wide %>%
    pull(dCT_name)
  
  data_dCT_ddCT_FC_wide[[new_col_name]] <- eff_i^-(gene_i)
  
}

#####################################
####   Calculation of the ddFC   ####
#####################################

# Same principle as for the FC, accounting for the efficiency of each primer set. 

#Reusing d_FC created at the FC calculation step: 
#d_FC = data.frame(TargetList, dCT_nameList, Efficiencies) 

# Apply a loop similar to that of the FC calculation: 

for (i in TargetList) {
  #ddFC column is created (empty)
  new_col_name <- paste0("ddFC_", i)
  #efficiency is collected from d_FC
  eff_i = d_FC %>%
    filter(TargetList==i) %>% pull(Efficiencies)
  #Column to use as input is identified: 
  ddCT_name <- paste0("ddCT_", i)  
  
  #Selection of the values matching the gene: 
    gene_i = data_dCT_ddCT_FC_wide %>%
    pull(ddCT_name)
  
    data_dCT_ddCT_FC_wide[[new_col_name]] <- eff_i^-(gene_i)
  
}

data_dCT_ddCT_FC_wide -> data_dCT_ddCT_FC_ddFC_wide

# For the long format, specify the column range to pivot longer (here all columns that contain the name of a gene in GeneList: 
selected_cols_to_pivot <- names(data_dCT_ddCT_FC_wide)[grepl(paste(GeneList, collapse = "|"), names(data_dCT_ddCT_FC_wide), fixed = FALSE)]
# Then pivot with a kind of string split using the "_" in each colname (eg. "dCT_Def" will pivot on "dCT" and on "Def")
data_dCT_ddCT_FC_ddFC_long <- pivot_longer(data_dCT_ddCT_FC_wide, cols = selected_cols_to_pivot, 
                     names_sep = "\\_", names_to = c("stat", "gene"))


#################### 
##### PLOTTING #####
####################

df_plotting <- data_dCT_ddCT_FC_ddFC_long


df_plotting %>%
  mutate(gene = tolower(as.character(gene))) -> df_plotting_ready

#Reorder treatments so that UC comes first: 
df_plotting_ready$Treatment <- factor(df_plotting_ready$Treatment, levels = c("Unchallenged", "Sham", "Ecc15", "Mluteus", "Conidiospores"))
#Reorder genes for the facet:
df_plotting_ready$gene <- factor(df_plotting_ready$gene, levels = c("hym", "def1", "def2", "cru"))


######################### 
###### PLOT BY GENE #####
######################### 

df_plotting_ready %>%
  # Exclude PPO
  filter(stat %in% "ddFC") %>%
  ggplot() +
  aes(x = Treatment, y = value, fill = Treatment) +
  geom_boxplot(shape = "circle", outlier.shape = NA) +
  scale_fill_brewer(palette = "BuGn", direction = -1) +
  theme_bw() +
  facet_wrap(vars(gene), scales = "free", ncol = 4) +
  geom_point(position = position_jitter(width = 0.1), size = 0.7, shape = 16) +
  #Add a 20% blank space at the top of each graph
  scale_y_continuous(limits = function(x) c(0, max(x) * 1.2)) +
  guides(fill = "none") +
  labs(y = "Fold-change relative to Unchallenged", x = "") +
  theme(axis.text.x = element_text(angle = 50, hjust = 1)) +
  scale_x_discrete(labels = c("Unchallenged", "Sham", expression(paste(italic("Ecc15"))), 
                              expression(paste(italic("M. luteus"))), expression(paste(italic("M. brunneum"))))) +
  #mettre les titres des facets en italic: 
  theme(strip.text = element_text(face = "italic"))


#############################
###### PLOT BY PATHOGEN #####
#############################

custom_labels <- c("Sham" = "Sham",
                   "Ecc15" = "Ecc15",
                   "Mluteus" = "M. luteus",
                   "Conidiospores" = "M. brunneum")


# Create a data frame for the annotation
# ann_text1 <- data.frame(gene = "def", value = 7900, label = "***", Treatment = "Sham")
# ann_text2 <- data.frame(gene = "hym", value = 7900, label = "***", Treatment = "Sham")
# ann_text3 <- data.frame(gene = "hym", value = 7900, label = "***", Treatment = "Mluteus")
# ann_text4 <- data.frame(gene = "hym", value = 7900, label = "***", Treatment = "Ecc15")

df_plotting_ready %>%
  filter(stat %in% "ddFC", gene %in% c("hym", "def1", "def2", "cru"), !Treatment %in% "Unchallenged") %>%
  ggplot() +
  aes(x = gene, y = value, fill = gene) +
  geom_boxplot(shape = "circle", outlier.shape = NA) +
  scale_fill_brewer(palette = "BuGn", direction = -1) +
  theme_bw() +
  facet_wrap(vars(Treatment), scales = "fixed", ncol = 4, labeller = labeller(Treatment = custom_labels)) +
  geom_point(position = position_jitter(width = 0.1), size = 0.7, shape = 16) +
  # scale_y_continuous(limits = function(x) c(0, max(x) * 1.1)) +
  coord_cartesian(ylim = c(0, 18000)) +
  guides(fill = "none") +
  labs(y = "Fold-change relative to Unchallenged", x = "") +
  theme(axis.text.x = element_text(angle = 50, hjust = 1, face = "italic")) 



#################
##### STATS #####
#################

df_plotting %>%
  mutate(gene = tolower(as.character(gene))) -> df_stats


# HYMENOPTAECIN
df_stats %>%
  filter(stat %in% "ddCT", gene == "hym") -> df_hym

# DEFENSIN1
df_stats %>%
  filter(stat %in% "ddCT", gene == "def1") -> df_def1

# DEFENSIN2
df_stats %>%
  filter(stat %in% "ddCT", gene == "def2") -> df_def2

# CRUSTIN
df_stats %>%
  filter(stat %in% "ddCT", gene == "cru") -> df_cru



####### CONTRAST MATRIX #########
contrast_matrix <- rbind (
  "Sham - Unchallenged"=c(0,0,0,1,-1),
  "Conidiospores - Sham"=c(0,0,0,-1,0),
  "Ecc15 - Sham"=c(0,1,0,-1,0),
  "Mluteus - Sham"=c(0,0,1,-1,0),
  "Ecc15 - Mluteus"=c(0,1,-1,0,0)
)

#HYM

modelhym <- lmer(value ~ Treatment + (1|Replicate), data = df_hym)

residualshym <- residuals(modelhym)
shapiro.test(residualshym)

summary(modelhym)

Anova(modelhym, type="II") 

summary(glht (modelhym , linfct = contrast_matrix), test=adjusted("BH") )


#DEF1
modeldef <- lmer(value ~ Treatment + (1|Replicate), data = df_def1)

residualsdef <- residuals(modeldef)
shapiro.test(residualsdef)

summary(modeldef)

Anova(modeldef, type="II") 

summary(glht (modeldef , linfct = contrast_matrix), test=adjusted("BH") )


#DEF2
modeldef2 <- lmer(value ~ Treatment + (1|Replicate), data = df_def2)

residualsdef2 <- residuals(modeldef2)
shapiro.test(residualsdef2)

summary(modeldef2)

Anova(modeldef2, type="II") 

summary(glht (modeldef2 , linfct = contrast_matrix), test=adjusted("BH") )

#CRU

modelcru <- lmer(value ~ Treatment + (1|Replicate), data = df_cru)

residualscru <- residuals(modelcru)
shapiro.test(residualscru)

summary(modelcru)

Anova(modelcru, type="II")