diff --git a/scripts/plotting/corr_plots.R b/scripts/plotting/corr_plots.R new file mode 100755 index 0000000..8ab78ef --- /dev/null +++ b/scripts/plotting/corr_plots.R @@ -0,0 +1,242 @@ +#!/usr/bin/env Rscript +getwd() +setwd("~/git/LSHTM_analysis/scripts/plotting") +getwd() +source("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 +######################################################################= \ No newline at end of file