#======================================================================= # Task: To generate a logo plot or bar plot but coloured # aa properties. # step1: read mcsm file and OR file # step2: plot wild type positions # step3: plot mutants per position coloured by aa properties # step4: make the size of the letters/bars prop to OR if you can! # useful links # https://stackoverflow.com/questions/5438474/plotting-a-sequence-logo-using-ggplot2 # https://omarwagih.github.io/ggseqlogo/ # https://kkdey.github.io/Logolas-pages/workflow.html # A new sequence logo plot to highlight enrichment and depletion. # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6288878/ #very good: http://www.cbs.dtu.dk/biotools/Seq2Logo-2.0/ #======================================================================= #%% specify curr dir getwd() setwd('~/git/LSHTM_analysis/plotting_test/') getwd() #======================================================================= #%% load packages # header file header_dir = '~/git/LSHTM_analysis/' source(paste0(header_dir, '/', 'my_header.R')) #======================================================================= #%% variable assignment: input and output paths & filenames drug = 'pyrazinamide' gene = 'pncA' gene_match = paste0(gene,'_p.') cat(gene_match) #=========== # data dir #=========== datadir = paste0('~/git/Data') #=========== # input #=========== # source R script 'combining_two_df.R' #indir = paste0(datadir, '/', drug, '/', 'output') # reading files indir = '../meta_data_analysis' # sourcing R script in_filename = 'combining_df_ps.R' infile = paste0(indir, '/', in_filename) cat(paste0('Input is a R script: ', '\'', infile, '\'') , '\n========================================================') #=========== # output #=========== # 1) lineage dist of all muts outdir = paste0('~/git/Data', '/', drug, '/', 'output/plots') #same as indir #cat('Output dir: ', outdir, '\n') #file_type = '.svg' #out_filename1 = paste0(tolower(gene), '_lineage_dist_ps', file_type) #outfile1 = paste0(outdir, '/', out_filename1) #cat(paste0('Output plot1 :', outfile1) # , '\n========================================================') #%% end of variable assignment for input and output files #======================================================================= ##%% read input file cat('Reading input file(sourcing R script):', in_filename) source(infile) #========================== # 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 # since this is one-many relationship # i.e the same SNP can belong to multiple lineages # using the _comp dataset means # we lose some muts and at this level, we should use # as much info as available, hence use df with NA # This will the first plotting df # Then subset this to extract dr muts only (second plottig df) ########################### #%%%%%%%%%%%%%%%%%%%%%%%%% # uncomment as necessary # REASSIGNMENT #my_data = merged_df2 #my_data = merged_df2_comp #my_data = merged_df3 my_data = merged_df3_comp #%%%%%%%%%%%%%%%%%%%%%%%%%% # delete variables not required rm(merged_df2, merged_df2_comp, merged_df3, merged_df3_comp) # quick checks colnames(my_data) str(my_data) c1 = unique(my_data$Position) nrow(my_data) cat('No. of rows in my_data:', nrow(my_data) , '\nDistinct positions corresponding to snps:', length(c1) , '\n===========================================================') #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # FIXME: Think and decide what you want to remove # mut_pos_occurence < 1 or sample_pos_occurrence <1 # get freq count of positions so you can subset freq<1 require(data.table) #setDT(my_data)[, mut_pos_occurrence := .N, by = .(Position)] #265, 14 #extract freq_pos>1 #my_data_snp = my_data[my_data$occurrence!=1,] #u = unique(my_data_snp$Position) #73 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # REASSIGNMENT to prevent changing code my_data_snp = my_data #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #======================================================================= #%% logo plots from dataframe ############# # PLOTS: ggseqlogo with custom height # https://omarwagih.github.io/ggseqlogo/ ############# #require(ggplot2) #require(tidyverse) library(ggseqlogo) foo = my_data_snp[, c("Position", "Mutant_type","ratioDUET", "OR" , "mut_prop_polarity", "mut_prop_water") ] # log10OR # FIXME: at the source script (when calculating AFandOR) my_data_snp$log10or = log10(my_data_snp$OR) bar = my_data_snp[, c('Position', 'Mutant_type', 'OR', 'logor', 'log10or')] bar_or = my_data_snp[, c('Position', 'Mutant_type', 'OR')] wide_df_or <- bar_or %>% spread(Position, OR, fill = 0) wide_df_or = as.matrix(wide_df_or) rownames(wide_df_or) = wide_df_or[,1] wide_df_or = wide_df_or[,-1] # custom height (OR) logo plot: yayy works ggseqlogo(wide_df_or, method='custom', seq_type='aa') + ylab('my custom height') + theme(legend.position = "bottom" , axis.text.x = element_text(size = 11 , angle = 90 , hjust = 1 , vjust = 0.4) , axis.text.y = element_text(size = 15 , angle = 0 , hjust = 1 , vjust = 0))+ labs(title = "AA logo plot" , x = "Wild-type Position" , y = "OR") #%% end of logo plot with OR as height #======================================================================= # extracting data with log10OR bar_logor = my_data_snp[, c('Position', 'Mutant_type', 'log10or')] wide_df_logor <- bar_logor %>% spread(Position, log10or, fill = 0) wide_df_logor = as.matrix(wide_df_logor) rownames(wide_df_logor) = wide_df_logor[,1] wide_df_logor = wide_df_logor[,-1] # custom height (log10OR) logo plot: yayy works ggseqlogo(wide_df_logor, method='custom', seq_type='aa') + ylab('my custom height') + theme(legend.position = "bottom" , axis.text.x = element_text(size = 11 , angle = 90 , hjust = 1 , vjust = 0.4) , axis.text.y = element_text(size = 15 , angle = 0 , hjust = 1 , vjust = 0))+ labs(title = "AA logo plot" , x = "Wild-type Position" , y = "Log10(OR)") #======================================================================= #%% logo plot from sequence ################# # Plot: LOGOLAS (ED plots) # link: https://github.com/kkdey/Logolas # on all pncA samples: output of mutate.py ################ library(Logolas) 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') #%% end of script #=======================================================================