Skip to content

GitLab

Commit e67ee873 authored by Carine Rey's avatar Carine Rey
Browse files

fix errors if external usage

parent 6bce0845
Hide whitespace changes
Inline Side-by-side
Showing with 107 additions and 105 deletions
lib/scripts/calc_t_per_meth.R
View file @ e67ee873
......@@ -26,7 +26,7 @@ option_list = list(
make_option(c("--H0_NeG5_NeC_div2_res" ) , type="character", default=NA, help="merged_results H0_NeG5_NeC_div2_res" , metavar="character"),
make_option(c("--H0_NeG5_NeC_x2_res" ) , type="character", default=NA, help="merged_results H0_NeG5_NeC_x2_res" , metavar="character"),
make_option(c("--HaPCOC" ) , type="character", default=NA, help="merged_results HaPCOC" , metavar="character"),
make_option(c("--tree_prefix" ) , type="character", default="", help="tree_prefix" , metavar="character"),
......@@ -99,10 +99,10 @@ build_df_dist_couple = function (df_h0,df_ha,name) {
if ((! is.null(df_h0)) & (! is.null(df_ha))) {
id_vars = if ("P_JSD" %in% colnames(df_h0)) {c("Sites","P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2")} else {c("Sites")}
df_h0 = melt(df_h0, id.vars = id_vars, variable.name = "methode")
df_ha = melt(df_ha, id.vars = id_vars, variable.name = "methode")
df_ha = melt(df_ha, id.vars = id_vars, variable.name = "methode")
print(head(df_h0))
df = merge(df_h0, df_ha, by = c("Sites","P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2", "methode"), suffix = c("_H0","_Ha"))
df_melt = melt(df, id.vars = c("Sites","P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2", "methode"), variable.name = "val_H0Ha")
df = merge(df_h0, df_ha, by = c(id_vars, "methode"), suffix = c("_H0","_Ha"))
df_melt = melt(df, id.vars = c(id_vars, "methode"), variable.name = "val_H0Ha")
df_melt$couple = name
} else {
df_melt = NULL
......@@ -167,7 +167,7 @@ calc_TN_FP_TP_FN = function(meth, df_H0, df_Ha){
t_list =c(0,1,df_H0[,meth],df_Ha[,meth])
t_list = sort(t_list)
t_list = unique(t_list)
df_TN_FP = do.call(rbind.data.frame,lapply(t_list, calc_TN_FP, vals=df_H0[,meth]))
df_TP_FN = do.call(rbind.data.frame,lapply(t_list, calc_TP_FN, vals=df_Ha[,meth]))
......@@ -181,7 +181,7 @@ calc_TN_FP_TP_FN = function(meth, df_H0, df_Ha){
build_df_couple = function (df_h0,df_ha,name) {
print(name)
if ((! is.null(df_h0)) & (! is.null(df_ha))) {
id_vars = if ("P_JSD" %in% colnames(df_h0)) {c("Sites","P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2")} else {c("Sites")}
id_vars = if ("P_JSD" %in% colnames(df_h0)) {c("Sites","P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2")} else {c("Sites")}
meth_l = colnames(df_h0)[ ! colnames(df_h0) %in% id_vars]
df = do.call(rbind.data.frame,lapply(unique(meth_l) , calc_TN_FP_TP_FN, df_H0 = df_h0, df_Ha = df_ha))
df$couple = name
......@@ -190,7 +190,7 @@ build_df_couple = function (df_h0,df_ha,name) {
df = NULL
print(df)
}
return(df)
}
......@@ -349,7 +349,7 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
limitsize = FALSE
)
print("plot recall_precision")
if (length(meths) <= 7) {
plot = ggplot(df_out, aes(x=sensitivity, y=precision98_02, col = methode))
plot = plot + theme_bw()
......@@ -372,13 +372,13 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
colour ="black")) +
guides( colour = guide_legend(override.aes = list(alpha = 1,size = 3), nrow=2))
)
plot_recall_precision_papier = plot_grid(legend_PR, plot + theme(legend.position="none"), ncol = 1, scale = 1,
labels = c("",""),
align = "h",
rel_heights = c( 0.3, 0.7),
hjust = 0, vjust = 1)
save_plot(paste0(opt$out,suffix,".recall_precision.pdf"),
......@@ -391,7 +391,7 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
} else {
plot_recall_precision_papier = NULL
}
print("plot FPR")
#dcast
x_labs = "Threshold"
......@@ -400,7 +400,7 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
df_out_melt2$couple = gsub("/HaPC ","",df_out_melt2$couple)
df_out_melt2 = subset(df_out_melt2, variable == "specificity")
df_out_melt2$FPR = 1 - df_out_melt2$value
if (length(df_out_melt2$methode) > 0){
df_recall_sup09_per_meth2 = subset(df_recall_sup09_per_meth, couple %in%c("H0/HaPC NeG5", "H0/HaPC NeG5_NeC_div2" , "H0/HaPC NeG5_NeC_x2"))
df_recall_sup09_per_meth2$couple = gsub("/HaPC ","",df_recall_sup09_per_meth2$couple)
......@@ -428,21 +428,21 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
base_aspect_ratio = 1.5,
limitsize = FALSE
)
if (length(unique(df_out_melt2$couple)) == 3){
print(head(df_out_melt2))
df_out_melt2_dcast = dcast(df_out_melt2, methode + threshold ~ couple, value.var="FPR" )
print(head(df_out_melt2_dcast))
print(df_out_melt2_dcast)
# complete values for not existing threshold between couple
complete_value=function(df, col_name) {
print(colnames(df))
print(col_name)
df = df [,c("methode","threshold",col_name)]
colnames(df) = c("methode","threshold","value")
fix_val = function(r,df_no_na) {
x = r["value"]
t = r["threshold"]
......@@ -469,19 +469,19 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
df_out_melt2_dcast$d_H0NeG5_NeC_x2_H0NeG5 = df_out_melt2_dcast$H0NeG5_NeC_x2 - df_out_melt2_dcast$H0NeG5
df_out_melt2_dcast$d_H0NeG5_NeC_div2_H0NeG5 = df_out_melt2_dcast$H0NeG5_NeC_div2 - df_out_melt2_dcast$H0NeG5
print(df_out_melt2_dcast)
df_out_melt2_dcast_melt = melt(df_out_melt2_dcast[,c("methode", "threshold","d_H0NeG5_NeC_x2_H0NeG5","d_H0NeG5_NeC_div2_H0NeG5")], id.vars= c("methode", "threshold"))
print(head(df_out_melt2_dcast_melt))
df_out_melt2_dcast_melt$variable = factor(df_out_melt2_dcast_melt$variable, levels = c("H0NeG5","H0NeG5_NeC_x2","H0NeG5_NeC_div2","d_H0NeG5_NeC_x2_H0NeG5","d_H0NeG5_NeC_div2_H0NeG5"))
df_recall_sup09_per_meth3 = subset(df_recall_sup09_per_meth2, couple %in%c("H0NeG5"))
df_recall_sup09_per_meth3$variable = df_recall_sup09_per_meth3$couple
df_recall_sup09_per_meth3$variable = factor(df_recall_sup09_per_meth3$variable, levels = c("H0NeG5","H0NeG5_NeC_x2","H0NeG5_NeC_div2","d_H0NeG5_NeC_x2_H0NeG5","d_H0NeG5_NeC_div2_H0NeG5"))
x_labs = "Threshold"
y_labs = "Delta FP Rate"
......@@ -516,9 +516,8 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
plot = ggplot(df_out_melt, aes(x=threshold, y=value, col=couple))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top",
plot.margin = unit(x = c(0, 0, 0, 0), units = "cm"),
legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="black")) +
guides(colour = guide_legend(override.aes = list(alpha = 1), nrow=2))
legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="black")) +
guides(colour = guide_legend(override.aes = list(alpha = 1), nrow=2))
plot = plot + guides(fill=FALSE) + theme(legend.position="top")
plot = plot + labs(x=x_labs, y=y_labs) #+ ylim(y_lim)
......@@ -533,101 +532,104 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
plot_max_MCC = plot
save_plot(paste0(opt$out,suffix,".indicator_per_meth.pdf"),
plot_max_MCC,
ncol = 0.4 * length(unique(df_out_melt$methode)),
nrow = 1.7,
base_aspect_ratio = 1.5,
limitsize = FALSE
)
plot_max_MCC,
ncol = 0.4 * length(unique(df_out_melt$methode)),
nrow = 1.7,
base_aspect_ratio = 1.5,
limitsize = FALSE
)
print("plot value/distance")
if ( all((df_out_melt) %in% c("P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2"))){
print("plot value/distance")
plot = ggplot(df_d, aes(y=P_JSD, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="P JSD", x="Value")
plot = plot + xlim(c(0,1)) + ylim(c(0, max(max(df_d$P_JSD), 1)))
plot = plot + guides(col=FALSE)
plot = plot + scale_color_manual(values=colors)
plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 1, show.legend = NA, linetype="dashed")
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_value_JSD = plot
plot = ggplot(df_d, aes(y=P_ED, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="P ED", x="Value")
plot = plot + xlim(c(0,1)) + ylim(c(0, max(max(df_d$P_ED), 1)))
plot = plot + guides(col=FALSE)
plot = plot + scale_color_manual(values=colors)
plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 1, show.legend = NA, linetype="dashed")
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_value_ED = plot
plot = ggplot(subset(df_d,val_H0Ha=="value_Ha"), aes(y=entropy_C2, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="Entropy C2", x="Value")
plot = plot + xlim(c(0,1)) + ylim(c(0, max(max(df_d$entropy_C2), 1)))
plot = plot + guides(col=FALSE)
plot = plot + scale_color_manual(values=colors)
plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 1, show.legend = NA, linetype="dashed")
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_entropy_C2 = plot
plot = ggplot(df_d, aes(y=entropy_C1, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="Entropy C1", x="Value")
plot = plot + xlim(c(0,1)) + ylim(c(0, max(max(df_d$entropy_C2), 1)))
plot = plot + guides(col=FALSE)
plot = plot + scale_color_manual(values=colors)
plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 1, show.legend = NA, linetype="dashed")
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_entropy_C1 = plot
plot = ggplot(subset(df_d,val_H0Ha=="value_Ha") , aes(y=entropy_C2, x=entropy_C1, shape = val_H0Ha, col = value))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="Entropy C2", x="Entropy C1")
plot = plot + xlim(c(0, max(max(df_d$entropy_C1), 1))) + ylim(c(0, max(max(df_d$entropy_C2), 1)))
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + scale_colour_gradientn(colours = c("red","yellow","green","lightblue","darkblue"),
values=c(1.0,0.8,0.6,0.4,0.2,0))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_entropy_C2_C1 = plot
plot = ggplot(df_d, aes(y=P_JSD, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="P JSD", x="Value")
plot = plot + xlim(c(0,1)) + ylim(c(0, max(max(df_d$P_JSD), 1)))
plot = plot + guides(col=FALSE)
plot = plot + scale_color_manual(values=colors)
plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 1, show.legend = NA, linetype="dashed")
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_value_JSD = plot
plot = ggplot(df_d, aes(y=P_ED, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="P ED", x="Value")
plot = plot + xlim(c(0,1)) + ylim(c(0, max(max(df_d$P_ED), 1)))
plot = plot + guides(col=FALSE)
plot = plot + scale_color_manual(values=colors)
plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 1, show.legend = NA, linetype="dashed")
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_value_ED = plot
plot = ggplot(subset(df_d,val_H0Ha=="value_Ha"), aes(y=entropy_C2, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="Entropy C2", x="Value")
plot = plot + xlim(c(0,1)) + ylim(c(0, max(max(df_d$entropy_C2), 1)))
plot = plot + guides(col=FALSE)
plot = plot + scale_color_manual(values=colors)
plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 1, show.legend = NA, linetype="dashed")
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_entropy_C2 = plot
plot = ggplot(df_d, aes(y=entropy_C1, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="Entropy C1", x="Value")
plot = plot + xlim(c(0,1)) + ylim(c(0, max(max(df_d$entropy_C2), 1)))
plot = plot + guides(col=FALSE)
plot = plot + scale_color_manual(values=colors)
plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 1, show.legend = NA, linetype="dashed")
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_entropy_C1 = plot
plot = ggplot(subset(df_d,val_H0Ha=="value_Ha") , aes(y=entropy_C2, x=entropy_C1, shape = val_H0Ha, col = value))
plot = plot + theme_bw()
plot = plot + theme(legend.position="top")
plot = plot + labs(y="Entropy C2", x="Entropy C1")
plot = plot + xlim(c(0, max(max(df_d$entropy_C1), 1))) + ylim(c(0, max(max(df_d$entropy_C2), 1)))
plot = plot + geom_point(size=1.2, alpha=0.7, stroke=0.01)
plot = plot + scale_shape_manual(values = c(16,17))
plot = plot + scale_colour_gradientn(colours = c("red","yellow","green","lightblue","darkblue"),
values=c(1.0,0.8,0.6,0.4,0.2,0))
plot = plot + facet_grid(couple ~ methode)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_entropy_C2_C1 = plot
save_plot(paste0(opt$out,suffix,".value_distance_per_meth.pdf"),
plot_value_JSD,
ncol = 0.4* length(unique(df_d$methode)),
nrow = 0.4* length(unique(df_d$couple)),
base_aspect_ratio = 1,
limitsize = FALSE
)
save_plot(paste0(opt$out,suffix,".value_distance_per_meth.pdf"),
plot_value_JSD,
ncol = 0.4* length(unique(df_d$methode)),
nrow = 0.4* length(unique(df_d$couple)),
base_aspect_ratio = 1,
limitsize = FALSE
)
#plot = plot_grid(plot_recall_precision,plot_max_MCC,plot_value_JSD,plot_value_ED,plot_entropy_C2_C1,plot_entropy_C1,plot_entropy_C2,
# labels = c("A", "B","C","D","E","F"),
# rel_heights = c(length(unique(df_out$couple))*0.8,3,length(unique(df_out$couple))*0.8,length(unique(df_out$couple))*0.8,length(unique(df_out$couple))*0.8,length(unique(df_out$couple))*0.8,length(unique(df_out$couple))*0.8),
# nrow=7
# )
}
plot = plot_grid(plot_recall_precision,plot_recall_precision_papier,plot_max_MCC,plot_FPR,plot_FPR2,
labels = c("A", "A'", "B","B'","B''"),
rel_heights = c(length(unique(df_out$couple))*0.8,2,4,2,2),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment