345 lines
13 KiB
R
345 lines
13 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"
|
|
, "position"
|
|
, "X5uhc_position"
|
|
, "X5uhc_offset"
|
|
, "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
|
|
, scaled_cols_affinity)]
|
|
|
|
df3_plot = df3[, cols_to_extract]
|
|
##############################################################
|
|
# 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")
|
|
cols_mcsm_lig
|
|
df3_lig_ens = df3[, cols_mcsm_lig]
|
|
|
|
cols_to_numeric = c("ligand_outcome","mmcsm_lig_outcome")
|
|
df3_lig_ens[, cols_to_numeric] <- sapply(df3_lig_ens[, cols_to_numeric]
|
|
, function(x){ifelse(x == "Destabilising", 0, 1)})
|
|
|
|
cols_to_average_lig = which(colnames(df3_lig_ens)%in%cols_to_numeric)
|
|
cols_to_average_lig
|
|
|
|
# ensemble average across predictors
|
|
df3_lig_ens$ens_lig = rowMeans(df3_lig_ens[,cols_to_average_lig])
|
|
|
|
head(df3_lig_ens$position); head(df3_lig_ens$mutationinformation)
|
|
head(df3_lig_ens$ens_lig)
|
|
table(df3_lig_ens$ens_lig)
|
|
|
|
#===============================
|
|
# Filter ligand distance <10
|
|
# from globals else uncomment
|
|
#LigDist_cutoff = 10
|
|
#LigDist_colname = "ligand_distance"
|
|
#===============================
|
|
table(df3_lig_ens[LigDist_colname]<LigDist_cutoff)
|
|
expected_npos = table(df3_lig_ens$position[df3_lig_ens[LigDist_colname]<10])
|
|
expected_npos = length(expected_npos)
|
|
sum(table(df3_lig_ens$position[df3_lig_ens[LigDist_colname]<10]))
|
|
|
|
df3_lig_ens_filtered = df3_lig_ens[df3_lig_ens[LigDist_colname]<10,]
|
|
|
|
if ( nrow(df3_lig_ens_filtered) == sum(table(df3_lig_ens$position[df3_lig_ens[LigDist_colname]<10])) ){
|
|
cat(paste0("\nPASS:", LigDist_colname, "filtered according to criteria:", LigDist_cutoff, angstroms_symbol ))
|
|
}else{
|
|
stop(paste0("\nAbort:", LigDist_colname, "could not be filtered according to criteria:", LigDist_cutoff, angstroms_symbol))
|
|
}
|
|
|
|
# ensemble average of predictors by position
|
|
mean_ens_lig_by_position <- df3_lig_ens_filtered %>%
|
|
dplyr::group_by(position) %>%
|
|
dplyr::summarize(avg_ens_lig = mean(ens_lig))
|
|
|
|
class(mean_ens_lig_by_position)
|
|
|
|
# convert to a df
|
|
mean_ens_lig_by_position = as.data.frame(mean_ens_lig_by_position)
|
|
|
|
table(mean_ens_lig_by_position$avg_ens_lig)
|
|
|
|
# REscale b/w -1 and 1
|
|
#en_aff_min = min(mean_ens_affinity_by_position['ens_affinity'])
|
|
#en_aff_max = max(mean_ens_affinity_by_position['ens_affinity'])
|
|
|
|
# scale the average affintiy value between -1 and 1
|
|
# mean_ens_affinity_by_position['avg_ens_affinity_scaled'] = lapply(mean_ens_affinity_by_position['avg_ens_affinity']
|
|
# , function(x) ifelse(x < 0, x/abs(en_aff_min), x/en_aff_max))
|
|
|
|
mean_ens_lig_by_position['avg_ens_lig_scaled'] = lapply(mean_ens_lig_by_position['avg_ens_lig']
|
|
, function(x) {
|
|
scales::rescale(x, to = c(-1,1)
|
|
#, from = c(en_aff_min,en_aff_max))
|
|
, from = c(0,1))
|
|
})
|
|
|
|
cat(paste0('Average (mcsm-lig+mmcsm-lig) scores:\n'
|
|
, head(mean_ens_lig_by_position['avg_ens_lig'])
|
|
, '\n---------------------------------------------------------------'
|
|
, '\nAverage (mcsm-lig+mmcsm-lig) scaled scores:\n'
|
|
, head(mean_ens_lig_by_position['avg_ens_lig_scaled'])))
|
|
|
|
if ( nrow(mean_ens_lig_by_position) == expected_npos ){
|
|
cat("\nPASS: Generated ensemble average values for ligand affinity" )
|
|
}else{
|
|
stop(paste0("\nAbort: length mismatch for ligand affinity data"))
|
|
}
|
|
|
|
#################################################################
|
|
|
|
#=====================================
|
|
# Affintiy (mCSM-ppi2):
|
|
#D1148G for rpob DOES NOT EXIST for 5UHC
|
|
#=====================================
|
|
cols_mcsm_ppi2 = c("mutationinformation"
|
|
, "position"
|
|
, "X5uhc_position"
|
|
, "X5uhc_offset"
|
|
, "sensitivity"
|
|
, "interface_dist"
|
|
#, "mcsm_ppi2_affinity"
|
|
#, "mcsm_ppi2_scaled"
|
|
, "mcsm_ppi2_outcome"
|
|
)
|
|
cols_mcsm_ppi2
|
|
df3_ppi2_raw = df3[, c(cols_mcsm_ppi2, "mcsm_ppi2_affinity", "mcsm_ppi2_scaled") ]
|
|
|
|
table(df3_ppi2_raw$mcsm_ppi2_outcome)
|
|
|
|
df3_ppi2 = df3[, cols_mcsm_ppi2]
|
|
|
|
cols_to_numeric_ppi2 = c("mcsm_ppi2_outcome")
|
|
df3_ppi2[, cols_to_numeric_ppi2] <- sapply(df3_ppi2[, cols_to_numeric_ppi2]
|
|
, function(x){ifelse(x == "Descreasing", 0, 1)})
|
|
|
|
cols_to_average_ppi2 = which(colnames(df3_ppi2)%in%cols_to_numeric_ppi2)
|
|
cols_to_average_ppi2
|
|
|
|
#===============================
|
|
# Filter interface <10
|
|
Dist_cutoff = 10
|
|
ppi2Dist_colname = "interface_dist"
|
|
#===============================
|
|
table(df3_ppi2[ppi2Dist_colname]<Dist_cutoff)
|
|
expected_npos = table(df3_ppi2$position[df3_ppi2[ppi2Dist_colname]<10])
|
|
expected_npos = length(expected_npos)
|
|
sum(table(df3_ppi2$position[df3_ppi2[ppi2Dist_colname]<10]))
|
|
|
|
df3_ppi2_filtered = df3_ppi2[df3_ppi2[ppi2Dist_colname]<10,]
|
|
|
|
if (tolower(gene)== "rpob"){
|
|
check = nrow(df3_ppi2_filtered) == ( sum(table(df3_ppi2$position[df3_ppi2[ppi2Dist_colname]<10])) + 1)
|
|
}else{
|
|
check = nrow(df3_ppi2_filtered) == sum(table(df3_ppi2$position[df3_ppi2[ppi2Dist_colname]<10]))
|
|
}
|
|
|
|
if (check ){
|
|
cat(paste0("\nPASS:", ppi2Dist_colname
|
|
,"filtered according to criteria:"
|
|
, Dist_cutoff
|
|
, angstroms_symbol ))
|
|
}else{
|
|
stop(paste0("\nAbort:", ppi2Dist_colname
|
|
, "could not be filtered according to criteria:"
|
|
, Dist_cutoff, angstroms_symbol))
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
######################################################################
|
|
##################
|
|
# 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
|
|
#===============================================================
|