---

title: "5. Alpha Diversity Estimates"

description: |

Reproducible workflow to assess alpha diversity across temperature treatments using Hill numbers.

format:

html:

toc: true

toc-depth: 3

---

```{r}

#| message: false

#| results: hide

#| code-fold: true

#| code-summary: "Click here for setup information"

knitr::opts_chunk$set(echo = TRUE, eval = FALSE)

set.seed(119)

#library(conflicted)

library(phyloseq); packageVersion("phyloseq")

library(Biostrings); packageVersion("Biostrings")

pacman::p_load(tidyverse, hilldiv,

microbiome, phytools, phangorn,

pairwiseAdonis, naniar,

labdsv, patchwork, agricolae, ggpubr,

reactable, downloadthis, captioner,

install = FALSE, update = FALSE)

source("hack_code/div_test_plot_jjs.R")

source("hack_code/div_test_jjs.R")

options(scipen=999)

knitr::opts_current$get(c(

"cache",

"cache.path",

"cache.rebuild",

"dependson",

"autodep"

))

```

```{r}

#| include: false

#| eval: true

## Load to build page only #2

remove(list = ls())

load("page_build/alpha_wf.rdata")

```

```{r}

#| include: false

#| eval: true

# Create the caption(s) with captioner

caption_tab_ssu <- captioner(prefix = "(16S rRNA) Table", suffix = " |", style = "b")

caption_fig_ssu <- captioner(prefix = "(16S rRNA) Figure", suffix = " |", style = "b")

caption_tab_its <- captioner(prefix = "(ITS) Table", suffix = " |", style = "b")

caption_fig_its <- captioner(prefix = "(ITS) Figure", suffix = " |", style = "b")

# Create a function for referring to the tables in text

ref <- function(x) str_extract(x, "[^|]*") %>%

trimws(which = "right", whitespace = "[ ]")

```

```{r}

#| include: false

#| eval: true

source(file.path("assets", "functions.R"))

source("assets/captions/captions_alpha.R")

```

# Synopsis

This workflow contains diversity assessments for the 2018 high temperature data sets. In order to run the workflow, you either need to first run the [DADA2 Workflow](dada2.html) **and** the [Data Preparation workflow](data-prep.html) **or** download the files linked below. See the [Data Availability](data-availability.html) page for complete details.

## Workflow Input

Files needed to run this workflow can be downloaded from figshare.

<iframe src="https://widgets.figshare.com/articles/14690739/embed?show_title=1" width="100%" height="251" allowfullscreen frameborder="0"></iframe>

<iframe src="https://widgets.figshare.com/articles/14701440/embed?show_title=1" width="100%" height="351" allowfullscreen frameborder="0"></iframe>

# 16s rRNA

```{r}

#| include: false

## Initial Load for ANALYSIS #1

set.seed(119)

## ASV FULL AND PIME AND PERfect

ssu18_ps_work <- readRDS("files/data-prep/rdata/ssu18_ps_work.rds")

ssu18_ps_filt <- readRDS("files/filtering/arbitrary/rdata/ssu18_ps_filt.rds")

ssu18_ps_perfect <- readRDS("files/filtering/perfect/rdata/ssu18_ps_perfect.rds")

ssu18_ps_pime <- readRDS("files/filtering/pime/rdata/ssu18_ps_pime.rds")

```

```{r}

#| echo: false

swel_col <- c("#2271B2", "#71B222", "#B22271")

```

To account for presence of rare sequence variants caused by sequencing errors or other technical artifacts, we use Hill numbers [@alberdi2019guide]. Hill numbers allow the weight put on rare versus abundant sequence variants to be scaled while providing intuitive comparisons of diversity levels using “effective number of ASVs” as a measuring unit. This approach allows for balancing the over representation of rare ASVs that might be inflated due to sequencing errors.

We will then use Shapiro-Wilk tests to test for normalcy and then, depending on the results, either use parametric ANOVA or non-parametric Kruskal-Wallis to compare alpha diversity among treatments.

## Calculate Hill Numbers

To calculate Hill numbers, we use the R package `hilldiv` [@alberdi2019hilldiv]. We calculate three metrics that put more or less weight on common species:

1) Observed richness, where `q-value = 0`.

2) Shannon exponential, which weighs ASVs by their frequency, where `q-value = 1`.

3) Simpson multiplicative inverse, which over weighs abundant ASVs, where `q-value = 2`.

We perform each analysis against the Full (unfiltered) data set as well as the Arbitrary, PERfect and PIME, filtered data sets using the function `hill_div`.

The command is as follows:

```

hill_div(count = x, qvalue = i, tree = ultrametric_tree)

```

where `x` is the sample by ASV table, `i` is the q-value corresponding to the metric of interest and `tree` is an ultrametric formatted phylogenetic tree if you choose to look at lineage, rather than ASV, diversity.

```{r}

#| echo: false

## this is for subsettig samples..

rm(ssu18_asv_pime_sum)

ssu18_alpha_ds <- c("ssu18_ps_work", "ssu18_ps_filt", "ssu18_ps_perfect", "ssu18_ps_pime")

rm(list = ls(pattern = "_0"))

rm(list = ls(pattern = "_3"))

rm(list = ls(pattern = "_8"))

for (i in ssu18_alpha_ds) {

tmp_get <- get(i)

tmp_0 <- prune_samples(c("P10_D00_010_C0E", "P02_D00_010_C0A", "P04_D00_010_C0B", "P06_D00_010_C0C", "P08_D00_010_C0D"), tmp_get)

tmp_3 <- prune_samples(c("P01_D00_010_W3A", "P03_D00_010_W3B", "P05_D00_010_W3C", "P07_D00_010_W3D", "P09_D00_010_W3E"), tmp_get)

tmp_8 <- prune_samples(c("P01_D00_010_W8A", "P03_D00_010_W8B", "P05_D00_010_W8C", "P07_D00_010_W8D", "P09_D00_010_W8E"), tmp_get)

tmp_0 <- prune_taxa(taxa_sums(tmp_0) > 0, tmp_0)

tmp_3 <- prune_taxa(taxa_sums(tmp_3) > 0, tmp_3)

tmp_8 <- prune_taxa(taxa_sums(tmp_8) > 0, tmp_8)

tmp_ps_name_0 <- purrr::map_chr(i, ~ paste0(., "_0"))

assign(tmp_ps_name_0, tmp_0)

tmp_ps_name_3 <- purrr::map_chr(i, ~ paste0(., "_3"))

assign(tmp_ps_name_3, tmp_3)

tmp_ps_name_8 <- purrr::map_chr(i, ~ paste0(., "_8"))

assign(tmp_ps_name_8, tmp_8)

rm(list = ls(pattern = "tmp_"))

}

tmp_objects <- data.frame(c("FULL data set (0C)",

"FULL data set (3C)",

"FULL data set (8C)",

"Arbitrary filter data (0C)",

"Arbitrary filter data (3C)",

"Arbitrary filter data (8C)",

"PERfect filtered data (0C)",

"PERfect filtered data (3C)",

"PERfect filtered data (8C)",

"PIME (0C)",

"PIME (3C)",

"PIME (8C)"))

tmp_samples <- c("ssu18_ps_work_0", "ssu18_ps_work_3", "ssu18_ps_work_8",

"ssu18_ps_filt_0", "ssu18_ps_filt_3", "ssu18_ps_filt_8",

"ssu18_ps_perfect_0", "ssu18_ps_perfect_3", "ssu18_ps_perfect_8",

"ssu18_ps_pime_0", "ssu18_ps_pime_3", "ssu18_ps_pime_8")

tmp_no_samp <- c()

for (i in tmp_samples) {

tmp_get <- nsamples(get(i))

tmp_no_samp <- c(append(tmp_no_samp, tmp_get))

}

tmp_no_samp <- data.frame(tmp_no_samp)

tmp_rc <- c()

for (i in tmp_samples) {

tmp_get <- sum(readcount(get(i)))

tmp_rc <- c(append(tmp_rc, tmp_get))

}

tmp_rc <- data.frame(tmp_rc)

tmp_asv <- c()

for (i in tmp_samples) {

tmp_get <- ntaxa(get(i))

tmp_asv <- c(append(tmp_asv, tmp_get))

}

tmp_asv <- data.frame(tmp_asv)

ssu18_asv_pime_sum <- dplyr::bind_cols(tmp_objects, tmp_samples) %>%

dplyr::bind_cols(., tmp_no_samp) %>%

dplyr::bind_cols(., tmp_rc) %>%

dplyr::bind_cols(., tmp_asv) %>%

dplyr::rename("Description" = 1, "object name" = 2, "no. samples" = 3,

"total reads" = 4, "total asvs" = 5)

rm(list = ls(pattern = "tmp_"))

write.table(ssu18_asv_pime_sum, "ssu_split_by_temp.txt", sep = "\t", quote = FALSE, row.names = FALSE)

```

```{r}

#| echo: false

ssu18_alpha_ds <- c("ssu18_ps_work", "ssu18_ps_filt", "ssu18_ps_perfect", "ssu18_ps_pime")

tmp_objects <- data.frame(c("FULL data set",

"Arbitrary filter",

"PERfect filtered",

"PIME filtered"))

tmp_samples <- c("ssu18_ps_work",

"ssu18_ps_filt",

"ssu18_ps_perfect",

"ssu18_ps_pime")

tmp_no_samp <- c()

for (i in tmp_samples) {

tmp_get <- nsamples(get(i))

tmp_no_samp <- c(append(tmp_no_samp, tmp_get))

}

tmp_no_samp <- data.frame(tmp_no_samp)

tmp_rc <- c()

for (i in tmp_samples) {

tmp_get <- sum(readcount(get(i)))

tmp_rc <- c(append(tmp_rc, tmp_get))

}

tmp_rc <- data.frame(tmp_rc)

tmp_asv <- c()

for (i in tmp_samples) {

tmp_get <- ntaxa(get(i))

tmp_asv <- c(append(tmp_asv, tmp_get))

}

tmp_asv <- data.frame(tmp_asv)

ssu18_asv_pime_sum <- dplyr::bind_cols(tmp_objects, tmp_samples) %>%

dplyr::bind_cols(., tmp_no_samp) %>%

dplyr::bind_cols(., tmp_rc) %>%

dplyr::bind_cols(., tmp_asv) %>%

dplyr::rename("Description" = 1, "object name" = 2, "no. samples" = 3,

"total reads" = 4, "total asvs" = 5)

rm(list = ls(pattern = "tmp_"))

write.table(ssu18_asv_pime_sum, "ssu_filter_summmary.txt", sep = "\t", quote = FALSE, row.names = FALSE)

```

```{r}

#| echo: false

mergedGP <- merge_samples(ssu18_ps_perfect_3, "TEMP")

SD <- merge_samples(sample_data(ssu18_ps_perfect_3), "TEMP")

mergedGP

```

```{r}

#| echo: false

objects()

ssu18_ps_work@phy_tree <- NULL

ssu18_ps_filt@phy_tree <- NULL

ssu18_ps_pime@phy_tree <- NULL

ssu18_ps_perfect@phy_tree <- NULL

```

We first transform all the data to relative abundance values, and compute new trees.

```{r}

#| code-fold: true

ssu18_alpha_ds <- c("ssu18_ps_work", "ssu18_ps_filt", "ssu18_ps_perfect", "ssu18_ps_pime")

for (i in ssu18_alpha_ds) {

tmp_ps <- transform_sample_counts(get(i), function(otu) otu/sum(otu))

tmp_ps@phy_tree <- NULL

tmp_tree <- rtree(ntaxa(tmp_ps), rooted = TRUE,

tip.label = taxa_names(tmp_ps))

tmp_ps_norm <- merge_phyloseq(tmp_ps, sample_data, tmp_tree)

tmp_asv <- data.frame(t(otu_table(tmp_ps_norm)))

tmp_ps_name <- purrr::map_chr(i, ~ paste0(., "_norm"))

assign(tmp_ps_name, tmp_ps_norm)

tmp_asv_name <- purrr::map_chr(i, ~ paste0(., "_tu"))

assign(tmp_asv_name, tmp_asv)

rm(list = ls(pattern = "tmp_"))

}

```

Next, we run the analysis for all three metrics on the data sets (without a tree).

```{r}

#| results: hide

#| code-fold: true

qvalue <- c(0,1,2)

for (i in qvalue) {

for (j in ssu18_alpha_ds) {

tmp_asv <- get(purrr::map_chr(j, ~ paste0(., "_tu")))

tmp_df <- data.frame(hill_div(tmp_asv, qvalue = i))

tmp_df <- tmp_df %>% dplyr::rename("tmp_name" = 1) %>%

tibble::rownames_to_column("SamName")

tmp_name <- purrr::map_chr(j, ~ paste0(., "_h", i))

print(tmp_name)

assign(tmp_name, tmp_df)

rm(list = ls(pattern = "tmp_"))

}

}

objects(pattern = "_h")

```

And make summary tables to add back into each `ps` object.

```{r}

#| code-fold: true

for (i in ssu18_alpha_ds) {

tmp_obs <- get(purrr::map_chr(i, ~ paste0(., "_h0")))

tmp_sha <- get(purrr::map_chr(i, ~ paste0(., "_h1")))

tmp_sim <- get(purrr::map_chr(i, ~ paste0(., "_h2")))

tmp_hill <- dplyr::left_join(tmp_obs, tmp_sha, by = "SamName") %>%

dplyr::left_join(., tmp_sim, by = "SamName") %>%

dplyr::rename("Observed" = 2, "Shannon_exp" = 3, "InvSimpson" = 4)

tmp_name <- purrr::map_chr(i, ~ paste0(., "_hill"))

assign(tmp_name, tmp_hill)

rm(list = ls(pattern = "tmp_"))

}

objects(pattern = "_hill")

objects()

```

And then create the new objects with the diversity data.

```{r}

#| code-fold: true

for (i in ssu18_alpha_ds) {

tmp_ps <- get(i)

tmp_tree <- rtree(ntaxa(tmp_ps), rooted = TRUE,

tip.label = taxa_names(tmp_ps))

tmp_ps <- merge_phyloseq(tmp_ps, sample_data, tmp_tree)

tmp_hill <- get(purrr::map_chr(i, ~ paste0(., "_hill")))

tmp_hill_samp <- dplyr::left_join(data.frame(sample_data(tmp_ps)),

tmp_hill, by = "SamName")

tmp_hill_samp$ID <- tmp_hill_samp$SamName

tmp_hill_samp <- tmp_hill_samp %>% tibble::column_to_rownames("ID")

tmp_ps2 <- merge_phyloseq(otu_table(tmp_ps),

sample_data(tmp_hill_samp),

tax_table(tmp_ps),

phy_tree(tmp_ps))

assign(i, tmp_ps2)

tmp_path <- file.path("files/alpha/rdata/")

saveRDS(tmp_ps2, paste(tmp_path, i, ".rds", sep = ""))

rm(list = ls(pattern = "tmp_"))

}

objects()

```

