---
title: "Kean et al R code"
author: "Iain RL Kean"
date: "25/05/2021"
output: pdf_document
editor_options: 
  chunk_output_type: console
---
```{r environment, include=FALSE}

#Load packages into R environment, if package is not installed, initiates installation.

if (!require("tidyverse")){  
  install.packages("tidyverse")
  library(tidyverse)
}
if (!require("ggfortify")){
  install.packages("ggfortify")
  library(ggfortify)
}
if (!require("reshape2")){
  install.packages("reshape2")
  library(reshape2)
}
if (!require("gridExtra")){
  install.packages("gridExtra")
  library(gridExtra)
}
if (!require("vegan")){
  install.packages("vegan")
  library(vegan)
}
if (!require("conover.test")){
  install.packages("conover.test")
  library(conover.test)
}
if (!require("ggpubr")){
  install.packages("ggpubr")
  library(ggpubr)
}


```

```{r Create output directory, include=FALSE}
if(!dir.exists("~/Documents/Kean_2021")){
  dir.create(file.path("~/Documents/", "Kean_2021"))
}

```

```{r statistical layers, include=FALSE}

#Creates a dataframe from conover.test output to use with ggplot2 
make_conover_layer <- function(df=NULL , group= NULL , method = "bh", step.increase = 0.1, text.adjust = 0.02, include.ns = F) {
  if (!require("conover.test")){
    install.packages("conover.test")
    library(conover.test)
  }
  if (!require("tidyverse")){
    install.packages("tidyverse")
    library("tidyverse")
  }
  con.df<- as.data.frame(conover.test::conover.test( x = df, g = group, method = method, altp = T)[2:5])
  con.df <- con.df %>%  separate(comparisons, c("x","xend")," - ")
  con.df$y <- max(df)
  lev <- levels(as.factor(group))
  if (include.ns == T){
    con.df$sig <- ifelse(con.df$altP.adjusted <= 0.05, "*",round(con.df$altP.adjusted, digits = 3))
  }
  else{
    con.df$sig <- ifelse(con.df$altP.adjusted <= 0.05, "*","NS")
  }
  
  con.df$sig <- ifelse(con.df$altP.adjusted <= 0.01, "**",con.df$sig)
  con.df$sig <- ifelse(con.df$altP.adjusted <= 0.001, "***",con.df$sig)
  con.df$tx <- match(con.df$x, lev) + (match(con.df$xend,lev) - match(con.df$x,lev))*0.5
  if (include.ns == F){
    con.df <- subset(con.df, as.numeric(con.df$altP.adjusted) <= 0.05 )
  }
  if (nrow(con.df) >0){
  for (i in 1:nrow(con.df)){
    con.df$y[i] <- max(df) + (max(df)* i * step.increase )
    con.df$ty[i] <- con.df$y[i] + (con.df$y[i] * text.adjust)
  }
  }
  con.df$yend <- con.df$y 
  return(con.df)
}


```

# Add the below line to your ggplot code

+ annotate("segment", x = df$x, xend = df$xend, y = df$y, yend=df$yend)  + annotate("text", x = df$tx, y= df$ty, label = df$sig)

# Example
data("iris")
iris.layer <- make_conover_layer(df=iris$Sepal.Length, group = iris$Species)
iris.layer
ggplot(iris, mapping = aes( x = Species, y = Sepal.Length)) + geom_boxplot() + annotate(geom = "segment", x=iris.layer$x, xend = iris.layer$xend, y= iris.layer$y, yend=iris.layer$yend) + annotate( geom = "text", y = iris.layer$ty, x= iris.layer$tx, label = iris.layer$sig )


```{r Import bile acid data, include=FALSE}
#Read file into dataframe
ghici_full_3.df <- read.delim('~/Documents/Kean_2021/Kean_2021_data_table.csv', sep = ",", header = TRUE)
#Add column for new date ranges
ghici_full_3.df$DateRange2 <- ifelse(grepl("C", ghici_full_3.df$Sample_ID, ignore.case = T), "Control", "X")
ghici_full_3.df$DateRange2 <- ifelse( as.numeric(ghici_full_3.df$DayspostPICU) <3.5 & ghici_full_3.df$DayspostPICU != "Control", "d1-3",  ghici_full_3.df$DateRange2)
ghici_full_3.df$DateRange2 <- ifelse( as.numeric(ghici_full_3.df$DayspostPICU) >=3.5 & as.numeric(ghici_full_3.df$DayspostPICU) <=7.0, "d4-7", ghici_full_3.df$DateRange2)
ghici_full_3.df$DateRange2 <- ifelse( as.numeric(ghici_full_3.df$DayspostPICU) >7.0 ,"d8-10", ghici_full_3.df$DateRange2)
ghici_full_3.df$DateRange2 <- ifelse( is.na(ghici_full_3.df$DateRange2) ,  "d1-3", ghici_full_3.df$DateRange2)
ghici_full_3.df$DateRange2 <- ifelse(grepl("C", ghici_full_3.df$Sample_ID, ignore.case = T), "Control",ghici_full_3.df$DateRange2)


ghici_full_3.df$secBA <- ghici_full_3.df$Deoxycholic_Acid + ghici_full_3.df$Lithocholic_acid + ghici_full_3.df$Isolithocholic_Acid + ghici_full_3.df$Ursodeoxycholic_acid + ghici_full_3.df$X3a.Hydroxy.12_Ketolithocholic_Acid + ghici_full_3.df$X3a.Hydroxy.12_Ketolithocholic_Acid.1 + ghici_full_3.df$Taurohyocholic_Acid + ghici_full_3.df$Glycoursodeoxycholic_Acid

ghici_full_3.df$priBA <- ghici_full_3.df$Cholic_acid + ghici_full_3.df$Chenodeoxycholic_Acid + ghici_full_3.df$Glycochenodeoxycholic_Acid + ghici_full_3.df$Taurocholic_Acid

#Subset of dataframe with bile acids
by.df<- subset(ghici_full_3.df, ghici_full_3.df$Cholic_acid >0)[,c(1:32,216:218)]
by.df <-droplevels(by.df)



```



