#!/usr/bin/env Rscript  
#########################################################
# TASK: producing boxplots for dr and other muts

#########################################################
#=======================================================================
# working dir and loading libraries
getwd()
setwd("~/git/LSHTM_analysis/scripts/plotting")
getwd()

source("~/git/LSHTM_analysis/scripts/Header_TT.R")
#library(ggplot2)
#library(data.table)
#library(dplyr)
source("other_plots_data.R")

#=======
# output
#=======
dr_other_combined = "dr_other_muts.svg"
plot_dr_other_combined  =  paste0(plotdir,"/", dr_other_combined)

dr_other_combined_labelled = "dr_other_muts_labelled.svg"
plot_dr_other_combined_labelled  =  paste0(plotdir,"/", dr_other_combined_labelled)

#dr_other_foldx = "dr_other_muts_foldx.svg"
#plot_dr_other_foldx  =  paste0(plotdir,"/", dr_other_foldx)

#point_legend = "point_legend.svg"
#plot_point_legend  =  paste0(plotdir,"/", point_legend)
########################################################################
#               end of data extraction and cleaning for plots          #
########################################################################
#my_comparisons <- list( c(dr_muts_col, other_muts_col) )
my_comparisons <- list( c("DM", "OM") )

my_ats = 22 # axis text size
my_als = 20 # axis label size
my_fls = 20 # facet label size
my_pts = 22 # plot title size

##########################################################################

#===========
# Plot1: PS
#===========
# my_stat_ps = compare_means(param_value~mutation_info
#                            , group.by = "param_type"
#                            , data = df_lf_ps
#                            , paired = FALSE
#                            , p.adjust.method = "BH")

y_value = "param_value"

p1 = ggplot(lf_duet, aes(x = mutation_info
                    , y = eval(parse(text=y_value)) ))  + 
  facet_wrap(~ param_type
             , nrow = 1
             , scales = "free_y") + 
  geom_boxplot(fill = "white", outlier.colour = NA
               #, position = position_dodge(width = 0.9)
               , width = 0.2) +
  geom_point(position = position_jitterdodge(dodge.width=0.01)
             , alpha = 0.5
             , show.legend = FALSE
             , aes(colour = duet_outcome)) +
  theme(axis.text.x = element_text(size = my_ats)
        , axis.text.y = element_text(size = my_ats
                                     , angle = 0
                                     , hjust = 1
                                     , vjust = 0)
        , axis.title.x = element_text(size = my_ats)
        , axis.title.y = element_text(size = my_ats)
        , plot.title = element_text(size = my_pts , hjust = 0.5, colour = "black", face = "bold")
        , strip.background = element_rect(fill = "khaki2")
        , strip.text.x = element_text(size = my_fls, colour = "black")
        , legend.title = element_text(color = "black", size = my_als)
        , legend.text = element_text(size = my_ats)
        , legend.direction = "vertical") +
  labs(title = "DUET"
       , x = ""
       , y = "")+ 
  stat_compare_means(comparisons = my_comparisons
                     , method = "wilcox.test"
                     , paired = FALSE
                     #, label = "p.format")
                     , label = "p.signif")
p1

##########################################################################
#=================
# Plot 2: Foldx
#=================
my_stat_foldx = compare_means(param_value~mutation_info, group.by = "param_type"
                              , data = df_lf_foldx,  paired = FALSE, p.adjust.method = "BH")

y_value = "param_value"

p2 = ggplot(df_lf_foldx, aes(x = mutation_info
                             , y = eval(parse(text=y_value)) ))  + 
  facet_wrap(~ param_type
             , nrow = 1
             , scales = "free_y") + 
  geom_boxplot(fill = "white", outlier.colour = NA
               #, position = position_dodge(width = 0.2)
               , width = 0.5) +
  geom_point(position = position_jitterdodge(dodge.width=0.01)
             , alpha = 0.5
             , show.legend = FALSE
             , aes(colour = factor(foldx_outcome))) +
  theme(axis.text.x = element_text(size = my_ats)
        , axis.text.y = element_text(size = my_ats
                                     , angle = 0
                                     , hjust = 1
                                     , vjust = 0)
        , axis.title.x = element_text(size = my_ats)
        , axis.title.y = element_text(size = my_ats)
        , plot.title = element_text(size = my_pts , hjust = 0.5, colour = "black", face = "bold")
        , strip.background = element_rect(fill = "palegreen1")
        , strip.text.x = element_text(size = my_fls, colour = "black")
        , legend.title = element_text(color = "black", size = my_als)
        , legend.text = element_text(size = my_ats)
        #, legend.title = element_blank()
        , legend.direction = "vertical") +
  labs(title = "Foldx"
       , x = ""
       , y = "") + 
  stat_compare_means(comparisons = my_comparisons
                     , method = "wilcox.test"
                     , paired = FALSE
                     #, label = "p.format")
                     , label = "p.signif")

