ancombc_example.Rmd

---
  title: "Exploratory analysis for Comparative Metataxonomy"
author: "Roberto Siani"
date: "24.06.21"

---

  # SETUP

  Set up working environment: packages and default graphic themes.

```{r, cache = TRUE, echo = FALSE, include = FALSE}

# pacman to install and load libraries

if (!require("pacman")) install.packages(
  "pacman",
  verbose = F)

# BiocManager for Bioconductor libraries

if (!require("BiocManager")) install.packages(
  "BiocManager",
  verbose = F)

# devtools for GitHub repositories

if (!require("devtools")) install.packages(
  "devtools",
  verbose = F)

# GitHub libraries

pacman::p_load_gh(
  "jbisanz/qiime2R",
  "benjjneb/decontam",
  "mikemc/speedyseq",
  "adw96/breakaway",
  "adw96/DivNet")

# install/load the remainder of libraries

pacman::p_load(
  tidyverse,
  vegan,
  ape,
  microbiome,
  ggpubr,
  patchwork,
  hrbrthemes,
  phangorn,
  DECIPHER,
  MicrobiotaProcess,
  hues)

# load sensible themes and palettes

source("~/Desktop/library/scripts/helpeR.R")

theme_set(my_theme)

select = dplyr::select
transform = microbiome::transform

```

# INPUT

Congratulations, now that you have your phyloseq object prepared, you can finally start to explore the microbial wonderland in your samples

```{r, echo = FALSE, include = FALSE}

# load your uBiome dataset

load("~/Desktop/ancombc/uBiome_181021_1511.RData")

# again, take a good look at your results. Look at the sparsity and the percentage of singletons, select your variables of interest

summarize_phyloseq(uBiome)

## check your samples and groups

table(meta(uBiome)$Description, meta(uBiome)$Bacteria)

## set up palette

uBiome = 
  uBiome %>% 
  subset_samples(Bacteria == "WT1" | Bacteria == "LjR176")

palette_list =
  list(
    Description = c("#cc7722", "#2277cc"))

map(palette_list,
    ~swatches::show_palette(.x))

```

### ANCOM-BC

```{r}

p_load(ANCOMBC)

res_ancombc =
  ancombc(uBiome,
          formula = "Description",
          p_adj_method = "fdr",
          conserve = T)


## Description

res_df =
  res_ancombc$res$beta %>%
  as.data.frame() %>%
  rownames_to_column("ASV") %>%
  left_join(
    res_ancombc$res$q_val %>%
      as.data.frame() %>%
      rownames_to_column("ASV"),
    by = "ASV",
    suffix = c("_beta", "_q")) %>%
  left_join(
    res_ancombc$res$W %>%
      as.data.frame() %>%
      rownames_to_column("ASV")) %>%
  rename(DescriptionInoculated_W = DescriptionInoculated) %>%
  mutate(DescriptionInoculated_DA = as.factor(case_when(
    DescriptionInoculated_q > 0.01 & DescriptionInoculated_q <= 0.05 ~ "*",
    DescriptionInoculated_q > 0.001 & DescriptionInoculated_q <= 0.01 ~ "**",
    DescriptionInoculated_q > 0.0001 & DescriptionInoculated_q <= 0.001 ~ "***",
    DescriptionInoculated_q <= 0.0001 ~ "****",
    TRUE ~ "ns"))) %>%
  pivot_longer(-ASV,
               names_to = c("var", ".value"),
               names_sep = "_") %>%
  mutate(var = factor(var, levels = c("DescriptionInoculated"),
                      labels = c("<- Control|Treatment ->")))

res_df$DA = relevel(res_df$DA, ref = "ns")

plot_7 =
  ggplot(res_df) +
  geom_point(
    aes(x = beta,
        y = -log10(q + 1e-5),
        color = DA)) +
  scale_color_manual(values = michelangelo$grad3) +
  facet_wrap(~var) +
  theme(legend.position = "right") +
  labs(y = expression(-log[10]~q-value))


relAbu_df =
  res_df %>%
  filter(q <= 0.05) %>%
  left_join(tax_table(uBiome) %>%
              as.data.frame() %>%
              rownames_to_column("ASV")) %>%
  mutate(Genus = vctrs::vec_as_names(Genus, repair = "unique") %>%
           sub("...", "_", fixed = T, .))

plot_8 =
  relAbu_df %>%
  ggplot() +
  geom_bar(aes(y = fct_reorder(Genus, W),
               x = W,
               fill = Phylum),
           stat = "identity",
           position = "dodge") +
  scale_fill_iwanthue(hmin = 0,
                      hmax = 360,
                      cmin = 0,
                      cmax = 30,
                      lmin = 50,
                      lmax = 75) +
  theme(axis.text.y = element_text(size = 12),
        axis.title.y = element_blank(),
        legend.position = "right",
        strip.text = element_blank()) +
  guides(fill = guide_legend(ncol = 1))


(plot_7 | plot_8) + plot_layout(widths = c(2, 1))

```