LSHTM_analysis/scripts/plotting/mcsm_affinity_data_only.R

243 lines
8.5 KiB
R

#source("~/git/LSHTM_analysis/config/pnca.R")
#source("~/git/LSHTM_analysis/config/alr.R")
#source("~/git/LSHTM_analysis/config/gid.R")
#source("~/git/LSHTM_analysis/config/embb.R")
#source("~/git/LSHTM_analysis/config/katg.R")
source("~/git/LSHTM_analysis/config/rpob.R")
source("/home/tanu/git/LSHTM_analysis/my_header.R")
#########################################################
# TASK: Generate averaged affinity values
# across all affinity tools for a given structure
# as applicable...
#########################################################
#=======
# output
#=======
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
#OutFile1
outfile_mean_aff = paste0(outdir_images, "/", tolower(gene)
, "_mean_affinity_all.csv")
print(paste0("Output file:", outfile_mean_aff))
#OutFile2
outfile_mean_aff_priorty = paste0(outdir_images, "/", tolower(gene)
, "_mean_affinity_priority.csv")
print(paste0("Output file:", outfile_mean_aff_priorty))
#%%===============================================================
#=============
# Input
#=============
df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv")
df3 = read.csv(df3_filename)
length(df3$mutationinformation)
# mut_info checks
table(df3$mutation_info)
table(df3$mutation_info_orig)
table(df3$mutation_info_labels_orig)
# used in plots and analyses
table(df3$mutation_info_labels) # different, and matches dst_mode
table(df3$dst_mode)
# create column based on dst mode with different colname
table(is.na(df3$dst))
table(is.na(df3$dst_mode))
#===============
# Create column: sensitivity mapped to dst_mode
#===============
df3$sensitivity = ifelse(df3$dst_mode == 1, "R", "S")
table(df3$sensitivity)
length(unique((df3$mutationinformation)))
all_colnames = as.data.frame(colnames(df3))
common_cols = c("mutationinformation"
, "X5uhc_position"
, "X5uhc_offset"
, "position"
, "dst_mode"
, "mutation_info_labels"
, "sensitivity"
, "ligand_distance"
, "interface_dist")
all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)]
all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)]
#===================
# stability cols
#===================
# raw_cols_stability = c("duet_stability_change"
# , "deepddg"
# , "ddg_dynamut2"
# , "ddg_foldx")
#
# scaled_cols_stability = c("duet_scaled"
# , "deepddg_scaled"
# , "ddg_dynamut2_scaled"
# , "foldx_scaled")
#
# outcome_cols_stability = c("duet_outcome"
# , "deepddg_outcome"
# , "ddg_dynamut2_outcome"
# , "foldx_outcome")
#===================
# affinity cols
#===================
raw_cols_affinity = c("ligand_affinity_change"
, "mmcsm_lig"
, "mcsm_ppi2_affinity"
, "mcsm_na_affinity")
scaled_cols_affinity = c("affinity_scaled"
, "mmcsm_lig_scaled"
, "mcsm_ppi2_scaled"
, "mcsm_na_scaled" )
outcome_cols_affinity = c( "ligand_outcome"
, "mmcsm_lig_outcome"
, "mcsm_ppi2_outcome"
, "mcsm_na_outcome")
#===================
# conservation cols
#===================
# raw_cols_conservation = c("consurf_score"
# , "snap2_score"
# , "provean_score")
#
# scaled_cols_conservation = c("consurf_scaled"
# , "snap2_scaled"
# , "provean_scaled")
#
# # CANNOT strictly be used, as categories are not identical with conssurf missing altogether
# outcome_cols_conservation = c("provean_outcome"
# , "snap2_outcome"
# #consurf outcome doesn't exist
# )
######################################################################
cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols
, raw_cols_affinity
, scaled_cols_affinity
, outcome_cols_affinity
# , raw_cols_stability
# , scaled_cols_stability
# , outcome_cols_stability
)]
cols_to_extract = cols_to_consider[cols_to_consider%in%c(common_cols
, raw_cols_affinity
, scaled_cols_affinity)]
df3_plot = df3[, cols_to_extract]
DistCutOff_colnames = c("ligand_distance", "interface_dist")
DistCutOff = 10
df3_affinity_filtered = df3_plot[ (df3_plot$ligand_distance<10 | df3_plot$interface_dist <10),]
c0u = unique(df3_affinity_filtered$position)
length(c0u)
foo = df3_affinity_filtered[df3_affinity_filtered$ligand_distance<10,]
bar = df3_affinity_filtered[df3_affinity_filtered$interface_dist<10,]
wilcox.test(foo$mmcsm_lig_scaled~foo$sensitivity)
wilcox.test(foo$mmcsm_lig~foo$sensitivity)
wilcox.test(foo$affinity_scaled~foo$sensitivity)
wilcox.test(foo$ligand_affinity_change~foo$sensitivity)
wilcox.test(bar$mcsm_na_scaled~bar$sensitivity)
wilcox.test(bar$mcsm_na_affinity~bar$sensitivity)
wilcox.test(bar$mcsm_ppi2_scaled~bar$sensitivity)
wilcox.test(bar$mcsm_ppi2_affinity~bar$sensitivity)
##############################################################
df = df3_affinity_filtered
sum(is.na(df))
df2 = na.omit(df) # Apply na.omit function
a = df2[df2$position==37,]
sel_cols = c("mutationinformation", "position", scaled_cols_affinity)
a = a[, sel_cols]
##############################################################
# FIXME: ADD distance to NA when SP replies
#####################
# Ensemble affinity: affinity_cols
# mcsm_lig, mmcsm_lig and mcsm_na
#####################
# extract outcome cols and map numeric values to the categories
# Destabilising == 0, and stabilising == 1 so rescaling can let -1 be destabilising
#########################################
#=====================================
# Affintiy (2 cols): average the scores
# across predictors ==> average by
# position ==> scale b/w -1 and 1
# column to average: ens_affinity
#=====================================
cols_mcsm_lig = c("mutationinformation"
, "position"
, "sensitivity"
, "X5uhc_position"
, "X5uhc_offset"
, "ligand_distance"
, "ligand_outcome"
, "mmcsm_lig_outcome")
######################################################################
##################
# merge: mean ensemble stability and affinity by_position
####################
# if ( class(mean_ens_stability_by_position) && class(mean_ens_affinity_by_position) != "data.frame"){
# cat("Y")
# }
common_cols = intersect(colnames(mean_ens_stability_by_position), colnames(mean_ens_affinity_by_position))
if (dim(mean_ens_stability_by_position) && dim(mean_ens_affinity_by_position)){
print(paste0("PASS: dim's match, mering dfs by column :", common_cols))
#combined = as.data.frame(cbind(mean_duet_by_position, mean_affinity_by_position ))
combined_df = as.data.frame(merge(mean_ens_stability_by_position
, mean_ens_affinity_by_position
, by = common_cols
, all = T))
cat(paste0("\nnrows combined_df:", nrow(combined_df)
, "\nnrows combined_df:", ncol(combined_df)))
}else{
cat(paste0("FAIL: dim's mismatch, aborting cbind!"
, "\nnrows df1:", nrow(mean_duet_by_position)
, "\nnrows df2:", nrow(mean_affinity_by_position)))
quit()
}
#%%============================================================
# output
write.csv(combined_df, outfile_mean_ens_st_aff
, row.names = F)
cat("Finished writing file:\n"
, outfile_mean_ens_st_aff
, "\nNo. of rows:", nrow(combined_df)
, "\nNo. of cols:", ncol(combined_df))
# end of script
#===============================================================