```{r}

#| echo: false

## THIS is legacy code for using trees. Not included in this analysis

#### With a phylogenetic tree

## We can also run the tests using the phylogenetic tree to assess **lineage** diversity rather than ASV diversity for the filtered data sets. We do a quick check to ensure the ASV names in the ASV table and the tip names in the new phylogenetic tree are identical.

ssu18_ps_pime_tree_ult <- readRDS("files/alpha/rdata/ssu18_ps_pime_tree_ult.rds")

ssu18_ps_pime_otu_tree_ult <- readRDS("files/alpha/rdata/ssu18_ps_pime_otu_tree_ult.rds")

ssu18_ps_perfect_tree_ult <- readRDS("files/alpha/rdata/ssu18_ps_perfect_tree_ult.rds")

ssu18_ps_perfect_otu_tree_ult <- readRDS("files/alpha/rdata/ssu18_ps_perfect_otu_tree_ult.rds")

ssu18_ps_filt_tree_ult <- readRDS("files/alpha/rdata/ssu18_ps_filt_tree_ult.rds")

is.ultrametric(ssu18_ps_pime_tree_ult)

is.ultrametric(ssu18_ps_pime_otu_tree_ult)

is.ultrametric(ssu18_ps_perfect_tree_ult)

is.ultrametric(ssu18_ps_perfect_otu_tree_ult)

is.ultrametric(ssu18_ps_filt_tree_ult)

## To generate trees, we used the R-package [`phytools`](https://github.com/liamrevell/phytools), specifically the function `force.ultrametric` since `hill_div` requires ultrametric trees.

rm(i)

ssu18_alpha_ds_tree <- c("ssu18_ps_filt", "ssu18_ps_perfect", "ssu18_ps_pime")

for (i in ssu18_alpha_ds_tree) {

tmp_ps <- get(i)

tmp_tree <- tmp_ps@phy_tree

tmp_tree_ult <- force.ultrametric(tmp_tree, method = c("nnls"))

print(is.ultrametric(tmp_tree_ult))

tmp_tree_name <- purrr::map_chr(i, ~ paste0(., "_tree_ult"))

assign(tmp_tree_name, tmp_tree_ult)

tmp_path <- file.path("files/alpha/rdata/")

saveRDS(tmp_tree_ult, paste(tmp_path, tmp_tree_name, ".rds", sep = ""))

tmp_tmp_asv <- get(purrr::map_chr(i, ~ paste0(., "_tu")))

print(identical(sort(rownames(tmp_tmp_asv)), sort(tmp_tree_ult$tip.label)))

rm(list = ls(pattern = "tmp_"))

}

objects()

rm(i, j)

qvalue <- c(0,1,2)

for (i in qvalue) {

for (j in ssu18_alpha_ds_tree) {

tmp_asv <- get(purrr::map_chr(j, ~ paste0(., "_tu")))

tmp_tree <- get(purrr::map_chr(j, ~ paste0(., "_tree_ult")))

tmp_df <- data.frame(hill_div(tmp_asv, qvalue = i, tree = tmp_tree))

tmp_df <- tmp_df %>% dplyr::rename("tmp_name" = 1) %>%

tibble::rownames_to_column("SamName")

tmp_name <- purrr::map_chr(j, ~ paste0(., "_ht", i))

print(tmp_name)

assign(tmp_name, tmp_df)

rm(list = ls(pattern = "tmp_"))

}

}

ssu18_alpha_ds_tree_filt <- c("ssu18_ps_filt")

for (i in ssu18_alpha_ds_tree_filt) {

tmp_obs <- get(purrr::map_chr(i, ~ paste0(., "_ht0")))

tmp_sha <- get(purrr::map_chr(i, ~ paste0(., "_ht1")))

tmp_sim <- get(purrr::map_chr(i, ~ paste0(., "_ht2")))

tmp_hill <- dplyr::left_join(tmp_obs, tmp_sha, by = "SamName") %>%

dplyr::left_join(., tmp_sim, by = "SamName") %>%

dplyr::rename("Observed_fit" = 2, "Shannon_exp_fit" = 3, "InvSimpson_fit" = 4)

tmp_name <- purrr::map_chr(i, ~ paste0(., "_hill_t"))

assign(tmp_name, tmp_hill)

rm(list = ls(pattern = "tmp_"))

}

ssu18_alpha_ds_tree_perfect <- c("ssu18_ps_perfect")

for (i in ssu18_alpha_ds_tree_perfect) {

tmp_obs <- get(purrr::map_chr(i, ~ paste0(., "_ht0")))

tmp_sha <- get(purrr::map_chr(i, ~ paste0(., "_ht1")))

tmp_sim <- get(purrr::map_chr(i, ~ paste0(., "_ht2")))

tmp_hill <- dplyr::left_join(tmp_obs, tmp_sha, by = "SamName") %>%

dplyr::left_join(., tmp_sim, by = "SamName") %>%

dplyr::rename("Observed_pet" = 2, "Shannon_exp_pet" = 3, "InvSimpson_pet" = 4)

tmp_name <- purrr::map_chr(i, ~ paste0(., "_hill_t"))

assign(tmp_name, tmp_hill)

rm(list = ls(pattern = "tmp_"))

}

ssu18_alpha_ds_tree_pime <- c("ssu18_ps_pime")

for (i in ssu18_alpha_ds_tree_pime) {

tmp_obs <- get(purrr::map_chr(i, ~ paste0(., "_ht0")))

tmp_sha <- get(purrr::map_chr(i, ~ paste0(., "_ht1")))

tmp_sim <- get(purrr::map_chr(i, ~ paste0(., "_ht2")))

tmp_hill <- dplyr::left_join(tmp_obs, tmp_sha, by = "SamName") %>%

dplyr::left_join(., tmp_sim, by = "SamName") %>%

dplyr::rename("Observed_pit" = 2, "Shannon_exp_pit" = 3, "InvSimpson_pit" = 4)

tmp_name <- purrr::map_chr(i, ~ paste0(., "_hill_t"))

assign(tmp_name, tmp_hill)

rm(list = ls(pattern = "tmp_"))

}

rm(i)

for (i in ssu18_alpha_ds_tree) {

tmp_ps <- get(i)

tmp_hill <- get(purrr::map_chr(i, ~ paste0(., "_hill_t")))

tmp_hill_samp <- dplyr::left_join(data.frame(sample_data(tmp_ps)),

tmp_hill, by = "SamName")

tmp_hill_samp$ID <- tmp_hill_samp$SamName

tmp_hill_samp <- tmp_hill_samp %>% tibble::column_to_rownames("ID")

tmp_ps2 <- merge_phyloseq(otu_table(tmp_ps),

sample_data(tmp_hill_samp),

tax_table(tmp_ps),

phy_tree(tmp_ps))

assign(i, tmp_ps2)

tmp_path <- file.path("files/alpha/rdata/")

saveRDS(tmp_ps2, paste(tmp_path, i, ".rds", sep = ""))

rm(list = ls(pattern = "tmp_"))

}

```

## Hill Numbers Summary

Now we summarize the data for each sample against all three metrics. The table contains the results of ASV diversity estimates from the full data set and the three filtered data sets.

The suffix `_fi` indicates metrics for the Arbitrary data set. The suffix `_pe` indicates metrics for the PERfect data set and the suffix `_pet` indicates the lineage diversity for the PERfect data set.

The suffix `_pi` indicates metrics for the PIME data set and the suffix `_pit` indicates the lineage diversity for the PIME data set.

