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

## Install and load related packages

```{r}
library(remotes)
library(microbiomeMarker)
library(htmltools)
library(phyloseq)
library(cli)
library(treeio)
library(microbiome)
library(tidyverse)
library(palmerpenguins)
library(here)
library(knitr)
library(kableExtra)
library(patchwork)
library(broom)
library(readxl)
library(vegan)
library(devtools)
library(microViz)
library(Biostrings)
library(ggplot2)
library(dplyr)
library(lubridate)
library(circlize)
library(dendextend)
library(RColorBrewer)
library(viridis)
library(dada2)
library(BiMiCo)
library(MicEco)
library(forcats)
library(gridExtra)
library(flextable)
library(picante)
library(ggpubr)
library(rstatix)
library("DESeq2")
library(ggrepel)
library(gridExtra)
library(cowplot)
library(ggh4x)
library(pheatmap)
library(pairwiseAdonis)
library(writexl)
```

Setting working directory to location of fasta file and reference silva_nr99_v138.2_toGenus_trainset.fa.gz run for revised taxonomy of lactic acid bacteria.

```{r}

# Load sequences into a DNAStringSet
asv_seqs <- readDNAStringSet("ASV-dna-sequences.fasta")

set.seed(100)
taxa <- assignTaxonomy(asv_seqs, "silva_nr99_v138.2_toGenus_trainset.fa.gz", multithread=TRUE)
unname(taxa)

asv <- names(asv_seqs)  # these are the headers in the FASTA
asv

# Convert to data.frame and save
tax_df <- as.data.frame(taxa)
tax_df

# rownames(taxa) are DNA sequences, so match them to FASTA sequences
tax_df$asv <- asv[match(rownames(tax_df), as.character(asv_seqs))]
tax_df


#Drop sequences completely by resetting rownames
rownames(tax_df) <- NULL

##Optional: reorder columns so ASV_ID is first
tax_df <- tax_df[, c("asv", colnames(taxa))]
tax_df

write_xlsx(tax_df, "tax_df.xlsx")

##########end##########
```

Import ASV table and metadata Excel files, merge with tax_df for downstream analysis

```{r}

otu_mat<- read_excel("sn_F038_seqIDLEfSeNT.xlsx", sheet = "asv")
tax_mat<- read_excel("tax_df.xlsx", sheet = 1) 
samples_df <- read_excel("sn_F038_seqIDLEfSeNT.xlsx", sheet = "met")
summary(samples_df)
```

#change variables from character to factor

```{r}

cols <- c("sample_ID", "substrate_type", "cycle_no", "propagation_phase", "site", "replicate")

samples_df %<>%
       mutate_each_(funs(factor(.)),cols)

str(samples_df)


```

#Re-organise the data into a matrix

```{r}
otu_mat <- otu_mat %>%
  tibble::column_to_rownames("asv")
tax_mat <- tax_mat %>% 
  tibble::column_to_rownames("asv")

samples_df <- samples_df %>% 
  tibble::column_to_rownames("sample_ID")

otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)

class(otu_mat) <- "numeric" #change OTU table to numeric, to correct OTU non-numeric object output.
OTU <- otu_table(otu_mat, taxa_are_rows = TRUE)
TAX <- tax_table(tax_mat)
samples <- sample_data(samples_df)
```

#Create a phyloseq object with tree and refseq

```{r}
rs <- readDNAStringSet(file = "ASV-dna-sequences.fasta",
      format = "fasta", nrec = -1L, skip = 0L, seek.first.rec = FALSE, use.names = TRUE)

tree <- read.tree(file = "tree.nwk")
sample_names()   
psrt <- phyloseq(OTU, TAX, samples, rs, tree)
psrt
```

#plot tree of top 40 taxa

```{r}
tree_40 <- prune_taxa(taxa_names(psrt)[1:40], psrt)

colnames(meta(psrt))

plot_tree(tree_40, nodelabf=nodeplotboot(), ladderize="right", 
          color="cycle_no", label.tips="Genus") #check potential of unwanted hits: Archaea, chloroplasts and mitochondria.
```

#culculate the total number of reads for each sample and sort them to ID problematic samples.

