---

title: "8. Metadata Analysis"

description: |

Workflow to assess the relationship between environmental parameters and microbial communities.

---

```{r}

#| message: false

#| results: hide

#| code-fold: true

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

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

set.seed(119)

#library(conflicted)

#pacman::p_depends(vegan, local = TRUE)

#pacman::p_depends_reverse(vegan, local = TRUE)

library(phyloseq); packageVersion("phyloseq")

pacman::p_load(ape, bestNormalize, tidyverse, gdata, mctoolsr,

ggpubr, ggvegan, ggfortify, vegan, reactable,

downloadthis, captioner, ggplot2,

install = FALSE, update = FALSE)

options(scipen=999)

knitr::opts_current$get(c(

"cache",

"cache.path",

"cache.rebuild",

"dependson",

"autodep"

))

```

```{r}

#| include: false

#| eval: true

#| comment: "Load for summary #3"

remove(list = ls())

load("page_build/metadata_summary.rdata")

objects()

```

```{r}

#| echo: false

#| eval: true

#| results: hide

#| comment: "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")

caption_tab <- captioner(prefix = "(common) Table", suffix = " |", style = "b")

caption_fig <- captioner(prefix = "(common) 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}

#| echo: false

#| eval: true

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

```

# Synopsis

In this workflow we compare the environmental metadata with the microbial community data. Part 1 involves the [16S rRNA](#s-rrna-1) community data and Part 2 deals with the [ITS](#its-1) community data. Each workflow contains the same major steps:

1) Metadata Normality Tests: Shapiro-Wilk Normality Test to test whether each matadata parameter is normally distributed.

2) Normalize Parameters: R package [`bestNormalize`](https://github.com/petersonR/bestNormalize) to find and execute the best normalizing transformation.

3) Split Metadata parameters into groups: a) Environmental and edaphic properties, b) Microbial functional responses, and c) Temperature adaptation properties.

4) Autocorrelation Tests: Test all possible pair-wise comparisons, on both normalized and non-normalized data sets, for each group.

5) Remove autocorrelated parameters from each group.

6) Dissimilarity Correlation Tests: Use Mantel Tests to see if any on the metadata groups are significantly correlated with the community data.

7) Best Subset of Variables: Determine which of the metadata parameters from each group are the most strongly correlated with the community data. For this we use the `bioenv` function from the `vegan` [@oksanen2013community] package.

8) Distance-based Redundancy Analysis: Ordination analysis of samples and metadata vector overlays using `capscale`, also from the `vegan` package.

## Workflow Input

Files needed to run this workflow can be downloaded from figshare. You will also need the metadata table, which can be accessed below.

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

## Metadata

In this section we provide access to the complete environmental metadata containing 61 parameters across all 15 samples.

<br/>

```{r}

#| echo: false

#| eval: true

metad <- read.table("files/metadata/tables/metadata.txt",

header = TRUE, sep = "\t")

seq_table <- metad

```

<small>`r caption_tab("metadata_results")`</small>

```{r}

#| echo: false

#| eval: true

reactable(seq_table,

defaultColDef = colDef(headerStyle = list(fontSize = "0.8em"),

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

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

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

footerStyle = list(fontWeight = "bold")

),

columns = list(

Sample_ID = colDef(name = "Sample ID",

sticky = "left",

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

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

align = "left",

minWidth = 150),

raw_reads = colDef(footer = function(values) sprintf("%.0f", sum(values)))

),

searchable = TRUE, defaultPageSize = 5,

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

showPageSizeOptions = TRUE, highlight = TRUE,

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

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

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

rm(seq_table)

```

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

# Multivariate Summary

Summary results for of `envfit` and `bioenv` tests for edaphic properties, soil functional response, and temperature adaptation for both 16S rRNA and ITS data sets.

## 16S rRNA

Metadata parameters removed based on autocorrelation results.

- edaphic properties: TEB + DON + Na + Al + Ca

- soil functional response: micN + micNP + enzCN + enzCP + BP~ase~ + CE~ase~ + LP~ase~ + N~ase~ + P~ase~

- temperature adaptation: NUE + PUE + SI

Metadata parameters removed from capscale analysis based on Degrees of Freedom.

