---
title: "Soil properties in agricultural systems affect microbial genomic traits"
output: html_document
date: "2024-10-04"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Load packages
```{r}
# dependencies
library(ggtree)
library(Biostrings)
library(phylofactor)
library(here)
library(tidyverse)
library(ggplot2)
library(ggpmisc)
library(devtools)
library(dplyr)
library(tidyr)
library(progress)
library(mgcv)
#for microeco dependencies
# First select the packages of interest
packages <- c("reshape2", "GUniFrac", "MASS", "ggpubr", "randomForest", "ggdendro", "ggrepel", 
"agricolae", "gridExtra", "picante", "pheatmap", "igraph", "rgexf", "RJSONIO", "ggalluvial", 
"tidyverse", "magrittr", "microeco", "ggtree")
# Now check or install
lapply(packages, function(x) {
  if(!require(x, character.only = TRUE)) {
    install.packages(x, dependencies = TRUE)
  }})
```

set common theme for plotting
```{r}
#set theme
theme_set(theme_bw() +
            theme(strip.background = element_rect(colour = "black", fill = "white"), 
                  strip.text = element_text(size=14)) +
            theme(axis.text.x = element_text(angle = 90, size = 10, hjust = 0.95),
                  axis.text.y = element_text(size = 14),
                  axis.title.y = element_text(size = 12),
                  axis.title.x=element_text(size = 12, face = "bold"),
                  axis.ticks.x=element_blank()) +
            theme(plot.title = element_text(hjust = 0.5)) +
            theme(legend.text = element_text(size = 10),
                  legend.title = element_text(size = 10, face = "bold"), 
                  legend.direction = "vertical", legend.box = "vertical"))
```


Set working directory and import files
```{r}
# working directory
here::i_am("Soil properties in agricultural systems affect microbial genomic traits.RMD")
here()

#read in ASV table
ASV.tab <- read.csv(here("Analysis211014", "seqtab_16s.csv"), header = TRUE, row.names = 1, check.names = F)
nrow(ASV.tab) #check row numbers
ncol(ASV.tab) #check column numbers

#read in taxonomy
ASV.taxonomy <- read.csv(here("Analysis211014", "new.new.tax.csv"), header = TRUE, row.names = 1, check.names = F)  

#metadata, set the row names to the col containing the sequencer ID
meta1  <- read.csv(here("Analysis211014", "AgriiRhizaEnv_v5_tidy.csv"), header = TRUE, row.names = 6, check.names = F)
colnames(meta1)
sample=as.data.frame(meta1)
#make sure all are data.frame format
class(ASV.tab)
class(ASV.taxonomy)
class(sample)
```

Filter out samples with low read numbers
```{r}
#filter out samples here
reads<-colSums(ASV.tab)#calculates number of reads per sample
hist(reads,breaks=50)
#filter out <8000 
ASV.tab.sub<-ASV.tab[,reads>8000]
# #rename
ASV.tab=ASV.tab.sub
```

Make microeco object, filter unwanted reads and calculate abundances at all taxa levels 
```{r}
#create the R6 object
dataset <- microtable$new(sample_table = sample, otu_table = ASV.tab, tax_table = ASV.taxonomy)
#make the species and sample information consistent 
dataset$tidy_dataset()
print(dataset)
#remove ASVs which are not assigned in the Kingdom "k__Archaea" or "k__Bacteria".
dataset$tax_table %<>% base::subset(domain == "Archaea" | domain == "Bacteria")
print(dataset)
#remove ASVs with the taxonomic assignments "mitochondria" or "chloroplast".
dataset$filter_pollution(taxa = c("mitochondria", "chloroplast"))
print(dataset)
#make the ASVs same in otu_table, tax_table and phylo_tree, we use tidy_dataset() again.
dataset$tidy_dataset()
print(dataset)
#check the sequence numbers in each sample.
dataset$sample_sums() %>% range
# rarefy to below minread number
set.seed(7)
dataset$rarefy_samples(sample.size = 7555)
dataset$sample_sums() %>% range
#save the rarefied table
#write.csv(dataset$otu_table, file = "microeco_rarefied_seqtab_updated_silva.csv")
#save the taxa for the rarefied reads
#write.csv(dataset$tax_table, file="microeco_rarefied_taxa_updated_silva.csv")
#calculate the taxa abundance at each taxonomic rank using cal_abund()
dataset$cal_abund()
# return dataset$taxa_abund
class(dataset$taxa_abund)


```

Plot pH v LOI
```{r}
#plot pH v LOI
meta2=dataset$sample_table

meta2 %>% 
  ggplot(aes(x=pH, y=Organic_Matter_LOI, size = Organic_Matter_LOI, colour= pH))+
  geom_point(alpha=0.5) +
  scale_colour_gradient(low = "red", high="blue")+
  scale_size(range = c(.01, 10), name="LOI")+
  ylab("Organic_Matter_LOI")  +
  xlab("pH")+
  stat_poly_eq(use_label(c("eq", "adj.R2", "f", "p", "n")))


```

Calculate alpha and beta diversity
```{r}
####

# 
#calculate the alpha diversity and return dataset$alpha_diversity
dataset$cal_alphadiv(PD = FALSE)
class(dataset$alpha_diversity)

#calculate distance matrix of beta diversity and return dataset$beta_diversity
dataset$cal_betadiv(unifrac = FALSE)
class(dataset$beta_diversity)



```

Make barplots at phyla level per sample, include ten most abundant taxa and facet by sample's pH classification
```{r}
# #make barplots
# #trans_abund class
t1 <- trans_abund$new(dataset = dataset, taxrank = "phylum", ntaxa =10)
#with pH bin facet
t1$plot_bar(others_color = "grey70", xtext_keep = FALSE, legend_text_italic = FALSE, facet = "pH2")

```


Do correlation heatmap between env variables and phyla abundance
```{r}
#do correlation analysis
#create trans_diff object
colnames(sample)
t1 <- trans_env$new(dataset = dataset, add_data = sample[, (13:20)])
t1$cal_cor(use_data = "phylum", p_adjust_method = "fdr", p_adjust_type = "Env")
#keep only phyla with sign. response
t1$plot_cor(filter_feature = c("", "*", "**"))

```


Plot simpsons diversity scores for each sample against pH
```{r}
#make df of simpsons div
simps.df=vegan::diversity(t(dataset$otu_table), "simpson") %>% as.data.frame() %>%   dplyr::rename("simps.div"=".") %>% rownames_to_column(.,var="rowname")
#add simpsons to meta
meta=left_join(meta4, simps.df, by="rowname") %>% dplyr::relocate("simps.div",.after = "rowname")
#plot simps div figures with LOI and pH and stats
meta %>% 
  ggplot(aes(x=pH, y=simps.div, size = Organic_Matter_LOI, colour= alpha.div))+
  geom_point(alpha=0.5) +
  scale_colour_gradient(low = "yellow", high="darkgreen")+
  scale_size(range = c(.01, 10), name="LOI")+
  ylab("simps.div")  +
  xlab("pH")+
  stat_poly_eq(use_label(c("eq", "adj.R2", "f", "p", "n")))
```

