calc_t_per_meth.R 41 KB
Newer Older
Carine Rey's avatar
Carine Rey committed
1 2 3 4 5 6 7 8 9
#!/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")
10
library("flux")
Carine Rey's avatar
Carine Rey committed
11

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

Carine Rey's avatar
Carine Rey committed
14
option_list = list(
15 16 17 18 19 20 21 22 23 24
  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"),
25 26 27 28
  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"),
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 35 36 37 38
  make_option(c("--HaPC_NeG4_NeC_1_res"       ) , type="character", default=NA, help="merged_results HaPC_NeG4_NeC_1_res"       , metavar="character"),
  make_option(c("--H0_NeG4_NeC_1_res"         ) , type="character", default=NA, help="merged_results H0_NeG4_NeC_1_res"         , metavar="character"),
  make_option(c("--HaPC_NeG1_NeC_4_res"         ) , type="character", default=NA, help="merged_results HaPC_NeG1_NeC_4_res"         , metavar="character"),
  make_option(c("--H0_NeG1_NeC_4_res"           ) , type="character", default=NA, help="merged_results H0_NeG1_NeC_4_res"           , metavar="character"),
  
Carine Rey's avatar
Carine Rey committed
39 40
  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"),
41
  make_option(c("--HaPCOC"          ) , type="character", default=NA, help="merged_results HaPCOC"          , metavar="character"),
Carine Rey's avatar
Carine Rey committed
42

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

45 46

  make_option(c("-o","--out"), type="character", default="out",
Carine Rey's avatar
Carine Rey committed
47 48 49 50 51 52
              help="output prefix [default= %default]", metavar="character")
);

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

53 54
print(opt)

55
if (is.na(opt$H0_NeG1)){
Carine Rey's avatar
Carine Rey committed
56
  print_help(opt_parser)
57
  stop("At least one argument must be supplied (H0_NeG1 or H0 input file)", call.=FALSE)
Carine Rey's avatar
Carine Rey committed
58
}
59
if (is.null(opt$HaPCOC)){
Carine Rey's avatar
Carine Rey committed
60
  print_help(opt_parser)
61
  stop("At least one argument must be supplied (HaPCOC input file)", call.=FALSE)
Carine Rey's avatar
Carine Rey committed
62 63
}

Carine Rey's avatar
Carine Rey committed
64 65 66 67 68 69 70


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

print("prep input files")

read_hyp = function(opt_name) {
71 72
  print(opt_name)
  if (! is.na(opt_name)){
Carine Rey's avatar
Carine Rey committed
73 74
    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
75 76 77
    #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
78 79 80 81 82
  } else {
    return(NULL)
  }
}

Carine Rey's avatar
Carine Rey committed
83 84 85 86 87 88 89 90 91 92
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 )
93

Carine Rey's avatar
Carine Rey committed
94 95 96 97
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    )
98

99 100 101 102
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
103

104 105 106 107 108
df_HaPC_NeG4_NeC_1      = read_hyp(opt$HaPC_NeG4_NeC_1)
df_H0_NeG4_NeC_1      = read_hyp(opt$H0_NeG4_NeC_1  )
df_HaPC_NeG1_NeC_4        = read_hyp(opt$HaPC_NeG1_NeC_4  )
df_H0_NeG1_NeC_4        = read_hyp(opt$H0_NeG1_NeC_4    )

Carine Rey's avatar
Carine Rey committed
109 110
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
111

Carine Rey's avatar
Carine Rey committed
112
df_HaPCOC = read_hyp(opt$HaPCOC)
113

Carine Rey's avatar
Carine Rey committed
114

115 116 117



Carine Rey's avatar
Carine Rey committed
118 119 120 121 122 123
########################################################################

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
124 125
    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
126
    df_ha = melt(df_ha, id.vars = id_vars, variable.name = "methode")
Carine Rey's avatar
Carine Rey committed
127
    print(name)
Carine Rey's avatar
Carine Rey committed
128
    print(head(df_h0))
Carine Rey's avatar
Carine Rey committed
129 130
    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
131 132 133 134 135 136 137
    df_melt$couple = name
  } else {
    df_melt = NULL
  }
  return(df_melt)
}

