---
title: "Species sorting shapes the divergence of a natural bacterial community with repeatable functionality during propagation with alternative substrates"
author: "Shepherd"
format: docx
editor: visual
---

## Install and load packages

```{r}
library(pairwiseAdonis)
library(pheatmap)
library(tidyverse)
library(readxl)
library(dplyr)
library(data.table)
library(factoextra)
library(RColorBrewer)
library(ggpubr)
library(pheatmap)
library(tidyr)
library(ggplot2)
library(viridis)
library(vegan)
library(permute)
library(lattice)
library(gridExtra)
library(ggdendro)
library(grid)
library(FactoMineR)
library(dunn.test) 
library(grid)
```

Set working directory, read files, organize data frame, plot and analyse data

Volatile organic compounds (VOC) plots before and after fermentation

```{r}
#read files

VOC_B <- read_excel("raw_substrate_GCMS_UptoC10LEfSe.xlsx", sheet = "F", col_names = T)
View(VOC_B)

VOC_A <- read_excel("sn_F038_GCMS_manuscriptLEfSe.xlsx", sheet = "M", col_names = T)
View(VOC_A)


#define character and numeric data in files
VOC_B[, 1:5] <- lapply(VOC_B[, 1:5], as.character)
VOC_A[, 1:5] <- lapply(VOC_A[, 1:5], as.character)

VOC_B[, 6:ncol(VOC_B)] <- lapply(VOC_B[, 6:ncol(VOC_B)], as.numeric)
VOC_A[, 6:ncol(VOC_A)] <- lapply(VOC_A[, 6:ncol(VOC_A)], as.numeric)

#merge files
VOC_combined <- bind_rows(VOC_B,VOC_A)

View(VOC_combined)

print(sapply(VOC_combined, class)) # to view data class: numeric and character

########## normaliseby log trnasformation recommended for conference by Prof. Smid ############################
VOC_combinedN <- VOC_combined %>%
  select(-sample_ID, -substrate_type, -cycle_no, -propagation_phase, -site) %>% mutate_all(~ log(. + 0.001)) %>% mutate_all(function(x) log(x + 1)) %>% mutate_all(function(x) {
    median_val <- median(x, na.rm = TRUE)
        x / median_val
  }) %>% mutate_all(~ ifelse(is.infinite(.) | is.nan(.) | is.na(.), 0.001, .))

print(VOC_combinedN)
summary(VOC_combinedN)
#######################################################

#option: impute dataVOC_imput <- as.data.frame(lapply(VOC_n, function(x) ifelse(is.na(x), median(x, na.rm = TRUE), x)))str(VOC_imput)

#Combine meta to normalised data

VOC_NC <- cbind(VOC_combined[, 1:5], VOC_combinedN)
View(VOC_NC)

########################## reorder starter to be first, then pheatmap with both raw and mabisi samples ##############

VOC_NC$site <- factor(VOC_NC$site, levels = c("starter", "raw", "farm_1", "farm_2", "farm_3"))
VOC_NC$substrate_type <- factor(VOC_NC$substrate_type, levels = c("starter", "RCM", "F100", "S26", "LFM", "FCM"))
VOC_NC$propagation_phase <- factor(VOC_NC$propagation_phase, levels = c("starter", "early_phase", "late_phase"))

VOC_NC <- VOC_NC %>%
  arrange(site)

dim(as.matrix(VOC_NC))
mat <- as.matrix(VOC_NC %>% select(-sample_ID, -substrate_type, -cycle_no, -propagation_phase, -site))
mat[!is.finite(mat)] <- 0.001
rownames(mat) <- VOC_NC$sample_ID
colnames(mat) <- colnames(VOC_NC %>% select(-sample_ID, -substrate_type, -cycle_no, -propagation_phase, -site))
site <- VOC_NC$site
substrate_type <- VOC_NC$substrate_type
propagation_phase <- VOC_NC$propagation_phase

annotation_data <- data.frame(site=site, substrate_type=substrate_type,propagation_phase=propagation_phase, row.names = VOC_NC$sample_ID)

my_color = list(
    site = c(starter = "#A6CEE3", raw = "#1F78B4", farm_1 = "#B2DF8A", farm_2 = "#33A02C", farm_3 = "#FB9A99"),
    substrate_type = c(starter = "#A6CEE3", RCM = "#FDBF6F", F100 = "#FF7F00", S26 = "#CAB2D6", LFM = "#6A3D9A", FCM = "#FFFF99"),
    propagation_phase = c(starter = "#A6CEE3", early_phase = "#E7298A", late_phase = "#666666")
)

heatVOC <- pheatmap::pheatmap(t(mat),
                           annotation_col =  annotation_data, 
                           annotation_colors = my_color,
                           cluster_cols = FALSE, # or TRUE if you want to cluster columns
                           fontsize = 15, 
                           angle_col = 90, 
                           fontsize_row = 12, 
                           fontsize_col = 5, 
                           width = 10,
                           height = 12,
                           main = "", 
                           fontface_row = "bold",
                           fontface_col = "bold")
heatVOC
ggsave(filename = "heatmapVOC.png", plot = heatVOC, 
       width = 20, height = 9)

###################################### filter out raw samples due to low VOC and photo oxidation ########

VOC_NC_afterFermentation <- VOC_NC %>%
  filter(site != "raw")
VOC_NC_afterFermentation

dim(as.matrix(VOC_NC_afterFermentation))
mat <- as.matrix(VOC_NC_afterFermentation[,6:ncol(VOC_NC_afterFermentation)])
mat[!is.finite(mat)] <- 0.001
rownames(mat) <- VOC_NC_afterFermentation$sample_ID
colnames(mat) <- colnames(VOC_NC_afterFermentation[,6:ncol(VOC_NC_afterFermentation)])
site <- VOC_NC_afterFermentation$site
substrate_type <- VOC_NC_afterFermentation$substrate_type
propagation_phase <- VOC_NC_afterFermentation$propagation_phase

annotation_data <- data.frame(site=site, substrate_type=substrate_type,propagation_phase=propagation_phase, row.names = VOC_NC_afterFermentation$sample_ID)

heatVOCAF <- pheatmap::pheatmap(t(mat),
                           annotation_col =  annotation_data, 
                           annotation_colors = my_color,
                           cluster_cols = FALSE, # or TRUE if you want to cluster columns
                           fontsize = 15, 
                           angle_col = 90, 
                           fontsize_row = 12, 
                           fontsize_col = 6, 
                           width = 10,
                           height = 12,
                           main = "",
                           fontface_row = "bold",
                           fontface_col = "bold")
heatVOCAF
ggsave(filename = "heatmapVOC_AF.png", plot = heatVOCAF, 
       width = 20, height = 9)

##############################

```