```{r}

#| echo: false

######################### FULL ASV

ssu18_work_hill_samp <- ssu18_ps_work_hill %>% dplyr::rename("Observed" = 2,

"Shannon_exp" = 3,

"InvSimpson" = 4)

ssu18_work_hill_samp <- dplyr::left_join(data.frame(sample_data(ssu18_ps_work)), ssu18_work_hill_samp)

ssu18_work_hill_samp <- ssu18_work_hill_samp[order(ssu18_work_hill_samp$PLOT, ssu18_work_hill_samp$TEMP),]

######################### Filt ASV

ssu18_filt_hill_samp <- ssu18_ps_filt_hill %>% dplyr::rename("Observed_fi" = 2,

"Shannon_exp_fi" = 3,

"InvSimpson_fi" = 4)

######################### PERfect ASV

ssu18_perfect_hill_samp <- ssu18_ps_perfect_hill %>% dplyr::rename("Observed_pe" = 2,

"Shannon_exp_pe" = 3,

"InvSimpson_pe" = 4)

######################### PIME ASV

ssu18_pime_hill_samp <- ssu18_ps_pime_hill %>% dplyr::rename("Observed_pi" = 2,

"Shannon_exp_pi" = 3,

"InvSimpson_pi" = 4)

```

```{r}

#| echo: false

######################### ASV TABLE

ssu_tab_alpha_div <- dplyr::left_join(ssu18_work_hill_samp, ssu18_filt_hill_samp) %>%

dplyr::left_join(., ssu18_perfect_hill_samp) %>%

dplyr::left_join(., ssu18_pime_hill_samp) %>%

dplyr::relocate(c("Observed", "Observed_fi", "Observed_pe", "Observed_pi"), .after = "PAIR") %>%

dplyr::relocate(c("Shannon_exp", "Shannon_exp_fi", "Shannon_exp_pe", "Shannon_exp_pi"), .after = "Observed_pi") %>%

dplyr::rename("Sample_ID" = "SamName")

ssu_tab_alpha_div <- ssu_tab_alpha_div%>% dplyr::mutate(across(.cols = c(11:ncol(ssu_tab_alpha_div)),

round, digits = 1))

```

### Hill Summary

<small>`r caption_tab_ssu("ssu_alpha_div")`</small>

```{r}

#| echo: false

#| eval: true

seq_table <- ssu_tab_alpha_div

seq_table[2:6] <- NULL

write.table(seq_table, "files/alpha/tables/ssu_tab_alpha_div.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

columns = list(

Sample_ID = colDef(name = "Sample_ID",

sticky = "left",

style = list(borderRight = "5px solid #eee"),

headerStyle = list(borderRight = "5px solid #eee"),

align = "left",

minWidth = 150),

Observed = colDef(name = "FULL"),

Observed_fi = colDef(name = "Arbitrary"),

Observed_pe = colDef(name = "PERfect"),

Observed_pi = colDef(name = "PIME",

style = list(borderRight = "5px solid #eee"),

headerStyle = list(borderRight = "5px solid #eee")),

Shannon_exp = colDef(name = "FULL"),

Shannon_exp_fi = colDef(name = "Arbitrary"),

Shannon_exp_pe = colDef(name = "PERfect"),

Shannon_exp_pi = colDef(name = "PIME",

style = list(borderRight = "5px solid #eee"),

headerStyle = list(borderRight = "5px solid #eee")),

InvSimpson = colDef(name = "FULL"),

InvSimpson_fi = colDef(name = "Arbitrary"),

InvSimpson_pe = colDef(name = "PERfect"),

InvSimpson_pi = colDef(name = "PIME")

),

columnGroups = list(

colGroup(name = "Observed richness",

columns = c("Observed", "Observed_fi", "Observed_pe", "Observed_pi"),

headerStyle = list(borderRight = "5px solid #eee", fontSize = "1.1em")),

colGroup(name = "Shannon exponential",

columns = c("Shannon_exp", "Shannon_exp_fi", "Shannon_exp_pe", "Shannon_exp_pi"),

headerStyle = list(fontSize = "1.1em")),

colGroup(name = "Simpson multiplicative inverse",

columns = c("InvSimpson", "InvSimpson_fi", "InvSimpson_pe", "InvSimpson_pi"),

headerStyle = list(fontSize = "1.1em"))

),

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = TRUE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100

),

searchable = TRUE, defaultPageSize = 5, showPagination = TRUE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = TRUE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/alpha/tables/ssu_tab_alpha_div.txt dname="ssu_tab_alpha_div" label="Download table as text file" icon="table" type="primary" >}}

## Normality Tests

Before running significance tests, we need to know if data is normally distributed, which will tell use whether to use a parametric or non-parametric test. To test if the data are normally distributed, we use the Shapiro-Wilk Normality test and the Bartlett Test of Homogeneity of Variances.

If the p-values are both **not significant** (p > 0.05) from the tests, we accept the null hypothesis (that the results are normally distributed) and test for significance between samples using an ANOVA. If the p-values are both **significant** (p < 0.05), we reject the null hypothesis (that the results are normally distributed) and test for significance between samples using Kruskal-Wallis (non-parametric equivalent of ANOVA).

The commands are as follows:

`shapiro.test(x)`, where `x` is a numeric vector of alpha diversity values from the sample data table.

`bartlett.test(Value ~ Group, data = df)` Where `Value` is the metric of interest, `Group` in the treatment to compare, and `df` is the data frame.

First the Shapiro-Wilk Normality test.

```{r}

#| code-fold: true

ssu_tab_alpha_div

ssu18_div_tab <- ssu_tab_alpha_div

ssu18_shap_tests_asv <- c()

for (i in colnames(ssu18_div_tab[,7:ncol(ssu18_div_tab)])) {

tmp_name <- purrr::map_chr(i, ~ paste0("ssu18_shap_asv_", .))

ssu18_shap_tests_asv <- append(ssu18_shap_tests_asv, tmp_name)

tmp_test <- eval(shapiro.test(ssu18_div_tab[[i]]))

tmp_test$data.name <- tmp_name

assign(tmp_name, tmp_test)

rm(list = ls(pattern = "tmp_"))

}