- edaphic properties: Mg + Mn

- soil functional response: NONE

- temperature adaptation: NONE

<br/>

```{r}

#| echo: false

#| eval: true

seq_table <- read.table("files/metadata/tables/ssu18_sig_envfit_bioenv_summary.txt", header = TRUE, sep = "\t")

```

```{r}

#| echo: false

#| eval: true

#| warning: false

#| message: false

seq_table <- ssu18_sig_summary %>% dplyr::rename("envfit_r2" = 3,

"envfit_pval" = 4,

"bioenv_r" = 5,

"bioenv_pval" = 6)

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

seq_table[is.na(seq_table)] <- ""

```

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

```{r}

#| echo: false

#| eval: true

reactable(seq_table, groupBy = "metadata_set", defaultExpanded = TRUE,

defaultColDef = colDef(

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

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

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

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

columns = list(

metadata_set = colDef(name = "Metadata group",

sticky = "left",

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

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

align = "left",

minWidth = 90),

parameter = colDef(name = "parameter",

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

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

envfit_r2 = colDef(name = "r2"),

envfit_pval = colDef(name = "P(>r)",

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",

borderRight = "5px solid #eee")},

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

bioenv_r = colDef(name = "r", minWidth = 40),

bioenv_pval = colDef(name = "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")})),

columnGroups = list(

colGroup(name = "envfit",

columns = c("envfit_r2", "envfit_pval"), headerStyle = list(fontSize = "1.2em")),

colGroup(name = "bioenv",

columns = c("bioenv_r", "bioenv_pval"), headerStyle = list(fontSize = "1.2em"))),

searchable = FALSE, defaultPageSize = 15, 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.8em")))

rm(seq_table)

```

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

## ITS

Metadata parameters removed based on autocorrelation results.

- edaphic properties: TEB + DON + Na + Al + Ca

- soil functional response: micN + micNP + enzCN + enzCP + BP~ase~ + CE~ase~ + LP~ase~ + N~ase~ + P~ase~

- temperature adaptation: NUE + PUE + P~Q10~ + SI

Metadata parameters removed from capscale analysis based on Degrees of Freedom.

- edaphic properties: Mg + Mn + Na + Al + Fe + K

- soil functional response: NONE

- temperature adaptation: S~Q10~

<br/>

```{r}

#| echo: false

#| eval: true

seq_table <- read.table("files/metadata/tables/its18_sig_envfit_bioenv_summary.txt",

header = TRUE, sep = "\t")

```

```{r}

#| echo: false

#| eval: true

seq_table <- its18_sig_summary %>% dplyr::rename("envfit_r2" = 3,

"envfit_pval" = 4,

"bioenv_r" = 5,

"bioenv_pval" = 6)

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

seq_table[is.na(seq_table)] <- ""

```

<small>`r caption_tab_its("its_multi_summ")`</small>

```{r}

#| echo: false

#| eval: true

reactable(seq_table, groupBy = "metadata_set", defaultExpanded = TRUE,

defaultColDef = colDef(

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

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

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

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

columns = list(

metadata_set = colDef(name = "Metadata group",

sticky = "left",

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

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

align = "left",

minWidth = 90),

parameter = colDef(name = "parameter",

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

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

envfit_r2 = colDef(name = "r2"),

envfit_pval = colDef(name = "P(>r)",

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",

borderRight = "5px solid #eee")},

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

bioenv_r = colDef(name = "r", minWidth = 40),

bioenv_pval = colDef(name = "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")})),

columnGroups = list(

colGroup(name = "envfit",

columns = c("envfit_r2", "envfit_pval"), headerStyle = list(fontSize = "1.2em")),

colGroup(name = "bioenv",

columns = c("bioenv_r", "bioenv_pval"), headerStyle = list(fontSize = "1.2em"))),

searchable = FALSE, defaultPageSize = 15, 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.8em")))

rm(seq_table)

```

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

```{r}

#| include: false

#| eval: true

#| comment: "Load to build page only #2"

load("page_build/metadata_ssu18_wf.rdata")

objects()

```

