library(tidyverse)
library(cowplot)
library(PMCMRplus)
data <- read.delim("Biochemical markers.txt")
data$Group <- as.factor(data$Group)
data$Group <- factor(data$Group, levels(data$Group)[c(4,2,3,1)])

#Descriptive statistics data frame
descriptive <- data.frame()
for(i in 1:4){
  descriptive[i,1] <- levels(data$Group)[i]
  for(j in 1:(ncol(data)-1)){
    descriptive[i,2*j] <- mean(data[data$Group == levels(data$Group)[i], j+1])
    names(descriptive)[2*j] <- names(data)[j+1]
    descriptive[i,2*j+1] <- sd(data[data$Group == levels(data$Group)[i], j+1])
    names(descriptive)[2*j+1] <- paste(names(data)[j+1], "sd", sep = "_")
  }
}
names(descriptive)[1] <- "Group"
descriptive$Group <- as.factor(descriptive$Group)
descriptive$Group <- factor(descriptive$Group, levels(descriptive$Group)[c(4,2,3,1)])

#Expression of Hmox1
RTqPCR <- read.delim("RT-qPCR.txt")
RTqPCR$Group <- as.factor(RTqPCR$Group)
RTqPCR$Group <- factor(RTqPCR$Group, levels(RTqPCR$Group)[c(4,2,3,1)])

for(i in 1:nrow(RTqPCR)){
  RTqPCR$CT_bact[i] <- mean(c(RTqPCR$CT1_bact[i], RTqPCR$CT2_bact[i]))
  RTqPCR$CT_Hmox1[i] <- mean(c(RTqPCR$CT1_Hmox1[i], RTqPCR$CT2_Hmox1[i]))
  RTqPCR$dCT[i] <- RTqPCR$CT_Hmox1[i] - RTqPCR$CT_bact[i]
}

for(i in 1:nrow(RTqPCR)){
  RTqPCR$ddCT[i] <- RTqPCR$dCT[i] - mean(RTqPCR[RTqPCR$Group == "SH", "dCT"])
  RTqPCR$expression[i] <- 2^-(RTqPCR$ddCT[i])
  RTqPCR$log2foldChange[i] <- log2(RTqPCR$expression[i])
}

for(i in 1:nrow(RTqPCR)){
  data$Hmox1[i] <- RTqPCR$log2foldChange[i]
}

Gene_expression <- data.frame()

for(i in 1:nrow(RTqPCR)){
  Gene_expression[i,1] <- data$Group[i]
  Gene_expression[i,2] <- RTqPCR$ddCT[i]*(-1)
}

#Vegfa and Pdk1 expression
RTqPCR2 <- read.delim("Vegfa_Pdk1.txt")
RTqPCR2$Group <- as.factor(RTqPCR2$Group)
RTqPCR2$Group <- factor(RTqPCR2$Group, levels(RTqPCR2$Group)[c(4,2,3,1)])

for(i in 1:nrow(RTqPCR2)){
  RTqPCR2$CT_bact[i] <- mean(c(RTqPCR2$CT1_bact[i], RTqPCR2$CT2_bact[i]), na.rm = TRUE)
  RTqPCR2$CT_Vegfa[i] <- mean(c(RTqPCR2$CT1_Vegfa[i], RTqPCR2$CT2_Vegfa[i]), na.rm = TRUE)
  RTqPCR2$CT_Pdk1[i] <- mean(c(RTqPCR2$CT1_Pdk1[i], RTqPCR2$CT2_Pdk1[i]), na.rm = TRUE)
  RTqPCR2$dCT_Vegfa[i] <- RTqPCR2$CT_Vegfa[i] - RTqPCR2$CT_bact[i]
  RTqPCR2$dCT_Pdk1[i] <- RTqPCR2$CT_Pdk1[i] - RTqPCR2$CT_bact[i]
}

for(i in 1:nrow(RTqPCR2)){
  RTqPCR2$ddCT_Vegfa[i] <- RTqPCR2$dCT_Vegfa[i] - mean(RTqPCR2[RTqPCR2$Group == "SH", "dCT_Vegfa"])
  RTqPCR2$ddCT_Pdk1[i] <- RTqPCR2$dCT_Pdk1[i] - mean(RTqPCR2[RTqPCR2$Group == "SH", "dCT_Pdk1"])
  RTqPCR2$Vegfa_expression[i] <- 2^-(RTqPCR2$ddCT_Vegfa[i])
  RTqPCR2$Pdk1_expression[i] <- 2^-(RTqPCR2$ddCT_Pdk1[i])
}

