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

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

c_detail <- c_detail[,17:40]
s_detail <- s_detail[,17:40]

reaction <- as.data.frame(c(rep('Coral-associated bacteria',6),rep('Seawater bacteria',6)))
colnames(reaction) <- 'Models'
reaction[,2] <- as.data.frame(rep(c('CarveMe/Consensus','KBase/Consensus','gapseq/Consensus',
                                    'CarveMe/gapseq','CarveMe/KBase','gapseq/KBase'),2))
colnames(reaction)[2] <- 'Model'

for(i in 1:6){
  reaction[i,3] <- mean(c_detail[,i])
  reaction[i+6,3] <- mean(s_detail[,i])
}
colnames(reaction)[3] <- 'Mean'

for(i in 1:6){
  reaction[i,4] <- sd(c_detail[,i])
  reaction[i+6,4] <- sd(s_detail[,i])
}
colnames(reaction)[4] <- 'SD'

reaction$Model <- factor(reaction$Model, levels=c('gapseq/Consensus', 'CarveMe/Consensus', 
                                                  'KBase/Consensus', 'gapseq/KBase',
                                                  'CarveMe/KBase','CarveMe/gapseq'))

library(ggplot2)
rxns <- ggplot(reaction, aes(x = Model, y = Mean))+
  ylim(0,0.8)+
  theme_bw()+
  ggtitle('a. Reactions')+
  labs(y = 'Jaccard index')+
  theme(text = element_text(size = 10),
        legend.title = element_blank(),
        legend.text = element_text(size=13),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold"),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
        panel.grid = element_blank())+
  geom_point(aes(color = Models), position = position_dodge(0.3), size = 3)+
  geom_errorbar(aes(ymin = Mean-SD, ymax=Mean+SD, color = Models), position = position_dodge(.3),width = 0.2)


######### metabolites #################
metabolite <- reaction

for(i in 1:6){
  metabolite[i,3] <- mean(c_detail[,i+6])
  metabolite[i+6,3] <- mean(s_detail[,i+6])
}

for(i in 1:6){
  metabolite[i,4] <- sd(c_detail[,i+6])
  metabolite[i+6,4] <- sd(s_detail[,i+6])
}

mets <- ggplot(metabolite, aes(x = Model, y = Mean))+
  ylim(0,0.8)+
  theme_bw()+
  ggtitle('b. Metabolites')+
  labs(y = 'Jaccard index')+
  theme(text = element_text(size = 10),
        legend.title = element_blank(),
        legend.text = element_text(size=13),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold"),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
        panel.grid = element_blank())+
  geom_point(aes(color = Models), position = position_dodge(0.3), size = 3)+
  geom_errorbar(aes(ymin = Mean-SD, ymax=Mean+SD, color = Models), position = position_dodge(.3),width = 0.2)

############ dead-end metabolite ##############
dead <- reaction

for(i in 1:6){
  dead[i,3] <- mean(c_detail[,i+12])
  dead[i+6,3] <- mean(s_detail[,i+12])
}

for(i in 1:6){
  dead[i,4] <- sd(c_detail[,i+12])
  dead[i+6,4] <- sd(s_detail[,i+12])
}

dead_end <- ggplot(dead, aes(x = Model, y = Mean))+
  ylim(0,0.8)+
  theme_bw()+
  ggtitle('c. Dead_end_mets')+
  labs(y = 'Jaccard index')+
  theme(text = element_text(size = 10),
        legend.title = element_blank(),
        legend.text = element_text(size=13),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold"),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
        panel.grid = element_blank())+
  geom_point(aes(color = Models), position = position_dodge(0.3), size = 3)+
  geom_errorbar(aes(ymin = Mean-SD, ymax=Mean+SD, color = Models), position = position_dodge(.3),width = 0.2)

############ genes ##############
genes <- reaction

for(i in 1:6){
  genes[i,3] <- mean(c_detail[,i+18])
  genes[i+6,3] <- mean(s_detail[,i+18])
}

for(i in 1:6){
  genes[i,4] <- sd(c_detail[,i+18])
  genes[i+6,4] <- sd(s_detail[,i+18])
}

gen <- ggplot(genes, aes(x = Model, y = Mean))+
  ylim(0,0.8)+
  theme_bw()+
  ggtitle('d. Genes')+
  labs(y = 'Jaccard index')+
  theme(text = element_text(size = 10),
        legend.title = element_blank(),
        legend.text = element_text(size=13),
        plot.title = element_text(face = "bold"),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold"),
        axis.text.y = element_text(size = 13),
        axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
        panel.grid = element_blank())+
  geom_point(aes(color = Models), position = position_dodge(0.3), size = 3)+
  geom_errorbar(aes(ymin = Mean-SD, ymax=Mean+SD, color = Models), position = position_dodge(.3),width = 0.2)

library(ggpubr)
ggarrange(rxns, mets, dead_end, gen, 
          common.legend = TRUE,
          legend = 'top',
          ncol = 2,nrow = 2)


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