```{r}

#| echo: false

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

tmp_1 <- ssu18_ps_work_mc$data_loaded

tmp_2 <- ssu18_ps_pime_mc$data_loaded

tmp_1 <- tmp_1 %>% tibble::rownames_to_column("ID")

tmp_2 <- tmp_2 %>% tibble::rownames_to_column("ID")

tmp_3 <- left_join(tmp_2, tmp_1, by = "ID")

tmp_3 <- tmp_3[ , order(names(tmp_3))]

write.table(tmp_3, "tmp_3.txt", sep = "\t", quote = FALSE)

```

# 16s rRNA

## Data Formating

```{r}

#| include: false

#| comment: "Initial Load for ANALYSIS #1"

set.seed(119)

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

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

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

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

```

1. Transform `ps` Objects

The first step is to transform the phyloseq objects.

```{r}

#| code-fold: true

samp_ps <- c("ssu18_ps_work", "ssu18_ps_pime", "ssu18_ps_perfect",

"ssu18_ps_work_otu", "ssu18_ps_pime_otu", "ssu18_ps_perfect_otu")

for (i in samp_ps) {

tmp_get <- get(i)

tmp_ps <- transform_sample_counts(tmp_get, function(otu) 1e5 * otu/sum(otu))

tmp_ps@phy_tree <- NULL

tmp_ps <- prune_samples(sample_sums(tmp_ps) > 0, tmp_ps)

tmp_tree <- rtree(ntaxa(tmp_ps), rooted = TRUE, tip.label = taxa_names(tmp_ps))

tmp_samp <- data.frame(sample_data(tmp_get))

tmp_samp <- tmp_samp %>% dplyr::rename("TREAT_T" = "TEMP")

tmp_ps <- merge_phyloseq(tmp_ps, sample_data, tmp_tree)

sample_data(tmp_ps) <- tmp_samp

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

print(tmp_name)

assign(tmp_name, tmp_ps)

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

}

```

2. Format for `mctoolsr`

Then a bit of formatting to make the data compatible with `mctoolsr`.

```{r}

#| code-fold: true

for (i in samp_ps) {

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

tmp_tax <- data.frame(tax_table(tmp_get))

tmp_tax$ASV_SEQ <- NULL

tmp_col_names <- colnames(tmp_tax)

tmp_tax_merge <- tmp_tax %>% tidyr::unite(taxonomy,

all_of(tmp_col_names),

sep = ";")

tmp_tax_merge <- tmp_tax_merge %>% tibble::rownames_to_column("#OTU_ID")

tmp_otu <- data.frame(t(otu_table(tmp_get)))

tmp_otu <- tmp_otu %>% tibble::rownames_to_column("#OTU_ID")

tmp_otu_tax <- dplyr::left_join(tmp_otu, tmp_tax_merge, by = "#OTU_ID")

tmp_samp <- data.frame(sample_data(tmp_get))

tmp_samp[,c(1,3)] <- NULL

tmp_samp <- tmp_samp %>% tibble::rownames_to_column("#SampleID")

tmp_metad <- metad %>% dplyr::rename("#SampleID" = "id")

tmp_metad[,2:5] <- NULL

tmp_md <- dplyr::left_join(tmp_samp, tmp_metad, by = "#SampleID")

tmp_otu_name <- purrr::map_chr(i, ~ paste0(., "_otu_tax"))

print(tmp_otu_name)

assign(tmp_otu_name, tmp_otu_tax)

tmp_md_name <- purrr::map_chr(i, ~ paste0(., "_md"))

print(tmp_md_name)

assign(tmp_md_name, tmp_md)

tmp_path <- file.path("files/metadata/tables/")

write_delim(tmp_otu_tax, paste(tmp_path, tmp_otu_name, ".txt", sep = ""), delim = "\t")

write_delim(tmp_md, paste(tmp_path, tmp_md_name, ".txt", sep = ""), delim = "\t")

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

}

```

```{r}

#| code-fold: true

for (i in samp_ps) {

tmp_path <- file.path("files/metadata/tables/")

tmp_otu_name <- purrr::map_chr(i, ~ paste0(., "_otu_tax"))

tmp_md_name <- purrr::map_chr(i, ~ paste0(., "_md"))

tmp_tax_table_fp <- paste(tmp_path, tmp_otu_name, ".txt", sep = "")

tmp_map_fp <- paste(tmp_path, tmp_md_name, ".txt", sep = "")

tmp_input <- load_taxa_table(tmp_tax_table_fp, tmp_map_fp)

tmp_input_name <- purrr::map_chr(i, ~ paste0(., "_mc"))

print(tmp_input_name)

assign(tmp_input_name, tmp_input)

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

}

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

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

```