Carine Rey's avatar
Carine Rey committed
138 139 140 141 142 143
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")
144

Carine Rey's avatar
Carine Rey committed
145 146
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")
147

148 149
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 NeG1_NeC_5")
Carine Rey's avatar
Carine Rey committed
150

151 152 153
df_d_H0HaPC_NeG4_NeC_1   = build_df_dist_couple(df_H0_NeG4_NeC_1, df_HaPC_NeG4_NeC_1, "H0/HaPC NeG4_NeC_1")
df_d_H0HaPC_NeG1_NeC_4   = build_df_dist_couple(df_H0_NeG1_NeC_4, df_HaPC_NeG1_NeC_4, "H0/HaPC NeG1_NeC_4")

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

Carine Rey's avatar
Carine Rey committed
156 157 158
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")

159 160 161
df_d_H0H0_NeG1NeG1_NeC4   = build_df_dist_couple(df_H0_NeG1, df_H0_NeG1_NeC_4, "H0/H0 NeG1/NeG1_NeC_4")
df_d_H0H0_NeG4NeG4_NeC1   = build_df_dist_couple(df_H0_NeG4, df_H0_NeG4_NeC_1, "H0/H0 NeG4/NeG4_NeC_1")

162
df_d = rbind.data.frame(
163 164 165 166 167
                      df_d_H0HaPCOC_NeG1,
                      df_d_H0HaPC_NeG1,
                      df_d_H0HaPC_NeG2,
                      df_d_H0HaPC_NeG3,
                      df_d_H0HaPC_NeG4,
168 169
                      df_d_H0HaPC_NeG5,
                      df_d_H0HaPC_NeG5_NeC_div2,
Carine Rey's avatar
Carine Rey committed
170
                      df_d_H0HaPC_NeG5_NeC_x2,
171 172
                      df_d_H0HaPC_NeG5_NeC_1,
                      df_d_H0HaPC_NeG1_NeC_5,
Carine Rey's avatar
Carine Rey committed
173 174
                      df_d_H0H0_NeG1NeG1_NeC5,
                      df_d_H0H0_NeG5NeG5_NeC1,
175 176
                      df_d_H0HaPC_NeG4_NeC_1,
                      df_d_H0HaPC_NeG1_NeC_4,
177 178
                      df_d_H0H0_NeG4NeG4_NeC1,
                      df_d_H0H0_NeG1NeG1_NeC4,
Carine Rey's avatar
Carine Rey committed
179 180
                      df_d_H0HaPC_NeG5_indel
                      )
Carine Rey's avatar
Carine Rey committed
181 182 183 184 185 186

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

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

187 188


Carine Rey's avatar
Carine Rey committed
189 190 191 192
########################################################################

print("prep df")

Carine Rey's avatar
Carine Rey committed
193 194
## fun...

Carine Rey's avatar
Carine Rey committed
195
calc_TN_FP = function(t, vals){
Carine Rey's avatar
Carine Rey committed
196 197
    TN = 0
    FP = 0
198
    vals = na.omit(vals)
Carine Rey's avatar
Carine Rey committed
199
    if (length(vals) > 0) {
Carine Rey's avatar
Carine Rey committed
200 201
        TN = sum(vals < t)
        FP = sum(vals >= t)
Carine Rey's avatar
Carine Rey committed
202
    }
Carine Rey's avatar
Carine Rey committed
203
    return(data.frame(t=t, TN = TN, FP = FP))
Carine Rey's avatar
Carine Rey committed
204 205
}

Carine Rey's avatar
Carine Rey committed
206
calc_TP_FN = function(t, vals){
Carine Rey's avatar
Carine Rey committed
207 208 209 210
    FN = 0
    TP = 0
    vals = na.omit(vals)
    if (length(vals) > 0) {
Carine Rey's avatar
Carine Rey committed
211 212
        FN = sum(vals < t)
        TP = sum(vals >= t)
Carine Rey's avatar
Carine Rey committed
213
    }
Carine Rey's avatar
Carine Rey committed
214
    return(data.frame(t=t, TP = TP, FN = FN))
Carine Rey's avatar
Carine Rey committed
215 216
}

Carine Rey's avatar
Carine Rey committed
217 218 219 220
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
221

