library(readxl)
library(dplyr)



#Figure S1a --------------------------------------
# Figure S1a. genes ####
dat.genes <- read_excel("Fig S1 Source Data.xlsx", sheet = "Fig S1a rarefaction genes")
dat.KO <- read_excel("Fig S1 Source Data.xlsx", sheet = "Fig S1a rarefaction KOs")


data1<-dat.genes[(dat.genes$Group=="COPD"),]
data2<-dat.genes[(dat.genes$Group=="Health"),]
dat_combined <- merge(data1,data2, by = "Samples",all = T)


library(ggplot2)
p1<-ggplot(dat_combined,aes(x=factor(Samples),y=Genes.x/1000000))+
  geom_boxplot(colour="#f0999f",outlier.colour=NULL,outlier.shape=19,outlier.size=1)+
  #scale_y_continuous(limits = c(0,2))  +
  theme_bw() + theme(panel.grid.minor =element_blank(), panel.grid.major = element_blank()) 
p1

library(ggpubr)
library(cowplot)
p2<-ggplot(dat_combined,aes(x=factor(Samples),y=Genes.y/1000000)) + 
  geom_boxplot(colour="#46bbc0",outlier.colour=NULL,outlier.shape=19,outlier.size=1)+
  scale_y_continuous(limits = c(0,1.25), position = "right",
                     breaks = seq(0,1.25, by=0.25))  +
  theme_half_open(11, rel_small = 1) + 
  rremove("x.axis")+
  rremove("xlab") +
  rremove("x.text") +
  rremove("x.ticks") +
  rremove("legend") 

p2

aligned_plots <- align_plots(p1, p2, align="hv", axis="tblr")
FigS1a.genes <- ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])
FigS1a.genes


# Figure S1a. KO ####

data1<-dat.KO[(dat.KO$Group=="COPD"),]
data2<-dat.KO[(dat.KO$Group=="Health"),]
dat_combined <- merge(data1,data2, by = "Samples",all = T)

p1<-ggplot(dat_combined,aes(x=factor(Samples),y=Genes.x/1000))+
  geom_boxplot(colour="#f0999f",outlier.colour=NULL,outlier.shape=19,outlier.size=1)+
  scale_y_continuous(limits = c(2.5,7))  +
  theme_bw() + theme(panel.grid.minor =element_blank(), panel.grid.major = element_blank()) 
p1


p2<-ggplot(dat_combined,aes(x=factor(Samples),y=Genes.y/1000)) + 
  geom_boxplot(colour="#46bbc0",outlier.colour=NULL,outlier.shape=19,outlier.size=1)+
  scale_y_continuous(limits =  c(2.5,7), position = "right",
                     breaks = seq(3,7, by=1))  +
  theme_half_open(11, rel_small = 1) + 
  rremove("x.axis")+
  rremove("xlab") +
  rremove("x.text") +
  rremove("x.ticks") +
  rremove("legend") 

p2

aligned_plots <- align_plots(p1, p2, align="hv", axis="tblr")
FigS1a.KO <- ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])
FigS1a.KO



# Figure S1b. -----------------------------------------------------
rm(list = ls())
# Fig S1b.stats #####
data<-read_excel("Fig S1 Source Data.xlsx" ,sheet =  "Fig S1b bin stats")
plotDat<-melt(data,id.vals="NAME")
group<-factor(plotDat$Group,levels=c("High","Medium","Low"))
FigS1b.stats <- 
  ggplot(plotDat,aes(x=group,y=value))+
  geom_boxplot(outlier.shape=NA,aes(fill=group))+
  geom_jitter(alpha=0.1,shape=20,size=3,position=position_jitter(0.2))+
  facet_wrap(~variable,scale="free")+
  scale_fill_manual(values=c("#3081bc","#9fc8e2","#dfebf5")) +
  theme_bw()+theme(axis.line = element_line(colour = "black"),
                   panel.grid.major = element_blank(),
                   panel.grid.minor = element_blank(),
                   panel.background = element_blank())
FigS1b.stats

# Fig S1b.taxa #####
data <- read_excel("Fig S1 Source Data.xlsx", sheet = "Fig S1b bin taxa")
plotDat <- data %>% reshape2::melt() %>% 
  group_by(variable, Genus) %>% 
  summarise(sumed = sum(value), Phylum=unique(Phylum)) %>%
  mutate(freq=sumed/sum(sumed))

plotDat$Phylum <- factor(plotDat$Phylum, 
                         levels = c("Bacteroidetes","Proteobacteria","Actinobacteria",
                                    "Firmicutes","Spirochaetes","TM7"))
genus.level <- c("Prevotella","Porphyromonadaceae","Bacteroidales",
                 "Neisseria","Haemophilus","Moraxella","Campylobacter","Ralstonia","Pseudomonas","Lautropia",
                 "Rothia","Schaalia","Actinomyces",
                 "Streptococcus","Veillonella","Lachnospiraceae","Mogibacterium",
                 "Megasphaera","Selenomonas","Lactobacillus",
                 "Treponema","Saccharibacteria","Others","Unclassified")
plotDat$Genus[!plotDat$Genus  %in% genus.level] <- "Others"

plotDat$Genus <- factor(plotDat$Genus,levels = rev(genus.level))

mypalette <- c("#E25554","#E88D8D","#F2C3C3",
               "#EEA020","#F0AD4D","#F4BF76","#F7CF9E","#F6D9B9","#FBE9D8","#FCF2EA",
               "#49A9CE","#80CCE0","#BFE4EE",
               "#33A863","#4EB676","#77C392","#96CFA9",
               "#B4DCC1","#C5E3D0","#E2F1E6",
               "#D0BB8A","#968CBC","#B7B6B6","#E5E3E2")

FigS1b.taxa <- ggplot(plotDat) +
  geom_col(aes(x=variable, y=freq, fill=Genus))+
  scale_fill_manual(values = rev(mypalette))+
  theme_classic()
FigS1b.taxa