tanu/LSHTM_analysis
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

added all scripts

Browse source
This commit is contained in:
Tanushree Tunstall 2022-08-27 23:12:00 +01:00
parent da8a1069a8
commit 7bed6a1e22
8 changed files with 1840 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

177
scripts/plotting/plotting_thesis/alr/DONOTUSEdm_om_plots_layout_alr.R Normal file
View file

@ -0,0 +1,177 @@
# # source dm_om_plots.R
# source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/dm_om_plots.R")
#
# ##### plots to combine ####
# duetP
# foldxP
# deepddgP
# dynamut2P
# genomicsP
# consurfP
# proveanP
# snap2P
# mcsmligP
# mcsmlig2P
# mcsmppi2P
#
#
# # Plot labels
# tit1 = "Stability changes"
# tit2 = "Genomic measure"
# tit3 = "Distance to partners"
# tit4 = "Evolutionary Conservation"
# tit5 = "Affinity changes"
# pt_size = 30
#
# theme_georgia <- function(...) {
# theme_gray(base_family = "sans", ...) +
# theme(plot.title = element_text(face = "bold"))
# }
#
#
# title_theme <- calc_element("plot.title", theme_georgia())
#
# pt1 = ggdraw() +
# draw_label(
# tit1,
# fontfamily = title_theme$family,
# fontface = title_theme$face,
# #size = title_theme$size
# size = pt_size
# )
#
# pt2 = ggdraw() +
# draw_label(
# tit2,
# fontfamily = title_theme$family,
# fontface = title_theme$face,
# size = pt_size
# )
#
# pt3 = ggdraw() +
# draw_label(
# tit3,
# fontfamily = title_theme$family,
# fontface = title_theme$face,
# size = pt_size
# )
#
# pt4 = ggdraw() +
# draw_label(
# tit4,
# fontfamily = title_theme$family,
# fontface = title_theme$face,
# size = pt_size
# )
#
#
# pt5 = ggdraw() +
# draw_label(
# tit5,
# fontfamily = title_theme$family,
# fontface = title_theme$face,
# size = pt_size
# )
#
# #======================
# # Output plot function
# #======================
# OutPlot_dm_om = function(x){
#
# # dist b/w plot title and plot
# relH_tp = c(0.08, 0.92)
#
# my_label_size = 25
# #----------------
# # Top panel
# #----------------
# top_panel = cowplot::plot_grid(
# cowplot::plot_grid(pt1,
# cowplot::plot_grid(duetP, foldxP, deepddgP, dynamut2P
# , nrow = 1
# , labels = c("A", "B", "C", "D")
# , label_size = my_label_size)
# , ncol = 1
# , rel_heights = relH_tp
# ),
# NULL,
# cowplot::plot_grid(pt2,
# cowplot::plot_grid(genomicsP
# , nrow = 1
# , labels = c("E")
# , label_size = my_label_size)
# , ncol = 1
# , rel_heights = relH_tp
# ),
# NULL,
# cowplot::plot_grid(pt3,
# cowplot::plot_grid( #distanceP
# distanceP_lig
# , distanceP_ppi2
# , nrow = 1
# , labels = c("F", "G")
# , label_size = my_label_size)
# , ncol = 1
# , rel_heights = relH_tp
# ),
# nrow = 1,
# rel_widths = c(2/7, 0.1/7, 0.5/7, 0.1/7, 1/7)
# )
#
# #----------------
# # Bottom panel
# #----------------
# bottom_panel = cowplot::plot_grid(
# cowplot::plot_grid(pt4,
# cowplot::plot_grid(consurfP, proveanP, snap2P
# , nrow = 1
# , labels = c("H", "I", "J")
# , label_size = my_label_size)
# , ncol = 1
# , rel_heights =relH_tp
# ),NULL,
# cowplot::plot_grid(pt5,
# cowplot::plot_grid(mcsmligP
# , mcsmlig2P
# , mcsmppi2P
# , nrow = 1
# , labels = c("K", "L", "M")
# , label_size = my_label_size)
# , ncol = 1
# , rel_heights = relH_tp
# ),NULL,
# nrow = 1,
# rel_widths = c(3/6,0.1/6,3/6, 0.1/6 )
# )
#
# #-------------------------------
# # combine: Top and Bottom panel
# #-------------------------------
# cowplot::plot_grid (top_panel, bottom_panel
# , nrow =2
# , rel_widths = c(1, 1)
# , align = "hv")
# }
#
# #=====================
# # OutPlot: svg and png
# #======================
# dm_om_combinedP = paste0(outdir_images
# ,tolower(gene)
# ,"_dm_om_all.svg")
#
# cat("DM OM plots with stats:", dm_om_combinedP)
# svg(dm_om_combinedP, width = 32, height = 18)
#
# OutPlot_dm_om()
# dev.off()
#
#
# dm_om_combinedP_png = paste0(outdir_images
# ,tolower(gene)
# ,"_dm_om_all.png")
# cat("DM OM plots with stats:", dm_om_combinedP_png)
# png(dm_om_combinedP_png, width = 32, height = 18, units = "in", res = 300)
#
# OutPlot_dm_om()
# dev.off()

198
scripts/plotting/plotting_thesis/alr/alr_ORandSNP_results.R Normal file
View file

