rm(list=ls())
setwd('~/comparison of tools/data')

###### Number of exchange reactions ######
rxns <- read.table('figures/exchange_rxns.csv',header = T,sep = ',')

t_rxns <- as.data.frame(c(rxns$c_carve,rxns$c_gapseq,rxns$c_kbase,rxns$c_merge,rxns$s_carve,rxns$s_gapseq,rxns$s_kbase,rxns$s_merge))
t_rxns[2] <- as.data.frame(c(rep('CarveMe',55),
                             rep('gapseq',55),
                             rep('KBase',55),
                             rep('Consensus',55),
                             rep('CarveMe',55),
                             rep('gapseq',55),
                             rep('KBase',55),
                             rep('Consensus',55)))
t_rxns[3] <- as.data.frame(c(rep('Coral',220),
                             rep('Seawater',220)))
colnames(t_rxns) <- c('Number', 'Methods','Type') 

t_rxns$Methods <- factor(t_rxns$Methods, levels=c('CarveMe', 'gapseq', 'KBase', 'Consensus'))

coral <- t_rxns[which(t_rxns$Type == 'Coral'),]
sea <- t_rxns[which(t_rxns$Type == 'Seawater'),]

library(ggplot2)
library(ggpubr)
coral_ex <- ggplot(coral, aes(Methods, Number)) +
  geom_boxplot(aes(color = Methods))+
  ylim(50,250)+
  theme_bw()+
  ggtitle('b.')+
  labs(y = 'Number of exchange reactions')+
  stat_compare_means(comparisons = list(c('CarveMe','gapseq'),c('CarveMe','KBase'),c('CarveMe','Consensus'),
                                        c('gapseq','KBase'),c('gapseq','Consensus'),c('KBase','Consensus')),
                     label.y = c(170,185,200,215,225,235))+
  theme(text = element_text(size = 10),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold", size = 13),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1,face = "bold"),
        panel.grid = element_blank())

sea_ex <- ggplot(sea, aes(Methods, Number)) +
  geom_boxplot(aes(color = Methods))+
  ylim(50,250)+
  theme_bw()+
  ggtitle('b.')+
  labs(y = 'Number of exchange reactions')+
  stat_compare_means(comparisons = list(c('CarveMe','gapseq'),c('CarveMe','KBase'),c('CarveMe','Consensus'),
                                        c('gapseq','KBase'),c('gapseq','Consensus'),c('KBase','Consensus')),
                     label.y = c(170,185,200,215,225,235))+
  theme(text = element_text(size = 10),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold", size = 13),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1,face = "bold"),
        panel.grid = element_blank())

###### Number of added reactions ######
habitat <- c('coral', 'seawater')
methods <- c('carveme','gapseq', 'kbase','consensus')

coral_add <- data_frame()
sea_add <- data_frame()

for (i in 1:length(habitat)){
  for (j in 1:length(methods)){
    fil_dir <- paste0(habitat[i], '_', methods[j], '_gf_vs_abundance.txt')
    data <- read.table(file = fil_dir, header = T, sep = "\t")
    if (i == 1){
      coral_add <- as.data.frame(c(coral_add,data$added))
    }
    else if (i == 2){
      sea_add <- as.data.frame(c(sea_add,data$added))
    }
  }
}

coral_add = as.data.frame(t(coral_add))
coral_add[2] <- as.data.frame(c(rep('CarveMe',50),
                             rep('gapseq',50),
                             rep('KBase',50),
                             rep('Consensus',50)))
colnames(coral_add) <- c('Number', 'Methods')
coral_add$Methods <- factor(coral_add$Methods, levels=c('CarveMe', 'gapseq', 'KBase', 'Consensus'))

sea_add = as.data.frame(t(sea_add))
sea_add[2] <- as.data.frame(c(rep('CarveMe',55),
                                rep('gapseq',55),
                                rep('KBase',55),
                                rep('Consensus',55)))
colnames(sea_add) <- c('Number', 'Methods') 
sea_add$Methods <- factor(sea_add$Methods, levels=c('CarveMe', 'gapseq', 'KBase', 'Consensus'))

coral_add_p <- ggplot(coral_add, aes(Methods, Number)) +
  geom_boxplot(aes(color = Methods))+
  ylim(0,250)+
  theme_bw()+
  ggtitle('a.')+
  labs(y = 'Number of added reactions')+
  stat_compare_means(comparisons = list(c('CarveMe','gapseq'),c('CarveMe','KBase'),c('CarveMe','Consensus'),
                                        c('gapseq','KBase'),c('gapseq','Consensus'),c('KBase','Consensus')),
                     label.y = c(170,185,200,215,225,235))+
  theme(text = element_text(size = 10),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold", size = 13),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1,face = "bold"),
        panel.grid = element_blank())

sea_add_p <- ggplot(sea_add, aes(Methods, Number)) +
  geom_boxplot(aes(color = Methods))+
  ylim(0,250)+
  theme_bw()+
  ggtitle('a.')+
  labs(y = 'Number of added reactions')+
  stat_compare_means(comparisons = list(c('CarveMe','gapseq'),c('CarveMe','KBase'),c('CarveMe','Consensus'),
                                        c('gapseq','KBase'),c('gapseq','Consensus'),c('KBase','Consensus')),
                     label.y = c(170,185,200,215,225,235))+
  theme(text = element_text(size = 10),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold", size = 13),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1,face = "bold"),
        panel.grid = element_blank())

ggarrange(coral_add_p, coral_ex, 
          common.legend = TRUE,
          legend = 'right')

ggsave(filename = "~/comparison of tools/Figure_4.svg",
       width = 25,height = 15,units = "cm")

ggarrange(sea_add_p, sea_ex, 
          common.legend = TRUE,
          legend = 'right')

ggsave(filename = "~/comparison of tools/Figure_S4.svg",
       width = 25,height = 15,units = "cm")