Plot NMDS of beta-diversity
```{r}
#NMDS ordination 
#create an object and select NMDS for ordination,
t1 <- trans_beta$new(dataset = dataset, measure = "bray")
t1$cal_ordination(method = "NMDS")
# t1$res_ordination is the ordination result list
class(t1$res_ordination)
# plot the  result
nmds=t1$plot_ordination(point_size = 0.1)
nmds+
  geom_point(aes(color = dataset$sample_table$pH, size = dataset$sample_table$Organic_Matter_LOI, alpha=0.8, shape = factor(dataset$sample_table$Crop_rotation)))+
  scale_colour_gradient(low = "red", high="blue", name = "pH")+
  scale_size(range = c(.01, 10), name="LOI")

```
**Use agrii_taxa_conversion_20240920.md script to 
- convert DADA2 assigned taxa to NCBI's taxaID; 
- run output through genomesizeR; 
- run output through SNB calculator; 
- run output though ANNA16 rRNA copy number estimator
**

now import the individual taxa/ASVs genome size, rRNA and SNB for calculating each sample's 
weighted mean values of gc (GC content), gs (genome size), rrna (rRNA copy), and snb (social niche breadth)

```{r}
rename=dplyr::rename
#import taxa's genome sizes 
img=read.delim(here("Analysis211014", "genome_size_new","IMG_taxa_genomesize.tsv")) %>%   select(-c("taxon_oid","Sequencing.Status", "Genome.Name...Sample.Name", 
            "Sequencing.Center", "Study.Name" ))
#import the taxa table, post rarefication
tax=read.csv(here("Analysis211014", "microeco_rarefied_taxa_updated_silva.csv")) 

#Code to run through the img table, taxonomic level by taxonomic level to extract 

# Step 1: Prepare the img table by selecting relevant columns
img_table_selected <- img %>%
  select(NCBI.Domain, NCBI.Phylum, NCBI.Class, NCBI.Order, NCBI.Family, NCBI.Genus, NCBI.Species,
         GTDB.Domain, GTDB.Phylum, GTDB.Class, GTDB.Order, GTDB.Family, GTDB.Genus, GTDB.Species,
         NCBI.Taxon.ID, `Genome.Size....assembled`, `Gene.Count....assembled`, 
         `CDS.Count....assembled`, `CDS......assembled`, `rRNA.Count....assembled`, 
         `X16S.rRNA.Count....assembled`, `GC....assembled`)

# Step 2: Create a function to lookup values from the img table and print query & result
lookup_value <- function(row, row_index, total_rows) {
  # Print progress for each row
  message(sprintf("Processing row %d of %d", row_index, total_rows))
  
  # Print the query (the current row being processed)
  print("Query:")
  print(row)
  
  # Perform the lookup
  match_row <- img_table_selected %>%
    filter(
        (NCBI.Genus == row$genus | GTDB.Genus == row$genus) |
        (NCBI.Family == row$family | GTDB.Family == row$family) |
        (NCBI.Order == row$order | GTDB.Order == row$order) |
        (NCBI.Class == row$class | GTDB.Class == row$class) |
        (NCBI.Phylum == row$phylum | GTDB.Phylum == row$phylum)
    ) %>%
    select(NCBI.Taxon.ID, `Genome.Size....assembled`, `Gene.Count....assembled`, 
           `CDS.Count....assembled`, `CDS......assembled`, `rRNA.Count....assembled`, 
           `X16S.rRNA.Count....assembled`, `GC....assembled`) %>%
    head(1)  # Only take the first match (if any)
  
  # If a match is found, print the match and return the original value + match data
  if (nrow(match_row) > 0) {
    print("Matched Result:")
    print(match_row)
    return(data.frame(row, match_row))  # Return as data frame
  } else {
    print("No match found.")
    return(data.frame(row, NCBI.Taxon.ID = NA, `Genome.Size....assembled` = NA, 
                      `Gene.Count....assembled` = NA, `CDS.Count....assembled` = NA, 
                      `CDS......assembled` = NA, `rRNA.Count....assembled` = NA, 
                      `X16S.rRNA.Count....assembled` = NA, `GC....assembled` =NA))  # Return NAs as data frame
  }
}

# Step 3: Apply the lookup function to the tax table with a progress bar
total_rows <- nrow(tax)
pb <- progress_bar$new(total = total_rows)

result_list <- list()  # Store results as a list to combine later

for (i in 1:total_rows) {
  pb$tick()  # Update progress bar
  result_list[[i]] <- lookup_value(tax[i, ], i, total_rows)  # Call the function with row index and total count
}

# Combine the results into a final table
result_table <- bind_rows(result_list)

# Step 4: View the result
fix(result_table)

#read in rarefied sequence table 
seq=read.csv(here("Analysis211014", "microeco_rarefied_seqtab_updated_silva.csv"))

#merge 
traits.seq=left_join(result_table, seq, by="X") %>% 
  left_join(.,snb, by="X") %>% 
  relocate(.,SNB.score,.after="X")
 
#write.csv(traits.seq,here("genome_size_new", "trait_seq_with_seqtab.csv"))

#reading in additional genomic traits and social niche breadth

#for rRNA file
rrna=read.csv(here("Analysis211014", "rRNA", "agri_copy_number.csv")) %>% 
  mutate(ASV = str_replace(index, ">","ASV_")) %>% 
  select(-index)
head(rrna)

#for SNB file
snb=read.delim(here("Analysis211014", "social_niche_bredth_new", "agrii_raw_snb_output.txt")) %>% 
  mutate(X = str_remove(taxonomic.lineage, ".*species\\.")) %>% 
  select(-`taxonomic.lineage`)
head(snb)

#now do a crazy merge of all traits into one mega file
genome.traits=left_join(traits.seq, snb, by = "X") %>% 
  left_join(., rrna, by= "ASV") #%>% 
  relocate(c(`number.of.samples`:predicted_copy_number), .after = "X16S.rRNA.Count....assembled")


#write.csv(genome.traits,here("genome_size_new", "genome_traits_with_seqtab.csv"))
#genome.traits=read.csv(here("Analysis211014", "genome_size_new", "genome_traits_with_seqtab.csv"))

###Calculations for weighted means.....
#first create long format table 
genome.traits.long=genome.traits %>% 
  pivot_longer(., cols = "AgriiRhiza.1.A.H.06.535502.2mm.16S.1":ncol(.),
               names_to = "samples", values_to = "counts")
select=dplyr::select

#create weighted snb averages
snb.weight=genome.traits.long %>% drop_na(SNB.score) %>% 
  group_by(ASV, samples) %>% 
  mutate(snb.weight=SNB.score*counts)

snb.final=snb.weight %>% dplyr::select(counts, samples, snb.weight) %>% 
  group_by(samples) %>% 
  mutate(ave.snb=sum(snb.weight)/sum(counts)) %>% 
  dplyr::select(samples, ave.snb) %>% 
  distinct()


#create weighted genome size averages
gs.weight=genome.traits.long %>% drop_na("Genome.Size....assembled") %>% 
  group_by(ASV, samples) %>% 
  mutate(gs.weight=`Genome.Size....assembled`*counts)

gs.final=gs.weight %>% select(counts, samples, gs.weight) %>% 
  group_by(samples) %>% 
  mutate(ave.gs=sum(gs.weight)/sum(counts)) %>% 
  select(samples, ave.gs) %>% 
  distinct()


#create weighted rrna averages
rrna.weight=genome.traits.long %>% drop_na(predicted_copy_number) %>% 
  group_by(ASV, samples) %>% 
  mutate(rrna.weight=predicted_copy_number*counts)

rrna.final=rrna.weight %>% select(counts, samples, rrna.weight) %>% 
  group_by(samples) %>% 
  mutate(ave.rrna=sum(rrna.weight)/sum(counts)) %>% 
  select(samples, ave.rrna) %>% 
  distinct()

#create weighted GC averages
gc.weight=genome.traits.long %>% drop_na(`GC....assembled`) %>% 
  group_by(ASV, samples) %>% 
  mutate(gc.weight=`GC....assembled`*counts)

gc.final=gc.weight %>% select(counts, samples, gc.weight) %>% 
  group_by(samples) %>% 
  mutate(ave.gc=sum(gc.weight)/sum(counts)) %>% 
  select(samples, ave.gc) %>% 
  distinct()

#merging weighted mean traits with metadata
#read in meta data and change indexing column to match those in traits files
# meta2=read.csv(here("Analysis211014", "AgriiRhizaEnv_v5_tidy.csv"), header = TRUE, check.names = F) %>% 
#   dplyr::rename(samples=vlookup_16s) %>% 
#   mutate(samples = str_replace_all(samples, "-","\\."))

meta6=meta5 %>% 
  rename(samples=rowname) %>% 
  mutate(samples = str_replace_all(samples, "-","\\."))


#combine all sample traits and meta data
mean.traits=left_join(snb.final, rrna.final, by = "samples") %>% 
  left_join(., gs.final, by = "samples") %>% 
  left_join(., gc.final, by = "samples") %>% 
  left_join(., meta6, by = "samples")
  

#write.csv(mean.traits,here("genome_size_new", "mean_traits2.csv"))

```


