5.microbiome_characterization.Rmd

---
title: "microbiome_characterization"
author: "Erica Ryu"
date: "3/15/2022"
output: pdf_document
---

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

# 5. Microbiome Characterization

The purpose of this script is to examine standard microbiome characteristics (alpha and beta diversity)

## load packages
```{r}
library(phyloseq)
library(ggplot2)
library(ggsignif)
library(FSA)
library(vegan)
library(dplyr)
library(rbiom)
library(ape)
library(pairwiseAdonis)
```

## load data
```{r}
phyloseq_complete <- readRDS("output/ps_complete.rds")
rfxn_micro <- read.csv("output/rarefaction_edit03292024_final.csv", header = TRUE)

# subset based on extraction kit
qiagen <- subset_samples(phyloseq_complete, Condition == "Qiagen")
psoil <- subset_samples(phyloseq_complete, Condition == "Psoil")
```

## load functions
```{r}
# calculate beta diversity
beta_ordinate <- function(physeq, beta_dist){
  # set up data
  ps.prop <- transform_sample_counts(physeq, function(otu) otu/sum(otu))
  ps.prop@otu_table <- na.omit(ps.prop@otu_table)
  beta_inp <- t(as(otu_table(ps.prop), "matrix"))
  ## calculate distance and ordinate
  if (beta_dist == "unweighted unifrac") {
    dist <- unifrac(beta_inp, weighted = FALSE, tree = ps.prop@phy_tree)
    return(pcoa(dist))
  } else if (beta_dist == "weighted unifrac") {
    dist <- unifrac(beta_inp, weighted = TRUE, tree = ps.prop@phy_tree)
    return(pcoa(dist))
  } else {
    ord.pcoa <- ordinate(ps.prop, method = "PCoA", distance = "bray", weighted = FALSE)
  }
}

# set up beta diversity data for plotting
plot_beta <- function(beta, physeq){
  # extract axes
  PCOAaxes <- beta$vectors[,c(1,2,3,4)]
  # extract lifestyle column from metadata and add to vectors
  lifestyle <- physeq@sam_data$Lifestyle
  PCOAaxes_meta <- cbind(PCOAaxes, lifestyle)
  df_PCOA<- as.data.frame(as.matrix(PCOAaxes_meta))
  # change industrial to American industrial
  df_PCOA$lifestyle<- gsub("Industrial", "American Industrial", df_PCOA$lifestyle)
  # add space to Recently Settled
  df_PCOA$lifestyle <- gsub("RecentlySettled", "Recently Settled", df_PCOA$lifestyle)
  df_PCOA$lifestyle <- factor(df_PCOA$lifestyle, levels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"))
  return(df_PCOA)
}

# calculate permanova
permanova <- function(physeq, beta_dist){
  # set up data
  ps.prop <- transform_sample_counts(physeq, function(otu) otu/sum(otu))
  ps.prop@otu_table <- na.omit(ps.prop@otu_table)
  df <- data.frame(sample_data(ps.prop))
  beta_inp <- t(as(otu_table(ps.prop), "matrix"))
  if (beta_dist == "unweighted unifrac") {
    dist <- unifrac(beta_inp, weighted = FALSE, tree = ps.prop@phy_tree)
  } else if (beta_dist == "weighted unifrac") {
    dist <- unifrac(beta_inp, weighted = TRUE, tree = ps.prop@phy_tree)
  } else {
    dist = phyloseq::distance(ps.prop, method="bray", normalized=TRUE, parallel=FALSE, fast=TRUE) 
  }
   perm <- adonis2(dist ~ Lifestyle, data = df, permutations = 99999)
   print(perm)
}
 
```

## set colors
```{r}
fivecolors <- c("darkslateblue", "deepskyblue", "lightblue3", "lightsalmon" , "firebrick")
```

