calc_t_per_meth.R 35.1 KB
Newer Older
Carine Rey's avatar
Carine Rey committed
1 2 3 4 5 6 7 8 9 10
#!/usr/bin/env Rscript

#https://www.r-bloggers.com/passing-arguments-to-an-r-script-from-command-lines/


library("optparse")
library("reshape2")
library("ggplot2")
library("cowplot")

Carine Rey's avatar
Carine Rey committed
11 12
date = format(Sys.time(), format="%Y-%m-%d %X")

Carine Rey's avatar
Carine Rey committed
13
option_list = list(
14 15 16 17 18 19 20 21 22 23
  make_option(c("--H0_NeG1"         ) , type="character", default=NA, help="merged_results H0_NeG1)"        , metavar="character"),
  make_option(c("--H0_NeG2"         ) , type="character", default=NA, help="merged_results H0_NeG2)"        , metavar="character"),
  make_option(c("--H0_NeG3"         ) , type="character", default=NA, help="merged_results H0_NeG3)"        , metavar="character"),
  make_option(c("--H0_NeG4"         ) , type="character", default=NA, help="merged_results H0_NeG4)"        , metavar="character"),
  make_option(c("--H0_NeG5"         ) , type="character", default=NA, help="merged_results H0_NeG5)"        , metavar="character"),
  make_option(c("--HaPC_NeG1"       ) , type="character", default=NA, help="merged_results HaPC_NeG1"       , metavar="character"),
  make_option(c("--HaPC_NeG2"       ) , type="character", default=NA, help="merged_results HaPC_NeG2"       , metavar="character"),
  make_option(c("--HaPC_NeG3"       ) , type="character", default=NA, help="merged_results HaPC_NeG3"       , metavar="character"),
  make_option(c("--HaPC_NeG4"       ) , type="character", default=NA, help="merged_results HaPC_NeG4"       , metavar="character"),
  make_option(c("--HaPC_NeG5"       ) , type="character", default=NA, help="merged_results HaPC_NeG5"       , metavar="character"),
24 25 26 27
  make_option(c("--HaPC_NeG5_NeC_div2_res"       ) , type="character", default=NA, help="merged_results HaPC_NeG5_NeC_div2_res"       , metavar="character"),
  make_option(c("--HaPC_NeG5_NeC_x2_res"         ) , type="character", default=NA, help="merged_results HaPC_NeG5_NeC_x2_res"         , metavar="character"),
  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"),
Carine Rey's avatar
Carine Rey committed
28 29 30 31 32
  make_option(c("--HaPC_NeG5_NeC_1_res"       ) , type="character", default=NA, help="merged_results HaPC_NeG5_NeC_1_res"       , metavar="character"),
  make_option(c("--H0_NeG5_NeC_1_res"         ) , type="character", default=NA, help="merged_results H0_NeG5_NeC_1_res"         , metavar="character"),
  make_option(c("--HaPC_NeG1_NeC_5_res"         ) , type="character", default=NA, help="merged_results HaPC_NeG1_NeC_5_res"         , metavar="character"),
  make_option(c("--H0_NeG1_NeC_5_res"           ) , type="character", default=NA, help="merged_results H0_NeG1_NeC_5_res"           , metavar="character"),
  
Carine Rey's avatar
Carine Rey committed
33 34
  make_option(c("--H0_NeG5_indel_res"         ) , type="character", default=NA, help="merged_results H0_NeG5_indel_res"         , metavar="character"),
  make_option(c("--HaPC_NeG5_indel_res"           ) , type="character", default=NA, help="merged_results HaPC_NeG5_indel_res"           , metavar="character"),
35
  make_option(c("--HaPCOC"          ) , type="character", default=NA, help="merged_results HaPCOC"          , metavar="character"),
Carine Rey's avatar
Carine Rey committed
36

Carine Rey's avatar
Carine Rey committed
37
  make_option(c("--tree_prefix" ) , type="character", default="", help="tree_prefix"          , metavar="character"),
38

39 40

  make_option(c("-o","--out"), type="character", default="out",
Carine Rey's avatar
Carine Rey committed
41 42 43 44 45 46
              help="output prefix [default= %default]", metavar="character")
);

opt_parser = OptionParser(option_list=option_list);
opt = parse_args(opt_parser);

47 48
print(opt)

49
if (is.na(opt$H0_NeG1)){
Carine Rey's avatar
Carine Rey committed
50
  print_help(opt_parser)
51
  stop("At least one argument must be supplied (H0_NeG1 or H0 input file)", call.=FALSE)
Carine Rey's avatar
Carine Rey committed
52
}
53
if (is.null(opt$HaPCOC)){
Carine Rey's avatar
Carine Rey committed
54
  print_help(opt_parser)
55
  stop("At least one argument must be supplied (HaPCOC input file)", call.=FALSE)
Carine Rey's avatar
Carine Rey committed
56 57
}