plot scatter graphs with lm equations
```{r}

#GS v rRNA
mean.traits %>% 
  ggplot(aes(x=ave.rrna, y=ave.gs, size = Organic_Matter_LOI, colour= pH))+
  geom_point(alpha=0.5) +
  scale_colour_gradient(low = "red", high="blue")+
  scale_size(range = c(.01, 10), name="LOI")+
  stat_poly_eq(use_label(c("eq", "adj.R2", "f", "p", "n")))


#GS v GC
mean.traits %>% 
  ggplot(aes(x=ave.gc, y=ave.gs, size = Organic_Matter_LOI, colour= pH))+
  geom_point(alpha=0.5) +
  scale_colour_gradient(low = "red", high="blue")+
  scale_size(range = c(.01, 10), name="LOI")+
  stat_poly_eq(use_label(c("eq", "adj.R2", "f", "p", "n")))

#GS v snb
mean.traits %>% 
  ggplot(aes(x=ave.snb, y=ave.gs, size = Organic_Matter_LOI, colour= pH))+
  geom_point(alpha=0.5) +
  scale_colour_gradient(low = "red", high="blue")+
  scale_size(range = c(.01, 10), name="LOI")+
  stat_poly_eq(use_label(c("eq", "adj.R2", "f", "p", "n")))

#GC v snb
mean.traits %>% 
  ggplot(aes(x=ave.snb, y=ave.gc, size = Organic_Matter_LOI, colour= pH))+
  geom_point(alpha=0.5) +
  scale_colour_gradient(low = "red", high="blue")+
  scale_size(range = c(.01, 10), name="LOI")+
  stat_poly_eq(use_label(c("eq", "adj.R2", "f", "p", "n")))


```


run generalised additive model to examine relationships between variables and weighted means
```{r}

#GAM for pH associations - are genome size and rRNA linked to pH?
Gam.pH=gam(pH~s(Organic_Matter_LOI)+s(ave.rrna)+s(ave.gs)+s(ave.gc),data = mean.traits)
summary(Gam.pH)


#GAM for snb - are pH, gs and rrna linked to SNB score?
Gam.snb=gam(ave.snb~s(Organic_Matter_LOI)+s(pH)+s(ave.gs)+s(ave.rrna)+s(ave.gc),data = mean.traits)
summary(Gam.snb)

```


phylofactorisation - examining the phylogenetic tree to check for relationships with SNB 
and to make sure weighted mean SNB scores are not being skewed by increased occurrence of 'rare' taxa in higher pH soils

```{r}
# loading files
asv = read.csv(here("Analysis211014", "trait_seq_with_seqtab.csv")) %>% 
  select(-`X.1`) %>% 
  relocate(ASV)
```

Creating a multifasta file to build tree
```{r}
multifasta <- asv %>%
  # Create a new column combining the seq_id and sequences into FASTA format
  mutate(fasta_format = paste0(">", ASV, "\n", X)) %>%
  
  # Pull the fasta_format column as a vector
  pull(fasta_format) %>%
  
  # Write the FASTA formatted sequences to a file
  write_lines(here("agrii_multifasta.fasta"))
```

Building a tree
```{r}
library(msa)
library(phangorn)
library(DECIPHER)

# align using 'msa' package - taking too long so using DECIPHER - see below.
# alignment <- msa(type = "dna", asv$X, method = "ClustalOmega", verbose = TRUE)  # Assuming the sequences are in column X

# Read the FASTA file
dna_sequences <- readDNAStringSet(here("Analysis211014", "agrii_multifasta.fasta"), format = "fasta")
dna_sequences 

# Perform the alignment using AlignSeqs from 'DECIPHER'
alignment <- AlignSeqs(dna_sequences, processors = NULL, verbose = TRUE) # processors = NULL, detects all and uses them
```

```{r}
# constructing a distance matrix and tree
dist_matrix <- dist.ml(as.DNAbin(alignment))
tree <- NJ(dist_matrix)  # Neighbor-Joining Tree
plot(tree)
#write.tree(tree, file = here("agrii_rarefied_asv_tree.nwk"))
```

Maximum-likelihood tree
```{r}
# Create a phyDat object for DNA sequences
aa_phyDat <- phyDat(as.DNAbin(alignment), type = "DNA")

# Estimate an ML tree using the GTR model
fit <- pml(tree, aa_phyDat)
fit_opt <- optim.pml(fit, model = "GTR")  # Optimize the ML tree

# Plot the optimized ML tree
plot(fit_opt$tree, main = "Maximum Likelihood Tree")
```

```{r}
# Visualize the tree with tip labels
plot(fit_opt$tree, main = "agrii rarefied ASVs", cex = 0.6, edge.width = 2)
#write.tree(fit_opt$tree, file = here("agrii_optimised_ML_tree.nwk"))
```

Preparing files for `Phylofactor`
```{r}
# setting up the taxonomy file
taxonomy <- asv %>% select(ASV, domain:genus) %>%
  unite(col = "taxonomy", domain:genus, sep = ";", na.rm = TRUE) %>% 
  rename("OTU_ID" = "ASV")

# generating the OTU table
otu <- asv %>% select(ASV, contains("AgriiRhiza")) %>% 
  column_to_rownames(., var = "ASV")

# reading in the metadata
meta <- read.csv(here("Analysis211014","metadata_agrii.csv")) %>% 
  select(vlookup_16s, pH, Organic_Matter_LOI) %>% 
  rename("LOI" = "Organic_Matter_LOI") %>% 
  mutate(samples = str_replace_all(vlookup_16s, "-", "\\.")) %>% 
  select(-vlookup_16s) %>% 
  mutate(group = case_when(
    pH < 7 ~ "low-neutral",    # Rename "neutral" to "low" for consistency
    pH > 7 ~ "high"  # Handle pH > 7
  )) %>% 
  column_to_rownames(., var = "samples")

# Filter meta to include only rows with samples matching the column names in otu
meta_filtered <- meta %>%
  rownames_to_column(., var = "samples") %>% 
  filter(samples %in% colnames(otu)) %>% 
  column_to_rownames(., var = "samples")
```