## Choose Data Set

At this point in the code we need to choose a data set to use, formatted with `mctoolsr`. Remember, there are four choices:

1) FULL, unfiltered data set.

2) Arbitrary filtered data set.

3) PERfect filtered data set.

4) PIME filtered data set.

This way, if we want to test other data sets we only need to change the name here.

```{r}

objects(pattern = "_mc")

ssu18_select_mc <- ssu18_ps_pime_mc

```

## Normality Tests

Before proceeding, we need to test each parameter in the metadata to see which ones are and are not normally distributed. For that, we use the Shapiro-Wilk Normality Test. Here we only need one of the metadata files.

```{r}

#| code-fold: true

temp_md <- ssu18_select_mc$map_loaded

temp_md[,1:9] <- NULL

shap_results <- NULL

for (i in colnames(temp_md)) {

tmp_shap <- shapiro.test(temp_md[[i]])

tmp_p <- round(tmp_shap$p.value, digits = 5)

tmp_res <- eval(isTRUE(tmp_shap$p.value < 0.05))

shap_results <- rbind(shap_results, data.frame(i, tmp_p, tmp_res))

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

}

colnames(shap_results) <- c("parameter", "p-value", "tranform")

shap_results

dplyr::filter(shap_results, tranform == "TRUE")

md_to_tranform <- shap_results$parameter[shap_results$tranform == TRUE]

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

```

<details>

<summary>**Show the results** of each normality test for metadata parameters</summary>

</br>

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

```{r}

#| echo: false

#| eval: true

seq_table <- shap_results

reactable(seq_table,

defaultColDef = colDef(headerStyle = list(fontSize = "0.8em"),

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

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

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

footerStyle = list(fontWeight = "bold")

),

columns = list(

parameter = colDef(name = "Parameter",

sticky = "left",

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

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

align = "left",

minWidth = 150),

`p-value` = colDef(name = "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")})

),

searchable = TRUE, defaultPageSize = 60,

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

showPageSizeOptions = FALSE, highlight = TRUE,

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

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

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

rm(seq_table)

```

</details>

Looks like we need to transform `r length(md_to_tranform)` metadata parameters.

## Normalize Parameters