```{r}
RN <- sample_sums(psrt)
sort(RN)

rank_names(psrt)
head(tax_table(psrt)) 

phy <- table(tax_table(psrt)[, 2])
phy

genus <- table(tax_table(psrt)[, 6])
genus
```

#filter out unwanted taxa: Archaea, chloroplasts, and mitochondria & to minimise clutter on plots, take OTU that represent atleast 0.25% of reads in atleast one sample

```{r}
#option 1: use subset taxa functionto remove non bacteria sequences 

psrt_filt <- subset_taxa(psrt, !Kingdom %in% c("Archaea","Eukaryota") | is.na(Kingdom)) #removing archaea,eukaryota, but keeping NA
psrt_filt <- subset_taxa(psrt_filt, !Order %in% c("Chloroplast") | is.na(Order)) #removing Chloroplast, but keeping NA
psrt_filt <- subset_taxa(psrt_filt, !Family %in% c("Mitochondria") | is.na(Family)) #removing Mitochondria, but keeping NA
psrt_filt

tree_2 <- prune_taxa(taxa_names(psrt_filt)[1:40], psrt_filt)

plot_tree(tree_2, nodelabf=nodeplotboot(), ladderize="right", 
          color="cycle_no", label.tips="Genus") #check potential of unwanted hits: Archaea, chloroplasts and mitochondria.

#Option 2: using rm nonbac function via BiMiCo package to removing non bacteria sequences, or spurious taxa: 

psrtB <-rm_nonbac(psrt)
psrtB
 
tree_3 <- prune_taxa(taxa_names(psrtB)[1:40], psrtB)

plot_tree(tree_3, nodelabf=nodeplotboot(), ladderize="right", 
          color="cycle_no", label.tips="Genus")

#option one data cleaning and tree ploting gives cleaner tree;


```

\# Normalisation:rarefaction & rarefaction curve with minimum library size, and data trimming by 0.25% minimum relative abundance

```{r}

#part A - rafraction:

tab <- otu_table(psrt_filt) #to build/access the OTU_table
tab

class(tab) <- "matrix" #set data frame as matrix

tab <- t(tab) #transpose observation
tab

raremin <- min(rowSums(tab)) #to ID minimum sample count
raremin

rare_psrtFR <- rarefy_even_depth(psrt_filt, rngseed=711, sample.size=19534, replace=F) #resamples an OTU table such that all samples have same library size - rarefaction/normalisation

rare_psrtFRc <- rarecurve(tab, step=10000, lwd=2, ylab="OTU",  label=F) #plots rarefaction curve

rare_psrtFRcm <- rarecurve(tab, step=50, lwd=1, ylab="Number of ASVs",  main = "Rarefaction curves", label=F, col= "blue", cex = 0.4)
 abline(v = raremin) #plots rarefaction curve with min library size position
 
#part B - Normalise by median sequencing depth and relative abundance (trim off low counts) of  0.25%

total = median(sample_sums(psrt_filt))
normalised0.25 <- filter_taxa(psrt_filt, function(x) sum(x > total*0.025) > 0, TRUE)
normalised0.25

#We select rarefaction (part A) as recommended in microbiome data analysis (Cameron et al., 2021)

```

\# Rarefy and prune 0.25% to exclude spurious taxa

```{r}

rare_psrtFR

sample_sums(psrt_filt) #checking the sequencing depth before rafaction

sample_sums(rare_psrtFR) #checking the sequencing depth after rafaction

#Combination of nonbacterial removal in option 1, rarefaction and & applying ps_prune 0.25% relative abundance cut off to prune spurious taxa (MicEco package).  

rarefyphy0.25 <- ps_prune(rare_psrtFR, min.samples = 0, min.reads = 4, min.abundance = 0.0025)

rarefyphy0.25

phy_df <- psmelt(rarefyphy0.25) #convert whole phyloseq to data frame

phy_df

```

Plot graphs for visualisation: by propagation cyphase, substrate type and site

