moved all test scripts for functions to tests/

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"
+              
+)

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")
@@ -141,3 +165,15 @@ if(!require(protr)){
 
 #BiocManager::install("Logolas")
 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!
+

scripts/plotting/get_plotting_dfs.R

@@ -39,8 +39,8 @@ import_dirs(drug, gene)
 #---------------------------
 # call: plotting_data()
 #---------------------------
-#if (!exists("infile_params") && exists("gene")){
+if (!exists("infile_params") && exists("gene")){
-if (!is.character(infile_params) && exists("gene")){ # when running as cmd
+#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 (!exists("infile_metadata") && exists("gene")){
-if (!is.character(infile_metadata) && exists("gene")){ # when running as cmd
+#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")

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)
+lin_lf_plot$sel_lineages = factor(lin_lf_plot$sel_lineages)
-levels(lin_lf_plot$sel_lineages_f)
+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)
+lin_wf_plot$sel_lineages = factor(lin_wf_plot$sel_lineages)
-levels(lin_wf_plot$sel_lineages_f)
+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"

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,12 +101,8 @@ 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
+# Important: Check factors so that x-axis categ appear as you want
-lin_wf$sel_lineages_f = gsub("lineage", "L", lin_wf$sel_lineages)
+lin_wf$sel_lineages =  factor(lin_wf$sel_lineages, c("L1"
-lin_wf
-
-# 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"
@@ -96,7 +123,7 @@ lin_wf$sel_lineages_f =  factor(lin_wf$sel_lineages_f, c("L1"
                                                      , "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,7 +156,7 @@ 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"
+lin_lf$sel_lineages =  factor(lin_lf$sel_lineages, c("L1"
                                                      , "L2"
                                                      , "L3"
                                                      , "L4"
@@ -150,4 +177,4 @@ lin_lf$sel_lineages_f =  factor(lin_lf$sel_lineages_f, c("L1"
                                                      , "L4;L7"
                                                      , ""))
 
-levels(lin_lf$sel_lineages_f)
+levels(lin_lf$sel_lineages)

scripts/plotting/lineage_dist_combined_PS.R

@@ -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()

scripts/plotting/lineage_dist_dm_om_combined_PS.R

@@ -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()
-
-##############################################################################