PCA and Multivariate analysis of VOC after fermentation only!!

```{r}

#perform PCA

rownames(VOC_NC_afterFermentation) <- VOC_NC_afterFermentation$sample_ID

pca.data <- PCA(VOC_NC_afterFermentation[,-c(1:5)])
pca.data

#Check for outliers through PCA###

PCA_VOCsuboutlier <- fviz_pca_biplot(
  pca.data,
  geom = "text",                               # keep points
  label = "ind",                                # show sample IDs
  repel = TRUE,                                 # avoid overlaps
  habillage = VOC_NC_afterFermentation$substrate_type, # color by substrate
  labelsize = 2,
  arrowsize = 0.1
)
PCA_VOCsuboutlier
ggsave("PCA_VOCsubOutlier.png", PCA_VOCsuboutlier, dpi = 300, bg = "white")

############################Exclude 9 outlier samples, due to technical error from GC-MS###################################################

VOC_no_outliers <- VOC_NC_afterFermentation[!VOC_NC_afterFermentation$sample_ID %in% c("S72", "S107", "S139", "S144", "S152", "S153", "S170", "S176", "S179"), ]

pca.data <- PCA(VOC_no_outliers[,-c(1:5)])
pca.data

#Check for outliers through PCA###

PCA_VOCNoutlier <- fviz_pca_biplot(
  pca.data,
  geom = "text",                               # keep points
  label = "ind",                                # show sample IDs
  repel = TRUE,                                 # avoid overlaps
  habillage = VOC_no_outliers$substrate_type, # color by substrate
  labelsize = 2,
  arrowsize = 0.1
)
PCA_VOCNoutlier
#################################################################################
# make color pellet a colorblind friendly

color_palette <- c(starter = "#E69F00", #orange
                   RCM = "#56B4E9", #skyblue
                   F100 = "#009E73", #green
                   S26 = "#CC79A7", #light pink
                   LFM = "#0072B2", #dark blue
                   FCM = "#D55E00")#dark orange

##############################################################################
#biplot, color by VOC representation and repel to avoid overlap

PCA_VOCsub <- fviz_pca_biplot(pca.data, 
                geom = "text",
                label = "var", # Show only variable labels 
                repel = TRUE,                   # Avoid overlapping for variables
                col.var = "#000000",            # Variables color
                habillage = VOC_no_outliers$substrate_type,  # Group by substrate type
                palette = color_palette,               # Color-blind friendly palette
                addEllipses = TRUE,             # Add ellipses for substrate types
                legend.title = "Substrate type", 
                select.var = list(contrib = 18),         
                labelsize = 2,# Adjust label size
                arrowsize = 0.1,
                title = "")                       
  
PCA_VOCsub
ggsave("PCA_VOCandsub.png", PCA_VOCsub, dpi = 300, bg = "white")

##############################################################################
# by substrate over propagation phase

PCA_MaS <- fviz_pca_biplot(pca.data, label = "none", geom.ind = "point", geom.var = "text") + 
  geom_point(aes (colour = (VOC_no_outliers$substrate_type)), size = 2) +
  stat_ellipse(aes(group = VOC_no_outliers$propagation_phase, linetype = VOC_no_outliers$propagation_phase), geom = "path", level = 0.95, alpha = 0.8, linewidth = 1.5) + 
  labs(linetype = "Propagation phase") + 
  guides(colour = guide_legend(title = "substrate_type")) +
  ggtitle("") +
  theme(plot.title = element_text(margin = margin(b = 20)),
        legend.margin = margin(t = 0, b = 0, l = 0, r = 0), 
        legend.text = element_text(size = 8), 
        legend.key.size = unit(0.5, "cm"),
        plot.margin = margin(10, 10, 10, 10))+
   scale_colour_manual(name = "substrate_type",
                      values = color_palette,
                      breaks = unique(VOC_no_outliers$substrate_type),  
                      labels = unique(VOC_no_outliers$substrate_type)) +
  scale_linetype_manual(values = c("solid","dotted", "dashed"))
PCA_MaS
ggsave("PCA_substrate and propagation.png", PCA_MaS, dpi = 300, bg = "white")

########################################################################

# by site over propagation phase
site_color_palette <- c(starter = "#E69F00",   # Orange
                        farm_1 = "#56B4E9",      # Sky Blue
                        farm_2 = "#009E73",   # Green
                        farm_3 = "#F0E442")   # Yellow

PCA_MaL <- fviz_pca_biplot(pca.data, label = "none", geom.ind = "point", geom.var = "text") + 
  geom_point(aes (colour = (VOC_no_outliers$site)), size = 2) +
  stat_ellipse(aes(group = VOC_no_outliers$propagation_phase, linetype = VOC_no_outliers$propagation_phase), geom = "path", level = 0.95, alpha = 0.8, linewidth = 1.5) + 
  labs(linetype = "Propagation phase") + 
  guides(colour = guide_legend(title = "site")) +
  ggtitle("") +
  theme(plot.title = element_text(margin = margin(b = 20)),
        legend.margin = margin(t = 0, b = 0, l = 0, r = 0), 
        legend.text = element_text(size = 8), 
        legend.key.size = unit(0.5, "cm"),
        plot.margin = margin(10, 10, 10, 10))+
   scale_colour_manual(name = "site",
                      values = site_color_palette,
                      breaks = unique(VOC_no_outliers$site),  
                      labels = unique(VOC_no_outliers$site)) +
  scale_linetype_manual(values = c("solid", "dotted","dashed"))
PCA_MaL
ggsave("PCA_site.png", PCA_MaL, dpi = 300, bg = "white")

###########################################################################

# MULTIVARIATE ANALYSIS for VOC

# Compute bray distance
VOC_after_dist <- vegdist(as.matrix(VOC_no_outliers[, -c(1:5)]), method = "bray")

# by substrate type
VOC_SubstrateType <- adonis2(VOC_after_dist ~ as.factor(VOC_no_outliers$substrate_type), data = VOC_no_outliers, permutations = 999) 
VOC_SubstrateType

VOC_SubstrateTypePW <- pairwise.adonis(VOC_after_dist,VOC_no_outliers$substrate_type) 
corrected_p_values <- p.adjust(VOC_SubstrateType$p.value, method = "BH")
VOC_SubstrateTypePW

#by propagation phase
VOC_after.propagation <- adonis2(VOC_after_dist ~ as.factor(VOC_no_outliers$propagation_phase), data = VOC_no_outliers, permutations = 999) #no pairwise, P > 0.05
VOC_after.propagation

VOC_after.propagationW <- pairwise.adonis(VOC_after_dist,VOC_no_outliers$propagation_phase) 
corrected_p_values <- p.adjust(VOC_after.propagation$p.value, method = "BH")
VOC_after.propagationW

##by site

VOC_after.site <- adonis2(VOC_after_dist ~ as.factor(VOC_no_outliers$site), data = VOC_no_outliers, permutations = 999) #no pairwise, P > 0.05
VOC_after.site
#################


```