objects(pattern = "shap_")

```

And then the Bartlett Test of Homogeneity of Variances.

```{r}

#| code-fold: true

ssu18_div_tab_asv <- ssu_tab_alpha_div

ssu18_bart_tests_asv <- c()

for (i in colnames(ssu18_div_tab_asv[,7:ncol(ssu18_div_tab_asv)])) {

tmp_name <- purrr::map_chr(i, ~ paste0("ssu18_bart_asv_", .))

ssu18_bart_tests_asv <- append(ssu18_bart_tests_asv, tmp_name)

tmp_test <- eval(bartlett.test(ssu18_div_tab_asv[[i]] ~ TEMP, data = ssu18_div_tab_asv))

tmp_test$data.name <- tmp_name

assign(tmp_name, tmp_test)

rm(list = ls(pattern = "tmp_"))

}

objects(pattern = "_bart_")

```

Here we see which Shapiro-Wilk Normality **and** Bartlett tests were significant and which were not. So wherever the value of both p-values in `> 0.05` we can use an ANOVA, otherwise we use Kruskal-Wallis.

```{r}

#| echo: false

#| eval: true

#| warning: false

#| message: false

ssu18_norm_res <- data.frame()

for (i in colnames(ssu18_div_tab_asv[,7:ncol(ssu18_div_tab_asv)])) {

tmp_get_shap <- get(purrr::map_chr(i, ~ paste0("ssu18_shap_asv_", .)))

tmp_shap_p <- round(tmp_get_shap[[2]], 4)

tmp_get_bart <- get(purrr::map_chr(i, ~ paste0("ssu18_bart_asv_", .)))

tmp_bart_p <- round(tmp_get_bart[[3]], 4)

tmp_test <- if(tmp_shap_p > 0.05 & tmp_bart_p > 0.05 ){

tmp_test <- "ANOVA"

} else {

tmp_test <- "Kruskal-Wallis"

}

tmp_df <- data.frame(i, tmp_shap_p, tmp_bart_p, tmp_test)

ssu18_norm_res <- dplyr::bind_rows(ssu18_norm_res, tmp_df)

rm(list = ls(pattern = "tmp_"))

}

ssu18_norm_res[,1] <- gsub('_exp', '', ssu18_norm_res[,1])

ssu18_norm_res <- ssu18_norm_res %>% separate(col = 1,

into = c("metric", "dataset"),

sep = "_", remove = TRUE)

ssu18_norm_res <- ssu18_norm_res %>% tidyr::replace_na(list(dataset = "FULL"))

ssu18_norm_res$dataset <- stringr::str_replace(ssu18_norm_res$dataset, "^fi$", "FILT") %>%

stringr::str_replace(., "^pi$", "PIME") %>%

stringr::str_replace(., "^pe$", "PERfect")

ssu18_norm_res <- ssu18_norm_res %>%

dplyr::rename(c("p-value (normality)" = 3, "p-value (homogeneity)" = 4, "method" = 5))

```

<small>`r caption_tab_ssu("ssu_norm_tests")`</small>

```{r}

#| echo: false

#| eval: true

seq_table <- ssu18_norm_res

write.table(ssu18_norm_res, "files/alpha/tables/ssu_norm_tests.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = FALSE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100, headerVAlign = "center"

),

columns = list(

metric = colDef(name = "Metric",

sticky = "left",

style = list(borderRight = "1px solid #eee"),

headerStyle = list(borderRight = "1px solid #eee"),

align = "left",

minWidth = 100),

`p-value (normality)` = colDef(name = "p-value (normality)", minWidth = 100, resizable = FALSE,

style = function(value) {

if (value > 0.05) {

color <- "#2271B2"

} else if (value < 0.05) {

color <- "#B22271"

} else {

color <- "#777"

}

list(color = color, fontWeight = "bold")

}),

`p-value (homogeneity)` = colDef(name = "p-value (homogeneity)", minWidth = 130, resizable = FALSE,

style = function(value) {

if (value > 0.05) {

color <- "#2271B2" #green "#", "#"

} else if (value < 0.05) {

color <- "#B22271" #RED

} else {

color <- "#777"

}

list(color = color, fontWeight = "bold")

}),

`dataset` = colDef(name = "dataset", minWidth = 80, resizable = FALSE),

`method` = colDef(name = "method", minWidth = 120, resizable = FALSE)

),

searchable = FALSE, defaultPageSize = 12, showPagination = FALSE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = TRUE, showSortable = FALSE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.8em")))

rm(seq_table)

```

{{< downloadthis files/alpha/tables/ssu_norm_tests.txt dname="ssu_norm_tests" label="Download table as text file" icon="table" type="primary" >}}

## Significance Tests

To begin, we need to create a hierarchy variable; a two-column matrix specifying the relationship between samples (first column) and groups (second column).

```{r}

#| code-fold: true

ssu18_hill_hier <- ssu_tab_alpha_div

ssu18_hill_hier <- ssu18_hill_hier %>% dplyr::select("Sample_ID", "TEMP") %>%

tibble::remove_rownames()

ssu18_hill_hier <- ssu18_hill_hier[order(ssu18_hill_hier$TEMP), ]

ssu18_hill_hier$TEMP = paste(ssu18_hill_hier$TEMP, 'C', sep='')

ssu18_hill_hier <- ssu18_hill_hier %>% tibble::remove_rownames()

saveRDS(ssu18_hill_hier, "files/alpha/rdata/ssu18_hill_hier.rds")

```

Again, we start by testing significance of ASV diversity for the data sets against each of the three metrics using the `div_test` function.

The command is as follows:

```

div_test(countable = x, qvalue = i, hierarchy = hier,

tree = ultrametric_tree, posthoc = TRUE)

