From 03031d2eb6ddd2d27cc142375bc3be90d2f735da Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Thu, 9 Sep 2021 13:12:07 +0100 Subject: [PATCH] moved all test scripts for functions to tests/ --- .../functions/{ => tests}/test_aa_prop_bp.R | 0 .../functions/{ => tests}/test_af_or_calcs.R | 0 scripts/functions/{ => tests}/test_bp.R | 0 .../functions/{ => tests}/test_bp_lineage.R | 0 .../{ => tests}/test_combining_dfs_plotting.R | 0 scripts/functions/{ => tests}/test_lf_bp.R | 0 .../{ => tests}/test_lf_unpaired_stats.R | 0 scripts/functions/tests/test_lineage_dist.R | 32 ++ .../{ => tests}/test_plotting_data.R | 0 scripts/plotting/Header_TT.R | 40 +- scripts/plotting/get_plotting_dfs.R | 8 +- .../lineage_basic_barplots_combined.R | 39 +- scripts/plotting/lineage_bp_data.R | 129 +++--- scripts/plotting/lineage_dist_combined_PS.R | 303 -------------- .../plotting/lineage_dist_dm_om_combined_PS.R | 387 ------------------ 15 files changed, 162 insertions(+), 776 deletions(-) rename scripts/functions/{ => tests}/test_aa_prop_bp.R (100%) rename scripts/functions/{ => tests}/test_af_or_calcs.R (100%) rename scripts/functions/{ => tests}/test_bp.R (100%) rename scripts/functions/{ => tests}/test_bp_lineage.R (100%) rename scripts/functions/{ => tests}/test_combining_dfs_plotting.R (100%) rename scripts/functions/{ => tests}/test_lf_bp.R (100%) rename scripts/functions/{ => tests}/test_lf_unpaired_stats.R (100%) create mode 100644 scripts/functions/tests/test_lineage_dist.R rename scripts/functions/{ => tests}/test_plotting_data.R (100%) delete mode 100755 scripts/plotting/lineage_dist_combined_PS.R delete mode 100755 scripts/plotting/lineage_dist_dm_om_combined_PS.R diff --git a/scripts/functions/test_aa_prop_bp.R b/scripts/functions/tests/test_aa_prop_bp.R similarity index 100% rename from scripts/functions/test_aa_prop_bp.R rename to scripts/functions/tests/test_aa_prop_bp.R diff --git a/scripts/functions/test_af_or_calcs.R b/scripts/functions/tests/test_af_or_calcs.R similarity index 100% rename from scripts/functions/test_af_or_calcs.R rename to scripts/functions/tests/test_af_or_calcs.R diff --git a/scripts/functions/test_bp.R b/scripts/functions/tests/test_bp.R similarity index 100% rename from scripts/functions/test_bp.R rename to scripts/functions/tests/test_bp.R diff --git a/scripts/functions/test_bp_lineage.R b/scripts/functions/tests/test_bp_lineage.R similarity index 100% rename from scripts/functions/test_bp_lineage.R rename to scripts/functions/tests/test_bp_lineage.R diff --git a/scripts/functions/test_combining_dfs_plotting.R b/scripts/functions/tests/test_combining_dfs_plotting.R similarity index 100% rename from scripts/functions/test_combining_dfs_plotting.R rename to scripts/functions/tests/test_combining_dfs_plotting.R diff --git a/scripts/functions/test_lf_bp.R b/scripts/functions/tests/test_lf_bp.R similarity index 100% rename from scripts/functions/test_lf_bp.R rename to scripts/functions/tests/test_lf_bp.R diff --git a/scripts/functions/test_lf_unpaired_stats.R b/scripts/functions/tests/test_lf_unpaired_stats.R similarity index 100% rename from scripts/functions/test_lf_unpaired_stats.R rename to scripts/functions/tests/test_lf_unpaired_stats.R diff --git a/scripts/functions/tests/test_lineage_dist.R b/scripts/functions/tests/test_lineage_dist.R new file mode 100644 index 0000000..1f40d16 --- /dev/null +++ b/scripts/functions/tests/test_lineage_dist.R @@ -0,0 +1,32 @@ +############################### +# TEST function lineage_dist.R +# to plot lineage +# dist plots with or without facet +############################## +getwd() +setwd("~/git/LSHTM_analysis/scripts/plotting/") +getwd() + +source("Header_TT.R") + +source("get_plotting_dfs.R") + +cat("cols imported:" + , mcsm_red2, mcsm_red1, mcsm_mid, mcsm_blue1, mcsm_blue2) + + +############################################################# + +lineage_distP(lin_dist_plot + , with_facet = F + , leg_label = "Mutation Class" +) + +lineage_distP(lin_dist_plot + , with_facet = T + , facet_wrap_var = "mutation_info_labels" + , leg_label = "Mutation Class" + , leg_pos_wf = "none" + , leg_dir_wf = "horizontal" + +) \ No newline at end of file diff --git a/scripts/functions/test_plotting_data.R b/scripts/functions/tests/test_plotting_data.R similarity index 100% rename from scripts/functions/test_plotting_data.R rename to scripts/functions/tests/test_plotting_data.R diff --git a/scripts/plotting/Header_TT.R b/scripts/plotting/Header_TT.R index e4593d0..47599d3 100755 --- a/scripts/plotting/Header_TT.R +++ b/scripts/plotting/Header_TT.R @@ -1,8 +1,12 @@ ######################################################### -### A) Installing and loading required packages +# A) Installing and loading required packages +# B) My functions +######################################################### + ######################################################### #lib_loc = "/usr/local/lib/R/site-library") + require("getopt", quietly = TRUE) # cmd parse arguments if (!require("tidyverse")) { @@ -10,6 +14,21 @@ if (!require("tidyverse")) { library(tidyverse) } +if (!require("shiny")) { + install.packages("shiny", dependencies = TRUE) + library(shiny) +} + +if (!require("gridExtra")) { + install.packages("gridExtra", dependencies = TRUE) + library(gridExtra) +} + +if (!require("ggridges")) { + install.packages("ggridges", dependencies = TRUE) + library(ggridges) +} + # if (!require("ggplot2")) { # install.packages("ggplot2", dependencies = TRUE) # library(ggplot2) @@ -20,6 +39,11 @@ if (!require("tidyverse")) { # library(dplyr) # } +if (!require ("plyr")){ + install.packages("plyr") + library(plyr) + } + # Install #if(!require(devtools)) install.packages("devtools") #devtools::install_github("kassambara/ggcorrplot") @@ -140,4 +164,16 @@ if(!require(protr)){ # install.packages("BiocManager") #BiocManager::install("Logolas") -library("Logolas") \ No newline at end of file +library("Logolas") + + +#################################### +# Load all my functions: +# only works if tidyverse is loaded +# hence included it here! +#################################### + +func_path = "~/git/LSHTM_analysis/scripts/functions/" +source_files <- list.files(func_path, "\\.R$") # locate all .R files +map(paste0(func_path, source_files), source) # source all your R scripts! + diff --git a/scripts/plotting/get_plotting_dfs.R b/scripts/plotting/get_plotting_dfs.R index 5876d8d..89b477c 100755 --- a/scripts/plotting/get_plotting_dfs.R +++ b/scripts/plotting/get_plotting_dfs.R @@ -39,8 +39,8 @@ import_dirs(drug, gene) #--------------------------- # call: plotting_data() #--------------------------- -#if (!exists("infile_params") && exists("gene")){ -if (!is.character(infile_params) && exists("gene")){ # when running as cmd +if (!exists("infile_params") && exists("gene")){ +#if (!is.character(infile_params) && exists("gene")){ # when running as cmd #in_filename_params = paste0(tolower(gene), "_all_params.csv") #for pncA in_filename_params = paste0(tolower(gene), "_comb_afor.csv") # part combined for gid infile_params = paste0(outdir, "/", in_filename_params) @@ -67,8 +67,8 @@ cat("\nLigand distance cut off, colname:", LigDist_colname #-------------------------------- # call: combining_dfs_plotting() #-------------------------------- -#if (!exists("infile_metadata") && exists("gene")){ -if (!is.character(infile_metadata) && exists("gene")){ # when running as cmd +if (!exists("infile_metadata") && exists("gene")){ +#if (!is.character(infile_metadata) && exists("gene")){ # when running as cmd in_filename_metadata = paste0(tolower(gene), "_metadata.csv") # part combined for gid infile_metadata = paste0(outdir, "/", in_filename_metadata) cat("\nInput file for gene metadata not specified, assuming filename: ", infile_metadata, "\n") diff --git a/scripts/plotting/lineage_basic_barplots_combined.R b/scripts/plotting/lineage_basic_barplots_combined.R index 94bcd4a..b6f25e6 100755 --- a/scripts/plotting/lineage_basic_barplots_combined.R +++ b/scripts/plotting/lineage_basic_barplots_combined.R @@ -12,7 +12,6 @@ getwd() # Installing and loading required packages ########################################################## source("Header_TT.R") -source("../functions/bp_lineage.R") #=========== # input @@ -40,24 +39,6 @@ if(is.null(drug)|is.null(gene)) { source ('get_plotting_dfs.R') -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 - , "\ndr_muts_col:", dr_muts_col - , "\nother_muts_col:", other_muts_col - , "\ndrtype_col:", resistance_col) - #======= # output #======= @@ -74,21 +55,21 @@ plot_basic_bp_lineage_cl # Data: All lineages or # selected few #------------------------ -sel_lineages = levels(lin_lf$sel_lineages_f)[1:4] +sel_lineages = levels(lin_lf$sel_lineages)[1:4] sel_lineages -lin_lf_plot = lin_lf[lin_lf$sel_lineages_f%in%sel_lineages,] +lin_lf_plot = lin_lf[lin_lf$sel_lineages%in%sel_lineages,] str(lin_lf_plot) # drop unused factor levels -lin_lf_plot$sel_lineages_f = factor(lin_lf_plot$sel_lineages_f) -levels(lin_lf_plot$sel_lineages_f) +lin_lf_plot$sel_lineages = factor(lin_lf_plot$sel_lineages) +levels(lin_lf_plot$sel_lineages) str(lin_lf_plot) #------------------------ # plot from my function: #------------------------ lin_countP = lin_count_bp(lin_lf_plot - , x_categ = "sel_lineages_f" + , x_categ = "sel_lineages" , y_count = "p_count" , bar_fill_categ = "count_categ" , display_label_col = "p_count" @@ -109,21 +90,21 @@ lin_countP # Data: All lineages or # selected few #------------------------ -sel_lineages = levels(lin_wf$sel_lineages_f)[1:4] +sel_lineages = levels(lin_wf$sel_lineages)[1:4] sel_lineages -lin_wf_plot = lin_wf[lin_wf$sel_lineages_f%in%sel_lineages,] +lin_wf_plot = lin_wf[lin_wf$sel_lineages%in%sel_lineages,] str(lin_wf_plot) # drop unused factor levels -lin_wf_plot$sel_lineages_f = factor(lin_wf_plot$sel_lineages_f) -levels(lin_wf_plot$sel_lineages_f) +lin_wf_plot$sel_lineages = factor(lin_wf_plot$sel_lineages) +levels(lin_wf_plot$sel_lineages) str(lin_wf_plot) #------------------------ # plot from my function: #------------------------ lin_diversityP = lin_count_bp(lin_wf_plot - , x_categ = "sel_lineages_f" + , x_categ = "sel_lineages" , y_count = "snp_diversity" , display_label_col = "snp_diversity_f" , bar_stat_stype = "identity" diff --git a/scripts/plotting/lineage_bp_data.R b/scripts/plotting/lineage_bp_data.R index 2d51ab0..e9ab929 100755 --- a/scripts/plotting/lineage_bp_data.R +++ b/scripts/plotting/lineage_bp_data.R @@ -27,13 +27,44 @@ cat("\nMissing samples with lineage classification:", table(merged_df2$lineage = } +# Add pretty lineage labels and mut_info_labels +class(merged_df2$lineage); table(merged_df2$lineage) +merged_df2$lineage_labels = gsub("lineage", "L", merged_df2$lineage) +table(merged_df2$lineage_labels) + +class(merged_df2$lineage_labels) + +merged_df2$lineage_labels = factor(merged_df2$lineage_labels, c("L1" + , "L2" + , "L3" + , "L4" + , "L5" + , "L6" + , "L7" + , "LBOV" + , "L1;L2" + , "L1;L3" + , "L1;L4" + , "L2;L3" + , "L2;L3;L4" + , "L2;L4" + , "L2;L6" + , "L2;LBOV" + , "L3;L4" + , "L4;L6" + , "L4;L7" + , "")) + +class(merged_df2$lineage_labels); nlevels(merged_df2$lineage_labels) + + ################################## # WF data: lineages with # snp count # total_samples # snp diversity (perc) ################################## -sel_lineages = levels(as.factor(merged_df2$lineage)) +sel_lineages = levels(merged_df2$lineage_labels) lin_wf = data.frame(sel_lineages) #4, 1 total_snps_u = NULL @@ -41,12 +72,12 @@ total_samples = NULL for (i in sel_lineages){ #print(i) - curr_total = length(unique(merged_df2$id)[merged_df2$lineage==i]) + curr_total = length(unique(merged_df2$id)[merged_df2$lineage_labels==i]) #print(curr_total) total_samples = c(total_samples, curr_total) print(total_samples) - foo = merged_df2[merged_df2$lineage==i,] + foo = merged_df2[merged_df2$lineage_labels==i,] print(paste0(i, "=======\n")) print(length(unique(foo$mutationinformation))) curr_count = length(unique(foo$mutationinformation)) @@ -70,33 +101,29 @@ lin_wf lin_wf$snp_diversity_f = round( (lin_wf$snp_diversity * 100), digits = 0) lin_wf$snp_diversity_f = paste0(lin_wf$snp_diversity_f, "%") -# Lineage names -lin_wf$sel_lineages_f = gsub("lineage", "L", lin_wf$sel_lineages) -lin_wf +# Important: Check factors so that x-axis categ appear as you want +lin_wf$sel_lineages = factor(lin_wf$sel_lineages, c("L1" + , "L2" + , "L3" + , "L4" + , "L5" + , "L6" + , "L7" + , "LBOV" + , "L1;L2" + , "L1;L3" + , "L1;L4" + , "L2;L3" + , "L2;L3;L4" + , "L2;L4" + , "L2;L6" + , "L2;LBOV" + , "L3;L4" + , "L4;L6" + , "L4;L7" + , "")) -# Important: Relevel factors so that x-axis categ appear as you want -lin_wf$sel_lineages_f = factor(lin_wf$sel_lineages_f, c("L1" - , "L2" - , "L3" - , "L4" - , "L5" - , "L6" - , "L7" - , "LBOV" - , "L1;L2" - , "L1;L3" - , "L1;L4" - , "L2;L3" - , "L2;L3;L4" - , "L2;L4" - , "L2;L6" - , "L2;LBOV" - , "L3;L4" - , "L4;L6" - , "L4;L7" - , "")) - -levels(lin_wf$sel_lineages_f) +levels(lin_wf$sel_lineages) ################################## # LF data: lineages with @@ -106,7 +133,7 @@ levels(lin_wf$sel_lineages_f) ################################## names(lin_wf) tot_cols = ncol(lin_wf) -pivot_cols = c("sel_lineages", "sel_lineages_f", "snp_diversity", "snp_diversity_f") +pivot_cols = c("sel_lineages", "snp_diversity", "snp_diversity_f") pivot_cols_n = length(pivot_cols) expected_rows = nrow(lin_wf) * ( length(lin_wf) - pivot_cols_n ) @@ -129,25 +156,25 @@ if ( nrow(lin_lf) == expected_rows ){ } # Important: Relevel factors so that x-axis categ appear as you want -lin_lf$sel_lineages_f = factor(lin_lf$sel_lineages_f, c("L1" - , "L2" - , "L3" - , "L4" - , "L5" - , "L6" - , "L7" - , "LBOV" - , "L1;L2" - , "L1;L3" - , "L1;L4" - , "L2;L3" - , "L2;L3;L4" - , "L2;L4" - , "L2;L6" - , "L2;LBOV" - , "L3;L4" - , "L4;L6" - , "L4;L7" - , "")) +lin_lf$sel_lineages = factor(lin_lf$sel_lineages, c("L1" + , "L2" + , "L3" + , "L4" + , "L5" + , "L6" + , "L7" + , "LBOV" + , "L1;L2" + , "L1;L3" + , "L1;L4" + , "L2;L3" + , "L2;L3;L4" + , "L2;L4" + , "L2;L6" + , "L2;LBOV" + , "L3;L4" + , "L4;L6" + , "L4;L7" + , "")) -levels(lin_lf$sel_lineages_f) +levels(lin_lf$sel_lineages) diff --git a/scripts/plotting/lineage_dist_combined_PS.R b/scripts/plotting/lineage_dist_combined_PS.R deleted file mode 100755 index bf1c75b..0000000 --- a/scripts/plotting/lineage_dist_combined_PS.R +++ /dev/null @@ -1,303 +0,0 @@ -#!/usr/bin/env Rscript -######################################################### -# TASK: Lineage dist plots: ggridges - -# Output: 2 SVGs for PS stability - -# 1) all muts -# 2) dr_muts - -########################################################## -# Installing and loading required packages -########################################################## -getwd() -setwd("~/git/LSHTM_analysis/scripts/plotting/") -getwd() - -source("Header_TT.R") -library(ggridges) -source("combining_dfs_plotting.R") -# 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) - -#======= -# output -#======= -lineage_dist_combined = "lineage_dist_combined_PS.svg" -plot_lineage_dist_combined = paste0(plotdir,"/", lineage_dist_combined) -#======================================================================== - -########################### -# Data for plots -# you need merged_df2 or merged_df2_comp -# since this is one-many relationship -# i.e the same SNP can belong to multiple lineages -# using the _comp dataset means -# we lose some muts and at this level, we should use -# as much info as available, hence use df with NA -########################### -# REASSIGNMENT -my_df = merged_df2 - -# delete variables not required -rm(my_df_u, merged_df2, 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) - -table(my_df$mutation_info) - -# subset df with dr muts only -my_df_dr = subset(my_df, mutation_info == "dr_mutations_pyrazinamide") -table(my_df_dr$mutation_info) - -######################################################################## -# end of data extraction and cleaning for plots # -######################################################################## - -#========================== -# Plot 1: ALL Muts -# x = mcsm_values, y = dist -# fill = stability -#============================ - -my_plot_name = 'lineage_dist_PS.svg' - -plot_lineage_duet = paste0(plotdir,"/", my_plot_name) - -#=================== -# Data for plots -#=================== -table(my_df$lineage); str(my_df$lineage) - -# subset only lineages1-4 -sel_lineages = c("lineage1" - , "lineage2" - , "lineage3" - , "lineage4" - #, "lineage5" - #, "lineage6" - #, "lineage7" - ) - -# uncomment as necessary -df_lin = subset(my_df, subset = lineage %in% sel_lineages ) -table(df_lin$lineage) - -# refactor -df_lin$lineage = factor(df_lin$lineage) - -sum(table(df_lin$lineage)) #{RESULT: Total number of samples for lineage} - -table(df_lin$lineage)#{RESULT: No of samples within lineage} - -length(unique(df_lin$mutationinformation))#{Result: No. of unique mutations the 4 lineages contribute to} - -length(df_lin$mutationinformation) - -u2 = unique(my_df$mutationinformation) -u = unique(df_lin$mutationinformation) -check = u2[!u2%in%u]; print(check) #{Muts not present within selected lineages} - -#%%%%%%%%%%%%%%%%%%%%%%%%% -# REASSIGNMENT -df <- df_lin -#%%%%%%%%%%%%%%%%%%%%%%%%% - -rm(df_lin) - -#****************** -# generate distribution plot of lineages -#****************** -# 2 : ggridges (good!) -my_ats = 15 # axis text size -my_als = 20 # axis label size - -my_labels = c('Lineage 1', 'Lineage 2', 'Lineage 3', 'Lineage 4' - #, 'Lineage 5', 'Lineage 6', 'Lineage 7' - ) -names(my_labels) = c('lineage1', 'lineage2', 'lineage3', 'lineage4' - # , 'lineage5', 'lineage6', 'lineage7' - ) -# check plot name -plot_lineage_duet - -# output svg -#svg(plot_lineage_duet) -p1 = ggplot(df, aes(x = duet_scaled - , y = duet_outcome))+ - - #printFile=geom_density_ridges_gradient( - geom_density_ridges_gradient(aes(fill = ..x..) - #, jittered_points = TRUE - , scale = 3 - , size = 0.3 ) + - facet_wrap( ~lineage - , scales = "free" - #, switch = 'x' - , labeller = labeller(lineage = my_labels) ) + - coord_cartesian( xlim = c(-1, 1)) + - scale_fill_gradientn(colours = c("#f8766d", "white", "#00bfc4") - , name = "DUET" ) + - theme(axis.text.x = element_text(size = my_ats - , angle = 90 - , hjust = 1 - , vjust = 0.4) - - , axis.text.y = element_blank() - , axis.title.x = element_blank() - , axis.title.y = element_blank() - , axis.ticks.y = element_blank() - , plot.title = element_blank() - , strip.text = element_text(size = my_als) - , legend.text = element_text(size = my_als-5) - , legend.title = element_text(size = my_als) -) - -print(p1) -#dev.off() - -####################################################################### -# lineage distribution plot for dr_muts -####################################################################### - -#========================== -# Plot 2: dr muts ONLY -# x = mcsm_values, y = dist -# fill = stability -#============================ - -my_plot_name_dr = 'lineage_dist_dr_muts_PS.svg' - -plot_lineage_dr_duet = paste0(plotdir,"/", my_plot_name_dr) - -#=================== -# Data for plots -#=================== -table(my_df_dr$lineage); str(my_df_dr$lineage) - -# uncomment as necessary -df_lin_dr = subset(my_df_dr, subset = lineage %in% sel_lineages) -table(df_lin_dr$lineage) - -# refactor -df_lin_dr$lineage = factor(df_lin_dr$lineage) - -sum(table(df_lin_dr$lineage)) #{RESULT: Total number of samples for lineage} - -table(df_lin_dr$lineage)#{RESULT: No of samples within lineage} - -length(unique(df_lin_dr$mutationinformation))#{Result: No. of unique mutations the 4 lineages contribute to} - -length(df_lin_dr$mutationinformation) - -u2 = unique(my_df_dr$mutationinformation) -u = unique(df_lin_dr$mutationinformation) -check = u2[!u2%in%u]; print(check) #{Muts not present within selected lineages} - -#%%%%%%%%%%%%%%%%%%%%%%%%% -# REASSIGNMENT -df_dr <- df_lin_dr -#%%%%%%%%%%%%%%%%%%%%%%%%% - -rm(df_lin_dr) - -#****************** -# generate distribution plot of lineages -#****************** -# 2 : ggridges (good!) -my_ats = 15 # axis text size -my_als = 20 # axis label size - - -# check plot name -plot_lineage_dr_duet - -# output svg -#svg(plot_lineage_dr_duet) -p2 = ggplot(df_dr, aes(x = duet_scaled - , y = duet_outcome))+ - - geom_density_ridges_gradient(aes(fill = ..x..) - #, jittered_points = TRUE - , scale = 3 - , size = 0.3) + - #geom_point(aes(size = or_mychisq))+ - facet_wrap( ~lineage - , scales = "free" - #, switch = 'x' - , labeller = labeller(lineage = my_labels) ) + - coord_cartesian( xlim = c(-1, 1) - #, ylim = c(0, 6) - #, clip = "off" - ) + - scale_fill_gradientn(colours = c("#f8766d", "white", "#00bfc4") - , name = "DUET" ) + - theme(axis.text.x = element_text(size = my_ats - , angle = 90 - , hjust = 1 - , vjust = 0.4) - , axis.text.y = element_blank() - , axis.title.x = element_blank() - , axis.title.y = element_blank() - , axis.ticks.y = element_blank() - , plot.title = element_blank() - , strip.text = element_text(size = my_als) - , legend.text = element_text(size = 10) - , legend.title = element_text(size = my_als) - #, legend.position = "none" - ) - -print(p2) -#dev.off() -######################################################################## -#============== -# combine plot -#=============== - -svg(plot_lineage_dist_combined, width = 12, height = 6) - -printFile = cowplot::plot_grid(p1, p2 - , label_size = my_als+10) - -print(printFile) -dev.off() diff --git a/scripts/plotting/lineage_dist_dm_om_combined_PS.R b/scripts/plotting/lineage_dist_dm_om_combined_PS.R deleted file mode 100755 index 07912ac..0000000 --- a/scripts/plotting/lineage_dist_dm_om_combined_PS.R +++ /dev/null @@ -1,387 +0,0 @@ -#!/usr/bin/env Rscript -######################################################### -# TASK: Lineage dist plots: ggridges - -# Output: 2 SVGs for PS stability - -# 1) all muts -# 2) dr_muts - -########################################################## -# Installing and loading required packages -########################################################## -getwd() -setwd("~/git/LSHTM_analysis/scripts/plotting/") -getwd() - -source("Header_TT.R") -library(ggridges) -library(plyr) -source("combining_dfs_plotting.R") -# 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) - -cat("cols imported:" - , mcsm_red2, mcsm_red1, mcsm_mid, mcsm_blue1, mcsm_blue2) - -#======= -# output -#======= -lineage_dist_combined_dm_om = "lineage_dist_combined_dm_om_PS.svg" -plot_lineage_dist_combined_dm_om = paste0(plotdir,"/", lineage_dist_combined_dm_om) - -lineage_dist_combined_dm_om_L = "lineage_dist_combined_dm_om_PS_labelled.svg" -plot_lineage_dist_combined_dm_om_L = paste0(plotdir,"/", lineage_dist_combined_dm_om_L) - -#======================================================================== - -########################### -# Data for plots -# you need merged_df2 or merged_df2_comp -# since this is one-many relationship -# i.e the same SNP can belong to multiple lineages -# using the _comp dataset means -# we lose some muts and at this level, we should use -# as much info as available, hence use df with NA -########################### -# REASSIGNMENT -my_df = merged_df2 - -# delete variables not required -rm(my_df_u, merged_df2, merged_df2_comp, merged_df3, merged_df3_comp - , merged_df2_lig, merged_df2_comp_lig, merged_df3_lig, merged_df3_comp_lig) - -# quick checks -colnames(my_df) -str(my_df) - -table(my_df$mutation_info) - -#=================== -# Data for plots -#=================== -table(my_df$lineage); str(my_df$lineage) - -# select lineages 1-4 -sel_lineages = c("lineage1" - , "lineage2" - , "lineage3" - , "lineage4") - #, "lineage5" - #, "lineage6" - #, "lineage7") - -# works nicely with facet wrap using labeller, but not otherwise -#my_labels = c('Lineage 1' -# , 'Lineage 2' -# , 'Lineage 3' -# , 'Lineage 4') -# #, 'Lineage 5' -# #, 'Lineage 6' -# #, 'Lineage 7') - -#names(my_labels) = c('lineage1' -# , 'lineage2' -# , 'lineage3' -# , 'lineage4') -# #, 'lineage5' -# #, 'lineage6' -# #, 'lineage7') - -#========================== -# subset selected lineages -#========================== -df_lin = subset(my_df, subset = lineage %in% sel_lineages) -table(df_lin$lineage) - -#{RESULT: Total number of samples for lineage} -sum(table(df_lin$lineage)) - -#{RESULT: No of samples within lineage} -table(df_lin$lineage) - -#{Result: No. of unique mutations the 4 lineages contribute to} -length(unique(df_lin$mutationinformation)) - -u2 = unique(my_df$mutationinformation) -u = unique(df_lin$mutationinformation) - -#{Result:Muts not present within selected lineages} -check = u2[!u2%in%u]; print(check) - -# workaround to make labels appear nicely for in otherwise cases -#================== -# lineage: labels -# from "plyr" -#================== -#{Result:No of samples in selected lineages} -table(df_lin$lineage) - -df_lin$lineage_labels = mapvalues(df_lin$lineage - , from = c("lineage1","lineage2", "lineage3", "lineage4") - , to = c("Lineage 1", "Lineage 2", "Lineage 3", "Lineage 4")) -table(df_lin$lineage_labels) - -table(df_lin$lineage_labels) == table(df_lin$lineage) - -#======================== -# mutation_info: labels -#======================== -#{Result:No of DM and OM muts in selected lineages} -table(df_lin$mutation_info) - -df_lin$mutation_info_labels = ifelse(df_lin$mutation_info == dr_muts_col, "DM", "OM") -table(df_lin$mutation_info_labels) - -table(df_lin$mutation_info) == table(df_lin$mutation_info_labels) - - -#======================== -# duet_outcome: labels -#======================== -#{Result: No. of D and S mutations in selected lineages} -table(df_lin$duet_outcome) - -df_lin$duet_outcome_labels = ifelse(df_lin$duet_outcome == "Destabilising", "D", "S") -table(df_lin$duet_outcome_labels) - -table(df_lin$duet_outcome) == table(df_lin$duet_outcome_labels) - - -#======================= -# subset dr muts only -#======================= -#my_df_dr = subset(df_lin, mutation_info == dr_muts_col) -#table(my_df_dr$mutation_info) -#table(my_df_dr$lineage) - -#========================= -# subset other muts only -#========================= -#my_df_other = subset(df_lin, mutation_info == other_muts_col) -#table(my_df_other$mutation_info) -#table(my_df_other$lineage) - -######################################################################## -# end of data extraction and cleaning for plots # -######################################################################## - -#========================== -# Distribution plots -#============================ - -#%%%%%%%%%%%%%%%%%%%%%%%%% -# REASSIGNMENT -df <- df_lin -#%%%%%%%%%%%%%%%%%%%%%%%%% - -rm(df_lin) - -#****************** -# generate distribution plot of lineages -#****************** -# 2 : ggridges (good!) -my_ats = 15 # axis text size -my_als = 20 # axis label size -n_colours = length(unique(df$duet_scaled)) -my_palette <- colorRampPalette(c(mcsm_red2, mcsm_red1, mcsm_mid, mcsm_blue1, mcsm_blue2))(n = n_colours+1) - -#======================================= -# Plot 1: lineage dist: geom_density_ridges_gradient (allows aesthetics to vary along ridgeline, no alpha setting!) -# else same as geom_density_ridges) -# x = duet_scaled -# y = duet_outcome -# fill = duet_scaled -# Facet: Lineage -#======================================= -# output individual svg -#plot_lineage_dist_duet_f paste0(plotdir,"/", "lineage_dist_duet_f.svg") -#plot_lineage_dist_duet_f -#svg(plot_lineage_dist_duet_f) - -p1 = ggplot(df, aes(x = duet_scaled - , y = duet_outcome))+ - geom_density_ridges_gradient(aes(fill = ..x..) - #, jittered_points = TRUE - , scale = 3 - , size = 0.3 ) + - facet_wrap( ~lineage_labels - # , scales = "free" - # , labeller = labeller(lineage = my_labels) - ) + - coord_cartesian( xlim = c(-1, 1)) + - scale_fill_gradientn(colours = my_palette - , name = "DUET" - #, breaks = c(-1, 0, 1) - #, labels = c(-1,0,1) - #, limits = c(-1,1) - ) + - theme(axis.text.x = element_text(size = my_ats - , angle = 90 - , hjust = 1 - , vjust = 0.4) - #, axis.text.y = element_blank() - , axis.text.y = element_text(size = my_ats) - , axis.title.x = element_text(size = my_ats) - , axis.title.y = element_blank() - , axis.ticks.y = element_blank() - , plot.title = element_blank() - , strip.text = element_text(size = my_als) - , legend.text = element_text(size = my_als-10) - #, legend.title = element_text(size = my_als-6) - , legend.title = element_blank() - , legend.position = c(-0.08, 0.41) - #, legend.direction = "horizontal" - #, legend.position = "left" -)+ - labs(x = "DUET") - -p1 - - -#p1_with_legend = p1 + guides(fill = guide_colourbar(label = FALSE)) - -#======================================= -# Plot 2: lineage dist: geom_density_ridges, allows alpha to be set -# x = duet_scaled -# y = lineage_labels -# fill = mutation_info -# NO FACET -#======================================= -# output svg -#plot_lineage_dist_duet_dm_om = paste0(plotdir,"/", "lineage_dist_duet_dm_om.svg") -#plot_lineage_dist_duet_dm_om -#svg(plot_lineage_dist_duet_dm_om) - -p2 = ggplot(df, aes(x = duet_scaled - , y = lineage_labels))+ - geom_density_ridges(aes(fill = factor(mutation_info_labels)) - , scale = 3 - , size = 0.3 - , alpha = 0.8) + - coord_cartesian( xlim = c(-1, 1)) + - scale_fill_manual(values = c("#E69F00", "#999999")) + - theme(axis.text.x = element_text(size = my_ats - , angle = 90 - , hjust = 1 - , vjust = 0.4) - , axis.text.y = element_text(size = my_ats) - , axis.title.x = element_text(size = my_ats) - , axis.title.y = element_blank() - , axis.ticks.y = element_blank() - , plot.title = element_blank() - , strip.text = element_text(size = my_als) - , legend.text = element_text(size = my_als-4) - , legend.title = element_text(size = my_als-4) - , legend.position = c(0.8, 0.9)) + - labs(x = "DUET" - , fill = "Mutation class") # legend title - -p2 - -#======================================= -# Plot 3: lineage dist: geom_density_ridges_gradient (allows aesthetics to vary along ridgeline, no alpha setting!) -# else same as geom_density_ridges) -# x = duet_scaled -# y = lineage_labels -# fill = duet_scaled -# NO FACET (nf) -#======================================= -# output individual svg -#plot_lineage_dist_duet_nf = paste0(plotdir,"/", "lineage_dist_duet_nf.svg") -#plot_lineage_dist_duet_nf -#svg(plot_lineage_dist_duet_nf) - -p3 = ggplot(df, aes(x = duet_scaled - , y = lineage_labels))+ - geom_density_ridges_gradient(aes(fill = ..x..) - #, jittered_points = TRUE - , scale = 3 - , size = 0.3 ) + - coord_cartesian( xlim = c(-1, 1)) + - scale_fill_gradientn(colours = my_palette, name = "DUET") + - theme(axis.text.x = element_text(size = my_ats - , angle = 90 - , hjust = 1 - , vjust = 0.4) - - , axis.text.y = element_text(size = my_ats) - , axis.title.x = element_text(size = my_ats) - , axis.title.y = element_blank() - , axis.ticks.y = element_blank() - , plot.title = element_blank() - , strip.text = element_text(size = my_als) - , legend.text = element_text(size = my_als-10) - , legend.title = element_text(size = my_als-3) - , legend.position = c(0.8, 0.8)) + - #, legend.direction = "horizontal")+ - #, legend.position = "top")+ - labs(x = "DUET") - -p3 - -######################################################################## -#============== -# combine plots -#=============== -# 1) without labels -plot_lineage_dist_combined_dm_om -svg(plot_lineage_dist_combined_dm_om, width = 12, height = 6) - -OutPlot1 = cowplot::plot_grid(p1, p2 - , rel_widths = c(0.5/2, 0.5/2)) - -print(OutPlot1) -dev.off() - - -# 2) with labels -plot_lineage_dist_combined_dm_om_L -svg(plot_lineage_dist_combined_dm_om_L, width = 12, height = 6) - -OutPlot2 = cowplot::plot_grid(p1, p2 - #, labels = c("(a)", "(b)") - , labels = "AUTO" - #, label_x = -0.045, label_y = 0.92 - #, hjust = -0.7, vjust = -0.5 - #, align = "h" - , rel_widths = c(0.5/2, 0.5/2) - , label_size = my_als) - -print(OutPlot2) -dev.off() - -##############################################################################