## examine alpha diversity
```{r}
# order by lifestyle
rfxn_micro$Lifestyle <- factor(rfxn_micro$Lifestyle, ordered = TRUE, levels=c("Foragers", "RecentlySettled", "Agriculturalists", "Expats", "Industrial", "control"))

qiagen_rfxn <- subset(rfxn_micro, Condition == "Qiagen")
psoil_rfxn <- subset(rfxn_micro, Condition == "Psoil")

# fix names
qiagen_rfxn$Lifestyle <- gsub("RecentlySettled", "Recently Settled", qiagen_rfxn$Lifestyle)
psoil_rfxn$Lifestyle <- gsub("RecentlySettled", "Recently Settled", psoil_rfxn$Lifestyle)

qiagen_rfxn$Lifestyle <- gsub("Industrial", "American Industrial", qiagen_rfxn$Lifestyle)
psoil_rfxn$Lifestyle <- gsub("Industrial", "American Industrial", psoil_rfxn$Lifestyle)

# filter by alpha diversity measure
shannon_qiagen <- subset(qiagen_rfxn, measure == "Shannon")
shannon_psoil <- subset(psoil_rfxn, measure == "Shannon")

faiths_qiagen <- subset(qiagen_rfxn, measure == "Faiths")
faiths_psoil <- subset(psoil_rfxn, measure == "Faiths")

rich_qiagen <- subset(qiagen_rfxn, measure == "Richness")
rich_psoil <- subset(psoil_rfxn, measure == "Richness")

simpson_qiagen <- subset(qiagen_rfxn, measure == "Simpson")
simpson_psoil <- subset(psoil_rfxn, measure == "Simpson")

fisher_qiagen <- subset(qiagen_rfxn, measure == "Fisher")
fisher_psoil <- subset(psoil_rfxn, measure == "Fisher")

# overall alpha diversity
qiagen_alpha <- ggplot(qiagen_rfxn, aes(x = Lifestyle, y = mean, group = Lifestyle), color = black) +
  geom_violin(aes(fill=Lifestyle), alpha=0.8) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) +
  geom_signif(comparisons = list(c("Expats", "American Industrial"), c("Agriculturalists", "American Industrial"), c("Recently Settled", "American Industrial"), c("American Industrial", "Foragers")), map_signif_level = TRUE, step_increase = 0.1) +
  scale_x_discrete(limits=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats",  "American Industrial")) +
  geom_jitter(size = 1, width=0.15, col="darkgreen") + facet_wrap(~measure, ncol = 3, scales = "free") +
  theme(axis.title.x = element_blank(),axis.title.y = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 0.95, vjust=0.95)) +
  theme(legend.position = "none") + xlab("")

ggsave(file = "figures/qiagen_alpha.pdf", width = 6, height = 8, plot = qiagen_alpha)

# subset to specific metrics for visualization
qiagen_alpha_plot <- subset(qiagen_rfxn, measure == "Shannon" | measure == "Faiths")

qiagen_alpha_plot_obj <- ggplot(qiagen_alpha_plot, aes(x = Lifestyle, y = mean, group = Lifestyle), color = black) +
  geom_violin(aes(fill=Lifestyle), alpha=0.8) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) +
  geom_signif(comparisons = list(c("Expats", "American Industrial"), c("Agriculturalists", "American Industrial"), c("Recently Settled", "American Industrial"), c("American Industrial", "Foragers")), map_signif_level = TRUE, step_increase = 0.1) +
  scale_x_discrete(limits=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats",  "American Industrial")) +
  geom_jitter(size = 1, width=0.15, col="darkgreen") + facet_wrap(~measure, ncol = 3, scales = "free") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 0.95, vjust=0.95)) +
  theme(legend.position = "none") + xlab("")

ggsave(file = "figures/qiagen_alpha_plot_obj.pdf", width = 6, height = 4, plot = qiagen_alpha_plot_obj)

psoil_alpha <- ggplot(psoil_rfxn, aes(x = Lifestyle, y = mean, group = Lifestyle), color = black) +
  geom_violin(aes(fill=Lifestyle), alpha=0.8) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) +
  scale_x_discrete(limits=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats",  "American Industrial")) +
  geom_jitter(size = 1, width=0.15, col="darkgreen") + facet_wrap(~measure, ncol = 3, scales = "free") +
  theme(axis.title.x = element_blank(),axis.title.y = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 0.95, vjust=0.95)) +
  theme(legend.position = "none") + xlab("")

ggsave(file = "figures/psoil_alpha.pdf", width = 6, height = 7, plot = psoil_alpha)

## calculate significance
# KW - qiagen
kruskal.test(mean ~ Lifestyle, data = shannon_qiagen)
kruskal.test(mean ~ Lifestyle, data = faiths_qiagen)
kruskal.test(mean ~ Lifestyle, data = rich_qiagen)
kruskal.test(mean ~ Lifestyle, data = simpson_qiagen)
kruskal.test(mean ~ Lifestyle, data = fisher_qiagen)

# KW - psoil
kruskal.test(mean ~ Lifestyle, data = shannon_psoil)
kruskal.test(mean ~ Lifestyle, data = faiths_psoil)
kruskal.test(mean ~ Lifestyle, data = rich_psoil)
kruskal.test(mean ~ Lifestyle, data = simpson_psoil)
kruskal.test(mean ~ Lifestyle, data = fisher_psoil)

## dunn's test
faiths_dunn <- dunnTest(mean ~ Lifestyle, data = faiths_qiagen, method="bh")
rich_dunn <- dunnTest(mean ~ Lifestyle, data = rich_qiagen, method="bh")
fisher_dunn <- dunnTest(mean ~ Lifestyle, data = fisher_qiagen, method="bh")

# extract as df
df_faiths_dunn <- faiths_dunn[["res"]]
df_rich_dunn <- rich_dunn[["res"]]
df_fisher_dunn <- fisher_dunn[["res"]]

# order by p value
df_faiths_dunn <- df_faiths_dunn[order(df_faiths_dunn$P.adj),]
df_rich_dunn <- df_rich_dunn[order(df_rich_dunn$P.adj),]
df_fisher_dunn <- df_fisher_dunn[order(df_fisher_dunn$P.adj),]

# extract sig p values
df_faiths_dunn_sig <- df_faiths_dunn[(df_faiths_dunn$P.adj < 0.05),]
df_rich_dunn_sig <- df_rich_dunn[(df_rich_dunn$P.adj < 0.05),]
df_fisher_dunn_sig <- df_fisher_dunn[(df_fisher_dunn$P.adj < 0.05),]

# add column for diversity metric
df_faiths_dunn_sig$metric <- rep(c("faiths"), times = length(df_faiths_dunn_sig$P.adj))
df_rich_dunn_sig$metric <- rep(c("rich"), times = length(df_rich_dunn_sig$P.adj))
df_fisher_dunn_sig$metric <- rep(c("fisher"), times = length(df_fisher_dunn_sig$P.adj))

# combine all sig
df_dunn_sig <- rbind(df_faiths_dunn_sig, df_rich_dunn_sig, df_fisher_dunn_sig)
df_dunn_sig

## alpha diversity for ONLY nepali individuals
faiths_qiagen_nepal <- filter(faiths_qiagen, Lifestyle != "American Industrial")
rich_qiagen_nepal <- filter(rich_qiagen, Lifestyle != "American Industrial")
fisher_qiagen_nepal <- filter(fisher_qiagen, Lifestyle != "American Industrial")
# KW
kruskal.test(mean ~ Lifestyle, data = faiths_qiagen_nepal)
kruskal.test(mean ~ Lifestyle, data = rich_qiagen_nepal)
kruskal.test(mean ~ Lifestyle, data = fisher_qiagen_nepal)
```