Carine Rey's avatar
Carine Rey committed
58 59 60 61 62 63 64


########################################################################

print("prep input files")

read_hyp = function(opt_name) {
65 66
  print(opt_name)
  if (! is.na(opt_name)){
Carine Rey's avatar
Carine Rey committed
67 68
    df   = read.table(opt_name, header=TRUE, sep = '\t',na.strings = "NA")
    df = df[, ! colnames(df) %in% c("Indel_prop", "Indel_prop.ConvLeaves.")]
Carine Rey's avatar
Carine Rey committed
69 70 71
    #id_vars = if ("P_JSD" %in% colnames(df)) {c("Sites","P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2")} else {c("Sites")}
    #df_melt = melt(df, id.vars = id_vars, variable.name = "methode")
    return(df)
Carine Rey's avatar
Carine Rey committed
72 73 74 75 76
  } else {
    return(NULL)
  }
}

Carine Rey's avatar
Carine Rey committed
77 78 79 80 81 82 83 84 85 86
df_H0_NeG1        = read_hyp(opt$H0_NeG1   )
df_H0_NeG2        = read_hyp(opt$H0_NeG2   )
df_H0_NeG3        = read_hyp(opt$H0_NeG3   )
df_H0_NeG4        = read_hyp(opt$H0_NeG4   )
df_H0_NeG5        = read_hyp(opt$H0_NeG5   )
df_HaPC_NeG1      = read_hyp(opt$HaPC_NeG1 )
df_HaPC_NeG2      = read_hyp(opt$HaPC_NeG2 )
df_HaPC_NeG3      = read_hyp(opt$HaPC_NeG3 )
df_HaPC_NeG4      = read_hyp(opt$HaPC_NeG4 )
df_HaPC_NeG5      = read_hyp(opt$HaPC_NeG5 )
87

Carine Rey's avatar
Carine Rey committed
88 89 90 91
df_HaPC_NeG5_NeC_div2      = read_hyp(opt$HaPC_NeG5_NeC_div2)
df_HaPC_NeG5_NeC_x2        = read_hyp(opt$HaPC_NeG5_NeC_x2  )
df_H0_NeG5_NeC_div2      = read_hyp(opt$H0_NeG5_NeC_div2  )
df_H0_NeG5_NeC_x2        = read_hyp(opt$H0_NeG5_NeC_x2    )
92

Carine Rey's avatar
Carine Rey committed
93 94 95 96 97
df_HaPC_NeG5_NeC_1      = read_hyp(opt$HaPC_NeG5_NeC_1)
df_H0_NeG5_NeC_1      = read_hyp(opt$H0_NeG5_NeC_1  )
df_HaPC_NeG1_NeC_5        = read_hyp(opt$HaPC_NeG1_NeC_5  )
df_H0_NeG1_NeC_5        = read_hyp(opt$H0_NeG1_NeC_5    )

Carine Rey's avatar
Carine Rey committed
98 99
df_HaPC_NeG5_indel        = read_hyp(opt$HaPC_NeG5_indel_res  )
df_H0_NeG5_indel      = read_hyp(opt$H0_NeG5_indel_res  )
Carine Rey's avatar
Carine Rey committed
100

Carine Rey's avatar
Carine Rey committed
101
df_HaPCOC = read_hyp(opt$HaPCOC)
102

Carine Rey's avatar
Carine Rey committed
103

104 105 106



Carine Rey's avatar
Carine Rey committed
107 108 109 110 111 112
########################################################################

print("prep df_d")

build_df_dist_couple = function (df_h0,df_ha,name) {
  if ((! is.null(df_h0)) & (! is.null(df_ha))) {
Carine Rey's avatar
Carine Rey committed
113 114
    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")
Carine Rey's avatar
Carine Rey committed
115
    df_ha = melt(df_ha, id.vars = id_vars, variable.name = "methode")
Carine Rey's avatar
Carine Rey committed
116
    print(head(df_h0))
Carine Rey's avatar
Carine Rey committed
117 118
    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")
Carine Rey's avatar
Carine Rey committed
119 120 121 122 123 124 125
    df_melt$couple = name
  } else {
    df_melt = NULL
  }
  return(df_melt)
}