for(i in 1:nrow(RTqPCR)){
  RTqPCR2$Vegfalog2foldChange[i] <- log2(RTqPCR2$Vegfa_expression[i])
  RTqPCR2$Pdk1log2foldChange[i] <- log2(RTqPCR2$Pdk1_expression[i])
}

for(i in 1:nrow(RTqPCR)){
  data$Vegfa[i] <- RTqPCR2$Vegfalog2foldChange[i]
  data$Pdk1[i] <- RTqPCR2$Pdk1log2foldChange[i]
}

for(i in 1:nrow(RTqPCR2)){
  Gene_expression[i,3] <- RTqPCR2$ddCT_Vegfa[i]*(-1)
  Gene_expression[i,4] <- RTqPCR2$ddCT_Pdk1[i]*(-1)
}

names(Gene_expression) <- c("Group", "Hmox1", "Vegfa", "Pdk1")

#Histology scores
histo <- read.delim("Histology.txt")
histo$Group <- as.factor(histo$Group)
histo$Group <- factor(histo$Group, levels(histo$Group)[c(4,2,3,1)])
##The data from the histology evaluation is by definition non-parametric; hence, 
##it will not be subjected to the normality test, since datasets in which all 
##the values of a variable are equal trigger an error.

#Histology descriptive statistics (Median and Interquartile range)
descriptive_histo <- data.frame()
for(i in 1:4){
  descriptive_histo[i,1] <- levels(histo$Group)[i]
  for(j in 1:(ncol(histo)-1)){
    descriptive_histo[i,4*j-2] <- median(histo[histo$Group == levels(histo$Group)[i], j+1])
    names(descriptive_histo)[4*j-2] <- names(histo)[j+1]
    descriptive_histo[i,4*j-1] <- quantile(histo[histo$Group == levels(histo$Group)[i], j+1], 1/4)
    names(descriptive_histo)[4*j-1] <- paste(names(histo)[j+1], "1st", sep = "_")
    descriptive_histo[i,4*j] <- quantile(histo[histo$Group == levels(histo$Group)[i], j+1], 3/4)
    names(descriptive_histo)[4*j] <- paste(names(histo)[j+1], "3rd", sep = "_")
    descriptive_histo[i,4*j+1] <- IQR(histo[histo$Group == levels(histo$Group)[i], j+1])
    names(descriptive_histo)[4*j+1] <- paste(names(histo)[j+1], "IQR", sep = "_")
  }
}
names(descriptive_histo)[1] <- "Group"
descriptive_histo$Group <- as.factor(descriptive_histo$Group)
descriptive_histo$Group <- factor(descriptive_histo$Group, levels(descriptive_histo$Group)[c(4,2,3,1)])

#Normality test
shapiro <- data.frame()
for(i in 1:4){
  shapiro[i,1] <- levels(data$Group)[i]
  for(j in 2:ncol(data)){
    shapiro[i,j] <-  shapiro.test(data[data$Group == levels(data$Group)[i], j])$p
    names(shapiro)[j] <- names(data)[j]
  }
}
names(shapiro)[1] <- "Group"

#ANOVA
p_values <- data.frame()
for(i in 2:ncol(data)){
  ANOVA <- aov(data[,i] ~ Group, data)
  Tukey <- as.data.frame(TukeyHSD(ANOVA)$Group)
  p_values[1:6,1] <- rownames(Tukey)
  p_values[,i] <- Tukey$"p adj"
  names(p_values)[i] <- names(data)[i]
}
names(p_values)[1] <- "Comparisons"

#Kruskal-Wallis
## The total bilirubin dataset did not pass the normality test. Hence, a Kruskal-Wallis test with Dunn's post hoc was performed.
tb <- kwAllPairsDunnTest(TB ~ Group, data = data, p.adjust.method = "bonferroni")
tb <- tb$p.value

##Histological non-parametric analysis
Congestion <- kwAllPairsDunnTest(Congestion ~ Group, data = histo, p.adjust.method = "bonferroni")
Congestion <- Congestion$p.value

Vacuolization <- kwAllPairsDunnTest(Vacuolization ~ Group, data = histo, p.adjust.method = "bonferroni")
Vacuolization <- Vacuolization$p.value

Necrosis <- kwAllPairsDunnTest(Necrosis ~ Group, data = histo, p.adjust.method = "bonferroni")
Necrosis <- Necrosis$p.value