```{r}
# Create color group #########
num_colors <- 11
palette_colors <- colorRampPalette(brewer.pal(8, "Dark2"))(num_colors)
palette.colors()
barplot(rep(1, num_colors), col = palette_colors, space = 0, axes = FALSE)
print(palette_colors)
############## histoghraph with early and late propagation phase###########


#reorder factor levels for propagation phase

phy_df <- phy_df %>%
  mutate(propagation_phase = factor(propagation_phase, levels = c('O', 'early_phase', 'late_phase')))
phy_df

# Create the plot
Phy_histp_NT <- phy_df %>%
  group_by(Sample) %>%
  mutate(percent_abundance = Abundance/sum(Abundance) * 100) %>%
  ggplot(aes(x = Sample, y = percent_abundance, fill = Phylum)) +
  theme_bw() + 
  geom_bar(stat = "identity") +
  facet_nested(~ fct_relevel(propagation_phase, "AC", "short_term", "long_term") + 
      fct_relevel(substrate_type, "RCM", "F100", "S26", "LFM", "FCM"), 
    scales = "free", space = "free") +
  theme(axis.text.x = element_text(angle = 90, size = 4, hjust = 1, vjust = 0.5,face = "bold"), 
        strip.text = element_text(size = 5, face = "bold")) + 
  scale_fill_manual(values = palette_colors, 
    labels = c("Bacillota","Pseudomonadota","not assigned")) + 
    labs(x = "sample ID", y = "relative abundance")
Phy_histp_NT

ggsave("hist_phylum_genus.png", plot = Phy_histp_NT, dpi = 300, width = 12, height = 4)

################################# propagation by substrate by genus##############


prop_histNT2 <- phy_df %>%
  group_by(Sample) %>%
  mutate(percent_abundance = Abundance/sum(Abundance) * 100) %>%
  ggplot(aes(x = Sample, y = percent_abundance, fill = Genus)) +
  theme_bw() +
  geom_bar(stat = "identity") +
  facet_nested(~ fct_relevel(propagation_phase, "AC", "short_term", "long_term") + 
      fct_relevel(substrate_type, "RCM", "F100", "S26", "LFM", "FCM"), 
    scales = "free", space = "free") +
  theme(axis.text.x = element_text(angle=90, size=4, hjust = 1, vjust = 0.5,face = "bold"), 
        strip.text = element_text(size = 5, face = "bold")) + 
  guides(fill=guide_legend(title="Genus")) + 
  scale_fill_manual(values = palette_colors, labels = c("Acetobacter", "Clostridium", "Enterobacter","Enterococcus","Klebsiella","Lacticaseibacillus", "Lactiplantibacillus","Lactococcus","Leuconostoc","Paucilactobacillus", "not assigned")) +
  labs(x = "sample ID", y = "relative abundance")
prop_histNT2

ggsave("hist_SubstPropagation_genus.png", plot = prop_histNT2, dpi = 300, width = 12, height = 4)

##########################propagation by farm graph ###########

subs_PropFarmNT <- phy_df %>%
  group_by(Sample) %>%
  mutate(percent_abundance = Abundance/sum(Abundance) * 100) %>%
  ggplot(aes(x = Sample, percent_abundance, fill = Genus)) +
  theme_bw() + 
  geom_bar(stat = "identity") +
  facet_nested(~ fct_relevel(propagation_phase, "AC", "short_term", "long_term") + 
      fct_relevel(site, "farm_1", "farm_2", "farm_3"), 
    scales = "free", space = "free")  +
  theme(axis.text.x = element_text(angle=90, size=4, face = "bold")) + guides(fill=guide_legend(title="Genus")) + 
  scale_fill_manual(values = palette_colors, labels = c("Acetobacter", "Clostridium", "Enterobacter","Enterococcus","Klebsiella","Lacticaseibacillus", "Lactiplantibacillus","Lactococcus","Leuconostoc","Paucilactobacillus", "not assigned")) +
  labs(x = "sample ID", y = "relative abundance")
subs_PropFarmNT

ggsave("hist_FarmSubstrate_genus.png", plot = subs_PropFarmNT, dpi = 300, width = 12, height = 4)

```

#Alpha diversity - box plots & alpha diversity indices