Here we use the R package [`bestNormalize`](https://cran.r-project.org/web/packages/bestNormalize/vignettes/bestNormalize.html) to find and execute the best normalizing transformation. The function will test the following normalizing transformations:

- `arcsinh_x` performs an arcsinh transformation.

- `boxcox` Perform a Box-Cox transformation and center/scale a vector to attempt normalization. `boxcox` estimates the optimal value of lambda for the Box-Cox transformation. The function will return an error if a user attempt to transform nonpositive data.

- `yeojohnson` Perform a Yeo-Johnson Transformation and center/scale a vector to attempt normalization. `yeojohnson` estimates the optimal value of lambda for the Yeo-Johnson transformation. The Yeo-Johnson is similar to the Box-Cox method, however it allows for the transformation of nonpositive data as well.

- `orderNorm` The Ordered Quantile (ORQ) normalization transformation, `orderNorm()`, is a rank-based procedure by which the values of a vector are mapped to their percentile, which is then mapped to the same percentile of the normal distribution. Without the presence of ties, this essentially guarantees that the transformation leads to a uniform distribution.

- `log_x` performs a simple log transformation. The parameter a is essentially estimated by the training set by default (estimated as the minimum possible to some extent epsilon), while the base must be specified beforehand. The default base of the log is 10.

- `sqrt_x` performs a simple square-root transformation. The parameter a is essentially estimated by the training set by default (estimated as the minimum possible), while the base must be specified beforehand.

- `exp_x` performs a simple exponential transformation.

See this GitHub issue ([#5](https://github.com/petersonR/bestNormalize/issues/5)) for a description on getting reproducible results. Apparently, you can get different results because the `bestNormalize()` function uses repeated cross-validation (and doesn't automatically set the seed), so the results will be slightly different each time the function is executed.

```{r}

#| code-fold: true

set.seed(119)

for (i in md_to_tranform) {

tmp_md <- ssu18_select_mc$map_loaded

tmp_best_norm <- bestNormalize(tmp_md[[i]], r = 1, k = 5, loo = TRUE)

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

assign(tmp_name, tmp_best_norm)

print(tmp_name)

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

}

```

<details>

<summary>Show the chosen transformations</summary>

```{r}

#| echo: false

#| eval: true

for (i in md_to_tranform) {

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

tmp_print <- c("bestNormalize Chosen transformation of", i)

cat("\n")

cat("##", tmp_print, "##", "\n")

print(tmp_get$chosen_transform)

cat("_____________________________________")

cat("\n")

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

}

```

</details>

<details>

<summary>Show the complete `bestNormalize` results</summary>

```{r}

#| echo: false

#| eval: true

for (i in md_to_tranform) {

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

tmp_print <- c("Results of bestNormalize for", i)

cat("\n")

cat("##", tmp_print, "##", "\n")

print(tmp_get)

cat("_____________________________________")

cat("\n")

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

}

```

</details>

Great, now we can add the normalized transformed data back to our `mctoolsr` metadata file.

```{r}

#| code-fold: true

ssu18_select_mc_norm <- ssu18_select_mc

for (i in md_to_tranform) {

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

tmp_new_data <- tmp_get$x.t

ssu18_select_mc_norm$map_loaded[[i]] <- tmp_new_data

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

}

```

And rerun the Shapiro Tests.

```{r}

#| code-fold: true

temp_md_norm <- ssu18_select_mc_norm$map_loaded

temp_md_norm[,1:9] <- NULL

shap_results_norm <- NULL

for (i in colnames(temp_md_norm)) {

tmp_shap <- shapiro.test(temp_md_norm[[i]])

tmp_p <- round(tmp_shap$p.value, digits = 5)

tmp_res <- eval(isTRUE(tmp_shap$p.value < 0.05))

shap_results_norm <- rbind(shap_results_norm, data.frame(i, tmp_p, tmp_res))

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

}

colnames(shap_results_norm) <- c("parameter", "p-value", "tranform")

shap_results_norm

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

```

And check if there are any parameters that are still significant for the normality test.

```{r}

#| eval: true

shap_results$parameter[shap_results_norm$tranform == TRUE]

```

Ok. Looks like `bestNormalize` was unable to find a suitable transformation for `Al` and `Fe`. This is likely because there is very little variation in these metadata and/or there are too few significant digits.

## Normalized Metadata

Finally, here is a new summary table that includes all of the normalized data.

</br>

```{r}

#| echo: false

#| eval: true

metad_norm <- ssu18_select_mc_norm$map_loaded

metad_norm <- metad_norm %>% tibble::rownames_to_column("Sample_ID")

metad_norm <- metad_norm[order(match(metad_norm$Sample_ID, metad$Sample_ID)), ]

metad_norm[,2:9] <- NULL

seq_table <- metad_norm

write.table(metad_norm, "files/metadata/tables/ssu18_normalized_metadata.txt",

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

```

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

```{r}

#| echo: false

#| eval: true

seq_table <- seq_table %>% dplyr::mutate(across(.cols = c(3:5, 7:9, 14:15,19:20,

24,26,30,34:38,40, 42:43,

49, 52:53,55,57,60:62),

round, digits = 3))

reactable(seq_table,

defaultColDef = colDef(headerStyle = list(fontSize = "0.8em"),

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

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

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

footerStyle = list(fontWeight = "bold")

),

columns = list(

Sample_ID = colDef(name = "Sample ID",

sticky = "left",

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

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

align = "left",

minWidth = 150),

raw_reads = colDef(footer = function(values) sprintf("%.0f", sum(values)))

),

searchable = TRUE, defaultPageSize = 5,

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

showPageSizeOptions = TRUE, highlight = TRUE,

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

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

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

rm(seq_table)

```

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

```{r}

#| echo: false

#| eval: false

#| comment: "Code to run individual transformations for each metadata and generate histograms"

##for (i in md_to_tranform) {

## tmp_data <- ssu18_select_mc$map_loaded

## tmp_no_transform_obj <- no_transform(tmp_data[[i]])

## tmp_arcsinh_obj <- arcsinh_x(tmp_data[[i]])

## #tmp_boxcox_obj <- boxcox(tmp_data[[i]])

## tmp_yeojohnson_obj <- yeojohnson(tmp_data[[i]])

## tmp_orderNorm_obj <- orderNorm(tmp_data[[i]])

## tmp_log_obj <- log_x(tmp_data[[i]])

## tmp_sqrt_x_obj <- sqrt_x(tmp_data[[i]])

## tmp_exp_x_obj <- exp_x(tmp_data[[i]])

##

## tmp_plot_no_transform <- hist(tmp_no_transform_obj$x.t, main = i, breaks = 6, xlab = "no_transform")

## tmp_plot_name <- purrr::map_chr(i, ~ paste0(., "_hist_plot_no_transform"))

## assign(tmp_plot_name, tmp_plot_no_transform)

## print(tmp_plot_name)

##

## tmp_plot_arcsinh <- hist(tmp_arcsinh_obj$x.t, main = i, breaks = 6, xlab = "arcsinh")

## tmp_plot_name <- purrr::map_chr(i, ~ paste0(., "_hist_plot_arcsinh"))

## assign(tmp_plot_name, tmp_plot_arcsinh)

## print(tmp_plot_name)

## # cannot run because of Al and Fe

## #tmp_plot_boxcox <- hist(tmp_boxcox_obj$x.t, main = i, breaks = 6, xlab = "boxcox")

## #tmp_plot_name <- purrr::map_chr(i, ~ paste0(., "_hist_plot_boxcox"))

## #assign(tmp_plot_name, tmp_plot_boxcox)

## #print(tmp_plot_name)

##

## tmp_plot_yeojohnson <- hist(tmp_yeojohnson_obj$x.t, main = i, breaks = 6, xlab = "yeojohnson")

## tmp_plot_name <- purrr::map_chr(i, ~ paste0(., "_hist_plot_yeojohnson"))

## assign(tmp_plot_name, tmp_plot_yeojohnson)

## print(tmp_plot_name)

##

## tmp_plot_orderNorm <- hist(tmp_orderNorm_obj$x.t, main = i, breaks = 6, xlab = "orderNorm")

## tmp_plot_name <- purrr::map_chr(i, ~ paste0(., "_hist_plot_orderNorm"))

## assign(tmp_plot_name, tmp_plot_orderNorm)

## print(tmp_plot_name)

##

## tmp_plot_log <- hist(tmp_log_obj$x.t, main = i, breaks = 6, xlab = "log")

## tmp_plot_name <- purrr::map_chr(i, ~ paste0(., "_hist_plot_log"))

## assign(tmp_plot_name, tmp_plot_log)

## print(tmp_plot_name)

##

## tmp_plot_sqrt <- hist(tmp_sqrt_x_obj$x.t, main = i, breaks = 6, xlab = "sqrt")

## tmp_plot_name <- purrr::map_chr(i, ~ paste0(., "_hist_plot_sqrt"))

## assign(tmp_plot_name, tmp_plot_sqrt)

## print(tmp_plot_name)

##

## tmp_plot_exp <- hist(tmp_exp_x_obj$x.t, main = i, breaks = 6, xlab = "exp")

## tmp_plot_name <- purrr::map_chr(i, ~ paste0(., "_hist_plot_exp"))

## assign(tmp_plot_name, tmp_plot_exp)

## print(tmp_plot_name)

##

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

##}

```

## Autocorrelation Tests

Next, we test the metadata for autocorrelations. Do we do this on the original data or the transformed data? No idea, so let's do both.

### Split Metadata

We need to split the data into different groups.

A) Environmental and edaphic properties

B) Microbial functional responses

C) Temperature adaptation properties