#Graphs
data_summary <- function(x) {
  m <- mean(x)
  ymin <- m-sd(x)
  ymax <- m+sd(x)
  return(c(y=m,ymin=ymin,ymax=ymax))
}

graph <- function(biomarker, label, binwidth){
  data %>% ggplot(aes(x = Group, y = {{biomarker}}, fill = Group)) +
    stat_summary(fun = "mean", geom = "bar", color = "black", width = 0.7, lwd = 1) +
    stat_summary(fun.data = data_summary, geom = "errorbar", width = 0.5, lwd = 1) +
    geom_dotplot(binaxis ="y", stackdir = "center", method = "dotdensity", binwidth = {{binwidth}}, dotsize = 1, stroke = 2) +
    scale_fill_manual(values = c("#FFFFFF", "#55A0FB", "#FF0000", "#106c97")) + ylab({{label}}) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
    theme(axis.title.x = element_blank(), legend.position="none", axis.line = 
            element_line(colour = "black", size = 1), panel.background = element_blank(),
          axis.text.x = element_text(size = 16, colour = "black", face = "bold"), 
          axis.text.y = element_text(size = 16, colour = "black", face = "bold"),
          axis.title.y = element_text(size = 18, colour = "black", face = "bold"),
          axis.ticks = element_line(size = 1, colour = "black"))
}

graph_mol <- function(biomarker, label, binwidth){
  Gene_expression %>% ggplot(aes(x = Group, y = {{biomarker}}, fill = Group)) + geom_hline(yintercept = 0, size = 1) +
    stat_summary(fun = "mean", geom = "bar", color = "black", width = 0.7, lwd = 1) +
    stat_summary(fun.data = data_summary, geom = "errorbar", width = 0.5, lwd = 1) +
    geom_dotplot(binaxis ="y", stackdir = "center", method = "dotdensity", binwidth = {{binwidth}}, dotsize = 1, stroke = 2) +
    scale_fill_manual(values = c("#FFFFFF", "#55A0FB", "#FF0000", "#106c97")) + ylab({{label}}) +
    theme(axis.title.x = element_blank(), legend.position="none", axis.line = 
            element_line(colour = "black", size = 1), panel.background = element_blank(), 
          axis.text.x = element_text(size = 16, colour = "black", face = "bold"), 
          axis.text.y = element_text(size = 16, colour = "black", face = "bold"),
          axis.title.y = element_text(size = 18, colour = "black", face = "bold"),
          axis.ticks = element_line(size = 1, colour = "black"))
}

#Binwidths are determined as 1/30 of the highest value for each biomarker (or rank in graphs spanning to negative values)

add_p_value <- function(graph, biomarker_as_string, position_IR, position_HGIR, position_HGTox){
  if(p_values[which(p_values$Comparisons == "IR-SH"), biomarker_as_string] < 0.05){
    graph <- graph + geom_text(aes(x = "IR", y = {{position_IR}}), label = "*", size = 8)}
  else{}
  if(p_values[which(p_values$Comparisons == "HGIR-SH"), biomarker_as_string] < 0.05 & 
     p_values[which(p_values$Comparisons == "HGIR-IR"), biomarker_as_string] < 0.05 &
     p_values[which(p_values$Comparisons == "IR-SH"), biomarker_as_string] >= 0.05){
    graph <- graph + geom_text(aes(x = "HGIR", y = {{position_HGIR}}), label = "#,‡", size = 6)}
  else{}
  if(p_values[which(p_values$Comparisons == "HGIR-IR"), biomarker_as_string] < 0.05 &
     p_values[which(p_values$Comparisons == "IR-SH"), biomarker_as_string] < 0.05){
    graph <- graph + geom_text(aes(x = "HGIR", y = {{position_HGIR}}), label = "#", size = 6)}
  else{}
  if(p_values[which(p_values$Comparisons == "HGTox-SH"), biomarker_as_string] < 0.05){
    graph <- graph + geom_text(aes(x = "HGTox", y = {{position_HGTox}}), label = "†", size = 6)}
  else{}
  if(p_values[which(p_values$Comparisons == "HGIR-SH"), biomarker_as_string] < 0.05 & 
     p_values[which(p_values$Comparisons == "HGIR-IR"), biomarker_as_string] >= 0.05 &
     p_values[which(p_values$Comparisons == "IR-SH"), biomarker_as_string] >= 0.05){
    graph <- graph + geom_text(aes(x = "HGIR", y = {{position_HGIR}}), label = "‡", size = 6)}
  else{graph}
}