Carine Rey's avatar
Carine Rey committed
126 127 128 129 130 131
df_d_H0HaPCOC_NeG1 = build_df_dist_couple(df_H0_NeG1, df_HaPCOC, "H0/HaPCOC NeG1")
df_d_H0HaPC_NeG1   = build_df_dist_couple(df_H0_NeG1, df_HaPC_NeG1, "H0/HaPC NeG1")
df_d_H0HaPC_NeG2   = build_df_dist_couple(df_H0_NeG2, df_HaPC_NeG2, "H0/HaPC NeG2")
df_d_H0HaPC_NeG3   = build_df_dist_couple(df_H0_NeG3, df_HaPC_NeG3, "H0/HaPC NeG3")
df_d_H0HaPC_NeG4   = build_df_dist_couple(df_H0_NeG4, df_HaPC_NeG4, "H0/HaPC NeG4")
df_d_H0HaPC_NeG5   = build_df_dist_couple(df_H0_NeG5, df_HaPC_NeG5, "H0/HaPC NeG5")
132

Carine Rey's avatar
Carine Rey committed
133 134
df_d_H0HaPC_NeG5_NeC_div2   = build_df_dist_couple(df_H0_NeG5_NeC_div2, df_HaPC_NeG5_NeC_div2, "H0/HaPC NeG5_NeC_div2")
df_d_H0HaPC_NeG5_NeC_x2   = build_df_dist_couple(df_H0_NeG5_NeC_x2, df_HaPC_NeG5_NeC_x2, "H0/HaPC NeG5_NeC_x2")
135

Carine Rey's avatar
Carine Rey committed
136 137 138
df_d_H0HaPC_NeG5_NeC_1   = build_df_dist_couple(df_H0_NeG5_NeC_1, df_HaPC_NeG5_NeC_1, "H0/HaPC NeG5_NeC_1")
df_d_H0HaPC_NeG1_NeC_5   = build_df_dist_couple(df_H0_NeG1_NeC_5, df_HaPC_NeG1_NeC_5, "H0/HaPC NeG5_NeC_5")

Carine Rey's avatar
Carine Rey committed
139 140
df_d_H0HaPC_NeG5_indel   = build_df_dist_couple(df_H0_NeG5_indel, df_HaPC_NeG5_indel, "H0/HaPC NeG5_indel")

141
df_d = rbind.data.frame(
142 143 144 145 146
                      df_d_H0HaPCOC_NeG1,
                      df_d_H0HaPC_NeG1,
                      df_d_H0HaPC_NeG2,
                      df_d_H0HaPC_NeG3,
                      df_d_H0HaPC_NeG4,
147 148
                      df_d_H0HaPC_NeG5,
                      df_d_H0HaPC_NeG5_NeC_div2,
Carine Rey's avatar
Carine Rey committed
149
                      df_d_H0HaPC_NeG5_NeC_x2,
Carine Rey's avatar
Carine Rey committed
150 151
                      df_d_H0HaPC_NeG5_NeC_1,
                      df_d_H0HaPC_NeG1_NeC_5,
Carine Rey's avatar
Carine Rey committed
152 153
                      df_d_H0HaPC_NeG5_indel
                      )
Carine Rey's avatar
Carine Rey committed
154 155 156 157 158 159 160 161 162 163

df_d = df_d[order(df_d$methode),]

print(head(df_d))
print(tail(df_d))

########################################################################

print("prep df")

Carine Rey's avatar
Carine Rey committed
164 165
## fun...

Carine Rey's avatar
Carine Rey committed
166
calc_TN_FP = function(t, vals){
Carine Rey's avatar
Carine Rey committed
167 168
    TN = 0
    FP = 0
169
    vals = na.omit(vals)
Carine Rey's avatar
Carine Rey committed
170
    if (length(vals) > 0) {
Carine Rey's avatar
Carine Rey committed
171 172
        TN = sum(vals < t)
        FP = sum(vals >= t)
Carine Rey's avatar
Carine Rey committed
173
    }
Carine Rey's avatar
Carine Rey committed
174
    return(data.frame(t=t, TN = TN, FP = FP))
Carine Rey's avatar
Carine Rey committed
175 176
}

Carine Rey's avatar
Carine Rey committed
177
calc_TP_FN = function(t, vals){
Carine Rey's avatar
Carine Rey committed
178 179 180 181
    FN = 0
    TP = 0
    vals = na.omit(vals)
    if (length(vals) > 0) {
Carine Rey's avatar
Carine Rey committed
182 183
        FN = sum(vals < t)
        TP = sum(vals >= t)
Carine Rey's avatar
Carine Rey committed
184
    }
Carine Rey's avatar
Carine Rey committed
185
    return(data.frame(t=t, TP = TP, FN = FN))
Carine Rey's avatar
Carine Rey committed
186 187
}

Carine Rey's avatar
Carine Rey committed
188 189 190 191
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)
Carine Rey's avatar
Carine Rey committed
192

Carine Rey's avatar
Carine Rey committed
193 194 195 196 197
    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]))

    df_TN_FP_TP_FN = merge(df_TN_FP, df_TP_FN, by="t")
    df_TN_FP_TP_FN$methode = meth