## examine beta diversity
```{r}
bray_ordinate <- beta_ordinate(qiagen, beta_dist = "bray")
bray_plot <- plot_beta(bray_ordinate, qiagen)

bray <- ggplot(bray_plot, aes(x = as.numeric(Axis.1), y = as.numeric(Axis.2))) +
  geom_point(shape = 21, color = "black", size = 3, aes(fill = lifestyle)) +
  scale_x_continuous(breaks=seq(-1,1,0.1)) +
  scale_y_continuous(breaks=seq(-1,1,0.1)) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"))+
  labs(x = "PCoA1 (28.62%)",
       y = "PCoA2 (14.99%)")

ggsave(file = "figures/bray_small.pdf", plot = bray, width = 6, height = 3)

# permanova
set.seed(105)
bray_perm <- permanova(qiagen, beta_dist = "bray")

# pairwise adonis
ps_RA <- transform_sample_counts(qiagen, function(otu) otu/sum(otu))
ps_RA@otu_table <- na.omit(ps_RA@otu_table)
dist_bray <- phyloseq::distance(ps_RA, method="bray", normalized=TRUE, parallel=FALSE, fast=TRUE) 
pw_adonis <-  pairwise.adonis(dist_bray, phyloseq::sample_data(ps_RA)$Lifestyle)

pairwise_adonis(x, factors, sim_method = "bray",
  p_adjust_m = "bonferroni", reduce = NULL)

## PCoA1
bray_pcoa1 <- ggplot(bray_plot, aes(x = lifestyle, y = as.numeric(Axis.1), group = lifestyle), color = black) +
  geom_violin(alpha = 0.8, aes(fill=lifestyle)) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) + 
  scale_y_continuous(breaks=seq(-1,1,0.1), labels = scales::comma) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  geom_jitter(size = 1, col="darkgreen", width = 0.1) +
  labs(x = "Lifestyle", y = "PCoA1 (28.62%)") + 
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5, vjust=0.5)) + 
  theme(legend.position = "none")

# load package for trend test
library(DescTools) # this MUST be loaded after running permanova and all other vegan functions, because loading both DescTools and vegan causes vegan to stop working

# trend test for axis 1
JonckheereTerpstraTest(as.numeric(Axis.1) ~ lifestyle, data = bray_plot)

## PCoA2
bray_pcoa2 <- ggplot(bray_plot, aes(x = lifestyle, y = as.numeric(Axis.2), group = lifestyle), color = black) +
  geom_violin(alpha = 0.8, aes(fill=lifestyle)) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) +
  scale_y_continuous(breaks=seq(-1,1,0.1), labels = scales::comma) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  geom_jitter(size = 1, col="darkgreen", width = 0.1) +
  labs(x = "Lifestyle", y = "PCoA2 (14.99%)") +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5, vjust=0.5)) +
  theme(legend.position = "none")

# trend test for axis 2
JonckheereTerpstraTest(as.numeric(Axis.2) ~ lifestyle, data = bray_plot)
ggsave(file = "figures/bray_pcoa1_small.pdf", plot = bray_pcoa1, width = 6, height = 3)
ggsave(file = "figures/bray_pcoa2_small.pdf", plot = bray_pcoa2, width = 6, height = 3)

## all powersoil metrics
bray_ordinate_psoil <- beta_ordinate(psoil, beta_dist = "bray")
bray_plot_psoil <- plot_beta(bray_ordinate_psoil, psoil)

bray_psoil <- ggplot(bray_plot_psoil, aes(x = as.numeric(Axis.1), y = as.numeric(Axis.2))) +
  geom_point(shape = 21, color = "black", size = 3, aes(fill = lifestyle)) +
  scale_x_continuous(breaks=seq(-1,1,0.1)) +
  scale_y_continuous(breaks=seq(-1,1,0.1)) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"))+
  labs(x = "PCoA1 (25.99%)",
       y = "PCoA2 (13.83%)")

ggsave(file = "figures/bray_psoil_small.pdf", plot = bray_psoil, width = 6, height = 3)

bray_perm_psoil <- permanova(psoil, beta_dist = "bray")

bray_pcoa1_psoil <- ggplot(bray_plot_psoil, aes(x = lifestyle, y = as.numeric(Axis.1)*-1, group = lifestyle), color = black) +
  geom_violin(alpha = 0.8, aes(fill=lifestyle)) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) + 
  scale_y_continuous(breaks=seq(-1,1,0.1), labels = scales::comma) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  geom_jitter(size = 1, col="darkgreen", width = 0.1) +
  labs(x = "Lifestyle", y = "PCoA1 (25.99%)") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust=0.5)) + 
  theme(legend.position = "none")

ggsave(file = "figures/bray_pcoa1_psoil_small.pdf", plot = bray_pcoa1_psoil, width = 4, height = 3)

# trend test for axis 1
JonckheereTerpstraTest(as.numeric(Axis.1) ~ lifestyle, data = bray_plot_psoil)

bray_pcoa2_psoil <- ggplot(bray_plot_psoil, aes(x = lifestyle, y = as.numeric(Axis.2)*-1, group = lifestyle), color = black) +
  geom_violin(alpha = 0.8, aes(fill=lifestyle)) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) + 
  scale_y_continuous(breaks=seq(-1,1,0.1), labels = scales::comma) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  geom_jitter(size = 1, col="darkgreen", width = 0.1) +
  labs(x = "Lifestyle", y = "PCoA2 (13.83%)") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust=0.5)) + 
  theme(legend.position = "none")

ggsave(file = "figures/bray_pcoa2_psoil_small.pdf", plot = bray_pcoa2_psoil, width = 4, height = 3)

# trend test for axis 2
JonckheereTerpstraTest(as.numeric(Axis.2) ~ lifestyle, data = bray_plot_psoil)
```