Testing phylogenetic signal with pH
To test whether `pH` values from the metadata are driving the separation in your phylogenetic tree, you can use methods like **phylogenetic signal tests** or **ordination techniques** (e.g., PCoA) that integrate phylogenetic information. Below is an outline of steps you can follow:

1: Phylogenetic Signal (Blomberg’s K or Pagel’s Lambda)
These tests assess whether the `pH` values exhibit a phylogenetic signal, i.e., whether closely related ASVs tend to have similar `pH` values.
```{r}
library(phytools)
library(phylosignal)

# setting the tree variable
# tree <- fit_opt$tree  # OR 
tree <- read.tree(here("Analysis211014", "agrii_optimised_ML_tree.nwk"))

# reading in the mean_traits per csv
snb_meta <- asv %>% select(ASV, SNB.score) %>% rename("SNB" = "SNB.score")

# Ensure that `metadata` and the `tree$tip.label` have the same ASV names:
metadata <- snb_meta[match(tree$tip.label, snb_meta$ASV), ]
rownames(metadata) <- snb_meta$ASV

# Extract SNB values corresponding to the tree tips
SNB_values <- metadata$SNB

# Test for SNB phylogenetic signal (Blomberg’s K or Pagel’s Lambda)
signal_test <- phylosig(tree, SNB_values, method="lambda", test=TRUE)
print(signal_test)
```

```{r}
# Alternatively, for Blomberg’s K:
library(picante)
library(phylobase)
# convert tree to phylo4 format
phylo_tree <- as(tree, "phylo4")
names(SNB_values) <- tree$tip.label

# Remove NAs from SNB_values and keep only matching tree tips
non_na_indices <- !is.na(SNB_values)
SNB_values_filtered <- SNB_values[non_na_indices]

# Subset the tree to match SNB_values (only the tips present in SNB_values_filtered)
tree_filtered <- drop.tip(tree, tree$tip.label[!non_na_indices])

# Check if the tree is rooted
if (!is.rooted(tree_filtered)) {
  # Root the tree (using midpoint rooting if no outgroup is available)
  tree_filtered <- midpoint.root(tree_filtered)
}

# Ensure the tree is fully dichotomous
if (!is.binary.tree(tree_filtered)) {
  # Resolve polytomies
  tree_filtered <- multi2di(tree_filtered)
}

# Run phylosignal again
K <- phylosignal(SNB_values_filtered, tree_filtered)

# View the results
print(K$K)
print(K)
```

2: Phylogenetic Comparative Methods with `phylolm`
Using **phylogenetic generalized least squares (PGLS)** to model the relationship between `SNB` and the tree while accounting for phylogenetic relationships.
```{r}
library(phylolm)

# Run a PGLS regression:
model <- phylolm(SNB_values ~ 1, data = metadata, phy = tree)
summary(model)
```

3: Phylogenetic Ordination (PCoA on Unifrac Distances)
Calculate **Unifrac distances** (which account for phylogenetic relationships) and perform **PCoA** to see how the tree is separating based on `SNB`.
```{r}
# Calculate UniFrac distances:
library(phyloseq)

# creating phyloseq object
ps_otu <- otu_table(as.matrix(otu), taxa_are_rows = TRUE)  # 'otu' is your OTU data
ps_meta <- sample_data(meta_filtered)

taxonomy_edited <- taxonomy %>% column_to_rownames(., var = "OTU_ID")
ps_taxonomy <- tax_table(as.matrix(taxonomy_edited))

physeq <- phyloseq(ps_otu, ps_meta, tree, ps_taxonomy)

# Assuming your tree and ASV abundance table are already imported in phyloseq
unifrac_dist <- UniFrac(physeq, weighted = TRUE)

# PCA
pcoa_results <- ape::pcoa(unifrac_dist)
plot(pcoa_results$vectors[, 1], pcoa_results$vectors[, 2],
        col = meta_filtered$pH, pch = 19, 
        xlab = "PCoA1", ylab = "PCoA2")

# Test correlation between SNB and ordination:
cor.test(pcoa_results$vectors[, 1], meta_filtered$pH)
```

```{r}
# Create a data frame with the PCoA results and metadata
pcoa_df <- data.frame(PC1 = pcoa_results$vectors[, 1], 
                      PC2 = pcoa_results$vectors[, 2],
                      pH = meta_filtered$pH)

ggplot(pcoa_df, aes(x = PC1, y = PC2, color = pH)) +
  geom_point(size = 3, alpha = 0.7) +  # Plot points with transparency
  scale_color_gradient(low = "red", high = "blue") +  # Color gradient based on pH values
  labs(x = "PCoA1", y = "PCoA2", color = "pH") +  # Labels
  theme_minimal() +  # Minimalistic theme
  theme(legend.position = "right") +  # Legend position
  ggtitle("PCoA Plot Based on Weighted UniFrac Distances")
```

```{r}
# Correlation plot with ggplot2
cor_plot <- ggplot(pcoa_df, aes(x = pH, y = PC1)) +
  geom_point(size = 3, alpha = 0.7, color = "darkblue") +  # Scatter plot
  geom_smooth(method = "lm", color = "red") +  # Add linear regression line
  labs(x = "pH", y = "PCoA1", title = "Correlation between PCoA1 and pH") +
  theme_minimal() +  # Minimal theme
  theme(plot.title = element_text(hjust = 0.5))

# Display the correlation plot
print(cor_plot)

# Optional: Show correlation value
cor_test <- cor.test(pcoa_df$PC1, pcoa_df$pH)
print(cor_test)
```

```{r}
# Create a data frame with the PCoA results and metadata
pcoa_df <- data.frame(PC1 = pcoa_results$vectors[, 1], 
                      PC2 = pcoa_results$vectors[, 2],
                      pH = meta_filtered$pH)

# Fit a linear model for PCoA1 vs pH
lm_model <- lm(PC1 ~ pH, data = pcoa_df)

# Extract R^2 and p-value
summary_lm <- summary(lm_model)
r_squared <- summary_lm$r.squared
p_value <- summary_lm$coefficients[2, 4]

# Create correlation plot with R^2 and p-value annotated
cor_plot <- ggplot(pcoa_df, aes(x = pH, y = PC1)) +
  geom_point(size = 3, alpha = 0.7, color = "darkblue") +  # Scatter plot
  geom_smooth(method = "lm", color = "red") +  # Add linear regression line
  labs(x = "pH", y = "PCoA1", title = "Correlation between PCoA1 and pH") +
  theme_minimal() +  # Minimal theme
  theme(plot.title = element_text(hjust = 0.5)) +
  annotate("text", x = min(pcoa_df$pH), y = max(pcoa_df$PC1), 
           label = paste0("R² = ", round(r_squared, 3), "\n", 
                          "p = ", format(p_value, digits = 3, scientific = TRUE)),
           hjust = 0, size = 5, color = "black")

# Display the correlation plot
print(cor_plot)
```