@ -0,0 +1,198 @@
#!/usr/bin/env Rscript
#source("~/git/LSHTM_analysis/config/gid.R")
#source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
#=======
# output
#=======
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
###################################################################
# FIXME: ADD distance to NA when SP replies
geneL_normal = c("pnca")
geneL_na = c("gid", "rpob")
geneL_ppi2 = c("alr", "embb", "katg", "rpob")
# LigDist_colname # from globals used
# ppi2Dist_colname #from globals used
# naDist_colname #from globals used
common_cols = c("mutationinformation"
, "X5uhc_position"
, "X5uhc_offset"
, "position"
, "dst_mode"
, "mutation_info_labels"
, "sensitivity", dist_columns )
########################################
categ_cols_to_factor = grep( "_outcome|_info", colnames(merged_df3) )
fact_cols = colnames(merged_df3)[categ_cols_to_factor]
if (any(lapply(merged_df3[, fact_cols], class) == "character")){
cat("\nChanging", length(categ_cols_to_factor), "cols to factor")
merged_df3[, fact_cols] <- lapply(merged_df3[, fact_cols], as.factor)
if (all(lapply(merged_df3[, fact_cols], class) == "factor")){
cat("\nSuccessful: cols changed to factor")
}
}else{
cat("\nRequested cols aready factors")
}
cat("\ncols changed to factor are:\n", colnames(merged_df3)[categ_cols_to_factor] )
####################################
# merged_df3: NECESSARY pre-processing
###################################
#df3 = merged_df3
plot_cols = c("mutationinformation", "mutation_info_labels", "position", "dst_mode"
, all_cols)
all_cols = c(common_cols
, all_stability_cols
, all_affinity_cols
, all_conserv_cols)
# counting
foo = merged_df3[, c("mutationinformation"
, "wild_pos"
, "position"
, "sensitivity"
, "avg_lig_affinity"
, "avg_lig_affinity_scaled"
, "avg_lig_affinity_outcome"
, "ligand_distance"
, "ligand_affinity_change"
, "affinity_scaled"
, "ligand_outcome"
, "consurf_colour_rev"
, "consurf_outcome")]
table(foo$consurf_outcome)
foo2 = foo[foo$ligand_distance<10,]
table(foo2$ligand_outcome)
#############################
# wide plots SNP
# DRUG
length(aa_pos_drug); aa_pos_drug
drug = foo[foo$position%in%aa_pos_drug,]
drug$wild_pos
#CA
###############################################
# OR plot
bar = merged_df3[, c("mutationinformation"
, "wild_pos"
, "position"
, "sensitivity"
, affinity_dist_colnames
, "or_mychisq"
, "pval_fisher"
#, "pval_chisq"
, "neglog_pval_fisher"
, "log10_or_mychisq")]
# bar$p_adj_bonferroni = p.adjust(bar$pval_fisher, method = "bonferroni")
# bar$signif_bon = bar$p_adj_bonferroni
# bar = dplyr::mutate(bar
# , signif_bon = case_when(signif_bon == 0.05 ~ "."
# , signif_bon <=0.0001 ~ '****'
# , signif_bon <=0.001 ~ '***'
# , signif_bon <=0.01 ~ '**'
# , signif_bon <0.05 ~ '*'
# , TRUE ~ 'ns'))
bar$p_adj_fdr = p.adjust(bar$pval_fisher, method = "BH")
bar$signif_fdr = bar$p_adj_fdr
bar = dplyr::mutate(bar
, 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'))
# sort df
bar = bar[order(bar$or_mychisq, decreasing = T), ]
bar = bar[, c("mutationinformation"
, "wild_pos"
, "position"
, "sensitivity"
, affinity_dist_colnames
, "or_mychisq"
#, "pval_fisher"
#, "pval_chisq"
#, "neglog_pval_fisher"
#, "log10_or_mychisq"
#, "signif_bon"
, "p_adj_fdr"
, "signif_fdr")]
table(bar$sensitivity)
str(bar)
sen = bar[bar$or_mychisq<1,]
sen = na.omit(sen)
res = bar[bar$or_mychisq>1,]
res = na.omit(res)
# comp
bar_or = bar[!is.na(bar$or_mychisq),]
table(bar_or$sensitivity)
sen1 = bar_or[bar_or$or_mychisq<1,] # sen and res ~OR
res1 = bar_or[bar_or$or_mychisq>1,] # sen and res ~OR
# sanity check
if (nrow(bar_or) == nrow(sen1) + nrow(res1) ){
cat("\nPASS: df with or successfully sourced"
, "\nCalculating % of muts with OR>1")
}else{
stop("Abort: df with or numbers mimatch")
}
# percent for OR muts
pc_orR = nrow(res1)/(nrow(sen1) + nrow(res1)); pc_orR
cat("Number of R muts with OR>1:", nrow(res1)
, "\nPercentage of muts with OR>1 i.e resistant:"
, pc_orR *100 )
# muts with highest OR
head(bar_or$mutationinformation, 10)
# sort df
bar_or = bar_or[order(bar_or$or_mychisq
, bar_or$ligand_distance
#, bar_or$nca_dist
, bar_or$interface_dist
, decreasing = T), ]
nrow(bar_or)
bar_or$drug_site = ifelse(bar_or$position%in%aa_pos_drug, "drug", "no")
table(bar_or$drug_site)
bar_or$plp_site = ifelse(bar_or$position%in%aa_pos_plp, "plp", "no")
table(bar_or$plp_site)
top10_or = bar_or[1:10,]
# are these active sites
top10_or$position[top10_or$position%in%active_aa_pos]
# clostest most sig
bar_or_lig = bar_or[bar_or$ligand_distance<10,]
bar_or_lig = bar_or_lig[order(bar_or_lig$ligand_distance, -bar_or_lig$or_mychisq), ]
table(bar_or_lig$signif_fdr)
bar_or_ppi = bar_or[bar_or$interface_dist<10,]
bar_or_ppi = bar_or_ppi[order(bar_or_ppi$interface_dist, -bar_or_ppi$or_mychisq), ]
table(bar_or_ppi$signif_fdr)

459
scripts/plotting/plotting_thesis/alr/alr_appendix_tables.R Normal file
View file

@ -0,0 +1,459 @@
#!/usr/bin/env Rscript
#source("~/git/LSHTM_analysis/config/alr.R")
#source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.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)
###################################################################
geneL_normal = c("pnca")
#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
ppi2Dist_colname
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("position"
, "sensitivity"
, "mutationinformation"
, affinity_dist_colnames[1]
, "ligand_affinity_change"
, "ligand_outcome"
, "mmcsm_lig"
, "mmcsm_lig_outcome"
, affinity_dist_colnames[2]
# #, affinity_dist_colnames[3]
# , "mcsm_na_affinity"
# , "mcsm_na_outcome"
# #, "mcsm_nca_affinity"
# #, "mcsm_nca_outcome"
, gene_colnames
, "maf"
, "or_mychisq"
, "pval_fisher")
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
#==============================================
#-----------------------------
# adjusted P-values
#-----------------------------
# add cols: p_adj_fdr and signif_fdr
df3_output$p_adj_fdr = p.adjust(df3_output$pval_fisher, method = "fdr")
df3_output$signif_fdr = df3_output$p_adj_fdr
df3_output = dplyr::mutate(df3_output
, 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_output$or_mychisq = round(df3_output$or_mychisq,2)
df3_output$p_adj_fdr = round(df3_output$p_adj_fdr,2)
head(df3_output)
#----------
# MAF (%)
#----------
# add col maf_percent
df3_output$maf_percent = df3_output$maf*100
# rounding
df3_output$maf_percent = round(df3_output$maf_percent,2)
head(df3_output$af); head(df3_output$maf);head(df3_output$maf_percent)
#----------
# P-value
#----------
df3_output$pval_fisher = round(df3_output$pval_fisher,2)
class(df3_output)
head(df3_output)
####################################
# Appendix: ligand affinity
####################################
df_lig = df3_output[df3_output[[LigDist_colname]]<DistCutOff,]
cols_to_output_lig = c("position"
, "sensitivity"
, "mutationinformation"
, LigDist_colname
, "ligand_affinity_change"
, "ligand_outcome"
, "mmcsm_lig"
, "mmcsm_lig_outcome"
, "maf_percent"
, "or_mychisq"
, "pval_fisher"
, "p_adj_fdr"
, "signif_fdr")
# select cols
Out_df_lig = df_lig[, cols_to_output_lig]
# sort df by OR and then MAF: highest OR and highest MAF
#Out_df_ligS1 = Out_df_lig[order(Out_df_lig$or_mychisq, decreasing = T), ]
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("position"
, "sensitivity"
, "Mutation"
, paste0("Lig-Dist (", angstroms_symbol, ")")
, "mCSM-ligand affinity"
, "mCSM ligand_outcome"
, "mmCSM-ligand affinity"
, "mmCSM ligand_outcome"
, paste0("MAF ","(%)")
, "Odds Ratio"
, "P-value"
, "Adjusted P-value"
, "Adjusted P-value significance")
colnames(Out_df_ligS) = colsNames_to_output_lig
head(Out_df_ligS)
# ADD: active site annot
nrow(Out_df_ligS)
Out_df_ligS$drug_site = ifelse(Out_df_ligS$position%in%aa_pos_drug, "drug", "no")
table(Out_df_ligS$drug_site)
Out_df_ligS$plp_site = ifelse(Out_df_ligS$position%in%aa_pos_plp, "plp", "no")
table(Out_df_ligS$plp_site)
#--------------------
# write output file: KS test within grpup
#----------------------
Out_ligT = paste0(outdir_stats
, tolower(gene)
, "_lig_muts.csv")
cat("Output of Ligand muts:", Out_ligT )
write.csv(Out_df_ligS, Out_ligT, row.names = FALSE)
########################################################################
####################################
# Appendix: NA affinity or PPI2
# naDist_colname
####################################
# Filtered data
#df_nca = df3_output[df3_output[[naDist_colname]]<DistCutOff,]
df_nca = df3_output[df3_output[[ppi2Dist_colname]]<DistCutOff,]
# select cols
cols_to_output_nca = c("position"
, "sensitivity"
, "mutationinformation"
#, naDist_colname
, ppi2Dist_colname
, gene_colnames
, "maf_percent"
, "or_mychisq"
, "pval_fisher"
, "p_adj_fdr"
, "signif_fdr")
# extract output cols
Out_df_nca = df_nca[, cols_to_output_nca]
# sort df by OR and then MAF: Highest OR and Highest MAF
#Out_df_ncaS = Out_df_nca[order(Out_df_nca$or_mychisq, decreasing = T), ]
Out_df_ncaS = Out_df_nca[order(-Out_df_nca$or_mychisq, Out_df_nca$maf_percent), ]
colsNames_to_output_nca = c("position"
, "sensitivity"
, "Mutation"
#, 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")
colnames(Out_df_ncaS) = colsNames_to_output_nca
Out_df_ncaS
# ADD: active site annot
nrow(Out_df_ncaS)
Out_df_ncaS$drug_site = ifelse(Out_df_ncaS$position%in%aa_pos_drug, "drug", "no")
table(Out_df_ncaS$drug_site)
Out_df_ncaS$plp_site = ifelse(Out_df_ncaS$position%in%aa_pos_plp, "plp", "no")
table(Out_df_ncaS$plp_site)
#--------------------
# write output file: KS test within grpup
#----------------------
Out_ncaT = paste0(outdir_stats
, tolower(gene)
, "_na_muts.csv")
cat("Output of NA muts:", Out_ncaT )
write.csv(Out_df_ncaS, Out_ncaT, 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
# FIXME DUBIOUS as min and max can be both negative
#-------------------
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 = "CAUTION: Most DE/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_ss$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 1: colnames ppi2 mismatch")
}
if (identical(colnames(mut_h_ppi2_dd), colnames(mut_h_lig_dd)) ){
cat("\nPass 2 : ppi2")
}else{
quit("Abort 2: colnames ppi2 mismatch")
}
#combine
gene_aff_combined = rbind(mut_h_ppi2_dd, mut_h_ppi2_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")
}
nrow(combined_table)
combined_table$drug_site = ifelse(combined_table$position%in%aa_pos_drug, "drug", "no")
table(combined_table$drug_site)
combined_table$plp_site = ifelse(combined_table$position%in%aa_pos_plp, "plp", "no")
table(combined_table$plp_site)
#--------------------
# 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)

545
scripts/plotting/plotting_thesis/alr/alr_ks_test_all_PS.R Executable file
View file

@ -0,0 +1,545 @@
#!/usr/bin/env Rscript
#########################################################
# TASK: KS test for PS/DUET lineage distributions
#=======================================================================
#!/usr/bin/env Rscript
#source("~/git/LSHTM_analysis/config/alr.R")
# get plottting dfs
#source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
#=============
# Output
#=============
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
outdir_stats = paste0(outdir_images,"stats/")
# ks test by lineage
#ks_lineage = paste0(outdir, "/KS_lineage_all_muts.csv")
###########################
# Data for stats
# you need df2 or df2_comp
# since this is one-many relationship
###########################
# REASSIGNMENT
df2 = merged_df2[, colnames(merged_df2)%in%plotting_cols]
# quick checks
colnames(df2)
str(df2)
########################################################################
table(df2$lineage); str(df2$lineage)
# subset only lineages1-4
sel_lineages = c("L1"
, "L2"
, "L3"
, "L4")
# subset selected lineages
df_lin = subset(df2, subset = lineage %in% sel_lineages)
table(df_lin$lineage)
table(df_lin$sensitivity)
table(df_lin$lineage, df_lin$sensitivity)
#ensure lineage is a factor
#str(df_lin$lineage); str(df_lin$lineage_labels)
#df_lin$lineage = as.factor(df_lin$lineage)
#df_lin$lineage_labels = as.factor(df_lin$lineage)
table(df_lin$lineage); table(df_lin$lineage_labels)
#==============================
# Stats for average stability
#=============================
# individual: CHECKS
lin1 = df_lin[df_lin$lineage == "L1",]$avg_stability_scaled
lin2 = df_lin[df_lin$lineage == "L2",]$avg_stability_scaled
lin3 = df_lin[df_lin$lineage == "L3",]$avg_stability_scaled
lin4 = df_lin[df_lin$lineage == "L4",]$avg_stability_scaled
ks.test(lin1, lin4)
ks.test(df_lin$avg_stability_scaled[df_lin$lineage == "L1"]
, df_lin$avg_stability_scaled[df_lin$lineage == "L4"])
#=======================================================================
my_lineages = levels(factor(df_lin$lineage)); my_lineages
#=======================================================================
# Loop
#0 : < 2.2e-16
#=====================
# Lineage 1 comparisons
#=====================
my_lin1 = "L1"
#my_lineages_comp_l1 = c("L2", "L3", "L4")
my_lineages_comp_l1 = my_lineages[-match(my_lin1, my_lineages)]
ks_df_l1 = data.frame()
for (i in my_lineages_comp_l1){
cat(i)
l1_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
, n_samples = NA
, n_samples_total = NA)
lineage_comp = paste0(my_lin1, " vs ", i)
n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin1,])
n_samples_i = nrow(df_lin[df_lin$lineage == i,])
n_samples_all = paste0(my_lin1, "(", n_samples_lin, ")"
, ", "
, i, "(", n_samples_i, ")")
n_samples_total = n_samples_lin + n_samples_i
ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$method
ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic
ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
l1_df$comparison = lineage_comp
l1_df$method = ks_method
l1_df$ks_statistic = ks_statistic[[1]]
l1_df$ks_pvalue = ks_pvalue
l1_df$n_samples = n_samples_all
l1_df$n_samples_total= n_samples_total
print(l1_df)
ks_df_l1 = rbind(ks_df_l1,l1_df)
}
ks_df_l1
# adjusted p-value: bonferroni
ks_df_l1$p_adj_bonferroni = p.adjust(ks_df_l1$ks_pvalue, method = "bonferroni")
ks_df_l1$signif_bon = ks_df_l1$p_adj_bonferroni
ks_df_l1 = dplyr::mutate(ks_df_l1
, signif_bon = case_when(signif_bon == 0.05 ~ "."
, signif_bon <=0.0001 ~ '****'
, signif_bon <=0.001 ~ '***'
, signif_bon <=0.01 ~ '**'
, signif_bon <0.05 ~ '*'
, TRUE ~ 'ns'))
# adjusted p-value:fdr
ks_df_l1$p_adj_fdr = p.adjust(ks_df_l1$ks_pvalue, method = "fdr")
ks_df_l1$signif_fdr = ks_df_l1$p_adj_fdr
ks_df_l1 = dplyr::mutate(ks_df_l1
, 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'))
ks_df_l1
#####################################################################
#=====================
# Lineage 2 comparisons
#=====================
my_lin2 = "L2"
#my_lineages_comp_l2 = c("L1", lineage3", "L4")
my_lineages_comp_l2 = my_lineages[-match(my_lin2, my_lineages)]
ks_df_l2 = data.frame()
for (i in my_lineages_comp_l2){
cat(i)
l2_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
, n_samples = NA
, n_samples_total = NA)
lineage_comp = paste0(my_lin2, " vs ", i)
n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin2,])
n_samples_i = nrow(df_lin[df_lin$lineage == i,])
n_samples_all = paste0(my_lin2, "(", n_samples_lin, ")"
, ", "
, i, "(", n_samples_i, ")")
n_samples_total = n_samples_lin + n_samples_i
ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$method
ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic
ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
l2_df$comparison = lineage_comp
l2_df$method = ks_method
l2_df$ks_statistic = ks_statistic[[1]]
l2_df$ks_pvalue = ks_pvalue
l2_df$n_samples = n_samples_all
l2_df$n_samples_total = n_samples_total
print(l2_df)
ks_df_l2 = rbind(ks_df_l2, l2_df)
}
# adjusted p-value: bonferroni
ks_df_l2$p_adj_bonferroni = p.adjust(ks_df_l2$ks_pvalue, method = "bonferroni")
ks_df_l2$signif_bon = ks_df_l2$p_adj_bonferroni
ks_df_l2 = dplyr::mutate(ks_df_l2
, signif_bon = case_when(signif_bon == 0.05 ~ "."
, signif_bon <=0.0001 ~ '****'
, signif_bon <=0.001 ~ '***'
, signif_bon <=0.01 ~ '**'
, signif_bon <0.05 ~ '*'
, TRUE ~ 'ns'))
# adjusted p-value:fdr
ks_df_l2$p_adj_fdr = p.adjust(ks_df_l2$ks_pvalue, method = "fdr")
ks_df_l2$signif_fdr = ks_df_l2$p_adj_fdr
ks_df_l2 = dplyr::mutate(ks_df_l2
, 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'))
ks_df_l2
#=====================
# Lineage 3 comparisons
#=====================
my_lin3 = "L3"
#my_lineages_comp_l3 = c("L1", lineage2", "L4")
my_lineages_comp_l3 = my_lineages[-match(my_lin3, my_lineages)]
ks_df_l3 = data.frame()
for (i in my_lineages_comp_l3){
cat(i)
l3_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
, n_samples = NA
, n_samples_total = NA)
lineage_comp = paste0(my_lin3, " vs ", i)
n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin3,])
n_samples_i = nrow(df_lin[df_lin$lineage == i,])
n_samples_all = paste0(my_lin3, "(", n_samples_lin, ")"
, ", "
, i, "(", n_samples_i, ")")
n_samples_total = n_samples_lin + n_samples_i
ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$method
ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic
ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3]
, df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
l3_df$comparison = lineage_comp
l3_df$method = ks_method
l3_df$ks_statistic = ks_statistic[[1]]
l3_df$ks_pvalue = ks_pvalue
l3_df$n_samples = n_samples_all
l3_df$n_samples_total = n_samples_total
print(l3_df)
ks_df_l3 = rbind(ks_df_l3, l3_df)
}
# adjusted p-value: bonferroni
ks_df_l3$p_adj_bonferroni = p.adjust(ks_df_l3$ks_pvalue, method = "bonferroni")
ks_df_l3$signif_bon = ks_df_l3$p_adj_bonferroni
ks_df_l3 = dplyr::mutate(ks_df_l3
, signif_bon = case_when(signif_bon == 0.05 ~ "."
, signif_bon <=0.0001 ~ '****'
, signif_bon <=0.001 ~ '***'
, signif_bon <=0.01 ~ '**'
, signif_bon <0.05 ~ '*'
, TRUE ~ 'ns'))
# adjusted p-value:fdr
ks_df_l3$p_adj_fdr = p.adjust(ks_df_l3$ks_pvalue, method = "fdr")
ks_df_l3$signif_fdr = ks_df_l3$p_adj_fdr
ks_df_l3 = dplyr::mutate(ks_df_l3
, 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'))
ks_df_l3
####################################################################
# combine all 3 ks_dfs
n_dfs = 3
if ( all.equal(nrow(ks_df_l1), nrow(ks_df_l2), nrow(ks_df_l3)) &&
all.equal(ncol(ks_df_l1), ncol(ks_df_l2), ncol(ks_df_l3)) ){
cat("\nPASS: Calculating expected rows and cols for sanity checks on combined_dfs")
expected_rows = nrow(ks_df_l1) * n_dfs
expected_cols = ncol(ks_df_l1)
ks_df_combined = rbind(ks_df_l1, ks_df_l2, ks_df_l3)
if ( nrow(ks_df_combined) == expected_rows && ncol(ks_df_combined) == expected_cols ){
cat("\nPASS: combined df successfully created"
, "\nnrow combined_df:", nrow(ks_df_combined)
, "\nncol combined_df:", ncol(ks_df_combined))
}
else{
cat("\nFAIL: Dim mismatch"
, "\nExpected rows:", expected_rows
, "\nGot:", nrow(ks_df_combined)
, "\nExpected cols:", expected_cols
, "\nGot:", ncol(ks_df_combined))
}
}else{
cat("\nFAIL: Could not generate expected_rows and/or expected_cols"
, "\nCheck hardcoded value of n_dfs")
}
#----------------------------
# ADD extra cols: formatting
#----------------------------
# adding pvalue_signif
ks_df_combined$pvalue_signif = ks_df_combined$ks_pvalue
str(ks_df_combined$pvalue_signif)
ks_df_combined = dplyr::mutate(ks_df_combined
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
, pvalue_signif <=0.0001 ~ '****'
, pvalue_signif <=0.001 ~ '***'
, pvalue_signif <=0.01 ~ '**'
, pvalue_signif <0.05 ~ '*'
, TRUE ~ 'ns'))
# Remove duplicates
rows_to_remove = c("L2 vs L1", "L3 vs L1", "L3 vs L2")
ks_df_combined_f = ks_df_combined[-match(rows_to_remove, ks_df_combined$comparison),]
#################################################################################
#================================
# R vs S distribution: Overall
#=================================
df_lin_R = df_lin[df_lin$sensitivity == "R",]
df_lin_S = df_lin[df_lin$sensitivity == "S",]
overall_RS_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
, n_samples = NA
, n_samples_total = NA)
comp_type = "overall R vs S distributions"
n_samples_R = nrow(df_lin_R)
n_samples_S = nrow(df_lin_S)
n_samples_all = paste0("R(n=", n_samples_R, ")" , " ", "S(n=", n_samples_S, ")")
n_samples_total = n_samples_R+n_samples_S
ks.test(df_lin_R$avg_stability_scaled
, df_lin_S$avg_stability_scaled)
ks_method = ks.test(df_lin_R$avg_stability_scaled
, df_lin_S$avg_stability_scaled)$method
ks_statistic = ks.test(df_lin_R$avg_stability_scaled
, df_lin_S$avg_stability_scaled)$statistic
ks_pvalue = ks.test(df_lin_R$avg_stability_scaled
, df_lin_S$avg_stability_scaled)$p.value
overall_RS_df$comparison = comp_type
overall_RS_df$method = ks_method
overall_RS_df$ks_statistic = ks_statistic
overall_RS_df$ks_pvalue = ks_pvalue
overall_RS_df$n_samples = n_samples_all
overall_RS_df$n_samples_total= n_samples_total
#----------------------------
# ADD extra cols
#----------------------------
# adjusted p-value: bonferroni
overall_RS_df$p_adj_bonferroni = p.adjust(overall_RS_df$ks_pvalue, method = "bonferroni")
overall_RS_df$signif_bon = overall_RS_df$p_adj_bonferroni
overall_RS_df = dplyr::mutate(overall_RS_df
, signif_bon = case_when(signif_bon == 0.05 ~ "."
, signif_bon <=0.0001 ~ '****'
, signif_bon <=0.001 ~ '***'
, signif_bon <=0.01 ~ '**'
, signif_bon <0.05 ~ '*'
, TRUE ~ 'ns'))
# adjusted p-value:fdr
overall_RS_df$p_adj_fdr = p.adjust(overall_RS_df$ks_pvalue, method = "fdr")
overall_RS_df$signif_fdr = overall_RS_df$p_adj_fdr
overall_RS_df = dplyr::mutate(overall_RS_df
, 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'))
# unadjusted p-values
overall_RS_df$pvalue_signif = overall_RS_df$ks_pvalue
overall_RS_df = dplyr::mutate(overall_RS_df
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
, pvalue_signif <=0.0001 ~ '****'
, pvalue_signif <=0.001 ~ '***'
, pvalue_signif <=0.01 ~ '**'
, pvalue_signif <0.05 ~ '*'
, TRUE ~ 'ns'))
#####################################################################
if (all(colnames(ks_df_combined_f) == colnames(overall_RS_df))){
cat("\nPASS:combining KS test results")
}else{
stop("\nAbort: Cannot combine results for KS test")
}
ks_df_combined_f2 = rbind(ks_df_combined_f, overall_RS_df)
# ADD extra cols
ks_df_combined_f2$ks_comp_type = "between_lineages"
ks_df_combined_f2$gene_name = tolower(gene)
ks_df_combined_f2
###########################################################################
#=================================
# Within lineage R vs S: MANUAL
#=================================
lin1 = df_lin[df_lin$lineage == my_lin1,]
ks.test(lin1$avg_stability_scaled[lin1$sensitivity == "R"]
, lin1$avg_stability_scaled[lin1$sensitivity == "S"])
lin2 = df_lin[df_lin$lineage == my_lin2,]
ks.test(lin2$avg_stability_scaled[lin2$sensitivity == "R"]
, lin2$avg_stability_scaled[lin2$sensitivity == "S"])
lin3 = df_lin[df_lin$lineage == "L3",]
ks.test(lin3$avg_stability_scaled[lin3$sensitivity == "R"]
, lin3$avg_stability_scaled[lin3$sensitivity == "S"])
lin4 = df_lin[df_lin$lineage == "L4",]
ks.test(lin4$avg_stability_scaled[lin4$sensitivity == "R"]
, lin4$avg_stability_scaled[lin4$sensitivity == "S"])
###################################################################
#=======================
# Within lineage R vs S: LOOP
#=======================
all_within_lin_df = data.frame()
for (i in c(my_lin1
, my_lin2
, my_lin3
)){
#cat(i)
within_lin_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
, n_samples = NA
, n_samples_total = NA)
within_lineage_comp = paste0(i, ": R vs S")
my_lin_df = df_lin[df_lin$lineage == i,]
lin_R_n = nrow(my_lin_df[my_lin_df$sensitivity == "R",])
lin_S_n = nrow(my_lin_df[my_lin_df$sensitivity == "S",])
lin_total_n = lin_R_n+lin_S_n
n_samples_all = paste0("R(n=", lin_R_n, ")" , ", ", "S(n=", lin_S_n, ")")
ks_method = ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"]
, my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$method
ks_statistic = ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"]
, my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$statistic
ks_pvalue =ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"]
, my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$p.value
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
within_lin_df$comparison = within_lineage_comp
within_lin_df$method = ks_method
within_lin_df$ks_statistic = ks_statistic[[1]]
within_lin_df$ks_pvalue = ks_pvalue
within_lin_df$n_samples = n_samples_all
within_lin_df$n_samples_total=lin_total_n
print(my_lin_df)
all_within_lin_df = rbind(all_within_lin_df, within_lin_df)
}
all_within_lin_df
#----------------------------
# ADD extra cols: formatting
#----------------------------
# adjusted p-value: bonferroni
all_within_lin_df$p_adj_bonferroni = p.adjust(all_within_lin_df$ks_pvalue, method = "bonferroni")
all_within_lin_df$signif_bon = all_within_lin_df$p_adj_bonferroni
all_within_lin_df = dplyr::mutate(all_within_lin_df
, signif_bon = case_when(signif_bon == 0.05 ~ "."
, signif_bon <=0.0001 ~ '****'
, signif_bon <=0.001 ~ '***'
, signif_bon <=0.01 ~ '**'
, signif_bon <0.05 ~ '*'
, TRUE ~ 'ns'))
# adjusted p-value:fdr
all_within_lin_df$p_adj_fdr = p.adjust(all_within_lin_df$ks_pvalue, method = "fdr")
all_within_lin_df$signif_fdr = all_within_lin_df$p_adj_fdr
all_within_lin_df = dplyr::mutate(all_within_lin_df
, 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'))
# unadjusted p-value
all_within_lin_df$pvalue_signif = all_within_lin_df$ks_pvalue
str(all_within_lin_df$pvalue_signif)
all_within_lin_df = dplyr::mutate(all_within_lin_df
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
, pvalue_signif <=0.0001 ~ '****'
, pvalue_signif <=0.001 ~ '***'
, pvalue_signif <=0.01 ~ '**'
, pvalue_signif <0.05 ~ '*'
, TRUE ~ 'ns'))
# ADD info cols
all_within_lin_df$ks_comp_type = "within_lineages"
all_within_lin_df$gene_name = tolower(gene)
##################################################################
if (all(colnames(ks_df_combined_f2) == colnames(all_within_lin_df))){
cat("\nPASS:combining KS test results")
}else{
stop("\nAbort: Cannot combine results for KS test")
}
ks_df_combined_all = rbind(ks_df_combined_f2, all_within_lin_df)
#--------------------
# write output file: KS test within grpup
#----------------------
Out_ks_all = paste0(outdir_stats
, tolower(gene)
, "_ks_lineage_all_comp.csv")
cat("Output of KS test all comparisons:", Out_ks_all )
write.csv(ks_df_combined_all, Out_ks_all, row.names = FALSE)