pH

```{r}


#read files, assign character and numeric to variables

pHCons <- read_excel("Shepherd_F038_pH&Consistency_ManuscriptLEfSe.xlsx", sheet = "Sheet1")
View(pHCons) 

pHCons[, 1:5] <- lapply(pHCons[, 1:5], as.character)
pHCons[, 6:9] <- lapply(pHCons[, 6:9], as.numeric)

str(pHCons)

#### boxplot of pH before fermentation (T0)
pHCons_filtered <- pHCons %>% filter(!is.na(pH_T0) & site != "starter") #drop starter, has no pH values at T0 and for consistency.

############################# plot pH T0, with RCM ref#################

RCM_pHT0 <- median(pHCons$pH_T0[pHCons$substrate_type == "RCM"], na.rm = TRUE) # calculatiing median for RCM as reference

box_phT0 <- ggplot(data = pHCons_filtered, aes(x = propagation_phase, y = pH_T0, fill = substrate_type)) +
  geom_violin(alpha = 0.3) +  
  geom_hline(yintercept = RCM_pHT0, linetype = "dashed", color = "black", size = 0.8) +
  theme_classic() +  
  scale_x_discrete(limits = c("early_phase", "late_phase")) + 
  scale_fill_brewer(palette = "Set3") +  
  labs(title = "", x = "Propagation phase", y = "pH") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = -90),  # Rotate x-axis labels
    strip.background = element_blank()
  )
box_phT0
ggsave("boxPlot_pHT0.png", box_phT0, dpi = 300)

#normality

qqnorm(pHCons$Mabisi_pH_T48)
qqline(pHCons$Mabisi_pH_T48, col = "red") #qq plot reveals not normaly distributed, data points away from reference line

shapiro_test_result <- shapiro.test(pHCons$Mabisi_pH_T48)
print(shapiro_test_result) # p < 0.05, reject null hypothesis that data is normaly distributed, aggrement with qq plot.

### box plot of pH at T48: substrate over propagation phases##
RCM_pH <- median(pHCons$Mabisi_pH_T48[pHCons$substrate_type == "RCM"], na.rm = TRUE) # calculatiing median for starter as reference

pHCons_filtered$propagation_phase <- factor(pHCons_filtered$propagation_phase, levels = c("early_phase", "late_phase"))

box_phT48 <- ggplot(data = pHCons_filtered, aes(x = propagation_phase, y = Mabisi_pH_T48, fill = substrate_type)) + 
  geom_violin(alpha = 0.3) +  
  geom_hline(yintercept = RCM_pH, linetype = "dashed", color = "black", size = 0.8) +  
  theme_classic() +  
  scale_x_discrete(limits = c("early_phase", "late_phase")) + 
  scale_fill_brewer(palette = "Set3") +  
  labs(title = "", x = "Propagation phase", y = "pH") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = -90),  # Rotate x-axis labels
    strip.background = element_blank()
  )
box_phT48
ggsave("boxPlot_pHT48.png", box_phT48, dpi = 300)

############################ pH at different farm sites T48##################

box_phT48FS <- ggplot(data = pHCons_filtered, aes(x = site, y = Mabisi_pH_T48, fill = substrate_type)) +
  geom_violin(alpha = 0.3) +  
  geom_hline(yintercept = RCM_pH, linetype = "dashed", color = "black", size = 0.8) +  
  theme_classic() +  
  scale_x_discrete(limits = c("Farm_1", "Farm_2", "Farm_3")) + 
  scale_fill_brewer(palette = "Set3") +  
  labs(title = "", x = "Farm site", y = "pH") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = -90),  # Rotate x-axis labels
    strip.background = element_blank()
  )
box_phT48FS
ggsave("boxPlot_FarmSitepHT48.png", box_phT48FS, dpi = 300)

######### non parametric Multivariate analysis ###############################



#T48 and propagation phase
Wcox_pHT48_s <- wilcox.test(Mabisi_pH_T48 ~ propagation_phase, data = pHCons_filtered)
Wcox_pHT48_s


#T48 and substrate type
Wcox_pHSub <- kruskal.test(Mabisi_pH_T48 ~ substrate_type, data = pHCons_filtered)
Wcox_pHSub
Wcox_pHSub_dunntest <- dunn.test(pHCons_filtered$Mabisi_pH_T48, g = pHCons_filtered$substrate_type, method = "BH")
Wcox_pHSub_dunntest

# location

Wcox_pHSit <- kruskal.test(Mabisi_pH_T48 ~ site, data = pHCons_filtered)
Wcox_pHSit
Wcox_pHSit_dunntest <- dunn.test(pHCons_filtered$Mabisi_pH_T48, g = pHCons_filtered$site, method = "BH")
Wcox_pHSit_dunntest

```