Carine Rey's avatar
Carine Rey committed
198 199 200 201 202

    return(df_TN_FP_TP_FN)
}


203
build_df_couple = function (df_h0,df_ha,name) {
204
  print(name)
205
  if ((! is.null(df_h0)) & (! is.null(df_ha))) {
Carine Rey's avatar
Carine Rey committed
206
    id_vars = if ("P_JSD" %in% colnames(df_h0)) {c("Sites","P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2")} else {c("Sites")}
Carine Rey's avatar
Carine Rey committed
207 208
    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))
209
    df$couple = name
210
    print(head(df))
211 212
  } else {
    df = NULL
213
    print(df)
214
  }
Carine Rey's avatar
Carine Rey committed
215

216 217
  return(df)
}
Carine Rey's avatar
Carine Rey committed
218

Carine Rey's avatar
Carine Rey committed
219 220 221 222 223 224
df_H0HaPCOC_NeG1 = build_df_couple(df_H0_NeG1, df_HaPCOC , "H0/HaPCOC NeG1")
df_H0HaPC_NeG1   = build_df_couple(df_H0_NeG1, df_HaPC_NeG1, "H0/HaPC NeG1")
df_H0HaPC_NeG2   = build_df_couple(df_H0_NeG2, df_HaPC_NeG2, "H0/HaPC NeG2")
df_H0HaPC_NeG3   = build_df_couple(df_H0_NeG3, df_HaPC_NeG3, "H0/HaPC NeG3")
df_H0HaPC_NeG4   = build_df_couple(df_H0_NeG4, df_HaPC_NeG4, "H0/HaPC NeG4")
df_H0HaPC_NeG5   = build_df_couple(df_H0_NeG5, df_HaPC_NeG5, "H0/HaPC NeG5")
225 226


Carine Rey's avatar
Carine Rey committed
227 228
df_H0HaPC_NeG5_NeC_div2   = build_df_couple(df_H0_NeG5_NeC_div2, df_HaPC_NeG5_NeC_div2, "H0/HaPC NeG5_NeC_div2")
df_H0HaPC_NeG5_NeC_x2   = build_df_couple(df_H0_NeG5_NeC_x2, df_HaPC_NeG5_NeC_x2, "H0/HaPC NeG5_NeC_x2")
229

Carine Rey's avatar
Carine Rey committed
230 231 232
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")

Carine Rey's avatar
Carine Rey committed
233 234 235
df_H0HaPC_NeG5_indel   = build_df_couple(df_H0_NeG5_indel, df_HaPC_NeG5_indel, "H0/HaPC NeG5_indel")


Carine Rey's avatar
Carine Rey committed
236 237
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,
Carine Rey's avatar
Carine Rey committed
238
            df_H0HaPC_NeG5_NeC_1,df_H0HaPC_NeG1_NeC_5,
Carine Rey's avatar
Carine Rey committed
239
            df_H0HaPC_NeG5_NeC_x2,df_H0HaPC_NeG5_indel)
240 241 242

df_l = df_l[-which(sapply(df_l, is.null))]
df =  do.call("rbind",df_l)
243 244 245 246

print(head(df))
print(tail(df))

Carine Rey's avatar
Carine Rey committed
247 248

## Sensitivity (= recall)
Carine Rey's avatar
Carine Rey committed
249 250 251
df$sens = df$TP / (df$TP+df$FN)
df$sens[is.na(df$sens)] = 0

252
## Specificity
Carine Rey's avatar
Carine Rey committed
253 254 255
df$spe =  df$TN / (df$FP+df$TN)
df$spe[is.na(df$spe)] = 0

256
## MCC
Carine Rey's avatar
Carine Rey committed
257 258 259
n_sites = sum(df[1,c("TP","FN")])
p = 140 / n_sites
n = (6000 -140)/ n_sites
Carine Rey's avatar
Carine Rey committed
260
df$FP_2 = df$FP * n
261
df$TP_2 = df$TP * p
Carine Rey's avatar
Carine Rey committed
262
df$FN_2 = df$FN * p
263
df$TN_2 = df$TN * n
Carine Rey's avatar
Carine Rey committed
264 265 266 267

df$mcc = (df$TP_2 * df$TN_2 - df$FP_2 * df$FN_2) / sqrt((df$TP_2+df$FP_2)*(df$TP_2+df$FN_2)*(df$TN_2+df$FP_2)*(df$TN_2+df$FN_2))
df$mcc[is.na(df$mcc)] = 0

268 269

## Precision
270 271
df$precision98_02 = (df$sens * 0.02) / (df$sens * 0.02 + (1-df$spe) * 0.98)
df$precision50_50 = (df$sens * 0.5) / (df$sens * 0.5 + (1-df$spe) * 0.5)
272 273 274 275 276