Carine Rey's avatar
Carine Rey committed
222 223 224 225 226
    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
227 228 229 230 231

    return(df_TN_FP_TP_FN)
}


232
build_df_couple = function (df_h0,df_ha,name) {
233
  print(name)
234
  if ((! is.null(df_h0)) & (! is.null(df_ha))) {
Carine Rey's avatar
Carine Rey committed
235
    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
236 237
    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))
238
    df$couple = name
239
    print(head(df))
240 241
  } else {
    df = NULL
242
    print(df)
243
  }
Carine Rey's avatar
Carine Rey committed
244

245 246
  return(df)
}
Carine Rey's avatar
Carine Rey committed
247

Carine Rey's avatar
Carine Rey committed
248 249 250 251 252 253
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")
254 255


Carine Rey's avatar
Carine Rey committed
256 257
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")
258

259 260
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
261

Carine Rey's avatar
Carine Rey committed
262 263 264
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")

265 266 267 268 269 270
df_H0HaPC_NeG4_NeC_1   = build_df_couple(df_H0_NeG4_NeC_1, df_HaPC_NeG4_NeC_1, "H0/HaPC NeG4_NeC_1")
df_H0HaPC_NeG1_NeC_4   = build_df_couple(df_H0_NeG1_NeC_4, df_HaPC_NeG1_NeC_4, "H0/HaPC NeG1_NeC_4")

df_H0H0_NeG4NeG4_NeC_1   = build_df_couple(df_H0_NeG4, df_H0_NeG4_NeC_1, "H0/H0 NeG4/NeG4_NeC_1")
df_H0H0_NeG1NeG1_NeC_4   = build_df_couple(df_H0_NeG1, df_H0_NeG1_NeC_4, "H0/H0 NeG1/NeG1_NeC_4")

Carine Rey's avatar
Carine Rey committed
271 272 273
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
274 275
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,
276
            df_H0HaPC_NeG5_NeC_1,df_H0HaPC_NeG1_NeC_5,
Carine Rey's avatar
Carine Rey committed
277
            df_H0H0_NeG5NeG5_NeC_1, df_H0H0_NeG1NeG1_NeC_5,
278 279
            df_H0HaPC_NeG4_NeC_1,df_H0HaPC_NeG1_NeC_4,
            df_H0H0_NeG4NeG4_NeC_1, df_H0H0_NeG1NeG1_NeC_4,
Carine Rey's avatar
Carine Rey committed
280
            df_H0HaPC_NeG5_NeC_x2,df_H0HaPC_NeG5_indel)
281 282 283

df_l = df_l[-which(sapply(df_l, is.null))]
df =  do.call("rbind",df_l)
284 285 286 287

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

Carine Rey's avatar
Carine Rey committed
288 289

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

293
## Specificity
Carine Rey's avatar
Carine Rey committed
294 295 296
df$spe =  df$TN / (df$FP+df$TN)
df$spe[is.na(df$spe)] = 0

297
## MCC
Carine Rey's avatar
Carine Rey committed
298 299 300
n_sites = sum(df[1,c("TP","FN")])
p = 140 / n_sites
n = (6000 -140)/ n_sites
Carine Rey's avatar
Carine Rey committed
301
df$FP_2 = df$FP * n
302
df$TP_2 = df$TP * p
Carine Rey's avatar
Carine Rey committed
303
df$FN_2 = df$FN * p
304
df$TN_2 = df$TN * n
Carine Rey's avatar
Carine Rey committed
305 306 307 308

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

309 310

## Precision
311 312
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)
313 314 315 316 317

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


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

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

323 324
print(summary(df_out))

Carine Rey's avatar
Carine Rey committed
325
########################################################################
326 327 328
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
329
                                            function(x) {
330
                                                return(x[which.max(x$MCC),c("couple","methode", "threshold", "MCC","sensitivity","specificity","precision98_02","precision50_50")])
Carine Rey's avatar
Carine Rey committed
331
                                                }))
Carine Rey's avatar
Carine Rey committed
332 333 334 335 336 337 338

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
339
########################################################################
340 341
print("prep plot recall_precision_per_meth")

Carine Rey's avatar
Carine Rey committed
342 343

