calc_t_per_meth.R 43.3 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)]
Carine Rey's avatar
Carine Rey committed
369 370 371
                                                P = P[order(R)]    
                                                R = R[order(R)]
                                                
372
                                                auc = 0
Carine Rey's avatar
Carine Rey committed
373 374 375
                                                auc2 = 0
                                                best_R_90_P = NA
                                               
376
                                                if (length(P) > 1) {
Carine Rey's avatar
Carine Rey committed
377 378 379 380 381 382 383 384 385 386 387 388 389
                                                    # add point for step curves
                                                    if (length(P) < 20) {
                                                    P_tmp = P    
                                                    R_tmp = R
                                                    for (i in seq(1,length(P)-1)) {
                                                        R_i = R[i] + (R[i+1] - R[i]) * 0.001
                                                        P_i = P[i+1]
                                                        P_tmp = c(P_tmp, P_i)
                                                        R_tmp = c(R_tmp, R_i)
                                                    }
                                                    P = P_tmp[order(R_tmp)]
                                                    R = R_tmp[order(R_tmp)]
                                                    }
390
                                                auc = auc(x = R, y = P, dens = 1)
Carine Rey's avatar
Carine Rey committed
391 392 393 394
                                                best_R_90_P = max(0,R[P>=0.9])
                                                P = c(P[1],P)  
                                                R = c(0,R)
                                                auc2 = auc(x = R, y = P, dens = 1)
395
                                                }
Carine Rey's avatar
Carine Rey committed
396 397 398
                                                
                                                
                                                df= data.frame(methode =x$methode[1] ,couple =x$couple[1] ,auc=auc,best_R_90_P=best_R_90_P,auc2=auc2)
399 400 401 402
                                                return(df)
                                                }))

df_auc = df_auc[order(df_auc$couple),]
Carine Rey's avatar
Carine Rey committed
403 404 405 406 407
condensed_meths = c("PCOC","Diffsel_mean","Identical_LG08","Mutinomial_1MinusLRT","Tdg09_1MinusLRT","Msd_1MinusP","Msd_0.05_1MinusP","Topological_LG08")
df_auc = subset(df_auc, methode %in% condensed_meths)
df_auc$auc_rank = unlist(tapply(df_auc$auc,df_auc$couple, function(x){rank(-x)}))
df_auc$auc2_rank = unlist(tapply(df_auc$auc2,df_auc$couple, function(x){rank(-x)}))
df_auc$best_R_90_P_rank = unlist(tapply(df_auc$best_R_90_P,df_auc$couple,function(x){rank(-x)}))
408 409
print(df_auc)

Carine Rey's avatar
Carine Rey committed
410
########################################################################
411 412


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

Carine Rey's avatar
Carine Rey committed
415
    nb_c = length(unique(df_out$couple))
Carine Rey's avatar
Carine Rey committed
416
    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
417
    colors2 = c("#984EA3","#4AA947","#377EB8","#E41A1C","#F5BE5B","#90EE90","#8B6914","#BFBFBF","#7F7F7F","#ADD8E6")
418

Carine Rey's avatar
Carine Rey committed
419 420 421 422
    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
423
    meths = meths[meths %in% df_out$methode]
Carine Rey's avatar
Carine Rey committed
424 425 426
    } else {
    meths = sort(unique(df_out$methode))
    }
427

Carine Rey's avatar
Carine Rey committed
428
    alpha = 0.7
Carine Rey's avatar
Carine Rey committed
429 430 431 432 433 434 435


    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)
436
    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
437 438 439 440 441 442 443

    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)
444
    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
445 446 447 448
    df_out_melt$couple =  factor(df_out_melt$couple, levels = couple_levels)

    print("plot recall_precision_per_meth")

449
    plot = ggplot(df_out, aes(x=sensitivity, y=precision98_02, col = couple))
Carine Rey's avatar
Carine Rey committed
450
    plot = plot + theme_bw()
451
    plot = plot + labs(x="Sensitivity (= Recall)", y="Precision (98/2)")
Carine Rey's avatar
Carine Rey committed
452 453 454 455
    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
456 457 458
    #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
459 460
    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
461
    plot = plot + ggtitle(paste0("Tree: ", opt$tree_prefix , "\n Date: " , date ))
Carine Rey's avatar
Carine Rey committed
462 463 464 465 466 467 468 469 470 471
    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
472
    print("plot recall_precision")
Carine Rey's avatar
Carine Rey committed
473

Carine Rey's avatar
Carine Rey committed
474
    if (length(meths) <= 9) {
475
        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
476
        plot = ggplot(subset(df_out, couple %in% couple_i), aes(x=sensitivity, y=precision98_02, col = methode))
477 478 479 480 481 482 483 484 485
        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
486
        plot = plot + facet_wrap(~ couple, ncol=4)
487 488 489 490 491 492 493 494 495
        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
496

497 498 499 500 501
        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
502

503 504 505 506 507 508 509


        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,
510 511
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
512 513
        
        plot_PR_c = function(couple_l) {
514 515 516 517 518 519
        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))
