-
Notifications
You must be signed in to change notification settings - Fork 1
/
Balanced_optimal_evalue_MCC.R
111 lines (87 loc) · 3.09 KB
/
Balanced_optimal_evalue_MCC.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
library(dplyr)
library(ggplot2)
library(mltools)
library(gridExtra)
library(reshape2)
setwd("~/Documents/research/Vakirlis_Carvunis_McLysaght_2019/")
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
#
# load data
#
df_fn <- read.csv("Figure_3-source_data_1.csv")
df_fn$tag <- factor(df_fn$tag, levels=c("yeast", "fruitfly", "human"), ordered = TRUE)
#
df_fn <- filter(df_fn, species!="Caenorhabditis_elegans")
#
cer_df_tbl <- filter(df_fn, tag=="yeast")
dros_df_tbl <- filter(df_fn, tag=="fruitfly")
ver_df_tbl <- filter(df_fn, tag=="human")
df_fp <- read.csv("Figure_3-source_data_2.csv")
df_fp$tag <- factor(df_fp$tag, levels=c("yeast", "fruitfly", "human"), ordered = TRUE)
#
df_fp <- filter(df_fp, species!="Caenorhabditis_elegans")
#
cer_fp_df_tbl <- filter(df_fp, tag=="yeast")
dros_fp_df_tbl <- filter(df_fp, tag=="fruitfly")
ver_fp_df_tbl <- filter(df_fp, tag=="human")
#
# calculate FP, FN, TP, TN in each dataset
#
cer_f1 <- mutate(cer_fp_df_tbl,
FPR = found/total,
TPR = 1-(cer_df_tbl$not_found/cer_df_tbl$total),
TP = cer_df_tbl$total - cer_df_tbl$not_found,
FP = found,
TN = total-found,
FN = cer_df_tbl$not_found)
dros_df_tbl <- filter(dros_df_tbl, evalue %in% dros_fp_df_tbl$evalue)
dros_f1 <- mutate(dros_fp_df_tbl,
FPR = found/total,
TPR = 1-(dros_df_tbl$not_found/dros_df_tbl$total),
TP = dros_df_tbl$total - dros_df_tbl$not_found,
FP = found,
TN = total-found,
FN = dros_df_tbl$not_found)
ver_f1 <- mutate(ver_fp_df_tbl,
FPR = found/total,
TPR = 1-(ver_df_tbl$not_found/ver_df_tbl$total),
TP = ver_df_tbl$total - ver_df_tbl$not_found,
FP = found,
TN = total-found,
FN = ver_df_tbl$not_found)
#
# calculate F1 scores and MCC for each dataset
#
f1_score <- function(tp, fp, fn, tn)
{
precision = tp/(tp+fp)
recall = tp/(tp+fn)
return((2*(precision*recall))/(precision+recall))
}
cer_f1_final = cer_f1 %>%
group_by(evalue, species) %>%
dplyr::summarize(f1score = f1_score(TP, FP, FN, TN),
mcc = round(mcc(TP=TP, FP=FP, FN=FN, TN=TN),2),
FPR=FPR) %>%
group_by(species) %>%
slice(which(mcc == max(mcc))) %>%
group_by(species) %>%
slice(which.max(evalue))
dros_f1_final = dros_f1 %>%
group_by(evalue, species) %>%
dplyr::summarize(f1score = f1_score(TP, FP, FN, TN),
mcc = round(mcc(TP=TP, FP=FP, FN=FN, TN=TN),2),
FPR=FPR) %>%
group_by(species) %>%
slice(which(mcc == max(mcc))) %>%
group_by(species) %>%
slice(which.max(evalue))
ver_f1_final = ver_f1 %>%
group_by(evalue, species) %>%
dplyr::summarize(f1score = f1_score(TP, FP, FN, TN),
mcc = round(mcc(TP=TP, FP=FP, FN=FN, TN=TN),2),
FPR=FPR) %>%
group_by(species) %>%
slice(which(mcc == max(mcc))) %>%
group_by(species) %>%
slice(which.max(evalue))