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

attempting affintiy stuff

Browse source
This commit is contained in:
Tanushree Tunstall 2022-08-01 21:41:02 +01:00
parent 0d8979dfcb
commit ccc877e811
2 changed files with 151 additions and 47 deletions
Download patch file Download diff file Expand all files Collapse all files

189
scripts/plotting/mcsm_mean_affinity_ensemble.R
View file

@ -59,10 +59,14 @@ 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")
, "ligand_distance"
, "interface_dist")
all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)]
all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)]
@ -122,25 +126,28 @@ outcome_cols_affinity = c( "ligand_outcome"
######################################################################
cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols
, raw_cols
, scaled_cols
, outcome_cols_affinity)]
, 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
# , outcome_cols_affinity)]
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
# df3_plot = df3[, cols_to_extract]
#
# df3_plot[, outcome_cols_affinity] <- sapply(df3_plot[, outcome_cols_affinity]
# , function(x){ifelse(x == "Destabilising", 0, 1)})
df3_plot = df3[, c(common_cols, scaled_cols)]
#########################################
#=====================================
# Affintiy (2 cols): average the scores
# across predictors ==> average by
@ -148,20 +155,61 @@ df3_plot = df3[, c(common_cols, scaled_cols)]
# column to average: ens_affinity
#=====================================
cols_to_average_affinity = which(colnames(df3_plot)%in%outcome_cols_affinity)
cols_to_average_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_plot_affinity$ens_affinity = rowMeans(df3_plot_affinity[,cols_to_average_affinity])
df3_lig_ens$ens_lig = rowMeans(df3_lig_ens[,cols_to_average_lig])
head(df3_plot_affinity$position); head(df3_plot_affinity$mutationinformation)
head(df3_plot_affinity$ens_affinity)
table(df3_plot_affinity$ens_affinity)
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_affinity_by_position <- df3_plot_affinity %>%
mean_ens_lig_by_position <- df3_lig_ens_filtered %>%
dplyr::group_by(position) %>%
dplyr::summarize(avg_ens_affinity = mean(ens_affinity))
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'])
@ -171,38 +219,91 @@ mean_ens_affinity_by_position <- df3_plot_affinity %>%
# 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_affinity_by_position['avg_ens_affinity_scaled'] = lapply(mean_ens_affinity_by_position['avg_ens_affinity']
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 affintiy scores:\n'
, head(mean_ens_affinity_by_position['avg_ens_affinity'])
cat(paste0('Average (mcsm-lig+mmcsm-lig) scores:\n'
, head(mean_ens_lig_by_position['avg_ens_lig'])
, '\n---------------------------------------------------------------'
, '\nAverage affintiy scaled scores:\n'
, head(mean_ens_affinity_by_position['avg_ens_affinity_scaled'])))
, '\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))
}
#convert to a df
mean_ens_affinity_by_position = as.data.frame(mean_ens_affinity_by_position)
#FIXME: sanity checks
# TODO: predetermine the bounds
# l_bound_ens_aff = min(mean_ens_affintiy_by_position['avg_ens_affinity_scaled'])
# u_bound_ens_aff = max(mean_ens_affintiy_by_position['avg_ens_affinity_scaled'])
#
# if ( (l_bound_ens_aff == -1) && (u_bound_ens_aff == 1) ){
# cat(paste0("PASS: ensemble affinity scores averaged by position and then scaled"
# , "\nmin ensemble averaged affinity: ", l_bound_ens_aff
# , "\nmax ensemble averaged affinity: ", u_bound_ens_aff))
# }else{
# cat(paste0("FAIL: ensemble affinity scores could not be scaled b/w -1 and 1"
# , "\nmin ensemble averaged affinity: ", l_bound_ens_aff
# , "\nmax ensemble averaged affinity: ", u_bound_ens_aff))
# quit()
# }
######################################################################

9
scripts/plotting/mcsm_mean_stability_ensemble.R
View file

@ -55,7 +55,8 @@ common_cols = c("mutationinformation"
, "dst_mode"
, "mutation_info_labels"
, "sensitivity"
, "ligand_distance")
, "ligand_distance"
, "interface_dist")
all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)]
all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)]
@ -117,11 +118,13 @@ outcome_cols_affinity = c( "ligand_outcome"
cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols
, raw_cols_stability
, scaled_cols_stability
, outcome_cols_stability)]
, outcome_cols_stability
, raw_cols_affinity
, scaled_cols_affinity
, outcome_cols_affinity)]
cols_to_extract = cols_to_consider[cols_to_consider%in%c(common_cols
, outcome_cols_stability)]
##############################################################
#####################
# Ensemble stability: outcome_cols_stability