From 17081949123a6a2d0759af6a29cb811c81e63a71 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Thu, 10 Sep 2020 14:09:12 +0100 Subject: [PATCH] added boxplots and stats for other numerical params --- scripts/plotting/other_plots.R | 173 +++++++++++++++++++ scripts/plotting/other_plots_data.R | 252 ++++++++++++++++++++++++++++ 2 files changed, 425 insertions(+) create mode 100644 scripts/plotting/other_plots.R create mode 100644 scripts/plotting/other_plots_data.R diff --git a/scripts/plotting/other_plots.R b/scripts/plotting/other_plots.R new file mode 100644 index 0000000..9fccb46 --- /dev/null +++ b/scripts/plotting/other_plots.R @@ -0,0 +1,173 @@ +#!/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("Header_TT.R") +#library(ggplot2) +#library(data.table) +#library(dplyr) +source("other_plots_data.R") + +#======= +# output +#======= +#dr_other_plots_combined = "dr_other_combined.svg" +#plot_dr_other_plots_combined = paste0(plotdir,"/", dr_other_plots_combined) + + + +######################################################################## +# end of data extraction and cleaning for plots # +######################################################################## +#my_comparisons <- list( c(dr_muts_col, other_muts_col) ) +my_comparisons <- list( c("DM", "OM") ) + +#=========== +# 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(df_lf_ps, 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.5) + + geom_point(position = position_jitterdodge(dodge.width=0.01) + , show.legend = FALSE + , aes(colour = factor(duet_outcome))) + + theme(axis.text.x = element_text(size = 15) + , axis.text.y = element_text(size = 15 + , angle = 0 + , hjust = 1 + , vjust = 0) + , axis.title.x = element_text(size = 15) + , axis.title.y = element_text(size = 15) + , plot.title = element_text(size = 20, hjust = 0.5) + , strip.text.x = element_text(size = 15, colour = "black") + , legend.title = element_text(color = "black", size = 20) + , legend.text = element_text(size = 15) + , 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.9) + , width = 0.5) + + geom_point(position = position_jitterdodge(dodge.width=0.01) + , show.legend = FALSE + , aes(colour = factor(foldx_outcome))) + + theme(axis.text.x = element_text(size = 15) + , axis.text.y = element_text(size = 15 + , angle = 0 + , hjust = 1 + , vjust = 0) + , axis.title.x = element_text(size = 15) + , axis.title.y = element_text(size = 15) + , plot.title = element_text(size = 20, hjust = 0.5) + , strip.text.x = element_text(size = 15, colour = "black") + , legend.title = element_text(color = "black", size = 20) + , legend.text = element_text(size = 15) + , 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.5) + + geom_point(position = position_jitterdodge(dodge.width=0.01) + , show.legend = FALSE + , aes(colour = factor(ligand_outcome))) + + theme(axis.text.x = element_text(size = 15) + , axis.text.y = element_text(size = 15 + , angle = 0 + , hjust = 1 + , vjust = 0) + , axis.title.x = element_text(size = 15) + , axis.title.y = element_text(size = 15) + , plot.title = element_text(size = 20, hjust = 0.5) + , strip.text.x = element_text(size = 15, colour = "black") + , legend.title = element_text(color = "black", size = 20) + , legend.text = element_text(size = 15) + , 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 + +#=========================== +# combine +#=========================== +#svg(plot_or_combined, width = 32, height = 12) + +theme_set(theme_gray()) # to preserve default theme + +printFile = cowplot::plot_grid(plot_grid(p1, p2, p3 + , nrow = 3 + , align = 'h' + , labels = c("", "", "") + , label_size = 20)) +print(printFile) +dev.off() + diff --git a/scripts/plotting/other_plots_data.R b/scripts/plotting/other_plots_data.R new file mode 100644 index 0000000..8f6836a --- /dev/null +++ b/scripts/plotting/other_plots_data.R @@ -0,0 +1,252 @@ +#!/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("Header_TT.R") +library(ggplot2) +library(data.table) +library(dplyr) +source("combining_dfs_plotting.R") + +#======= +# output +#======= +#lineage_dist_combined = "lineage_dist_combined_PS.svg" +#plot_lineage_dist_combined = paste0(plotdir,"/", lineage_dist_combined) + + +rm(merged_df2, merged_df2_comp, merged_df2_lig, merged_df2_comp_lig + , merged_df3_comp, merged_df3_comp_lig + , my_df_u, my_df_u_lig) +#======================================================================== +#%%%%%%%%%%%%%%%%%%% +# REASSIGNMENT: PS +#%%%%%%%%%%%%%%%%%%%% +df_ps = merged_df3 + +# name tidying +df_ps$mutation_info = as.factor(df_ps$mutation_info) +df_ps$duet_outcome = as.factor(df_ps$duet_outcome) +df_ps$foldx_outcome = as.factor(df_ps$foldx_outcome) +df_ps$ligand_outcome = as.factor(df_ps$ligand_outcome) + +# check +table(df_ps$mutation_info) + +#%%%%%%%%%%%%%%%%%%% +# REASSIGNMENT: LIG +#%%%%%%%%%%%%%%%%%%%% + +df_lig = merged_df3_lig + +# name tidying +df_lig$mutation_info = as.factor(df_lig$mutation_info) +df_lig$duet_outcome = as.factor(df_lig$duet_outcome) +df_lig$ligand_outcome = as.factor(df_lig$ligand_outcome) + +# check +table(df_lig$mutation_info) + +#======================================================================== +n = which(colnames(df_ps) == "ddg"); n # 10 + +my_min = min(df_ps[,n]); my_min +my_max = max(df_ps[,n]); my_max + +df_ps$foldx_scaled = ifelse(df_ps[,n] < 0 + , df_ps[,n]/abs(my_min) + , df_ps[,n]/my_max) #335, 15 +#sanity check +my_min = min(df_ps$foldx_scaled); my_min +my_max = max(df_ps$foldx_scaled); my_max + +# assign foldx +#ddg<0 = "Stabilising" (-ve) +c1 = table(df_ps$ddg < 0) +df_ps$foldx_outcome = ifelse(df_ps$ddg < 0, "Stabilising", "Destabilising") +c2 = table(df_ps$ddg < 0) + +if ( all(c1 == c2) ){ + cat("PASS: foldx outcome successfully created") +}else{ + cat("FAIL: foldx outcome could not be created. Aborting!") + exit() +} + +#=========== +# Data: ps +#=========== +# keep similar dtypes cols together +cols_to_select_ps = c("mutationinformation", "mutation", "position", "mutation_info" + , "duet_outcome" + + , "duet_scaled" + , "ligand_distance" + , "asa" + , "rsa" + , "rd_values" + , "kd_values") + +df_wf_ps = df_ps[, cols_to_select_ps] + +pivot_cols_ps = cols_to_select_ps[1:5]; pivot_cols_ps + +expected_rows_lf_ps = nrow(df_wf_ps) * (length(df_wf_ps) - length(pivot_cols_ps)) +expected_rows_lf_ps + +# LF data: duet +df_lf_ps = gather(df_wf_ps, param_type, param_value, duet_scaled:kd_values, factor_key=TRUE) + +if (nrow(df_lf_ps) == expected_rows_lf_ps){ + cat("PASS: long format data created for duet") +}else{ + cat("FAIL: long format data could not be created for duet") + exit() +} + +str(df_wf_ps) +str(df_lf_ps) + +# assign pretty labels: param_type +levels(df_lf_ps$param_type); table(df_lf_ps$param_type) + +levels(df_lf_ps$param_type)[levels(df_lf_ps$param_type)=="duet_scaled"] <- "DUET" +levels(df_lf_ps$param_type)[levels(df_lf_ps$param_type)=="ligand_distance"] <- "Ligand Distance" +levels(df_lf_ps$param_type)[levels(df_lf_ps$param_type)=="asa"] <- "ASA" +levels(df_lf_ps$param_type)[levels(df_lf_ps$param_type)=="rsa"] <- "RSA" +levels(df_lf_ps$param_type)[levels(df_lf_ps$param_type)=="rd_values"] <- "RD" +levels(df_lf_ps$param_type)[levels(df_lf_ps$param_type)=="kd_values"] <- "KD" +# check +levels(df_lf_ps$param_type); table(df_lf_ps$param_type) + +# assign pretty labels: mutation_info +levels(df_lf_ps$mutation_info); table(df_lf_ps$mutation_info) +sum(table(df_lf_ps$mutation_info)) == nrow(df_lf_ps) + +levels(df_lf_ps$mutation_info)[levels(df_lf_ps$mutation_info)==dr_muts_col] <- "DM" +levels(df_lf_ps$mutation_info)[levels(df_lf_ps$mutation_info)==other_muts_col] <- "OM" +# check +levels(df_lf_ps$mutation_info); table(df_lf_ps$mutation_info) + +############################################################################ +#=========== +# Data: foldx +#=========== +# keep similar dtypes cols together +cols_to_select_foldx = c("mutationinformation", "mutation", "position", "mutation_info" + , "foldx_outcome" + + , "foldx_scaled" + , "ligand_distance" + , "asa" + , "rsa" + , "rd_values" + , "kd_values") + + +df_wf_foldx = df_ps[, cols_to_select_foldx] + +pivot_cols_foldx = cols_to_select_foldx[1:5]; pivot_cols_foldx + +expected_rows_lf_foldx = nrow(df_wf_foldx) * (length(df_wf_foldx) - length(pivot_cols_foldx)) +expected_rows_lf_foldx + +# LF data: foldx +df_lf_foldx = gather(df_wf_foldx, param_type, param_value, foldx_scaled:kd_values, factor_key=TRUE) + +if (nrow(df_lf_foldx) == expected_rows_lf_foldx){ + cat("PASS: long format data created for foldx") +}else{ + cat("FAIL: long format data could not be created for foldx") + exit() +} + +# assign pretty labels: param type +levels(df_lf_foldx$param_type); table(df_lf_foldx$param_type) + +levels(df_lf_foldx$param_type)[levels(df_lf_foldx$param_type)=="foldx_scaled"] <- "Foldx" +levels(df_lf_foldx$param_type)[levels(df_lf_foldx$param_type)=="ligand_distance"] <- "Ligand Distance" +levels(df_lf_foldx$param_type)[levels(df_lf_foldx$param_type)=="asa"] <- "ASA" +levels(df_lf_foldx$param_type)[levels(df_lf_foldx$param_type)=="rsa"] <- "RSA" +levels(df_lf_foldx$param_type)[levels(df_lf_foldx$param_type)=="rd_values"] <- "RD" +levels(df_lf_foldx$param_type)[levels(df_lf_foldx$param_type)=="kd_values"] <- "KD" +# check +levels(df_lf_foldx$param_type); table(df_lf_foldx$param_type) + +# assign pretty labels: mutation_info +levels(df_lf_foldx$mutation_info); table(df_lf_foldx$mutation_info) +sum(table(df_lf_foldx$mutation_info)) == nrow(df_lf_foldx) + +levels(df_lf_foldx$mutation_info)[levels(df_lf_foldx$mutation_info)==dr_muts_col] <- "DM" +levels(df_lf_foldx$mutation_info)[levels(df_lf_foldx$mutation_info)==other_muts_col] <- "OM" +# check +levels(df_lf_foldx$mutation_info); table(df_lf_foldx$mutation_info) + +############################################################################ +#=========== +# LF data: LIG +#=========== +# keep similar dtypes cols together +cols_to_select_lig = c("mutationinformation", "mutation", "position", "mutation_info" + , "ligand_outcome" + + , "affinity_scaled" + , "ligand_distance" + , "asa" + , "rsa" + , "rd_values" + , "kd_values") + +df_wf_lig = df_lig[, cols_to_select_lig] + +pivot_cols_lig = cols_to_select_lig[1:5]; pivot_cols_lig + +expected_rows_lf_lig = nrow(df_wf_lig) * (length(df_wf_lig) - length(pivot_cols_lig)) +expected_rows_lf_lig + +# LF data: foldx +df_lf_lig = gather(df_wf_lig, param_type, param_value, affinity_scaled:kd_values, factor_key=TRUE) + +if (nrow(df_lf_lig) == expected_rows_lf_lig){ + cat("PASS: long format data created for foldx") +}else{ + cat("FAIL: long format data could not be created for foldx") + exit() +} + +# assign pretty labels: param_type +levels(df_lf_lig$param_type); table(df_lf_lig$param_type) + +levels(df_lf_lig$param_type)[levels(df_lf_lig$param_type)=="affinity_scaled"] <- "Ligand Affinity" +levels(df_lf_lig$param_type)[levels(df_lf_lig$param_type)=="ligand_distance"] <- "Ligand Distance" +levels(df_lf_lig$param_type)[levels(df_lf_lig$param_type)=="asa"] <- "ASA" +levels(df_lf_lig$param_type)[levels(df_lf_lig$param_type)=="rsa"] <- "RSA" +levels(df_lf_lig$param_type)[levels(df_lf_lig$param_type)=="rd_values"] <- "RD" +levels(df_lf_lig$param_type)[levels(df_lf_lig$param_type)=="kd_values"] <- "KD" +#check +levels(df_lf_lig$param_type); table(df_lf_lig$param_type) + +# assign pretty labels: mutation_info +levels(df_lf_lig$mutation_info); table(df_lf_lig$mutation_info) +sum(table(df_lf_lig$mutation_info)) == nrow(df_lf_lig) + +levels(df_lf_lig$mutation_info)[levels(df_lf_lig$mutation_info)==dr_muts_col] <- "DM" +levels(df_lf_lig$mutation_info)[levels(df_lf_lig$mutation_info)==other_muts_col] <- "OM" +# check +levels(df_lf_lig$mutation_info); table(df_lf_lig$mutation_info) + +############################################################################### + +# clear excess variables +rm(cols_to_select_ps, cols_to_select_foldx, cols_to_select_lig + , pivot_cols_ps, pivot_cols_foldx, pivot_cols_lig + , expected_rows_lf_ps, expected_rows_lf_foldx, expected_rows_lf_lig + , my_max, my_min, na_count, na_count_df2, na_count_df3, dup_muts_nu + , c1, c2, n)