```{r Prepare taxanomic dataframes , include=FALSE}}
genus.df <-  subset(ghici_full_3.df, rowSums(ghici_full_3.df[,57:215]) > 0)[,c(34:216)]

colnames(ghici_full_3.df)[216]

nrow(genus.df)

genus.df$Shannon <- diversity(genus.df[,c(24:182)], index = "shannon")

genus.df$first_sample <- ifelse(genus.df$first_sample =="first_sample ","First",genus.df$first_sample)
genus.df$first_sample <- ifelse(genus.df$first_sample =="second_sample","Second",genus.df$first_sample)
genus.df$first_sample <- ifelse(genus.df$first_sample =="third_sample","Third",genus.df$first_sample)
genus.df$first_sample <- ifelse(genus.df$first_sample =="control","Control",genus.df$first_sample)

genus.df$infection <- ifelse(grepl(x = genus.df$diagnosis_group1 , pattern = "infection",ignore.case = T),"Infection",ifelse(grepl(x = genus.df$diagnosis_group1 , pattern = "control",ignore.case = T), "Control","Other"))

genus.df$diagnosis_group1 <- ifelse(genus.df$diagnosis_group1 =="resp_infection","Resp Infection",genus.df$diagnosis_group1)
genus.df$diagnosis_group1 <- ifelse(genus.df$diagnosis_group1 =="control","Control",genus.df$diagnosis_group1)
genus.df$diagnosis_group2 <- ifelse(genus.df$diagnosis_group2 =="lung","Resp Infection" , ifelse(genus.df$diagnosis_group2 =="control","Control","No Resp Infection"))

genus.df$amo <- ifelse( grepl("amo" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#penicillins
genus.df$aug <- ifelse( grepl("aug" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#penicillins #TAZOBACTAM
genus.df$azi <- ifelse( grepl("azi" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#macorlide
genus.df$benz <- ifelse( grepl("benz" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0) #penicillins
genus.df$cef <- ifelse( grepl("cef" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0) #cephalosporin
genus.df$cet <- ifelse( grepl("cet" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#cephalosporin
genus.df$chl <- ifelse( grepl("chl" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0) #Amphenicols
genus.df$cip <- ifelse( grepl("cip" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0) #quinolone
genus.df$clar <- ifelse( grepl("clar" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#macorlide
genus.df$clin <- ifelse( grepl("clin" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#Lincosamide
genus.df$coa <- ifelse( grepl("coa" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#penicillins #TAZOBACTAM
genus.df$ery <- ifelse( grepl("ery" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#macorlide
genus.df$flu <- ifelse( grepl("flu" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)
genus.df$gent <- ifelse( grepl("gent" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#Aminoglycoside
genus.df$met <- ifelse( grepl("met" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#penicillins
genus.df$pen <- ifelse( grepl("pen" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#penicillins
genus.df$pip <- ifelse( grepl("pip" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#penicillins
genus.df$tax <- ifelse( grepl("tax" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#TAZOBACTAM
genus.df$taz <- ifelse( grepl("taz" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#TAZOBACTAM
genus.df$tei <- ifelse( grepl("tei" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#glycopeptide
genus.df$tob <- ifelse( grepl("tob" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#Aminoglycoside
genus.df$van <- ifelse( grepl("van" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)#glycopeptide
genus.df$unknown <- ifelse( grepl("unknown" , genus.df$antibiotic_type, ignore.case = TRUE),  1, 0)
genus.df$abNum <- rowSums(genus.df[, 186:208])

genus.df$RAB2 <- genus.df$Alistipes+genus.df$Bacteroides+ genus.df$Bifidobacterium+  genus.df$Coprococcus+ genus.df$Desulfovibrio+  genus.df$Faecalibacterium+ genus.df$Parabacteroides+ genus.df$Roseburia+ genus.df$Ruminococcus+ genus.df$Subdoligranulum

# Infer location from ventilation free days
genus.df$Location <- ifelse( genus.df$group_patient =="patient" , ifelse( ((28 - as.numeric(genus.df$PICU_free_days) )- (as.numeric(genus.df$day_sample))) >=0, "PICU","General"), "Control")

genus.df$Location[6] <- "PICU"

genus.df$Location <- factor(genus.df$Location,levels = c("Control","PICU","General") )

droplevels(genus.df$Location)


#Create distance matrix of taxanomic data
tc.mat <- as.matrix.data.frame(genus.df[,24:182])
tc.dist <- vegdist(tc.mat)
tc.mds <- metaMDS(tc.dist, scale.=T)

#Add identifying data to NMDS dataframe
tc2.dat <- as.data.frame(tc.mds$points)
tc2.dat$Samlpe_num <- genus.df$Sample_num
tc2.dat$GHICI_number<- genus.df$GHICI_num
tc2.dat$DateRange2 <- genus.df$DateRange2



```

```{r  patient data}

weight.df <- subset(ghici_full_3.df, !is.na(ghici_full_3.df$weight_kg))
conover.test(x=weight.df$weight_kg, g = weight.df$Patient_0no1yes)
summary(subset(weight.df, weight.df$Patient_0no1yes ==1)$weight_kg)
summary(subset(weight.df, weight.df$Patient_0no1yes ==0)$weight_kg)
#ggplot(weight.df, mapping = aes(x = as.factor(Patient_0no1yes), y = weight_kg)) + geom_boxplot(notch = T)+ xlab("Patient No:0, Yes:1") + ylab("Weight (kg)")

age.df <- subset(ghici_full_3.df, !is.na(ghici_full_3.df$age_months))

conover.test(x=age.df$age_months, g = age.df$Patient_0no1yes)
summary(subset(age.df, age.df$Patient_0no1yes ==1)$age_months/12)
summary(subset(age.df, age.df$Patient_0no1yes ==0)$age_months/12)
#ggplot(age.df, mapping = aes(x = as.factor(Patient_0no1yes), y = age_months/12)) + geom_boxplot(notch = T) + xlab("Patient No:0, Yes:1")+ ylab("Age (years)")

summary(as.factor(ghici_full_3.df$Primary_Organ_failure))


```