print(head(df))
print(tail(df))


Carine Rey's avatar
Carine Rey committed
277
df_out = df[,c("methode", "t", "sens", "spe", "mcc", "precision98_02","precision50_50","couple")]
278
colnames(df_out) = c("methode",  "threshold", "sensitivity", "specificity", "MCC", "precision98_02","precision50_50","couple")
279

Carine Rey's avatar
Carine Rey committed
280 281
df_out = df_out[order(df_out$methode),]

282 283
print(summary(df_out))

Carine Rey's avatar
Carine Rey committed
284
########################################################################
285 286 287
print("prep plot max mcc")

df_max_mcc_per_method = do.call(rbind, lapply(split(df_out,paste0(df_out$methode,df_out$couple)),
Carine Rey's avatar
Carine Rey committed
288
                                            function(x) {
289
                                                return(x[which.max(x$MCC),c("couple","methode", "threshold", "MCC","sensitivity","specificity","precision98_02","precision50_50")])
Carine Rey's avatar
Carine Rey committed
290
                                                }))
Carine Rey's avatar
Carine Rey committed
291 292 293 294 295 296 297

print(df_max_mcc_per_method)


df_max_mcc_per_method_2 = df_max_mcc_per_method
df_max_mcc_per_method_2$variable="MCC"

Carine Rey's avatar
Carine Rey committed
298
########################################################################
299 300
print("prep plot recall_precision_per_meth")

Carine Rey's avatar
Carine Rey committed
301 302

df_recall_sup09_per_meth = do.call(rbind, lapply(split(df_out,paste0(df_out$methode," ", df_out$couple)),
303
                                            function(x) {
304 305
                                                x$precision98_02[is.na(x$precision98_02)] = 0
                                                x2 = x[x$precision98_02 > 0.9,]
Carine Rey's avatar
Carine Rey committed
306
                                                if (nrow(x2) > 0) {
307
                                                    res = x2[which.max(x2$sensitivity),c("couple","methode", "threshold", "MCC","sensitivity","specificity","precision98_02","precision50_50")]
308
                                                } else {
309 310
                                                   x = x[1,c("couple","methode", "threshold", "MCC","sensitivity","specificity","precision98_02","precision50_50")]
                                                   x[,c("threshold", "MCC","sensitivity","specificity","precision98_02","precision50_50")] = NA
Carine Rey's avatar
Carine Rey committed
311
                                                   res = x
312
                                                }
Carine Rey's avatar
Carine Rey committed
313
                                                return(res)
314 315 316 317
                                            }))

print(df_recall_sup09_per_meth)

Carine Rey's avatar
Carine Rey committed
318
########################################################################
319 320


Carine Rey's avatar
Carine Rey committed
321
plot_out = function(df_out, df_d , df_recall_sup09_per_meth, meths = NULL, suffix="") {
322

Carine Rey's avatar
Carine Rey committed
323
    nb_c = length(unique(df_out$couple))
Carine Rey's avatar
Carine Rey committed
324 325
    colors = c(c("#984EA3","#4AA947","#377EB8","#E41A1C","#F5BE5B","#90EE90","#8B6914","#BFBFBF","#0000FF","#FFC0CB")[1:nb_c], c("#7F7F7F","#ADD8E6"))
    colors2 = c("#984EA3","#4AA947","#377EB8","#E41A1C","#F5BE5B","#90EE90","#8B6914","#BFBFBF","#7F7F7F","#ADD8E6")
326

Carine Rey's avatar
Carine Rey committed
327 328 329 330
    if (! is.null(meths)) {
    df_out = df_out[df_out$methode %in%meths,]
    df_d   = df_d[df_d$methode %in%meths,]
    df_recall_sup09_per_meth = df_recall_sup09_per_meth[df_recall_sup09_per_meth$methode %in%meths,]
Carine Rey's avatar
Carine Rey committed
331
    meths = meths[meths %in% df_out$methode]
Carine Rey's avatar
Carine Rey committed
332 333 334
    } else {
    meths = sort(unique(df_out$methode))
    }
335

Carine Rey's avatar
Carine Rey committed
336
    alpha = 0.7
Carine Rey's avatar
Carine Rey committed
337 338 339 340 341 342 343


    df_out$methode =  factor(df_out$methode, levels = meths)
    df_recall_sup09_per_meth$methode =  factor(df_recall_sup09_per_meth$methode, levels = meths)
    df_d$methode =  factor(df_d$methode, levels = meths)

    couple_levels = unique(df_out$couple)
344
    couple_levels = c(couple_levels[! couple_levels %in% c("H0/HaPCOC","H0/HaPCOC NeG1")] , "H0/HaPCOC", "H0/HaPCOC NeG1")
Carine Rey's avatar
Carine Rey committed
345 346 347 348 349 350 351

    df_out$couple =  factor(df_out$couple, levels = couple_levels)
    df_d$couple =  factor(df_d$couple, levels = couple_levels)
    df_recall_sup09_per_meth$couple =  factor(df_recall_sup09_per_meth$couple, levels = couple_levels )

    df_out_melt = melt(df_out, id.vars = c("couple", "methode","threshold"))
    df_out_melt$methode =  factor(df_out_melt$methode, levels = meths)
352
    df_out_melt$variable =  factor(df_out_melt$variable, levels = c("sensitivity", "specificity", "precision98_02","precision50_50", "MCC"))
Carine Rey's avatar
Carine Rey committed
353 354 355 356
    df_out_melt$couple =  factor(df_out_melt$couple, levels = couple_levels)

    print("plot recall_precision_per_meth")

357
    plot = ggplot(df_out, aes(x=sensitivity, y=precision98_02, col = couple))
Carine Rey's avatar
Carine Rey committed
358
    plot = plot + theme_bw()
359
    plot = plot + labs(x="Sensitivity (= Recall)", y="Precision (98/2)")
Carine Rey's avatar
Carine Rey committed
360 361 362 363
    plot = plot + theme(legend.position="top")
    plot = plot + ylim(c(0,1)) + xlim(c(0,1))
    plot = plot + guides(fill=FALSE) + guides(col=FALSE)
    plot = plot + scale_color_manual(values=colors)
Carine Rey's avatar
Carine Rey committed
364 365 366
    #plot = plot + geom_point(size=1, alpha=alpha)
    plot = plot + geom_step(direction="vh", size=1, alpha=alpha)
    plot = plot + geom_hline( aes(yintercept = 0.9), col="black" ,  size = 0.5, show.legend = NA,linetype="dashed")
Carine Rey's avatar
Carine Rey committed
367 368
    plot = plot + facet_grid(couple ~ methode)
    plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))
