---
# System prepended metadata

title: Diversity biz (26-May-2021)
tags: [Aquaponics]

---

---
tags: Aquaponics
title: Diversity biz (26-May-2021)
---

# Diversity biz
# Summary

---
> **Nothing significantly different based on estimated shannon diversity at the family-level taxonomic classifications we have (that's the lowest taxonomic resolution we have).**

---

[toc]

# Setting up

```R
library(phyloseq)
library(DivNet)
library(tidyverse)
library(magrittr)

theme_set(theme_bw)

# helper functions for mass running things
make_div_obj <- function(phy_obj, wanted_formula, cores = 1) {

    # running divnet
    div_obj <- phy_obj %>% divnet(formula = wanted_formula, ncores = cores)

    return(div_obj)

}

run_div_test <- function(phy_obj, div_obj, target_ref = NULL) {

    # getting covariate
    target_covariate <- attributes(div_obj$X)$contrasts %>% names

    # getting single estimate and se for each treatment to give to betta function
        # getting sample info tab alone
    sample_tab <- phy_obj %>% sample_data %>% data.frame
        # getting unique treatment types
    unique_treatment_types <- sample_tab %>% pull(target_covariate) %>% unique %>% as.character

    unique_samples_needed <- c()

    for ( type in unique_treatment_types ) {

        sample_to_add <- sample_tab %>% filter(!!as.symbol(target_covariate) == !!type) %>% head(1) %>% row.names()

        unique_samples_needed <- c(unique_samples_needed, sample_to_add)
    }

        # subsetting divnet object to one per unique type
    sub_div_obj_shannon <- div_obj$shannon[unique_samples_needed]
    sub_div_obj_shannon_variance <- div_obj$`shannon-variance`[unique_samples_needed]

        # getting estimates and standard errors for one entry for each type
    estimates <- c()
    ses <- c()

    for ( i in 1:length(unique_samples_needed) ) {
        estimates <- c(estimates, unlist(sub_div_obj_shannon[i])[1] %>% as.numeric())
        ses <- c(ses, sqrt(sub_div_obj_shannon_variance[i]))
    }


    # making subset sample tab
    sub_sample_tab <- sample_tab[unique_samples_needed, ]

    # setting up matrix
    predictors_ <- sub_sample_tab[, names(sub_sample_tab) %in% target_covariate]

        # there is no built-in way to change the baseline in the contrasts, so modifying the factor levels as suggested on divnet issue here: https://github.com/adw96/DivNet/issues/37
    if ( is.null(target_ref) ) {
        X <- model.matrix(~predictors_)
        baseline <- predictors_ %>% levels %>% head(1)
    } else {
        predictors_ <- predictors_ %>% relevel(ref = target_ref)
        baseline <- predictors_ %>% levels %>% head(1)
        X <- model.matrix(~predictors_)
    }



    res <- betta(estimates, ses, X)

    cat("With '", baseline, "' as the baseline:\n\n", sep = "")

    tab <- res$table %>% data.frame

    # p-values are calculated as: 2*(1-pnorm(abs(Estimate/SE))) (see betta function)
    # I think getting 0s is just due to how many digits R is working with. It seems
    # pnorm(5.32672)  = 0.9999999 and pnorm(5.32673) = 1. So if our abs(Estimate/SE)
    # is greater than 5.32672, we’re going to get a p-value of 0. So setting 0's to
    # < 1e-8 for what is printed out

    new_p_vec <- vector()

    for ( value in tab$p.values ) {
        if ( value == 0 ) {
            new_p_vec <- c(new_p_vec, "< 1e-8")
        } else {
            new_p_vec <- c(new_p_vec, value)
        }
    }

    tab$p.values <- new_p_vec

    print(tab)

    return(res)

}

plot_DivNet_diversity <- function(phy_obj, div_obj, rank = "Family") {

    # getting covariate
    target_covariate <- attributes(div_obj$X)$contrasts %>% names

    # getting sample table from phyloseq object
    sample_tab <- phy_obj %>% sample_data %>% data.frame

    # getting unique treatment types
    unique_treatment_types <- sample_tab %>% pull(target_covariate) %>% unique %>% as.character

    unique_samples_needed <- c()

    for ( type in unique_treatment_types ) {

        sample_to_add <- sample_tab %>% filter(!!as.symbol(target_covariate) == type) %>% head(1) %>% row.names()

        unique_samples_needed <- c(unique_samples_needed, sample_to_add)
    }

    sub_div_obj_shannon <- div_obj$shannon[c(unique_samples_needed)]

    div_est_vec <- c()
    ci_vec <- c()

    for ( i in 1:length(unique_samples_needed) ) {
        div_est_vec <- c(div_est_vec, unlist(sub_div_obj_shannon[i])[1] %>% as.numeric())
        ci_vec <- c(ci_vec, unlist(sub_div_obj_shannon[i])[2] %>% as.numeric() * 2)
    }

    sub_plot_tab <- data.frame(target_covariate = unique_treatment_types, shannon = div_est_vec, ci = ci_vec)

    # getting max value for plot y-axis
    # max_y <- max(sub_plot_tab$shannon + sub_plot_tab$ci) %>% ceiling

    y_axis_label = paste0("Shannon diversity estimate\n(", rank, " level)")
    plot <- ggplot(sub_plot_tab, aes(x = target_covariate, y = shannon, color = target_covariate)) +
        geom_errorbar(aes(ymin = shannon - ci, ymax = shannon + ci), width = 0) +
        geom_point() + theme_bw() +
        xlab(target_covariate) + ylab(y_axis_label) +
        theme(legend.position = "none") # + coord_cartesian(ylim = c(0,max_y))

    return(plot)

}
```

```
# reading in our data (excluding EM1 and Lambda)
otu_counts_tab <- read.table("combined_otu_counts_tab.txt", header=TRUE, row.names=1, sep="\t")[,c(1:30)]
taxonomy_tab <- read.table("combined_taxonomy_tab.txt", header=TRUE, row.names=1, sep="\t", na.strings="<NA>")
sample_info_tab <- read.table("sample_info_tab.txt", header=T, check.names=F, row.names=1, sep="\t")[c(1:30),]

# making phyloseq object
counts_phy <- otu_table(otu_counts_tab, taxa_are_rows=TRUE)
tax_phy <- tax_table(as.matrix(taxonomy_tab))
sample_phy <- sample_data(sample_info_tab)
aqua_phy <- phyloseq(counts_phy, tax_phy, sample_phy)

# grouping at family rank
aqua_family_phy <- tax_glom(aqua_phy, taxrank = "Family")
```

# Based on timepoint

```
time_point_div_obj <- make_div_obj(aqua_family_phy, ~Time_point, cores = 4)

time_point_T1_base_div_res <- run_div_test(aqua_family_phy, time_point_div_obj)
```

## Not Sig (p = 0.307)

```
With 'T1' as the baseline:

                Estimates Standard.Errors p.values
(Intercept)    3.07915048       0.1057571   < 1e-8
predictors_T2  0.01151054       0.1735481    0.947
predictors_T3 -0.23486316       0.2413842    0.331
```


```
time_point_T2_base_div_res <- run_div_test(aqua_family_phy, time_point_div_obj, target_ref = "T2")
```

```
With 'T2' as the baseline:

                Estimates Standard.Errors p.values
(Intercept)    3.09066102       0.1057571   < 1e-8
predictors_T1 -0.01151054       0.1600364    0.943
predictors_T3 -0.24637370       0.2413842    0.307
```

## Plotting
```
plot_DivNet_diversity(aqua_family_phy, time_point_div_obj)
```

<a href="https://i.imgur.com/1HVkXIv.png"><img src="https://i.imgur.com/1HVkXIv.png"></a>

---

# Based on lighting
```
lighting_div_obj <- make_div_obj(aqua_family_phy, ~Lighting, cores = 4)
lighting_div_res <- run_div_test(aqua_family_phy, lighting_div_obj)
```

## Not Sig (p = 0.208)

```
With 'Artificial' as the baseline:

                   Estimates Standard.Errors p.values
(Intercept)        2.9169832       0.1190661   < 1e-8
predictors_Natural 0.1772417       0.1407935    0.208
```


## Plotting
```
plot_DivNet_diversity(aqua_family_phy, lighting_div_obj)
```

<a href="https://i.imgur.com/C4LXsD9.png"><img src="https://i.imgur.com/C4LXsD9.png"></a>


---

# Based on plants

```
plants_div_obj <- make_div_obj(aqua_family_phy, ~Plants, cores = 4)
plants_div_res <- run_div_test(aqua_phy, plants_div_obj)
```

## Tomatos lower diversity than Peas (0.068)

```
With 'Basil' as the baseline:

                     Estimates Standard.Errors p.values
(Intercept)          3.0236533       0.1563061   < 1e-8
predictors_Peas      0.1597980       0.2672525     0.55
predictors_Tomatoes -0.3516465       0.2807357     0.21
```


```
plants_div_res_peas_baseline <- run_div_test(aqua_phy, plants_div_obj, target_ref = "Peas")
```

```
With 'Peas' as the baseline:

                     Estimates Standard.Errors p.values
(Intercept)          3.1834513       0.1563061   < 1e-8
predictors_Basil    -0.1597980       0.2649870    0.546
predictors_Tomatoes -0.5114445       0.2807357    0.068
```

## Plotting
```
plot_DivNet_diversity(aqua_family_phy, plants_div_obj)
```

<a href="https://i.imgur.com/944FNda.png"><img src="https://i.imgur.com/944FNda.png"></a>

---

# Based on water

```
water_div_obj <- make_div_obj(aqua_family_phy, ~Water_circulation, cores = 4)
water_div_res <- run_div_test(aqua_family_phy, water_div_obj)
```

## Not Sig (0.26)

```With 'Pumps' as the baseline:

                  Estimates Standard.Errors p.values
(Intercept)       3.1727239       0.1948949   < 1e-8
predictors_Wicks -0.3268373       0.2880062    0.256
```

## Plotting
```
plot_DivNet_diversity(aqua_family_phy, water_div_obj)
```

<a href="https://i.imgur.com/VZezc88.png"><img src="https://i.imgur.com/VZezc88.png"></a>

---

# Based on nitrifying components

```
N_div_obj <- make_div_obj(aqua_family_phy, ~Nitrifying_components, cores = 4)
N_div_res <- run_div_test(aqua_family_phy, N_div_obj)
```

## Not Sig (0.172)

```
With 'Collected_gravel' as the baseline:

                      Estimates Standard.Errors p.values
(Intercept)            2.775672       0.2106208   < 1e-8
predictors_Commercial  0.347640       0.2547854    0.172
```

## Plotting
```
plot_DivNet_diversity(aqua_family_phy, N_div_obj)
```

<a href="https://i.imgur.com/QXejwuE.png"><img src="https://i.imgur.com/QXejwuE.png"></a>

---

# Splitting plants

```
aqua_peas_family_phy <- subset_samples(aqua_family_phy, Plants == "Peas")
aqua_basil_family_phy <- subset_samples(aqua_family_phy, Plants == "Basil")
aqua_tomatoes_family_phy <- subset_samples(aqua_family_phy, Plants == "Tomatoes")

```

## Peas

### Timepoint (not sig: p = 0.45)

```
peas_time_point_div_obj <- make_div_obj(aqua_peas_family_phy, ~Time_point, cores = 4)

peas_time_point_T1_base_div_res <- run_div_test(aqua_peas_family_phy, peas_time_point_div_obj)
```

```
With 'T1' as the baseline:

               Estimates Standard.Errors p.values
(Intercept)    3.0053968       0.2185220   < 1e-8
predictors_T2  0.2398329       0.3193968    0.453
predictors_T3 -0.2985799       0.4024406    0.458
```


```
peas_time_point_T2_base_div_res <- run_div_test(aqua_peas_family_phy, peas_time_point_div_obj, target_ref = "T2")
```

```
With 'T2' as the baseline:

               Estimates Standard.Errors p.values
(Intercept)    3.2452297       0.2185220   < 1e-8
predictors_T1 -0.2398329       0.4488045    0.593
predictors_T3 -0.5384127       0.4024406    0.181
```

### Plotting
```
plot_DivNet_diversity(aqua_peas_family_phy, peas_time_point_div_obj)
```

<a href="https://i.imgur.com/b5fyCwo.png"><img src="https://i.imgur.com/b5fyCwo.png"></a>


---

### Lighting (not sig: p = 0.555)

```
peas_lighting_div_obj <- make_div_obj(aqua_peas_family_phy, ~Lighting, cores = 4)
peas_lighting_div_res <- run_div_test(aqua_peas_family_phy, peas_lighting_div_obj)
```

```
With 'Artificial' as the baseline:

                   Estimates Standard.Errors p.values
(Intercept)        3.0981946       0.2239162   < 1e-8
predictors_Natural 0.1983297       0.3363078    0.555
```

### Plotting
```
plot_DivNet_diversity(aqua_peas_family_phy, peas_lighting_div_obj)
```

<a href="https://i.imgur.com/gNJwu5t.png"><img src="https://i.imgur.com/gNJwu5t.png"></a>

---

### Water (not sig: p = 0.304)

```
peas_water_div_obj <- make_div_obj(aqua_peas_family_phy, ~Water_circulation, cores = 4)
peas_water_div_res <- run_div_test(aqua_peas_family_phy, peas_water_div_obj)
```
```
With 'Pumps' as the baseline:

                  Estimates Standard.Errors p.values
(Intercept)       3.2475024       0.2249006   < 1e-8
predictors_Wicks -0.3565429       0.3470048    0.304
```


### Plotting
```
plot_DivNet_diversity(aqua_peas_family_phy, peas_water_div_obj)
```

<a href="https://i.imgur.com/dD34bvX.png"><img src="https://i.imgur.com/dD34bvX.png"></a>

---

### Nitrifying components (not sig: p = 0.159)

```
peas_N_div_obj <- make_div_obj(aqua_peas_family_phy, ~Nitrifying_components, cores = 4)
peas_N_div_res <- run_div_test(aqua_peas_family_phy, peas_N_div_obj)
```
```
With 'Collected_gravel' as the baseline:

                      Estimates Standard.Errors p.values
(Intercept)           2.8891087       0.2149975   < 1e-8
predictors_Commercial 0.3563273       0.2528028    0.159
```


### Plotting
```
plot_DivNet_diversity(aqua_peas_family_phy, peas_water_div_obj)
```

<a href="https://i.imgur.com/awDr1eQ.png"><img src="https://i.imgur.com/awDr1eQ.png"></a>

---

## Basil

### Timepoint (not sig: p = 0.239)

```
basil_time_point_div_obj <- make_div_obj(aqua_basil_family_phy, ~Time_point, cores = 4)

basil_time_point_T1_base_div_res <- run_div_test(aqua_basil_family_phy, basil_time_point_div_obj)
```

```
With 'T1' as the baseline:

               Estimates Standard.Errors p.values
(Intercept)   2.79147025       0.1281569   < 1e-8
predictors_T2 0.25188425       0.4396179    0.567
predictors_T3 0.05699726       0.1654831    0.731
```


```
basil_time_point_T2_base_div_res <- run_div_test(aqua_basil_family_phy, basil_time_point_div_obj, target_ref = "T2")
```

```
               Estimates Standard.Errors p.values
(Intercept)    3.0433545       0.1281569   < 1e-8
predictors_T1 -0.2518843       0.2282494     0.27
predictors_T3 -0.1948870       0.1654831    0.239
```

### Plotting
```
plot_DivNet_diversity(aqua_basil_family_phy, basil_time_point_div_obj)
```

<a href="https://i.imgur.com/1zpPvFg.png"><img src="https://i.imgur.com/1zpPvFg.png"></a>


---

### Lighting (not sig: p = 0.166)

```
basil_lighting_div_obj <- make_div_obj(aqua_basil_family_phy, ~Lighting, cores = 4)
basil_lighting_div_res <- run_div_test(aqua_basil_family_phy, basil_lighting_div_obj)
```

```
With 'Artificial' as the baseline:

                   Estimates Standard.Errors p.values
(Intercept)        2.7769250       0.1732599   < 1e-8
predictors_Natural 0.2877119       0.2078168    0.166
```

### Plotting
```
plot_DivNet_diversity(aqua_basil_family_phy, basil_lighting_div_obj)
```

<a href="https://i.imgur.com/Gy97pUP.png"><img src="https://i.imgur.com/Gy97pUP.png"></a>

---

### Water (not sig: p = 0.288)

```
basil_water_div_obj <- make_div_obj(aqua_basil_family_phy, ~Water_circulation, cores = 4)
basil_water_div_res <- run_div_test(aqua_basil_family_phy, basil_water_div_obj)
```

```
With 'Pumps' as the baseline:

                  Estimates Standard.Errors p.values
(Intercept)       3.0317425        0.115106   < 1e-8
predictors_Wicks -0.1771017        0.166659    0.288
```


### Plotting
```
plot_DivNet_diversity(aqua_basil_family_phy, basil_water_div_obj)
```

<a href="https://i.imgur.com/ZKDzqEt.png"><img src="https://i.imgur.com/ZKDzqEt.png"></a>

---

### Nitrifying components (not sig: p = 0.170)

```
basil_N_div_obj <- make_div_obj(aqua_basil_family_phy, ~Nitrifying_components, cores = 4)
basil_N_div_res <- run_div_test(aqua_basil_family_phy, basil_N_div_obj)
```
```
With 'Collected_gravel' as the baseline:

                      Estimates Standard.Errors p.values
(Intercept)           2.7015675       0.2525301   < 1e-8
predictors_Commercial 0.3786621       0.2756966     0.17
```


### Plotting
```
plot_DivNet_diversity(aqua_basil_family_phy, basil_water_div_obj)
```

<a href="https://i.imgur.com/RWbrAWV.png"><img src="https://i.imgur.com/RWbrAWV.png"></a>

---

## Tomatoes

### Timepoint (not sig: p = 0.146)

```
tomatoes_time_point_div_obj <- make_div_obj(aqua_tomatoes_family_phy, ~Time_point, cores = 4)

tomatoes_time_point_T1_base_div_res <- run_div_test(aqua_tomatoes_family_phy, tomatoes_time_point_div_obj)
```

```
With 'T1' as the baseline:

               Estimates Standard.Errors p.values
(Intercept)    2.9060773       0.2620180   < 1e-8
predictors_T2 -0.0650556       0.4363785    0.881
predictors_T3 -0.7539750       0.5183790    0.146
```


```
tomatoes_time_point_T2_base_div_res <- run_div_test(aqua_tomatoes_family_phy, tomatoes_time_point_div_obj, target_ref = "T2")
```

```
With 'T2' as the baseline:

               Estimates Standard.Errors p.values
(Intercept)    2.8410217        0.262018   < 1e-8
predictors_T1  0.0650556        0.422836    0.878
predictors_T3 -0.6889194        0.518379    0.184
```

### Plotting
```
plot_DivNet_diversity(aqua_tomatoes_family_phy, tomatoes_time_point_div_obj)
```

<a href="https://i.imgur.com/f3Qx8XU.png"><img src="https://i.imgur.com/f3Qx8XU.png"></a>


---

### Lighting (not sig: p = 0.329)

```
tomatoes_lighting_div_obj <- make_div_obj(aqua_tomatoes_family_phy, ~Lighting, cores = 4)
tomatoes_lighting_div_res <- run_div_test(aqua_tomatoes_family_phy, tomatoes_lighting_div_obj)
```

```
With 'Artificial' as the baseline:

                    Estimates Standard.Errors p.values
(Intercept)         2.7193033       0.1493807   < 1e-8
predictors_Natural -0.2234033       0.2288393    0.329
```

### Plotting
```
plot_DivNet_diversity(aqua_tomatoes_family_phy, tomatoes_lighting_div_obj)
```

<a href="https://i.imgur.com/X8tUhGx.png"><img src="https://i.imgur.com/X8tUhGx.png"></a>
---

### Water (not sig: p = 0.878)

```
tomatoes_water_div_obj <- make_div_obj(aqua_tomatoes_family_phy, ~Water_circulation, cores = 4)
tomatoes_water_div_res <- run_div_test(aqua_tomatoes_family_phy, tomatoes_water_div_obj)
```
```
With 'Pumps' as the baseline:

                  Estimates Standard.Errors p.values
(Intercept)      2.68439913       0.2048484   < 1e-8
predictors_Wicks 0.06047993       0.3926903    0.878
```


### Plotting
```
plot_DivNet_diversity(aqua_tomatoes_family_phy, tomatoes_water_div_obj)
```

<a href="https://i.imgur.com/XyrNsOE.png"><img src="https://i.imgur.com/XyrNsOE.png"></a>
---

### Nitrifying components (not sig: p = 0.293)

```
tomatoes_N_div_obj <- make_div_obj(aqua_tomatoes_family_phy, ~Nitrifying_components, cores = 4)
tomatoes_N_div_res <- run_div_test(aqua_tomatoes_family_phy, tomatoes_N_div_obj)
```
```
With 'Collected_gravel' as the baseline:

                      Estimates Standard.Errors p.values
(Intercept)            2.493579       0.1500111   < 1e-8
predictors_Commercial  0.225985       0.2147059    0.293
```


### Plotting
```
plot_DivNet_diversity(aqua_tomatoes_family_phy, tomatoes_water_div_obj)
```

<a href="https://i.imgur.com/KXqOmYY.png"><img src="https://i.imgur.com/KXqOmYY.png"></a>


---