## check unweighted unifrac
```{r}
bray_ordinate_uw <- beta_ordinate(qiagen, beta_dist = "unweighted unifrac")
bray_plot_uw <- plot_beta(bray_ordinate_uw, qiagen)

bray_uw <- ggplot(bray_plot_uw, aes(x = as.numeric(Axis.1), y = as.numeric(Axis.2))) +
  geom_point(shape = 21, color = "black", size = 3, aes(fill = lifestyle)) +
  scale_x_continuous(breaks=seq(-1,1,0.1)) +
  scale_y_continuous(breaks=seq(-1,1,0.1)) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"))+
  labs(x = "PCoA1 (21.90%)",
       y = "PCoA2 (11.77%)")

ggsave(file = "figures/bray_uw_small.pdf", plot = bray_uw, width = 6, height = 3)

# permanova
set.seed(105)
bray_perm_uw <- permanova(qiagen, beta_dist = "unweighted unifrac")

## PCoA1
bray_uw_pcoa1 <- ggplot(bray_plot_uw, aes(x = lifestyle, y = as.numeric(Axis.1), group = lifestyle), color = black) +
  geom_violin(alpha = 0.8, aes(fill=lifestyle)) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) + 
  scale_y_continuous(breaks=seq(-1,1,0.1), labels = scales::comma) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  geom_jitter(size = 1, col="darkgreen", width = 0.1) +
  labs(x = "Lifestyle", y = "PCoA1 (21.90%)") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust=0.5)) + 
  theme(legend.position = "none")

# load package for trend test
library(DescTools) # this MUST be loaded after running permanova and all other vegan functions, because loading both DescTools and vegan causes vegan to stop working

# trend test for axis 1
JonckheereTerpstraTest(as.numeric(Axis.1) ~ lifestyle, data = bray_plot_uw)

ggsave(file = "figures/bray_uw_pcoa1.pdf", plot = bray_uw_pcoa1, width = 4, height = 3)

## PCoA2
bray_uw_pcoa2 <- ggplot(bray_plot_uw, aes(x = lifestyle, y = as.numeric(Axis.2), group = lifestyle), color = black) +
  geom_violin(alpha = 0.8, aes(fill=lifestyle)) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) +
  scale_y_continuous(breaks=seq(-1,1,0.1), labels = scales::comma) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  geom_jitter(size = 1, col="darkgreen", width = 0.1) +
  labs(x = "Lifestyle", y = "PCoA2 (11.77%)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust=0.5)) +
  theme(legend.position = "none")

# trend test for axis 2
JonckheereTerpstraTest(as.numeric(Axis.2) ~ lifestyle, data = bray_plot_uw)

ggsave(file = "figures/bray_uw_pcoa2.pdf", plot = bray_uw_pcoa2, width = 4, height = 3)

```