```{r Figure 1}


pric <- make_conover_layer(df = (by.df$Cholic_acid+by.df$Chenodeoxycholic_Acid+by.df$Glycochenodeoxycholic_Acid+by.df$Taurocholic_Acid)/rowSums(by.df[,14:32]), group = by.df$DateRange2, method = "bh")

pri.plot <- ggplot(by.df, aes(x= DateRange2 , y = ((Cholic_acid+Chenodeoxycholic_Acid+Glycochenodeoxycholic_Acid+Taurocholic_Acid)/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  +  annotate("segment", x = pric$x, xend = pric$xend, y = pric$y, yend = pric$yend) + annotate( "text", x= pric$tx, y = pric$ty, label = pric$sig) + scale_y_continuous()  + ylab('RC of Primary Bile Acids') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none" ,axis.title.x=element_blank() , axis.text.y.left = element_text(angle = 90)) +labs(tag="a")#, axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.title.y = element_text(size = 15) ,  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) 

secc <- make_conover_layer(df=((by.df$Deoxycholic_Acid+by.df$Lithocholic_acid+by.df$Isolithocholic_Acid+by.df$Ursodeoxycholic_acid+by.df$`X3a.Hydroxy.12_Ketolithocholic_Acid`+by.df$`X3a.Hydroxy.12_Ketolithocholic_Acid.1`+by.df$Taurohyocholic_Acid+by.df$Glycoursodeoxycholic_Acid)/rowSums(by.df[,14:32])),group=by.df$DateRange2, method = "bh" )

sec.plot <- ggplot(by.df, aes(x= DateRange2 , y = ((Deoxycholic_Acid+Lithocholic_acid+Isolithocholic_Acid+Ursodeoxycholic_acid +`X3a.Hydroxy.12_Ketolithocholic_Acid`+`X3a.Hydroxy.12_Ketolithocholic_Acid.1`+Taurohyocholic_Acid+Glycoursodeoxycholic_Acid)/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  +  annotate("segment", x = secc$x, xend = secc$xend, y = secc$y, yend = secc$yend) + annotate( "text", x= secc$tx, y = secc$ty, label = secc$sig) + scale_y_continuous()  + ylab('RC of Secondary Bile Acids') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none" ,axis.title.x=element_blank() , axis.text.y.left = element_text(angle = 90)) +labs(tag="b") # ,axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.title.y = element_text(size = 15) , axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) 

crc <- make_conover_layer((by.df$Cholic_acid)/(by.df$Chenodeoxycholic_Acid), by.df$DateRange2, method = "bh")

chen_ratio.plot <- ggplot(by.df, aes(x= DateRange2 , y = (Cholic_acid/Chenodeoxycholic_Acid), fill=DateRange2)) + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE) + theme_bw() +  annotate("segment", x = crc$x, xend = crc$xend, y = crc$y, yend = crc$yend) + annotate( "text", x= crc$tx, y = crc$ty, label = crc$sig)+ scale_y_continuous()  + ylab('CA:CDCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none" ,axis.title.x=element_blank() , axis.text.y.left = element_text(angle = 90)) +labs(tag="c") #,axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title.y = element_text(size = 15)  ,axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))

ratioc <-make_conover_layer((by.df$Cholic_acid+by.df$Chenodeoxycholic_Acid+by.df$Glycochenodeoxycholic_Acid+by.df$Taurocholic_Acid)/(by.df$Deoxycholic_Acid+by.df$Lithocholic_acid+by.df$Isolithocholic_Acid+by.df$Ursodeoxycholic_acid+by.df$`X3a.Hydroxy.12_Ketolithocholic_Acid`+by.df$`X3a.Hydroxy.12_Ketolithocholic_Acid.1`+by.df$Taurohyocholic_Acid+by.df$Glycoursodeoxycholic_Acid),by.df$DateRange2, method = "bh")

ratio.plot <- ggplot(by.df, aes(x= DateRange2 , y = (Cholic_acid+Chenodeoxycholic_Acid+Glycochenodeoxycholic_Acid+ Taurocholic_Acid)/(Deoxycholic_Acid+ Lithocholic_acid+ Isolithocholic_Acid+Ursodeoxycholic_acid+`X3a.Hydroxy.12_Ketolithocholic_Acid`+`X3a.Hydroxy.12_Ketolithocholic_Acid.1`+Taurohyocholic_Acid+Glycoursodeoxycholic_Acid), fill=DateRange2)) + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE) + theme_bw()  +  annotate("segment", x = ratioc$x, xend = ratioc$xend, y = ratioc$y, yend = ratioc$yend) + annotate( "text", x= ratioc$tx, y = ratioc$ty, label = ratioc$sig)+ scale_y_continuous()  + ylab(expression(paste("Ratio of 1",degree," to 2",degree," Bile Acids"))) + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none" ,axis.title.x=element_blank() , axis.text.y.left = element_text(angle = 90))  +labs(tag="d") #,axis.title.x=element_blank(),  size = 15, hjust = 0.5),  axis.title.y = element_text(size = 15) ,   axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))

ratio.plot
cdc <- make_conover_layer((by.df$Cholic_acid)/(by.df$Deoxycholic_Acid), by.df$DateRange2, method ="bh")

chol_deox_ratio.plot <- c<- ggplot(by.df, aes(x= DateRange2 , y = (Cholic_acid/Deoxycholic_Acid), fill=DateRange2)) + theme_bw()+ geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)   +  annotate("segment", x = cdc$x, xend = cdc$xend, y = cdc$y, yend = cdc$yend) + annotate( "text", x= cdc$tx, y = cdc$ty, label = cdc$sig) + scale_y_continuous()  + ylab('CA:DCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none" ,axis.title.x=element_blank() , axis.text.y.left = element_text(angle = 90)) +labs(tag="e") #,axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.title.y = element_text(size = 15) ,   axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) 

clc <- make_conover_layer((by.df$Chenodeoxycholic_Acid)/(by.df$Lithocholic_acid+by.df$Isolithocholic_Acid), by.df$DateRange2, method = "by")

chen_lith_ratio.plot <- ggplot(by.df, aes(x= DateRange2 , y = (Chenodeoxycholic_Acid/(Lithocholic_acid+Isolithocholic_Acid)), fill=DateRange2))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  +  annotate("segment", x = clc$x, xend = clc$xend, y = clc$y, yend = clc$yend) + annotate( "text", x= clc$tx, y = clc$ty, label = clc$sig) + scale_y_continuous()  + ylab('CDCA:LCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none" ,axis.title.x=element_blank() , axis.text.y.left = element_text(angle = 90)) +labs(tag="f") # ,axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.title.y = element_text(size = 15) , axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="f")

fig1 <- grid.arrange(pri.plot, sec.plot ,chen_ratio.plot, ratio.plot,  chol_deox_ratio.plot, chen_lith_ratio.plot, nrow = 2)

ggsave( "~/Documents/Kean_2021/Kean_2021_Figure_1.png",plot = fig1,  height = 120 , width = 170, units = "mm", dpi = 300)

```



