--- a
+++ b/6-Figure scripts/Fig E2.r
@@ -0,0 +1,67 @@
+# Figure E2a. Total bacterial load ####
+dat <- read_excel("Fig E2 Source Data.xlsx", sheet = "Bacterial load")
+
+library(ggpubr)
+FigE2a <- ggplot(dat) +
+  geom_boxplot(aes(x=Group_Site, y=log10cpn, fill=Group),outlier.shape = NA, alpha=0.5) +
+  geom_jitter(aes(x=Group_Site, y=log10cpn, fill=Group),shape=21,size=2, width = 0.2) +
+  theme_bw() + theme(panel.grid = element_blank()) +
+  ylab("log10 (16S copy number)")
+
+FigE2a
+
+wilcox.test(log10cpn~Group, data = dat %>% as.data.frame() %>% dplyr::filter(Site=="Guangzhou")) 
+wilcox.test(log10cpn~Group, data = dat %>% as.data.frame() %>% dplyr::filter(Site=="Shenzhen")) 
+
+
+# Figure E2.b --------------------
+# Figure E2b. taxonomy #####
+dat <- read_excel("Fig E2 Source Data.xlsx", sheet = "Taxonomy")
+
+FigE2b.taxonomy <- ggplot(dat) +
+  geom_point(aes(x=NUE,y=EOS), size=2)+
+  geom_text_repel(aes(x=NUE,y=EOS, label = Label),size=3) +
+  theme_bw()+ theme(panel.grid = element_blank()) +
+  geom_hline(yintercept = 0, linetype="twodash") +
+  geom_hline(yintercept = 1, linetype="dashed") +
+  geom_vline(xintercept = -1, linetype="dashed") +
+  geom_vline(xintercept = 0, linetype="twodash") +
+  geom_vline(xintercept = 1, linetype="dashed") +
+  xlab("Correlation NEU: Directionality x - log(P)")+
+  ylab("Correlation EOS: Directionality x - log(P)")+
+  ggtitle("Differential species-level taxa")
+FigE2b.taxonomy
+
+# Figure E2b. KEGG modules #####
+dat <- read_excel("Fig E2 Source Data.xlsx", sheet = "Metagenome")
+
+FigE2b.KEGG <- ggplot(dat) +
+  geom_point(aes(x=NEU,y=EOS), size=2) +
+  theme_bw()+ theme(panel.grid = element_blank()) +
+  geom_hline(yintercept = 0, linetype="twodash") +
+  geom_hline(yintercept = 1, linetype="dashed") +
+  geom_vline(xintercept = -1, linetype="dashed") +
+  geom_vline(xintercept = 0, linetype="twodash") +
+  geom_vline(xintercept = 1, linetype="dashed") +
+  xlab("Correlation NEU: Directionality x - log(P)")+
+  ylab("Correlation EOS: Directionality x - log(P)")+
+  ggtitle("Differential KEGG modules")
+FigE2b.KEGG
+
+
+
+# Figure E2c. correlation #####
+dat <- read_excel("Fig E2 Source Data.xlsx", sheet = "Correlation")
+
+FigE2c <- ggplot(dat,aes(x=bin_based_FC, y=reads_based_FC)) +
+  geom_point( size=2 )+
+  theme_classic() + theme(panel.grid = element_blank()) +
+  geom_hline(yintercept = 0, linetype="dashed") +
+  geom_vline(xintercept = 0, linetype="dashed") +
+  geom_smooth(method = "lm", se = F, linetype="dashed", color="black", size=0.5) +
+  xlim(c(-20, 20)) + ylim(c(-3, 10)) +
+  scale_y_continuous(breaks = seq(-2, 10, 2))+
+  xlab("FC (Reads-based)")+
+  ylab("FC (Bins-based)")
+FigE2c
+