```{r}
# Load vegan for ordisurf
library(vegan)
library(viridis)

# Prepare data for ordisurf
ordisurf_plot <- ordiplot(pcoa_results$vectors, choices = c(1, 2), display = "sites", type = "n")

# Add points for PCoA with colors based on pH values
points(ordisurf_plot, "sites", pch = 19, col = viridis(20)[cut(pcoa_df$pH, breaks = 20)])

# Fit pH as a surface using ordisurf and add contour lines
ordisurf(ordisurf_plot, pcoa_df$pH, add = TRUE, col = "blue", labcex = 0.8)

# Add a color legend
legend("topright", legend = round(seq(min(pcoa_df$pH), max(pcoa_df$pH), length.out = 5), 2),
       fill = viridis(5), title = "pH", bty = "n")

# Add contour lines for pH using ordisurf
ordisurf(ordisurf_plot, pcoa_df$pH, add = TRUE, col = "red", labcex = 1.1)

# Add legend
legend("topright", legend = paste("pH:", round(range(pcoa_df$pH), 2)), col = "red", pch = "-")

library(cowplot)

# Combine the plots into one figure
combined_plot <- plot_grid(cor_plot, ordisurf_plot, labels = c("A", "B"), ncol = 2)

# Display the combined plot
print(combined_plot)
```

4: Mantel Test (Correlation Between Phylogeny and pH)
**Mantel test** to see if there is a significant correlation between the phylogenetic distances and pH.
```{r}
# Create a distance matrix of SNB values:
pH_dist <- dist(meta_filtered$pH)

# Perform Mantel test:
mantel_result <- vegan::mantel(unifrac_dist, pH_dist)
print(mantel_result)
```



## Phylofactorisation
```{r}
# Step 1: Load data and required packages
library(phylofactor)
library(phyloseq)
library(ape)
library(ggtree)

# Step 1: Replace 0 with 0.65 * minimum non-zero counts for that sample
otu_no_zero <- otu %>%
  mutate(across(everything(), ~ ifelse(. == 0, 0.65, .)))

meta_pf <- read.csv(here("Analysis211014", "mean_traits.csv")) %>% 
  select(X, samples, ave.snb, ave.rrna, ave.gs, ave.gc,
         simps.div, alpha.div, pH, Organic_Matter_LOI) %>%
  column_to_rownames(., var = "samples") %>%
  select(-X) %>%
  mutate(snb_group = case_when(          # creating SNB groups since median SNB value is 0.9228
    ave.snb < 0.923 ~ "low_snb", #based on median sample SNB value
    ave.snb >= 0.923 ~ "high_snb"
  ))

# Step 2: Create a phyloseq object
ps_otu <- otu_table(as.matrix(otu_no_zero), taxa_are_rows = TRUE)
ps_meta <- sample_data(meta_pf)
ps_taxonomy <- tax_table(as.matrix(taxonomy_edited))

physeq <- phyloseq(ps_otu, ps_meta, tree, ps_taxonomy)

# Step 3: Extract response variable (e.g., ave.snb)
response_variable <- sample_data(physeq)$ave.snb    # sample_data(physeq)$ave.snb

# Step 4: Ensure the OTU table is in the same order as the tree's tip labels
otu_no_zero_ordered <- otu_no_zero[tree$tip.label, , drop = FALSE]  # Reorder OTU table

# Step 5: Perform phylofactorization
pf <- PhyloFactor(otu_no_zero_ordered, tree, ps_meta, 
                  frmla = Data~`ave.snb`, 
                  transform.fcn = log,
                  nfactors = 3, 
                  ncores=12)
saveRDS(pf, here("pf.rds"))
pf

# Visualize the clades associated with the response variable
plot(pf)

# Summarize the results
summary(pf)
```

```{r}
# Extract important clades and display their names
clades <- pf$groups
print(clades)

# plot
library(ggtree)

# Create an empty list to store results
linked_asvs <- list(Group1 = character(), Group2 = character())

# Link ASVs for Group1
for (node in clades[[1]]$Group1) {
    # Check if the node number corresponds to an ASV
    asv_name <- names(node_numbers)[node_numbers == node]
    if (length(asv_name) > 0) {
        linked_asvs$Group1 <- c(linked_asvs$Group1, asv_name)
    }
}

# Link ASVs for Group2
for (node in clades[[1]]$Group2) {
    # Check if the node number corresponds to an ASV
    asv_name <- names(node_numbers)[node_numbers == node]
    if (length(asv_name) > 0) {
        linked_asvs$Group2 <- c(linked_asvs$Group2, asv_name)
    }
}

# Print linked ASVs for both groups
print(linked_asvs)

# Retrieve node numbers for Group1 and Group2 as numeric
group1_nodes <- as.numeric(sapply(clades[[1]]$Group1, function(node) {
    which(node_numbers == node)
}))

group2_nodes <- as.numeric(sapply(clades[[1]]$Group2, function(node) {
    which(node_numbers == node)
}))

# Check for valid node numbers
group1_nodes <- group1_nodes[!is.na(group1_nodes)]
group2_nodes <- group2_nodes[!is.na(group2_nodes)]

# Plot the tree with highlighted nodes
ggtree(tree, aes(color = "grey", linetype = "dashed")) + 
  geom_tree(color = "grey") +  # Set the edges to grey
  geom_hilight(node = group1_nodes, fill = "blue") + 
  geom_hilight(node = group2_nodes, fill = "red") +
  theme(legend.position = "none")  # Optionally remove the legend if it's not needed
```

RUNNING WITH LOW/HIGH SNB
````{r}
# Step 5: Perform phylofactorization with 'snb_group'
pf_snb <- PhyloFactor(otu_no_zero_ordered, tree, ps_meta, 
                  frmla = Data~snb_group, 
                  transform.fcn = log,
                  nfactors = 3, 
                  ncores=12)
saveRDS(pf_snb, here("pf_snb.rds"))
pf_snb

# Visualize the clades associated with the response variable
plot(pf_snb)

# Summarize the results
summary(pf_snb)
```

```{r}
# Extract important clades and display their names
clades_pf_snb <- pf_snb$groups
print(clades_pf_snb)

# Create an empty list to store results
linked_asvs_pf_snb <- list(Group1 = character(), Group2 = character())

# Link ASVs for Group1
for (node in clades_pf_snb[[1]]$Group1) {
    # Check if the node number corresponds to an ASV
    asv_name <- names(node_numbers)[node_numbers == node]
    if (length(asv_name) > 0) {
        linked_asvs_pf_snb$Group1 <- c(linked_asvs_pf_snb$Group1, asv_name)
    }
}

# Link ASVs for Group2
for (node in clades_pf_snb[[1]]$Group2) {
    # Check if the node number corresponds to an ASV
    asv_name <- names(node_numbers)[node_numbers == node]
    if (length(asv_name) > 0) {
        linked_asvs_pf_snb$Group2 <- c(linked_asvs_pf_snb$Group2, asv_name)
    }
}

# Print linked ASVs for both groups
print(linked_asvs_pf_snb)

# Retrieve node numbers for Group1 and Group2 as numeric
group1_nodes_pf_snb <- as.numeric(sapply(clades_pf_snb[[1]]$Group1, function(node) {
    which(node_numbers == node)
}))