```{r Figure S1}

ba1c <- make_conover_layer(df =by.df$`X3a.Hydroxy.6.7.DiketoCholanic_Acid`/rowSums(by.df[,14:32]), group= by.df$DateRange2 ,method= "bh")

figS1.1 <- ggplot(by.df, aes(x= DateRange2 , y = (`X3a.Hydroxy.6.7.DiketoCholanic_Acid`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba1c$x, xend = ba1c$xend, y = ba1c$y, yend=ba1c$yend)  + annotate("text", x = ba1c$tx, y= ba1c$ty, label = ba1c$sig) + scale_y_continuous()  + ylab(expression(paste('RC of 3',alpha,'-OH,6,7-diKClA'))) + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="a")

ba2c <- make_conover_layer(df =by.df$Taurohyocholic_Acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.2 <- ggplot(by.df, aes(x= DateRange2 , y = (Taurohyocholic_Acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba2c$x, xend = ba2c$xend, y = ba2c$y, yend=ba2c$yend)  + annotate("text", x = ba2c$tx, y= ba2c$ty, label = ba2c$sig)  + scale_y_continuous()  + ylab('RC of THCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="b")

ba3c <- make_conover_layer(df =by.df$Glycoursodeoxycholic_Acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.3 <- ggplot(by.df, aes(x= DateRange2 , y = (Glycoursodeoxycholic_Acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  +  annotate("segment", x = ba3c$x, xend = ba3c$xend, y = ba3c$y, yend=ba3c$yend)  + annotate("text", x = ba3c$tx, y= ba3c$ty, label = ba3c$sig) + scale_y_continuous()  + ylab('RC of GUDCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="c")

ba4c <- make_conover_layer(df =by.df$`X3_Dehydrocholic_Acid`/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.4 <- ggplot(by.df, aes(x= DateRange2 , y = (`X3_Dehydrocholic_Acid`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw()+ geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba4c$x, xend = ba4c$xend, y = ba4c$y, yend=ba4c$yend)  + annotate("text", x = ba4c$tx, y= ba4c$ty, label = ba4c$sig)  + scale_y_continuous()  + ylab('RC of 3-DHCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="d")

ba5c <- make_conover_layer(df =by.df$`X3.12.Diketocholanic_Acid`/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.5 <- ggplot(by.df, aes(x= DateRange2 , y = (`X3.12.Diketocholanic_Acid`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba5c$x, xend = ba5c$xend, y = ba5c$y, yend=ba5c$yend)  + annotate("text", x = ba5c$tx, y= ba5c$ty, label = ba5c$sig) + scale_y_continuous()  + ylab('RC of 3,12-diKClA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="e")

ba6c <- make_conover_layer(df =by.df$Ursodeoxycholic_acid/by.df$Chenodeoxycholic_Acid, group= by.df$DateRange2 , method ="bh")

 ggplot(by.df, aes(x= DateRange2 , y = (Ursodeoxycholic_acid/Chenodeoxycholic_Acid) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE) + annotate("segment", x = ba6c$x, xend = ba6c$xend, y = ba6c$y, yend=ba6c$yend)  + annotate("text", x = ba6c$tx, y= ba6c$ty, label = ba6c$sig) + scale_y_continuous()  + ylab('UDCA:CDCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="f")

ba6c <- make_conover_layer(df =by.df$Ursodeoxycholic_acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.6 <- ggplot(by.df, aes(x= DateRange2 , y = (Ursodeoxycholic_acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE) + annotate("segment", x = ba6c$x, xend = ba6c$xend, y = ba6c$y, yend=ba6c$yend)  + annotate("text", x = ba6c$tx, y= ba6c$ty, label = ba6c$sig) + scale_y_continuous()  + ylab('RC of UDCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="f")

ba7c <- make_conover_layer(df =by.df$`X12_Dehydrocholic_Acid`/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.7 <- ggplot(by.df, aes(x= DateRange2 , y = (`X12_Dehydrocholic_Acid`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba7c$x, xend = ba7c$xend, y = ba7c$y, yend=ba7c$yend)  + annotate("text", x = ba7c$tx, y= ba7c$ty, label = ba7c$sig)  + scale_y_continuous()  + ylab('RC of 12-DHCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="g")

ba8c <- make_conover_layer(df =by.df$Glycochenodeoxycholic_Acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.8 <- ggplot(by.df, aes(x= DateRange2 , y = (Glycochenodeoxycholic_Acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba8c$x, xend = ba8c$xend, y = ba8c$y, yend=ba8c$yend)  + annotate("text", x = ba8c$tx, y= ba8c$ty, label = ba8c$sig) + scale_y_continuous()  + ylab('RC of GCDCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="h")

ba9c <- make_conover_layer(df =by.df$`X3a.Hydroxy.12_Ketolithocholic_Acid`/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.9 <- ggplot(by.df, aes(x= DateRange2 , y = (`X3a.Hydroxy.12_Ketolithocholic_Acid`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE) + annotate("segment", x = ba9c$x, xend = ba9c$xend, y = ba9c$y, yend=ba9c$yend)  + annotate("text", x = ba9c$tx, y= ba9c$ty, label = ba9c$sig)  + scale_y_continuous()  + ylab(expression(paste('RC of 3',alpha,'-H,12-KLCA'))) + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="i")

ba10c <- make_conover_layer(df =by.df$Isolithocholic_Acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.10 <- ggplot(by.df, aes(x= DateRange2 , y = (Isolithocholic_Acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba10c$x, xend = ba10c$xend, y = ba10c$y, yend=ba10c$yend)  + annotate("text", x = ba10c$tx, y= ba10c$ty, label = ba10c$sig) + scale_y_continuous()  + ylab('RC of ILCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="j")

ba11c <- make_conover_layer(df =by.df$Chenodeoxycholic_Acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.11 <- ggplot(by.df, aes(x= DateRange2 ,y = (Chenodeoxycholic_Acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba11c$x, xend = ba11c$xend, y = ba11c$y, yend=ba11c$yend)  + annotate("text", x = ba11c$tx, y= ba11c$ty, label = ba11c$sig) + scale_y_continuous()  + ylab('RC of CDCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 10, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="k")

ba12c <- make_conover_layer(df =by.df$`X5ß.Cholanic_Acid.3ß._12a.diol`/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.12 <- ggplot(by.df, aes(x= DateRange2 , y = (`X5ß.Cholanic_Acid.3ß._12a.diol`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba12c$x, xend = ba12c$xend, y = ba12c$y, yend=ba12c$yend)  + annotate("text", x = ba12c$tx, y= ba12c$ty, label = ba12c$sig)+ ylab(expression(paste('RC of 5',beta,'-ClA,3',beta,',12',alpha,'-diol')))  + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 10, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="l")

ba13c <- make_conover_layer(df =by.df$`X3.Ketocholanic_Acid`/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.13 <- ggplot(by.df, aes(x= DateRange2 , y = (`X3.Ketocholanic_Acid`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba13c$x, xend = ba13c$xend, y = ba13c$y, yend=ba13c$yend)  + annotate("text", x = ba13c$tx, y= ba13c$ty, label = ba13c$sig)  + scale_y_continuous() + ylab('RC of 3-KClA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="m")

ba14c <- make_conover_layer(df =by.df$Lithocholic_acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.14 <- ggplot(by.df, aes(x= DateRange2 , y = (Lithocholic_acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba14c$x, xend = ba14c$xend, y = ba14c$y, yend=ba14c$yend)  + annotate("text", x = ba14c$tx, y= ba14c$ty, label = ba14c$sig) + scale_y_continuous()  + ylab('RC of LCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="n")

ba15c <- make_conover_layer(df =by.df$Deoxycholic_Acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.15 <- ggplot(by.df, aes(x= DateRange2 , y = (Deoxycholic_Acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba15c$x, xend = ba15c$xend, y = ba15c$y, yend=ba15c$yend)  + annotate("text", x = ba15c$tx, y= ba15c$ty, label = ba15c$sig)  + scale_y_continuous()  + ylab('RC of DCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="o")

ba16c <- make_conover_layer(df =by.df$Cholic_acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.16 <- ggplot(by.df, aes(x= DateRange2 , y = (Cholic_acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba16c$x, xend = ba16c$xend, y = ba16c$y, yend=ba16c$yend)  + annotate("text", x = ba16c$tx, y= ba16c$ty, label = ba16c$sig) + scale_y_continuous()  + ylab('RC of CA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="p")
 
ba17c <- make_conover_layer(df =by.df$`X23.nor.5b.Cholanic_Acid.3a._12a.diol`/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.17 <- ggplot(by.df, aes(x= DateRange2 , y = (`X23.nor.5b.Cholanic_Acid.3a._12a.diol`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw()+ geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba17c$x, xend = ba17c$xend, y = ba17c$y, yend=ba17c$yend)  + annotate("text", x = ba17c$tx, y= ba17c$ty, label = ba17c$sig) + scale_y_continuous()  + ylab(expression(paste('RC of 23 nor-,5',beta,'-ClA,3',alpha,',12',alpha,'-diol'))) + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="q")

ba18c <- make_conover_layer(df =by.df$`X3a.Hydroxy.12_Ketolithocholic_Acid.1`/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.18 <-ggplot(by.df, aes(x= DateRange2 , y = (`X3a.Hydroxy.12_Ketolithocholic_Acid.1`/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba18c$x, xend = ba18c$xend, y = ba18c$y, yend=ba18c$yend)  + annotate("text", x = ba18c$tx, y= ba18c$ty, label = ba18c$sig) + scale_y_continuous()  + ylab(expression(paste('RC of 3',alpha,'-H,12-KLCA 2'))) + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 10, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="r") 

ba19c <- make_conover_layer(df =by.df$Taurocholic_Acid/rowSums(by.df[,14:32]), group= by.df$DateRange2 , method ="bh")

figS1.19 <- ggplot(by.df, aes(x= DateRange2 , y = (Taurocholic_Acid/rowSums(by.df[,14:32], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = ba19c$x, xend = ba19c$xend, y = ba19c$y, yend=ba19c$yend)  + annotate("text", x = ba19c$tx, y= ba19c$ty, label = ba19c$sig) + scale_y_continuous()  + ylab('RC of TCA') + geom_jitter (aes( colour = DateRange2), width =0.2) + scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2',"forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange',"darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5),  axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="s")


figS1 <- grid.arrange(figS1.1,figS1.2, figS1.3, figS1.4, figS1.5, figS1.6, figS1.7, figS1.8,figS1.9, figS1.10, figS1.11, figS1.12, figS1.13, figS1.14, figS1.15, figS1.16, figS1.17, figS1.18, figS1.19, nrow=4)


ggsave("~/Documents/Kean_2021/Kean_2021_Figure_S1.eps", plot = figS1, width =470, height = 280 ,units = "mm", dpi = 300)

```