```{r}

#perform alpha diversity Chao1 richness and Shannon diversity estimators (top 40 ASV), refer to Yan Hui tutorial for microbiome analysis in R.
manual_colors <- brewer.pal(6, "Set2")
print(manual_colors)
###################### substrate variation ########################
alpha_ps1 <- plot_richness(rarefyphy0.25, x = "propagation_phase", measures = c("Chao1", "Shannon")) +
  geom_boxplot(aes(fill = substrate_type), alpha = 0.7) +  # Remove color to prevent duplicate legend
  theme_classic() + 
  scale_x_discrete(limits = c("O", "early_phase", "late_phase")) +
  scale_fill_manual(values = manual_colors) +  
  theme(strip.background = element_blank(), 
        axis.text.x = element_text(angle = -90))

alpha_ps1
ggsave("boxplot_alpha SubstrateF038.png", plot = alpha_ps1, dpi = 300, width = 12, height = 4)

################## farm location ########################################

alpha_ps2 <- plot_richness(rarefyphy0.25, x = "propagation_phase", measures = c("Chao1", "Shannon")) +
  geom_boxplot(aes(fill = site)) + 
  theme_classic() + 
  scale_x_discrete(limits = c("O", "early_phase", "late_phase"))+
  scale_fill_manual(values = manual_colors) +  # Set the order of the x-axis
  theme(strip.background = element_blank(), axis.text.x = element_text(angle = -90))
alpha_ps2
ggsave("boxplot_alpha FarmSiteF038.png", plot = alpha_ps2, dpi = 300, width = 12, height = 4)

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

#perform alpha diversity estimates, Wilcoxon rank-sum test: 

rich = estimate_richness(rarefyphy0.25)
rich

# alpha diversity pairwise Wilcox test - propagation cycle: (i) propagation phase Chao1

P_chao_psrt <- pairwise.wilcox.test(rich$Chao1, sample_data(rarefyphy0.25)$propagation_phase, p.adjust.method = "BH")
P_chao_psrt
P_tab_Cha <- P_chao_psrt$p.value %>% as.data.frame() %>% tibble::rownames_to_column(var = "group1") %>%
  gather(key="group2", value="p.adj", -group1) %>%
  na.omit()
P_tab_Cha
P_tab_Cha %>% flextable() %>% autofit()

# alpha diversity pairwise Wilcox test - propagation cycle: (ii) propagation phase Shannon

P_shan_psrt <- pairwise.wilcox.test(rich$Shannon, sample_data(rarefyphy0.25)$propagation_phase, p.adjust.method = "BH")
P_shan_psrt
P_tab_Sha <- P_shan_psrt$p.value %>% as.data.frame() %>% tibble::rownames_to_column(var = "group1") %>%
  gather(key="group2", value="p.adj", -group1) %>%
  na.omit()
P_tab_Sha %>% flextable() %>% autofit()

# alpha diversity pairwise Wilcox test - substrate variation: (iii) substrate Chao1
#Chao1
s_chao1_psrt <- pairwise.wilcox.test(rich$Chao1, sample_data(rarefyphy0.25)$substrate_type, p.adjust.method = "BH")
s_chao1_psrt
s_tab_chao1 <- s_chao1_psrt$p.value %>% as.data.frame() %>% tibble::rownames_to_column(var = "group1") %>%
  gather(key="group2", value="p.adj", -group1) %>%
  na.omit()
s_tab_chao1 %>% flextable() %>% autofit()
##substrate shannon
s_chao1_psrt <- pairwise.wilcox.test(rich$Shannon, sample_data(rarefyphy0.25)$substrate_type, p.adjust.method = "BH")
s_chao1_psrt
s_tab_chao1 <- s_chao1_psrt$p.value %>% as.data.frame() %>% tibble::rownames_to_column(var = "group1") %>%
  gather(key="group2", value="p.adj", -group1) %>%
  na.omit()
s_tab_chao1 %>% flextable() %>% autofit()

# diversity index: Chao1 & Shannon based on farm site

f_chao1_psrt <- pairwise.wilcox.test(rich$Chao1, sample_data(rarefyphy0.25)$site, p.adjust.method = "BH")
f_chao1_psrt
s_tab_chao1F <- f_chao1_psrt$p.value %>% as.data.frame() %>% tibble::rownames_to_column(var = "group1") %>%
  gather(key="group2", value="p.adj", -group1) %>%
  na.omit()
s_tab_chao1F %>% flextable() %>% autofit()

f_shan_psrt <- pairwise.wilcox.test(rich$Shannon, sample_data(rarefyphy0.25)$site, p.adjust.method = "BH")
f_shan_psrt
s_tab_shanF <- f_shan_psrt$p.value %>% as.data.frame() %>% tibble::rownames_to_column(var = "group1") %>%
  gather(key="group2", value="p.adj", -group1) %>%
  na.omit()
s_tab_shanF %>% flextable() %>% autofit()

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

```