We first create lists of metadata parameters.

```{r}

#| code-fold: true

div <- c("PLOT", "TREAT", "TREAT_T", "PAIR", "Observed", "Shannon_exp",

"InvSimpson", "ATAP")

edaphic <- c("AST", "H2O", "N", "P", "Al",

"Ca", "Fe", "K", "Mg", "Mn", "Na", "TEB", "ECEC", "pH",

"NH4", "NO3", "PO4", "DOC", "DON", "DOCN")

soil_funct <- c("micC", "micN", "micP", "micCN", "micCP", "micNP",

"AG_ase", "BG_ase", "BP_ase", "CE_ase", "P_ase", "N_ase",

"S_ase", "XY_ase", "LP_ase", "PX_ase", "CO2",

"enzCN", "enzCP", "enzNP")

temp_adapt <- c("AG_Q10", "BG_Q10", "BP_Q10", "CE_Q10", "P_Q10", "N_Q10",

"S_Q10", "XY_Q10", "LP_Q10", "PX_Q10", "CUEcn", "CUEcp",

"NUE","PUE", "Tmin", "SI")

md_groups <- c("edaphic", "soil_funct", "temp_adapt")

# NoT uSed: minPO4, minNH4, minNO3, minTIN

```

And then use the lists to split the data sets by metadata group. Here, we do this for the original metadata and the metadata after normalization.