```{r Figure 2}

#Shanon index
fig2a.con <- make_conover_layer(group = genus.df$DateRange2, df = genus.df$Shannon, method = "bh")

fig2.a <- ggplot( genus.df, mapping = aes(x = DateRange2, y = Shannon, fill = DateRange2))+ theme_bw()+ geom_violin(draw_quantiles = c(0.25,0.5,0.75))+ scale_fill_manual( values = c('lightblue','firebrick','goldenrod',"forestgreen",'firebrick','goldenrod',"forestgreen"))+ annotate("segment", x= fig2a.con$x, xend = fig2a.con$xend, y= fig2a.con$y, yend = fig2a.con$yend) + annotate("text" , x = fig2a.con$tx, y = fig2a.con$ty, label = fig2a.con $sig)+labs(tag="a")+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90, hjust = 0.5) , axis.text.x.bottom = element_text( hjust = 0.5))

# Diversity
fig2.b <-ggplot( tc2.dat, mapping = aes( x= MDS1, y = MDS2, colour = DateRange2) ) +theme_bw()+ geom_point() +stat_ellipse(level =0.6) + scale_colour_manual( values = c('darkblue','firebrick','goldenrod3',"forestgreen")) + theme(legend.position="bottom", axis.text.y.left = element_text(angle = 90, hjust = 0.5), axis.text.x= element_text( hjust = 0.7), legend.text = element_text(size = 7), legend.title = element_blank()  ) +labs(tag="b", colour ="Group")

fig2 <- grid.arrange(fig2.a, fig2.b ,nrow= 2)

ggsave(filename = "~/Documents/Kean_2021/Kean_2021_Figure_2.eps", plot = fig2, width = 85, height = 170, units = "mm", dpi = 300)

```

```{r Figure 3}

#Truncated family level taxonomy derived from SILVA138.1 database.

family.df <- read.delim('~/Documents/Kean_2021/Kean_2021_family_taxonomy.csv', sep=",")
family.df$DateRange2 <- ifelse( family.df$day_sample >=1 & family.df$day_sample <=3, "d1-3", "x")
family.df$DateRange2 <- ifelse( family.df$day_sample >=4 & family.df$day_sample <=7, "d4-7", family.df$DateRange2)
family.df$DateRange2 <- ifelse( family.df$day_sample >7 ,"d8-10", family.df$DateRange2)
family.df$DateRange2 <- ifelse( family.df$day_sample <=0, "Control", family.df$DateRange2)

family.df$infection <- genus.df$infection

#melt faily dataframe for plotting
fam.melt = melt(family.df[,c(1:2,27,3:4,6:26)])

fig3 <- ggplot( data = fam.melt, mapping = aes( x = GHICI_num, y = value, fill = variable)) + geom_col( position="fill", colour = "grey50" ) +theme_bw()+ theme(legend.title = element_text( face= "bold", size = 10), axis.title.x=element_blank() ,legend.text = element_text(size = 6,  face="italic"), axis.text.x=element_text( angle = 300, hjust = 0, vjust = 0.5, size = 5.5), axis.text.y = element_text(angle = 90, hjust = 0.5), legend.key.size = unit(0.5,"line")) + facet_grid(~ DateRange2, scales= "free_x", space = "free_x", )+ ylab("Relative Proportion") + scale_fill_manual(name = "Family" ,  values=c('#e31a1c','#fb9a99', '#a6cee3','#1f78b4', '#fdbf6f','#ff7f00', '#33a02c','#b2df8a', '#cab2d6' ,'#6a3d9a', "grey25",'#6a3d9a','#cab2d6' ,'#ff7f00','#fdbf6f', '#e31a1c','#fb9a99' ,'#33a02c','#b2df8a','#1f78b4','#a6cee3')) + scale_y_continuous(expand = c(0,0)) + guides(fill = guide_legend(ncol = 1)) 

fig3
ggsave(filename = "~/Documents/Kean_2021/Kean_2021_Figure_3.eps", plot = fig3, width = 250, height = 150, units = "mm", dpi = 300)

```

