# 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