426 lines
16 KiB
R
426 lines
16 KiB
R
#!/usr/bin/env Rscript
|
|
#########################################################
|
|
# TASK: prediction
|
|
|
|
#=======================================================================
|
|
# working dir and loading libraries
|
|
getwd()
|
|
setwd("~/git/LSHTM_analysis/scripts/")
|
|
getwd()
|
|
|
|
source("plotting/combining_dfs_plotting.R")
|
|
|
|
####################################################################
|
|
# end of loading libraries and functions
|
|
####################################################################
|
|
|
|
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
# ps
|
|
table(merged_df2$mutation_info)
|
|
merged_df2$mutation_info_labels = ifelse(merged_df2$mutation_info == dr_muts_col, 1, 0)
|
|
table(merged_df2$mutation_info_labels)
|
|
|
|
|
|
table(merged_df3$mutation_info)
|
|
merged_df3$mutation_info_labels = ifelse(merged_df3$mutation_info == dr_muts_col, 1, 0)
|
|
table(merged_df3$mutation_info_labels)
|
|
|
|
|
|
# lig
|
|
table(merged_df2_lig$mutation_info)
|
|
merged_df2_lig$mutation_info_labels = ifelse(merged_df2_lig$mutation_info == dr_muts_col, 1, 0)
|
|
table(merged_df2_lig$mutation_info_labels)
|
|
|
|
|
|
table(merged_df3_lig$mutation_info)
|
|
merged_df3_lig$mutation_info_labels = ifelse(merged_df3_lig$mutation_info == dr_muts_col, 1, 0)
|
|
table(merged_df3_lig$mutation_info_labels)
|
|
|
|
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
###############################################################################
|
|
model_ind = glm(mutation_info_labels ~ or_mychisq
|
|
, data = merged_df2
|
|
, family = "binomial")
|
|
summary(model_ind)
|
|
nobs(model_ind)
|
|
|
|
#=============
|
|
# try loop
|
|
#=============
|
|
my_ivs = c("or_mychisq", "or_kin", "pval_fisher", "af", "duet_stability_change", "duet_scaled")
|
|
|
|
for( i in my_ivs){
|
|
cat ("===============================\n")
|
|
cat(i)
|
|
cat ("\n===============================\n")
|
|
print(summary(glm(mutation_info_labels ~ eval(parse(text=i))
|
|
, data = merged_df2
|
|
, family = "binomial")))
|
|
}
|
|
###############################################################################
|
|
#========================================
|
|
# merged_df2: UNadjusted,loop
|
|
#========================================
|
|
my_ivs = c("or_mychisq", "or_kin", "pval_fisher", "af"
|
|
, "ligand_distance"
|
|
, "rsa"
|
|
, "rd_values"
|
|
, "kd_values"
|
|
, "duet_stability_change"
|
|
, "duet_scaled"
|
|
, "duet_outcome"
|
|
, "ddg"
|
|
, "foldx_scaled"
|
|
, "foldx_outcome")
|
|
|
|
ps_logistic_df2 = data.frame()
|
|
|
|
for( i in my_ivs){
|
|
print(i)
|
|
|
|
df = data.frame(var_name = NA
|
|
, number_samples = NA
|
|
, beta = NA
|
|
, odds_ratio = NA
|
|
, pvalue = NA
|
|
, se = NA
|
|
, zvalue = NA
|
|
, ci_lower = NA
|
|
, ci_upper = NA)
|
|
|
|
model = glm(mutation_info_labels ~ eval(parse(text=i))
|
|
, data = merged_df2
|
|
, family = "binomial")
|
|
|
|
var_name = i
|
|
|
|
number_samples = nobs(model)
|
|
|
|
beta_logistic = summary(model)$coefficients[2,1]; beta_logistic
|
|
|
|
or_logistic = exp(summary(model)$coefficients[2,1])
|
|
|
|
pval_logistic = summary(model)$coefficients[2,4]
|
|
|
|
se_logistic = summary(model)$coefficients[2,2]
|
|
|
|
zval_logistic = summary(model)$coefficients[2,3]
|
|
|
|
ci_mod = exp(confint(model))[2,]
|
|
|
|
ci_lower_logistic = ci_mod[["2.5 %"]]
|
|
ci_upper_logistic = ci_mod[["97.5 %"]]
|
|
|
|
print(c(var_name, beta_logistic, or_logistic, pval_logistic, se_logistic, zval_logistic, ci_mod))
|
|
|
|
df$var_name = var_name
|
|
df$number_samples = number_samples
|
|
df$beta = beta_logistic
|
|
df$odds_ratio = or_logistic
|
|
df$pvalue = pval_logistic
|
|
df$se = se_logistic
|
|
df$zvalue = zval_logistic
|
|
df$ci_lower = ci_lower_logistic
|
|
df$ci_upper = ci_upper_logistic
|
|
|
|
print(df)
|
|
ps_logistic_df2 = rbind(ps_logistic_df2, df)
|
|
|
|
}
|
|
#--------------------
|
|
# formatting df
|
|
#--------------------
|
|
ps_logistic_df2$data_source = "df2"
|
|
ps_logistic_df2$model = "unadjusted"
|
|
|
|
ps_logistic_df2$odds_ratio = round(ps_logistic_df2$odds_ratio, 2)
|
|
ps_logistic_df2$ci_lower = round(ps_logistic_df2$ci_lower, 2)
|
|
ps_logistic_df2$ci_upper = round(ps_logistic_df2$ci_upper, 2)
|
|
|
|
# adding pvalue_signif
|
|
ps_logistic_df2$pvalue_signif = ps_logistic_df2$pvalue
|
|
str(ps_logistic_df2$pvalue_signif)
|
|
|
|
ps_logistic_df2 = dplyr::mutate(ps_logistic_df2
|
|
, 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'))
|
|
# rearranging columns
|
|
colnames(ps_logistic_df2)
|
|
ps_logistic_df2_o = ps_logistic_df2 [c("var_name"
|
|
, "number_samples"
|
|
, "model"
|
|
, "odds_ratio"
|
|
, "pvalue"
|
|
, "pvalue_signif"
|
|
, "beta"
|
|
, "se"
|
|
, "zvalue"
|
|
, "ci_lower"
|
|
, "ci_upper"
|
|
, "data_source")]
|
|
|
|
###############################################################################
|
|
#========================================
|
|
# merged_df2: adjusted, loop
|
|
#========================================
|
|
#model_adjusted = glm(mutation_info_labels ~ or_mychisq + rsa + rd_values + kd_values +
|
|
# duet_stability_change + duet_scaled + duet_outcome + ddg + foldx_scaled + foldx_outcome
|
|
# , data = merged_df2
|
|
# , family = "binomial")
|
|
|
|
model_adjusted_df2 = glm(mutation_info_labels ~ or_mychisq + or_kin + rd_values + kd_values +
|
|
ligand_distance + duet_stability_change
|
|
, data = merged_df2
|
|
, family = "binomial");summary(model_adjusted_df2)
|
|
|
|
var_names_df = as.data.frame(names(model_adjusted_df2$coefficients))
|
|
names(var_names_df) = c("var_name")
|
|
|
|
ci_mod = exp(confint(model_adjusted_df2))
|
|
ci_mod_df = as.data.frame(ci_mod)
|
|
names(ci_mod_df) = c("ci_lower", "ci_upper")
|
|
ci_mod_df$ci_lower = round(ci_mod_df$ci_lower, 2)
|
|
ci_mod_df$ci_upper = round(ci_mod_df$ci_upper, 2)
|
|
|
|
|
|
estimates_df = as.data.frame(summary(model_adjusted_df2)$coefficients)
|
|
names(estimates_df) = c("beta", "se", "zvalue", "pvalue")
|
|
estimates_df$odds_ratio = round(exp(estimates_df$beta), 2)
|
|
number_samples = nobs(model_adjusted_df2)
|
|
estimates_df$number_samples = number_samples
|
|
estimates_df$data_source = "df2"
|
|
estimates_df$model = "adjusted"
|
|
|
|
names(ps_logistic_adjusted_df2)
|
|
|
|
if ( all(rownames(estimates_df) == rownames(ci_mod_df)) ){
|
|
cat("PASS: rownames match. Preparing to merge...")
|
|
ps_logistic_adjusted_df2 = merge(estimates_df, ci_mod_df, by = "row.names", all = T)
|
|
}
|
|
|
|
|
|
colnames(ps_logistic_adjusted_df2)[1] <- c("var_name")
|
|
d1 = which(ps_logistic_adjusted_df2$var_name == "(Intercept)")
|
|
ps_logistic_adjusted_df2 = ps_logistic_adjusted_df2[-d1,]
|
|
|
|
# adding pvalue_signif
|
|
ps_logistic_adjusted_df2$pvalue_signif = ps_logistic_adjusted_df2$pvalue
|
|
str(ps_logistic_adjusted_df2$pvalue_signif)
|
|
|
|
ps_logistic_adjusted_df2 = dplyr::mutate(ps_logistic_adjusted_df2
|
|
, 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'))
|
|
# rearranging columns
|
|
colnames(ps_logistic_adjusted_df2)
|
|
ps_logistic_adjusted_df2_o = ps_logistic_adjusted_df2[c("var_name"
|
|
, "number_samples"
|
|
, "model"
|
|
, "odds_ratio"
|
|
, "pvalue"
|
|
, "pvalue_signif"
|
|
, "beta"
|
|
, "se"
|
|
, "zvalue"
|
|
, "ci_lower"
|
|
, "ci_upper"
|
|
, "data_source")]
|
|
|
|
###############################################################################
|
|
|
|
###############################################################################
|
|
#========================================
|
|
# merged_df3: UNadjusted,loop
|
|
#========================================
|
|
|
|
my_ivs = c("or_mychisq", "or_kin", "pval_fisher", "af"
|
|
, "ligand_distance"
|
|
, "rsa"
|
|
, "rd_values"
|
|
, "kd_values"
|
|
, "duet_stability_change"
|
|
, "duet_scaled"
|
|
, "duet_outcome"
|
|
, "ddg"
|
|
, "foldx_scaled"
|
|
, "foldx_outcome")
|
|
|
|
ps_logistic_df3 = data.frame()
|
|
|
|
for( i in my_ivs){
|
|
print(i)
|
|
|
|
df = data.frame(var_name = NA
|
|
, number_samples = NA
|
|
, beta = NA
|
|
, odds_ratio = NA
|
|
, pvalue = NA
|
|
, se = NA
|
|
, zvalue = NA
|
|
, ci_lower = NA
|
|
, ci_upper = NA)
|
|
|
|
model = glm(mutation_info_labels ~ eval(parse(text=i))
|
|
, data = merged_df3
|
|
, family = "binomial")
|
|
|
|
var_name = i
|
|
|
|
number_samples = nobs(model)
|
|
|
|
beta_logistic = summary(model)$coefficients[2,1]; beta_logistic
|
|
|
|
or_logistic = exp(summary(model)$coefficients[2,1])
|
|
|
|
pval_logistic = summary(model)$coefficients[2,4]
|
|
|
|
se_logistic = summary(model)$coefficients[2,2]
|
|
|
|
zval_logistic = summary(model)$coefficients[2,3]
|
|
|
|
ci_mod = exp(confint(model))[2,]
|
|
|
|
ci_lower_logistic = ci_mod[["2.5 %"]]
|
|
ci_upper_logistic = ci_mod[["97.5 %"]]
|
|
|
|
print(c(var_name, beta_logistic, or_logistic, pval_logistic, se_logistic, zval_logistic, ci_mod))
|
|
|
|
df$var_name = var_name
|
|
df$number_samples = number_samples
|
|
df$beta = beta_logistic
|
|
df$odds_ratio = or_logistic
|
|
df$pvalue = pval_logistic
|
|
df$se = se_logistic
|
|
df$zvalue = zval_logistic
|
|
df$ci_lower = ci_lower_logistic
|
|
df$ci_upper = ci_upper_logistic
|
|
|
|
print(df)
|
|
ps_logistic_df3 = rbind(ps_logistic_df3, df)
|
|
|
|
}
|
|
#--------------------
|
|
# formatting df
|
|
#--------------------
|
|
ps_logistic_df3$data_source = "df3"
|
|
ps_logistic_df3$model = "unadjusted"
|
|
|
|
ps_logistic_df3$odds_ratio = round(ps_logistic_df3$odds_ratio, 2)
|
|
ps_logistic_df3$ci_lower = round(ps_logistic_df3$ci_lower, 2)
|
|
ps_logistic_df3$ci_upper = round(ps_logistic_df3$ci_upper, 2)
|
|
|
|
# adding pvalue_signif
|
|
ps_logistic_df3$pvalue_signif = ps_logistic_df3$pvalue
|
|
str(ps_logistic_df3$pvalue_signif)
|
|
|
|
ps_logistic_df3 = dplyr::mutate(ps_logistic_df3
|
|
, 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'))
|
|
# rearranging columns
|
|
ps_logistic_df3_o = ps_logistic_df3[c("var_name"
|
|
, "number_samples"
|
|
, "model"
|
|
, "odds_ratio"
|
|
, "pvalue"
|
|
, "pvalue_signif"
|
|
, "beta"
|
|
, "se"
|
|
, "zvalue"
|
|
, "ci_lower"
|
|
, "ci_upper"
|
|
, "data_source")]
|
|
|
|
#========================================
|
|
# merged_df3: adjusted, loop
|
|
#========================================
|
|
#model_adjusted_df3 = glm(mutation_info_labels ~ or_mychisq + rsa + rd_values + kd_values +
|
|
# duet_stability_change + duet_scaled + duet_outcome + ddg + foldx_scaled + foldx_outcome
|
|
# , data = merged_df3
|
|
# , family = "binomial")
|
|
|
|
model_adjusted_df3 = glm(mutation_info_labels ~ rd_values +
|
|
ligand_distance + duet_stability_change
|
|
, data = merged_df3
|
|
, family = "binomial");summary(model_adjusted_df3)
|
|
|
|
var_names_df = as.data.frame(names(model_adjusted_df3$coefficients))
|
|
names(var_names_df) = c("var_name")
|
|
|
|
ci_mod = exp(confint(model_adjusted_df3))
|
|
ci_mod_df = as.data.frame(ci_mod)
|
|
names(ci_mod_df) = c("ci_lower", "ci_upper")
|
|
ci_mod_df$ci_lower = round(ci_mod_df$ci_lower, 2)
|
|
ci_mod_df$ci_upper = round(ci_mod_df$ci_upper, 2)
|
|
|
|
|
|
estimates_df = as.data.frame(summary(model_adjusted_df3)$coefficients)
|
|
names(estimates_df) = c("beta", "se", "zvalue", "pvalue")
|
|
estimates_df$odds_ratio = round(exp(estimates_df$beta), 2)
|
|
number_samples = nobs(model_adjusted_df3)
|
|
estimates_df$number_samples = number_samples
|
|
estimates_df$data_source = "df3"
|
|
estimates_df$model = "adjusted"
|
|
|
|
names(ps_logistic_adjusted_df3)
|
|
|
|
if ( all(rownames(estimates_df) == rownames(ci_mod_df)) ){
|
|
cat("PASS: rownames match. Preparing to merge...")
|
|
ps_logistic_adjusted_df3 = merge(estimates_df, ci_mod_df, by = "row.names", all = T)
|
|
}
|
|
|
|
colnames(ps_logistic_adjusted_df3)[1] <- c("var_name")
|
|
d2 = which(ps_logistic_adjusted_df3$var_name == "(Intercept)")
|
|
ps_logistic_adjusted_df3 = ps_logistic_adjusted_df3[-d2,]
|
|
|
|
# adding pvalue_signif
|
|
ps_logistic_adjusted_df3$pvalue_signif = ps_logistic_adjusted_df3$pvalue
|
|
str(ps_logistic_adjusted_df3$pvalue_signif)
|
|
|
|
ps_logistic_adjusted_df3 = dplyr::mutate(ps_logistic_adjusted_df3
|
|
, 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'))
|
|
# rearranging columns
|
|
colnames(ps_logistic_adjusted_df3)
|
|
ps_logistic_adjusted_df3_o = ps_logistic_adjusted_df3[c("var_name"
|
|
, "number_samples"
|
|
, "model"
|
|
, "odds_ratio"
|
|
, "pvalue"
|
|
, "pvalue_signif"
|
|
, "beta"
|
|
, "se"
|
|
, "zvalue"
|
|
, "ci_lower"
|
|
, "ci_upper"
|
|
,"data_source")]
|
|
|
|
#-------------
|
|
# lm
|
|
#-------------
|
|
|
|
model_lm = lm(or_kin ~ rsa + rd_values + duet_stability_change + ddg + mutation_info_labels
|
|
, data = merged_df3)
|
|
|
|
summary(model_lm)
|
|
|
|
|
|
model_lm1 = lm(or_mychisq ~ mutation_info_labels
|
|
, data = merged_df2)
|
|
|
|
summary(model_lm1)
|