calc_t_per_meth.R 44.4 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)}))
Carine Rey's avatar
Carine Rey committed
408 409 410 411 412 413 414 415 416 417 418 419 420 421

df_auc$string_auc = paste0(df_auc$auc_rank, " (", round(df_auc$auc,3),")")

df_auc$string_auc2 = paste0(df_auc$auc_rank, " (", round(df_auc$auc2,3),")")
df_auc$string_best_R_90_P = paste0(df_auc$best_R_90_P_rank, " (", round(df_auc$best_R_90_P,3),")")

df_only_auc = dcast(df_auc, methode  ~ couple, value.var="string_auc" )
colnames(df_only_auc)[1] = "methode\auc"
df_only_auc2 = dcast(df_auc, methode  ~ couple, value.var="string_auc2" )
colnames(df_only_auc2)[1] = "methode\auc2"
df_only_string_best_R_90_P = dcast(df_auc, methode  ~ couple, value.var="string_best_R_90_P" )
colnames(df_only_string_best_R_90_P)[1] = "methode\best_R_90_P"

print(df_only_auc)
422

Carine Rey's avatar
Carine Rey committed
423
########################################################################
424 425


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

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

Carine Rey's avatar
Carine Rey committed
432 433 434 435
    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
436
    meths = meths[meths %in% df_out$methode]
Carine Rey's avatar
Carine Rey committed
437 438 439
    } else {
    meths = sort(unique(df_out$methode))
    }
440

Carine Rey's avatar
Carine Rey committed
441
    alpha = 0.7
Carine Rey's avatar
Carine Rey committed
442 443 444 445 446 447 448


    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)
449
    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
450 451 452 453 454 455 456

    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)
457
    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
458 459 460 461
    df_out_melt$couple =  factor(df_out_melt$couple, levels = couple_levels)

    print("plot recall_precision_per_meth")

462
    plot = ggplot(df_out, aes(x=sensitivity, y=precision98_02, col = couple))
Carine Rey's avatar
Carine Rey committed
463
    plot = plot + theme_bw()
464
    plot = plot + labs(x="Sensitivity (= Recall)", y="Precision (98/2)")
Carine Rey's avatar
Carine Rey committed
465 466 467 468
    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
469 470 471
    #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
472 473
    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
474
    plot = plot + ggtitle(paste0("Tree: ", opt$tree_prefix , "\n Date: " , date ))
Carine Rey's avatar
Carine Rey committed
475 476 477 478 479 480 481 482 483 484
    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
485
    print("plot recall_precision")
Carine Rey's avatar
Carine Rey committed
486

Carine Rey's avatar
Carine Rey committed
487
    if (length(meths) <= 9 & all(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") %in% df_out$couple)) {
488
        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
489
        plot = ggplot(subset(df_out, couple %in% couple_i), aes(x=sensitivity, y=precision98_02, col = methode))
490 491 492 493 494 495 496 497 498
        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
499
        plot = plot + facet_wrap(~ couple, ncol=4)
500 501 502 503 504 505 506 507 508
        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
509

510 511 512 513 514
        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
515

516 517 518 519 520 521 522


        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,
523 524
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
525 526
        
        plot_PR_c = function(couple_l) {
527 528 529 530 531 532
        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))
533 534 535 536 537 538 539 540 541 542 543
        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
544 545
        plot}
        
546
        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"))
547 548 549
        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"))
550 551 552 553 554 555 556 557 558

        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
559 560 561 562 563
        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,
564
                  labels = c("","Only Pressure Shift","Only Change in Selection Efficacy","Pressure Shift & Change in Selection Efficacy"),
565
                  align = "h",
Carine Rey's avatar
Carine Rey committed
566 567 568
                  
                  rel_heights = c( 0.3, 0.7,0.7,0.7),
                  hjust = 0, vjust = 0)
