Chapter 14 Supp table 2

Beta diversity plots of samples using the preprocessed rarefied abundance values.

14.1 Libraries

library("phyloseq")
library("microbiome")
#Will use Unifrac() & adonis() function from vegan package
library(vegan)
library(tidyverse)

14.1.1 Statistical test

#Carry out PERMANOVA test
#Extract metadata
meta <- microbiome::meta(physeq)
#Create new temp column
meta$Storagetempcategory <- meta$Storagetemp
meta$Storagetempcategory <- gsub(pattern = "-80", replacement = "Baseline", x = meta$Storagetempcategory)
meta$Storagetempcategory <- gsub(pattern = "4", replacement = "Fridge", x = meta$Storagetempcategory)
meta$Storagetempcategory <- gsub(pattern = "20", replacement = "RT", x = meta$Storagetempcategory)
meta$Storageconditions <- paste0(meta$Storagetempcategory,"_",meta$DaysstoredpriortoDNAextraction,"D")
#adonis doesn't have weighted unifrac distance so need to calculate it
unifrac.dist <- as.matrix(UniFrac(physeq, weighted = TRUE,normalized = TRUE,  
                        parallel = FALSE,fast = TRUE))

#Pairwise comparison function
paiwise_function <- function(dist_df, meta_df, meta_column) {
  #Combination
  cbn <- combn(x = unique(meta_df[,meta_column]), m = 2)
  #Create empty final data frame with 4 columns
  # and a number of rows equal to the the number of combinations
  pairwise_permanova_df <- as.data.frame(matrix(data = NA, nrow = ncol(cbn), ncol = 4))
  #Add column names
  colnames(pairwise_permanova_df) <- c("X","Y","p","p.adj")
  #Loop through the combinations
  for(i in 1:ncol(cbn)){
    #Subset metadf to contain only our samples of interest
    metadf_subset <- meta[meta[,meta_column]%in% cbn[,i],]
    #Subset distance matrix
    samples_subset <- row.names(metadf_subset)
    dist.subset <- dist_df[samples_subset,samples_subset]
    #PERMANOVA/ADONIS
    f <- paste0("dist.subset~",meta_column)
    pairwise_adonis <- vegan::adonis2(formula = formula(f), permutations = 999,
                                               data = metadf_subset, by = "margin")
    #Add the group names and p-value to the main data frame
    pairwise_permanova_df[i,1:2] <- cbn[,i]
    pairwise_permanova_df[i,3] <- pairwise_adonis[1,"Pr(>F)"]
  }
  #Add adjusted P-values
  pairwise_permanova_df$p.adj <- p.adjust(pairwise_permanova_df$p, method = "BH")
  #Convert to wide table with adjusted p-values
  p_adjust_long <- pairwise_permanova_df[,c(1,2,4)]
  p_adjust_wide <- spread(data = p_adjust_long, Y, p.adj, fill=NA)
  #Reformat
  row.names(p_adjust_wide) <- p_adjust_wide[,1]
  p_adjust_wide <- p_adjust_wide[,-1]
  return(p_adjust_wide)
}

#PERMANOVA tests
#Patient
patient_permanova <- 
  vegan::adonis2(unifrac.dist ~ Patientnumber,
               data = meta, permutations=999, by = "margin")
