#****************************************************************************************
## Comparison of IRS distributions for PGM-Set and 4-way set
#****************************************************************************************
library("ggpubr")
library("ggthemes")
# Load IES table
# Load Table
IES_tab <- read.table(' path to IES table "Supplemental Table S2" ', h = T, stringsAsFactors = F)
IES_tab <- IES_tab[(IES_tab$IES_PLUS_25F0 + IES_tab$IES_MINUS_25F0) > 20 & 
                     (IES_tab$IES_PLUS_25F1 + IES_tab$IES_MINUS_25F1) > 20 &
                     (IES_tab$IES_PLUS_18F1 + IES_tab$IES_MINUS_18F1) > 20 &
                     (IES_tab$IES_PLUS_32F1 + IES_tab$IES_MINUS_32F1) > 20, ]

# Create combined dataframe
F025_all <- as.data.frame(IES_tab$IRS_25F0[IES_tab$IRS_25F0 > 0 &
                                             !is.na(IES_tab$IES_ID)])
colnames(F025_all) <- "IRS"
F025_all$Set <- "All_25F0"

F125_all <- as.data.frame(IES_tab$IRS_25F1[IES_tab$IRS_25F1 > 0 &
                                             !is.na(IES_tab$IES_ID)])
colnames(F125_all) <- "IRS"
F125_all$Set <- "All_25F1"

F118_all <- as.data.frame(IES_tab$IRS_18F1[IES_tab$IRS_18F1 > 0 & 
                                             !is.na(IES_tab$IES_ID)])
colnames(F118_all) <- "IRS"
F118_all$Set <- "All_18F1"

F132_all <- as.data.frame(IES_tab$IRS_32F1[IES_tab$IRS_32F1 > 0 &
                                             !is.na(IES_tab$IES_ID)])
colnames(F132_all) <- "IRS"
F132_all$Set <- "All_32F1"

F025_4way <- as.data.frame(na.omit(IES_tab$IRS_25F0[IES_tab$IRS_25F0 > 0 &
                                                      IES_tab$IRS_25F1 > 0 &
                                                      IES_tab$IRS_18F1 > 0 &
                                                      IES_tab$IRS_32F1 > 0]))
colnames(F025_4way) <- "IRS"
F025_4way$Set <- "4way_25F0"

F125_4way <- as.data.frame(na.omit(IES_tab$IRS_25F1[IES_tab$IRS_25F0 > 0 &
                                                      IES_tab$IRS_25F1 > 0 &
                                                      IES_tab$IRS_18F1 > 0 &
                                                      IES_tab$IRS_32F1 > 0]))
colnames(F125_4way) <- "IRS"
F125_4way$Set <- "4way_25F1"

F118_4way <- as.data.frame(na.omit(IES_tab$IRS_18F1[IES_tab$IRS_25F0 > 0 &
                                                      IES_tab$IRS_25F1 > 0 &
                                                      IES_tab$IRS_18F1 > 0 &
                                                      IES_tab$IRS_32F1 > 0]))
colnames(F118_4way) <- "IRS"
F118_4way$Set <- "4way_18F1"

F132_4way <- as.data.frame(na.omit(IES_tab$IRS_32F1[IES_tab$IRS_25F0 > 0 &
                                                      IES_tab$IRS_25F1 > 0 &
                                                      IES_tab$IRS_18F1 > 0 &
                                                      IES_tab$IRS_32F1 > 0]))
colnames(F132_4way) <- "IRS"
F132_4way$Set <- "4way_32F1"

df <- rbind(F025_all, F125_all, F118_all, F132_all, F025_4way, F125_4way, F118_4way, F132_4way)
df$IRS <- as.numeric(df$IRS)
#*******************************************************************************************************
# Box plot with ggboxplot and multiple testing
p <- ggboxplot(df, x = "Set", y = "IRS",
               fill = "Set",
               outline = F, outlier.shape=NA, notch = T, varwidth = T, 
               xlab = "", ylab = "IRS")

my_comparisons <- list( c("All_25F0", "4way_25F0"), c("All_25F1", "4way_25F1"), c("All_18F1", "4way_18F1"), c("All_32F1", "4way_32F1"), 
                        c("4way_25F0", "4way_25F1"), c("4way_25F0", "4way_18F1"), c("4way_25F0", "4way_32F1"))

p + theme_tufte()  +
  theme(axis.text.y=element_text(size=20), axis.title.y=element_text(size=20, hjust=0.3, vjust = +100), 
        axis.text.x=element_blank(), 
        legend.background = element_rect(fill= NULL, size=.5, linetype="dotted"),
        legend.text=element_text(size=10)) + 
  theme(plot.margin=unit(c(1,1,1.5,1.2),"cm")) +
  guides(fill=guide_legend(title=NULL)) +
  scale_y_continuous(breaks=seq(0,1.0,0.1)) +
  scale_fill_manual(values = c(rep("white", 4), c("medium blue", "sky blue","medium sea green", "fire brick"))) +
  scale_color_manual(values = c(c("medium blue", "sky blue","medium sea green", "fire brick"), rep("black", 4))) +
  stat_compare_means(label.y = 1.6, label.x = 1.2, size = 5) +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif", p.adjust.method = "BH", size = 0.6)