group2_nodes_pf_snb <- as.numeric(sapply(clades_pf_snb[[1]]$Group2, function(node) {
    which(node_numbers == node)
}))

# Check for valid node numbers
group1_nodes_pf_snb <- group1_nodes_pf_snb[!is.na(group1_nodes_pf_snb)]
group2_nodes_pf_snb <- group2_nodes_pf_snb[!is.na(group2_nodes_pf_snb)]

# Plot the tree with highlighted nodes
ggtree(tree, aes(color = "yellow", linetype = "dashed")) + 
  geom_tree(color = "yellow") +  # Set the edges to grey
  geom_hilight(node = group1_nodes_pf_snb, fill = "blue") + 
  geom_hilight(node = group2_nodes_pf_snb, fill = "red") +
  theme(legend.position = "none")  # Optionally remove the legend if it's not needed
```
TESTING WITH UNIFRAC
```{r}
# Assuming your tree and ASV abundance table are already imported in phyloseq
unifrac_dist <- UniFrac(physeq, weighted = TRUE)

# PCA
pcoa_results <- ape::pcoa(unifrac_dist)
plot(pcoa_results$vectors[, 1], pcoa_results$vectors[, 2],
        col = meta_pf$ave.snb, pch = 19, 
        xlab = "PCoA1", ylab = "PCoA2")

# Test correlation between SNB and ordination:
cor.test(pcoa_results$vectors[, 1], meta_pf$ave.snb)

# Create a data frame with the PCoA results and metadata
pcoa_df <- data.frame(PC1 = pcoa_results$vectors[, 1], 
                      PC2 = pcoa_results$vectors[, 2],
                      snb = meta_pf$ave.snb)

ggplot(pcoa_df, aes(x = PC1, y = PC2, color = snb)) +
  geom_point(size = 3, alpha = 0.7) +  # Plot points with transparency
  scale_color_gradient(low = "red", high = "blue") +  # Color gradient based on pH values
  labs(x = "PCoA1", y = "PCoA2", color = "snb") +  # Labels
  theme_minimal() +  # Minimalistic theme
  theme(legend.position = "right") +  # Legend position
  ggtitle("PCoA Plot Based on Weighted UniFrac Distances")
```

```{r}
# Correlation plot with ggplot2
cor_plot <- ggplot(pcoa_df, aes(x = snb, y = PC1)) +
  geom_point(size = 3, alpha = 0.7, color = "darkblue") +  # Scatter plot
  geom_smooth(method = "lm", color = "red") +  # Add linear regression line
  labs(x = "snb", y = "PCoA1", title = "Correlation between PCoA1 and pH") +
  theme_minimal() +  # Minimal theme
  theme(plot.title = element_text(hjust = 0.5))

# Display the correlation plot
print(cor_plot)

# Optional: Show correlation value
cor_test <- cor.test(pcoa_df$PC1, pcoa_df$snb)
print(cor_test)

# from here: https://chrismitbiz.github.io/ABlab-workflows/phylofactor.html#phylofactor
# create/extract metadata object from phyloseq object
meta.df <- data.frame(sample_data(ps.phylo))

# Extract rownames (Sample IDs)
rnames <- rownames(meta.df)

# Explanatory variable for model
X <- as.factor(meta.df$snb_group)

# check explanatory
summary(X)

# name the vector data with sample IDs 
names(X) <- rnames

# create an ASV table for phylofactor 
Data = as(otu_table(ps.phylo), "matrix")
Data = as.data.frame(Data)

# create tree object from phyloseq object for phylofactor
tree <- phy_tree(ps.phylo) 

# unroot the tree, this is critical for phylofactor to work
tree_ur <- ape::unroot(tree)

# check if unrooted
ape::is.rooted(tree_ur) # FALSE IS GOOD

# check that all in order
all(rownames(Data) == row.names(tax_table(ps.phylo))) 
```

```{r}
# Run phylofactor
PF_new <- PhyloFactor(Data, tree_ur, X, 
        frmla=Data~ X,    # Data (ASV abundances) as a response to X (snb_group)
        ncores = 12,      # number of cores 
        stop.early = NULL, 
        transform.fcn=log,    # log-ratio transform
        KS.Pthreshold = 0.01, # Kolmogorov Smirnov significance level 
        nfactors = 15, # the first 13 phylofactors only - in interest of time 
        choice = "var")  # see notes above
saveRDS(PF_new,'pf_new.rds')  

# Overview
PF_new # showing the output of the PF_new
PF_new$factors
```

```{r}
## Taxa summaries
# To create a taxa summary a taxonomy dataframe is needed first, 
# whose first column contains the species in the phylogeny and whose
# second column contains a semicolon-delimited taxonomy. 
# Here, we extract that out of the phyloseq object as following.
taxonomy_pf_new <- data.frame(tax_table(ps.phylo)) %>%
  unite(Taxon, sep=";", c("domain","phylum","class","order","family","genus"))  %>%
  rownames_to_column("Feature.ID")  
 
# running pf.taxa function from phylofactor
phylofactor::pf.taxa(PF_new, taxonomy_pf_new, factor = 6)
```

```{r}
# create a summary for a factor, this time for factor 6
s <- summary(PF_new, taxonomy_pf_new, factor=6)
s
# List all taxa for each group
s$taxa.split
```

```{r}
## Summary tables
# Following a custom function to provide a dataframe of all ASVs and their phylofactors 
# It requires as input the phylofactor and phyloseq objects. 
# It needs to be exactly the same pholoseq object (contain the ASVs and taxa) 
# that was used prior to running PhyloFactor. 
# The number (5 or 4) indicates which `choice` ('var' or 'F') was provided 
# to the PhyloFactor function. Either 'var' (5) or 'F' (4). Just my poor coding here.   
PF_summary <- phyfacsummary(PF_new, ps.phylo, 5)  
head(PF_summary)
```

Visualising factors
```{r}
## ILR boxplots
# If needed add another element to this plot. The boxplots of ILRs. 
# Use the custom addILRtophyloseq function and create a new phyloseq object with 
# added ILR results in the sample_data (for plotting)
ps.ILR <- addILRtophyloseq(ps.phylo, PF_new)
# Quick Check that they have 'arrived'.
head(data.frame(sample_data(ps.ILR)))
```

```{r}
## Create plotdata
# Some vectors to filter and change the appearance of the plot
factorvector <- paste0('Phylofactor.', 1:15)  # This creates a vector of "Phylofactor.1" to "Phylofactor.15"
labelvec <- paste0('Phylofactor ', 1:15)  # Creates labels like "Phylofactor 1", "Phylofactor 2", etc.

# rearrange the sample data for plotting
plotdata <- data.frame(sample_data(ps.ILR)) %>% 
  pivot_longer(cols = starts_with("Phylofactor"), 
    names_to = "ILR")

# Select specific factor vectors for plotting using the factorvector to filter
plotdata_sub <- plotdata %>% 
  filter(ILR %in% factorvector)

# create factors of the ILR column 
plotdata_sub$ILR <- factor(plotdata_sub$ILR, levels = plotdata_sub[1:length(unique(plotdata_sub$ILR)), ]$ILR)

