From ebca0f42210ec43b74ed6f76e5f9483cbaf55538 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Tue, 7 Sep 2021 10:52:26 +0100 Subject: [PATCH] moved old lineage_basic_barplot.R to redundant --- .../redundant/lineage_basic_barplot.R | 214 +++++++++++++ .../redundant/other_plots_data_SAFEGUARD.R | 301 ++++++++++++++++++ 2 files changed, 515 insertions(+) create mode 100644 scripts/plotting/redundant/lineage_basic_barplot.R create mode 100644 scripts/plotting/redundant/other_plots_data_SAFEGUARD.R diff --git a/scripts/plotting/redundant/lineage_basic_barplot.R b/scripts/plotting/redundant/lineage_basic_barplot.R new file mode 100644 index 0000000..e4503d1 --- /dev/null +++ b/scripts/plotting/redundant/lineage_basic_barplot.R @@ -0,0 +1,214 @@ +#!/usr/bin/env Rscript +getwd() +setwd("~/git/LSHTM_analysis/scripts/plotting/") +getwd() +######################################################### +# TASK: Basic lineage barplot showing numbers + +# Output: Basic barplot with lineage samples and mut count + +########################################################## +# Installing and loading required packages +########################################################## +source("Header_TT.R") +require(data.table) +source("combining_dfs_plotting.R") +# should return the following dfs, directories and variables + +# PS combined: +# 1) merged_df2 +# 2) merged_df2_comp +# 3) merged_df3 +# 4) merged_df3_comp + +# LIG combined: +# 5) merged_df2_lig +# 6) merged_df2_comp_lig +# 7) merged_df3_lig +# 8) merged_df3_comp_lig + +# 9) my_df_u +# 10) my_df_u_lig + +cat("Directories imported:" + , "\n====================" + , "\ndatadir:", datadir + , "\nindir:", indir + , "\noutdir:", outdir + , "\nplotdir:", plotdir) + +cat("Variables imported:" + , "\n=====================" + , "\ndrug:", drug + , "\ngene:", gene + , "\ngene_match:", gene_match + , "\nAngstrom symbol:", angstroms_symbol + , "\nNo. of duplicated muts:", dup_muts_nu + , "\nNA count for ORs:", na_count + , "\nNA count in df2:", na_count_df2 + , "\nNA count in df3:", na_count_df3 + , "\ndr_muts_col:", dr_muts_col + , "\nother_muts_col:", other_muts_col + , "\ndrtype_col:", resistance_col) + + +#=========== +# input +#=========== +# output of combining_dfs_plotting.R + +#======= +# output +#======= +# plot 1 +basic_bp_lineage = "basic_lineage_barplot.svg" +plot_basic_bp_lineage = paste0(plotdir,"/", basic_bp_lineage) + +#======================================================================= +#================ +# Data for plots: +# you need merged_df2, comprehensive one +# since this has one-many relationship +# i.e the same SNP can belong to multiple lineages +#================ +# REASSIGNMENT as necessary +my_df = merged_df2 + +# clear excess variable +rm(merged_df2_comp, merged_df3, merged_df3_comp) + +# quick checks +colnames(my_df) +str(my_df) + +# Ensure correct data type in columns to plot: need to be factor +is.factor(my_df$lineage) +my_df$lineage = as.factor(my_df$lineage) +is.factor(my_df$lineage) + +#========================== +# Plot: Lineage barplot +# x = lineage y = No. of samples +# col = Lineage +# fill = lineage +#============================ +table(my_df$lineage) +as.data.frame(table(my_df$lineage)) + +#============= +# Data for plots +#============= +# REASSIGNMENT +df <- my_df + +rm(my_df) + +# get freq count of positions so you can subset freq<1 +#setDT(df)[, lineage_count := .N, by = .(lineage)] + +#****************** +# generate plot: barplot of mutation by lineage +#****************** +sel_lineages = c("lineage1" + , "lineage2" + , "lineage3" + , "lineage4" + #, "lineage5" + #, "lineage6" + #, "lineage7" + ) + +df_lin = subset(df, subset = lineage %in% sel_lineages) + +# Create df with lineage inform & no. of unique mutations +# per lineage and total samples within lineage +# this is essentially barplot with two y axis + +bar = bar = as.data.frame(sel_lineages) #4, 1 +total_snps_u = NULL +total_samples = NULL + +for (i in sel_lineages){ + #print(i) + curr_total = length(unique(df$id)[df$lineage==i]) + total_samples = c(total_samples, curr_total) + print(total_samples) + + foo = df[df$lineage==i,] + print(paste0(i, "=======")) + print(length(unique(foo$mutationinformation))) + curr_count = length(unique(foo$mutationinformation)) + + total_snps_u = c(total_snps_u, curr_count) +} + +print(total_snps_u) +bar$num_snps_u = total_snps_u +bar$total_samples = total_samples +bar + +#***************** +# generate plot: lineage barplot with two y-axis +#https://stackoverflow.com/questions/13035295/overlay-bar-graphs-in-ggplot2 +#***************** + +y1 = bar$num_snps_u +y2 = bar$total_samples +x = sel_lineages + +to_plot = data.frame(x = x + , y1 = y1 + , y2 = y2) +to_plot + +# FIXME later: will be depricated! +melted = melt(to_plot, id = "x") +melted + + +svg(plot_basic_bp_lineage) + +my_ats = 20 # axis text size +my_als = 22 # axis label size + +g = ggplot(melted, aes(x = x + , y = value + , fill = variable)) + +printFile = g + geom_bar(stat = "identity" + , position = position_stack(reverse = TRUE) + , alpha=.75 + , colour='grey75') + + theme(axis.text.x = element_text(size = my_ats) + , axis.text.y = element_text(size = my_ats + #, angle = 30 + , hjust = 1 + , vjust = 0) + , axis.title.x = element_text(size = my_als + , colour = 'black') + , axis.title.y = element_text(size = my_als + , colour = 'black') + , legend.position = "top" + , legend.text = element_text(size = my_als)) + + #geom_text() + + geom_label(aes(label = value) + , size = 5 + , hjust = 0.5 + , vjust = 0.5 + , colour = 'black' + , show.legend = FALSE + #, check_overlap = TRUE + , position = position_stack(reverse = T)) + + labs(title = '' + , x = '' + , y = "Number" + , fill = 'Variable' + , colour = 'black') + + scale_fill_manual(values = c('grey50', 'gray75') + , name='' + , labels=c('Mutations', 'Total Samples')) + + scale_x_discrete(breaks = c('lineage1', 'lineage2', 'lineage3', 'lineage4') + , labels = c('Lineage 1', 'Lineage 2', 'Lineage 3', 'Lineage 4')) + +print(printFile) +dev.off() diff --git a/scripts/plotting/redundant/other_plots_data_SAFEGUARD.R b/scripts/plotting/redundant/other_plots_data_SAFEGUARD.R new file mode 100644 index 0000000..df5c1e3 --- /dev/null +++ b/scripts/plotting/redundant/other_plots_data_SAFEGUARD.R @@ -0,0 +1,301 @@ +#!/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) +library(tidyverse) +source("combining_dfs_plotting.R") + +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) + + +cols_to_select = c("mutation", "mutationinformation" + , "wild_type", "position", "mutant_type" + , "mutation_info") + +merged_df3_short = merged_df3[, cols_to_select] + +# write merged_df3 to generate structural figure +write.csv(merged_df3_short, "merged_df3_short.csv") + +#======================================================================== +#%%%%%%%%%%%%%%%%%%% +# REASSIGNMENT: PS +#%%%%%%%%%%%%%%%%%%%% +df_ps = merged_df3 + +#============================ +# adding foldx scaled values +# scale data b/w -1 and 1 +#============================ +n = which(colnames(df_ps) == "ddg"); n + +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) +# sanity check +my_min = min(df_ps$foldx_scaled); my_min +my_max = max(df_ps$foldx_scaled); my_max + +if (my_min == -1 && my_max == 1){ + cat("PASS: foldx ddg successfully scaled b/w -1 and 1" + , "\nProceeding with assigning foldx outcome category") +}else{ + cat("FAIL: could not scale foldx ddg values" + , "Aborting!") +} + +#================================ +# adding foldx outcome category +# 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() +} +#======================================================================= +# 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) + + # further checks to make sure dr and other muts are indeed unique +dr_muts = df_ps[df_ps$mutation_info == dr_muts_col,] +dr_muts_names = unique(dr_muts$mutation) + +other_muts = df_ps[df_ps$mutation_info == other_muts_col,] +other_muts_names = unique(other_muts$mutation) + +if ( table(dr_muts_names%in%other_muts_names)[[1]] == length(dr_muts_names) && + table(other_muts_names%in%dr_muts_names)[[1]] == length(other_muts_names) ){ + cat("PASS: dr and other muts are indeed unique") +}else{ + cat("FAIL: dr adn others muts are NOT unique!") + quit() +} + + +#%%%%%%%%%%%%%%%%%%% +# 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) + +#======================================================================== +#=========== +# 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) + +ligand_dist_colname = paste0("Distance to ligand (", angstroms_symbol, ")") +ligand_dist_colname + +duet_stability_name = paste0(delta_symbol, delta_symbol, "G") +duet_stability_name + +#levels(df_lf_ps$param_type)[levels(df_lf_ps$param_type)=="duet_scaled"] <- "Stability" +levels(df_lf_ps$param_type)[levels(df_lf_ps$param_type)=="duet_scaled"] <- duet_stability_name +#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)=="ligand_distance"] <- ligand_dist_colname +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) + +############################################################################ + +#=========== +# 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) + +############################################################################# +#=========== +# 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, 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() +} + +foldx_stability_name = paste0(delta_symbol, delta_symbol, "G") +foldx_stability_name + +# 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"] <- "Stability" +levels(df_lf_foldx$param_type)[levels(df_lf_foldx$param_type)=="foldx_scaled"] <- foldx_stability_name +#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) + +############################################################################ + +# 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)