569 570 571 572 573



        save_plot(paste0(opt$out,suffix,".recall_precision_ok.pdf"),
                  plot_recall_precision_papier_ok,
574 575
                  ncol = 2.5,
                  nrow = 2.5,
576
                  base_aspect_ratio = 1,
577 578
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
579
    } else {
580
        plot_recall_precision_papier = NULL
Carine Rey's avatar
Carine Rey committed
581 582 583 584 585 586 587 588
        plot_recall_precision_papier_ok = NULL
        save_plot(paste0(opt$out,suffix,".recall_precision_ok.pdf"),
                  plot_recall_precision_papier_ok,
                  ncol = 2.5,
                  nrow = 2.5,
                  base_aspect_ratio = 1,
                  limitsize = FALSE
                  )
Carine Rey's avatar
Carine Rey committed
589
    }
Carine Rey's avatar
Carine Rey committed
590

Carine Rey's avatar
Carine Rey committed
591
    print("plot FPR")
592
    #dcast
Carine Rey's avatar
Carine Rey committed
593 594
    x_labs = "Threshold"
    y_labs = "FP Rate (1 -Spe)"
Carine Rey's avatar
Carine Rey committed
595
    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
596 597 598
    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
599

Carine Rey's avatar
Carine Rey committed
600
    if (length(df_out_melt2$methode) > 0){
Carine Rey's avatar
Carine Rey committed
601
        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"))
602 603 604 605 606 607 608 609 610 611
        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
612
        plot = plot + geom_step(size=1,direction="vh", alpha=alpha)
613 614 615 616 617 618 619 620 621 622 623 624 625 626
        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
627

Carine Rey's avatar
Carine Rey committed
628
        if (length(unique(df_out_melt2$couple)) == 4){
629 630 631
            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
632

633
            # complete values for not existing threshold between couple
Carine Rey's avatar
Carine Rey committed
634

635 636 637 638 639
            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
640

641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
                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
661 662 663 664 665 666 667 668
            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")

669

Carine Rey's avatar
Carine Rey committed
670 671
            #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
672

Carine Rey's avatar
Carine Rey committed
673 674
            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
675

Carine Rey's avatar
Carine Rey committed
676 677
            #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"))
678
            print(head(df_out_melt2_dcast_melt))
Carine Rey's avatar
Carine Rey committed
679

Carine Rey's avatar
Carine Rey committed
680 681
            #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"))
682

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

685
            df_recall_sup09_per_meth3$variable = df_recall_sup09_per_meth3$couple
Carine Rey's avatar
Carine Rey committed
686
            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
687

688 689 690 691 692 693 694 695 696 697 698 699 700

            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
701
            plot = plot + scale_color_manual(values=colors[c(2,4,1,3)])
702 703 704 705 706 707 708
            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
709 710
    } else {
    plot_FPR = NULL
711
    plot_FPR2 = NULL
Carine Rey's avatar
Carine Rey committed
712
    }
Carine Rey's avatar
Carine Rey committed
713 714


Carine Rey's avatar
Carine Rey committed
715 716 717 718 719 720
    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
721
    plot = plot + theme(legend.position="top",
Carine Rey's avatar
Carine Rey committed
722 723
    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
724

Carine Rey's avatar
Carine Rey committed
725 726
    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
727
    #plot = plot + geom_point(size=0.5, alpha=alpha)
Carine Rey's avatar
Carine Rey committed
728
    plot = plot + geom_step(size=1, direction="vh",alpha=alpha)
Carine Rey's avatar
Carine Rey committed
729 730
    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
731
    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
732 733 734 735 736
    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
737
    save_plot(paste0(opt$out,suffix,".indicator_per_meth.pdf"),
Carine Rey's avatar
Carine Rey committed
738 739 740 741 742 743
    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
744

Carine Rey's avatar
Carine Rey committed
745
    plot_value_JSD = NULL
Carine Rey's avatar
Carine Rey committed
746
    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
747

Carine Rey's avatar
Carine Rey committed
748
        print("plot value/distance")
Carine Rey's avatar
Carine Rey committed
749 750


Carine Rey's avatar
Carine Rey committed
751

Carine Rey's avatar
Carine Rey committed
752
        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
753 754 755 756 757 758 759 760 761 762 763 764
        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
765

Carine Rey's avatar
Carine Rey committed
766 767 768 769 770

        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
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 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823
        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