updated logo_plot.R with functions

scripts/plotting/logo_plot.R

@@ -2,43 +2,80 @@
 #########################################################
 # TASK: producing logoplot 
 # from data and/or from sequence
-
 #########################################################
-#=======================================================================
 # working dir and loading libraries
 getwd()
 setwd("~/git/LSHTM_analysis/scripts/plotting")
 getwd()
 
 source("Header_TT.R")
-#library(ggplot2)
+source("../functions/plotting_globals.R")
-#library(data.table)
+source("../functions/plotting_data.R")
-#library(dplyr)
+source("../functions/combining_dfs_plotting.R")
+###########################################################
+# command line args
+#********************
+drug = 'streptomycin'
+gene = 'gid'
 
 #===========
 # input
 #===========
-source("combining_dfs_plotting.R")
+#---------------------
+# call: import_dirs()
+#---------------------
+import_dirs(drug, gene)
+
+#---------------------------
+# call: plotting_data()
+#---------------------------
+if (!exists("infile_params") && exists("gene")){
+  #if (!is.character(infile_params) && exists("gene")){
+  #in_filename_params = paste0(tolower(gene), "_all_params.csv") 
+  in_filename_params = paste0(tolower(gene), "_comb_afor.csv") # part combined for gid
+  infile_params = paste0(outdir, "/", in_filename_params)
+  cat("\nInput file for mcsm comb data not specified, assuming filename: ", infile_params, "\n")
+}
+
+# Input 1: read <gene>_comb_afor.csv
+pd_df = plotting_data(infile_params)
+my_df_u       = pd_df[[1]] # this forms one of the input for combining_dfs_plotting()
+
+#--------------------------------
+# call: combining_dfs_plotting()
+#--------------------------------
+if (!exists("infile_metadata") && exists("gene")){
+  #if (!is.character(infile_params) && exists("gene")){{
+  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")
+}
+
+# Input 2: read <gene>_meta data.csv
+cat("\nReading meta data file:", infile_metadata)
+
+gene_metadata <- read.csv(infile_metadata
+                          , stringsAsFactors = F
+                          , header = T)
+
+all_plot_dfs = combining_dfs_plotting(my_df_u
+                                      , gene_metadata
+                                      , lig_dist_colname = 'ligand_distance'
+                                      , lig_dist_cutoff = 10)
+
+merged_df2        = all_plot_dfs[[1]]
+merged_df3        = all_plot_dfs[[2]]
+#merged_df2_comp   = all_plot_dfs[[3]]
+#merged_df3_comp   = all_plot_dfs[[4]]
+#merged_df2_lig    = all_plot_dfs[[5]]
+#merged_df3_lig    = all_plot_dfs[[6]]
 
 #===========
 # output
 #===========
-
 logo_plot = "logo_plot.svg"
 plot_logo_plot = paste0(plotdir,"/", logo_plot)
 
-#==========================
-# This will return:
-
-# df with NA for pyrazinamide:
-# merged_df2
-# merged_df3
-
-# df without NA for pyrazinamide:
-# merged_df2_comp
-# merged_df3_comp
-#===========================
-
 ###########################
 # Data for plots
 # you need merged_df2 or merged_df2_comp
@@ -58,13 +95,9 @@ plot_logo_plot = paste0(plotdir,"/", logo_plot)
 #my_data = merged_df2
 #my_data =  merged_df2_comp
 my_data = merged_df3
-my_data = merged_df3_comp
+#my_data = merged_df3_comp
 #%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-# delete variables not required
-rm(merged_df2, merged_df2_comp)
-#rm(merged_df3, merged_df3_comp)
-
 # quick checks
 colnames(my_data)
 str(my_data)
