#!/usr/bin/env Rscript
getwd()
setwd("~/git/LSHTM_analysis/scripts/plotting")
getwd()
source("~/git/LSHTM_analysis/scripts/Header_TT.R")

spec = matrix(c(
  "drug"       , "d",  1, "character",
  "gene"       , "g",  1, "character",
  "data_file1" , "fa", 2, "character",
  "data_file2" , "fb", 2, "character" 
), byrow = TRUE, ncol = 4)

opt = getopt(spec)

drug            = opt$drug
gene            = opt$gene
infile_params   = opt$data_file1
infile_metadata = opt$data_file2

if(is.null(drug)|is.null(gene)) {
  stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)")
}

#===========
# Input
#===========

source("get_plotting_dfs.R")

#===========
# output
#===========
# PS
corr_ps = "corr_PS.svg"
plot_corr_ps =  paste0(plotdir,"/", corr_ps)

corr_ps_all = "corr_PS_all.svg"
plot_corr_ps_all =  paste0(plotdir,"/", corr_ps_all)


# LIG
corr_lig = "corr_LIG.svg"
plot_corr_lig =  paste0(plotdir,"/", corr_lig)

corr_lig_all = "corr_LIG_all.svg"
plot_corr_lig_all =  paste0(plotdir,"/", corr_lig_all)

##############################################################################
foo = corr_ps_df3
#foo2 = corr_ps_df2

bar = corr_lig_df3
#bar2 = corr_lig_df2

#================================
# Data for Correlation plots: PS
#================================
# subset data to generate pairwise correlations
cols_to_select = c("DUET"
                     , "Foldx"
                     , "Log (OR)"
                     , "-Log (P)"
                     , "MAF"
                     , "duet_outcome"
                     , drug)
corr_data_ps = foo[names(foo)%in%cols_to_select]
length(cols_to_select)

colnames(corr_data_ps)

start = 1
end = which(colnames(corr_data_ps) == drug); end # should be the last column
offset = 1

my_corr_ps = corr_data_ps[start:(end - offset)]
head(my_corr_ps)

#---------------------
# Corr plot PS: short
# data: corr_ps_df3
# cols: 7
#---------------------
cat("Corr plot PS DUET with coloured dots:", plot_corr_ps)
svg(plot_corr_ps, width = 15, height = 15)

pairs.panels(my_corr_ps[1:(length(my_corr_ps)-1)]
             , method = "spearman" # correlation method
             , hist.col = "grey" ##00AFBB
             , density = TRUE  # show density plots
             , ellipses = F # show correlation ellipses
             , stars = T
             , rug = F
             , breaks = "Sturges"
             , show.points = T
             , bg = c("#f8766d", "#00bfc4")[unclass(factor(my_corr_ps$duet_outcome))] # foldx colours are reveresed
             , pch = 21 # for bg
             , jitter = T
             , alpha = 1
             , cex = 1.8
             , cex.axis = 2
             , cex.labels = 4
             , cex.cor = 1
             , smooth = F
)
dev.off()

corr_ps_rho = corr.test(my_corr_ps[1:5], method = "spearman")$r
corr_ps_p = corr.test(my_corr_ps[1:5], method = "spearman")$p

#---------------------
# Corr plot PS: ALL
# data: corr_ps_df3
# cols: 10
#---------------------
end_ps_all = which(colnames(foo) == drug); end_ps_all # should be the last column

my_corr_ps_all = foo[start:(end_ps_all - offset)]
cols_to_drop = "Mutation"
my_corr_ps_all = my_corr_ps_all[, !(names(my_corr_ps_all)%in%cols_to_drop)]
head(my_corr_ps_all)
length(colnames(my_corr_ps_all))

cat("Corr plot PS DUET with coloured dots:", plot_corr_ps_all)
svg(plot_corr_ps_all, width = 15, height = 15)

pairs.panels(my_corr_ps_all[1:(length(my_corr_ps_all)-1)]
             , method = "spearman" # correlation method
             , hist.col = "grey" ##00AFBB
             , density = TRUE  # show density plots
             , ellipses = F # show correlation ellipses
             , stars = T
             , rug = F
             , breaks = "Sturges"
             , show.points = T
             , bg = c("#f8766d", "#00bfc4")[unclass(factor(my_corr_ps_all$duet_outcome))] # foldx colours are reveresed
             , pch = 21 # for bg
             , jitter = T
             , alpha = 1
             , cex = 1.5
             , cex.axis = 2
             , cex.labels = 2.5
             , cex.cor = 1
             , smooth = F
)
dev.off()

#==================================
# Data for Correlation plots: LIG
#==================================
cols_to_select_lig = c("Ligand Affinity"
                     , "Log (OR)"
                     , "-Log (P)"
                     , "MAF"
                     , "ligand_outcome"
                     , drug)

corr_data_lig = bar[names(bar)%in%cols_to_select_lig]
length(cols_to_select_lig)

colnames(corr_data_lig)

start_lig = 1
end_lig = which(colnames(corr_data_lig) == drug); end_lig # should be the last column
offset_lig = 1

my_corr_lig = corr_data_lig[start_lig:(end_lig-offset_lig)]
head(my_corr_lig)

#---------------------
# Corr plot LIG: short
# data: corr_lig_df3
# cols: 7
#---------------------
cat("Corr LIG plot with coloured dots:", plot_corr_lig)
svg(plot_corr_lig, width = 15, height = 15)

pairs.panels(my_corr_lig[1:(length(my_corr_lig)-1)]
             , method = "spearman" # correlation method
             , hist.col = "grey" ##00AFBB
             , density = TRUE  # show density plots
             , ellipses = F # show correlation ellipses
             , stars = T
             , rug = F
             , breaks = "Sturges"
             , show.points = T
             , bg = c("#f8766d", "#00bfc4")[unclass(factor(my_corr_lig$ligand_outcome))] 
             , pch = 21 # for bg
             , jitter = T
             , cex = 2
             , cex.axis = 2
             , cex.labels = 4
             , cex.cor = 1
             , smooth = F
)

dev.off()

corr_lig_rho = corr.test(my_corr_lig[1:4], method = "spearman")$r
corr_lig_p = corr.test(my_corr_lig[1:4], method = "spearman")$p

#---------------------
# Corr plot LIG: ALL
# data: corr_lig_df3
# cols: 9
#---------------------
end_lig_all = which(colnames(bar) == drug); end_lig_all # should be the last column

my_corr_lig_all = bar[start_lig:(end_lig_all - offset_lig)]
cols_to_drop = "Mutation"
my_corr_lig_all = my_corr_lig_all[, !(names(my_corr_lig_all)%in%cols_to_drop)]
head(my_corr_lig_all)
length(colnames(my_corr_lig_all))

cat("Corr plot LIG with coloured dots:", plot_corr_lig_all)
svg(plot_corr_lig_all, width = 15, height = 15)

pairs.panels(my_corr_lig_all[1:(length(my_corr_lig_all)-1)]
             , method = "spearman" # correlation method
             , hist.col = "grey" ##00AFBB
             , density = TRUE  # show density plots
             , ellipses = F # show correlation ellipses
             , stars = T
             , rug = F
             , breaks = "Sturges"
             , show.points = T
             , bg = c("#f8766d", "#00bfc4")[unclass(factor(my_corr_lig_all$ligand_outcome))] # foldx colours are reveresed
             , pch = 21 # for bg
             , jitter = T
             , alpha = 1
             , cex = 1.5
             , cex.axis = 2
             , cex.labels = 2.2
             , cex.cor = 1
             , smooth = F
)
dev.off()


######################################################################=
#                             End of script
######################################################################=