165
scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist.R Normal file
View file

@ -0,0 +1,165 @@
#!/usr/bin/env Rscript
#########################################################
# TASK: Lineage plots [merged_df2]
# Count
# Diversity
# Average stability dist
# Avergae affinity dist: optional
#########################################################
#=======
# output
#=======
# outdir_images = paste0("~/git/Writing/thesis/images/results/"
# , tolower(gene), "/")
# cat("plots will output to:", outdir_images)
#########################################################
#===============
#Quick numbers checks
#===============
nsample_lin = merged_df2[merged_df2$lineage%in%c("L1", "L2", "L3", "L4"),]
if ( all(table(nsample_lin$sensitivity)== table(nsample_lin$mutation_info_labels)) ){
cat("\nTotal no. of samples belonging to L1-l4 for", gene,":", nrow(nsample_lin)
, "\nCounting R and S samples")
if( sum(table(nsample_lin$sensitivity)) == nrow(nsample_lin) ){
cat("\nPASSNumbers cross checked:")
print(table(nsample_lin$sensitivity))
}
}else{
stop("Abort: Numbers mismatch. Please check")
}
########################################################################
###################################################
# Lineage barplots #
###################################################
my_xats = 13 # x axis text size # were 25
my_yats = 13# y axis text sized_lab_size
my_xals = 13 # x axis label size
my_yals = 13 # y axis label size
my_lls = 13 # legend label size
d_lab_size = 5
#===============================
# lineage sample and SNP count
#===============================
lin_countP = lin_count_bp(lf_data = lineage_dfL[['lin_lf']]
, all_lineages = F
, x_categ = "sel_lineages"
, y_count = "p_count"
, use_lineages = c("L1", "L2", "L3", "L4")
, bar_fill_categ = "count_categ"
, display_label_col = "p_count"
, bar_stat_stype = "identity"
, d_lab_size = d_lab_size
, d_lab_col = "black"
, my_xats = my_xats # x axis text size
, my_yats = my_yats # y axis text sized_lab_size
, my_xals = my_xals # x axis label size
, my_yals = my_yals # y axis label size
, my_lls = my_lls # legend label size
, bar_col_labels = c("nsSNPs", "Total Samples")
, bar_col_values = c("grey50", "gray75")
, bar_leg_name = ""
, leg_location = "top"
, y_log10 = F
, y_scale_percent = FALSE
, y_label = c("Count")
)
lin_countP
#===============================
# lineage SNP diversity count
#===============================
lin_diversityP = lin_count_bp_diversity(lf_data = lineage_dfL[['lin_wf']]
, x_categ = "sel_lineages"
, y_count = "snp_diversity"
#, all_lineages = F
, use_lineages = c("L1", "L2", "L3", "L4")
, display_label_col = "snp_diversity_f"
, bar_stat_stype = "identity"
, x_lab_angle = 90
, d_lab_size = d_lab_size
, my_xats = my_xats # x axis text size
, my_yats = my_yats # y axis text sized_lab_size
, my_xals = my_xals # x axis label size
, my_yals = my_yals # y axis label size
, my_lls = my_lls # legend label size
, y_log10 = F
, y_scale_percent = F
, leg_location = "top"
, y_label = "Percent" #"SNP diversity"
, bp_plot_title = "nsSNP diversity"
, title_colour = "black" #"chocolate4"
, subtitle_text = NULL
, sts = 10
, subtitle_colour = "#350E20FF")
lin_diversityP
###################################################
# Stability dist #
###################################################
# scaled_cols_stability = c("duet_scaled"
# , "deepddg_scaled"
# , "ddg_dynamut2_scaled"
# , "foldx_scaled"
# , "avg_stability_scaled")
my_ats = 13 # x axis text size # were 25
my_als = 13# y axis text sized_lab_size
my_leg_ts = 13 # x axis label size
my_leg_title =11 # y axis label size
my_strip_ts = 13 #
my_xlabel = paste0("Average stability ", "(", stability_suffix, ")"); my_xlabel
#plotdf = merged_df2[merged_df2$lineage%in%c("L1", "L2", "L3", "L4"),]
linP_dm_om = lineage_distP(merged_df2
, with_facet = F
, x_axis = "avg_stability_scaled"
, y_axis = "lineage_labels"
, x_lab = my_xlabel
, use_lineages = c("L1", "L2", "L3", "L4")
#, fill_categ = "mutation_info_orig", fill_categ_cols = c("#E69F00", "#999999")
, fill_categ = "sensitivity"
, fill_categ_cols = c("red", "blue")
, label_categories = c("Resistant", "Sensitive")
, leg_label = "Mutation group"
, my_ats = my_ats # axis text size
, my_als = my_als # axis label size
, my_leg_ts = my_leg_ts
, my_leg_title = my_leg_title
, my_strip_ts = my_strip_ts
, alpha = 0.56
)
linP_dm_om
###################################################
# Affinity dist [OPTIONAL] #
###################################################
# scaled_cols_affinity = c("affinity_scaled"
# , "mmcsm_lig_scaled"
# , "mcsm_ppi2_scaled"
# , "mcsm_na_scaled"
# , "avg_lig_affinity_scaled")
# lineage_distP(merged_df2
# , with_facet = F
# , x_axis = "avg_lig_affinity_scaled"
# , y_axis = "lineage_labels"
# , x_lab = my_xlabel
# , use_lineages = c("L1", "L2", "L3", "L4")
# #, fill_categ = "mutation_info_orig", fill_categ_cols = c("#E69F00", "#999999")
# , fill_categ = "sensitivity"
# , fill_categ_cols = c("red", "blue")
# , label_categories = c("Resistant", "Sensitive")
# , leg_label = "Mutation group"
# , my_ats = 22 # axis text size
# , my_als = 22 # axis label size
# , my_leg_ts = 22
# , my_leg_title = 22
# , my_strip_ts = 22
# , alpha = 0.56
# )

