--- a
+++ b/6-Figure scripts/Fig S6.R
@@ -0,0 +1,113 @@
+library(data.table)
+library(dplyr)
+library(reshape2)
+library(ggplot2)
+library(readxl)
+
+
+list.files()
+
+excel_sheets("Fig S6 Source Data.xlsx")
+
+# Fig S6 ------------------------------------------
+# Fig S6. LOSO density plot #######
+rm(list = ls())
+dat.hist <- read_excel("Fig S6 Source Data.xlsx", sheet = "LOSO")
+colnames(dat.hist) <- c("Figs_map", "Figs_color","species","X1")
+dat.hist <- dat.hist %>% mutate(X2=X1, Y1=-0.002, Y2=-0.009)
+
+taxonomy_df <- read_excel("Fig S6 Source Data.xlsx", sheet = "taxonomy")
+
+genus <- sapply(strsplit(dat.hist$species," ", fixed = T),"[[", 1)
+genus[!genus %in% taxonomy_df$Genus]
+
+phylum <- sapply(genus, 
+                 function(x){
+                   if(x %in% taxonomy_df$Genus) {
+                     taxonomy_df$Phylum[which(taxonomy_df$Genus == x)[1]]
+                   }else "Unclassified"
+                 })
+phylum_color_df <- cbind.data.frame(phylum=c("Bacteroidetes","Actinobacteria","TM7","Proteobacteria","Firmicutes","other"),
+                                    colors=c("#EF5656","#47B3DA","#9A8FC3","#F7A415","#2BB065","#BABABA"),
+                                    stringsAsFactors=F)
+dat.hist$phylum <- phylum
+dat.hist$phylum_other <- sapply(dat.hist$phylum, function(x)if(x %in% phylum_color_df$phylum) x else "other")
+
+Fig.map_ymax_df <- cbind.data.frame(Fig.map = unique(dat.hist$Figs_map),
+                                    # ymax=c(0.04,0.16,0.03,0.015,0.035,0.038),
+                                    ymax = c(0.06,0.19, 0.06, 0.05),
+                                    stringsAsFactors=F)
+
+
+for(Fig in unique(dat.hist$Figs_map)){
+  
+  #Fig=unique(dat.hist$Figs_map)[1]
+  
+  ymax = Fig.map_ymax_df$ymax[which(Fig.map_ymax_df$Fig.map == Fig)]
+  
+  dat_sub <-  dat.hist %>%  filter(Figs_map %in% Fig)
+  
+  phylum_ABC <- unique(dat_sub$phylum_other)[order(unique(dat_sub$phylum_other))]
+  phylum_rank <- c(phylum_ABC[phylum_ABC != "other"],"other")
+  
+  dat_sub$phylum_other <- factor(dat_sub$phylum_other, 
+                                 levels = phylum_rank)
+  Colors <- sapply(phylum_rank, function(x) phylum_color_df$colors[which(phylum_color_df$phylum ==x)])
+  
+  histP <- ggplot(dat_sub) +
+    geom_density(aes(x=X1,y=(..count..)/sum(..count..)),fill = "gray") + ylab("Density") + 
+    xlab("") +
+    theme_bw()+theme(panel.grid = element_blank()) +
+    ylim(-0.01,ymax) +
+    geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2, color=phylum_other )) +
+    scale_color_manual(values = Colors) + 
+    theme(legend.position = "none") + ggtitle(Fig) 
+  
+  
+  histP
+  
+  assign(paste("histP_", Fig, sep = ""), histP, envir = .GlobalEnv)
+}
+
+histP_MetaG_M00044
+histP_MetaG_P00250
+histP_MetaG_P00380
+histP_MetaG_P00564
+
+
+library(ggpubr)
+FigS6.LOSO <- ggarrange(histP_MetaG_P00380, histP_MetaG_M00044, histP_MetaG_P00250, histP_MetaG_P00564, ncol = 1)
+
+
+# Fig S6. KO contribution #######
+
+zscore_df <- read_excel("Fig S6 Source Data.xlsx", sheet = "KO contrib" )
+colnames(zscore_df) <- c("K","sp","score")
+
+tmp <- zscore.top3_df <- zscore_df %>%
+  group_by(K) %>%
+  top_n(n = 3, wt = score) %>% as.data.frame()
+tmp <- tmp %>% arrange(desc(score)) %>% arrange(K)
+
+zscore.top3_df <- tmp %>% #keep the top 3 with highest score
+  mutate(X = rep(c("X1","X2","X3"), length(unique(zscore_df$K)))) %>%
+  mutate(genus = sapply(strsplit(sp," ", fixed = T),"[[", 1)) %>%
+  mutate(phylum = sapply(genus,
+                         function(x){
+                           if(x %in% taxonomy_df$Genus){
+                             taxonomy_df$Phylum[which(taxonomy_df$Genus == x)[1]]
+                           }else "Unclassified"
+                         } ) ) %>%
+  mutate(phylum_other = sapply(phylum, function(x) if(x %in% phylum_color_df$phylum) x else "other"))
+
+zscore.top3_df$phylum_other <- factor(zscore.top3_df$phylum_other, levels = phylum_rank)
+
+
+FigS6.KOcontri <- ggplot(zscore.top3_df) +
+  geom_point(aes(x=X, y=K, size=score, fill=phylum_other), shape=21) +
+  scale_fill_manual(values = Colors)  + 
+  scale_size (range = c (3, 6)) +
+  theme(panel.grid = element_blank(), axis.title = element_blank(), panel.background = element_blank(),
+        axis.text.x = element_blank(), legend.position = "none", axis.ticks = element_blank())
+
+FigS6.KOcontri