#Graphs are rendered in a PDF format, with a size of 16x10 inches (panel of 6 graphs) or 16x5 inches (panel of 3 graphs)

#Animal weight
Weight <- graph(Weight, "Body weight (g)", 9.33)
Weight_p <- add_p_value(Weight, "Weight", 275, 275, 275)

Weight_liver <- graph(Weight_liver, "Liver weight (g)", 0.30)
Weight_liver_p <- add_p_value(Weight_liver, "Weight_liver", 8, 7.7, 8)

Liver_rat <- graph(Liver_rat, "Liver/Body weight", 0.0013)
Liver_rat_p <- add_p_value(Liver_rat, "Liver_rat", 0.035, 0.035, 0.040)

plot_grid(Weight_p, Weight_liver_p, Liver_rat_p, nrow = 1, ncol = 3, labels = "AUTO", label_size = 24, scale = 0.9)

#Biochemical markers
ALP <- graph(ALP, "ALP (U/L)", 12.66)
ALP_p <- add_p_value(ALP, "ALP", 350, 250, 200)

ALT <- graph(ALT, "ALT (U/L)", 75.33)
ALT_p <- add_p_value(ALT, "ALT", 2500, 1600, 250)

AST <- graph(AST, "AST (U/L)", 113.33)
AST_p <- add_p_value(AST, "AST", 3700, 3200, 500)

GLU <- graph(GLU, "Glucose (mg/dL)", 15.33)
GLU_p <- add_p_value(GLU, "GLU", 370, 300, 500)

LDH <- graph(LDH, "LDH (U/L)", 621.33)
LDH_p <- add_p_value(LDH, "LDH", 20000, 14000, 5000)

TB <- graph(TB, "Total bilirubin (mg/dL)", 0.0167)
## K-W test and post hoc test did not show significant differences. No symbols were added to the graph.

plot_grid(ALP_p, ALT_p, AST_p, GLU_p, LDH_p, TB, nrow = 2, ncol = 3, labels = "AUTO", label_size = 24, scale = 0.9)

#Proinflammatory cytokines
IL_1b <- graph(IL_1b, "IL-1\u03b2 (pg/mg of protein)", 0.6087)
IL_1b_p <- add_p_value(IL_1b, "IL_1b", 20, 19, 12)

IL_6 <- graph(IL_6, "IL-6 (pg/mg of protein)", 1.306)
IL_6_p <- add_p_value(IL_6, "IL_6", 35, 23, 32)

TNF_a <- graph(TNF_a, "TNF (pg/mg of protein)", 0.4163)
TNF_a_p <- add_p_value(TNF_a, "TNF_a", 10, 9, 9)

plot_grid(IL_1b_p, IL_6_p, TNF_a_p, nrow = 1, ncol = 3, labels = "AUTO", label_size = 24, scale = 0.9)

#Oxidative stress biomarkers
MDA <- graph(MDA, "MDA (nmol/mg of protein)", 31.19)
MDA_p <- add_p_value(MDA, "MDA", 750, 800, 1000)

SOD <- graph(SOD, "SOD (IU/mg of protein)", 1.3567)
SOD_p <- add_p_value(SOD, "SOD", 35, 40, 40)

GPx <- graph(GPx, "GPx (nmol/min/mg of protein)", 33.58)
GPx_p <- add_p_value(GPx, "GPx", 800, 900, 100)

plot_grid(MDA_p, SOD_p, GPx_p, nrow = 1, ncol = 3, labels = "AUTO", label_size = 24, scale = 0.9)

#RTqPCR
Hmox1 <- graph_mol(Hmox1, expression(italic("Hmox1")~"(-"*Delta*Delta*"C"[T]*" from SH)"), 0.06887)
Hmox1_p <- add_p_value(Hmox1, "Hmox1", 0.2, 0.4, 0.6)

Vegfa <- graph_mol(Vegfa, expression(italic("Vegfa")~"(-"*Delta*Delta*"C"[T]*" from SH)"), 0.1136)
Vegfa_p <- add_p_value(Vegfa, "Vegfa", 0.5, 0.5, 0.5)

Pdk1 <- graph_mol(Pdk1, expression(italic("Pdk1")~"(-"*Delta*Delta*"C"[T]*" from SH)"), 0.1423)
Pdk1_p <- add_p_value(Pdk1, "Pdk1", 0.5, 0.5, 0.5)

plot_grid(Hmox1_p, Vegfa_p, Pdk1_p, nrow = 1, ncol = 3, labels = "AUTO", label_size = 24, scale = 0.9)