Carine Rey's avatar
Carine Rey committed
369
    plot = plot + ggtitle(paste0("Tree: ", opt$tree_prefix , "\n Date: " , date ))
Carine Rey's avatar
Carine Rey committed
370 371 372 373 374 375 376 377 378 379
    plot_recall_precision = plot


    save_plot(paste0(opt$out,suffix,".recall_precision_per_meth.pdf"),
              plot_recall_precision,
              ncol = 0.4* length(unique(df_out_melt$methode)),
              nrow = 0.4* length(unique(df_out_melt$couple)),
              base_aspect_ratio = 1,
              limitsize = FALSE
              )
Carine Rey's avatar
Carine Rey committed
380
    print("plot recall_precision")
Carine Rey's avatar
Carine Rey committed
381

Carine Rey's avatar
Carine Rey committed
382
    if (length(meths) <= 9) {
383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
        plot = ggplot(df_out, 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")
        plot = plot + ylim(c(0,1)) + xlim(c(0,1))
        plot = plot + guides(fill=FALSE)
        plot = plot + scale_color_manual(values=colors2)
        #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 + 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))
                                )
Carine Rey's avatar
Carine Rey committed
404

405 406 407 408 409
        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)
Carine Rey's avatar
Carine Rey committed
410

411 412 413 414 415 416 417


        save_plot(paste0(opt$out,suffix,".recall_precision.pdf"),
                  plot_recall_precision_papier,
                  ncol = 0.6* length(unique(df_out_melt$couple)),
                  nrow = 1,
                  base_aspect_ratio = 1,
418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455
                  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 = plot + theme_bw()
        plot = plot + labs(x="Sensitivity (= Recall)", y="Precision (98/2)")
        plot = plot + theme(legend.position="top")
        plot = plot + ylim(c(0,1)) + xlim(c(0,1))
        plot = plot + guides(fill=FALSE)
        plot = plot + scale_color_manual(values=colors2)
        #plot = plot + geom_point(size=1, alpha=alpha)
        plot = plot + geom_step(direction="vh", size=0.7, 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 + 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))
                                )

        plot_recall_precision_papier_ok = 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_ok.pdf"),
                  plot_recall_precision_papier_ok,
                  ncol = 0.6* 3,
                  nrow = 1,
                  base_aspect_ratio = 1,
456 457
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
458
    } else {
459
        plot_recall_precision_papier = NULL
Carine Rey's avatar
Carine Rey committed
460
    }
Carine Rey's avatar
Carine Rey committed
461

Carine Rey's avatar
Carine Rey committed
462
    print("plot FPR")
463
    #dcast
Carine Rey's avatar
Carine Rey committed
464 465 466 467 468 469
    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$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
Carine Rey's avatar
Carine Rey committed
470

