In cas-bioinf/metagenboot: Diagnostics for Metagenomic Data via Model-Based Bootstrapping

From the paper "Abiotic rather than biotic filtering shapes the arbuscular mycorrhizal fungal communities of European seminatural grasslands"

library(tidyverse)
library(here)
library(readxl)
library(vegan)
library(cowplot)
theme_set(theme_cowplot())

devtools::load_all()

The data has to downloaded manually at https://nph.onlinelibrary.wiley.com/doi/abs/10.1111/nph.14947 - TableS2 and saved to local_data/vanGeel2017.xlsx

main_data_file <- here("local_data","vanGeel2017.xlsx")
if(!file.exists(main_data_file)) {
  stop(paste0("The file `", main_data_file, "` not found, you need to download it manually."))
}

The important part - read and format the data so that rows are observations (biological samples) and columns are OTUs. Also gather basic information about the sampples (I call this "mapping" but maybe not a great name).

otus_raw <- read_excel(main_data_file, sheet = "Sample OTU matrix")
otu_tab <- otus_raw %>% select(OTU_16:OTU_5257) %>% as.matrix()

mapping <- otus_raw %>% select(species_id:plot) %>% mutate(Sample = paste0("S_",1:n()), type = if_else(species_id == "soil", "Soil", "Root"))

rownames(otu_tab) <- mapping$Sample

A custom plot function that wraps plot_mds and fixes some aesthetics that are reused - completely unnecessary, just convenient.

plot_mds_van_geel <- function(base_mds, mapping, aligned_samples = NULL, show_paths = "none", ...) {
  plot_mds(base_mds, mapping, aligned_samples = aligned_samples, color_aes = grassland_type, show_paths = show_paths, sample_point_alpha = 0.1)
}

set.seed(201905103)
base_mds <- otu_tab %>% metaMDS(trymax = 100, parallel = metagenboot_options('cores'))
plot_mds_van_geel(base_mds, mapping) + facet_wrap(~type)
plot_mds_van_geel(base_mds, mapping)

Running the Dirichlet-multinomial bootstrap

#Here, group_column needs to be filled, it is needed for the connectivity metric (not directly using right now, but it is of interest). group_column should be the main factor in the data.
sens <- mds_sensitivity_check(20, otu_tab, mapping, group_column = grassland_type, trace = 0)
plot(sens,color_aes = grassland_type, show_paths = "none", sample_point_alpha = 0.1)
sens$connectivity_stats
sens$consistency_stats

cas-bioinf/metagenboot documentation built on Feb. 25, 2021, 3:58 p.m.