# create ggplot object
pilr <- plotdata_sub %>% 
  ggboxplot(x = "snb_group", y = "value", 
    combine = TRUE, 
    facet.by = "ILR", 
    fill = "ILR", 
    ncol = 2, 
    alpha = 0.5, 
    ylab = ("Isometric log ratio (ILR)"), 
    panel.labs.font = list(size = 11), 
    panel.labs = list(ILR = labelvec)) + 
    geom_point(position = "jitter", alpha = 0.5) +  # add jitter 
    theme(legend.position = "none") +               # do not show legend
    theme(axis.text.x = element_text(angle = 90, hjust = 1, 
    vjust = 0.5), axis.title.x = element_blank(), axis.text = element_text(size = 12))

# Choose custom colors with + scale_fill_manual(values = yourcolvector)
pilr 
```

Basic tree with factors highlighted in colors
```{r}
gtree_new <- phylofactor::pf.tree(PF_new)
ggtree::rotate_tree(gtree_new$ggplot,-45)
```

```{r}
## Create color vector
# use custom prevalence table function and select phyloseq object and phylum-level
taxa.df.phylum <- prevalencedf(ps.phylo, phylum)
# set seed differently will change the random sampling of the colours
set.seed(20240929)
# randomly sample colours from this colour scheme
#treecols <- sample(colorspace::rainbow_hcl(nrow(unique(taxa.df.phylum))))
# Other color options
treecols <- sample(viridis::turbo(nrow(unique(taxa.df.phylum)), 
                                  direction = -1, begin = 0, end = 1))
# treecols <- viridis::rainbow(nrow(unique(taxa.df.phylum)), direction = -1)
names(treecols) <- unique(taxa.df.phylum$phylum)
rm(taxa.df.phylum)
treecols
```

```{r}
# Skip this and do the next code chunk if you just want to use the custom function

# Use custom prevalence table function and select phyloseq object and phylum-level
taxa.df.phylum <- prevalencedf(ps.phylo, phylum)

# Choose the phyla with highest ASV counts
# This vector will determine the length and order of the phyla shown in the legend
# change this to your needs. 
# In this case we fiter n > 3 and create a vector of the Phylum column
ordervec <- (taxa.df.phylum %>% 
            dplyr::filter(n > 3))$phylum
# You could also create a custom vector. For example: 
#ordervec <- c("Proteobacteria" ,  "Firmicutes" , "Actinobacteriota",
#              "Bacteroidota",  "Chloroflexi",  "Acidobacteriota")
  
## Plot the tree
# setting the variables
layout = "circular"
branch.length = "none"
metadatacolumn <- 'snb_group'
offset.text = 5
pwidth = 0.05
factorvector = c(1:15)  # choose which factors to highlight

# 'Copy' tree object from the phylofactor object into the phyloseq object for plotting
phyloseq::phy_tree(ps.phylo) <- PF_new$tree
# create long dataframe of metadata and select the columns OTU and Abundance
# then filter only unique rows 
melt_simple <- psmelt(ps.phylo) %>% 
  dplyr::select(OTU, Abundance) %>% 
  unique()

# Create a phyla list to include into the tree for colouring
taxa.df.phylum <- phyloseq::tax_table(ps.phylo)@.Data %>% 
  as.data.frame %>% 
  rownames_to_column("OTUID") %>% 
  dplyr::select(OTUID, phylum) %>% 
  column_to_rownames("OTUID")

phyla.list <- split(rownames(taxa.df.phylum), 
                    taxa.df.phylum$phylum, 
                    drop = TRUE)

tree.df <- phyloseq::phy_tree(PF_new$tree)
ggtree_gps <- ggtree::groupOTU(tree.df, phyla.list, "phylum")

# saving the previous plot
saveRDS(p, "plot_pf.rds")

# gg plotting
p <- ggtree::ggtree(ggtree_gps, aes(color = phylum), layout = layout, 
        branch.length = branch.length) + 
        theme(legend.position = "right") + 
        scale_colour_manual(values = treecols, limits = ordervec, 
        breaks = ordervec, 
        guide = guide_legend(override.aes = list(size = 5), 
        keywidth = 0.2, keyheight = 0.2, order = 1))

# add barplots on the outside, representing relative abundances
p <- p + ggnewscale::new_scale_fill() + 
         geom_fruit(data = melt_simple, 
          geom = geom_bar, 
          mapping = aes(y = OTU, x = Abundance, group = label, 
          fill = phylum), pwidth = 0.38, orientation = "y", 
          stat = "identity", show.legend = FALSE) + 
          scale_fill_manual(values = treecols, 
          limits = ordervec, 
          breaks = ordervec)

## add phylofactor highlights in grey 
factor.map = data.frame(Phylofactor = 1:PF_new$nfactors, group = rep(1, PF_new$nfactors))
m <- nrow(factor.map)
n = ape::Ntip(PF_new$tree)
nd <- NULL
for (i in 1:m) {
    grp <- PF_new$tree$tip.label[PF_new$groups[[factor.map[i,1]]][[factor.map[i, 2]]]]
    grp <- intersect(grp, PF_new$tree$tip.label)
    nd <- c(nd, tidytree::MRCA(PF_new$tree, grp))
}

df <- data.frame(factor.map, node = nd)
df$Phylofactor <- factor(df$Phylofactor, levels = df$Phylofactor)

p <- p + ggtree::geom_hilight(data = df, aes(node = node), 
          fill = "black", alpha = 0.3)

# Add a text label to identify the number and location of the phylofactor
  for (i in factorvector) {
p <- p + geom_cladelabel(node = 
                 (df %>% dplyr::filter(Phylofactor ==  i))$node, label = paste("PF", 
                 (df %>% dplyr::filter(Phylofactor == i))$Phylofactor), 
                 offset.text = offset.text, 
                 hjust = "center")
  }
rm(taxa.df.phylum, df,factor.map)

# plot the PF tree
p
```

Highlight the taxa
```{r}
# check which phyla to highlight in tree
s$taxa.split$group2 # this is the summary for factor 6 create above

# create ordervec
df <- prevalencedf(ps.phylo, phylum)
ordervec <- (df %>% 
  dplyr::filter(phylum %in% c("Proteobacteria", "Actinobacteriota", "Firmicutes", 
                              "Bacteroidota",  "Verrucomicrobiota", "Planctomycetota",
                              "Myxococcota", "Acidobacteriota",
                              "Chloroflexi","Cyanobacteria")))
# # plot tree
# # Requires a vector called 'treecols'
# p <- mytreefunctionwithbarplots(PF_new, ps.phylo, 
#          layout = "circular",
#          branch.length = "none",
#          factorvector = c(1,2,3,4,6,11,13),
#          offset.text = 5,
#          pwidth = 0.05)
```

```{r}
# A small df to map the factors and groups
factor.map = data.frame(Phylofactor = 1:PF_new$nfactors, group = rep(1, PF_new$nfactors))
# Then get tip labels of selected factors and groups and use as index to the # taxonomy to get taxon names
# Create a df of taxa names for annotation - make sure that each taxon appears only once
# here the aim is to show the first common taxon level for each phylofactor
tax.names <- data.frame( Taxa = c(data.frame(tax_table(ps.phylo))[PF_new$tree$tip.label
                           [PF_new$groups[[factor.map[1,1]]][[factor.map[1, 2]]]],]$genus,
                           
                           data.frame(tax_table(ps.phylo))[PF_new$tree$tip.label
                           [PF_new$groups[[factor.map[2,1]]][[factor.map[2, 2]]]],]$family,
                           
                           data.frame(tax_table(ps.phylo))[PF_new$tree$tip.label
                           [PF_new$groups[[factor.map[4,1]]][[factor.map[4, 2]]]],]$family,
                           
                           data.frame(tax_table(ps.phylo))[PF_new$tree$tip.label
                           [PF_new$groups[[factor.map[6,1]]][[factor.map[5, 2]]]],]$genus,
                           
                           data.frame(tax_table(ps.phylo))[PF_new$tree$tip.label
                           [PF_new$groups[[factor.map[11,1]]][[factor.map[11, 2]]]],]$family) ) %>% 
                      unique()
