modified bp with option for adding stats and boxplplots. Moved old one to redundant

scripts/functions/lf_bp.R

@@ -0,0 +1,193 @@
+#############################
+# Barplots: ggplot
+# stats +/-
+# violin +/-
+# barplot +/
+# beeswarm
+#############################
+
+lf_bp <- function(lf_df
+                  , p_title = ""
+                  , colour_categ = ""
+                  , x_grp = "mutation_info"
+                  , y_var = "param_value"
+                  , facet_var = "param_type"
+                  , n_facet_row = 1
+                  , y_scales = "free_y"
+                  , colour_bp_strip = "khaki2"
+                  , dot_size = 3
+                  , dot_transparency = 0.3
+                  , violin_quantiles = c(0.25, 0.5, 0.75) # can be NULL
+                  , my_ats = 22 # axis text size
+                  , my_als = 20 # axis label size
+                  , my_fls = 20 # facet label size
+                  , my_pts = 22 # plot title size)
+                  , make_boxplot = FALSE
+                  , bp_width = c("auto", 0.5)
+                  , add_stats = FALSE
+                  , stat_grp_comp = c("DM", "OM")
+                  , stat_method = "wilcox.test"
+                  , my_paired = FALSE
+                  , stat_label = c("p.format", "p.signif") ){
+
+  p1 <- ggplot(lf_df, aes(x = eval(parse(text = x_grp))
+                    , y = eval(parse(text = y_var)) ))  + 
+    
+    facet_wrap(~ eval(parse(text = facet_var))
+               , nrow = n_facet_row
+               , scales = y_scales) +
+    
+    geom_violin(trim = T
+                , scale = "width"
+                #, position = position_dodge(width = 0.9)
+                , draw_quantiles = violin_quantiles)
+    
+    if (make_boxplot){
+      
+      if (bp_width == "auto"){
+        bp_width = 0.5/length(unique(lf_df[[x_grp]]))
+        cat("\nAutomatically calculated boxplot width, using bp_width:\n", bp_width, "\n")
+      }else{
+        cat("\nBoxplot width value provided, using:",  bp_width, "\n")
+        bp_width = bp_width} 
+    
+    p2 = p1 + geom_boxplot(fill = "white"
+                  , outlier.colour = NA
+                  #, position = position_dodge(width = 0.9)
+                  , width = bp_width) + 
+      geom_beeswarm(priority = "density"
+                    #, shape = 21
+                    , size = dot_size
+                    , alpha = dot_transparency
+                    , show.legend = FALSE
+                    , cex = 0.8
+                    , aes(colour = factor(eval(parse(text = colour_categ))) ))  
+      
+    } else {
+    # ggbeeswarm (better than geom_point)
+    p2 = p1 +  geom_beeswarm(priority = "density"
+                    #, shape = 21
+                    , size = dot_size
+                    , alpha = dot_transparency
+                    , show.legend = FALSE
+                    , cex = 0.8
+                    , aes(colour = factor(eval(parse(text = colour_categ))) )) 
+    }
+  
+   # Add foramtting to graph
+   OutPlot = p2 + theme(axis.text.x = element_text(size = my_ats)
+                   , axis.text.y = element_text(size = my_ats
+                                       , angle = 0
+                                       , hjust = 1
+                                       , vjust = 0)
+                   , axis.title.x = element_text(size = my_ats)
+                   , axis.title.y = element_text(size = my_ats)
+                   , plot.title = element_text(size = my_pts
+                                      , hjust = 0.5
+                                      , colour = "black"
+                                      , face = "bold")
+                   , strip.background = element_rect(fill = colour_bp_strip)
+                   , strip.text.x = element_text(size = my_fls
+                                                  , colour = "black")
+                   , legend.title = element_text(color = "black"
+                                                  , size = my_als)
+                   , legend.text = element_text(size = my_ats)
+                   , legend.direction = "vertical") +
+    
+    labs(title = p_title
+         , x = ""
+         , y = "") 
+   
+    if (add_stats){
+      my_comparisonsL <- list( stat_grp_comp )
+      
+    OutPlot = OutPlot + stat_compare_means(comparisons = my_comparisonsL
+                       , method = stat_method
+                       , paired = my_paired
+                       , label = stat_label[1])
+    } 
+   return(OutPlot)
+}
+
+#############################
+# Barplot NO stats: plotly
+# violin +/-
+# barplot +/
+# beeswarm
+
+# TODO: plot_ly()
+#############################
+lf_bp_plotly <- function(lf_df
+                         , p_title = ""
+                         , colour_categ = ""
+                         , x_grp = mutation_info
+                         , y_var = param_value
+                         , facet_var = param_type
+                         , n_facet_row = 1
+                         , y_scales = "free_y"
+                         , colour_bp_strip = "khaki2"
+                         , dot_size = 3
+                         , dot_transparency = 0.3
+                         , violin_quantiles = c(0.25, 0.5, 0.75) # can be NULL
+                         , my_ats = 20 # axis text size
+                         , my_als = 18 # axis label size
+                         , my_fls = 18 # facet label size
+                         , my_pts = 22 # plot title size)
+                         #, make_boxplot = FALSE
+                         , bp_width = c("auto", 0.5)
+                         #, add_stats = FALSE
+                         #, stat_grp_comp = c("DM", "OM")
+                         #, stat_method = "wilcox.test"
+                         #, my_paired = FALSE
+                         #, stat_label = c("p.format", "p.signif") 
+                         ){
+  
+  OutPlotly = ggplot(lf_df, aes(x = eval(parse(text = x_grp))
+                              , y = eval(parse(text = y_var))
+                              , label1 = x_grp
+                              , label2 = y_var
+                              , lable3 = colour_categ) )  +
+
+      facet_wrap(~ eval(parse(text = facet_var))
+               , nrow = n_facet_row
+               , scales = y_scales) +
+    
+    geom_violin(trim = T
+                , scale = "width"
+                , draw_quantiles = violin_quantiles) +
+    
+   geom_beeswarm(priority = "density"
+                  , size = dot_size
+                  , alpha = dot_transparency
+                  , show.legend = FALSE
+                  , cex = 0.8
+                  , aes(colour = factor(eval(parse(text = colour_categ) ) ) ) ) +
+    theme(axis.text.x = element_text(size = my_ats)
+          , axis.text.y = element_text(size = my_ats
+                                       , angle = 0
+                                       , hjust = 1
+                                       , vjust = 0)
+          , axis.title.x = element_text(size = my_ats)
+          , axis.title.y = element_text(size = my_ats)
+          , plot.title = element_text(size = my_pts
+                             , hjust = 0.5
+                             , colour = "black"
+                             , face = "bold")
+          , strip.background = element_rect(fill = colour_bp_strip)
+          , strip.text.x = element_text(size = my_fls
+                                         , colour = "black")
+          , legend.title = element_text(color = "black"
+                                         , size = my_als)
+          , legend.text = element_text(size = my_ats)
+          , legend.position = "none")+
+    
+    labs(title = p_title
+         , x = ""
+         , y = "") 
+  
+  OutPlotly = ggplotly(OutPlotly
+                       #, tooltip = c("label")
+                       )
+  return(OutPlotly)
+  
+}