75
scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist_layout.R Normal file
View file

@ -0,0 +1,75 @@
#!/usr/bin/env Rscript
source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist.R")
###########################################
# TASK: generate plots for lineage
# Individual plots in
#lineage_bp_both.R
#linage_dist_ens_stability.R
###########################################
#### just A+B ####
# png
my_label_size = 12
linPlots_combined = paste0(outdir_images
, tolower(gene)
,"_linP_AB_combined.png")
cat("\nOutput plot:", linPlots_combined)
png(linPlots_combined, width = 9, height = 6, units = "in", res = 300)
cowplot::plot_grid(lin_countP, lin_diversityP
, ncol = 2
, labels = "AUTO"
, label_size = my_label_size)
dev.off()
# svg
# linPlots_combined = paste0(outdir_images
# , tolower(gene)
# ,"_linP_combined.svg")
#
# cat("\nOutput plot:", linPlots_combined)
# svg(linPlots_combined, width = 18, height = 12)
#
# cowplot::plot_grid(
# cowplot::plot_grid(lin_countP, lin_diversityP
# , nrow = 2
# , rel_heights = c(1.2,1)
# , labels = "AUTO"
# , label_size = my_label_size),
# NULL,
# linP_dm_om,
# nrow = 1,
# labels = c("", "", "C"),
# label_size = my_label_size,
# rel_widths = c(35, 3, 52)
# )
# dev.off()
# png
# my_label_size = 12
# linPlots_combined = paste0(outdir_images
# , tolower(gene)
# ,"_linP_combined.png")
#
# cat("\nOutput plot:", linPlots_combined)
# png(linPlots_combined, width = 9, height = 6, units = "in" ,res = 300)
#
# cowplot::plot_grid(
# cowplot::plot_grid(lin_countP, lin_diversityP
# , nrow = 2
# , rel_heights = c(1.2,1)
# , labels = "AUTO"
# , label_size = my_label_size),
# NULL,
# linP_dm_om,
# nrow = 1,
# labels = c("", "", "C"),
# label_size = my_label_size,
# rel_widths = c(35, 3, 52)
# )
# dev.off()