520 521 522 523 524 525 526 527 528 529 530
        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
531 532
        plot}
        
533
        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"))
534 535 536
        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"))
537 538 539 540 541 542 543 544 545

        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
546 547 548 549 550
        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,
551
                  labels = c("","Only Pressure Shift","Only Change in Selection Efficacy","Pressure Shift & Change in Selection Efficacy"),
552
                  align = "h",
Carine Rey's avatar
Carine Rey committed
553 554 555
                  
                  rel_heights = c( 0.3, 0.7,0.7,0.7),
                  hjust = 0, vjust = 0)
556 557 558 559 560



        save_plot(paste0(opt$out,suffix,".recall_precision_ok.pdf"),
                  plot_recall_precision_papier_ok,
561 562
                  ncol = 2.5,
                  nrow = 2.5,
563
                  base_aspect_ratio = 1,
564 565
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
566
    } else {
567
        plot_recall_precision_papier = NULL
Carine Rey's avatar
Carine Rey committed
568
    }
Carine Rey's avatar
Carine Rey committed
569

Carine Rey's avatar
Carine Rey committed
570
    print("plot FPR")
571
    #dcast
Carine Rey's avatar
Carine Rey committed
572 573
    x_labs = "Threshold"
    y_labs = "FP Rate (1 -Spe)"
Carine Rey's avatar
Carine Rey committed
574
    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
575 576 577
    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
578

Carine Rey's avatar
Carine Rey committed
579
    if (length(df_out_melt2$methode) > 0){
Carine Rey's avatar
Carine Rey committed
580
        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"))
581 582 583 584 585 586 587 588 589 590
        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))
Carine Rey's avatar
Carine Rey committed
591
        plot = plot + geom_step(size=1,direction="vh", alpha=alpha)
592 593 594 595 596 597 598 599 600 601 602 603 604 605
        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
606

Carine Rey's avatar
Carine Rey committed
607
        if (length(unique(df_out_melt2$couple)) == 4){
608 609 610
            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
611

612
            # complete values for not existing threshold between couple
Carine Rey's avatar
Carine Rey committed
613

614 615 616 617 618
            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
619

620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639
                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
640 641 642 643 644 645 646 647
            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")

648

Carine Rey's avatar
Carine Rey committed
649 650
            #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
651

Carine Rey's avatar
Carine Rey committed
652 653
            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
654

Carine Rey's avatar
Carine Rey committed
655 656
            #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"))
657
            print(head(df_out_melt2_dcast_melt))
Carine Rey's avatar
Carine Rey committed
658

Carine Rey's avatar
Carine Rey committed
659 660
            #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"))
661

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

664
            df_recall_sup09_per_meth3$variable = df_recall_sup09_per_meth3$couple
Carine Rey's avatar
Carine Rey committed
665
            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
666

667 668 669 670 671 672 673 674 675 676 677 678 679

            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
680
            plot = plot + scale_color_manual(values=colors[c(2,4,1,3)])
681 682 683 684 685 686 687
            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
688 689
    } else {
    plot_FPR = NULL
690
    plot_FPR2 = NULL
Carine Rey's avatar
Carine Rey committed
691
    }
Carine Rey's avatar
Carine Rey committed
692 693


Carine Rey's avatar
Carine Rey committed
694 695 696 697 698 699
    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
700
    plot = plot + theme(legend.position="top",
Carine Rey's avatar
Carine Rey committed
701 702
    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
703

Carine Rey's avatar
Carine Rey committed
704 705
    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
706
    #plot = plot + geom_point(size=0.5, alpha=alpha)
Carine Rey's avatar
Carine Rey committed
707
    plot = plot + geom_step(size=1, direction="vh",alpha=alpha)
Carine Rey's avatar
Carine Rey committed
708 709
    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
710
    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
711 712 713 714 715
    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
716
    save_plot(paste0(opt$out,suffix,".indicator_per_meth.pdf"),
Carine Rey's avatar
Carine Rey committed
717 718 719 720 721 722
    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
723

Carine Rey's avatar
Carine Rey committed
724
    plot_value_JSD = NULL
Carine Rey's avatar
Carine Rey committed
725
    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
726

Carine Rey's avatar
Carine Rey committed
727
        print("plot value/distance")
Carine Rey's avatar
Carine Rey committed
728 729


Carine Rey's avatar
Carine Rey committed
730

Carine Rey's avatar
Carine Rey committed
731
        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
732 733 734 735 736 737 738 739 740 741 742 743
        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
744

Carine Rey's avatar
Carine Rey committed
745 746 747 748 749

        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
750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802
        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
803

Carine Rey's avatar
Carine Rey committed
804 805 806 807 808 809 810 811

        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
812

Carine Rey's avatar
Carine Rey committed
813

Carine Rey's avatar
Carine Rey committed
814
    }
Carine Rey's avatar
Carine Rey committed
815 816 817 818 819

    #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