df_recall_sup09_per_meth = do.call(rbind, lapply(split(df_out,paste0(df_out$methode," ", df_out$couple)),
344
                                            function(x) {
345 346
                                                x$precision98_02[is.na(x$precision98_02)] = 0
                                                x2 = x[x$precision98_02 > 0.9,]
Carine Rey's avatar
Carine Rey committed
347
                                                if (nrow(x2) > 0) {
348
                                                    res = x2[which.max(x2$sensitivity),c("couple","methode", "threshold", "MCC","sensitivity","specificity","precision98_02","precision50_50")]
349
                                                } else {
350 351
                                                   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
352
                                                   res = x
353
                                                }
Carine Rey's avatar
Carine Rey committed
354
                                                return(res)
355 356 357 358
                                            }))

print(df_recall_sup09_per_meth)

359 360 361 362 363 364 365 366 367 368

print("calc auc")

df_auc = do.call(rbind, lapply(split(df_out,paste0(df_out$methode,df_out$couple)),
                                            function(x) {
                                                print(x)
                                                R = x$sensitivity
                                                P = x$precision98_02
                                                R = R[! is.na(P)]
                                                P = P[! is.na(P)]
369 370 371 372
                                                auc = 0
                                                if (length(P) > 1) {
                                                auc = auc(x = R, y = P, dens = 1)
                                                }
373 374 375 376 377 378 379
                                                df= data.frame(methode =x$methode[1] ,couple =x$couple[1] ,auc=auc)
                                                return(df)
                                                }))

df_auc = df_auc[order(df_auc$couple),]
print(df_auc)

Carine Rey's avatar
Carine Rey committed
380
########################################################################
381 382


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

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

Carine Rey's avatar
Carine Rey committed
389 390 391 392
    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
393
    meths = meths[meths %in% df_out$methode]
Carine Rey's avatar
Carine Rey committed
394 395 396
    } else {
    meths = sort(unique(df_out$methode))
    }
397

Carine Rey's avatar
Carine Rey committed
398
    alpha = 0.7
Carine Rey's avatar
Carine Rey committed
399 400 401 402 403 404 405


    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)
406
    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
407 408 409 410 411 412 413

    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)
414
    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
415 416 417 418
    df_out_melt$couple =  factor(df_out_melt$couple, levels = couple_levels)

    print("plot recall_precision_per_meth")

419
    plot = ggplot(df_out, aes(x=sensitivity, y=precision98_02, col = couple))
Carine Rey's avatar
Carine Rey committed
420
    plot = plot + theme_bw()
421
    plot = plot + labs(x="Sensitivity (= Recall)", y="Precision (98/2)")
Carine Rey's avatar
Carine Rey committed
422 423 424 425
    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
426 427 428
    #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
429 430
    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
431
    plot = plot + ggtitle(paste0("Tree: ", opt$tree_prefix , "\n Date: " , date ))
Carine Rey's avatar
Carine Rey committed
432 433 434 435 436 437 438 439 440 441
    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
442
    print("plot recall_precision")
Carine Rey's avatar
Carine Rey committed
443

Carine Rey's avatar
Carine Rey committed
444
    if (length(meths) <= 9) {
445
        couple_i = c("H0/HaPC NeG1","H0/HaPC NeG2","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","H0/H0 NeG5/NeG4_NeC_1","H0/H0 NeG3/NeG1_NeC_4", "H0/HaPC NeG4_NeC_1","H0/HaPC NeG1_NeC_4")
Carine Rey's avatar
Carine Rey committed
446
        plot = ggplot(subset(df_out, couple %in% couple_i), aes(x=sensitivity, y=precision98_02, col = methode))
447 448 449 450 451 452 453 454 455
        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")
Carine Rey's avatar
Carine Rey committed
456
        plot = plot + facet_wrap(~ couple, ncol=4)
457 458 459 460 461 462 463 464 465
        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),
                                            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
466

467 468 469 470 471
        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
472

473 474 475 476 477 478 479


        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,