## check weighted unifrac
```{r}
bray_ordinate_w <- beta_ordinate(qiagen, beta_dist = "weighted unifrac")
bray_plot_w <- plot_beta(bray_ordinate_w, qiagen)

bray_w <- ggplot(bray_plot_w, aes(x = as.numeric(Axis.1), y = as.numeric(Axis.2))) +
  geom_point(shape = 21, color = "black", size = 3, aes(fill = lifestyle)) +
  scale_x_continuous(breaks=seq(-1,1,0.2)) +
  scale_y_continuous(breaks=seq(-1,1,0.2)) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"))+
  labs(x = "PCoA1 (53.61%)",
       y = "PCoA2 (21.16%)")

ggsave(file = "figures/bray_w_small.pdf", plot = bray_w, width = 6, height = 3)

# permanova
set.seed(105)
bray_perm_w <- permanova(qiagen, beta_dist = "weighted unifrac")

## PCoA1
bray_w_pcoa1 <- ggplot(bray_plot_w, aes(x = lifestyle, y = as.numeric(Axis.1), group = lifestyle), color = black) +
  geom_violin(alpha = 0.8, aes(fill=lifestyle)) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) + 
  scale_y_continuous(breaks=seq(-1,1,0.2), labels = scales::comma) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  geom_jitter(size = 1, col="darkgreen", width = 0.1) +
  labs(x = "Lifestyle", y = "PCoA1 (53.61%)") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust=0.5)) + 
  theme(legend.position = "none")

# load package for trend test
library(DescTools) # this MUST be loaded after running permanova and all other vegan functions, because loading both DescTools and vegan causes vegan to stop working

# trend test for axis 1
JonckheereTerpstraTest(as.numeric(Axis.1) ~ lifestyle, data = bray_plot_w)

ggsave(file = "figures/bray_w_pcoa1.pdf", plot = bray_w_pcoa1, width = 4, height = 3)

## PCoA2
bray_w_pcoa2 <- ggplot(bray_plot_w, aes(x = lifestyle, y = as.numeric(Axis.2), group = lifestyle), color = black) +
  geom_violin(alpha = 0.8, aes(fill=lifestyle)) +
  geom_boxplot(width=0.1, color="black", alpha=0.2, outlier.shape = NA) +
  scale_y_continuous(breaks=seq(-1,1,0.2), labels = scales::comma) +
  scale_fill_manual(name=NULL,
                    values=fivecolors,
                    breaks=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial"),
                    labels=c("Foragers", "Recently Settled", "Agriculturalists", "Expats", "American Industrial")) +
  geom_jitter(size = 1, col="darkgreen", width = 0.1) +
  labs(x = "Lifestyle", y = "PCoA2 (21.16%)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust=0.5)) +
  theme(legend.position = "none")

# trend test for axis 2
JonckheereTerpstraTest(as.numeric(Axis.2) ~ lifestyle, data = bray_plot_w)

ggsave(file = "figures/bray_w_pcoa2.pdf", plot = bray_w_pcoa2, width = 4, height = 3)
```