Consistency

```{r}
#### Consistency after fermentation: substrate over time
view(pHCons_filtered)

############################seting reference consistency from RCM substrate###

### box plot of pH at T48: substrate over propagation phases##
RCM_cons <- median(pHCons_filtered$consistency_cm[pHCons_filtered$substrate_type == "RCM"], na.rm = TRUE) # calculatiing median for starter as reference

############## box plot of consistency cm by propagation phase###############

pHCons_filtered$propagation_phase <- factor(pHCons_filtered$propagation_phase, levels = c("early_phase", "late_phase"))

box_Consist <- ggplot(data = pHCons_filtered, aes(x = propagation_phase, y = consistency_cm, fill = substrate_type)) + 
  geom_violin(alpha = 0.3) +  
  geom_hline(yintercept = RCM_cons, linetype = "dashed", color = "black", size = 0.8) +  
  theme_classic() +  
  scale_x_discrete(limits = c("early_phase", "late_phase")) + 
  scale_fill_brewer(palette = "Set3") +  
  labs(title = "", x = "Propagation phase", y = "Consistency (cm)") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = -90),  # Rotate x-axis labels
    strip.background = element_blank()
  )
box_Consist 
ggsave("ConsPlot_phases.png", box_Consist, dpi = 300) 

   
#####################################################################

#### Consistency after fermentation by substrate and farm site####

box_consistencySite <- ggplot(data = pHCons_filtered, aes(x = site, y = consistency_cm, fill = substrate_type)) + 
  geom_violin(alpha = 0.3) +  # Violin plot to show distribution
  geom_hline(yintercept = RCM_cons, linetype = "dashed", color = "black", size = 0.8) +
 theme_classic()  +  
  scale_x_discrete(limits = c("Farm_1", "Farm_2", "Farm_3")) + 
  scale_fill_brewer(palette = "Set3") + 
  labs(title = "", x = "Farm site", y = "Consistency (cm)") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = -90),  # Rotate x-axis labels
    strip.background = element_blank()
  )
box_consistencySite
ggsave("ConsPlot_site.png", box_consistencySite, dpi = 300)
#############################################################################
### PERMANOVA of consistency outcome#####
KW_Cons_Substrate <- kruskal.test(consistency_cm ~ substrate_type, data = pHCons_filtered)
KW_Cons_Substrate
KW_Cons_SubstratePW <- dunn.test(pHCons_filtered$consistency_cm, g = pHCons_filtered$substrate_type, method = "BH")
KW_Cons_SubstratePW

KW_ConsProp <- wilcox.test(consistency_cm ~ propagation_phase, data = pHCons_filtered)
KW_ConsProp # no differencein consistency between propagation phase, no need for pairwise analysis.

KW_ConsSite <- kruskal.test(consistency_cm ~ site, data = pHCons_filtered)
KW_ConsSite #site did not differ in consistency, no need for pairwise analysis.
#######################################################################

```

Combine consistency and pH T48 plots

```{r}
pHCons_plots <- (box_phT48 / box_CT48)+
  plot_annotation(tag_levels = 'a')
pHCons_plots

ggsave("pH&Consplots combined.png", plot = pHCons_plots, dpi = 300, width = 7, height = 8)

################################### combined pH before propagation, after propagation by farm site and consistency by farm site############

pH_cons_farm <- (box_phT0/box_phT48FS/box_consistencySite) +
  plot_annotation(tag_levels = 'a')
pH_cons_farm
ggsave("pH_Cons_farmplots combined.png", plot = pH_cons_farm, dpi = 300, width = 7, height = 8)
```