```{r Figure 4}

raballc <-make_conover_layer(df = genus.df$RAB2/rowSums(genus.df[,24:182]), group= genus.df$DateRange2 , method = "bh" )

fig4.a <- ggplot(genus.df, aes(x= DateRange2 , y = ((RAB2)/rowSums(genus.df[,24:182], na.rm= TRUE) ), fill = DateRange2 ) ) + theme_bw()+ geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = raballc$x, xend = raballc$xend, y = raballc$y, yend=raballc$yend)  + annotate("text", x = raballc$tx, y= raballc$ty, label = raballc$sig) + scale_y_continuous()  + ylab('Proportion of RAB') + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen','firebrick', 'goldenrod2', "forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen','red','darkorange', "darkgreen"))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90, hjust = 0.5), axis.text.x.bottom = element_text( hjust = 0.5))  +labs(tag="a") 


c.l <-make_conover_layer(df = family.df$Lachnospiraceae/rowSums(family.df[,6:26]), group= family.df$DateRange2  , method = "bh" )

fig4.b <- ggplot(family.df, aes(x= DateRange2 , y = (`Lachnospiraceae`/rowSums(family.df[,6:26], na.rm= TRUE)) , fill = DateRange2 ))+ theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE) + annotate("segment", x = c.l$x, xend = c.l$xend, y = c.l$y, yend=c.l$yend)  + annotate("text", x = c.l$tx, y= c.l$ty, label = c.l$sig) + scale_y_continuous() + xlab('Days Post PICU Admission') + ylab(expression('Prop. of '~italic('Lachnospiraceae'))) + geom_jitter (aes( colour = DateRange2))+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2', "forestgreen")) + scale_colour_manual(values=c('darkblue','red','darkorange', "darkgreen")) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  hjust = 0.5) , axis.text.x.bottom = element_text( hjust = 0.5)) +labs(tag="b")

fig4 <- grid.arrange(fig4.a , fig4.b, nrow =2)

ggsave("~/Documents/Kean_2021/Kean_2021_Figure_4.eps", plot= fig4, width = 85, height = 120, units = "mm", dpi = 300)

```

```{r Figur S2}

loc.con <- make_conover_layer(group = genus.df$first_sample, df = (genus.df$RAB2/rowSums(genus.df[,24:182])), method = "bh")

figS2.a <-ggplot(genus.df, aes(x= first_sample , y = (`RAB2`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = first_sample )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = loc.con$x, xend = loc.con$xend, y = loc.con$y, yend=loc.con$yend)  + annotate("text", x = loc.con$tx, y= loc.con$ty, label = loc.con$sig) + scale_y_continuous() + xlab('Sample number') + ylab(expression('Proportion of RAB')) + geom_jitter (aes( colour = DateRange2, shape = Location))+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 13, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="a") 

loc.con <- make_conover_layer(group = genus.df$diagnosis_group1, df = (genus.df$RAB2/rowSums(genus.df[,24:182])), method = "bh")

figS2.b <- ggplot(genus.df, aes(x= diagnosis_group1 , y = (`RAB2`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = diagnosis_group1 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = loc.con$x, xend = loc.con$xend, y = loc.con$y, yend=loc.con$yend)  + annotate("text", x = loc.con$tx, y= loc.con$ty, label = loc.con$sig) + scale_y_continuous() + xlab('Primary Diagnosis') + ylab(expression('Proportion of RAB')) + geom_jitter (aes( colour = diagnosis_group1))+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen',"darkblue","pink")) + scale_colour_manual(values=c('darkblue','red',"darkorange",'darkgreen',"royalblue","purple"))+ theme(legend.position = "none", axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 10, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="b")

loc.con <- make_conover_layer(group = genus.df$diagnosis_group1, df = (family.df$Ruminococcaceae/rowSums(family.df[,24:44])), method = "bh")

figS2.c <- ggplot(family.df, aes(x= diagnosis_group1 , y = (`Ruminococcaceae`/rowSums(family.df[,24:44], na.rm= TRUE)) , fill = diagnosis_group1 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = loc.con$x, xend = loc.con$xend, y = loc.con$y, yend=loc.con$yend)  + annotate("text", x = loc.con$tx, y= loc.con$ty, label = loc.con$sig) + scale_y_continuous() + xlab('Primary Diagnosis') + ylab(expression('Proportion of'~italic('Lachnospiraceae'))) + geom_jitter (aes( colour = diagnosis_group1))+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen','darkblue','pink')) + scale_colour_manual(values=c('darkblue','red',"darkorange",'darkgreen','royalblue','purple'))+ theme(legend.position = "none", axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 10, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="c") 

loc.con <- make_conover_layer(group = family.df$diagnosis_group2, df = (family.df$Lachnospiraceae/rowSums(family.df[,24:44])), method = "bh")

figS2.d <- ggplot(family.df, aes(x=diagnosis_group2 , y = (`Lachnospiraceae`/rowSums(family.df[,24:44], na.rm= TRUE))  , fill = diagnosis_group2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = loc.con$x, xend = loc.con$xend, y = loc.con$y, yend=loc.con$yend)  + annotate("text", x = loc.con$tx, y= loc.con$ty, label = loc.con$sig) + scale_y_continuous() + xlab('Respiratory Infection') + ylab(expression('Proportion of'~italic('Lachnospiraceae'))) + geom_jitter (aes( colour = diagnosis_group2))+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen','goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','orange','green'))+ theme(legend.position = "none", axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 13, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="d") 

loc.con <- make_conover_layer(group = family.df$infection, df = (family.df$Lachnospiraceae/rowSums(family.df[,24:44])), method = "bh")

figS2.e <- ggplot(family.df, aes(x= infection , y = (Lachnospiraceae/rowSums(family.df[,24:44], na.rm= TRUE)) , fill = infection )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = loc.con$x, xend = loc.con$xend, y = loc.con$y, yend=loc.con$yend)  + annotate("text", x = loc.con$tx, y= loc.con$ty, label = loc.con$sig) + scale_y_continuous() + xlab('Infection / Non-Infection') + ylab(expression('Proportion of'~italic('Lachnospiraceae'))) + geom_jitter (aes( colour = infection))+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen','goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','orange','green'))+ theme(legend.position = "none", axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 13, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="e") 

figS2 <- grid.arrange(figS2.a, figS2.b , figS2.c, figS2.d, figS2.e, nrow=2)

ggsave("~/Documents/Kean_2021/Kean_2021_Figure_S2.eps", plot = figS2 ,width = 500,height =283, units = "mm", dpi=300 )


```