Carine Rey's avatar
Carine Rey committed
471
    if (length(df_out_melt2$methode) > 0){
472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497
        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)

        plot = ggplot(df_out_melt2, aes(x=threshold, y=FPR, col=couple))
        plot = plot + theme_bw()
        plot = plot + theme(legend.position="top",
                            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)
        plot = plot + ylim(c(0,1)) + xlim(c(0,1))
        plot = plot + geom_step(size=1, alpha=alpha)
        plot = plot + scale_color_manual(values=colors)
        plot = plot + geom_vline( data = df_recall_sup09_per_meth2 , aes(xintercept = threshold, col=couple), size = 0.5, show.legend = NA,linetype="dashed")
        plot = plot + facet_grid(. ~ methode, scales="free")
        plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))

        plot_FPR = plot

        save_plot(paste0(opt$out,suffix,".FPR.pdf"),
                  plot_FPR,
                  ncol = 0.4 * length(unique(df_out_melt$methode)),
                  nrow = 1,
                  base_aspect_ratio = 1.5,
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
498

499 500 501 502 503
        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)
Carine Rey's avatar
Carine Rey committed
504

505
            # complete values for not existing threshold between couple
Carine Rey's avatar
Carine Rey committed
506

507 508 509 510 511
            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")
Carine Rey's avatar
Carine Rey committed
512

513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538
                fix_val = function(r,df_no_na) {
                                           x = r["value"]
                                           t = r["threshold"]
                                           if (is.na(x)) {
                                               t_n  = df_no_na$t [ df_no_na$t < t ]
                                               t_n = t_n [length(t_n)]
                                               x = df_no_na$value[df_no_na$t == t_n]
                                           }
                                           return(x)
                                       }
                df_ok = do.call(rbind, lapply(split(df,df$methode),
                                            function(df_meth) {
                                                df_no_na = df_meth[! is.na(df_meth$value),]
                                                df_meth$value = apply(df_meth[,c(2,3)],1,fix_val,df_no_na=df_no_na)
                                                return(df_meth)
                                                }))
                return(df_ok[,3])
            }

            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$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
Carine Rey's avatar
Carine Rey committed
539

540
            print(df_out_melt2_dcast)
Carine Rey's avatar
Carine Rey committed
541

542 543
            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))
Carine Rey's avatar
Carine Rey committed
544

545 546 547
            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"))
Carine Rey's avatar
Carine Rey committed
548

549 550
            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"))
Carine Rey's avatar
Carine Rey committed
551

552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572

            x_labs = "Threshold"
            y_labs = "Delta FP Rate"
            plot = ggplot()
            plot = plot + theme_bw()
            plot = plot + theme(legend.position="top",
                                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)
            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 + facet_grid(. ~ methode, scales="free")
            plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))

            plot_FPR2 = plot
        } else {
            plot_FPR2 = NULL
        }
Carine Rey's avatar
Carine Rey committed
573 574
    } else {
    plot_FPR = NULL
575
    plot_FPR2 = NULL
Carine Rey's avatar
Carine Rey committed
576
    }
Carine Rey's avatar
Carine Rey committed
577 578


Carine Rey's avatar
Carine Rey committed
579 580 581 582 583 584
    print("plot per indicator")

    x_labs = "Threshold"
    y_labs = ""
    plot = ggplot(df_out_melt, aes(x=threshold, y=value, col=couple))
    plot = plot + theme_bw()
Carine Rey's avatar
Carine Rey committed
585
    plot = plot + theme(legend.position="top",
Carine Rey's avatar
Carine Rey committed
586 587
    legend.background = element_rect(fill="white", size=0.5, linetype="solid", colour ="black")) +
    guides(colour = guide_legend(override.aes = list(alpha = 1), nrow=2))
Carine Rey's avatar
Carine Rey committed
588

Carine Rey's avatar
Carine Rey committed
589 590
    plot = plot + guides(fill=FALSE) + theme(legend.position="top")
    plot = plot + labs(x=x_labs, y=y_labs) #+ ylim(y_lim)
Carine Rey's avatar
Carine Rey committed
591 592
    #plot = plot + geom_point(size=0.5, alpha=alpha)
    plot = plot + geom_step(size=1, alpha=alpha)
Carine Rey's avatar
Carine Rey committed
593 594
    plot = plot + scale_color_manual(values=colors)
    #plot = plot + geom_hline( data = df_max_mcc_per_method_2, aes(yintercept = MCC), alpha=alpha, show.legend = NA)
Carine Rey's avatar
Carine Rey committed
595
    plot = plot + geom_vline( data = df_recall_sup09_per_meth, aes(xintercept = threshold, col=couple), size = 0.5, show.legend = NA,linetype="dashed")
