108 lines
3.4 KiB
R
Executable file
108 lines
3.4 KiB
R
Executable file
#!/usr/bin/env Rscript
|
|
#########################################################
|
|
# TASK: formatting data that will be used for various plots
|
|
|
|
# useful links
|
|
#https://stackoverflow.com/questions/38851592/r-append-column-in-a-dataframe-with-frequency-count-based-on-two-columns
|
|
#########################################################
|
|
# working dir and loading libraries
|
|
getwd()
|
|
setwd("~/git/LSHTM_analysis/scripts/plotting")
|
|
getwd()
|
|
|
|
#source("Header_TT.R")
|
|
library(ggplot2)
|
|
library(data.table)
|
|
library(dplyr)
|
|
|
|
require("getopt", quietly = TRUE) #cmd parse arguments
|
|
#========================================================
|
|
# command line args
|
|
#spec = matrix(c(
|
|
# "drug" , "d", 1, "character",
|
|
# "gene" , "g", 1, "character"
|
|
#), byrow = TRUE, ncol = 4)
|
|
|
|
#opt = getopt(spec)
|
|
|
|
#drug = opt$druggene = opt$gene
|
|
|
|
#if(is.null(drug)|is.null(gene)) {
|
|
# stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)")
|
|
#}
|
|
#========================================================
|
|
#%% variable assignment: input and output paths & filenames
|
|
drug = "pyrazinamide"
|
|
gene = "pncA"
|
|
gene_match = paste0(gene,"_p.")
|
|
cat(gene_match)
|
|
|
|
#=============
|
|
# directories
|
|
#=============
|
|
datadir = paste0("~/git/Data")
|
|
indir = paste0(datadir, "/", drug, "/input")
|
|
outdir = paste0("~/git/Data", "/", drug, "/output")
|
|
plotdir = paste0("~/git/Data", "/", drug, "/output/plots")
|
|
#======
|
|
# input
|
|
#======
|
|
#in_filename = "mcsm_complex1_normalised.csv"
|
|
in_filename_params = paste0(tolower(gene), "_all_params.csv")
|
|
infile_params = paste0(outdir, "/", in_filename_params)
|
|
cat(paste0("Input file 1:", infile_params) )
|
|
|
|
#%%===============================================================
|
|
###########################
|
|
# Read file: struct params
|
|
###########################
|
|
cat("Reading struct params including mcsm:", in_filename_params)
|
|
|
|
my_df = read.csv(infile_params, header = T)
|
|
|
|
cat("\nInput dimensions:", dim(my_df))
|
|
|
|
# quick checks
|
|
#colnames(my_df)
|
|
#str(my_df)
|
|
|
|
###########################
|
|
# 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..."))
|
|
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"))
|
|
my_df_u = my_df[!duplicated(my_df$mutationinformation),]
|
|
}else{
|
|
cat(paste0("\nNo duplicate mutations detected"))
|
|
my_df_u = my_df
|
|
}
|
|
|
|
upos = unique(my_df_u$position)
|
|
cat("\nDim of clean df:"); cat(dim(my_df_u))
|
|
cat("\nNo. of unique mutational positions:"); cat(length(upos), "\n")
|
|
|
|
|
|
###########################
|
|
# extract mutations <10Angstroms
|
|
###########################
|
|
table(my_df_u$ligand_distance<10)
|
|
|
|
my_df_u_lig = my_df_u[my_df_u$ligand_distance <10,]
|
|
angstroms_symbol = "\u212b"
|
|
cat(paste0("There are ", nrow(my_df_u_lig), " sites lying within 10", angstroms_symbol, " of the ligand"))
|
|
|
|
########################################################################
|
|
# end of data extraction and cleaning for plots #
|
|
########################################################################
|
|
# clear variables
|
|
rm(opt, spec)
|
|
|
|
|