```{r Figure S3}

rb1c <-make_conover_layer(df = genus.df$Alistipes/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )

figS3.a <- ggplot(genus.df, aes(x= DateRange2 , y = (`Alistipes`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = rb1c$x, xend = rb1c$xend, y = rb1c$y, yend=rb1c$yend)  + annotate("text", x = rb1c$tx, y= rb1c$ty, label = rb1c$sig) + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Alistipes'))) + geom_jitter (aes( colour = DateRange2))+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen')) + theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title = element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="a") 

rb3c <-make_conover_layer(df = genus.df$Bacteroides/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )

figS3.b <- ggplot(genus.df, aes(x= DateRange2 , y = (`Bacteroides`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = rb3c$x, xend = rb3c$xend, y = rb3c$y, yend=rb3c$yend)  + annotate("text", x = rb3c$tx, y= rb3c$ty, label = rb3c$sig)  + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Bacteroides'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title.y = element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="b")

rb4c <-make_conover_layer(df = genus.df$Bifidobacterium/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )

figS3.c <- ggplot(genus.df, aes(x= DateRange2 , y = (`Bifidobacterium`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = rb4c$x, xend = rb4c$xend, y = rb4c$y, yend=rb4c$yend)  + annotate("text", x = rb4c$tx, y= rb4c$ty, label = rb4c$sig) + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Bifidobacterium'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title.y = element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="c")

rb7c <-make_conover_layer(df = genus.df$Coprococcus/rowSums(genus.df[,24:182]), group= genus.df$DateRange2 , method ="bh" )

figS3.d <- ggplot(genus.df, aes(x= DateRange2 , y = (`Coprococcus`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  +  annotate("segment", x = rb7c$x, xend = rb7c$xend, y = rb7c$y, yend=rb7c$yend)  + annotate("text", x = rb7c$tx, y= rb7c$ty, label = rb7c$sig) + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Coprococcus'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title.y = element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="d")

rb8c <-make_conover_layer(df = genus.df$Desulfovibrio/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )

figS3.e <- ggplot(genus.df, aes(x= DateRange2 , y = (`Desulfovibrio`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = rb8c$x, xend = rb8c$xend, y = rb8c$y, yend=rb8c$yend)  + annotate("text", x = rb8c$tx, y= rb8c$ty, label = rb8c$sig) + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Desulfovibrio'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="e")

rb11c <-make_conover_layer(df = genus.df$Faecalibacterium/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )

figS3.f <- ggplot(genus.df, aes(x= DateRange2 , y = (Faecalibacterium/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = rb11c$x, xend = rb11c$xend, y = rb11c$y, yend=rb11c$yend)  + annotate("text", x = rb11c$tx, y= rb11c$ty, label = rb11c$sig) + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Faecalibacterium'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15)) +labs(tag="f") #g


rb14c <-make_conover_layer(df = genus.df$Parabacteroides/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )
figS3.g <- ggplot(genus.df, aes(x= DateRange2 , y = (`Parabacteroides`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = rb14c$x, xend = rb14c$xend, y = rb14c$y, yend=rb14c$yend)  + annotate("text", x = rb14c$tx, y= rb14c$ty, label = rb14c$sig) + scale_y_continuous() + ylab(expression('Proportion of '~italic('Parabacteroides'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="g")

rb16c <-make_conover_layer(df = genus.df$Roseburia/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )
figS3.h <- ggplot(genus.df, aes(x= DateRange2 , y = (`Roseburia`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = rb16c$x, xend = rb16c$xend, y = rb16c$y, yend=rb16c$yend)  + annotate("text", x = rb16c$tx, y= rb16c$ty, label = rb16c$sig) + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Roseburia'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="h") #h l

rb17c <-make_conover_layer(df = (genus.df$Ruminococcus)/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )
figS3.i <- ggplot(genus.df, aes(x= DateRange2 , y = ((Ruminococcus) /rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE)  + annotate("segment", x = rb17c$x, xend = rb17c$xend, y = rb17c$y, yend=rb17c$yend)  + annotate("text", x = rb17c$tx, y= rb17c$ty, label = rb17c$sig)  + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Ruminococcus'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="o") +labs(tag="i")

rb18c <-make_conover_layer(df = genus.df$Subdoligranulum/rowSums(genus.df[,24:182]), group= genus.df$DateRange2  )

figS3.j <- ggplot(genus.df, aes(x= DateRange2 , y = (`Subdoligranulum`/rowSums(genus.df[,24:182], na.rm= TRUE)) , fill = DateRange2 )) + theme_bw() + geom_boxplot( outlier.colour='black', outlier.shape = 'X', notch = TRUE) + annotate("segment", x = rb18c$x, xend = rb18c$xend, y = rb18c$y, yend=rb18c$yend)  + annotate("text", x = rb18c$tx, y= rb18c$ty, label = rb18c$sig) + scale_y_continuous()  + ylab(expression('Proportion of '~italic('Subdoligranulum'))) + geom_jitter (aes( colour = DateRange2), width = 0.2)+ scale_fill_manual(values=c('lightblue','firebrick', 'goldenrod2','forestgreen')) + scale_colour_manual(values=c('darkblue','red','darkorange','darkgreen'))+ theme(legend.position = "none",axis.title.x=element_blank(), axis.text.y.left = element_text(angle = 90,  size = 15, hjust = 0.5) , axis.title= element_text(size = 15), axis.text.x.bottom = element_text( size = 15, hjust = 0.5), plot.tag = element_text(size = 15))  +labs(tag="j")


figS3 <- grid.arrange(figS3.a, figS3.b, figS3.c, figS3.d, figS3.e, figS3.f, figS3.g, figS3.h, figS3.i, figS3.j, nrow = 3)

ggsave("~/Documents/Kean_2021/Kean_2021_Figure_S3.eps", plot = figS3, width = 450, height = 250, units = "mm", dpi = 300)

```

