#### Purpose of this analysis The goal of this analysis is to check the sensitivity of the biodiversity-ecosystem function relationships to rare OTUs. #### More specifically: We remove OTUs that have X number of sequences throughout the entire dataset and then check the relationship with diversity versus community-wide and per-capita heterotrophic production. This we call "singletons", "doubletons", "5-tons", "10-tons", etc all the way up to "300-tons". For context, removing 10-tons will be removing any OTUs that have a count of less than 10 throughout the entire dataset. ##### If you have any questions, please be welcome to email the corresponding author at marschmi@umich.edu or tweet her at [micro_marian](https://twitter.com/micro_marian?lang=en). ```{r setup, include=FALSE} # For width of code chunks and scroll bar options(width=250) knitr::opts_chunk$set(eval = TRUE, echo = TRUE, cache = TRUE, include = TRUE, warning = FALSE, collapse = FALSE, message = FALSE, engine = "R", # Chunks will always have R code, unless noted error = TRUE, fig.path="OTU_Removal_Analysis_Figs/", # Set the figure options fig.align = "center", fig.width = 7, fig.height = 7) ``` # Load Libraries ```{r load-libraries, message = FALSE, warning = FALSE} library(ggplot2) library(devtools) library(phyloseq) library(kableExtra) library(tidyr) library(dplyr) library(cowplot) library(forcats) library(picante) # Will also include ape and vegan library(car) # For residual analysis library(sandwich) # for vcovHC function in post-hoc test library(MASS) # For studres in plot_residuals function library(boot) # For cross validation library(DT) # Pretty HTML Table Output source("../code/Muskegon_functions.R") source("../code/set_colors.R") # Set the ggplot theme theme_set(theme_cowplot()) ``` # Prepare the data ## Load Data ```{r load-data, eval = TRUE, warning = FALSE} # Loads a phyloseq object named otu_merged_musk_pruned) load("../data/otu_merged_musk_pruned.RData") # The name of the phyloseq object is: otu_merged_musk_pruned # Productivity measurements are reliable only up to 1 decimal df1 <- sample_data(otu_merged_musk_pruned) %>% dplyr::mutate(tot_bacprod = round(tot_bacprod, digits = 1), SD_tot_bacprod = round(SD_tot_bacprod, digits = 1), frac_bacprod = round(frac_bacprod, digits = 1), SD_frac_bacprod = round(SD_frac_bacprod, digits = 1), fraction_bac_abund = as.numeric(fraction_bac_abund), fracprod_per_cell = frac_bacprod/(1000*fraction_bac_abund), fracprod_per_cell_noinf = ifelse(fracprod_per_cell == Inf, NA, fracprod_per_cell)) %>% dplyr::select(norep_filter_name, lakesite, limnion, fraction, year, season, tot_bacprod, SD_tot_bacprod, frac_bacprod, SD_frac_bacprod, fraction_bac_abund, fracprod_per_cell, fracprod_per_cell_noinf) row.names(df1) = df1$norep_filter_name # Add new sample data back into phyloseq object sample_data(otu_merged_musk_pruned) <- df1 # Remove MOTHJ715 and MBRHP715 because of low sequencing depth otu_merged_musk_pruned_noMOTHJ715_MBRHP715 <- subset_samples(otu_merged_musk_pruned, norep_filter_name != "MOTHJ715" & norep_filter_name != "MBRHP715") # Subset only the surface samples for the current study!! musk_surface <- subset_samples(otu_merged_musk_pruned_noMOTHJ715_MBRHP715, limnion == "Top" & year == "2015" & fraction %in% c("WholePart","WholeFree")) # Surface samples, 2015, and WholePart/WholeFree samples only! musk_surface_pruned <- prune_taxa(taxa_sums(musk_surface) > 0, musk_surface) # Remove tree notree_musk_surface_pruned <- phyloseq(tax_table(musk_surface_pruned), otu_table(musk_surface_pruned), sample_data(musk_surface_pruned)) # Remove singletons! musk_surface_pruned_rm1 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 1, notree_musk_surface_pruned) # Remove the tree for less computationally intensive steps notree_musk_surface_pruned_rm1 <- phyloseq(tax_table(musk_surface_pruned_rm1), otu_table(musk_surface_pruned_rm1), sample_data(musk_surface_pruned_rm1)) # If taxa with 2 counts are removed prune_taxa(taxa_sums(notree_musk_surface_pruned_rm1) > 2, notree_musk_surface_pruned_rm1) # If taxa with 5 counts are removed notree_musk_surface_pruned_rm5 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 5, notree_musk_surface_pruned) # If taxa with 10 counts are removed notree_musk_surface_pruned_rm10 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 10, notree_musk_surface_pruned) # If taxa with 20 counts are removed notree_musk_surface_pruned_rm20 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 20, notree_musk_surface_pruned) # If taxa with 30 counts are removed notree_musk_surface_pruned_rm30 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 30, notree_musk_surface_pruned) # If taxa with 60 counts are removed notree_musk_surface_pruned_rm60 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 60, notree_musk_surface_pruned) # If taxa with 90 counts are removed notree_musk_surface_pruned_rm90 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 90, notree_musk_surface_pruned) # If taxa with 150 counts are removed notree_musk_surface_pruned_rm150 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 150, notree_musk_surface_pruned) # If taxa with 300 counts are removed notree_musk_surface_pruned_rm225 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 225, notree_musk_surface_pruned) # If taxa with 300 counts are removed notree_musk_surface_pruned_rm300 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 300, notree_musk_surface_pruned) ``` ## Calculate Diversity ```{r calc-div} set.seed(777) ################## Remove singltons alpha_rm1 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm1) otu_alphadiv_rm1 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm1, richness_df = alpha_rm1$Richness, simpson_df = alpha_rm1$Inverse_Simpson, exp_shannon_df = alpha_rm1$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "1-tons") ################## Remove 5-tons alpha_rm5 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm5) otu_alphadiv_rm5 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm5, richness_df = alpha_rm5$Richness, simpson_df = alpha_rm5$Inverse_Simpson, exp_shannon_df = alpha_rm5$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "5-tons") ################## Remove 10-tons alpha_rm10 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm10) otu_alphadiv_rm10 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm10, richness_df = alpha_rm10$Richness, simpson_df = alpha_rm10$Inverse_Simpson, exp_shannon_df = alpha_rm10$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "10-tons") ################## Remove 20-tons alpha_rm20 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm20) otu_alphadiv_rm20 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm20, richness_df = alpha_rm20$Richness, simpson_df = alpha_rm20$Inverse_Simpson, exp_shannon_df = alpha_rm20$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "20-tons") ################## Remove 30-tons alpha_rm30 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm30) otu_alphadiv_rm30 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm30, richness_df = alpha_rm30$Richness, simpson_df = alpha_rm30$Inverse_Simpson, exp_shannon_df = alpha_rm30$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "30-tons") ################## Remove 60-tons alpha_rm60 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm60) otu_alphadiv_rm60 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm60, richness_df = alpha_rm60$Richness, simpson_df = alpha_rm60$Inverse_Simpson, exp_shannon_df = alpha_rm60$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "60-tons") ################## Remove 90-tons alpha_rm90 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm90) otu_alphadiv_rm90 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm90, richness_df = alpha_rm90$Richness, simpson_df = alpha_rm90$Inverse_Simpson, exp_shannon_df = alpha_rm90$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "90-tons") ################## Remove 90-tons alpha_rm150 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm150) otu_alphadiv_rm150 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm150, richness_df = alpha_rm150$Richness, simpson_df = alpha_rm150$Inverse_Simpson, exp_shannon_df = alpha_rm150$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "150-tons") ################## Remove 300-tons alpha_rm225 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm225) otu_alphadiv_rm225 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm225, richness_df = alpha_rm225$Richness, simpson_df = alpha_rm225$Inverse_Simpson, exp_shannon_df = alpha_rm225$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "225-tons") ################## Remove 300-tons alpha_rm300 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm300) otu_alphadiv_rm300 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm300, richness_df = alpha_rm300$Richness, simpson_df = alpha_rm300$Inverse_Simpson, exp_shannon_df = alpha_rm300$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "300-tons") ``` ## Minimum Sequences Plot ```{r min-seqs-plots} min_seqs <- c(min(sample_sums(notree_musk_surface_pruned_rm1)) - 1, min(sample_sums(notree_musk_surface_pruned_rm5)) - 1, min(sample_sums(notree_musk_surface_pruned_rm10)) - 1, min(sample_sums(notree_musk_surface_pruned_rm20)) - 1, min(sample_sums(notree_musk_surface_pruned_rm30)) - 1, min(sample_sums(notree_musk_surface_pruned_rm60)) - 1, min(sample_sums(notree_musk_surface_pruned_rm90)) - 1, min(sample_sums(notree_musk_surface_pruned_rm150)) - 1, min(sample_sums(notree_musk_surface_pruned_rm225)) - 1, min(sample_sums(notree_musk_surface_pruned_rm300)) - 1) num_otus <- c(ncol(otu_table(notree_musk_surface_pruned_rm1)), ncol(otu_table(notree_musk_surface_pruned_rm5)), ncol(otu_table(notree_musk_surface_pruned_rm10)), ncol(otu_table(notree_musk_surface_pruned_rm20)), ncol(otu_table(notree_musk_surface_pruned_rm30)), ncol(otu_table(notree_musk_surface_pruned_rm60)), ncol(otu_table(notree_musk_surface_pruned_rm90)), ncol(otu_table(notree_musk_surface_pruned_rm150)), ncol(otu_table(notree_musk_surface_pruned_rm225)), ncol(otu_table(notree_musk_surface_pruned_rm300))) Removed <- c("1-tons","5-tons", "10-tons", "20-tons", "30-tons", "60-tons", "90-tons", "150-tons", "225-tons","300-tons") statz <- data.frame(cbind(as.numeric(min_seqs), as.numeric(num_otus), Removed)) %>% mutate(Removed = factor(Removed, levels = c("1-tons","5-tons", "10-tons", "20-tons", "30-tons", "60-tons", "90-tons", "150-tons", "225-tons","300-tons"))) %>% mutate(num_otus = as.numeric(num_otus)) p1 <- ggplot(statz, aes(x = Removed, y = min_seqs, fill = Removed)) + geom_bar(stat = "identity") + ylab("Minimum # of Sequences") + scale_y_continuous(expand = c(0,0), limits = c(0, 8000)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1), legend.position = "none") p2 <-ggplot(statz, aes(x = Removed, y = num_otus, fill = Removed)) + geom_bar(stat = "identity") + ylab("Richness") + scale_y_continuous(expand = c(0,0)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1), legend.position = c(0.8, 0.65), legend.title = element_blank()) plot_grid(p1, p2, align = "v", labels = c("A", "B"), nrow = 2, ncol =1) ``` ## 15 and 25-tons analysis ```{r} ### 15-tons analysis notree_musk_surface_pruned_rm15 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 15, notree_musk_surface_pruned) alpha_rm15 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm15) otu_alphadiv_rm15 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm15, richness_df = alpha_rm15$Richness, simpson_df = alpha_rm15$Inverse_Simpson, exp_shannon_df = alpha_rm15$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "15-tons") ### To test for per-capita production summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(otu_alphadiv_rm15, measure == "Richness" & fraction == "WholePart"))) summary(lm(log10(fracprod_per_cell_noinf) ~ mean, data = filter(otu_alphadiv_rm15, measure == "Inverse_Simpson" & fraction == "WholePart"))) ### 25-tons analysis notree_musk_surface_pruned_rm25 <- prune_taxa(taxa_sums(notree_musk_surface_pruned) > 25, notree_musk_surface_pruned) alpha_rm25 <- calc_alpha_diversity(physeq = notree_musk_surface_pruned_rm25) otu_alphadiv_rm25 <- calc_mean_alphadiv(physeq = notree_musk_surface_pruned_rm25, richness_df = alpha_rm25$Richness, simpson_df = alpha_rm25$Inverse_Simpson, exp_shannon_df = alpha_rm25$Shannon) %>% mutate(fraction = factor(fraction, levels = c("WholePart", "Particle", "WholeFree", "Free")), lakesite = factor(lakesite, levels = c("MOT", "MDP", "MBR", "MIN")), measure = factor(measure, levels = c("Richness", "Exponential_Shannon", "Inverse_Simpson")), Removed = "25-tons") summary(lm(frac_bacprod ~ mean, data = filter(otu_alphadiv_rm25, measure == "Richness" & fraction == "WholePart"))) summary(lm(frac_bacprod ~ mean, data = filter(otu_alphadiv_rm25, measure == "Inverse_Simpson" & fraction == "WholePart"))) ``` ## Combine all data ```{r all_divs} # Combine all div metrics all_divs <- bind_rows(otu_alphadiv_rm1, otu_alphadiv_rm5, otu_alphadiv_rm10, otu_alphadiv_rm20, otu_alphadiv_rm30, otu_alphadiv_rm60, otu_alphadiv_rm90, otu_alphadiv_rm150, otu_alphadiv_rm225, otu_alphadiv_rm300) %>% dplyr::filter(fraction %in% c("WholePart", "WholeFree") & year == "2015") %>% mutate(fraction = fct_recode(fraction, "Particle" = "WholePart", "Free" = "WholeFree")) %>% mutate(Removed = factor(Removed, levels = c("1-tons","5-tons", "10-tons", "20-tons", "30-tons", "60-tons", "90-tons", "150-tons", "225-tons","300-tons"))) ``` # Richness ```{r rich-plots1, fig.width=9, fig.height=4, dependson="all_divs"} ### PLOT ggplot(dplyr::filter(all_divs, measure == "Richness"), aes(y = mean, x = Removed, color = Removed, fill = Removed)) + geom_boxplot(alpha = 0.3, outlier.shape = NA) + geom_point(size = 3, position = position_jitter(w = 0.1)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(.~fraction) + ylab("Mean Richness") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1), axis.title.x = element_blank()) ggplot(dplyr::filter(all_divs, measure == "Richness"), aes(y = mean, x = fraction, color = Removed, fill = Removed)) + geom_boxplot(alpha = 0.3, outlier.shape = NA) + geom_point(size = 3, position = position_jitter(w = 0.1)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(.~Removed) + ylab("Mean Richness") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1), axis.title.x = element_blank()) ``` ```{r richness-plots, fig.width=10, fig.height=5, dependson="all_divs"} # Linear Model output richness_lm_results <- lm_fraction_output(dataframe = dplyr::filter(all_divs, measure == "Richness")) sig_rich_lms_df <- richness_lm_results %>% bind_rows() %>% mutate(diversity_metric = "Richness") %>% filter(pval < 0.05) # Display significant models in a dataframe datatable(sig_rich_lms_df, options = list(pageLength = 40)) ### Community-Wide Production vs Richness sig_rich_lms_comm <- sig_rich_lms_df %>% filter(test == "Community-Wide Production" & fraction == "Particle") %>% dplyr::select(Removed) %>% .$Removed ggplot(dplyr::filter(all_divs, measure == "Richness"), aes(y = frac_bacprod, x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Richness") + ylab("Community-Wide Production") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Richness" & fraction == "Particle" & Removed %in% sig_rich_lms_comm)) + scale_color_manual(values = tons_colors) +scale_fill_manual(values = tons_colors) + facet_grid(fraction~Removed, scales = "free") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)) ### Per-capita production vs Richness sig_rich_lms_percap <- sig_rich_lms_df %>% filter(test == "Per-Capita Production") %>% dplyr::select(Removed) %>% .$Removed ggplot(dplyr::filter(all_divs, measure == "Richness"), aes(y = log10(fracprod_per_cell_noinf), x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Richness") + ylab("log10(Per-Capita Production)") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Richness" & fraction == "Particle" & Removed %in% sig_rich_lms_percap)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(fraction~Removed, scales = "free") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)) ``` # Exp(Shannon Entropy) ```{r shannon-plots1, fig.width=9, fig.height=4, dependson="all_divs"} ### PLOT ggplot(dplyr::filter(all_divs, measure == "Exponential_Shannon"), aes(y = mean, x = Removed, color = Removed, fill = Removed)) + geom_boxplot(alpha = 0.3, outlier.shape = NA) + geom_point(size = 3, position = position_jitter(w = 0.1)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(.~fraction) + ylab("Mean Exp(Shannon Entropy)") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1), axis.title.x = element_blank()) ggplot(dplyr::filter(all_divs, measure == "Exponential_Shannon"), aes(y = mean, x = fraction, color = Removed, fill = Removed)) + geom_boxplot(alpha = 0.3, outlier.shape = NA) + geom_point(size = 3, position = position_jitter(w = 0.1)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(.~Removed) + ylab("Mean Exp(Shannon Entropy)") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1), axis.title.x = element_blank()) ``` ```{r shannon-plots, fig.width=10, fig.height=5, dependson="all_divs"} # Linear Model output shannon_lm_results <- lm_fraction_output(dataframe = dplyr::filter(all_divs, measure == "Exponential_Shannon")) sig_shannon_lms_df <- shannon_lm_results %>% bind_rows() %>% mutate(diversity_metric = "Exponential_Shannon") %>% filter(pval < 0.05) # Display significant models in a table datatable(sig_shannon_lms_df, options = list(pageLength = 40)) ### Community-Wide Production vs Shannon Entropy sig_shannon_lms_comm <- sig_shannon_lms_df %>% filter(test == "Community-Wide Production" & fraction == "Particle") %>% dplyr::select(Removed) %>% .$Removed ### Community-Wide production vs Shannon Entropy ggplot(dplyr::filter(all_divs, measure == "Exponential_Shannon"), aes(y = frac_bacprod, x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Shannon Entropy") + ylab("Bacterial Production by Fraction") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Exponential_Shannon" & fraction == "Particle" & Removed %in% sig_shannon_lms_comm)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(fraction~Removed, scales = "free") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)) ### Per-Capita Production vs Shannon Entropy sig_shannon_lms_percap <- sig_shannon_lms_df %>% filter(test == "Per-Capita Production" & fraction == "Particle") %>% dplyr::select(Removed) %>% .$Removed ### Per-capita production vs Shannon Entropy ggplot(dplyr::filter(all_divs, measure == "Exponential_Shannon"), aes(y = log10(fracprod_per_cell_noinf), x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Shannon Entropy") + ylab("log10(Per-Capita Production)") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Exponential_Shannon" & fraction == "Particle" & Removed %in% sig_shannon_lms_percap)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(fraction~Removed, scales = "free") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)) ``` # Inverse Simpson ```{r invsimps-plots1, fig.width=9, fig.height=4, dependson="all_divs"} ### PLOT ggplot(dplyr::filter(all_divs, measure == "Inverse_Simpson"), aes(y = mean, x = Removed, color = Removed, fill = Removed)) + geom_boxplot(alpha = 0.3, outlier.shape = NA) + geom_point(size = 3, position = position_jitter(w = 0.1)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(.~fraction) + ylab("Mean Inverse_Simpson") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1), axis.title.x = element_blank()) ggplot(dplyr::filter(all_divs, measure == "Inverse_Simpson"), aes(y = mean, x = fraction, color = Removed, fill = Removed)) + geom_boxplot(alpha = 0.3, outlier.shape = NA) + geom_point(size = 3, position = position_jitter(w = 0.1)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(.~Removed) + ylab("Mean Inverse_Simpson") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1), axis.title.x = element_blank()) ``` ```{r inverse-simpson-plots, fig.width=10, fig.height=5, dependson="all_divs"} # Linear Model output invsimps_lm_results <- lm_fraction_output(dataframe = dplyr::filter(all_divs, measure == "Inverse_Simpson")) sig_invsimps_lms_df <- invsimps_lm_results %>% bind_rows() %>% mutate(diversity_metric = "Inverse_Simpson") %>% filter(pval < 0.05) # Display significant models in a dataframe datatable(sig_invsimps_lms_df, options = list(pageLength = 40)) ### Community-Wide Production vs Inverse Simpson sig_invsimps_lms_comm <- sig_invsimps_lms_df %>% filter(test == "Community-Wide Production" & fraction == "Particle") %>% dplyr::select(Removed) %>% .$Removed ### Community Wide production vs Inverse Simpson ggplot(dplyr::filter(all_divs, measure == "Inverse_Simpson"), aes(y = frac_bacprod, x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Inverse Simpson") + ylab("Community-Wide Production") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Inverse_Simpson" & fraction == "Particle" & Removed %in% sig_invsimps_lms_comm)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(fraction~Removed, scales = "free") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)) ### Per-Capita Production vs Inverse Simpson sig_invsimps_lms_percap <- sig_invsimps_lms_df %>% filter(test == "Per-Capita Production" & fraction == "Particle") %>% dplyr::select(Removed) %>% .$Removed ### Per-capita production vs Inverse Simpson ggplot(dplyr::filter(all_divs, measure == "Inverse_Simpson"), aes(y = log10(fracprod_per_cell_noinf), x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Inverse Simpson") + ylab("log10(Per-Capita Production)") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Inverse_Simpson" & fraction == "Particle" & Removed %in% sig_invsimps_lms_percap)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(fraction~Removed, scales = "free") + theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)) ``` # Figure S9 ```{r Figure-S9, fig.height=14, fig.width=12} p1 <- ggplot(dplyr::filter(all_divs, measure == "Richness" & fraction == "Particle"), aes(y = frac_bacprod, x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Richness") + ylab("Community-Wide \nProduction") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Richness" & fraction == "Particle" & Removed %in% sig_rich_lms_comm)) + scale_color_manual(values = tons_colors) +scale_fill_manual(values = tons_colors) + facet_grid(~Removed, scales = "free") + theme(legend.position = "none", legend.title = element_blank(), axis.text.x = element_blank(),axis.title.x = element_blank()) p2 <- ggplot(dplyr::filter(all_divs, measure == "Richness" & fraction == "Particle"), aes(y = log10(fracprod_per_cell_noinf), x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Richness") + ylab("log10(Per-Capita \nProduction)") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Richness" & fraction == "Particle" & Removed %in% sig_rich_lms_percap)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(~Removed, scales = "free") + theme(legend.position = "none", legend.title = element_blank(), axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 10), axis.title.x = element_text(face = "bold")) p3 <- ggplot(dplyr::filter(all_divs, measure == "Inverse_Simpson" & fraction == "Particle"), aes(y = frac_bacprod, x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Inverse Simpson") + ylab("Community-Wide \nProduction") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Inverse_Simpson" & fraction == "Particle" & Removed %in% sig_invsimps_lms_comm)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(~Removed, scales = "free") + theme(legend.position = "none", legend.title = element_blank(), axis.text.x = element_blank(),axis.title.x = element_blank()) p4 <- ggplot(dplyr::filter(all_divs, measure == "Inverse_Simpson" & fraction == "Particle"), aes(y = log10(fracprod_per_cell_noinf), x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Inverse Simpson") + ylab("log10(Per-Capita \nProduction)") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Inverse_Simpson" & fraction == "Particle" & Removed %in% sig_invsimps_lms_percap)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(~Removed, scales = "free") + theme(legend.position = "none", legend.title = element_blank(), axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 10), axis.title.x = element_text(face = "bold")) p5 <- ggplot(dplyr::filter(all_divs, measure == "Exponential_Shannon" & fraction == "Particle"), aes(y = frac_bacprod, x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Exponential Shannon") + ylab("Community-Wide \nProduction") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Exponential_Shannon" & fraction == "Particle" & Removed %in% sig_shannon_lms_comm)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(~Removed, scales = "free") + theme(legend.position = "none", legend.title = element_blank(), axis.text.x = element_blank(),axis.title.x = element_blank()) p6 <- ggplot(dplyr::filter(all_divs, measure == "Exponential_Shannon" & fraction == "Particle"), aes(y = log10(fracprod_per_cell_noinf), x = mean, color = Removed, fill = Removed)) + geom_point(size = 3) + xlab("Exponential Shannon") + ylab("log10(Per-Capita \nProduction)") + geom_smooth(method = "lm", data = filter(all_divs, measure == "Exponential_Shannon" & fraction == "Particle" & Removed %in% sig_shannon_lms_percap)) + scale_color_manual(values = tons_colors) + scale_fill_manual(values = tons_colors) + facet_grid(~Removed, scales = "free") + theme(legend.position = "none", legend.title = element_blank(), axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 10), axis.title.x = element_text(face = "bold")) plot_grid(p1, p2, p5, p6, p3, p4, nrow = 6, ncol = 1, rel_heights = c(0.9, 1, 0.9, 1.2, 0.9, 1.2), labels = c("A", "B", "C", "D", "E", "F"), align = "v") ``` # Session Information ```{r session-info} devtools::session_info() # This will include session info with all R package version information ```