```{r}

#| code-fold: true

select_md <- c("ssu18_select_mc", "ssu18_select_mc_norm")

for (i in select_md) {

#tmp_get <- get(purrr::map_chr(i, ~ paste0(i, "_mc")))

tmp_get <- get(i)

tmp_md_all <- tmp_get$map_loaded

tmp_div <- tmp_md_all %>% dplyr::select(all_of(div))

tmp_div <- tmp_div %>% tibble::rownames_to_column("SampleID")

## edaphic

tmp_sub_edaphic <- tmp_md_all %>% dplyr::select(all_of(edaphic))

tmp_sub_edaphic <- tmp_sub_edaphic %>% tibble::rownames_to_column("SampleID")

tmp_edaphic <- dplyr::left_join(tmp_div, tmp_sub_edaphic, by = "SampleID")

tmp_edaphic <- tmp_edaphic %>% tibble::column_to_rownames("SampleID")

## soil_funct

tmp_sub_soil_funct <- tmp_md_all %>% dplyr::select(all_of(soil_funct))

tmp_sub_soil_funct <- tmp_sub_soil_funct %>% tibble::rownames_to_column("SampleID")

tmp_soil_funct <- dplyr::left_join(tmp_div, tmp_sub_soil_funct, by = "SampleID")

tmp_soil_funct <- tmp_soil_funct %>% tibble::column_to_rownames("SampleID")

## temp_adapt

tmp_sub_temp_adapt <- tmp_md_all %>% dplyr::select(all_of(temp_adapt))

tmp_sub_temp_adapt <- tmp_sub_temp_adapt %>% tibble::rownames_to_column("SampleID")

tmp_temp_adapt <- dplyr::left_join(tmp_div, tmp_sub_temp_adapt, by = "SampleID")

tmp_temp_adapt <- tmp_temp_adapt %>% tibble::column_to_rownames("SampleID")

## combine

tmp_list <- list(data_loaded = ssu18_select_mc$data_loaded,

map_loaded = ssu18_select_mc$map_loaded,

taxonomy_loaded = ssu18_select_mc$taxonomy_loaded,

edaphic = tmp_edaphic,

soil_funct = tmp_soil_funct,

temp_adapt = tmp_temp_adapt)

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

print(tmp_name)

assign(tmp_name, tmp_list)

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

}

```

### Generate Autocorrelation Plots

A little housekeeping to get rid of parameters we don't need (e.g., plot number, pair, etc.).

```{r}

#| code-fold: true

edaphic_cor <- ssu18_select_mc_split$edaphic

edaphic_cor[,1:8] <- NULL

edaphic_norm_cor <- ssu18_select_mc_norm_split$edaphic

edaphic_norm_cor[,1:8] <- NULL

soil_funct_cor <- ssu18_select_mc_split$soil_funct

soil_funct_cor[,1:8] <- NULL

soil_funct_norm_cor <- ssu18_select_mc_norm_split$soil_funct

soil_funct_norm_cor[,1:8] <- NULL

temp_adapt_cor <- ssu18_select_mc_split$temp_adapt

temp_adapt_cor[,1:8] <- NULL

temp_adapt_norm_cor <- ssu18_select_mc_norm_split$temp_adapt

temp_adapt_norm_cor[,1:8] <- NULL

```