```{r Table S3}
both.df <- subset(ghici_full_3.df, (ghici_full_3.df$Cholic_acid >0))
both.df <- subset(both.df, both.df$Bacteroides >0 )

both_pat.df <- subset(both.df, both.df$Patient_0no1yes ==1)

tableS3.df <- data.frame(Genus = character(), `p-value (tLCA)` = numeric(), `rho[15] (tLCA)` = numeric(), `p-value (DCA)` = numeric(), `rho[15] (DCA)` = numeric() )

for (i in 57:215) {
  tlca <- cor.test(x =((both_pat.df[,i])/rowSums(both_pat.df[,57:215])), y = ((both_pat.df$Lithocholic_acid + both_pat.df$Isolithocholic_Acid)/rowSums( both_pat.df[,14:32])), method = "spearman", exact = F)
  dca <- cor.test(x =((both_pat.df[,i])/rowSums(both_pat.df[,57:215])), y = ((both_pat.df$Deoxycholic_Acid)/rowSums(both_pat.df[,14:32])), method = "spearman", exact = F)
  if (( tlca$p.value <= 0.05 && !is.na(tlca$p.value)) || ( dca$p.value <= 0.05 && !is.na(dca$p.value)) ) {
    tableS3.df %>% add_row( Genus = colnames(both_pat.df)[i], p.value..tLCA. = tlca$p.value, rho.15...tLCA. =tlca$estimate , p.value..DCA. = dca$p.value, rho.15...DCA. = dca$estimate) -> tableS3.df
  }
}


write_delim(x = tableS3.df, file = "~/Documents/Kean_2021/Kean_2021_Table_S3.csv", delim = ",")

```



```{r Figure S4}

#Set up dataframe for lachnospiraceae identified in genus.df

genus.df$group_patient

lachno.df <- genus.df[,c(1,4)]

lachno.df$Agathobacter <- genus.df$Agathobacter
lachno.df$Anaerobutyricum.hallii.group <- genus.df$Eubacterium.hallii.group
lachno.df$Anaerocolumna <- genus.df$Anaerocolumna
lachno.df$Anaerostignum <- genus.df$Anaerostignum
lachno.df$Anaerostipes <- genus.df$Anaerostipes
lachno.df$Blautia <- genus.df$Blautia
lachno.df$Butyrivibrio <- genus.df$Butyrivibrio
lachno.df$Coprococcus <- genus.df$Coprococcus
lachno.df$Dorea <- genus.df$Dorea
lachno.df$Eisenbergiella <- genus.df$Eisenbergiella
lachno.df$Fusicatenibacter <- genus.df$Fusicatenibacter
genus.df$Howardella -> lachno.df$Howardella
lachno.df$Hungatella <- genus.df$Hungatella
genus.df$Lachnoanaerobaculum -> lachno.df$Lachnoanaerobaculum
lachno.df$Lachnoclostridium <- genus.df$Lachnoclostridium
genus.df$Lachnospira -> lachno.df$Lachnospira
lachno.df$Mediterranibacter.gnavus.group <- genus.df$Ruminococcus.gnavus.group
lachno.df$Mediterranibacter.torques.group <- genus.df$Ruminococcus.torques.group
genus.df$Oribacterium -> lachno.df$Oribacterium
genus.df$Pseudobutyrivibrio -> lachno.df$Pseudobutyrivibrio
genus.df$Roseburia -> lachno.df$Roseburia
genus.df$Stomatobaculum -> lachno.df$Stomatobaculum
lachno.df$Tyzzerella <- genus.df$Tyzzerella

lachno.df$DateRange2 <- genus.df$DateRange2

#Reshape lachno.df for graphing
f8b.m <- melt(lachno.df)

#Change formatting
f8b.m$variable <- as.character(f8b.m$variable)

f8b.m$variable <- ifelse(f8b.m$variable=="Anaerobutyricum.hallii.group", "Anaerobutyricum hallii group", f8b.m$variable)
f8b.m$variable <- ifelse(f8b.m$variable=="Ruminococcus.gnavus.group", "Mediterraneibacter gnavus group", f8b.m$variable)
f8b.m$variable <- ifelse(f8b.m$variable=="Ruminococcus.torques.group", "Mediterraneibacter torques group", f8b.m$variable)

figS4 <- ggplot(f8b.m, mappin = aes(x=Sample_num, y = value, fill = variable))+ geom_col() +theme_bw() + ylab("Raw Count")+ xlab("Sample")+theme(legend.text = element_text(face="italic"), axis.text.x = element_text( angle = 300, hjust = 0, face = "bold"))+ facet_grid(~ DateRange2, scales= "free_x", space = "free_x") + scale_fill_manual(name ="Genus", values=c('royalblue','#e31a1c','firebrick' ,'#fdbf6f','#cab2d6' ,'magenta','#ff7f00','#a6cee3' ,'#6a3d9a','#ffff99','#b15928',"pink","pink2",'#b2df8a', '#33a02c','#fb9a99' ,'#1f78b4','#fdbf6f','#ff7f00','#cab2d6','skyblue','#e31a1c', "pink", " purple"))+ guides(fill = guide_legend(ncol = 1))

ggsave(filename = "~/Documents/Kean_2021/Kean_2021_Figure_S4.eps", plot = figS4, width = 450, height = 250, units = "mm", dpi = 300)

```

```{r Table S4}

d8.df <- subset(genus.df, genus.df$DateRange2 =="d8-10")

lachno.list <- list("Agathobacter", "Anaerobutyricum.hallii.group", "Anaerocolumna", "Anaerostignum", "Anaerostipes", "Blautia", "Butyrivibrio", "Coprococcus", "Dorea", "Eisenbergiella",  "Fusicatenibacter", "Howardella", "Hungatella", "Lachnoanaerobaculum", "Lachnoclostridium", "Lachnospira", "Mediterraneibacter.gnavus.group", "Mediterraneibacter.torques.group", "Oribacterium", "Pseudobutyrivibrio", "Roseburia", "Stomatobaculum" ,"Tyzzerella")

colnames(d8.df)[49] <- "Anaerobutyricum.hallii.group"
colnames(d8.df)[40] <- "Mediterraneibacter.torques.group"
colnames(d8.df)[93] <- "Mediterraneibacter.gnavus.group"

hallii_conf.df <- data.frame(Lac = character(),Comp = character(),S_cor=double(),S_p_val=double())
for (i in lachno.list){

  if (sum(d8.df[,i]) >=10){
    for (x in  25:182) {
      if (median(d8.df[,x]) >=5){
        s <- cor.test(d8.df[,colnames(d8.df)[x]]/rowSums(d8.df[24:182]), d8.df[,i]/rowSums(d8.df[24:182]), method = "spearman", exact  =F)
        #print(s)
        if ( s$p.value <= 0.1 ) {
          hallii_conf.df %>% add_row(Lac = i, Comp = colnames(d8.df)[x],  S_cor = s$estimate, S_p_val=s$p.value) -> hallii_conf.df 
     }}}}}

write_delim(x = hallii_conf.df, file = "~/Documents/Kean_2021/Kean_2021_Table_S4.csv", delim = ",")

```
