moved functions/ in the scripts dir
This commit is contained in:
parent
2dc81e72f4
commit
29fa99b914
6 changed files with 4 additions and 4 deletions
126
scripts/functions/plotting_data.R
Executable file
126
scripts/functions/plotting_data.R
Executable file
|
|
@ -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 #
|
||||
########################################################################
|
||||
61
scripts/functions/plotting_globals.R
Normal file
61
scripts/functions/plotting_globals.R
Normal file
|
|
@ -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
|
||||
125
scripts/functions/position_count_bp.R
Executable file
125
scripts/functions/position_count_bp.R
Executable file
|
|
@ -0,0 +1,125 @@
|
|||
#!/usr/bin/env Rscript
|
||||
|
||||
#########################################################
|
||||
# TASK: function for barplot showing no. of sites with nsSNP
|
||||
# count
|
||||
#########################################################
|
||||
# load libraries and functions
|
||||
library(ggplot2)
|
||||
library(data.table)
|
||||
library(dplyr)
|
||||
|
||||
theme_set(theme_grey())
|
||||
#=================================================================
|
||||
# site_snp_count_bp(): barplots for no. of sites and nsSNP count
|
||||
# input args
|
||||
## df containing data to plot
|
||||
## df column name containing site/position numbers
|
||||
## ...opt args
|
||||
|
||||
# TODO:
|
||||
# legend params, currently it appears as title
|
||||
# visually might be nicer for it to be inside the plot
|
||||
#=================================================================
|
||||
|
||||
site_snp_count_bp <- function (plotdf
|
||||
, df_colname = "position"
|
||||
#, bp_plot_title = ""
|
||||
#, leg_title = "Legend title"
|
||||
, leg_text_size = 20
|
||||
, axis_text_size = 25
|
||||
, axis_label_size = 22
|
||||
, xaxis_title = "Number of nsSNPs"
|
||||
, yaxis_title = "Number of Sites"){
|
||||
|
||||
# dim of plotdf
|
||||
cat(paste0("\noriginal df dimensions:"
|
||||
, "\nNo. of rows:", nrow(plotdf)
|
||||
, "\nNo. of cols:", ncol(plotdf)
|
||||
, "\nNow adding column: frequency of mutational positions"))
|
||||
|
||||
# adding snpcount for each position
|
||||
setDT(plotdf)[, pos_count := .N, by = .(eval(parse(text = df_colname)))]
|
||||
|
||||
cat("\nCumulative nssnp count\n"
|
||||
, table(plotdf$pos_count))
|
||||
|
||||
# calculating total no. of mutations
|
||||
tot_muts = sum(table(plotdf$pos_count))
|
||||
|
||||
# sanity check
|
||||
if(tot_muts == nrow(plotdf)){
|
||||
cat("\nPASS: total number of mutations match"
|
||||
, "\nTotal no. of nsSNPs:", tot_muts)
|
||||
} else{
|
||||
cat("\nWARNING: total no. of muts = ", tot_muts
|
||||
, "\nExpected = ", nrow(plotdf))
|
||||
}
|
||||
|
||||
# revised dimensions
|
||||
cat(paste0("\nrevised df dimensions:"
|
||||
, "\nNo. of rows:", nrow(plotdf)
|
||||
, "\nNo. of cols:", ncol(plotdf)))
|
||||
|
||||
# use group by on pos_count
|
||||
snpsBYpos_df <- plotdf %>%
|
||||
group_by(eval(parse(text = df_colname))) %>%
|
||||
summarize(snpsBYpos = mean(pos_count))
|
||||
|
||||
cat("\nnssnp count\n"
|
||||
, table(snpsBYpos_df$snpsBYpos))
|
||||
|
||||
# calculating total no. of sites associated with nsSNPs
|
||||
tot_sites = sum(table(snpsBYpos_df$snpsBYpos))
|
||||
|
||||
# sanity check
|
||||
if(tot_sites == length(unique(plotdf$position))){
|
||||
cat("\nPASS: total number of mutation sites match"
|
||||
, "\nTotal no. of sites with nsSNPs:", tot_sites)
|
||||
} else{
|
||||
cat("WARNING: total no. of sites = ", tot_sites
|
||||
, "\nExpected = ", length(unique(plotdf$position)))
|
||||
}
|
||||
|
||||
# FIXME: should really be legend title
|
||||
# but atm being using as plot title
|
||||
my_leg_title = paste0("Total nsSNPs:", tot_muts
|
||||
, ", Total no. of nsSNPs sites:", tot_sites)
|
||||
bp_plot_title = my_leg_title
|
||||
|
||||
#-------------
|
||||
# start plot 2
|
||||
#--------------
|
||||
# to make x axis display all positions
|
||||
# not sure if to use with sort or directly
|
||||
my_x = sort(unique(snpsBYpos_df$snpsBYpos))
|
||||
|
||||
g = ggplot(snpsBYpos_df, aes(x = snpsBYpos))
|
||||
OutPlot_pos_count = g + geom_bar(aes (alpha = 0.5)
|
||||
, show.legend = FALSE) +
|
||||
scale_x_continuous(breaks = unique(snpsBYpos_df$snpsBYpos)) +
|
||||
geom_label(stat = "count", aes(label = ..count..)
|
||||
, color = "black"
|
||||
, size = 10) +
|
||||
theme(axis.text.x = element_text(size = axis_text_size
|
||||
, angle = 0)
|
||||
, axis.text.y = element_text(size = axis_text_size
|
||||
, angle = 0
|
||||
, hjust = 1)
|
||||
, axis.title.x = element_text(size = axis_label_size)
|
||||
, axis.title.y = element_text(size = axis_label_size)
|
||||
#, 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 = leg_text_size)) +
|
||||
|
||||
labs(title = bp_plot_title
|
||||
, x = xaxis_title
|
||||
, y = yaxis_title)
|
||||
|
||||
return(OutPlot_pos_count)
|
||||
}
|
||||
|
||||
########################################################################
|
||||
# end of function
|
||||
########################################################################
|
||||
50
scripts/functions/stability_count_bp.R
Normal file
50
scripts/functions/stability_count_bp.R
Normal file
|
|
@ -0,0 +1,50 @@
|
|||
#!/usr/bin/env Rscript
|
||||
|
||||
#########################################################
|
||||
# TASK: function for basic barplot returning stability counts
|
||||
#########################################################
|
||||
# load libraries and functions
|
||||
library(ggplot2)
|
||||
theme_set(theme_grey())
|
||||
#==========================================================
|
||||
# stability_count_bp(): basic barplots for stability counts
|
||||
# input args
|
||||
## df containing data to plot
|
||||
## df column name containing stability outcome
|
||||
## legend title
|
||||
## ...opt args
|
||||
#==========================================================
|
||||
stability_count_bp <- function(plotdf
|
||||
, df_colname
|
||||
, leg_title = "Legend title"
|
||||
, axis_text_size = 25
|
||||
, axis_label_size = 22
|
||||
, leg_text_size = 20
|
||||
, yaxis_title = "Number of nsSNPs"
|
||||
, bp_plot_title = ""){
|
||||
|
||||
OutPlot_count = ggplot(plotdf, aes(x = eval(parse(text = df_colname)))) +
|
||||
geom_bar(aes(fill = eval(parse(text = df_colname))), show.legend = TRUE) +
|
||||
geom_label(stat = "count"
|
||||
, aes(label = ..count..)
|
||||
, color = "black"
|
||||
, show.legend = FALSE
|
||||
, size = 10) +
|
||||
theme(axis.text.x = element_blank()
|
||||
, axis.title.x = element_blank()
|
||||
, axis.title.y = element_text(size = axis_label_size)
|
||||
, axis.text.y = element_text(size = axis_text_size)
|
||||
, 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)) +
|
||||
labs(title = bp_plot_title
|
||||
, y = yaxis_title) +
|
||||
scale_fill_discrete(name = leg_title
|
||||
, labels = c("Destabilising", "Stabilising"))
|
||||
|
||||
return(OutPlot_count)
|
||||
}
|
||||
#############################################################
|
||||
# end of function
|
||||
#############################################################
|
||||
103
scripts/functions/test_functions.R
Normal file
103
scripts/functions/test_functions.R
Normal file
|
|
@ -0,0 +1,103 @@
|
|||
setwd("~/git/LSHTM_analysis/scripts/plotting/functions")
|
||||
getwd()
|
||||
#############################################################
|
||||
#===========================================
|
||||
# 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)
|
||||
my_df = pd_df[[1]]
|
||||
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")
|
||||
|
||||
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()
|
||||
#===============================================================
|
||||
Loading…
Add table
Add a link
Reference in a new issue