added rpob plot scripts

This commit is contained in:
Tanushree Tunstall 2022-08-29 23:26:43 +01:00
parent 8f97ab7cc8
commit 7c2e4b898e
8 changed files with 109 additions and 48 deletions

View file

@ -1,7 +1,6 @@
#!/usr/bin/env Rscript
source("~/git/LSHTM_analysis/config/gid.R")
source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
#source("~/git/LSHTM_analysis/config/gid.R")
#source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
#=======
# output
@ -13,7 +12,7 @@ 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_ppi2 = c("alr", "embb", "katg", "rpob")
geneL_both = c("rpob")
@ -21,10 +20,15 @@ if (tolower(gene)%in%geneL_na_v2) {
gene_colnames = c("mcsm_na_affinity", "mcsm_na_outcome")
}
if (tolower(gene)%in%geneL_nca) {
if (tolower(gene)%in%geneL_na) {
gene_colnames = c("mcsm_nca_affinity", "mcsm_nca_outcome")
}
if (tolower(gene)%in%geneL_both) {
gene_colnames = c("mcsm_na_affinity", "mcsm_na_outcome"
, "mcsm_ppi2_affinity", "mcsm_ppi2_outcome")
}
#from plotting_globals()
LigDist_colname
@ -36,6 +40,9 @@ angstroms_symbol
cat("\nAffinity Distance colnames:", length(affinity_dist_colnames)
, "\nThese are:", affinity_dist_colnames)
affinity_dist_colnames = c("ligand_distance","nca_distance","interface_dist") # for consistency
#===========
# Data used
#===========
@ -50,7 +57,7 @@ cols_to_output = c("position"
, "mmcsm_lig"
, "mmcsm_lig_outcome"
, affinity_dist_colnames[2]
# #, affinity_dist_colnames[3]
, affinity_dist_colnames[3]
# , "mcsm_na_affinity"
# , "mcsm_na_outcome"
# #, "mcsm_nca_affinity"
@ -148,8 +155,13 @@ colsNames_to_output_lig = c("position"
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)
#--------------------
# write output file: KS test within grpup
# write output file: lig
#----------------------
Out_ligT = paste0(outdir_stats
, tolower(gene)
@ -172,7 +184,9 @@ cols_to_output_nca = c("position"
, "sensitivity"
, "mutationinformation"
, naDist_colname
, gene_colnames
, "mcsm_na_affinity"
, "mcsm_na_outcome"
# , gene_colnames
, "maf_percent"
, "or_mychisq"
, "pval_fisher"
@ -201,8 +215,14 @@ colsNames_to_output_nca = c("position"
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)
#--------------------
# write output file: KS test within grpup
# write output file: na
#----------------------
Out_ncaT = paste0(outdir_stats
, tolower(gene)
@ -223,9 +243,9 @@ cols_to_output_ppi2 = c("position"
, "sensitivity"
, "mutationinformation"
, ppi2Dist_colname
#, "mcsm_ppi2_affinity"
#, "mcsm_ppi2_outcome"
, gene_colnames
, "mcsm_ppi2_affinity"
, "mcsm_ppi2_outcome"
#, gene_colnames
, "maf_percent"
, "or_mychisq"
, "pval_fisher"
@ -254,8 +274,14 @@ colsNames_to_output_ppi2 = c("position"
colnames(Out_df_ppi2S) = colsNames_to_output_ppi2
Out_df_ppi2S
# ADD: active site annot
nrow(Out_df_ppi2S)
Out_df_ppi2S$drug_site = ifelse(Out_df_ppi2S$position%in%aa_pos_drug, "drug", "no")
table(Out_df_ppi2S$drug_site)
#--------------------
# write output file: KS test within grpup
# write output file: ppi2
#----------------------
Out_ppi2T = paste0(outdir_stats
, tolower(gene)
@ -270,7 +296,7 @@ write.csv(Out_df_ppi2S, Out_ppi2T, row.names = FALSE)
##########################################################
# higest or/maf and stability effects
###########################################################
# convert to percet
# convert to percent
df3$maf_percent = df3$maf*100
cols_to_output_effects = c("position"
@ -282,6 +308,7 @@ cols_to_output_effects = c("position"
, "avg_lig_affinity"
, "avg_lig_affinity_outcome"
, affinity_dist_colnames[2]
, affinity_dist_colnames[3]
, gene_colnames
, "maf_percent"
, "or_mychisq"
@ -402,7 +429,7 @@ if (tolower(gene)%in%geneL_both ){
# add cols: PPI2
mut_h_ppi2_dd$mutational_effect = "Most Destabilising for PPI affinity"
mut_h_ppi2_dd$mutational_effect = "Most Stabilising 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")
@ -417,7 +444,7 @@ if (tolower(gene)%in%geneL_both ){
}
# combine BOTH: NA and PPI2
gene_aff_combined = rbind(mut_h_na_dd, mut_h_na_ss)
gene_aff_combined = rbind(mut_h_na_dd, mut_h_na_ss, mut_h_ppi2_dd, mut_h_ppi2_ss)
}
#===============
@ -447,10 +474,11 @@ colsNames_combined_table = c("position"
, "Ligand affinity outcome"
, paste0("NA-Dist (", angstroms_symbol, ")")
, paste0("PPI-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"
@ -469,8 +497,13 @@ if ( length(colnames(combined_table)) == length(colsNames_combined_table) ) {
stop("\nAbort: No. of cols mismatch. Cannot assign pretty colnames for output")
}
# ADD: active site annot
nrow(combined_table)
combined_table$drug_site = ifelse(combined_table$position%in%aa_pos_drug, "drug", "no")
table(combined_table$drug_site)
#--------------------
# write output file: KS test within grpup
# write output file: extreme effects
#----------------------
Out_combined_effectsT = paste0(outdir_stats
, tolower(gene)