480 481
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
482 483
        
        plot_PR_c = function(couple_l) {
484 485 486 487 488 489
        tmp_df = subset(df_out,couple%in%couple_l)
        if ("NULL" %in% couple_l) {
            tmp_df[dim(tmp_df)[1]+1,] = tmp_df[dim(tmp_df)[1],]
            tmp_df$couple[dim(tmp_df)[1]] = "to_be_rm"
        }
        plot = ggplot(tmp_df, aes(x=sensitivity, y=precision98_02, col = methode))
490 491 492 493 494 495 496 497 498 499 500
        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))
Carine Rey's avatar
Carine Rey committed
501 502
        plot}
        
503
        plot = plot_PR_c(c("H0/HaPC NeG1","H0/HaPC NeG2","H0/HaPC NeG4","H0/H0 NeG4/NeG4_NeC_1","H0/H0 NeG1/NeG1_NeC_4", "H0/HaPC NeG4_NeC_1","H0/HaPC NeG1_NeC_4"))
504 505 506
        plot1 = plot_PR_c(c("H0/HaPC NeG1","H0/HaPC NeG2","H0/HaPC NeG4"))
        plot2 = plot_PR_c(c("H0/H0 NeG4/NeG4_NeC_1","H0/H0 NeG1/NeG1_NeC_4","NULL"))
        plot3 = plot_PR_c(c("H0/HaPC NeG4_NeC_1","H0/HaPC NeG1_NeC_4","NULL"))
507 508 509 510 511 512 513 514 515

        legend_PR = get_legend(plot + theme(legend.position="top",
                                            legend.text  = element_text(size=10),
                                            legend.title = element_text(size=0),
                                            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
516 517 518 519 520
        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,
521
                  labels = c("","Only Pressure Shift","Only Change in Selection Efficacy","Pressure Shift & Change in Selection Efficacy"),
522
                  align = "h",
Carine Rey's avatar
Carine Rey committed
523 524 525
                  
                  rel_heights = c( 0.3, 0.7,0.7,0.7),
                  hjust = 0, vjust = 0)
526 527 528 529 530



        save_plot(paste0(opt$out,suffix,".recall_precision_ok.pdf"),
                  plot_recall_precision_papier_ok,
531 532
                  ncol = 2.5,
                  nrow = 2.5,
533
                  base_aspect_ratio = 1,
534 535
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
536
    } else {
537
        plot_recall_precision_papier = NULL
Carine Rey's avatar
Carine Rey committed
538
    }
Carine Rey's avatar
Carine Rey committed
539

Carine Rey's avatar
Carine Rey committed
540
    print("plot FPR")
541
    #dcast
Carine Rey's avatar
Carine Rey committed
542 543
    x_labs = "Threshold"
    y_labs = "FP Rate (1 -Spe)"
Carine Rey's avatar
Carine Rey committed
544
    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"))
Carine Rey's avatar
Carine Rey committed
545 546 547
    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
548

Carine Rey's avatar
Carine Rey committed
549
    if (length(df_out_melt2$methode) > 0){
Carine Rey's avatar
Carine Rey committed
550
        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"))
551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575
        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
576

Carine Rey's avatar
Carine Rey committed
577
        if (length(unique(df_out_melt2$couple)) == 4){
578 579 580
            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))
Carine Rey's avatar
Carine Rey committed
581

582
            # complete values for not existing threshold between couple
Carine Rey's avatar
Carine Rey committed
583

584 585 586 587 588
            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
589

590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609
                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")
Carine Rey's avatar
Carine Rey committed
610 611 612 613 614 615 616 617
            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")

618

Carine Rey's avatar
Carine Rey committed
619 620
            #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
621

Carine Rey's avatar
Carine Rey committed
622 623
            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
624

Carine Rey's avatar
Carine Rey committed
625 626
            #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"))
627
            print(head(df_out_melt2_dcast_melt))
Carine Rey's avatar
Carine Rey committed
628

Carine Rey's avatar
Carine Rey committed
629 630
            #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"))
631

Carine Rey's avatar
Carine Rey committed
632
            df_recall_sup09_per_meth3 = subset(df_recall_sup09_per_meth2, couple %in%c("H0NeG5","H0NeG1"))
Carine Rey's avatar
Carine Rey committed
633

634
            df_recall_sup09_per_meth3$variable = df_recall_sup09_per_meth3$couple
Carine Rey's avatar
Carine Rey committed
635
            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"))
Carine Rey's avatar
Carine Rey committed
636

637 638 639 640 641 642 643 644 645 646 647 648 649

            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)
