Skip to content

GitLab

Commit 4bbbe597 authored by Carine Rey's avatar Carine Rey
Browse files

new plots

parent dd32b73a
Hide whitespace changes
Inline Side-by-side
Showing with 64 additions and 32 deletions
lib/scripts/calc_t_per_meth.R
View file @ 4bbbe597
......@@ -113,6 +113,7 @@ build_df_dist_couple = function (df_h0,df_ha,name) {
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")
print(name)
print(head(df_h0))
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")
......@@ -138,6 +139,9 @@ df_d_H0HaPC_NeG1_NeC_5 = build_df_dist_couple(df_H0_NeG1_NeC_5, df_HaPC_NeG1_N
df_d_H0HaPC_NeG5_indel = build_df_dist_couple(df_H0_NeG5_indel, df_HaPC_NeG5_indel, "H0/HaPC NeG5_indel")
df_d_H0H0_NeG1NeG1_NeC5 = build_df_dist_couple(df_H0_NeG1, df_H0_NeG1_NeC_5, "H0/H0 NeG1/NeG1_NeC_5")
df_d_H0H0_NeG5NeG5_NeC1 = build_df_dist_couple(df_H0_NeG5, df_H0_NeG5_NeC_1, "H0/H0 NeG5/NeG5_NeC_1")
df_d = rbind.data.frame(
df_d_H0HaPCOC_NeG1,
df_d_H0HaPC_NeG1,
......@@ -149,6 +153,8 @@ df_d = rbind.data.frame(
df_d_H0HaPC_NeG5_NeC_x2,
df_d_H0HaPC_NeG5_NeC_1,
df_d_H0HaPC_NeG1_NeC_5,
df_d_H0H0_NeG1NeG1_NeC5,
df_d_H0H0_NeG5NeG5_NeC1,
df_d_H0HaPC_NeG5_indel
)
......@@ -230,12 +236,16 @@ df_H0HaPC_NeG5_NeC_x2 = build_df_couple(df_H0_NeG5_NeC_x2, df_HaPC_NeG5_NeC_x2
df_H0HaPC_NeG5_NeC_1 = build_df_couple(df_H0_NeG5_NeC_1, df_HaPC_NeG5_NeC_1, "H0/HaPC NeG5_NeC_1")
df_H0HaPC_NeG1_NeC_5 = build_df_couple(df_H0_NeG1_NeC_5, df_HaPC_NeG1_NeC_5, "H0/HaPC NeG1_NeC_5")
df_H0H0_NeG5NeG5_NeC_1 = build_df_couple(df_H0_NeG5, df_H0_NeG5_NeC_1, "H0/H0 NeG5/NeG5_NeC_1")
df_H0H0_NeG1NeG1_NeC_5 = build_df_couple(df_H0_NeG1, df_H0_NeG1_NeC_5, "H0/H0 NeG1/NeG1_NeC_5")
df_H0HaPC_NeG5_indel = build_df_couple(df_H0_NeG5_indel, df_HaPC_NeG5_indel, "H0/HaPC NeG5_indel")
df_l = list(df_H0HaPCOC_NeG1, df_H0HaPC_NeG1, df_H0HaPC_NeG2, df_H0HaPC_NeG3,
df_H0HaPC_NeG4, df_H0HaPC_NeG5, df_H0HaPC_NeG5_NeC_div2,
df_H0HaPC_NeG5_NeC_1,df_H0HaPC_NeG1_NeC_5,
df_H0H0_NeG5NeG5_NeC_1, df_H0H0_NeG1NeG1_NeC_5,
df_H0HaPC_NeG5_NeC_x2,df_H0HaPC_NeG5_indel)
df_l = df_l[-which(sapply(df_l, is.null))]
......@@ -321,7 +331,7 @@ print(df_recall_sup09_per_meth)
plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffix="") {
nb_c = length(unique(df_out$couple))
colors = c(c("#984EA3","#4AA947","#377EB8","#E41A1C","#F5BE5B","#90EE90","#8B6914","#BFBFBF","#0000FF","#FFC0CB")[1:nb_c], c("#7F7F7F","#ADD8E6"))
colors = c(c("#984EA3","#4AA947","#377EB8","#E41A1C","#F5BE5B","#90EE90","#8B6914","#BFBFBF","#0000FF","#FFC0CB","#FFA500","#FFFF00")[1:nb_c], c("#7F7F7F","#ADD8E6"))
colors2 = c("#984EA3","#4AA947","#377EB8","#E41A1C","#F5BE5B","#90EE90","#8B6914","#BFBFBF","#7F7F7F","#ADD8E6")
if (! is.null(meths)) {
......@@ -380,7 +390,8 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
print("plot recall_precision")
if (length(meths) <= 9) {
plot = ggplot(df_out, aes(x=sensitivity, y=precision98_02, col = methode))
couple_i = c("H0/HaPC NeG1","H0/HaPC NeG3","H0/HaPC NeG4","H0/HaPC NeG5","H0/H0 NeG5/NeG5_NeC_1","H0/H0 NeG1/NeG1_NeC_5", "H0/HaPC NeG5_NeC_1","H0/HaPC NeG1_NeC_5")
plot = ggplot(subset(df_out, couple %in% couple_i), aes(x=sensitivity, y=precision98_02, col = methode))
plot = plot + theme_bw()
plot = plot + labs(x="Sensitivity (= Recall)", y="Precision (98/2)")
plot = plot + theme(legend.position="top")
......@@ -390,13 +401,12 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
#plot = plot + geom_point(size=1, alpha=alpha)
plot = plot + geom_step(direction="vh", size=0.5, alpha=alpha)
plot = plot + geom_hline( aes(yintercept = 0.9), col="black" , size = 0.5, show.legend = NA,linetype="dashed")
plot = plot + facet_grid(. ~ couple )
plot = plot + facet_wrap(~ couple, ncol=4)
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
legend_PR = get_legend(plot + theme(legend.position="top",
legend.text = element_text(size=10),
legend.title = element_text(size=0),
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,size = 3), nrow=2))
......@@ -417,8 +427,9 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
base_aspect_ratio = 1,
limitsize = FALSE
)
plot = ggplot(subset(df_out,couple%in%c("H0/HaPC NeG5", "H0/HaPC NeG5_NeC_div2" , "H0/HaPC NeG5_NeC_x2")), aes(x=sensitivity, y=precision98_02, col = methode))
plot_PR_c = function(couple_l) {
plot = ggplot(subset(df_out,couple%in%couple_l), aes(x=sensitivity, y=precision98_02, col = methode))
plot = plot + theme_bw()
plot = plot + labs(x="Sensitivity (= Recall)", y="Precision (98/2)")
plot = plot + theme(legend.position="top")
......@@ -430,28 +441,38 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
plot = plot + geom_hline( aes(yintercept = 0.9), col="black" , size = 0.5, show.legend = NA,linetype="dashed")
plot = plot + facet_grid(. ~ couple )
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot}
plot = plot_PR_c(c("H0/HaPC NeG1","H0/HaPC NeG3","H0/HaPC NeG4","H0/HaPC NeG5","H0/H0 NeG5/NeG5_NeC_1","H0/H0 NeG1/NeG1_NeC_5", "H0/HaPC NeG5_NeC_1","H0/HaPC NeG1_NeC_5"))
plot1 = plot_PR_c(c("H0/HaPC NeG1","H0/HaPC NeG3","H0/HaPC NeG4","H0/HaPC NeG5"))
plot2 = plot_PR_c(c("H0/H0 NeG5/NeG5_NeC_1","H0/H0 NeG1/NeG1_NeC_5"))
plot3 = plot_PR_c(c("H0/HaPC NeG5_NeC_1","H0/HaPC NeG1_NeC_5"))
legend_PR = get_legend(plot + theme(legend.position="top",
legend.text = element_text(size=10),
legend.title = element_text(size=0),
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,size = 3), nrow=2))
)
plot_recall_precision_papier_ok = plot_grid(legend_PR, plot + theme(legend.position="none"), ncol = 1, scale = 1,
labels = c("",""),
plot_recall_precision_papier_ok = plot_grid(legend_PR,
plot1 + theme(legend.position="none"),
plot2 + theme(legend.position="none"),
plot3 + theme(legend.position="none"),
ncol = 1, scale = 1,
labels = c("","Only Shift Pressure","Only Change in Selection Efficacity","Shift Pressure & Change in Selection Efficacity"),
align = "h",
rel_heights = c( 0.3, 0.7),
hjust = 0, vjust = 1)
rel_heights = c( 0.3, 0.7,0.7,0.7),
hjust = 0, vjust = 0)
save_plot(paste0(opt$out,suffix,".recall_precision_ok.pdf"),
plot_recall_precision_papier_ok,
ncol = 0.6* 3,
nrow = 1,
nrow = 2.1,
base_aspect_ratio = 1,
limitsize = FALSE
)
......@@ -463,13 +484,13 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
#dcast
x_labs = "Threshold"
y_labs = "FP Rate (1 -Spe)"
df_out_melt2 = subset(df_out_melt,couple%in%c("H0/HaPC NeG5", "H0/HaPC NeG5_NeC_div2" , "H0/HaPC NeG5_NeC_x2"))
df_out_melt2 = subset(df_out_melt,couple%in%c("H0/HaPC NeG5", "H0/HaPC NeG1", "H0/HaPC NeG5_NeC_1" , "H0/HaPC NeG1_NeC_5"))
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 = subset(df_recall_sup09_per_meth, couple %in%c("H0/HaPC NeG5", "H0/HaPC NeG1", "H0/HaPC NeG5_NeC_1" , "H0/HaPC NeG1_NeC_5"))
df_recall_sup09_per_meth2$couple = gsub("/HaPC ","",df_recall_sup09_per_meth2$couple)
plot = ggplot(df_out_melt2, aes(x=threshold, y=FPR, col=couple))
......@@ -496,11 +517,10 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
limitsize = FALSE
)
if (length(unique(df_out_melt2$couple)) == 3){
if (length(unique(df_out_melt2$couple)) == 4){
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
......@@ -530,24 +550,32 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
}
df_out_melt2_dcast$H0NeG5 = complete_value(df_out_melt2_dcast, "H0NeG5")
df_out_melt2_dcast$H0NeG5_NeC_x2 = complete_value(df_out_melt2_dcast, "H0NeG5_NeC_x2")
df_out_melt2_dcast$H0NeG5_NeC_div2 = complete_value(df_out_melt2_dcast, "H0NeG5_NeC_div2")
print(df_out_melt2_dcast)
df_out_melt2_dcast$H0NeG1 = complete_value(df_out_melt2_dcast, "H0NeG1")
#df_out_melt2_dcast$H0NeG5_NeC_x2 = complete_value(df_out_melt2_dcast, "H0NeG5_NeC_x2")
#df_out_melt2_dcast$H0NeG5_NeC_div2 = complete_value(df_out_melt2_dcast, "H0NeG5_NeC_div2")
df_out_melt2_dcast$H0NeG1_NeC_5 = complete_value(df_out_melt2_dcast, "H0NeG1_NeC_5")
df_out_melt2_dcast$H0NeG5_NeC_1 = complete_value(df_out_melt2_dcast, "H0NeG5_NeC_1")
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
#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$d_H0NeG1_NeC_5_H0NeG1 = df_out_melt2_dcast$H0NeG1_NeC_5 - df_out_melt2_dcast$H0NeG1
df_out_melt2_dcast$d_H0NeG5_NeC_1_H0NeG5 = df_out_melt2_dcast$H0NeG5_NeC_1 - df_out_melt2_dcast$H0NeG5
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"))
#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"))
df_out_melt2_dcast_melt = melt(df_out_melt2_dcast[,c("methode", "threshold","d_H0NeG5_NeC_1_H0NeG5","d_H0NeG1_NeC_5_H0NeG1")], 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_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_out_melt2_dcast_melt$variable = factor(df_out_melt2_dcast_melt$variable, levels = c("H0NeG5", "H0NeG1", "H0NeG1_NeC_5","H0NeG5_NeC_1","d_H0NeG5_NeC_1_H0NeG5","d_H0NeG1_NeC_5_H0NeG1"))
df_recall_sup09_per_meth3 = subset(df_recall_sup09_per_meth2, couple %in%c("H0NeG5"))
df_recall_sup09_per_meth3 = subset(df_recall_sup09_per_meth2, couple %in%c("H0NeG5","H0NeG1"))
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"))
df_recall_sup09_per_meth3$variable = factor(df_recall_sup09_per_meth3$variable, levels = c("H0NeG5", "H0NeG1", "H0NeG1_NeC_5","H0NeG5_NeC_1","d_H0NeG5_NeC_1_H0NeG5","d_H0NeG1_NeC_5_H0NeG1"))
x_labs = "Threshold"
......@@ -562,7 +590,7 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
plot = plot + ylim(c(-1,1)) + xlim(c(0,1))
plot = plot + geom_vline(data = df_recall_sup09_per_meth3 , aes(xintercept = threshold, col=variable), size = 0.5, show.legend = NA,linetype="dashed")
plot = plot + geom_step(data = df_out_melt2_dcast_melt, aes(x=threshold, y=value, col=variable), size=1, alpha=alpha)
plot = plot + scale_color_manual(values=colors[c(2,3,1)])
plot = plot + scale_color_manual(values=colors[c(2,4,1,3)])
plot = plot + facet_grid(. ~ methode, scales="free")
plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
......@@ -607,13 +635,13 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
)
plot_value_JSD = NULL
if ( all( c("P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2") %in% colnames(df_d))){
if ( all( c("P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2") %in% colnames(df_d)) & (length(meths) <= 9) ){
print("plot value/distance")
plot = ggplot(df_d, aes(y=P_JSD, x=value, shape = val_H0Ha, col = val_H0Ha))
plot = ggplot(subset(df_d, val_H0Ha == "value_Ha"), 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")
......@@ -627,7 +655,11 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
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))
print(head(df_d))
print(table(df_d$val_H0Ha,df_d$methode,df_d$couple))
plot = ggplot(subset(df_d, val_H0Ha == "value_Ha"), 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")
......@@ -701,7 +733,7 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
# )
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),
rel_heights = c(length(unique(df_out$couple))*0.8,4,4,2,2),
nrow=5)
if (is.null(plot_value_JSD)) {
......@@ -720,7 +752,7 @@ plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffi
save_plot(paste0(opt$out,suffix,".pdf"),
plot,
ncol = 0.4* length(unique(df_out_melt$methode)),
ncol = 0.5* length(unique(df_out_melt$methode)),
nrow = length(unique(df_out$couple)) * 0.5 + 2 + 1 +1 +1,
base_aspect_ratio = 1,
limitsize = FALSE
......
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