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

library(ggplot2)
library(ggpubr)

detail <- read.table('figures/details_coral.csv',header = T,sep = ',')
detail_SW <- read.table('figures/details_seawter.csv',header = T,sep = ',')

# extract the number of reactions
# coral models
data <- as.data.frame(c(detail$c_reactions,
                     detail$c_metabolites,
                     detail$c_dead_end_mets,
                     detail$c_genes,
                     detail$g_reactions,
                     detail$g_metabolites,
                     detail$g_dead_end_mets,
                     detail$g_genes,
                     detail$k_reactions,
                     detail$k_metabolites,
                     detail$k_dead_end_mets,
                     detail$k_genes,
                     detail$m_reactions,
                     detail$m_metabolites,
                     detail$m_dead_end_mets,
                     detail$m_genes))


vari = c(rep('Reactions',50),rep('Metabolites',50),rep('Dead_end_mets',50),rep('Genes',50))
data[2] <- as.data.frame(rep(vari,4))

data[3] <- as.data.frame(c(rep('CarveMe',200),rep('gapseq',200),rep('KBase',200),rep('Consensus',200)))

colnames(data) <- c('Number','Details','Methods')
data$Details <- factor(data$Details , levels=c('Reactions', 'Metabolites', 'Dead_end_mets', 'Genes'))
data$Methods <- factor(data$Methods , levels=c('CarveMe', 'gapseq', 'KBase', 'Consensus'))

# seawater models
data_SW <- as.data.frame(c(detail_SW$c_reactions,
                           detail_SW$c_metabolites,
                           detail_SW$c_dead_end_mets,
                           detail_SW$c_genes,
                           detail_SW$g_reactions,
                           detail_SW$g_metabolites,
                           detail_SW$g_dead_end_mets,
                           detail_SW$g_genes,
                           detail_SW$k_reactions,
                           detail_SW$k_metabolites,
                           detail_SW$k_dead_end_mets,
                           detail_SW$k_genes,
                           detail_SW$m_reactions,
                           detail_SW$m_metabolites,
                           detail_SW$m_dead_end_mets,
                           detail_SW$m_genes))


vari_SW = c(rep('Reactions',55),rep('Metabolites',55),rep('Dead_end_mets',55),rep('Genes',55))
data_SW[2] <- as.data.frame(rep(vari_SW,4))

data_SW[3] <- as.data.frame(c(rep('CarveMe',220),rep('gapseq',220),rep('KBase',220),rep('Consensus',220)))

colnames(data_SW) <- c('Number','Details','Methods')
data_SW$Details <- factor(data_SW$Details , levels=c('Reactions', 'Metabolites', 'Dead_end_mets', 'Genes'))
data_SW$Methods <- factor(data_SW$Methods , levels=c('CarveMe', 'gapseq', 'KBase', 'Consensus'))


coral <- ggplot(data, aes(Details, Number)) +
  geom_boxplot(aes(color = Methods))+
  ylim(0,5000)+
  theme_bw()+
  ggtitle('a. Coral-associated bacteria')+
  stat_compare_means(aes(group = Methods), label = "p.signif", label.y = 4500)+
  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),
        panel.grid = element_blank())

seawater <- ggplot(data_SW, aes(Details, Number)) +
  geom_boxplot(aes(color = Methods))+
  ylim(0,5000)+
  theme_bw()+
  ggtitle('b. Seawater bacteria')+
  stat_compare_means(aes(group = Methods), label = "p.signif", label.y = 4500)+
  theme(text = element_text(size = 10),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
        panel.grid = element_blank())

ggarrange(coral, seawater, 
          common.legend = TRUE,
          legend = 'right')

ggsave(filename = "~/comparison of tools/Figure_1.svg",
       width = 20,height = 10,units = "cm")