#Change a row name
row.names(patient_permanova)[1] <- "Participant number"
#pairwise
patient_pairwise <- paiwise_function(unifrac.dist, meta, "Patientnumber")
#Write out data
write.table("Participant number", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, col.names = FALSE, row.names = FALSE)
write.table(patient_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table(patient_pairwise, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(patient_permanova, patient_pairwise)

#Buffer
buffer_permanova <- 
  vegan::adonis2(unifrac.dist ~ Bufferused ,
               data = meta, permutations=999, by = "margin")
#pairwise
buffer_pairwise <- paiwise_function(unifrac.dist, meta, "Bufferused")
#Write out data
write.table("Bufferused", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append=TRUE,col.names = FALSE, row.names = FALSE)
write.table(buffer_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table(buffer_pairwise, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(buffer_permanova, buffer_pairwise)

#temp
temp_permanova <- 
  vegan::adonis2(unifrac.dist ~ Storagetempcategory,
               data = meta, permutations=999, by = "margin")
#Change a row name
row.names(temp_permanova)[1] <- "Temp"
#pairwise
temp_pairwise <- paiwise_function(unifrac.dist, meta, "Storagetempcategory")
#Write out data
write.table("Temp", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append=TRUE,col.names = FALSE, row.names = FALSE)
write.table(temp_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table(temp_pairwise, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(temp_permanova, temp_pairwise)

#Storage condition
condition_permanova <- 
  vegan::adonis2(unifrac.dist ~ Storageconditions,
               data = meta, permutations=999, by = "margin")
#pairwise
condition_pairwise <- paiwise_function(unifrac.dist, meta, "Storageconditions")
#Write out data
write.table("Storageconditions", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append=TRUE,col.names = FALSE, row.names = FALSE)
write.table(condition_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table(condition_pairwise, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(condition_permanova, condition_pairwise)

#Buffer and participant
buffer_patient_permanova <- 
  vegan::adonis2(unifrac.dist ~ Bufferused + Patientnumber,
               data = meta, permutations=999, by = "margin")
#Change row names
row.names(buffer_patient_permanova)[1:2] <- c("Bufferused","Participant number")
#pairwise
#create new metadata column
meta$Bufferused_and_Patientnumber <- paste0(meta$Bufferused,"_",meta$Patientnumber)
buffer_patient_pairwise <- paiwise_function(unifrac.dist, meta, "Bufferused_and_Patientnumber")
#Write out data
write.table("Bufferused and Participant number", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append=TRUE,col.names = FALSE, row.names = FALSE)
write.table(buffer_patient_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table(buffer_patient_pairwise, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(buffer_patient_permanova, buffer_patient_pairwise)

#Buffer and temp
buffer_temp_permanova <- 
  vegan::adonis2(unifrac.dist ~ Bufferused + Storagetempcategory,
               data = meta, permutations=999, by = "margin")
#Change row names
row.names(buffer_temp_permanova)[1:2] <- c("Bufferused","Temp")
#pairwise
#create new metadata column
meta$Bufferused_and_Temp <- paste0(meta$Bufferused,"_",meta$Storagetempcategory)
buffer_temp_pairwise <- paiwise_function(unifrac.dist, meta, "Bufferused_and_Temp")
#Write out data
write.table("Bufferused and Temp", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append=TRUE,col.names = FALSE, row.names = FALSE)
write.table(buffer_temp_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table(buffer_temp_pairwise, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(buffer_temp_permanova, buffer_temp_pairwise)

#Temp and Participant
temp_participant_permanova <- 
  vegan::adonis2(unifrac.dist ~ Storagetempcategory + Patientnumber,
               data = meta, permutations=999, by = "margin")
#Change row names
row.names(temp_participant_permanova)[1:2] <- c("Temp","Participant number")
#pairwise
#create new metadata column
meta$Temp_and_participant <- paste0(meta$Storagetempcategory,"_",meta$Patientnumber)
temp_participant_pairwise <- paiwise_function(unifrac.dist, meta, "Temp_and_participant")
#Write out data
write.table("Temp and Participant number", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append=TRUE,col.names = FALSE, row.names = FALSE)
write.table(temp_participant_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table(temp_participant_pairwise, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(temp_participant_permanova, temp_participant_pairwise)

#Storage condition and participant
condition_patient_permanova <- 
  vegan::adonis2(unifrac.dist ~ Storageconditions + Patientnumber,
               data = meta, permutations=999, by = "margin")
#Change row names
row.names(condition_patient_permanova)[1:2] <- c("Storageconditions","Participant number")
#Write out data
write.table("Storageconditions and Participant number", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append=TRUE,col.names = FALSE, row.names = FALSE)
write.table(condition_patient_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise PERMANOVA skipped as there would be too few samples per grouping", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(condition_patient_permanova)

#Buffer, Storage condition and participant
buffer_condition_patient_permanova <- 
  vegan::adonis2(unifrac.dist ~ Bufferused + Storageconditions + Patientnumber,
               data = meta, permutations=999, by = "margin")
#Change row names
row.names(buffer_condition_patient_permanova)[1:3] <- 
  c("Bufferused", "Storageconditions","Participant number")
#Write out data
write.table("Bufferused, Storageconditions and Participant number", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append=TRUE,col.names = FALSE, row.names = FALSE)
write.table(buffer_condition_patient_permanova, file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = NA)
write.table("Pairwise PERMANOVA skipped as there would be too few samples per grouping", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
write.table("", file = "./data/permanova_16s.tsv", 
            sep = "\t", quote=FALSE, append = TRUE, col.names = FALSE, row.names = FALSE)
#Remove results
remove(buffer_condition_patient_permanova)




#buffer stats
buffer_permanova <- 
  adonis(unifrac.dist ~ Bufferused ,
               data = meta, permutations=999, method="euclidean")
write.table(as.data.frame(buffer_permanova$aov.tab),
            file = "./data/16s_buffer_permanova.tsv",
            sep = "\t", quote=FALSE)
#supp fig 2A stats
buffer_patient_permanova <- 
  adonis(unifrac.dist ~ Bufferused + Patientnumber,
               data = meta, permutations=999, method="euclidean")
write.table(as.data.frame(buffer_patient_permanova$aov.tab),
            file = "./data/16s_buffer_patient_permanova.tsv",
            sep = "\t", quote=FALSE)
#supp fig 2b stats
temp_patient_permanova <- 
  adonis(unifrac.dist ~ Storagetemp + Patientnumber,
               data = meta, permutations=999, method="euclidean")
write.table(as.data.frame(temp_patient_permanova$aov.tab),
            file = "./data/16s_temp_patient_permanova.tsv",
            sep = "\t", quote=FALSE)
#Storage condition
condition_patient_permanova <- 
  adonis(unifrac.dist ~ Storageconditions,
               data = meta, permutations=999, method="euclidean")
write.table(as.data.frame(condition_patient_permanova$aov.tab),
            file = "./data/16s_condition_permanova.tsv",
            sep = "\t", quote=FALSE)
#supp fig 2c stats
condition_patient_permanova <- 
  adonis(unifrac.dist ~ Storageconditions + Patientnumber,
               data = meta, permutations=999, method="euclidean")
write.table(as.data.frame(condition_patient_permanova$aov.tab),
            file = "./data/16s_condition_patient_permanova.tsv",
            sep = "\t", quote=FALSE)
#supp fig 2d stats
buffer_condition_patient_permanova <- 
  adonis(unifrac.dist ~ Bufferused + Storageconditions + Patientnumber,
               data = meta, permutations=999, method="euclidean")
write.table(as.data.frame(buffer_condition_patient_permanova$aov.tab),
            file = "./data/16s_buffer_condition_patient_permanova.tsv",
            sep = "\t", quote=FALSE)

#Pairwise permanovas
#Buffers
buffer_pair_permanova <- 
  pairwise.perm.manova(resp=unifrac.dist,fact=meta$Bufferused,
                                              nperm = 999)
write.table(as.data.frame(buffer_pair_permanova$p.value),
            file = "./data/16s_buffer_pair_permanova.tsv",
            sep = "\t", quote=FALSE)
#Buffer & patients
meta$Bufferused_and_Patientnumber <- paste0(meta$Bufferused,"_",meta$Patientnumber)
buffer_patient_pair_permanova <- 
  pairwise.perm.manova(resp=unifrac.dist,fact=meta$Bufferused_and_Patientnumber,
                                              nperm = 999)
write.table(as.data.frame(buffer_patient_pair_permanova$p.value),
            file = "./data/16s_buffer_patient_pair_permanova.tsv",
            sep = "\t", quote=FALSE)
#Temp & patients
meta$Storagetemp_and_Patientnumber <- paste0(meta$Storagetempcategory,"_",meta$Patientnumber)
temp_patient_pair_permanova <- 
  pairwise.perm.manova(resp=unifrac.dist,fact=meta$Storagetemp_and_Patientnumber,
                                              nperm = 999)
write.table(as.data.frame(temp_patient_pair_permanova$p.value),
            file = "./data/16s_temp_patient_pair_permanova.tsv",
            sep = "\t", quote=FALSE)
#Condition
cond_pair_permanova <- 
  pairwise.perm.manova(resp=unifrac.dist,fact=meta$Storageconditions,
                                              nperm = 999)
write.table(as.data.frame(cond_pair_permanova$p.value),
            file = "./data/16s_cond_pair_permanova.tsv",
            sep = "\t", quote=FALSE)