Carine Rey's avatar
Carine Rey committed
650
            plot = plot + scale_color_manual(values=colors[c(2,4,1,3)])
651 652 653 654 655 656 657
            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
658 659
    } else {
    plot_FPR = NULL
660
    plot_FPR2 = NULL
Carine Rey's avatar
Carine Rey committed
661
    }
Carine Rey's avatar
Carine Rey committed
662 663


Carine Rey's avatar
Carine Rey committed
664 665 666 667 668 669
    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
670
    plot = plot + theme(legend.position="top",
Carine Rey's avatar
Carine Rey committed
671 672
    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
673

Carine Rey's avatar
Carine Rey committed
674 675
    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
676 677
    #plot = plot + geom_point(size=0.5, alpha=alpha)
    plot = plot + geom_step(size=1, alpha=alpha)
Carine Rey's avatar
Carine Rey committed
678 679
    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
680
    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
681 682 683 684 685
    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
686
    save_plot(paste0(opt$out,suffix,".indicator_per_meth.pdf"),
Carine Rey's avatar
Carine Rey committed
687 688 689 690 691 692
    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
693

Carine Rey's avatar
Carine Rey committed
694
    plot_value_JSD = NULL
Carine Rey's avatar
Carine Rey committed
695
    if ( all( c("P_JSD","P_ED", "C1", "C2", "entropy_C1","entropy_C2") %in% colnames(df_d)) & (length(meths) <= 9) ){
Carine Rey's avatar
Carine Rey committed
696

Carine Rey's avatar
Carine Rey committed
697
        print("plot value/distance")
Carine Rey's avatar
Carine Rey committed
698 699


Carine Rey's avatar
Carine Rey committed
700

Carine Rey's avatar
Carine Rey committed
701
        plot = ggplot(subset(df_d, val_H0Ha == "value_Ha"), aes(y=P_JSD, x=value, shape =  val_H0Ha, col = val_H0Ha))
Carine Rey's avatar
Carine Rey committed
702 703 704 705 706 707 708 709 710 711 712 713
        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
714

Carine Rey's avatar
Carine Rey committed
715 716 717 718 719

        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))
Carine Rey's avatar
Carine Rey committed
720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772
        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
773

Carine Rey's avatar
Carine Rey committed
774 775 776 777 778 779 780 781

        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
782

Carine Rey's avatar
Carine Rey committed
783

Carine Rey's avatar
Carine Rey committed
784
    }
Carine Rey's avatar
Carine Rey committed
785 786 787 788 789

    #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
790
    #                 )
Carine Rey's avatar
Carine Rey committed
791 792
    plot = plot_grid(plot_recall_precision,plot_recall_precision_papier,plot_max_MCC,plot_FPR,plot_FPR2,
                    labels = c("A", "A'", "B","B'","B''"),
Carine Rey's avatar
Carine Rey committed
793
                    rel_heights = c(length(unique(df_out$couple))*0.8,4,4,2,2),
Carine Rey's avatar
Carine Rey committed
794 795 796 797 798 799 800 801 802
                    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 {
Carine Rey's avatar
Carine Rey committed
803 804 805 806
        plot_png = plot_grid(plot_recall_precision,plot_max_MCC,plot_value_JSD,plot_value_ED,
                    labels = c("A", "B","C","D"),
                    rel_heights = c(0.5,0.5,1,1),
                    nrow=4
Carine Rey's avatar
Carine Rey committed
807 808
                     )
    }
Carine Rey's avatar
Carine Rey committed
809 810 811

    save_plot(paste0(opt$out,suffix,".pdf"),
              plot,
Carine Rey's avatar
Carine Rey committed
812
              ncol = 0.5* length(unique(df_out_melt$methode)),
Carine Rey's avatar
Carine Rey committed
813 814 815 816 817 818 819 820
              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
821 822 823 824
              base_aspect_ratio = 1,
              limitsize = FALSE
              )
}
825

Carine Rey's avatar
Carine Rey committed
826 827