@@ -80,7 +113,7 @@ cat("No. of rows in my_data:", nrow(my_data)
 # mut_pos_occurence < 1 or sample_pos_occurrence <1
 
 # get freq count of positions so you can subset freq<1
-require(data.table)
+#require(data.table)
 #setDT(my_data)[, mut_pos_occurrence := .N, by = .(position)] #265, 14
 
 #extract freq_pos>1
@@ -100,15 +133,13 @@ my_data_snp = my_data
 # PLOTS: ggseqlogo with custom height
 # https://omarwagih.github.io/ggseqlogo/
 #############
-#require(ggplot2)
+foo = my_data_snp[, c("position"
-#require(tidyverse)
+                      , "mutant_type","duet_scaled", "or_mychisq"
-library(ggseqlogo)
+                      , "mut_prop_polarity", "mut_prop_water")] 
-
-foo = my_data_snp[, c("position", "mutant_type","duet_scaled", "or_mychisq"
-                      , "mut_prop_polarity", "mut_prop_water") ] 
 
 my_data_snp$log10or = log10(my_data_snp$or_mychisq)
-logo_data = my_data_snp[, c("position", "mutant_type", "or_mychisq", "log10or")]
+logo_data = my_data_snp[, c("position"
+                            , "mutant_type", "or_mychisq", "log10or")]
 
 logo_data_or = my_data_snp[, c("position", "mutant_type", "or_mychisq")]
 wide_df_or <- logo_data_or %>% spread(position, or_mychisq, fill = 0.0)
@@ -124,12 +155,11 @@ position_or = as.numeric(colnames(wide_df_or))
 #custom height (OR) logo plot: CORRECT x-axis labelling
 #============================================  
 # custom height (OR) logo plot: yayy works
-
 cat("Logo plot with OR as y axis:", plot_logo_plot)
 svg(plot_logo_plot, width = 30 , height = 6)
 
 logo_or = ggseqlogo(wide_df_or, method="custom", seq_type="aa") + ylab("my custom height") +
-    theme( axis.text.x = element_text(size = 16
+    theme( axis.text.x = element_text(size = 12
                                        , angle = 90
                                        , hjust = 1
                                        , vjust = 0.4)
@@ -174,7 +204,7 @@ colnames(wide_df_logor_m)
 
 position_logor = as.numeric(colnames(wide_df_logor_m))
 
-#  custom height (log10OR) logo plot: yayy works
+# custom height (log10OR) logo plot: yayy works
 ggseqlogo(wide_df_logor_m, method="custom", seq_type="aa") + ylab("my custom height") +
   theme(legend.position = "bottom"
         , axis.text.x = element_text(size = 11
@@ -191,85 +221,6 @@ ggseqlogo(wide_df_logor_m, method="custom", seq_type="aa") + ylab("my custom hei
                    , limits = factor(1:length(position_logor)))+
   ylab("Log (Odds Ratio)") + 
   scale_y_continuous(limits = c(0, 9))
-
 #=======================================================================
-#%% logo plot from sequence
-
-#################
-# Plot: LOGOLAS (ED plots)
-# link: https://github.com/kkdey/Logolas
-# on all pncA samples: output of mutate.py
-################
-library(Logolas)
-
-# data was pnca_msa.txt
-#FIXME: generate this file 
-
-seqs = read.csv("~/git/Data/pyrazinamide/output/pnca_msa.txt"
-                , header = FALSE
-                , stringsAsFactors = FALSE)$V1
-
-
-# my_data: useful!
-logomaker(seqs, type = "EDLogo", color_type = "per_symbol"
-          , return_heights = TRUE)
-
-logomaker(seqs, type = "Logo", color_type = "per_symbol", return_heights = TRUE)
-
 #%% end of script
 #=======================================================================
-#==============
-# online logo:
-#http://www.cbs.dtu.dk/biotools/Seq2Logo/ # good for getting pssm and psfm matrices
-#https://weblogo.berkeley.edu/logo.cgi
-#http://weblogo.threeplusone.com/create.cgi   # weblogo3
-
-#===============
-# PSSM logos= example from logomaker
-#===============
-
-data(pssm)
-logo_pssm(pssm, control = list(round_off = 0))
-
-#=================
-# PSSM: output from http://www.cbs.dtu.dk/biotools/Seq2Logo/
-# of MSA: pnca_msa.txt
-#==================
-foo = read.csv("/home/tanu/git/Data/pyrazinamide/pssm_transpose.csv")
-rownames(foo) = foo[,1]
-df = subset(foo, select = -c(1) )
-colnames(df)
-colnames(df) = seq(1:length(df))
-
-df_m = as.matrix(df)
-logo_pssm(df_m, control = list(round_off = 0))
-
-#=================
-# # PSFM: output from http://www.cbs.dtu.dk/biotools/Seq2Logo/
-# of MSA: pnca_msa.txt
-#=================
-# not designed for PSFM!
-# may want to figure out how to do it!
-logo_data = read.csv("/home/tanu/git/Data/pyrazinamide/psfm_transpose.csv")
-rownames(logo_data) = logo_data[,1]
-df2 = subset(logo_data, select = -c(1) )
-colnames(df2)
-colnames(df2) = seq(1:length(df2))
-
-df2_m = as.matrix(df2)
-logo_pssm(df2_m, control = list(round_off = 0))
-
-
-#=================
-# ggplots:
-#https://stackoverflow.com/questions/5438474/plotting-a-sequence-logo-using-ggplot2
-#=================
-
-library(ggplot2)
-library(gglogo)
-ggplot(data = ggfortify(sequences, "peptide")) +      
-  geom_logo(aes(x=position, y=bits, group=element, 
-                label=element, fill=interaction(Polarity, Water)),
-            alpha = 0.6)  +
-  scale_fill_brewer(palette="Paired") +
-  theme(legend.position = "bottom")