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

From the paper "Characterization of the Skin Microbiota of the Cane Toad Rhinella cf. marina in Puerto Rico and Costa Rica"

The data can be downloaded at https://www.frontiersin.org/articles/10.3389/fmicb.2017.02624/full#supplementary-material

library(MASS)
library(tidyverse)
library(here)
library(vegan)
library(cowplot)
library(doParallel)
library(readxl)

options(mc.cores = parallel::detectCores())
rstan::rstan_options(auto_write = TRUE)
devtools::load_all()

registerDoParallel(parallel::detectCores())

otus_raw <- read_excel(here("private_data","abarca2018","table 3.xlsx"), sheet ="Supplementary-table7-otu50%",range = "A2:V5154")


otu_tab <- otus_raw %>% dplyr::select(-`#OTU ID`) %>% as.matrix() %>% t()
colnames(otu_tab) <- otus_raw$`#OTU ID`

mapping <- tibble(Country = "PR", Site = "Santa Ana", Sample = paste0("S",1:10)) %>% 
  rbind(tibble(Country = "CR", Site = "Sarapiqui", Sample = paste0("S", 1:7,"__1"))) %>%
  rbind(tibble(Country = "CR", Site = "Turrialba", Sample = paste0("T",1:4)))

plot_mds_abarca <- function(base_mds, mapping, aligned_samples = NULL) {
  plot_mds(base_mds, mapping, aligned_samples = aligned_samples, color_aes = Country, shape_aes = Site) +
    scale_color_manual(values = c(CR = "red", PR = "black")) +
    scale_shape_manual(values = c("Santa Ana" = 21, "Turrialba" = 3, "Sarapiqui" = 2)) +
    scale_y_continuous(trans = "reverse")
}

Can we get the original plot?

set.seed(20190510)
nmds_func <- function(x) {
  x %>% vegdist() %>% MASS::isoMDS()
}
base_mds <- otu_tab %>%   nmds_func()
plot_mds_abarca(base_mds, mapping)

base_meta_mds <- otu_tab %>% metaMDS()
plot_mds_abarca(base_meta_mds, mapping)

mds_samples_rarefy <- sample_posterior_rarefy(20, otu_tab) %>% apply_per_sample(nmds_func) %>% purrr::map(vegan::procrustes, X = base_mds)

plot_mds_abarca(base_mds, mapping, aligned_samples = mds_samples_rarefy)

mds_samples_dm <- sample_posterior_dm(20, otu_tab, prior = "ml",N_reads = "original") %>% apply_per_sample(nmds_func) %>% purrr::map(vegan::procrustes, X = base_mds)

plot_mds_abarca(base_mds, mapping, aligned_samples = mds_samples_dm)

mds_samples_dm_mds <- sample_posterior_dm(20, otu_tab, prior = "ml") %>% metaMDS_per_sample() %>% purrr::map(vegan::procrustes, X = base_meta_mds)


plot_mds_abarca(base_meta_mds, mapping, aligned_samples = mds_samples_dm_mds)

mds_samples_rarefy_mds <- sample_posterior_rarefy(20, otu_tab) %>% metaMDS_per_sample() %>% purrr::map(vegan::procrustes, X = base_meta_mds)


plot_mds_abarca(base_meta_mds,  mapping, aligned_samples = mds_samples_rarefy_mds)

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