165
scripts/plotting/plotting_thesis/lineage_bp_dist.R Normal file
View file

@ -0,0 +1,165 @@
#!/usr/bin/env Rscript
#########################################################
# TASK: Lineage plots [merged_df2]
# Count
# Diversity
# Average stability dist
# Avergae affinity dist: optional
#########################################################
#=======
# output
#=======
# outdir_images = paste0("~/git/Writing/thesis/images/results/"
# , tolower(gene), "/")
# cat("plots will output to:", outdir_images)
#########################################################
#===============
#Quick numbers checks
#===============
nsample_lin = merged_df2[merged_df2$lineage%in%c("L1", "L2", "L3", "L4"),]
if ( all(table(nsample_lin$sensitivity)== table(nsample_lin$mutation_info_labels)) ){
cat("\nTotal no. of samples belonging to L1-l4 for", gene,":", nrow(nsample_lin)
, "\nCounting R and S samples")
if( sum(table(nsample_lin$sensitivity)) == nrow(nsample_lin) ){
cat("\nPASSNumbers cross checked:")
print(table(nsample_lin$sensitivity))
}
}else{
stop("Abort: Numbers mismatch. Please check")
}
########################################################################
###################################################
# Lineage barplots #
###################################################
my_xats = 8 # x axis text size # were 25
my_yats = 8# y axis text sized_lab_size
my_xals = 8 # x axis label size
my_yals = 8 # y axis label size
my_lls = 8 # legend label size
d_lab_size = 2.3
#===============================
# lineage sample and SNP count
#===============================
lin_countP = lin_count_bp(lf_data = lineage_dfL[['lin_lf']]
, all_lineages = F
, x_categ = "sel_lineages"
, y_count = "p_count"
, use_lineages = c("L1", "L2", "L3", "L4")
, bar_fill_categ = "count_categ"
, display_label_col = "p_count"
, bar_stat_stype = "identity"
, d_lab_size = d_lab_size
, d_lab_col = "black"
, my_xats = my_xats # x axis text size
, my_yats = my_yats # y axis text sized_lab_size
, my_xals = my_xals # x axis label size
, my_yals = my_yals # y axis label size
, my_lls = my_lls # legend label size
, bar_col_labels = c("nsSNPs", "Total Samples")
, bar_col_values = c("grey50", "gray75")
, bar_leg_name = ""
, leg_location = "top"
, y_log10 = F
, y_scale_percent = FALSE
, y_label = c("Count")
)
lin_countP
#===============================
# lineage SNP diversity count
#===============================
lin_diversityP = lin_count_bp_diversity(lf_data = lineage_dfL[['lin_wf']]
, x_categ = "sel_lineages"
, y_count = "snp_diversity"
#, all_lineages = F
, use_lineages = c("L1", "L2", "L3", "L4")
, display_label_col = "snp_diversity_f"
, bar_stat_stype = "identity"
, x_lab_angle = 90
, d_lab_size = d_lab_size
, my_xats = my_xats # x axis text size
, my_yats = my_yats # y axis text sized_lab_size
, my_xals = my_xals # x axis label size
, my_yals = my_yals # y axis label size
, my_lls = my_lls # legend label size
, y_log10 = F
, y_scale_percent = F
, leg_location = "top"
, y_label = "Percent" #"SNP diversity"
, bp_plot_title = "nsSNP diversity"
, title_colour = "black" #"chocolate4"
, subtitle_text = NULL
, sts = 10
, subtitle_colour = "#350E20FF")
lin_diversityP
###################################################
# Stability dist #
###################################################
# scaled_cols_stability = c("duet_scaled"
# , "deepddg_scaled"
# , "ddg_dynamut2_scaled"
# , "foldx_scaled"
# , "avg_stability_scaled")
my_ats = 8 # x axis text size # were 25
my_als = 8# y axis text sized_lab_size
my_leg_ts = 8 # x axis label size
my_leg_title = 8 # y axis label size
my_strip_ts = 8 #
my_xlabel = paste0("Average stability ", "(", stability_suffix, ")"); my_xlabel
#plotdf = merged_df2[merged_df2$lineage%in%c("L1", "L2", "L3", "L4"),]
linP_dm_om = lineage_distP(merged_df2
, with_facet = F
, x_axis = "avg_stability_scaled"
, y_axis = "lineage_labels"
, x_lab = my_xlabel
, use_lineages = c("L1", "L2", "L3", "L4")
#, fill_categ = "mutation_info_orig", fill_categ_cols = c("#E69F00", "#999999")
, fill_categ = "sensitivity"
, fill_categ_cols = c("red", "blue")
, label_categories = c("Resistant", "Sensitive")
, leg_label = "Mutation group"
, my_ats = my_ats # axis text size
, my_als = my_als # axis label size
, my_leg_ts = my_leg_ts
, my_leg_title = my_leg_title
, my_strip_ts = my_strip_ts
, alpha = 0.56
)
linP_dm_om
###################################################
# Affinity dist [OPTIONAL] #
###################################################
# scaled_cols_affinity = c("affinity_scaled"
# , "mmcsm_lig_scaled"
# , "mcsm_ppi2_scaled"
# , "mcsm_na_scaled"
# , "avg_lig_affinity_scaled")
# lineage_distP(merged_df2
# , with_facet = F
# , x_axis = "avg_lig_affinity_scaled"
# , y_axis = "lineage_labels"
# , x_lab = my_xlabel
# , use_lineages = c("L1", "L2", "L3", "L4")
# #, fill_categ = "mutation_info_orig", fill_categ_cols = c("#E69F00", "#999999")
# , fill_categ = "sensitivity"
# , fill_categ_cols = c("red", "blue")
# , label_categories = c("Resistant", "Sensitive")
# , leg_label = "Mutation group"
# , my_ats = 22 # axis text size
# , my_als = 22 # axis label size
# , my_leg_ts = 22
# , my_leg_title = 22
# , my_strip_ts = 22
# , alpha = 0.56
# )