scripts/functions/redundant/lf_bp_with_stats.R

@@ -14,13 +14,23 @@ lf_bp_with_stats <- function(lf_df
                         , stat_grp_comp = c("DM", "OM")
                         , stat_method = "wilcox.test"
                         , my_paired = FALSE
-                        #, stat_label = "p.format")
+                        , bp_width = c("auto", 0.5)
+                        , dot_size = 3
+                        , dot_transparency = 0.3
                         , stat_label = c("p.format", "p.signif")
                         , my_ats = 22 # axis text size
                         , my_als = 20 # axis label size
                         , my_fls = 20 # facet label size
                         , my_pts = 22 # plot title size
 ) {
+  if (bp_width == "auto"){
+  bp_width = 0.5/length(unique(lf_df[[x_grp]]))
+  cat("\nAutomatically calculated boxplot width, using bp_width:\n", bp_width, "\n")
+  }else{
+  cat("\nBoxplot width value provided, using:",  bp_width, "\n")
+  bp_width = bp_width
+  }
+  
   my_comparisonsL <- list( stat_grp_comp )
 
   bp_statP <- ggplot(lf_df, aes(x = eval(parse(text = x_grp))
@@ -30,15 +40,30 @@ lf_bp_with_stats <- function(lf_df
                , nrow = n_facet_row
                , scales = y_scales) +
     
-    geom_boxplot(fill = "white", outlier.colour = NA
+    geom_violin(trim = T
-                 #, position = position_dodge(width = 0.9)
+                , scale = "width"
-                 , width = 0.2) +
+                #, position = position_dodge(width = 0.9)
-    
+                , draw_quantiles = c(0.25, 0.5, 0.75)) + 
-    geom_point(position = position_jitterdodge(dodge.width = 0.01)
-               , alpha = 0.5
-               , show.legend = FALSE
-               , aes(colour = factor(eval(parse(text = colour_categ))) )) +
     
+    # geom_boxplot(fill = "white"
+    #              , outlier.colour = NA
+    #              #, position = position_dodge(width = 0.9)
+    #              , width = bp_width) +
+
+    # geom_point(position = position_jitterdodge(dodge.width = 0.5)
+    #            , alpha = 0.5
+    #            , show.legend = FALSE
+    #            , aes(colour = factor(eval(parse(text = colour_categ))) )) +
+
+    # ggbeeswarm (better than geom_point)
+      geom_beeswarm(priority = "density"
+                    #, shape = 21
+                    , size = dot_size
+                    , alpha = dot_transparency
+                    , show.legend = FALSE
+                    , cex = 0.8
+                    , aes(colour = factor(eval(parse(text = colour_categ))) )) +
+
     theme(axis.text.x = element_text(size = my_ats)
           , axis.text.y = element_text(size = my_ats
                                        , angle = 0
@@ -46,10 +71,15 @@ lf_bp_with_stats <- function(lf_df
                                        , vjust = 0)
           , axis.title.x = element_text(size = my_ats)
           , axis.title.y = element_text(size = my_ats)
-          , plot.title = element_text(size = my_pts , hjust = 0.5, colour = "black", face = "bold")
+          , plot.title = element_text(size = my_pts
+                                      , hjust = 0.5
+                                      , colour = "black"
+                                      , face = "bold")
           , strip.background = element_rect(fill = colour_bp_strip)
-          , strip.text.x = element_text(size = my_fls, colour = "black")
+          , strip.text.x = element_text(size = my_fls
-          , legend.title = element_text(color = "black", size = my_als)
+                                        , colour = "black")
+          , legend.title = element_text(color = "black"
+                                        , size = my_als)
           , legend.text = element_text(size = my_ats)
           , legend.direction = "vertical") +
     
@@ -60,7 +90,6 @@ lf_bp_with_stats <- function(lf_df
     stat_compare_means(comparisons = my_comparisonsL
                        , method = stat_method
                        , paired = my_paired
-                       #, label = "p.format")
                        , label = stat_label[1])
   
   return(bp_statP)

scripts/functions/redundant/test_lf_bp_with_stats.R

@@ -0,0 +1,83 @@
+setwd("~/git/LSHTM_analysis/scripts/plotting/")
+
+source("../functions/lf_bp_with_stats.R")
+source("../functions/lf_bp.R")
+
+######################
+# Make  plot
+######################
+# Note: Data
+# run other_plots_data.R
+# to get the long format data to test this function 
+
+lf_bp(lf_df = lf_dynamut2
+                       , p_title = "Dynamut2"
+                       , colour_categ = "ddg_dynamut2_outcome"
+                       , x_grp = "mutation_info"
+                       , y_var = "param_value"
+                       , facet_var = "param_type"
+                       , n_facet_row = 1
+                       , y_scales = "free_y"
+                       , colour_bp_strip = "khaki2"
+                       , dot_size = 3
+                       , dot_transparency = 0.3
+                       , violin_quantiles = c(0.25, 0.5, 0.75)
+                       , my_ats = 22 # axis text size
+                       , my_als = 20 # axis label size
+                       , my_fls = 20 # facet label size
+                       , my_pts = 22 # plot title size 
+                       , make_boxplot = F
+                       , bp_width = "auto"
+                       , add_stats = T
+                       , stat_grp_comp = c("DM", "OM")
+                       , stat_method = "wilcox.test"
+                       , my_paired = FALSE
+                       , stat_label = c("p.format", "p.signif") )
+ 
+# foo = lf_dynamut2 %>%
+#   group_by(mutation_info, param_type) %>%
+#   summarise( Mean = mean(param_value, na.rm = T)
+#              , SD = sd(param_value, na.rm = T)
+#              , Median = median(param_value, na.rm = T)
+#              , IQR = IQR(param_value, na.rm = T) )
+
+# Quick tests
+plotdata_sel = subset(lf_dynamut2
+             , lf_dynamut2$param_type == "ASA")
+
+plot_sum = plotdata_sel %>%
+  group_by(mutation_info, param_type) %>%
+  summarise(n = n()
+             , Mean = mean(param_value, na.rm = T)
+             , SD = sd(param_value, na.rm = T)
+             , Min = min(param_value, na.rm = T)
+             , Q1 = quantile(param_value, na.rm = T, 0.25)
+             , Median = median(param_value, na.rm = T)
+             , Q3 = quantile(param_value, na.rm = T, 0.75)
+             , Max = max(param_value, na.rm = T) ) %>%
+  rename('Mutation Class' = mutation_info
+         , Parameter = param_type)
+plot_sum = as.data.frame(plot_sum, row.names = NULL)
+plot_sum
+
+bar = compare_means(param_value ~ mutation_info
+              , group.by = "param_type"
+              , data = plotdata_sel
+              , paired = FALSE
+              , p.adjust.method = "BH")
+bar2 = bar[c("param_type"
+              , "group1"
+              , "group2"
+              , "p.format"
+              , "p.signif"
+              , "p.adj")] %>%
+  rename(Parameter = param_type
+          , Group1 = group1
+          , Group2 = group2
+          , "P-value" = p.format
+          , "P-sig" = p.signif
+          , "P-adj" = p.adj)
+bar2 = data.frame(bar2); bar2
+
+library(Hmisc)
+describe(lf_dynamut2)

scripts/functions/test_lf_bp.R

@@ -0,0 +1,55 @@
+setwd("~/git/LSHTM_analysis/scripts/plotting/")
+source("Header_TT.R")
+source("../functions/lf_bp.R")
+# ================================================
+# Data: run get_plotting_data.R
+# to get the long format data to test this function 
+# drug = "streptomycin" 
+# gene = "gid"
+# source("get_plotting_dfs.R")
+# ==================================================
+
+######################
+# Make plot: ggplot
+######################
+lf_bp(lf_df = lf_dynamut2
+      , p_title = "Dynamut2"
+      , colour_categ = "ddg_dynamut2_outcome"
+      , x_grp = "mutation_info"
+      , y_var = "param_value"
+      , facet_var = "param_type"
+      , n_facet_row = 1
+      , y_scales = "free_y"
+      , colour_bp_strip = "khaki2"
+      , dot_size = 3
+      , dot_transparency = 0.3
+      , violin_quantiles = c(0.25, 0.5, 0.75)
+      , my_ats = 22 # axis text size
+      , my_als = 20 # axis label size
+      , my_fls = 20 # facet label size
+      , my_pts = 22 # plot title size 
+      , make_boxplot = F
+      , bp_width = "auto"
+      , add_stats = T
+      , stat_grp_comp = c("DM", "OM")
+      , stat_method = "wilcox.test"
+      , my_paired = FALSE
+      , stat_label = c("p.format", "p.signif") )
+
+######################
+# Make plot: plotly
+######################
+# FIXME: This labels are not working as I want!
+# lf_bp_plotly(lf_df = lf_deepddg
+#       , p_title = "DeepDDG"
+#       , colour_categ = "deepddg_outcome"
+#       , x_grp = "mutation_info"
+#       , y_var = "param_value"
+#       , facet_var = "param_type"
+#       , n_facet_row = 1
+#       , y_scales = "free_y"
+#       , colour_bp_strip = "khaki2"
+#       , dot_size = 3
+#       , dot_transparency = 0.3
+#       , violin_quantiles = c(0.25, 0.5, 0.75)
+#  )

scripts/functions/test_lf_bp_with_stats.R

@@ -1,28 +0,0 @@
-setwd("~/git/LSHTM_analysis/scripts/plotting/")
-
-source("../functions/lf_bp_with_stats.R")
-
-######################
-# call function
-######################
-# Note: Data
-# run other_plots_data.R
-# to get the long format data to test this function 
-
-lf_bp_with_stats(lf_df = lf_dynamut2
-                       , x_grp = "mutation_info"
-                       , y_var = "param_value"
-                       , facet_var = "param_type"
-                       , n_facet_row = 1
-                       , y_scales = "free_y"
-                       , p_title = "Dynamut2"
-                       , colour_categ = "ddg_dynamut2_outcome"
-                       , stat_grp_comp = c("DM", "OM")
-                       , stat_method = "wilcox.test"
-                       , my_paired = FALSE
-                       #, stat_label = "p.format")
-                       , stat_label = c("p.format", "p.signif")
-                       , my_ats = 22 # axis text size
-                       , my_als = 20 # axis label size
-                       , my_fls = 20 # facet label size
-                       , my_pts = 22 )# plot title size 

scripts/plotting/Header_TT.R

@@ -3,12 +3,6 @@
 #########################################################
 #lib_loc = "/usr/local/lib/R/site-library")
 
-#if (!require("gplots")) {
-#  install.packages("gplots", dependencies = TRUE)
-#  library(gplots)
-#}
-require(extrafont)
-
 require("getopt", quietly = TRUE) # cmd parse arguments
 
 if (!require("tidyverse")) {
@@ -16,9 +10,23 @@ if (!require("tidyverse")) {
   library(tidyverse)
 }
 
-if (!require("ggplot2")) {
+# if (!require("ggplot2")) {
-  install.packages("ggplot2", dependencies = TRUE)
+#   install.packages("ggplot2", dependencies = TRUE)
-  library(ggplot2)
+#   library(ggplot2)
+# }
+
+# if (!require ("dplyr")){
+#   install.packages("dplyr")
+#   library(dplyr)
+# }
+
+# Install
+#if(!require(devtools)) install.packages("devtools")
+#devtools::install_github("kassambara/ggcorrplot")
+
+if (!require ("ggbeeswarm")){
+   install.packages("ggbeeswarm")
+   library(ggbeeswarm)
 }
 
 if (!require("plotly")) {
@@ -101,11 +109,6 @@ if (!require ("psych")){
   library(psych)
 }
 
-if (!require ("dplyr")){
-  install.packages("dplyr")
-  library(dplyr)
-}
-
 if (!require ("compare")){
   install.packages("compare")
   library(compare)
@@ -116,31 +119,25 @@ if (!require ("arsenal")){
   library(arsenal)
 }
 
+if(!require(ggseqlogo)){
+  install.packages("ggseqlogo")
+  library(ggseqlogo)
+}
 
-#if (!requireNamespace("BiocManager", quietly = TRUE))
+# for PDB files
-#  install.packages("BiocManager")
-
-#BiocManager::install("Logolas")
-library("Logolas")
-
-#install.packages("ggseqlogo")
-library(ggseqlogo)
-
-
-####TIDYVERSE
-# Install
-#if(!require(devtools)) install.packages("devtools")
-#devtools::install_github("kassambara/ggcorrplot")
-
-library(ggcorrplot)
-
-
-###for PDB files
-#install.packages("bio3d") 
 if(!require(bio3d)){
   install.packages("bio3d")
   library(bio3d)
 }
 
-#install.packages("protr")
 library(protr)
+if(!require(protr)){
+  install.packages("protr")
+  library(protr)
+}
+
+#if (!requireNamespace("BiocManager", quietly = TRUE))
+#  install.packages("BiocManager")
+
+#BiocManager::install("Logolas")
+library("Logolas")

scripts/plotting/get_plotting_dfs.R

@@ -86,8 +86,10 @@ all_plot_dfs = combining_dfs_plotting(my_df_u
                                       , lig_dist_colname = LigDist_colname
                                       , lig_dist_cutoff = LigDist_cutoff)
 
-merged_df2 = all_plot_dfs[[1]]
+merged_df2      = all_plot_dfs[[1]]
-merged_df3 = all_plot_dfs[[2]]
+merged_df3      = all_plot_dfs[[2]]
+merged_df2_comp = all_plot_dfs[[3]]
+merged_df3_comp = all_plot_dfs[[4]]
 #======================================================================
 # read other files
 infilename_dynamut = paste0("~/git/Data/", drug, "/output/dynamut_results/", gene
@@ -98,10 +100,15 @@ infilename_dynamut2  = paste0("~/git/Data/", drug, "/output/dynamut_results/dyna
 
 infilename_mcsm_na = paste0("~/git/Data/", drug, "/output/mcsm_na_results/", gene
                             , "_complex_mcsm_na_norm.csv")
-
+                            
+infilename_mcsm_f_snps <- paste0("~/git/Data/", drug, "/output/", gene
+                      , "_mcsm_formatted_snps.csv")
+                      
 dynamut_df   = read.csv(infilename_dynamut)
 dynamut2_df  = read.csv(infilename_dynamut2)
 mcsm_na_df   = read.csv(infilename_mcsm_na)
+mcsm_f_snps  = read.csv(infilename_mcsm_f_snps, header = F)
+names(mcsm_f_snps) = "mutationinformation"
 
 ####################################################################
 #                        Data for subcols barplot (~heatmpa)
@@ -430,11 +437,17 @@ if (nrow(corr_ps_df3) == nrow(merged_df3) && nrow(merged_df3_comp) == check1) {
       , "\nGot: ", check1)
 }
 
+
+rm(foo)
+####################################################################
+#                        Data for DM OM Plots: Long format dfs
+####################################################################
+source("other_plots_data.R")
+
 ########################################################################
 #                           End of script
 ########################################################################
-rm(foo)
 
-cat("\n===================================================\n"
+cat("\n######################################################\n"
       , "\nSuccessful: get_plotting_dfs.R worked!"
-      , "\n====================================================")
+      , "\n###################################################\n")

scripts/plotting/other_plots_data.R

@@ -3,10 +3,9 @@
 # TASK: producing boxplots for dr and other muts
 
 #########################################################
-#=======================================================================
 # working dir and loading libraries
 # getwd()
-setwd("~/git/LSHTM_analysis/scripts/plotting")
+# setwd("~/git/LSHTM_analysis/scripts/plotting")
 # getwd()
 
 # make cmd
@@ -14,21 +13,21 @@ setwd("~/git/LSHTM_analysis/scripts/plotting")
 # drug = "streptomycin"
 # gene = "gid"
 
-source("get_plotting_dfs.R")
+# source("get_plotting_dfs.R")
 #=======================================================================
 # MOVE TO COMBINE or singular file for deepddg
+# 
+# cols_to_select = c("mutation", "mutationinformation"
+#                    , "wild_type", "position", "mutant_type"
+#                    , "mutation_info")
+# 
+# merged_df3_short = merged_df3[, cols_to_select]
 
-cols_to_select = c("mutation", "mutationinformation"
+# infilename_mcsm_f_snps <- paste0("~/git/Data/", drug, "/output/", gene
-                   , "wild_type", "position", "mutant_type"
+#                       , "_mcsm_formatted_snps.csv")
-                   , "mutation_info")
+# 
-
+# mcsm_f_snps<- read.csv(infilename_mcsm_f_snps, header = F)
-merged_df3_short = merged_df3[, cols_to_select]
+# names(mcsm_f_snps) <- "mutationinformation"
-
-infilename_mcsm_f_snps <- paste0("~/git/Data/", drug, "/output/", gene
-                      , "_mcsm_formatted_snps.csv")
-
-mcsm_f_snps<- read.csv(infilename_mcsm_f_snps, header = F)
-names(mcsm_f_snps) <- "mutationinformation"
 
 # write merged_df3 to generate structural figure on chimera
 #write.csv(merged_df3_short, "merged_df3_short.csv")
@@ -52,11 +51,11 @@ my_min = min(merged_df3$deepddg_scaled); my_min
 my_max = max(merged_df3$deepddg_scaled); my_max
 
 if (my_min == -1 && my_max == 1){
-  cat("PASS: DeepDDG successfully scaled b/w -1 and 1"
+  cat("\nPASS: DeepDDG successfully scaled b/w -1 and 1"
       #, "\nProceeding with assigning deep outcome category")
       , "\n")
 }else{
-  cat("FAIL: could not scale DeepDDG ddg values"
+  cat("\nFAIL: could not scale DeepDDG ddg values"
       , "Aborting!")
 }
 
@@ -100,7 +99,7 @@ if (merging_cols == "mutationinformation") {
   cols_check <- c(c1, c2, c3, c4)
   expected_cols = n_comb_cols - ( length(cols_check) - 1)
   if (all(cols_check)){
-    cat("\nStage 2:Proceeding with merging dfs:\n")
+    cat("\nStage 2: Proceeding with merging dfs:\n")
     comb_df <- Reduce(inner_join, list(cols_mcsm_df
                                        , cols_mcsm_na_df
                                        , dynamut_df
@@ -115,12 +114,13 @@ if (merging_cols == "mutationinformation") {
     
     }
 }
-names(comb_df_s)
+#names(comb_df_s)
+cat("\n!!!IT GOT TO HERE!!!!")
 #=======================================================================
 fact_cols = colnames(comb_df_s)[grepl( "_outcome|_info", colnames(comb_df_s) )]
 fact_cols
 lapply(comb_df_s[, fact_cols], class)
-comb_df_s[,fact_cols] <- lapply(comb_df_s[,cols],as.factor)
+comb_df_s[, fact_cols] <- lapply(comb_df_s[, fact_cols], as.factor)
 
 if (any(lapply(comb_df_s[, fact_cols], class) == "character")){
   cat("\nChanging cols to factor")
@@ -512,7 +512,6 @@ rm(all_plot_dfs
    , my_data_snp
    , my_df
    , my_df_u
-   , ols_mcsm_df
    , other_muts
    , pd_df
    , subcols_df_ps