```

where `x` is ASV by sample table, `i` is the q-value corresponding to the metric of interest, `hier` is the hierarchy matrix, `tree` is an ultrametric formatted phylogenetic tree (if you choose to explore lineage diversity), and `posthoc` indicates whether to run post hoc pairwise analyses.

```{r}

#| results: hide

#| code-fold: true

qvalue <- c(0,1,2)

for (i in ssu18_alpha_ds) {

for (j in qvalue) {

tmp_get <- get(purrr::map_chr(i, ~ paste0(., "_tu")))

tmp_test <- div_test(tmp_get, qvalue = j,

hierarchy = ssu18_hill_hier,

posthoc = TRUE)

tmp_name <- purrr::map_chr(i, ~ paste0(., "_q", j, "_adt"))

print(tmp_name)

assign(tmp_name, tmp_test)

rm(list = ls(pattern = "tmp_"))

}

}

```

```{r}

#| echo: false

## Legacy code for lineage diversity

## And then test significance of lineage diversity for the filtered data sets.

ssu18_alpha_ds_tree <- c("ssu18_ps_filt", "ssu18_ps_perfect", "ssu18_ps_pime")

qvalue <- c(0,1,2)

for (i in ssu18_alpha_ds_tree) {

for (j in qvalue) {

tmp_get <- get(purrr::map_chr(i, ~ paste0(., "_tu")))

tmp_tree <- get(purrr::map_chr(i, ~ paste0(., "_tree_ult")))

tmp_test <- div_test(tmp_get, qvalue = j,

hierarchy = ssu18_hill_hier,

posthoc = TRUE,

tree = tmp_tree)

tmp_name <- purrr::map_chr(i, ~ paste0(., "_t_q", j, "_adt"))

print(tmp_name)

assign(tmp_name, tmp_test)

rm(list = ls(pattern = "tmp_"))

}

}

```

```{r}

#| echo: false

tmp_objects <- c("ssu18_ps_work", "ssu18_ps_filt", "ssu18_ps_perfect", "ssu18_ps_pime")

tmp_metric <- data.frame(c("Observed", "Shannon exponential", "Inverse Simpson"))

tmp_qvalue <- data.frame(c("0", "1", "2"))

qvalue <- c(0,1,2)

for (i in tmp_objects) {

tmp_h_pvalue <- c()

for (j in qvalue) {

tmp_name <- purrr::map_chr(i, ~ paste0(., "_q", j, "_adt"))

tmp_get <- get(tmp_name)$homogeneity.pvalue

tmp_h_pvalue <- c(append(tmp_h_pvalue, tmp_get))

}

tmp_h_pvalue <- data.frame(tmp_h_pvalue)

tmp_n_pvalue <- c()

for (j in qvalue) {

tmp_name <- purrr::map_chr(i, ~ paste0(., "_q", j, "_adt"))

tmp_get <- get(tmp_name)$normality.pvalue

tmp_n_pvalue <- c(append(tmp_n_pvalue, tmp_get))

}

tmp_n_pvalue <- data.frame(tmp_n_pvalue)

tmp_method <- c()

for (j in qvalue) {

tmp_name <- purrr::map_chr(i, ~ paste0(., "_q", j, "_adt"))

tmp_get <- get(tmp_name)$method

tmp_method <- c(append(tmp_method, tmp_get))

}

tmp_method <- data.frame(tmp_method)

tmp_phoc_method <- c()

for (j in qvalue) {

tmp_name <- purrr::map_chr(i, ~ paste0(., "_q", j, "_adt"))

tmp_get <- get(tmp_name)$posthoc.method

tmp_phoc_method <- c(append(tmp_phoc_method, tmp_get))

}

tmp_phoc_method <- data.frame(tmp_phoc_method)

tmp_df <- dplyr::bind_cols(tmp_metric, tmp_qvalue) %>%

dplyr::bind_cols(., tmp_n_pvalue) %>%

dplyr::bind_cols(., tmp_h_pvalue) %>%

dplyr::bind_cols(., tmp_method) %>%

dplyr::bind_cols(., tmp_phoc_method) %>%

dplyr::rename("metric" = 1, "q-value" = 2, "normality p-value" = 3,

"homogeneity p-value" = 4, "method" = 5, "posthoc method" = 6)

tmp_name <- purrr::map_chr(i, ~ paste0(i, "_sig_tab"))

assign(tmp_name, tmp_df)

}

```

```{r}

#| echo: false

## Format as follows:

## FALSE Kruskal-Wallis Test = its18_pime_q1_adt$test[[3]]

## TRUE Tukey post-hoc = its18_ps_work_q0_adt$test[[1]][[1, 5]]

############################### WORK ##################################################

tmp_pvalue <- data.frame(c(ssu18_ps_work_q0_adt$test[[1]][[1, 5]],

ssu18_ps_work_q1_adt$test[[1]][[1, 5]],

ssu18_ps_work_q2_adt$test[[1]][[1, 5]]))

ssu18_ps_work_sig_tab <- dplyr::bind_cols(ssu18_ps_work_sig_tab, tmp_pvalue) %>%

dplyr::rename("posthoc p-value" = 7)

####################################Filt#############################################

tmp_pvalue <- data.frame(c(ssu18_ps_filt_q0_adt$test[[1]][[1, 5]],

ssu18_ps_filt_q1_adt$test[[1]][[1, 5]],

ssu18_ps_filt_q2_adt$test[[1]][[1, 5]]))

ssu18_ps_filt_sig_tab <- dplyr::bind_cols(ssu18_ps_filt_sig_tab, tmp_pvalue) %>%

dplyr::rename("posthoc p-value" = 7)

####################################PERfect##########################################