Carine Rey's avatar
Carine Rey committed
596 597 598 599 600
    plot = plot + facet_grid(variable ~ methode, scales="free")
    plot = plot + theme(axis.text.x = element_text(angle = 45, hjust = 1))

    plot_max_MCC = plot

Carine Rey's avatar
Carine Rey committed
601
    save_plot(paste0(opt$out,suffix,".indicator_per_meth.pdf"),
Carine Rey's avatar
Carine Rey committed
602 603 604 605 606 607
    plot_max_MCC,
    ncol = 0.4 * length(unique(df_out_melt$methode)),
    nrow = 1.7,
    base_aspect_ratio = 1.5,
    limitsize = FALSE
    )
Carine Rey's avatar
Carine Rey committed
608

Carine Rey's avatar
Carine Rey committed
609
    plot_value_JSD = NULL
Carine Rey's avatar
Carine Rey committed
610
    if ( all( c("P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2") %in% colnames(df_d))){
Carine Rey's avatar
Carine Rey committed
611

Carine Rey's avatar
Carine Rey committed
612
        print("plot value/distance")
Carine Rey's avatar
Carine Rey committed
613 614


Carine Rey's avatar
Carine Rey committed
615

Carine Rey's avatar
Carine Rey committed
616 617 618 619 620 621 622 623 624 625 626 627 628
        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
Carine Rey's avatar
Carine Rey committed
629

Carine Rey's avatar
Carine Rey committed
630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683
        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
Carine Rey's avatar
Carine Rey committed
684

Carine Rey's avatar
Carine Rey committed
685 686 687 688 689 690 691 692

        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
                  )
Carine Rey's avatar
Carine Rey committed
693

Carine Rey's avatar
Carine Rey committed
694

Carine Rey's avatar
Carine Rey committed
695
    }
Carine Rey's avatar
Carine Rey committed
696 697 698 699 700

    #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,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,length(unique(df_out$couple))*0.8),
    #                 nrow=7
Carine Rey's avatar
Carine Rey committed
701
    #                 )
Carine Rey's avatar
Carine Rey committed
702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719
    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),
                    nrow=5)

    if (is.null(plot_value_JSD)) {
        plot_png = plot_grid(plot_recall_precision,plot_max_MCC,
                    labels = c("A", "B"),
                    rel_heights = c(1,1),
                    nrow=2
                     )
    } else {
        plot_png = plot_grid(plot_recall_precision,plot_max_MCC,plot_value_JSD,
                    labels = c("A", "B","C"),
                    rel_heights = c(0.5,0.5,1),
                    nrow=3
                     )
    }
Carine Rey's avatar
Carine Rey committed
720 721 722 723

    save_plot(paste0(opt$out,suffix,".pdf"),
              plot,
              ncol = 0.4* length(unique(df_out_melt$methode)),
Carine Rey's avatar
Carine Rey committed
724 725 726 727 728 729 730 731
              nrow = length(unique(df_out$couple)) * 0.5  + 2 + 1 +1 +1,
              base_aspect_ratio = 1,
              limitsize = FALSE
              )
    save_plot(paste0(opt$out,suffix,".png"),
              plot_png,
              ncol = 0.4* length(unique(df_out_melt$methode)),
              nrow = length(unique(df_out$couple)) * 1.5,
Carine Rey's avatar
Carine Rey committed
732 733 734 735
              base_aspect_ratio = 1,
              limitsize = FALSE
              )
}
736

Carine Rey's avatar
Carine Rey committed
737 738


Carine Rey's avatar
Carine Rey committed
739
plot_out(df_out, df_d, df_recall_sup09_per_meth, suffix="")
Carine Rey's avatar
Carine Rey committed
740
condensed_meths = c("PCOC","Diffsel_mean","Identical_LG08","Mutinomial_1MinusLRT","Tdg09_1MinusLRT","Msd_1MinusP","Msd_0.05_1MinusP","Topological_LG08")
Carine Rey's avatar
Carine Rey committed
741
plot_out(df_out, df_d, df_recall_sup09_per_meth, meths=condensed_meths,  suffix="_condensed")
742 743


Carine Rey's avatar
Carine Rey committed
744
write.table(df_d, file=paste0(opt$out,".complete_d.tsv"), row.names=FALSE, quote=F, sep = "\t")
Carine Rey's avatar
Carine Rey committed
745
write.table(df_out, file=paste0(opt$out,".complete.tsv"), row.names=FALSE, quote=F, sep = "\t")
746 747
write.table(df_max_mcc_per_method, file=paste0(opt$out,".max_MCC_per_meth.tsv"), row.names=FALSE, quote=F, sep = "\t")
write.table(df_recall_sup09_per_meth, file=paste0(opt$out,".recall09_per_meth.tsv"), row.names=FALSE, quote=F, sep = "\t")