56
scripts/plotting/plotting_thesis/lineage_bp_dist_layout.R Normal file
View file

@ -0,0 +1,56 @@
#!/usr/bin/env Rscript
source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/linage_bp_dist.R")
###########################################
# TASK: generate plots for lineage
# Individual plots in
#lineage_bp_both.R
#linage_dist_ens_stability.R
###########################################
# svg
# linPlots_combined = paste0(outdir_images
# , tolower(gene)
# ,"_linP_combined.svg")
#
# cat("\nOutput plot:", linPlots_combined)
# svg(linPlots_combined, width = 18, height = 12)
#
# cowplot::plot_grid(
# cowplot::plot_grid(lin_countP, lin_diversityP
# , nrow = 2
# , rel_heights = c(1.2,1)
# , labels = "AUTO"
# , label_size = my_label_size),
# NULL,
# linP_dm_om,
# nrow = 1,
# labels = c("", "", "C"),
# label_size = my_label_size,
# rel_widths = c(35, 3, 52)
# )
# dev.off()
# png
my_label_size = 12
linPlots_combined = paste0(outdir_images
, tolower(gene)
,"_linP_combined.png")
cat("\nOutput plot:", linPlots_combined)
png(linPlots_combined, width = 9, height = 6, units = "in" ,res = 300)
cowplot::plot_grid(
cowplot::plot_grid(lin_countP, lin_diversityP
# , nrow = 2
# , rel_heights = c(1.2,1)
# , labels = "AUTO"
# , label_size = my_label_size),
# NULL,
# linP_dm_om,
# nrow = 1,
# labels = c("", "", "C"),
# label_size = my_label_size,
# rel_widths = c(35, 3, 52)
# )
# dev.off()