p2 

##########################################################################
#==============
# Plot 3: LIG
#==============

my_stat_lig = compare_means(param_value~mutation_info, group.by = "param_type"
                            , data = df_lf_lig,  paired = FALSE, p.adjust.method = "BH")

y_value = "param_value"

p3 = ggplot(df_lf_lig, aes(x = mutation_info
                           , y = eval(parse(text=y_value)) ))  + 
  facet_wrap(~ param_type
             , nrow = 1
             , scales = "free_y") + 
  geom_boxplot(fill = "white", outlier.colour = NA
               #, position = position_dodge(width = 0.9)
               , width = 0.2) +
  geom_point(position = position_jitterdodge(dodge.width=0.01)
             , alpha = 0.5
             , show.legend = FALSE
             , aes(colour = factor(ligand_outcome))) +
  theme(axis.text.x = element_text(size = my_ats)
        , axis.text.y = element_text(size = my_ats
                                     , angle = 0
                                     , hjust = 1
                                     , vjust = 0)
        , axis.title.x = element_text(size = my_ats)
        , axis.title.y = element_text(size = my_ats)
        , plot.title = element_text(size = my_pts , hjust = 0.5, colour = "black", face = "bold")
        , strip.background = element_rect(fill = "thistle3")
        , strip.text.x = element_text(size = my_fls, colour = "black")
        , legend.title = element_text(color = "black", size = my_als)
        , legend.text = element_text(size = my_ats)
        , legend.direction = "vertical") +
  labs(title = "Ligand Affinity"
       , x = ""
       , y = "")+ 
  stat_compare_means(comparisons = my_comparisons
                     , method = "wilcox.test"
                     , paired = FALSE
                     #, label = "p.format")
                     , label = "p.signif")

p3

##########################################################################
#===========
# legend
#===========
legend = ggplot(df_lf_foldx, aes(x = mutation_info
                             , y = eval(parse(text=y_value)) ))  + 
  facet_wrap(~ param_type
             , nrow = 1
             , scales = "free_y") + 
  geom_boxplot(fill = "white", outlier.colour = NA
               #, position = position_dodge(width = 0.2)
               , width = 0.5) +
  geom_point(position = position_jitterdodge(dodge.width=0.01)
             , aes(colour = factor(foldx_outcome))) +
  theme(legend.text = element_text(size = 25)
        , legend.direction = "vertical"
        , legend.title = element_blank()
        #, axis.text.y = element_text(size = 30, colour = "black", angle = 90)
        #, axis.text.x = element_text(size = 30, colour = "black", angle = 0)
        #, strip.background = element_rect(fill="orange")
        ) +
  guides(color = guide_legend(override.aes = list(size = 6)))

legend
##########################################################################
#===========================
# combine
#===========================
#---------
# plot 1
#---------
cat("Output plot:",plot_dr_other_combined)
svg(plot_dr_other_combined, width = 24, height = 12)

#theme_set(theme_gray()) # to preserve default theme


# first combine fold and lig plots (2 and 3)
c1 = cowplot::plot_grid(p2, p3
                        #, labels = c("(b)", "(c)")
                        , labels = c("B", "C")
                        , hjust = -2
                        , label_size = 25
                        , nrow = 1
                        , rel_widths = c(1/6, 5/6))
c1

OutPlot = cowplot::plot_grid(p1, c1
                        , nrow = 2)
print(OutPlot)
dev.off()

#==========================================
# labelled
cat("Output plot:", plot_dr_other_combined_labelled)
svg(plot_dr_other_combined_labelled, width = 24, height = 12)
plot_dr_other_combined_labelled

Outplot_labelled = cowplot::plot_grid(p1, c1
                        , nrow = 2
                        #, labels = c("(a)", "", "")
                        , labels = c("A", "", "")
                        , hjust = -2
                        #, vjust = 0
                        , label_size = 25)

print(Outplot_labelled)

dev.off()

#---------
# plot 3: legend
#---------
#svg(plot_point_legend, width = 6, height = 7)
#OutPlot3 =  legend
#print(OutPlot3)
#dev.off()