## comparison with covariates
```{r}
# sex
ps.prop_cov <- transform_sample_counts(qiagen, function(otu) otu/sum(otu))
ps.prop_cov@otu_table <- na.omit(ps.prop_cov@otu_table)
df_cov <- data.frame(sample_data(ps.prop_cov))
beta_inp_cov <- t(as(otu_table(ps.prop_cov), "matrix"))

set.seed(111)
dist_bray <- phyloseq::distance(ps.prop_cov, method="bray", normalized=TRUE, parallel=FALSE, fast=TRUE) 
perm_bray_sex <- adonis2(dist_bray ~ Lifestyle + SEX, data = df_cov, permutations = 99999)

# sex and age
no_am_ps <- subset_samples(qiagen, Lifestyle != "Industrial")
ps.prop_noam <- transform_sample_counts(no_am_ps, function(otu) otu/sum(otu))
ps.prop_noam@otu_table <- na.omit(ps.prop_noam@otu_table)
df_noam <- data.frame(sample_data(ps.prop_noam))
beta_inp_noam <- t(as(otu_table(ps.prop_noam), "matrix"))

set.seed(100)
dist_bray_noam <- phyloseq::distance(ps.prop_noam, method="bray", normalized=TRUE, parallel=FALSE, fast=TRUE) 
perm_bray_noam <- adonis2(dist_bray_noam ~ Lifestyle + SEX + AGE, data = df_noam, permutations = 99999)
```