moved plotting_func to functions and replaced 3 basic_barplots scripts with 1

scripts/plotting/functions/plotting_data.R

@@ -0,0 +1,126 @@
+#!/usr/bin/env Rscript             
+#########################################################
+# TASK: formatting data that will be used for various plots
+#########################################################
+# load libraries and functions
+library(data.table)
+library(dplyr)
+#========================================================
+# plotting_data(): formatting data for plots
+# input args: 
+ ## input csv file
+ ## lig cut off dist, default = 10 Ang
+# output: list of 4 dfs, that need to be decompressed
+  ## my_df
+  ## my_df_u
+  ## my_df_u_lig
+  ## dup_muts
+#========================================================
+plotting_data <- function(infile_params, mcsm_lig_cutoff = 10) {
+my_df       = data.frame()
+my_df_u     = data.frame()
+my_df_u_lig = data.frame()
+dup_muts    = data.frame()
+  
+cat(paste0("\nInput file to prepare for plotting:", infile_params, "\n") )
+
+#===========================
+# Read file: struct params
+#===========================
+my_df = read.csv(infile_params, header = T)
+
+cat("\nInput dimensions:", dim(my_df)) 
+
+#==================================
+# add foldx outcome category
+# and foldx scaled values 
+
+# This will enable to always have these variables available
+# when calling for plots
+#==================================
+
+#------------------------------
+# adding foldx scaled values
+# scale data b/w -1 and 1
+#------------------------------
+n = which(colnames(my_df) == "ddg"); n 
+
+my_min = min(my_df[,n]); my_min 
+my_max = max(my_df[,n]); my_max 
+
+my_df$foldx_scaled = ifelse(my_df[,n] < 0
+                            , my_df[,n]/abs(my_min)
+                            , my_df[,n]/my_max) 
+# sanity check
+my_min = min(my_df$foldx_scaled); my_min 
+my_max = max(my_df$foldx_scaled); my_max
+
+if (my_min == -1 && my_max == 1){
+  cat("\nPASS: foldx ddg successfully scaled b/w -1 and 1"
+      , "\nProceeding with assigning foldx outcome category")
+}else{
+  cat("\nFAIL: could not scale foldx ddg values"
+      , "Aborting!\n")
+}
+
+#------------------------------
+# adding foldx outcome category
+# ddg<0 = "Stabilising" (-ve)
+#------------------------------
+c1 = table(my_df$ddg < 0)
+my_df$foldx_outcome = ifelse(my_df$ddg < 0, "Stabilising", "Destabilising")
+c2 = table(my_df$ddg < 0)
+
+if ( all(c1 == c2) ){
+  cat("\nPASS: foldx outcome successfully created")
+}else{
+  cat("\nFAIL: foldx outcome could not be created. Aborting!\n")
+  exit()
+}
+
+#==================================
+# extract unique mutation entries
+#==================================
+
+# check for duplicate mutations
+if ( length(unique(my_df$mutationinformation)) != length(my_df$mutationinformation)){
+  cat(paste0("\nCAUTION:", " Duplicate mutations identified"
+             , "\nExtracting these...\n"))
+  #cat(my_df[duplicated(my_df$mutationinformation),])
+  dup_muts = my_df[duplicated(my_df$mutationinformation),]
+  dup_muts_nu = length(unique(dup_muts$mutationinformation))
+  cat(paste0("\nDim of duplicate mutation df:", nrow(dup_muts)
+             , "\nNo. of unique duplicate mutations:", dup_muts_nu
+             , "\n\nExtracting df with unique mutations only\n"))
+  my_df_u = my_df[!duplicated(my_df$mutationinformation),]
+}else{
+  cat(paste0("\nNo duplicate mutations detected\n"))
+  my_df_u = my_df
+}
+
+upos = unique(my_df_u$position)
+cat("\nDim of clean df:"); cat(dim(my_df_u), "\n")
+cat("\nNo. of unique mutational positions:"); cat(length(upos), "\n")
+  
+#===============================================
+# extract mutations <10 Angstroms and symbol
+#===============================================
+table(my_df_u$ligand_distance<mcsm_lig_cutoff)
+
+my_df_u_lig = my_df_u[my_df_u$ligand_distance <mcsm_lig_cutoff,]
+
+cat(paste0("There are ", nrow(my_df_u_lig), " sites lying within 10\u212b of the ligand\n"))
+
+# return list of DFs
+
+#return(list(my_df, my_df_u, my_df_u_lig, dup_muts))
+#df_names = c("my_df", "my_df_u", "my_df_u_lig", "dup_muts")
+all_df = list(my_df, my_df_u, my_df_u_lig, dup_muts)
+#all_df = Map(setNames, all_df, df_names)
+
+return(all_df)
+}
+
+########################################################################
+#               end of data extraction and cleaning for plots          #
+########################################################################

scripts/plotting/functions/plotting_globals.R