\# Beta diversity - plotting NMDS, PCoA and Beta diversity PERMANOVA tests

```{r}
#performing NMDS graphs by ordination & ANOSIM (NMDS is aligned with ANOSIM)
manual_colorsNMDS <- brewer.pal(8, "Dark2")
print(manual_colorsNMDS)
# Visualize the colors in manual_colorsNMDS

dist = phyloseq::distance(rarefyphy0.25, method="bray")

psrtNorm.o <- ordinate(rarefyphy0.25, "NMDS", "bray")


#################substrate variation, propagation phase and farm site combined###

D_NMDS_withfarmsite <- plot_ordination(rarefyphy0.25, psrtNorm.o, type = "sample", shape = "propagation_phase", color = "substrate_type", title = "") +
  geom_point(size = 3) +
  stat_ellipse(aes(group = site, linetype = site), level = 0.95, geom = "polygon", fill = NA, color = "black") + 
  scale_color_manual(values = manual_colorsNMDS) +
  scale_linetype_manual(name = "Farm site", values = c("solid", "dashed", "dotted", "dotdash")) +
  theme(legend.key.size = unit(0.25, "cm"),  
        legend.text = element_text(size = 10),  
        legend.title = element_text(size = 12),
        strip.background = element_blank()) +
  theme_classic()
D_NMDS_withfarmsite
ggsave("NMDS_withfamsites.png", D_NMDS_withfarmsite, width = 8, height = 6, dpi = 300)
####################################

metadata <- data.frame(sample_data(rarefyphy0.25))
metadata
anosim(dist, metadata$substrate_type)
anosim(dist, metadata$propagation_phase)
anosim(dist, metadata$site) 

####interaction: pairwise####: outcome similar

anosim_prop <- anosim(dist, interaction(metadata$substrate_type, metadata$propagation_phase))
print(anosim_prop)

anosim_site_prop <- anosim(dist, interaction(metadata$substrate_type, metadata$site))
print(anosim_site_prop)

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

```

LEfSE - indicator species per treatment variable

```{r}
# indicator species per treatment variable with LEfSe
# by substrate
lefse_result_ALL <- run_lefse(rare_psrtFR,
                              group = "substrate_type",
                              taxa = "Genus",
                              kw_cutoff = 0.05, lda_cutoff = 3)

lefse_result_ALL
marker_table(lefse_result_ALL)

# LDA bar plot
LDAsubp <- plot_ef_bar(
  lefse_result_ALL
)+ 
  scale_fill_manual(values = palette_colors)+
  theme(axis.text.y = element_text(size = 6, face = "bold"))
LDAsubp

##########################################
# by propagation phase
lefse_result_ALLp <- run_lefse(rare_psrtFR,
                              group = "propagation_phase",
                              taxa = "Genus",
                              kw_cutoff = 0.05, lda_cutoff = 3)
lefse_result_ALLp
marker_table(lefse_result_ALLp)

# LDA bar plot
LDAsubpr <- plot_ef_bar(
  lefse_result_ALLp
) + 
  scale_fill_manual(values = palette_colors)+
  theme(axis.text.y = element_text(size = 6, face = "bold"))
LDAsubpr

##########################################
# by site
lefse_result_ALLs <- run_lefse(rare_psrtFR,
                              group = "site",
                              taxa = "Genus",
                              kw_cutoff = 0.05, lda_cutoff = 3)
lefse_result_ALLs
marker_table(lefse_result_ALLs)

# LDA bar plot
LDAsubs <- plot_ef_bar(
  lefse_result_ALLs
) + 
  scale_fill_manual(values = palette_colors)+
  theme(axis.text.y = element_text(size = 8, face = "bold"))
LDAsubs
ggsave("LDAplotsite combined.png", plot = LDAsubs, dpi = 300, width = 7, height = 6)
```

Combine plots for differential abundance analysis by substrate type

```{r}
combined_plots <- (LDAsubp / LDAsubpr)+
  plot_annotation(tag_levels = 'a')
combined_plots

ggsave("LDAplots combined.png", plot = combined_plots, dpi = 300, width = 7, height = 6)

```
