added katg and rpob files
This commit is contained in:
parent
7c2e4b898e
commit
f39bbdcce7
9 changed files with 3145 additions and 0 deletions
215
scripts/plotting/plotting_thesis/katg/katg_ORandSNP_results.R
Normal file
215
scripts/plotting/plotting_thesis/katg/katg_ORandSNP_results.R
Normal file
|
|
@ -0,0 +1,215 @@
|
|||
#!/usr/bin/env Rscript
|
||||
#source("~/git/LSHTM_analysis/config/katg.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/")
|
||||
|
||||
###################################################################
|
||||
# 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
|
||||
###############################################
|
||||
# OR plot
|
||||
df3_or = merged_df3
|
||||
df3_or$maf_percent = df3_or$maf*100
|
||||
|
||||
bar = df3_or[, c("mutationinformation"
|
||||
, "wild_pos"
|
||||
, "position"
|
||||
, "sensitivity"
|
||||
, drug
|
||||
, affinity_dist_colnames
|
||||
, "maf_percent"
|
||||
, "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"
|
||||
, drug
|
||||
, affinity_dist_colnames
|
||||
, "maf_percent"
|
||||
, "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("\nNo.of DST muts:", nrow(bar_or)
|
||||
, "\nNo of DST (R):", table(bar_or$sensitivity)[[1]]
|
||||
, "\nNo of DST (S):", table(bar_or$sensitivity)[[2]]
|
||||
, "\nNumber of R muts with OR >1 (n = ", nrow(res1),")"
|
||||
, "\nPercentage of muts with OR>1 i.e resistant:" , pc_orR *100 )
|
||||
|
||||
table(bar_or$sensitivity)
|
||||
|
||||
# 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$heme_site = ifelse(bar_or$position%in%aa_pos_hem, "heme", "no")
|
||||
table(bar_or$heme_site)
|
||||
|
||||
top10_or = bar_or[1:10,]
|
||||
top10_or
|
||||
|
||||
write.csv(bar_or, paste0(outdir_stats, "katg_OR_10.csv"))
|
||||
|
||||
# are these active sites
|
||||
top10_or$position[top10_or$position%in%active_aa_pos]
|
||||
|
||||
|
||||
#########################################################
|
||||
# closest 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)
|
||||
460
scripts/plotting/plotting_thesis/katg/katg_appendix_tables.R
Normal file
460
scripts/plotting/plotting_thesis/katg/katg_appendix_tables.R
Normal file
|
|
@ -0,0 +1,460 @@
|
|||
#!/usr/bin/env Rscript
|
||||
#source("~/git/LSHTM_analysis/config/katg.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$heme_site = ifelse(Out_df_ligS$position%in%aa_pos_hem, "heme", "no")
|
||||
table(Out_df_ligS$heme_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/PPI2 affinity
|
||||
# naDist_colname
|
||||
# ppi2Dist_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$heme_site = ifelse(Out_df_ncaS$position%in%aa_pos_hem, "heme", "no")
|
||||
table(Out_df_ncaS$heme_site)
|
||||
|
||||
#--------------------
|
||||
# write output file: KS test within grpup
|
||||
#----------------------
|
||||
Out_ncaT = paste0(outdir_stats
|
||||
, tolower(gene)
|
||||
, "_ppi2_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: 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_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$heme_site = ifelse(combined_table$position%in%aa_pos_hem, "heme", "no")
|
||||
table(combined_table$heme_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)
|
||||
318
scripts/plotting/plotting_thesis/katg/katg_dm_om_plots.R
Normal file
318
scripts/plotting/plotting_thesis/katg/katg_dm_om_plots.R
Normal file
|
|
@ -0,0 +1,318 @@
|
|||
#################
|
||||
# Numbers
|
||||
##################
|
||||
all_dm_om_df = dm_om_wf_lf_data(df = merged_df3, gene = gene)
|
||||
#
|
||||
# lf_duet = all_dm_om_df[['lf_duet']]
|
||||
# table(lf_duet$param_type)
|
||||
################################################################
|
||||
|
||||
#======================
|
||||
# Data: Dist+Genomics
|
||||
#======================
|
||||
lf_dist_genP = all_dm_om_df[['lf_dist_gen']]
|
||||
wf_dist_genP = all_dm_om_df[['wf_dist_gen']]
|
||||
|
||||
levels(lf_dist_genP$param_type)
|
||||
#lf_dist_genP$param_type <- factor(lf_dist_genP$param_type, levels=c("Log10(MAF)", "Lig Dist(Å)", "PPI Dist(Å)"))
|
||||
table(lf_dist_genP$param_type)
|
||||
|
||||
genomics_param = c("Log10(MAF)")
|
||||
|
||||
dist_genP = lf_bp2(lf_dist_genP
|
||||
#, p_title
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
#dist_genP
|
||||
#-------------------
|
||||
# Genomics data plot
|
||||
#-------------------
|
||||
genomics_dataP = lf_dist_genP[lf_dist_genP$param_type%in%genomics_param,]
|
||||
genomics_dataP$param_type = factor(genomics_dataP$param_type)
|
||||
table(genomics_dataP$param_type)
|
||||
|
||||
genomicsP = lf_bp2(genomics_dataP
|
||||
#, p_title = ""
|
||||
, dot_transparency = 0.05 #0.3 default
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
genomicsP
|
||||
|
||||
# #check
|
||||
# wilcox.test(wf_dist_genP$`Log10(MAF)`[wf_dist_genP$mutation_info_labels=="R"]
|
||||
# , wf_dist_genP$`Log10(MAF)`[wf_dist_genP$mutation_info_labels=="S"], paired = FALSE)
|
||||
#
|
||||
# tapply(wf_dist_genP$`Log10(MAF)`, wf_dist_genP$mutation_info_labels, summary)
|
||||
|
||||
#-------------------
|
||||
# Distance data plot:
|
||||
#--------------------
|
||||
# not genomics
|
||||
dist_dataP = lf_dist_genP[!lf_dist_genP$param_type%in%genomics_param,]
|
||||
dist_dataP$param_type = factor(dist_dataP$param_type)
|
||||
table(dist_dataP$param_type)
|
||||
levels(dist_dataP$param_type)
|
||||
# relevel factor to control ordering of appearance of plot
|
||||
dist_dataP$param_type <-relevel(dist_dataP$param_type, "Lig Dist(Å)" )
|
||||
table(dist_dataP$param_type)
|
||||
levels(dist_dataP$param_type)
|
||||
|
||||
distanceP = lf_bp2(dist_dataP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
distanceP
|
||||
|
||||
# # check
|
||||
# wilcox.test(wf_dist_genP$`PPI Dist(Å)`[wf_dist_genP$mutation_info_labels=="R"]
|
||||
# , wf_dist_genP$`PPI Dist(Å)`[wf_dist_genP$mutation_info_labels=="S"], paired = FALSE)
|
||||
#
|
||||
# wilcox.test(wf_dist_genP$`Lig Dist(Å)`[wf_dist_genP$mutation_info_labels=="R"]
|
||||
# , wf_dist_genP$`Lig Dist(Å)`[wf_dist_genP$mutation_info_labels=="S"], paired = FALSE)
|
||||
#
|
||||
# tapply(wf_dist_genP$`PPI Dist(Å)`, wf_dist_genP$mutation_info_labels, summary)
|
||||
#
|
||||
# tapply(wf_dist_genP$`Lig Dist(Å)`, wf_dist_genP$mutation_info_labels, summary)
|
||||
|
||||
|
||||
#-------------------
|
||||
# Distance data plot: LigDist
|
||||
#--------------------
|
||||
levels(dist_dataP$param_type)[[1]]
|
||||
#Lig Dist(Å), PPI Dist(Å)
|
||||
dist_data_lig = dist_dataP[dist_dataP$param_type%in%c(levels(dist_dataP$param_type)[[1]]),]
|
||||
dist_data_lig$param_type = factor(dist_data_lig$param_type)
|
||||
table(dist_data_lig$param_type)
|
||||
levels(dist_data_lig$param_type)
|
||||
distanceP_lig = lf_bp2(dist_data_lig
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
distanceP_lig
|
||||
|
||||
if (tolower(gene)%in%geneL_ppi2){
|
||||
#-------------------
|
||||
# Distance data plot: LigDist
|
||||
#--------------------
|
||||
levels(dist_dataP$param_type)[[2]]
|
||||
#Lig Dist(Å), PPI Dist(Å)
|
||||
dist_data_ppi2 = dist_dataP[dist_dataP$param_type%in%c(levels(dist_dataP$param_type)[[2]]),]
|
||||
dist_data_ppi2$param_type = factor(dist_data_ppi2$param_type)
|
||||
table(dist_data_ppi2$param_type)
|
||||
levels(dist_data_ppi2$param_type)
|
||||
distanceP_ppi2 = lf_bp2(dist_data_ppi2
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5)
|
||||
, dot_transparency = 0.2
|
||||
, monochrome = F)
|
||||
|
||||
distanceP_ppi2
|
||||
}
|
||||
|
||||
if (tolower(gene)%in%geneL_na){
|
||||
#-------------------
|
||||
# Distance data plot: NADist
|
||||
#--------------------
|
||||
levels(dist_dataP$param_type)[[2]]
|
||||
#Lig Dist(Å), PPI Dist(Å)
|
||||
dist_data_na = dist_dataP[dist_dataP$param_type%in%c(levels(dist_dataP$param_type)[[2]]),]
|
||||
dist_data_na$param_type = factor(dist_data_na$param_type)
|
||||
table(dist_data_na$param_type)
|
||||
levels(dist_data_na$param_type)
|
||||
distanceP_na = lf_bp2(dist_data_na
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
distanceP_na
|
||||
}
|
||||
#==============
|
||||
# Plot:DUET
|
||||
#==============
|
||||
lf_duetP = all_dm_om_df[['lf_duet']]
|
||||
#lf_duetP = lf_duet[!lf_duet$param_type%in%c(static_colsP),]
|
||||
table(lf_duetP$param_type)
|
||||
lf_duetP$param_type = factor(lf_duetP$param_type)
|
||||
table(lf_duetP$param_type)
|
||||
|
||||
duetP = lf_bp2(lf_duetP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.2)
|
||||
|
||||
#==============
|
||||
# Plot:FoldX
|
||||
#==============
|
||||
lf_foldxP = all_dm_om_df[['lf_foldx']]
|
||||
#lf_foldxP = lf_foldx[!lf_foldx$param_type%in%c(static_colsP),]
|
||||
table(lf_foldxP$param_type)
|
||||
lf_foldxP$param_type = factor(lf_foldxP$param_type)
|
||||
table(lf_foldxP$param_type)
|
||||
|
||||
foldxP = lf_bp2(lf_foldxP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.1)
|
||||
|
||||
#==============
|
||||
# Plot:DeepDDG
|
||||
#==============
|
||||
lf_deepddgP = all_dm_om_df[['lf_deepddg']]
|
||||
#lf_deepddgP = lf_deepddg[!lf_deepddg$param_type%in%c(static_colsP),]
|
||||
table(lf_deepddgP$param_type)
|
||||
lf_deepddgP$param_type = factor(lf_deepddgP$param_type)
|
||||
table(lf_deepddgP$param_type)
|
||||
|
||||
deepddgP = lf_bp2(lf_deepddgP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.2)
|
||||
|
||||
deepddgP
|
||||
|
||||
#==============
|
||||
# Plot: Dynamut2
|
||||
#==============
|
||||
lf_dynamut2P = all_dm_om_df[['lf_dynamut2']]
|
||||
#lf_dynamut2P = lf_dynamut2[!lf_dynamut2$param_type%in%c(static_colsP),]
|
||||
table(lf_dynamut2P$param_type)
|
||||
lf_dynamut2P$param_type = factor(lf_dynamut2P$param_type)
|
||||
table(lf_dynamut2P$param_type)
|
||||
|
||||
dynamut2P = lf_bp2(lf_dynamut2P
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.2)
|
||||
|
||||
|
||||
#==============
|
||||
# Plot:ConSurf
|
||||
#==============
|
||||
lf_consurfP = all_dm_om_df[['lf_consurf']]
|
||||
#lf_consurfP = lf_consurf[!lf_consurf$param_type%in%c(static_colsP),]
|
||||
table(lf_consurfP$param_type)
|
||||
lf_consurfP$param_type = factor(lf_consurfP$param_type)
|
||||
table(lf_consurfP$param_type)
|
||||
|
||||
consurfP = lf_bp2(lf_consurfP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
#==============
|
||||
# Plot:PROVEAN
|
||||
#==============
|
||||
lf_proveanP = all_dm_om_df[['lf_provean']]
|
||||
#lf_proveanP = lf_provean[!lf_provean$param_type%in%c(static_colsP),]
|
||||
table(lf_proveanP$param_type)
|
||||
lf_proveanP$param_type = factor(lf_proveanP$param_type)
|
||||
table(lf_proveanP$param_type)
|
||||
|
||||
proveanP = lf_bp2(lf_proveanP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
#==============
|
||||
# Plot:SNAP2
|
||||
#==============
|
||||
lf_snap2P = all_dm_om_df[['lf_snap2']]
|
||||
#lf_snap2P = lf_snap2[!lf_snap2$param_type%in%c(static_colsP),]
|
||||
table(lf_snap2P$param_type)
|
||||
lf_snap2P$param_type = factor(lf_snap2P$param_type)
|
||||
table(lf_snap2P$param_type)
|
||||
|
||||
snap2P = lf_bp2(lf_snap2P
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
|
||||
############################################################################
|
||||
#================
|
||||
# Plot: mCSM-lig
|
||||
#================
|
||||
lf_mcsm_ligP = all_dm_om_df[['lf_mcsm_lig']]
|
||||
#lf_mcsm_ligP = lf_mcsm_lig[!lf_mcsm_lig$param_type%in%c(static_colsP),]
|
||||
table(lf_mcsm_ligP$param_type)
|
||||
lf_mcsm_ligP$param_type = factor(lf_mcsm_ligP$param_type)
|
||||
table(lf_mcsm_ligP$param_type)
|
||||
|
||||
mcsmligP = lf_bp2(lf_mcsm_ligP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.8)
|
||||
|
||||
mcsmligP
|
||||
#=================
|
||||
# Plot: mmCSM-lig2
|
||||
#=================
|
||||
lf_mmcsm_lig2P = all_dm_om_df[['lf_mmcsm_lig2']]
|
||||
#lf_mmcsm_lig2P = lf_mmcsm_lig2P[!lf_mmcsm_lig2P$param_type%in%c(static_colsP),]
|
||||
table(lf_mmcsm_lig2P$param_type)
|
||||
lf_mmcsm_lig2P$param_type = factor(lf_mmcsm_lig2P$param_type)
|
||||
table(lf_mmcsm_lig2P$param_type)
|
||||
|
||||
mcsmlig2P = lf_bp2(lf_mmcsm_lig2P
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.8)
|
||||
|
||||
mcsmlig2P
|
||||
|
||||
#================
|
||||
# Plot: mCSM-ppi2
|
||||
#================
|
||||
if (tolower(gene)%in%geneL_ppi2){
|
||||
lf_mcsm_ppi2P = all_dm_om_df[['lf_mcsm_ppi2']]
|
||||
#lf_mcsm_ppi2P = lf_mcsm_ppi2[!lf_mcsm_ppi2$param_type%in%c(static_colsP),]
|
||||
table(lf_mcsm_ppi2P$param_type)
|
||||
lf_mcsm_ppi2P$param_type = factor(lf_mcsm_ppi2P$param_type)
|
||||
table(lf_mcsm_ppi2P$param_type)
|
||||
|
||||
mcsmppi2P = lf_bp2(lf_mcsm_ppi2P
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.3)
|
||||
|
||||
}
|
||||
|
||||
#==============
|
||||
# Plot: mCSM-NA
|
||||
#==============
|
||||
if (tolower(gene)%in%geneL_na){
|
||||
lf_mcsm_naP = all_dm_om_df[['lf_mcsm_na']]
|
||||
#lf_mcsm_naP = lf_mcsm_na[!lf_mcsm_na$param_type%in%c(static_colsP),]
|
||||
table(lf_mcsm_naP$param_type)
|
||||
lf_mcsm_naP$param_type = factor(lf_mcsm_naP$param_type)
|
||||
table(lf_mcsm_naP$param_type)
|
||||
|
||||
mcsmnaP = lf_bp2(lf_mcsm_naP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.4)
|
||||
|
||||
}
|
||||
|
||||
######################################
|
||||
# Outplot with stats
|
||||
######################################
|
||||
# outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
|
||||
#
|
||||
# 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)
|
||||
# cowplot::plot_grid(
|
||||
# cowplot::plot_grid(duetP, foldxP, deepddgP, dynamut2P, genomicsP, distanceP
|
||||
# , nrow=1
|
||||
# , rel_widths = c(1/7, 1/7,1/7,1/7, 1/7, 1.75/7)),
|
||||
# #, rel_widths = c(1/8, 1/8,1/8,1/8, 1/8, 2.75/8)), # for 3 distances
|
||||
# cowplot::plot_grid(consurfP, proveanP, snap2P
|
||||
# , mcsmligP
|
||||
# , mcsmlig2P
|
||||
# , mcsmppi2P
|
||||
# #, mcsmnaP
|
||||
# , nrow=1),
|
||||
# nrow=2)
|
||||
#
|
||||
# dev.off()
|
||||
|
||||
|
||||
176
scripts/plotting/plotting_thesis/katg/katg_dm_om_plots_layout.R
Normal file
176
scripts/plotting/plotting_thesis/katg/katg_dm_om_plots_layout.R
Normal file
|
|
@ -0,0 +1,176 @@
|
|||
# source dm_om_plots.R
|
||||
source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/katg/katg_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()
|
||||
545
scripts/plotting/plotting_thesis/katg/katg_ks_test_all_PS.R
Executable file
545
scripts/plotting/plotting_thesis/katg/katg_ks_test_all_PS.R
Executable 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/katg.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)
|
||||
|
|
@ -0,0 +1,460 @@
|
|||
#!/usr/bin/env Rscript
|
||||
#source("~/git/LSHTM_analysis/config/katg.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$heme_site = ifelse(Out_df_ligS$position%in%aa_pos_hem, "heme", "no")
|
||||
table(Out_df_ligS$heme_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/PPI2 affinity
|
||||
# naDist_colname
|
||||
# ppi2Dist_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$heme_site = ifelse(Out_df_ncaS$position%in%aa_pos_hem, "heme", "no")
|
||||
table(Out_df_ncaS$heme_site)
|
||||
|
||||
#--------------------
|
||||
# write output file: KS test within grpup
|
||||
#----------------------
|
||||
Out_ncaT = paste0(outdir_stats
|
||||
, tolower(gene)
|
||||
, "_ppi2_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: 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_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$heme_site = ifelse(combined_table$position%in%aa_pos_hem, "heme", "no")
|
||||
table(combined_table$heme_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)
|
||||
441
scripts/plotting/plotting_thesis/rpob/TEMp_na_appendix_tables.R
Normal file
441
scripts/plotting/plotting_thesis/rpob/TEMp_na_appendix_tables.R
Normal file
|
|
@ -0,0 +1,441 @@
|
|||
#!/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), "/")
|
||||
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_nca = 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_nca) {
|
||||
gene_colnames = c("mcsm_nca_affinity", "mcsm_nca_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)
|
||||
|
||||
#--------------------
|
||||
# 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
|
||||
# naDist_colname
|
||||
####################################
|
||||
# Filtered data
|
||||
df_nca = df3_output[df3_output[[naDist_colname]]<DistCutOff,]
|
||||
|
||||
# select cols
|
||||
cols_to_output_nca = c("position"
|
||||
, "sensitivity"
|
||||
, "mutationinformation"
|
||||
, naDist_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("MAF ","(%)")
|
||||
, "Odds Ratio"
|
||||
, "P-value"
|
||||
, "Adjusted P-value"
|
||||
, "Adjusted P-value significance")
|
||||
|
||||
colnames(Out_df_ncaS) = colsNames_to_output_nca
|
||||
Out_df_ncaS
|
||||
#--------------------
|
||||
# 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: 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")
|
||||
}
|
||||
|
||||
|
||||
nrow(combined_table)
|
||||
combined_table$drug_site = ifelse(combined_table$position%in%aa_pos_drug, "drug", "no")
|
||||
table(combined_table$drug_site)
|
||||
|
||||
combined_table$rna_site = ifelse(combined_table$position%in%aa_pos_rna, "rna", "no")
|
||||
table(combined_table$rna_site)
|
||||
|
||||
combined_table$sam_site = ifelse(combined_table$position%in%aa_pos_sam, "sam", "no")
|
||||
table(combined_table$sam_site)
|
||||
|
||||
combined_table$amp_site = ifelse(combined_table$position%in%aa_pos_amp, "amp", "no")
|
||||
table(combined_table$amp_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)
|
||||
304
scripts/plotting/plotting_thesis/rpob/dm_om_plots_rpob.R
Normal file
304
scripts/plotting/plotting_thesis/rpob/dm_om_plots_rpob.R
Normal file
|
|
@ -0,0 +1,304 @@
|
|||
#################
|
||||
# Numbers
|
||||
##################
|
||||
all_dm_om_df = dm_om_wf_lf_data(df = merged_df3, gene = gene)
|
||||
#
|
||||
# lf_duet = all_dm_om_df[['lf_duet']]
|
||||
# table(lf_duet$param_type)
|
||||
################################################################
|
||||
|
||||
#======================
|
||||
# Data: Dist+Genomics
|
||||
#======================
|
||||
lf_dist_genP = all_dm_om_df[['lf_dist_gen']]
|
||||
wf_dist_genP = all_dm_om_df[['wf_dist_gen']]
|
||||
|
||||
levels(lf_dist_genP$param_type)
|
||||
#lf_dist_genP$param_type <- factor(lf_dist_genP$param_type, levels=c("Log10(MAF)", "Lig Dist(Å)", "PPI Dist(Å)"))
|
||||
table(lf_dist_genP$param_type)
|
||||
|
||||
genomics_param = c("Log10(MAF)")
|
||||
|
||||
dist_genP = lf_bp2(lf_dist_genP
|
||||
#, p_title
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
#dist_genP
|
||||
#-------------------
|
||||
# Genomics data plot
|
||||
#-------------------
|
||||
genomics_dataP = lf_dist_genP[lf_dist_genP$param_type%in%genomics_param,]
|
||||
genomics_dataP$param_type = factor(genomics_dataP$param_type)
|
||||
table(genomics_dataP$param_type)
|
||||
|
||||
genomicsP = lf_bp2(genomics_dataP
|
||||
#, p_title = ""
|
||||
, dot_transparency = 0.1 #0.3 default
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
genomicsP
|
||||
|
||||
#-------------------
|
||||
# Distance data plot:
|
||||
#--------------------
|
||||
# not genomics
|
||||
dist_dataP = lf_dist_genP[!lf_dist_genP$param_type%in%genomics_param,]
|
||||
dist_dataP$param_type = factor(dist_dataP$param_type)
|
||||
table(dist_dataP$param_type)
|
||||
levels(dist_dataP$param_type)
|
||||
# relevel factor to control ordering of appearance of plot
|
||||
dist_dataP$param_type <-relevel(dist_dataP$param_type, "Lig Dist(Å)" )
|
||||
table(dist_dataP$param_type)
|
||||
levels(dist_dataP$param_type)
|
||||
|
||||
distanceP = lf_bp2(dist_dataP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.4)
|
||||
|
||||
distanceP
|
||||
|
||||
#-----------------------
|
||||
# Distance data plot: LigDist
|
||||
#-----------------------
|
||||
levels(dist_dataP$param_type)[[1]]
|
||||
#Lig Dist(Å), PPI Dist(Å)
|
||||
dist_data_lig = dist_dataP[dist_dataP$param_type%in%c(levels(dist_dataP$param_type)[[1]]),]
|
||||
dist_data_lig$param_type = factor(dist_data_lig$param_type)
|
||||
table(dist_data_lig$param_type)
|
||||
levels(dist_data_lig$param_type)
|
||||
distanceP_lig = lf_bp2(dist_data_lig
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.4)
|
||||
|
||||
distanceP_lig
|
||||
|
||||
if (tolower(gene)%in%geneL_ppi2){
|
||||
#-------------------
|
||||
# Distance data plot: PPI2 dist
|
||||
#--------------------
|
||||
levels(dist_dataP$param_type)[[2]]
|
||||
#Lig Dist(Å), PPI Dist(Å)
|
||||
dist_data_ppi2 = dist_dataP[dist_dataP$param_type%in%c(levels(dist_dataP$param_type)[[2]]),]
|
||||
dist_data_ppi2$param_type = factor(dist_data_ppi2$param_type)
|
||||
table(dist_data_ppi2$param_type)
|
||||
levels(dist_data_ppi2$param_type)
|
||||
distanceP_ppi2 = lf_bp2(dist_data_ppi2
|
||||
#, p_title = "ppi"
|
||||
, violin_quantiles = c(0.5)
|
||||
, dot_transparency = 0.4
|
||||
, monochrome = F)
|
||||
|
||||
distanceP_ppi2
|
||||
}
|
||||
|
||||
if (tolower(gene)%in%geneL_na){
|
||||
#-------------------
|
||||
# Distance data plot: NADist
|
||||
#--------------------
|
||||
levels(dist_dataP$param_type)[[3]]
|
||||
#Lig Dist(Å), PPI Dist(Å)
|
||||
dist_data_na = dist_dataP[dist_dataP$param_type%in%c(levels(dist_dataP$param_type)[[3]]),]
|
||||
dist_data_na$param_type = factor(dist_data_na$param_type)
|
||||
table(dist_data_na$param_type)
|
||||
levels(dist_data_na$param_type)
|
||||
distanceP_na = lf_bp2(dist_data_na
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.4)
|
||||
|
||||
distanceP_na
|
||||
}
|
||||
|
||||
#==============
|
||||
# Plot:DUET
|
||||
#==============
|
||||
lf_duetP = all_dm_om_df[['lf_duet']]
|
||||
#lf_duetP = lf_duet[!lf_duet$param_type%in%c(static_colsP),]
|
||||
table(lf_duetP$param_type)
|
||||
lf_duetP$param_type = factor(lf_duetP$param_type)
|
||||
table(lf_duetP$param_type)
|
||||
|
||||
duetP = lf_bp2(lf_duetP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.08)
|
||||
|
||||
#==============
|
||||
# Plot:FoldX
|
||||
#==============
|
||||
lf_foldxP = all_dm_om_df[['lf_foldx']]
|
||||
#lf_foldxP = lf_foldx[!lf_foldx$param_type%in%c(static_colsP),]
|
||||
table(lf_foldxP$param_type)
|
||||
lf_foldxP$param_type = factor(lf_foldxP$param_type)
|
||||
table(lf_foldxP$param_type)
|
||||
|
||||
foldxP = lf_bp2(lf_foldxP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.05)
|
||||
|
||||
#==============
|
||||
# Plot:DeepDDG
|
||||
#==============
|
||||
lf_deepddgP = all_dm_om_df[['lf_deepddg']]
|
||||
#lf_deepddgP = lf_deepddg[!lf_deepddg$param_type%in%c(static_colsP),]
|
||||
table(lf_deepddgP$param_type)
|
||||
lf_deepddgP$param_type = factor(lf_deepddgP$param_type)
|
||||
table(lf_deepddgP$param_type)
|
||||
|
||||
deepddgP = lf_bp2(lf_deepddgP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.05)
|
||||
|
||||
deepddgP
|
||||
|
||||
#==============
|
||||
# Plot: Dynamut2
|
||||
#==============
|
||||
lf_dynamut2P = all_dm_om_df[['lf_dynamut2']]
|
||||
#lf_dynamut2P = lf_dynamut2[!lf_dynamut2$param_type%in%c(static_colsP),]
|
||||
table(lf_dynamut2P$param_type)
|
||||
lf_dynamut2P$param_type = factor(lf_dynamut2P$param_type)
|
||||
table(lf_dynamut2P$param_type)
|
||||
|
||||
dynamut2P = lf_bp2(lf_dynamut2P
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.08)
|
||||
|
||||
|
||||
#==============
|
||||
# Plot:ConSurf
|
||||
#==============
|
||||
lf_consurfP = all_dm_om_df[['lf_consurf']]
|
||||
#lf_consurfP = lf_consurf[!lf_consurf$param_type%in%c(static_colsP),]
|
||||
table(lf_consurfP$param_type)
|
||||
lf_consurfP$param_type = factor(lf_consurfP$param_type)
|
||||
table(lf_consurfP$param_type)
|
||||
|
||||
consurfP = lf_bp2(lf_consurfP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
#==============
|
||||
# Plot:PROVEAN
|
||||
#==============
|
||||
lf_proveanP = all_dm_om_df[['lf_provean']]
|
||||
#lf_proveanP = lf_provean[!lf_provean$param_type%in%c(static_colsP),]
|
||||
table(lf_proveanP$param_type)
|
||||
lf_proveanP$param_type = factor(lf_proveanP$param_type)
|
||||
table(lf_proveanP$param_type)
|
||||
|
||||
proveanP = lf_bp2(lf_proveanP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
#==============
|
||||
# Plot:SNAP2
|
||||
#==============
|
||||
lf_snap2P = all_dm_om_df[['lf_snap2']]
|
||||
#lf_snap2P = lf_snap2[!lf_snap2$param_type%in%c(static_colsP),]
|
||||
table(lf_snap2P$param_type)
|
||||
lf_snap2P$param_type = factor(lf_snap2P$param_type)
|
||||
table(lf_snap2P$param_type)
|
||||
|
||||
snap2P = lf_bp2(lf_snap2P
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F)
|
||||
|
||||
|
||||
############################################################################
|
||||
#================
|
||||
# Plot: mCSM-lig
|
||||
#================
|
||||
lf_mcsm_ligP = all_dm_om_df[['lf_mcsm_lig']]
|
||||
#lf_mcsm_ligP = lf_mcsm_lig[!lf_mcsm_lig$param_type%in%c(static_colsP),]
|
||||
table(lf_mcsm_ligP$param_type)
|
||||
lf_mcsm_ligP$param_type = factor(lf_mcsm_ligP$param_type)
|
||||
table(lf_mcsm_ligP$param_type)
|
||||
|
||||
mcsmligP = lf_bp2(lf_mcsm_ligP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.8)
|
||||
|
||||
mcsmligP
|
||||
#=================
|
||||
# Plot: mmCSM-lig2
|
||||
#=================
|
||||
lf_mmcsm_lig2P = all_dm_om_df[['lf_mmcsm_lig2']]
|
||||
#lf_mmcsm_lig2P = lf_mmcsm_lig2P[!lf_mmcsm_lig2P$param_type%in%c(static_colsP),]
|
||||
table(lf_mmcsm_lig2P$param_type)
|
||||
lf_mmcsm_lig2P$param_type = factor(lf_mmcsm_lig2P$param_type)
|
||||
table(lf_mmcsm_lig2P$param_type)
|
||||
|
||||
mcsmlig2P = lf_bp2(lf_mmcsm_lig2P
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.8)
|
||||
|
||||
mcsmlig2P
|
||||
|
||||
#================
|
||||
# Plot: mCSM-ppi2
|
||||
#================
|
||||
if (tolower(gene)%in%geneL_ppi2){
|
||||
lf_mcsm_ppi2P = all_dm_om_df[['lf_mcsm_ppi2']]
|
||||
#lf_mcsm_ppi2P = lf_mcsm_ppi2[!lf_mcsm_ppi2$param_type%in%c(static_colsP),]
|
||||
table(lf_mcsm_ppi2P$param_type)
|
||||
lf_mcsm_ppi2P$param_type = factor(lf_mcsm_ppi2P$param_type)
|
||||
table(lf_mcsm_ppi2P$param_type)
|
||||
|
||||
mcsmppi2P = lf_bp2(lf_mcsm_ppi2P
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.3)
|
||||
|
||||
}
|
||||
|
||||
#==============
|
||||
# Plot: mCSM-NA
|
||||
#==============
|
||||
if (tolower(gene)%in%geneL_na){
|
||||
lf_mcsm_naP = all_dm_om_df[['lf_mcsm_na']]
|
||||
#lf_mcsm_naP = lf_mcsm_na[!lf_mcsm_na$param_type%in%c(static_colsP),]
|
||||
table(lf_mcsm_naP$param_type)
|
||||
lf_mcsm_naP$param_type = factor(lf_mcsm_naP$param_type)
|
||||
table(lf_mcsm_naP$param_type)
|
||||
|
||||
mcsmnaP = lf_bp2(lf_mcsm_naP
|
||||
#, p_title = ""
|
||||
, violin_quantiles = c(0.5), monochrome = F
|
||||
, dot_transparency = 0.4)
|
||||
|
||||
}
|
||||
|
||||
######################################
|
||||
# Outplot with stats
|
||||
######################################
|
||||
# outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
|
||||
#
|
||||
# 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)
|
||||
# cowplot::plot_grid(
|
||||
# cowplot::plot_grid(duetP, foldxP, deepddgP, dynamut2P, genomicsP, distanceP
|
||||
# , nrow=1
|
||||
# , rel_widths = c(1/7, 1/7,1/7,1/7, 1/7, 1.75/7)),
|
||||
# #, rel_widths = c(1/8, 1/8,1/8,1/8, 1/8, 2.75/8)), # for 3 distances
|
||||
# cowplot::plot_grid(consurfP, proveanP, snap2P
|
||||
# , mcsmligP
|
||||
# , mcsmlig2P
|
||||
# , mcsmppi2P
|
||||
# #, mcsmnaP
|
||||
# , nrow=1),
|
||||
# nrow=2)
|
||||
#
|
||||
# dev.off()
|
||||
|
||||
|
||||
226
scripts/plotting/plotting_thesis/rpob/rpob_ORandSNP_results.R
Normal file
226
scripts/plotting/plotting_thesis/rpob/rpob_ORandSNP_results.R
Normal file
|
|
@ -0,0 +1,226 @@
|
|||
#!/usr/bin/env Rscript
|
||||
#source("~/git/LSHTM_analysis/config/katg.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/")
|
||||
|
||||
###################################################################
|
||||
# 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
|
||||
|
||||
dist_columns = c(LigDist_colname, ppi2Dist_colname, naDist_colname )
|
||||
|
||||
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
|
||||
all_cols = c(common_cols
|
||||
, all_stability_cols
|
||||
, all_affinity_cols
|
||||
, all_conserv_cols)
|
||||
|
||||
plot_cols = c("mutationinformation", "mutation_info_labels", "position", "dst_mode"
|
||||
#, all_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
|
||||
###############################################
|
||||
# OR plot
|
||||
df3_or = merged_df3
|
||||
df3_or$maf_percent = df3_or$maf*100
|
||||
|
||||
bar = df3_or[, c("mutationinformation"
|
||||
, "wild_pos"
|
||||
, "position"
|
||||
, "sensitivity"
|
||||
, drug
|
||||
, affinity_dist_colnames
|
||||
, "maf_percent"
|
||||
, "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"
|
||||
, drug
|
||||
, affinity_dist_colnames
|
||||
, "maf_percent"
|
||||
, "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("\nNo.of DST muts:", nrow(bar_or)
|
||||
, "\nNo of DST (R):", table(bar_or$sensitivity)[[1]]
|
||||
, "\nNo of DST (S):", table(bar_or$sensitivity)[[2]]
|
||||
, "\nNumber of R muts with OR >1 (n = ", nrow(res1),")"
|
||||
, "\nPercentage of muts with OR>1 i.e resistant:" , pc_orR *100 )
|
||||
|
||||
table(bar_or$sensitivity)
|
||||
|
||||
# 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_distance
|
||||
, 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)
|
||||
|
||||
top10_or = bar_or[1:10,]
|
||||
top10_or
|
||||
|
||||
write.csv(bar_or, paste0(outdir_stats, "rpob_OR_10.csv"))
|
||||
|
||||
# are these active sites
|
||||
top10_or$position[top10_or$position%in%active_aa_pos]
|
||||
|
||||
|
||||
# most frequent
|
||||
bar_maf = bar_or[order(bar_or$maf_percent
|
||||
, bar_or$ligand_distance
|
||||
, bar_or$nca_distance
|
||||
, bar_or$interface_dist
|
||||
, decreasing = T), ]
|
||||
bar_maf$drug_site = ifelse(bar_maf$position%in%aa_pos_drug, "drug", "no")
|
||||
table(bar_maf$drug_site)
|
||||
bar_maf[1:10,]
|
||||
#########################################################
|
||||
# closest 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_nca = bar_or[bar_or$nca_distance<10,]
|
||||
bar_or_nca = bar_or_nca[order(bar_or_nca$nca_distance, -bar_or_ppi$or_mychisq), ]
|
||||
table(bar_or_nca$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)
|
||||
Loading…
Add table
Add a link
Reference in a new issue