@@ -0,0 +1,61 @@
+#!/usr/bin/env Rscript     
+
+#########################################################
+# TASK: importing global variable for plotting
+# import_dirs()
+# other global variables
+#########################################################
+
+# import_dirs(): similar to mkdirs from python script in repo.
+# input args: 'drug' and 'gene'
+# output: dir names for input and output files
+
+import_dirs <- function(drug, gene) {
+  gene_match = paste0(gene,"_p.")
+  cat(gene_match)
+  
+  #============================
+  # directories and variables
+  #============================
+  
+  datadir        <<- paste0("~/git/Data/")
+  indir          <<- paste0(datadir, drug, "/input")
+  outdir         <<- paste0("~/git/Data/", drug, "/output")
+  plotdir        <<- paste0("~/git/Data/", drug, "/output/plots")
+  
+  dr_muts_col    <<- paste0('dr_mutations_', drug)
+  other_muts_col <<- paste0('other_mutations_', drug)
+  resistance_col <<- "drtype"
+  
+}
+
+# other globals
+#===============================
+# mcsm ligand distance cut off
+#===============================
+#mcsm_lig_cutoff <<- 10
+
+#==================
+# Angstroms symbol
+#==================
+
+angstroms_symbol <<- "\u212b"
+#cat(paste0("There are ", nrow(my_df_u_lig), " sites lying within 10", angstroms_symbol, " of the ligand\n"))
+
+#===============
+# Delta symbol
+#===============
+
+delta_symbol <<- "\u0394"; delta_symbol
+
+#==========
+# Colours
+#==========
+
+mcsm_red2 <<-  "#ae301e" # most negative
+mcsm_red1 <<-  "#f8766d"
+
+mcsm_mid <<- "white" # middle
+
+mcsm_blue1 <<- "#00bfc4"
+mcsm_blue2 <<- "#007d85" # most positive

scripts/plotting/functions/position_count_bp.R

@@ -26,7 +26,7 @@ site_snp_count_bp <- function (plotdf
                                , df_colname = "position"
                                #, bp_plot_title = ""
                                #, leg_title = "Legend title"
-                               #, leg_text_size = 20
+                               , leg_text_size = 20
                                , axis_text_size = 25
                                , axis_label_size = 22
                                , xaxis_title = "Number of nsSNPs"
@@ -111,7 +111,7 @@ site_snp_count_bp <- function (plotdf
           #, legend.position = c(0.73,0.8)
           #, legend.text = element_text(size = leg_text_size)
           #, legend.title = element_text(size =  axis_label_size)
-          , plot.title = element_text(size =  axis_label_size)) + 
+          , plot.title = element_text(size =  leg_text_size)) + 
     
     labs(title = bp_plot_title
       , x = xaxis_title

scripts/plotting/functions/test_functions.R

@@ -1,9 +1,10 @@
 setwd("~/git/LSHTM_analysis/scripts/plotting/functions")
 getwd()
-
-# =================
-# Test function
-# ==================
+#############################################################
+#===========================================
+# load functions, data, dirs, hardocded vars
+# that will be used in testing the functions
+#===========================================
 source("../plotting_data.R")
 infile = "/home/tanu/git/Data/streptomycin/output/gid_comb_stab_struc_params.csv"
 pd_df = plotting_data(infile)
@@ -12,30 +13,91 @@ my_df_u     = pd_df[[2]]
 my_df_u_lig = pd_df[[3]]
 dup_muts    = pd_df[[4]]
 
+source("../plotting_globals.R")
+drug = "streptomycin"
+gene = "gid"
+
+import_dirs(drug, gene)
+
+#=====================
+# functions to test 
+#=====================
+source("stability_count_bp.R")
+source("position_count_bp.R")
+
+##################################################################
 # ------------------------------
 # barplot for mscm stability
 # ------------------------------
+basic_bp_duet =  paste0(tolower(gene), "_basic_barplot_PS.svg")
+plot_basic_bp_duet  =  paste0(plotdir,"/", basic_bp_duet)
+
+svg(plot_basic_bp_duet)
+print(paste0("plot filename:", basic_bp_duet))
+
+# function only
 stability_count_bp(plotdf = my_df_u
                , df_colname = "duet_outcome"
                , leg_title = "DUET outcome")
 
+dev.off()
+
 # ------------------------------
 # barplot for ligand affinity
 # ------------------------------
+basic_bp_ligand =  paste0(tolower(gene), "_basic_barplot_LIG.svg")
+plot_basic_bp_ligand  =  paste0(plotdir, "/", basic_bp_ligand)
+
+svg(plot_basic_bp_ligand)
+print(paste0("plot filename:", basic_bp_ligand))
+
+# function only
 stability_count_bp(plotdf = my_df_u_lig
                , df_colname = "ligand_outcome"
                , leg_title = "Ligand outcome"
-               , bp_plot_title = "Sites < 10 Ang of ligand"
-)
+               , bp_plot_title = "Sites < 10 Ang of ligand")
 
+dev.off()
+# ------------------------------
+# barplot for foldX
+# ------------------------------
+basic_bp_foldx = paste0(tolower(gene), "_basic_barplot_foldx.svg")
+plot_basic_bp_foldx  =  paste0(plotdir,"/", basic_bp_foldx)
+
+svg(plot_basic_bp_foldx)
+print(paste0("plot filename:", plot_basic_bp_foldx))
+
+stability_count_bp(plotdf = my_df_u
+                   , df_colname = "foldx_outcome"
+                   , leg_title = "FoldX outcome")
+dev.off()
+#===============================================================
 # ------------------------------
 # barplot for nssnp site count: all
 # ------------------------------
+pos_count_duet =  paste0(tolower(gene), "_position_count_PS.svg")
+plot_pos_count_duet = paste0(plotdir, "/", pos_count_duet)
+
+svg(plot_pos_count_duet)
+print(paste0("plot filename:", plot_pos_count_duet))
+
+# function only
 site_snp_count_bp(plotdf = my_df_u
                   , df_colname = "position")
 
+dev.off()
 # ------------------------------
 # barplot for nssnp site count: within 10 Ang
 # ------------------------------
+pos_count_ligand =  paste0(tolower(gene), "_position_count_LIG.svg")
+plot_pos_count_ligand = paste0(plotdir, "/", pos_count_ligand)
+
+svg(plot_pos_count_ligand)
+print(paste0("plot filename:", plot_pos_count_ligand))
+
+# function only
 site_snp_count_bp(plotdf = my_df_u_lig
                   , df_colname = "position")
+
+dev.off()
+#===============================================================