From 9da808680dc5b19e63ff0fdee6fe18f673bfe94f Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Thu, 10 Sep 2020 15:28:10 +0100 Subject: [PATCH] re-adding deleted combining_dfs_plotting.R --- scripts/combining_dfs_plotting.R | 439 +++++++++++++++++++++++++++++++ scripts/plotting/other_plots.R | 10 +- 2 files changed, 443 insertions(+), 6 deletions(-) create mode 100644 scripts/combining_dfs_plotting.R diff --git a/scripts/combining_dfs_plotting.R b/scripts/combining_dfs_plotting.R new file mode 100644 index 0000000..0a4c303 --- /dev/null +++ b/scripts/combining_dfs_plotting.R @@ -0,0 +1,439 @@ +#!/usr/bin/env Rscript +######################################################### +# TASK: To combine struct params and meta data for plotting +# Input csv files: +# 1) _all_params.csv +# 2) _meta_data.csv + +# Output: +# 1) muts with opposite effects on stability +# 2) large combined df including NAs for AF, OR,etc +# Dim: same no. of rows as gene associated meta_data_with_AFandOR +# 3) small combined df including NAs for AF, OR, etc. +# Dim: same as mcsm data +# 4) large combined df excluding NAs +# Dim: dim(#1) - na_count_df2 +# 5) small combined df excluding NAs +# Dim: dim(#2) - na_count_df3 +# This script is sourced from other .R scripts for plotting +######################################################### +#======================================================================= +# working dir and loading libraries +getwd() +setwd("~/git/LSHTM_analysis/scripts/plotting/") +getwd() + +source("Header_TT.R") +#require(data.table) +#require(arsenal) +#require(compare) +#library(tidyverse) +source("plotting_data.R") + +# should return the following dfs, directories and variables +# my_df +# my_df_u +# my_df_u_lig +# dup_muts + +cat(paste0("Directories imported:" + , "\ndatadir:", datadir + , "\nindir:", indir + , "\noutdir:", outdir + , "\nplotdir:", plotdir)) + +cat(paste0("Variables imported:" + , "\ndrug:", drug + , "\ngene:", gene + , "\ngene_match:", gene_match + , "\nLength of upos:", length(upos) + , "\nAngstrom symbol:", angstroms_symbol)) + +# clear excess variable +rm(my_df, upos, dup_muts) +#======================================================== +#=========== +# input +#=========== +#in_file1: output of plotting_data.R +# my_df_u + +# quick checks +head(my_df_u[, c("mutation")]) + +cols_to_extract = c("mutationinformation", "mutation", "or_mychisq", "or_kin", "af", "af_kin") +foo = my_df_u[, colnames(my_df_u)%in% cols_to_extract] + + +table(which(is.na(my_df_u$af_kin)) == which(is.na(my_df_u$af))) + +baz = read.csv(file.choose()) + +baz = cbind(my_df_u$mutation, my_df_u$or_mychisq, bar$mutation, bar$or_mychisq) +baz = as.data.frame(baz) +colnames(baz) = c("my_df_u_muts", "my_df_u_or", "real_muts", "real_or") +sum(is.na(baz$my_df_u_or)) == sum(is.na(my_df_u$or_mychisq)) + +cat("\nNo. of with NA in or_mychisq:", sum(is.na(my_df_u$or_mychisq)) + ,"\nNo. of NA in or_kin:" , sum(is.na(my_df_u$or_kin))) + +# infile 2: gene associated meta data +#in_filename_gene_metadata = paste0(tolower(gene), "_meta_data_with_AFandOR.csv") +in_filename_gene_metadata = paste0(tolower(gene), "_metadata.csv") +infile_gene_metadata = paste0(outdir, "/", in_filename_gene_metadata) +cat(paste0("Input infile 2:", infile_gene_metadata)) + +#=========== +# output +#=========== +# other variables that you can write +# primarily called by other scripts for plotting + +# PS combined: +# 1) merged_df2 +# 2) merged_df2_comp +# 3) merged_df3 +# 4) merged_df3_comp + +# LIG combined: +# 5) merged_df2_lig +# 6) merged_df2_comp_lig +# 7) merged_df3_lig +# 8) merged_df3_comp_lig + +#%%=============================================================== + +########################### +# 2: Read file: _meta data.csv +########################### +cat("Reading meta data file:", infile_gene_metadata) + +gene_metadata <- read.csv(infile_gene_metadata + , stringsAsFactors = F + , header = T) +cat("Dim:", dim(gene_metadata)) + + +# counting NAs in AF, OR cols +# or_mychisq +if (identical(sum(is.na(my_df_u$or_mychisq)) + , sum(is.na(my_df_u$pval_fisher)) + , sum(is.na(my_df_u$af)))){ + cat("\nPASS: NA count match for OR, pvalue and AF\n") + na_count = sum(is.na(my_df_u$af)) + cat("\nNo. of NAs: ", sum(is.na(my_df_u$or_mychisq))) +} else{ + cat("\nFAIL: NA count mismatch" + , "\nNA in OR: ", sum(is.na(my_df_u$or_mychisq)) + , "\nNA in pvalue: ", sum(is.na(my_df_u$pval_fisher)) + , "\nNA in AF:", sum(is.na(my_df_u$af))) +} + +# or kin +if (identical(sum(is.na(my_df_u$or_kin)) + , sum(is.na(my_df_u$pwald_kin)) + , sum(is.na(my_df_u$af_kin)))){ + cat("\nPASS: NA count match for OR, pvalue and AF\n from Kinship matrix calculations") + na_count = sum(is.na(my_df_u$af_kin)) + cat("\nNo. of NAs: ", sum(is.na(my_df_u$or_kin))) +} else{ + cat("\nFAIL: NA count mismatch" + , "\nNA in OR: ", sum(is.na(my_df_u$or_kin)) + , "\nNA in pvalue: ", sum(is.na(my_df_u$pwald_kin)) + , "\nNA in AF:", sum(is.na(my_df_u$af_kin))) +} + +str(gene_metadata) + +# change category of ambiguos mutations +table(gene_metadata$mutation_info) + +cols_to_extract2 = c("mutationinformation", "mutation", "mutation_info") +foo2 = gene_metadata[, colnames(gene_metadata)%in% cols_to_extract2] + +dr_muts = foo2[foo2$mutation_info == dr_muts_col,] +other_muts = foo2[foo2$mutation_info == other_muts_col,] + +common_muts = dr_muts[dr_muts$mutation%in%other_muts$mutation,] +#write.csv(common_muts, 'common_muts.csv') +rm(common_muts) + +# FIXME read properly +# "ambiguous_mut_names.csv" +#"pnca_p.gly108arg", "pnca_p.gly132ala", "pnca_p.val180phe" +ambiguous_muts = read.csv(file.choose()) +ambiguous_muts_names = ambiguous_muts$mutation + +common_muts_all = gene_metadata[gene_metadata$mutation%in%ambiguous_muts_names,] + +if (gene_metadata$mutation_info[gene_metadata$mutation%in%ambiguous_muts_names] == other_muts_col){ + print('change me') +} + +# make a copy +gene_metadata2 = gene_metadata +table(gene_metadata$mutation_info) +count_check = as.data.frame(cbind(table(gene_metadata$mutationinformation, gene_metadata$mutation_info))) +#count_check$checks = ifelse(count_check$dr_mutations_pyrazinamide&&count_check$other_mutations_pyrazinamide>0, "ambi", "pass") +table(count_check$checks) + + +poo = c("V180F", "G132A", "D49G") +poo2 = count_check[rownames(count_check)%in%poo,] +poo2[[dr_muts_col]]&& poo2[[other_muts_col]]>0 +poo2$checks = ifelse(all(poo2$checkspoo2[[dr_muts_col]]&& poo2[[other_muts_col]])>0, "ambi", "pass") + +# remove common_muts_all +ids = gene_metadata$mutation%in%common_muts_all$mutation; table(ids) +gene_metadata_unambiguous = gene_metadata2[!ids,] + +# sanity checks: should be true +table(gene_metadata_unambiguous$mutation%in%common_muts_all$mutation)[[1]] == nrow(gene_metadata_unambiguous) +nrow(gene_metadata_unambiguous) + nrow(common_muts_all) == nrow(gene_metadata) + +# correct common muts +table(common_muts_all$mutation_info) +common_muts_all$mutation_info = as.factor(common_muts_all$mutation_info) + +# change the other_muts to dr_muts +common_muts_all$mutation_info[common_muts_all$mutation_info==other_muts_col] <- dr_muts_col + +table(common_muts_all$mutation_info) +common_muts_all$mutation_info = factor(common_muts_all$mutation_info) +table(common_muts_all$mutation_info) + +# add it back to +gene_metadata2 = rbind(gene_metadata_unambiguous, common_muts_all) +nrow(gene_metadata2) == nrow(gene_metadata) + +################################################################### +# combining: PS +################################################################### +# sort by position (same as my_df) +head(gene_metadata$position) +gene_metadata = gene_metadata[order(gene_metadata$position),] +head(gene_metadata$position) + +#========================= +# Merge 1: merged_df2 +# dfs with NAs in ORs +#========================= +head(my_df_u$mutationinformation) +head(gene_metadata$mutationinformation) + +# Find common columns b/w two df +merging_cols = intersect(colnames(my_df_u), colnames(gene_metadata)) + +cat(paste0("Merging dfs with NAs: big df (1-many relationship b/w id & mut)" + , "\nNo. of merging cols:", length(merging_cols) + , "\nMerging columns identified:")) +print(merging_cols) + +# important checks! +table(nchar(my_df_u$mutationinformation)) +table(nchar(my_df_u$wild_type)) +table(nchar(my_df_u$mutant_type)) +table(nchar(my_df_u$position)) + +# all.y because x might contain non-structural positions! +merged_df2 = merge(x = gene_metadata + , y = my_df_u + , by = merging_cols + , all.y = T) + +cat("Dim of merged_df2: ", dim(merged_df2)) +head(merged_df2$position) + +# sanity check +cat("Checking nrows in merged_df2") +if(nrow(gene_metadata) == nrow(merged_df2)){ + cat("PASS: nrow(merged_df2) = nrow (gene associated gene_metadata)" + ,"\nExpected no. of rows: ",nrow(gene_metadata) + ,"\nGot no. of rows: ", nrow(merged_df2)) +} else{ + cat("FAIL: nrow(merged_df2)!= nrow(gene associated gene_metadata)" + , "\nExpected no. of rows after merge: ", nrow(gene_metadata) + , "\nGot no. of rows: ", nrow(merged_df2) + , "\nFinding discrepancy") + merged_muts_u = unique(merged_df2$mutationinformation) + meta_muts_u = unique(gene_metadata$mutationinformation) + # find the index where it differs + unique(meta_muts_u[! meta_muts_u %in% merged_muts_u]) + quit() +} + +#========================= +# Merge 2: merged_df3 +# dfs with NAs in ORs +# +# Cannot trust lineage, country from this df as the same mutation +# can have many different lineages +# but this should be good for the numerical corr plots +#========================= +# remove duplicated mutations +cat("Merging dfs without NAs: small df (removing muts with no AF|OR associated)" + ,"\nCannot trust lineage info from this" + ,"\nlinking col: mutationinforamtion" + ,"\nfilename: merged_df3") + +merged_df3 = merged_df2[!duplicated(merged_df2$mutationinformation),] +head(merged_df3$position); tail(merged_df3$position) # should be sorted + +# sanity check +cat("Checking nrows in merged_df3") +if(nrow(my_df_u) == nrow(merged_df3)){ + cat("PASS: No. of rows match with my_df" + ,"\nExpected no. of rows: ", nrow(my_df_u) + ,"\nGot no. of rows: ", nrow(merged_df3)) +} else { + cat("FAIL: No. of rows mismatch" + , "\nNo. of rows my_df: ", nrow(my_df_u) + , "\nNo. of rows merged_df3: ", nrow(merged_df3)) + quit() +} + +# counting NAs in AF, OR cols in merged_df3 +# this is because mcsm has no AF, OR cols, +# so you cannot count NAs +if (identical(sum(is.na(merged_df3$or_kin)) + , sum(is.na(merged_df3$pwald_kin)) + , sum(is.na(merged_df3$af_kin)))){ + cat("PASS: NA count match for OR, pvalue and AF\n") + na_count_df3 = sum(is.na(merged_df3$af_kin)) + cat("No. of NAs: ", sum(is.na(merged_df3$or_kin))) +} else{ + cat("FAIL: NA count mismatch" + , "\nNA in OR: ", sum(is.na(merged_df3$or_kin)) + , "\nNA in pvalue: ", sum(is.na(merged_df3$pwald_kin)) + , "\nNA in AF:", sum(is.na(merged_df3$af_kin))) +} + +# check if the same or and afs are missing for +if ( identical( which(is.na(merged_df2$or_mychisq)), which(is.na(merged_df2$or_kin))) + && identical( which(is.na(merged_df2$af)), which(is.na(merged_df2$af_kin))) + && identical( which(is.na(merged_df2$pval_fisher)), which(is.na(merged_df2$pwald_kin))) ){ + cat("PASS: Indices match for mychisq and kin ors missing values") +} else{ + cat("Index mismatch: mychisq and kin ors missing indices match") + quit() +} + +#========================= +# Merge3: merged_df2_comp +# same as merge 1 but excluding NAs from ORs, etc. +#========================= +cat("Merging dfs without any NAs: big df (1-many relationship b/w id & mut)" + ,"\nlinking col: Mutationinforamtion" + ,"\nfilename: merged_df2_comp") + +if ( identical( which(is.na(merged_df2$af)), which(is.na(merged_df2$af_kin))) ){ + print("mychisq and kin ors missing indices match. Procedding with omitting NAs") + na_count_df2 = sum(is.na(merged_df2$af)) + merged_df2_comp = merged_df2[!is.na(merged_df2$af),] + # sanity check: no +-1 gymnastics + cat("Checking nrows in merged_df2_comp") + if(nrow(merged_df2_comp) == (nrow(merged_df2) - na_count_df2)){ + cat("\nPASS: No. of rows match" + ,"\nDim of merged_df2_comp: " + ,"\nExpected no. of rows: ", nrow(merged_df2) - na_count_df2 + , "\nNo. of rows: ", nrow(merged_df2_comp) + , "\nNo. of cols: ", ncol(merged_df2_comp)) + }else{ + cat("FAIL: No. of rows mismatch" + ,"\nExpected no. of rows: ", nrow(merged_df2) - na_count_df2 + ,"\nGot no. of rows: ", nrow(merged_df2_comp)) + } +}else{ + print("Index mismatch for mychisq and kin ors. Aborting NA ommission") +} + +#========================= +# Merge4: merged_df3_comp +# same as merge 2 but excluding NAs from ORs, etc or +# remove duplicate mutation information +#========================= + +if ( identical( which(is.na(merged_df3$af)), which(is.na(merged_df3$af_kin))) ){ + print("mychisq and kin ors missing indices match. Procedding with omitting NAs") + na_count_df3 = sum(is.na(merged_df3$af)) + #merged_df3_comp = merged_df3_comp[!duplicated(merged_df3_comp$mutationinformation),] # a way + merged_df3_comp = merged_df3[!is.na(merged_df3$af),] # another way + cat("Checking nrows in merged_df3_comp") + if(nrow(merged_df3_comp) == (nrow(merged_df3) - na_count_df3)){ + cat("\nPASS: No. of rows match" + ,"\nDim of merged_df3_comp: " + ,"\nExpected no. of rows: ", nrow(merged_df3) - na_count_df3 + , "\nNo. of rows: ", nrow(merged_df3_comp) + , "\nNo. of cols: ", ncol(merged_df3_comp)) + }else{ + cat("FAIL: No. of rows mismatch" + ,"\nExpected no. of rows: ", nrow(merged_df3) - na_count_df3 + ,"\nGot no. of rows: ", nrow(merged_df3_comp)) + } +} else{ + print("Index mismatch for mychisq and kin ors. Aborting NA ommission") +} + +# alternate way of deriving merged_df3_comp +foo = merged_df3[!is.na(merged_df3$af),] +bar = merged_df3_comp[!duplicated(merged_df3_comp$mutationinformation),] +# compare dfs: foo and merged_df3_com +all.equal(foo, bar) +#summary(comparedf(foo, bar)) + +#============================================================== +################# +# OPTIONAL: write ALL 4 output files +################# +#outvars = c("merged_df2" +# , "merged_df3" +# , "merged_df2_comp" +# , "merged_df3_comp") + +#cat("Writing output files: " +# , "\nPath:", outdir) + +#for (i in outvars){ +# out_filename = paste0(i, ".csv") +# outfile = paste0(outdir, "/", out_filename) +# cat("Writing output file:" +# ,"\nFilename: ", out_filename,"\n") +# write.csv(get(i), outfile, row.names = FALSE) +# cat("Finished writing: ", outfile +# , "\nNo. of rows: ", nrow(get(i)) +# , "\nNo. of cols: ", ncol(get(i)), "\n") +#} + +#************************* +# clear variables +rm(foo, bar, gene_metadata + , in_filename_params, infile_params, merging_cols + , in_filename_gene_metadata, infile_gene_metadata + , merged_df2v2, merged_df2v3) +#************************* +##################################################################### +# Combining: LIG +##################################################################### + +#========================= +# Merges 5-8 +#========================= + +merged_df2_lig = merged_df2[merged_df2$ligand_distance<10,] +merged_df2_comp_lig = merged_df2_comp[merged_df2_comp$ligand_distance<10,] + +merged_df3_lig = merged_df3[merged_df3$ligand_distance<10,] +merged_df3_comp_lig = merged_df3_comp[merged_df3_comp$ligand_distance<10,] + +# sanity check +if (nrow(merged_df3_lig) == nrow(my_df_u_lig)){ + print("PASS: verified merged_df3_lig") +}else{ + cat(paste0("FAIL: nrow mismatch for merged_df3_lig" + , "\nExpected:", nrow(my_df_u_lig) + , "\nGot:", nrow(merged_df3_lig))) +} + +#========================================================================== +# end of script +##========================================================================= diff --git a/scripts/plotting/other_plots.R b/scripts/plotting/other_plots.R index 9fccb46..2ca3868 100644 --- a/scripts/plotting/other_plots.R +++ b/scripts/plotting/other_plots.R @@ -18,10 +18,8 @@ source("other_plots_data.R") #======= # output #======= -#dr_other_plots_combined = "dr_other_combined.svg" -#plot_dr_other_plots_combined = paste0(plotdir,"/", dr_other_plots_combined) - - +dr_other_plots_combined = "dr_other_muts.svg" +plot_dr_other_plots_combined = paste0(plotdir,"/", dr_other_plots_combined) ######################################################################## # end of data extraction and cleaning for plots # @@ -160,8 +158,8 @@ p3 # combine #=========================== #svg(plot_or_combined, width = 32, height = 12) - -theme_set(theme_gray()) # to preserve default theme +svg("test.svg", width = 25, height = 12) +#theme_set(theme_gray()) # to preserve default theme printFile = cowplot::plot_grid(plot_grid(p1, p2, p3 , nrow = 3