# get node labels and add to factor.map
m <- nrow(factor.map)
nd <- NULL
for (i in 1:m) {
    grp <- PF_new$tree$tip.label[PF_new$groups[[factor.map[i,1]]][[factor.map[i, 2]]]]
    grp <- intersect(grp, PF_new$tree$tip.label)
      nd <- c(nd, tidytree::MRCA(PF_new$tree, grp))
}
factor.map <- data.frame(factor.map, node = nd)

# add taxon annotations using the factor map
p <- p + geom_cladelabel(node = 
                 (factor.map %>% dplyr::filter(Phylofactor ==  1))$node, 
                 label = paste('italic(', tax.names[1,],')'), 
                 offset.text = 20, 
                 hjust = "center",
                 fontsize = 3,
                 parse = TRUE)

p <- p + geom_cladelabel(node = 
                 (factor.map %>% dplyr::filter(Phylofactor == 2))$node, 
                 label = paste('italic(', tax.names[2,],')'), 
                 offset.text = 10, 
                 hjust = "center",
                 fontsize = 3,
                 parse = TRUE)

p <- p + geom_cladelabel(node = 
                 (factor.map %>% dplyr::filter(Phylofactor == 4))$node, 
                 label = paste('italic(', tax.names[3,],')'), 
                 offset.text = 10, 
                 hjust = "center",
                 fontsize = 3,
                 parse = TRUE)

p <- p + geom_cladelabel(node = 
                 (factor.map %>% dplyr::filter(Phylofactor == 6))$node, 
                 label = paste('italic(', tax.names[4,],')'), 
                 offset.text = 10, 
                 hjust = "center",
                 fontsize = 3,
                 parse = TRUE)

p <- p + geom_cladelabel(node = 
                 (factor.map %>% dplyr::filter(Phylofactor == 11))$node, 
                 label = paste('italic(', tax.names[5,],')'), 
                 offset.text = 10, 
                 hjust = "center",
                 fontsize = 3,
                 parse = TRUE)
p
```

Basic tree with colors
```{r}
layout = "circular"
factor.map = data.frame(PF=1:PF_new$nfactors, 
                        group = rep(1, PF_new$nfactors),
                        Phylofactor = paste0("Phylofactor.", 1:PF_new$nfactors))
# Use this colour vector for the factor highlights and the boxplots
cols <- viridis::viridis(nrow(factor.map), direction = 1)
names(cols) <- factor.map$Phylofactor

m <- nrow(factor.map)
n = ape::Ntip(PF_new$tree)
nd <- NULL
for (i in 1:m) {
  grp <- PF_new$tree$tip.label[PF_new$groups[[factor.map[i,1]]][[factor.map[i,2]]]]
  grp <- intersect(grp, PF_new$tree$tip.label)
  nd <- c(nd, tidytree::MRCA(PF_new$tree, grp))
}

df <- data.frame(factor.map, node = nd)
df$Phylofactor <- factor(df$Phylofactor, 
                         levels = df$Phylofactor)

p <- ggtree::ggtree(PF_new$tree,
                    layout=layout, 
                    aes(alpha = 0.3)) +
      ggtree::geom_hilight(data = df, 
                    aes(node = node, 
                    fill = Phylofactor)) + 
      scale_fill_manual(values = cols,
                    limits = c('Phylofactor.6', 'Phylofactor.11', 'Phylofactor.13'),
                    labels = c('Phylofactor 6','Phylofactor 11','Phylofactor 13') ) +
  theme(legend.title = element_blank(), 
        legend.position = "top") +
  guides(alpha = "none")
p
```

```{r}
# Use the custom addILRtophyloseq function and create a new phyloseq object with 
# added ILR results in the sample_data (for plotting)
ps.ILR <- addILRtophyloseq(ps.phylo, PF_new)
PF <- PF_new
p.ILR <- ILR_plotfunction(ps.ILR,
                factorvector = factorvector,  # choose factors 
                explanatory = "snb_group",  # metadata used in model
                labelvec = labelvec,
                ncol = 2)  + # of the facets (if more than one facet) 
          scale_fill_manual(values = cols)
# Combine plots
ggpubr::ggarrange(p, p.ILR, ncol = 2)
```

for all 15 factors
```{r}
# Generate the list of factors from Phylofactor.1 to Phylofactor.15
factorvector <- paste0('Phylofactor.', 1:15)  # This creates a vector of "Phylofactor.1" to "Phylofactor.15"
labelvec <- paste0('Phylofactor ', 1:15)  # Creates labels like "Phylofactor 1", "Phylofactor 2", etc.

# Plot using the ILR_plotfunction
p.ILR <- ILR_plotfunction(ps.ILR,
                          factorvector = factorvector,  # Factors from 1 to 15
                          explanatory = "snb_group",    # Metadata used in the model
                          labelvec = labelvec,          # Labels for each Phylofactor
                          ncol = 3) +                   # Number of columns for facetting
         scale_fill_manual(values = cols)               # Color scale for the fill
```

```{r}
# Create summary for factor 3
s <- summary(PF,taxonomy, factor=6) 

# extract ASV IDs for group 1
s.ids <- s$species.list$Group2

# subset ps object
my_subset <- subset(otu_table(ps.phylo), rownames(otu_table(ps.phylo)) %in% s.ids)
new_physeq <- merge_phyloseq(my_subset, tax_table(ps.phylo), sample_data(ps.phylo))

# Use custom function to get prevalance table on Phylum level
plotdata <- prevalencedf(new_physeq, phylum)
plotdata$Phylum <- factor(plotdata$phylum, levels = plotdata$phylum)

# create the order vectors, determining the Phyla names present in legend and their order
ordervec <- (plotdata %>% dplyr::filter(n > 2))$Phylum 

# Legend labels, has to be same length as ordervec
labelvec <- (plotdata %>% 
            dplyr::filter(n > 2) %>% 
            mutate(Phylum = paste0(Phylum, " (", n, ")")))$Phylum

# colors from the tree colours - to align with tree figure if necessary
# ploting HOW MANY GENERA are present in each phylum
p.pf3 <- plotdata %>% 
  ggpubr::ggpie("n", fill = "Phylum") +
  scale_fill_manual(values = treecols, 
                    breaks = ordervec,
                    labels = labelvec) +
  theme(axis.text.x=element_blank(),
        legend.position = "right") +
  guides(fill = guide_legend(ncol = 2, byrow = FALSE)) +
  labs(title = "Phylofactor 6 (Group 1)", subtitle = "Number of Genera") 

# Combine plots
ggpubr::ggarrange(ggarrange(p, p.ILR, nrow = 1),
                  p.pf3, nrow = 2)
```
