added colors for tiles for gid plots and appendix tables

This commit is contained in:
Tanushree Tunstall 2022-08-25 17:50:49 +01:00
parent ac72634b48
commit e4e8bd7278
3 changed files with 273 additions and 42 deletions

View file

@ -1,9 +1,10 @@
#!/usr/bin/env Rscript
#source("~/git/LSHTM_analysis/config/alr.R")
source("~/git/LSHTM_analysis/config/embb.R")
#source("~/git/LSHTM_analysis/config/katg.R")
#source("~/git/LSHTM_analysis/config/gid.R")
#!/usr/bin/env Rscript
#source("~/git/LSHTM_analysis/config/pnca.R")
#source("~/git/LSHTM_analysis/config/embb.R")
source("~/git/LSHTM_analysis/config/gid.R")
#source("~/git/LSHTM_analysis/config/alr.R")
#source("~/git/LSHTM_analysis/config/katg.R")
#source("~/git/LSHTM_analysis/config/rpob.R")
# get plottting dfs

View file

@ -1,27 +1,29 @@
#!/usr/bin/env Rscript
#source("~/git/LSHTM_analysis/config/alr.R")
source("~/git/LSHTM_analysis/config/embb.R")
#source("~/git/LSHTM_analysis/config/katg.R")
#source("~/git/LSHTM_analysis/config/gid.R")
#source("~/git/LSHTM_analysis/config/pnca.R")
#source("~/git/LSHTM_analysis/config/rpob.R")
# get plottting dfs
source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
source("~/git/LSHTM_analysis/scripts/plotting/plotting_colnames.R")
#=======
# output
#=======
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
outdir_stats = paste0(outdir_images,"stats/")
cat("\nOutput dir for stats:", outdir_stats)
###################################################################
# FIXME: ADD distance to NA when SP replies
geneL_normal = c("pnca")
geneL_na = c("gid", "rpob")
#geneL_na = c("gid", "rpob")
geneL_na_v2 = c("gid")
geneL_ppi2 = c("alr", "embb", "katg", "rpob")
geneL_both = c("rpob")
if (tolower(gene)%in%geneL_na_v2) {
gene_colnames = c("mcsm_na_affinity", "mcsm_na_outcome")
}
if (tolower(gene)%in%geneL_ppi2) {
gene_colnames = c("mcsm_ppi2_affinity", "mcsm_ppi2_outcome")
}
#from plotting_globals()
LigDist_colname
@ -31,22 +33,28 @@ naDist_colname
delta_symbol #delta_symbol = "\u0394"; delta_symbol
angstroms_symbol
cat("\nAffinity Distance colnames:", length(affinity_dist_colnames)
, "\nThese are:", affinity_dist_colnames)
#===========
# Data used
#===========
df3 = merged_df3
cols_to_output = c("mutationinformation"
, "position"
cols_to_output = c("position"
, "sensitivity"
, "mutationinformation"
, affinity_dist_colnames[1]
, "ligand_affinity_change"
, "ligand_outcome"
, "mmcsm_lig"
, "mmcsm_lig_outcome"
, affinity_dist_colnames[2]
, "mcsm_ppi2_affinity"
, "mcsm_ppi2_outcome"
# #, affinity_dist_colnames[3]
# , "mcsm_na_affinity"
# , "mcsm_na_outcome"
# #, "mcsm_ppi2_affinity"
# #, "mcsm_ppi2_outcome"
, gene_colnames
, "maf"
, "or_mychisq"
, "pval_fisher")
@ -54,7 +62,7 @@ cols_to_output = c("mutationinformation"
cols_to_output
df3_output = df3[, cols_to_output]
colnames(df3_output)
cat("\nSelecting columns:", length(colnames(df3_output)))
#===============================================
# Add COLS and rounding: adjusted P-values + MAF
#==============================================
@ -91,7 +99,6 @@ head(df3_output$af); head(df3_output$maf);head(df3_output$maf_percent)
#----------
df3_output$pval_fisher = round(df3_output$pval_fisher,2)
class(df3_output)
head(df3_output)
@ -100,8 +107,9 @@ head(df3_output)
####################################
df_lig = df3_output[df3_output[[LigDist_colname]]<DistCutOff,]
cols_to_output_lig = c("mutationinformation"
, "position"
cols_to_output_lig = c("position"
, "sensitivity"
, "mutationinformation"
, LigDist_colname
, "ligand_affinity_change"
, "ligand_outcome"
@ -122,8 +130,9 @@ Out_df_ligS = Out_df_lig[order(-Out_df_lig$or_mychisq, Out_df_lig$maf_percent),
#head(Out_df_ligS1); tail(Out_df_ligS1)
head(Out_df_ligS); tail(Out_df_ligS)
colsNames_to_output_lig = c("Mutation"
, "position"
colsNames_to_output_lig = c("position"
, "sensitivity"
, "Mutation"
, paste0("Lig-Dist (", angstroms_symbol, ")")
, "mCSM-ligand affinity"
, "mCSM ligand_outcome"
@ -150,18 +159,20 @@ write.csv(Out_df_ligS, Out_ligT, row.names = FALSE)
########################################################################
####################################
# Appendix: PPi2 affinity
# Appendix: ppi2 affinity
####################################
# Filtered data
df_ppi2 = df3_output[df3_output[[ppi2Dist_colname]]<DistCutOff,]
# select cols
cols_to_output_ppi2 = c("mutationinformation"
, "position"
cols_to_output_ppi2 = c("position"
, "sensitivity"
, "mutationinformation"
, ppi2Dist_colname
, "mcsm_ppi2_affinity"
, "mcsm_ppi2_outcome"
#, "mcsm_ppi2_affinity"
#, "mcsm_ppi2_outcome"
, gene_colnames
, "maf_percent"
, "or_mychisq"
, "pval_fisher"
@ -176,8 +187,9 @@ Out_df_ppi2 = df_ppi2[, cols_to_output_ppi2]
Out_df_ppi2S = Out_df_ppi2[order(-Out_df_ppi2$or_mychisq, Out_df_ppi2$maf_percent), ]
colsNames_to_output_ppi2 = c("Mutation"
, "position"
colsNames_to_output_ppi2 = c("position"
, "sensitivity"
, "Mutation"
, paste0("PPI2-Dist (", angstroms_symbol, ")")
, paste0("mCSM-PPI2 (", delta_symbol,delta_symbol,"G)")
, "mCSM-PPI2 outcome"
@ -199,3 +211,221 @@ Out_ppi2T = paste0(outdir_stats
cat("Output of PPI2 muts:", Out_ppi2T )
write.csv(Out_df_ppi2S, Out_ppi2T, row.names = FALSE)
#################################################################################
#################################################################################
#################################################################################
##########################################################
# higest or/maf and stability effects
###########################################################
# convert to percet
df3$maf_percent = df3$maf*100
cols_to_output_effects = c("position"
, "sensitivity"
, "mutationinformation"
, "avg_stability"
, "avg_stability_outcome"
, affinity_dist_colnames[1]
, "avg_lig_affinity"
, "avg_lig_affinity_outcome"
, affinity_dist_colnames[2]
, gene_colnames
, "maf_percent"
, "or_mychisq"
, "pval_fisher")
df3_effects = df3[, cols_to_output_effects]
nrow(df3_effects); ncol(df3_effects)
# add cols: p_adj_fdr and signif_fdr
df3_effects$p_adj_fdr = p.adjust(df3_effects$pval_fisher, method = "fdr")
df3_effects$signif_fdr = df3_effects$p_adj_fdr
df3_effects = dplyr::mutate(df3_effects
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
, signif_fdr <=0.0001 ~ '****'
, signif_fdr <=0.001 ~ '***'
, signif_fdr <=0.01 ~ '**'
, signif_fdr <0.05 ~ '*'
, TRUE ~ 'ns'))
# rounding
df3_effects$or_mychisq = round(df3_effects$or_mychisq,2)
df3_effects$p_adj_fdr = round(df3_effects$p_adj_fdr,2)
head(df3_effects)
#-------------------
# Highest OR and MAF
#-------------------
mut_hor = df3_effects[df3_effects$or_mychisq%in%(max(df3_effects$or_mychisq, na.rm = T)),]
mut_hmaf = df3_effects[df3_effects$maf_percent%in%(max(df3_effects$maf_percent, na.rm = T)),]
if (identical(colnames(mut_hor), colnames(mut_hmaf)) ){
# add cols
mut_hor$mutational_effect = "Mutation with highest OR"
mut_hmaf$mutational_effect = "Most frequent mutation"
cat("\nPass: or and maf")
}else{
quit("Abort: colnames or and maf mismatch")
}
#-------------------
# Avg stability
# most DD/SS: average stability
#-------------------
mut_h_avs_dd = df3_effects[df3_effects$avg_stability%in%(min(df3_effects$avg_stability, na.rm = T)),]
mut_h_avs_ss = df3_effects[df3_effects$avg_stability%in%(max(df3_effects$avg_stability, na.rm = T)),]
if (identical(colnames(mut_h_avs_dd), colnames(mut_h_avs_ss)) ){
# add cols
mut_h_avs_dd$mutational_effect = "Most Destabilising for protomer"
mut_h_avs_ss$mutational_effect = "Most Stabilising for protomer"
cat("\nPass : avg stability")
}else{
quit("Abort: colnames stability mismatch")
}
#-------------------
# Filtered columns
# most DD/SS: ligand
#-------------------
df3_effects_lig = df3_effects[df3_effects[[LigDist_colname]]<DistCutOff,]
nrow(df3_effects_lig)
mut_h_lig_dd = df3_effects_lig[df3_effects_lig$avg_lig_affinity%in%(min(df3_effects_lig$avg_lig_affinity, na.rm = T)),]
mut_h_lig_ss = df3_effects_lig[df3_effects_lig$avg_lig_affinity%in%(max(df3_effects_lig$avg_lig_affinity, na.rm = T)),]
if (identical(colnames(mut_h_lig_dd), colnames(mut_h_lig_ss)) ){
# add cols
mut_h_lig_dd$mutational_effect = "Most Destabilising for Ligand affinity"
mut_h_lig_ss$mutational_effect = "Most Stabilising for Ligand affinity"
cat("\nPass: avg ligand affinity")
}else{
quit("Abort: colnames lig mismatch")
}
#-------------------
# Filtered columns
# most DD/SS: NA
#-------------------
if (tolower(gene)%in%geneL_na_v2 ){
df3_effects_na = df3_effects[df3_effects[[naDist_colname]]<DistCutOff,]
nrow(df3_effects_na)
mut_h_na_dd = df3_effects_na[df3_effects_na$mcsm_na_affinity%in%(min(df3_effects_na$mcsm_na_affinity, na.rm = T)),]
mut_h_na_ss = df3_effects_na[df3_effects_na$mcsm_na_affinity%in%(max(df3_effects_na$mcsm_na_affinity, na.rm = T)),]
# add cols
mut_h_na_dd$mutational_effect = "Most Destabilising for NA affinity"
mut_h_na_ss$mutational_effect = "Most Stabilising for NA affinity"
if (identical(colnames(mut_h_na_dd), colnames(mut_h_na_ss)) ){
cat("\nPass 1: NCA")
}else{
quit("Abort: colnames nca mismatch")
}
if (identical(colnames(mut_h_na_dd), colnames(mut_h_lig_dd)) ){
cat("\nPass 2: NCA")
}else{
quit("Abort: colnames ppi2 mismatch")
}
#combine
gene_aff_combined = rbind(mut_h_na_dd, mut_h_na_ss)
}
#-------------------
# Filtered columns
# most DD/SS: ppi2
#-------------------
if (tolower(gene)%in%geneL_ppi2 ){
df3_effects_ppi2 = df3_effects[df3_effects[[ppi2Dist_colname]]<DistCutOff,]
nrow(df3_effects_ppi2)
mut_h_ppi2_dd = df3_effects_ppi2[df3_effects_ppi2$mcsm_ppi2_affinity%in%(min(df3_effects_ppi2$mcsm_ppi2_affinity, na.rm = T)),]
mut_h_ppi2_ss = df3_effects_ppi2[df3_effects_ppi2$mcsm_ppi2_affinity%in%(max(df3_effects_ppi2$mcsm_ppi2_affinity, na.rm = T)),]
# add cols
mut_h_ppi2_dd$mutational_effect = "Most Destabilising for PPI affinity"
mut_h_ppi2_dd$mutational_effect = "Most Stabilising for PPI affinity"
if (identical(colnames(mut_h_ppi2_dd), colnames(mut_h_ppi2_ss)) ){
cat("\nPass 1: ppi2")
}else{
quit("Abort: colnames ppi2 mismatch")
}
if (identical(colnames(mut_h_ppi2_dd), colnames(mut_h_lig_dd)) ){
cat("\nPass 2 : ppi2")
}else{
quit("Abort: colnames ppi2 mismatch")
}
#combine
gene_aff_combined = rbind(mut_h_na_dd, mut_h_na_ss)
}
#============
# final combine
#============
if ( identical(colnames(mut_hor), colnames(mut_h_lig_dd)) ){
cat("PASS: all")
combined_table = rbind(mut_hor, mut_hmaf,
mut_h_avs_dd, mut_h_avs_ss,
mut_h_lig_dd, mut_h_lig_ss,
gene_aff_combined)
cat("\nCombined table dim:", "\nnrow:", nrow(combined_table), "\nncol:", ncol(combined_table))
}else{
quit("Abort: colnames ppi2 mismatch")
}
# Assign pretty colnames
colnames(combined_table)
colsNames_combined_table = c("position"
, "sensitivity"
, "Mutation"
, paste0("Avg stability (", delta_symbol,delta_symbol,"G)")
, "avg stability outcome"
, paste0("Lig-Dist (", angstroms_symbol, ")")
, "Avg ligand affinity"
, "Ligand affinity outcome"
#, paste0("NA-Dist (", angstroms_symbol, ")")
#, paste0("mCSM-NA (", delta_symbol,delta_symbol,"G)")
#, "mCSM-NA outcome"
, paste0("PPI-Dist (", angstroms_symbol, ")")
, paste0("mCSM-PPI (", delta_symbol,delta_symbol,"G)")
, "mCSM-PPI outcome"
, paste0("MAF ","(%)")
, "Odds Ratio"
, "P-value"
, "Adjusted P-value"
, "Adjusted P-value significance"
, "Mutational effect")
if ( length(colnames(combined_table)) == length(colsNames_combined_table) ) {
cat("Assiging pretty colnames for output")
colnames(combined_table) <- colsNames_combined_table
#colnames(combined_table)
}else{
stop("\nAbort: No. of cols mismatch. Cannot assign pretty colnames for output")
}
#--------------------
# write output file: KS test within grpup
#----------------------
Out_combined_effectsT = paste0(outdir_stats
, tolower(gene)
, "_mut_effects.csv")
cat("Output of effects:", Out_combined_effectsT )
write.csv(combined_table, Out_combined_effectsT, row.names = FALSE)