And finally the code to create the plots.

```{r}

#| code-fold: true

for (i in objects(pattern = "_cor$")) {

tmp_get <- get(i)

tmp_cormat <- round(cor(tmp_get), 2)

tmp_melted_cormat <- reshape2::melt(tmp_cormat)

tmp_get_lower_tri <- function(tmp_cormat){

tmp_cormat[upper.tri(tmp_cormat)] <- NA

return(tmp_cormat)

}

# Get upper triangle of the correlation matrix

tmp_get_upper_tri <- function(tmp_cormat){

tmp_cormat[lower.tri(tmp_cormat)] <- NA

return(tmp_cormat)

}

tmp_upper_tri <- tmp_get_upper_tri(tmp_cormat)

tmp_melted_cormat <- reshape2::melt(tmp_upper_tri, na.rm = TRUE)

ggplot(data = tmp_melted_cormat, aes(x = Var1, y = Var2, fill = value)) +

geom_tile()

tmp_ggheatmap <- ggplot(data = tmp_melted_cormat, aes(Var2, Var1, fill = value)) +

geom_tile(color = "white") +

scale_fill_gradient2(low = "blue", high = "red", mid = "white",

midpoint = 0, limit = c(-1,1), space = "Lab",

name="Pearson\nCorrelation") +

theme_minimal() +

theme(axis.text.x = element_text(angle = 45, vjust = 1,

size = 7, hjust = 1),

axis.text.y = element_text(vjust = 1, size = 7, hjust = 1)) +

coord_fixed() +

geom_text(aes(Var2, Var1, label = value), color = "black", size = 1.75) +

theme(

axis.title.x = element_blank(),

axis.title.y = element_blank(),

panel.grid.major = element_blank(),

panel.border = element_blank(),

panel.background = element_blank(),

axis.ticks = element_blank(),

legend.justification = c(1, 0),

legend.position = c(0.6, 0.7),

legend.direction = "horizontal") +

guides(fill = guide_colorbar(barwidth = 7, barheight = 1,

title.position = "top",

title.hjust = 0.5))

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

assign(tmp_name, tmp_ggheatmap)

print(tmp_name)

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

}

objects(pattern = "_ggheatmap")

```

## Autocorrelation Plots

<br/>

::: {.column-body}

::: {.panel-tabset}

#### Edaphic properties

```{r}

#| echo: false

#| eval: true

#| warning: false

#| fig-height: 4.5

edaphic_cor_ggheatmap

```

<small>`r caption_fig_ssu("ssu_autocor_edaph_plots")`</small>

#### Edaphic properties (normalized)

```{r}

#| echo: false

#| eval: true

#| warning: false

#| fig-height: 4.5

edaphic_norm_cor_ggheatmap

```

<small>`r caption_fig_ssu("ssu_autocor_edaph_norm_plots")`</small>

:::

::: {.panel-tabset}

#### Functional responses

```{r}

#| echo: false

#| eval: true

#| warning: false

#| fig-height: 4.5

soil_funct_cor_ggheatmap

```

<small>`r caption_fig_ssu("ssu_autocor_funct_plots")`</small>

#### Functional responses (normalized)

```{r}

#| echo: false

#| eval: true

#| warning: false

#| fig-height: 4.5

soil_funct_norm_cor_ggheatmap

```

<small>`r caption_fig_ssu("ssu_autocor_funct_norm_plots")`</small>

:::

::: {.panel-tabset}

#### Temperature adaptation

```{r}

#| echo: false

#| eval: true

#| warning: false

#| fig-height: 4.5

temp_adapt_cor_ggheatmap

```

<small>`r caption_fig_ssu("ssu_autocor_temp_plots")`</small>

#### Temperature adaptation (normalized)

```{r}

#| echo: false

#| eval: true

#| warning: false

#| fig-height: 4.5

temp_adapt_norm_cor_ggheatmap

```

<small>`r caption_fig_ssu("ssu_autocor_temp_norm_plots")`</small>

:::

:::

```{r}

#| echo: false