# 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 GRUBBS TEST

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

data_raw <- fread('Fig3_data.csv') 
data_raw$unique_number <- c(1:nrow(data_raw))

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

# The dataframe MUST contain a column "Replicate" containing the replicate or colony ID, and a column "Treatment" (sham, exposed, etc.). 
# It can contain other columns. 

#####################################
#### 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   #### MODIF FOR QUEENS -> NOT ACCOUNTING FOR REPLICATE 
##################################### BC EACH SAMPLE IS A BIOLOGICAL REPLICATE

# 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 <- c(unique(data_dCT_clean_long_noOutliers$Replicate))

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

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

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

data_dCT_ddCT_FC_ddFC_long -> df_plotting

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

# # /!\ first plotting revealed a completely absurd point, unique_number 9 for hym and def, and an absurd value for unique_number 28 for Hym specifically.
# # These points are removed prior to plotting and statistical analysis. 

df_plotting$value[df_plotting$unique_number == 9] <- NA
df_plotting$value[df_plotting$gene == "Hym" & df_plotting$Treatment == "male" & df_plotting$unique_number == 28] <- NA


# correct gene names : 
df_plotting %>%
  mutate(gene = tolower(as.character(gene))) -> df_plotting_ready



#Reorder treatments so that virgin comes first: 
df_plotting_ready$Treatment <- factor(df_plotting_ready$Treatment, levels = c("virgin", "claustral", "colonial", "male"))
#Reorder genes for the facet:
df_plotting_ready$gene <- factor(df_plotting_ready$gene, levels = c("hym", "def1", "def2", "cru", "lysc1", "PPO"))

################# 
###### PLOT #####
################# 

df_plotting_ready %>%
  # Exclude PPO
  filter(stat %in% "FC", !gene %in% c("PPO", "lysc1")) %>%
  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 = "Expression relative to EF1a", x = "") +
  theme(axis.text.x = element_text(angle = 50, hjust = 1)) +
  # scale_x_discrete(labels = c("Unchallenged", "Sham", expression(paste(italic("M. luteus"))), 
  # expression(paste(italic("Ecc15"))), expression(paste(italic("M. brunneum"))))) +
  #mettre les titres des facets en italic: 
  theme(strip.text = element_text(face = "italic")) +
  geom_vline(xintercept = c(3.5), linetype = "dotted", color = "black") 

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

# HYMENOPTAECIN
df_plotting_ready %>%
  filter(stat %in% "dCT", gene == "hym") -> df_hym
# DEFENSIN1
df_plotting_ready %>%
  filter(stat %in% "dCT", gene == "def1") -> df_def1
# DEFENSIN2
df_plotting_ready %>%
  filter(stat %in% "dCT", gene == "def2") -> df_def2
# CRUSTIN
df_plotting_ready %>%
  filter(stat %in% "dCT", gene == "cru") -> df_cru


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

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

summary(modelhym)

Anova(modelhym, type="II") 


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

residualsDef1 <- residuals(modelDef1)
shapiro.test(residualsDef1)

summary(modelDef1)

Anova(modelDef1, type="II") 

post_hocDef1 <- summary(glht(modelDef1, linfct=mcp(Treatment="Tukey")), test=adjusted("BH"))
print(cld(post_hocDef1))

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

residualsDef2 <- residuals(modelDef2)
shapiro.test(residualsDef2)

summary(modelDef2)

Anova(modelDef2, type="II") 

post_hocDef2 <- summary(glht(modelDef2, linfct=mcp(Treatment="Tukey")), test=adjusted("BH"))
print(cld(post_hocDef2))


#CRU
modelCru <- glmer((value) ~ Treatment + (1|Replicate), data=df_cru)

residualsCru <- residuals(modelCru)
shapiro.test(residualsCru)

summary(modelCru)

Anova(modelCru, type="II") 

post_hocCru <- summary(glht(modelCru, linfct=mcp(Treatment="Tukey")), test=adjusted("BH"))
print(cld(post_hocCru))