tmp_pvalue <- data.frame(c(ssu18_ps_perfect_q0_adt$test[[1]][[1, 5]],

ssu18_ps_perfect_q1_adt$test[[1]][[1, 5]],

ssu18_ps_perfect_q2_adt$test[[1]][[1, 5]]))

ssu18_ps_perfect_sig_tab <- dplyr::bind_cols(ssu18_ps_perfect_sig_tab, tmp_pvalue) %>%

dplyr::rename("posthoc p-value" = 7)

################################ PIME ################################################

tmp_pvalue <- data.frame(c(ssu18_ps_pime_q0_adt$test[[1]][[1, 5]],

ssu18_ps_pime_q1_adt$test[[3]],

ssu18_ps_pime_q2_adt$test[[1]][[1, 5]]))

ssu18_ps_pime_sig_tab <- dplyr::bind_cols(ssu18_ps_pime_sig_tab, tmp_pvalue) %>%

dplyr::rename("posthoc p-value" = 7)

rm(list = ls(pattern = "tmp_"))

```

```{r}

#| echo: false

ssu18_ps_work_sig_tab$dataset <- "FULL"

ssu18_ps_work_sig_tab$type <- "ASV"

ssu18_ps_work_sig_tab$lineage <- "No"

ssu18_ps_work_sig_tab <- ssu18_ps_work_sig_tab[,c(8:10,1:7)]

ssu18_ps_pime_sig_tab$dataset <- "PIME"

ssu18_ps_pime_sig_tab$type <- "ASV"

ssu18_ps_pime_sig_tab$lineage <- "No"

ssu18_ps_pime_sig_tab <- ssu18_ps_pime_sig_tab[,c(8:10,1:7)]

ssu18_ps_perfect_sig_tab$dataset <- "PERfect"

ssu18_ps_perfect_sig_tab$type <- "ASV"

ssu18_ps_perfect_sig_tab$lineage <- "No"

ssu18_ps_perfect_sig_tab <- ssu18_ps_perfect_sig_tab[,c(8:10,1:7)]

ssu18_ps_filt_sig_tab$dataset <- "FILT"

ssu18_ps_filt_sig_tab$type <- "ASV"

ssu18_ps_filt_sig_tab$lineage <- "No"

ssu18_ps_filt_sig_tab <- ssu18_ps_filt_sig_tab[,c(8:10,1:7)]

ssu18_sig_tab_all <- rbind(ssu18_ps_work_sig_tab,

ssu18_ps_filt_sig_tab,

ssu18_ps_pime_sig_tab,

ssu18_ps_perfect_sig_tab

)

```

## Summary

<small>`r caption_tab_ssu("ssu_sign_tests")`</small>

```{r}

#| echo: false

#| eval: true

#| warning: false

#| message: false

seq_table <- ssu18_sig_tab_all

seq_table[,2:3] <- NULL

seq_table <- seq_table %>% dplyr::mutate(across(.cols = c(4:5), round, digits = 4))

write.table(seq_table, "files/alpha/tables/ssu_sign_tests.txt",

row.names = FALSE, sep = "\t", quote = FALSE)

```

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(

header = function(value) gsub("_", " ", value, fixed = TRUE),

cell = function(value) format(value, nsmall = 0),

align = "center", filterable = FALSE, sortable = TRUE, resizable = TRUE,

footerStyle = list(fontWeight = "bold"), minWidth = 100

),

columns = list(

dataset = colDef(name = "Dataset",

sticky = "left",

style = list(borderRight = "1px solid #eee"),

headerStyle = list(borderRight = "1px solid #eee"),

align = "left",

minWidth = 80),

`posthoc p-value` = colDef(name = "posthoc p-value",

style = function(value) {

if (value > 0.05) {

color <- "#2271B2" #green

} else if (value < 0.05) {

color <- "#B22271" #RED "#", "#"

} else {

color <- "#777"

}

list(color = color, fontWeight = "bold")

}, cell = function(value) format(value, digits = 3, scientific = TRUE)

),

metric = colDef(name = "metric", minWidth = 140, align = "left"),

`q-value` = colDef(name = "q-value", minWidth = 70, align = "center"),

`normality p-value` = colDef(name = "p-value (normality)", minWidth = 100,),

`homogeneity p-value` = colDef(name = "p-value (homogeneity)", minWidth = 130),

`dataset` = colDef(name = "dataset", minWidth = 80),

`method` = colDef(name = "method", minWidth = 120)

),

searchable = FALSE, defaultPageSize = 12, showPagination = FALSE,

pageSizeOptions = c(5, 10, nrow(seq_table)),

showPageSizeOptions = TRUE, highlight = TRUE,

bordered = TRUE, striped = TRUE, compact = TRUE,

wrap = FALSE, showSortable = FALSE, fullWidth = TRUE,

theme = reactableTheme(style = list(fontSize = "0.7em")))

#rm(seq_table)

```

{{< downloadthis files/alpha/tables/ssu_sign_tests.txt dname="ssu_sign_tests" label="Download table as text file" icon="table" type="primary" >}}

## PostHoc Analyses

First let's check the results of each posthoc analysis.

::: {.column-body}

::: {.panel-tabset}

### Observed (q-value = 0)

```{r}

#| echo: false

#| eval: true

#| results: hold

#| comment: ''

ssu18_asvf0_lab <- "FULL (Observed)"

ssu18_asvf0_lab

ssu18_ps_work_q0_adt$posthoc.method

data.frame(ssu18_ps_work_q0_adt$posthoc)

```

```{r}

#| echo: false

#| eval: true

#| results: hold

#| comment: ''

ssu18_asvl0_lab <- "FILT (Observed)"

ssu18_asvl0_lab

ssu18_ps_filt_q0_adt$posthoc.method

data.frame(ssu18_ps_filt_q0_adt$posthoc)

```