added untracked files in scripts/plotting

scripts/plotting/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()

scripts/plotting/lineage_dist_combined_PS.R

@@ -0,0 +1,303 @@
+#!/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()

scripts/plotting/lineage_dist_dm_om_combined_PS.R

@@ -0,0 +1,387 @@
+#!/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()
+
+##############################################################################

scripts/plotting/other_plots_data.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)