ANCOMBC2.Rmd

---
title: "ANCOM-BC2 Tutorial"
author: 
  - Huang Lin$^1$
  - $^1$NIEHS, Research Triangle Park, NC 27709, USA
date: '`r format(Sys.Date(), "%B %d, %Y")`'
output: rmarkdown::html_vignette
bibliography: bibliography.bib
vignette: >
  %\VignetteIndexEntry{ANCOM-BC2 Tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, message = FALSE, warning = FALSE, comment = NA}
knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = NA, 
                      fig.width = 6.25, fig.height = 5)
library(ANCOMBC)
library(tidyverse)
library(DT)
options(DT.options = list(
  initComplete = JS("function(settings, json) {",
  "$(this.api().table().header()).css({'background-color': 
  '#000', 'color': '#fff'});","}")))
```

# 1. Introduction

Analysis of Compositions of Microbiomes with Bias Correction 2 (ANCOM-BC2) is 
a methodology for performing differential abundance (DA) analysis of microbiome 
count data. This version extends and refines the previously published 
Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC) 
methodology [@lin2020analysis] in several ways as follows:

1. **Bias correction**: ANCOM-BC2 estimates and corrects both the 
sample-specific (sampling fraction) as well as 
taxon-specific (sequencing efficiency) biases. 

2. **Regularization of variance**: Inspired by Significance Analysis of 
Microarrays (SAM) [@tusher2001significance] methodology, a small positive 
constant is added to the denominator of ANCOM-BC2 test statistic corresponding 
to each taxon to avoid the significance due to extremely small standard errors, 
especially for rare taxa. By default, we used the 5-th percentile of the 
distribution of standard errors for each fixed effect as the regularization 
factor. 

3. **Sensitivity analysis for the pseudo-count addition**: Like other 
differential abundance analysis methods, ANCOM-BC2 applies a log transformation 
to the observed counts. However, the presence of zero counts poses a challenge, 
and researchers often consider adding a pseudo-count before the log 
transformation. However, it has been shown that the choice of pseudo-count can 
impact the results and lead to an inflated false positive 
rate [@costea2014fair; @paulson2014reply]. To address this issue, 
we conduct a sensitivity analysis to assess the impact of different 
pseudo-counts on zero counts for each taxon. This analysis involves adding a 
series of pseudo-counts (ranging from 0.01 to 0.5 in increments of 0.01) to 
the zero counts of each taxon. Linear regression models are then performed on 
the bias-corrected log abundance table using the different pseudo-counts. The 
sensitivity score for each taxon is calculated as the proportion of times 
that the p-value exceeds the specified significance level (alpha). If all 
p-values consistently show significance or nonsignificance across different 
pseudo-counts and are consistent with the results obtained without adding 
pseudo-counts to zero counts (using the default settings), then the taxon is 
considered not sensitive to the pseudo-count addition. 

4. **Multi-group comparisons and repeated measurements**: The ANCOM-BC2 
methodology extends ANCOM-BC for multiple groups and repeated measurements 
as follows:

    + Multiple pairwise comparisons: When performning multiple pairwise 
    comparisons, the mixed directional false discover rate (mdFDR) should be 
    taken into account. The mdFDR is the combination of false 
    discovery rate due to multiple testing, multiple pairwise comparisons, 
    and directional tests within each pairwise comparison. For example, suppose 
    we have five taxa and three experimental groups: g1, g2, and g3. 
    Thus, we are performing five tests corresponding to five taxa. 
    For each taxon, we are also conducting three pairwise comparisons
    (g1 vs. g2, g2 vs. g3, and g1 vs. g3). Within each pairwise comparison,
    we wish to determine if the abundance has increased or decreased or did not
    change (direction of the effect size). Errors could occur in each step.
    The overall false discovery rate is controlled by the mdFDR methodology we
    adopted from [@guo2010controlling; @grandhi2016multiple].
    
    + Multiple pairwise comparisons against a pre-specified group 
    (e.g., Dunnett's type of test): We use the same set-up as in the multiple 
    pairwise comparisons but use the Dunnett-type modification described 
    in [@grandhi2016multiple] to control the mdFDR which is more powerful.
    
    + Pattern analysis for ordered groups: In some instances, researchers are 
    interested in discovering abundance patterns of each taxon over the ordered 
    groups, for example, groups based on the health condition of subjects 
    (e.g., lean, overweight, obese). In such cases, in addition to pairwise 
    comparison, one may be interested in identifying taxa whose abundances are 
    increasing or decreasing or have other patterns over the groups.
    We adopted methodologies from [@jelsema2016clme] to 
    perform pattern analysis under the ANCOM-BC2 framework.

**A clarification regarding Structural zeros**: A taxon is considered to have 
structural zeros in some (>=1) groups if it is completely (or nearly completely) 
missing in these groups. For instance, suppose there are three groups: 
g1, g2, and g3. If the counts of taxon A in g1 are 0, but they are nonzero 
in g2 and g3, then taxon A will be considered to contain structural zeros in g1. 
In this example, taxon A is declared to be differentially abundant between 
g1 and g2, g1 and g3, and consequently, it is globally differentially abundant 
with respect to this group variable. Such taxa are not further analyzed using 
ANCOM-BC2, but the results are summarized in the overall summary. 

# 2. Installation

Download the package. 

```{r getPackage, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("ANCOMBC")
```

Load the package. 

```{r load, eval=FALSE}
library(ANCOMBC)
```

# 3. Run ANCOM-BC2 on a simulated dataset {.tabset}

## 3.1 Generate simulated data

1. The simulated data contain: one continuous covariate: cont_cov, and one
categorical covariate: cat_cov. The categorical covariate has three levels: "1", 
"2", and "3". Let's focus on the discussions on the group variable: cat_cov.

2. The true abundances are generated using the Poisson lognormal (PLN) 
model based on the mechanism described in the LDM paper [@hu2020testing]. The
PLN model relates the abundance vector $Y_i$ with a Gaussian latent vector
$Z_i$ for taxon $i$ as follows:
$$ \text{latent layer: } Z_i \sim N(\mu_i, \Sigma_i) \\ 
\text{observation layer: } Y_i|Z_i \sim POI(N \exp(Z_i)) $$
where $N$ is the scaling factor. Because of the presence of a latent layer, the 
PLN model displays a larger variance than the Poisson model (over-dispersion).
Also, the covariance (correlation) between abundances has the same sign as 
the covariance (correlation) between the corresponding latent variables. 
This property gives enormous flexibility in modeling the variance-covariance
structure of microbial abundances since it is easy to specify different
variance-covariance matrices in the multivariate Gaussian distribution.

3. Instead of specifying the variance-covariance matrix, we choose to estimate
the variance-covariance matrix from a real dataset: the Quantitative Microbiome 
Project (QMP) data [@vandeputte2017quantitative]. This dataset contains 
quantitative microbiome count data of 106 samples and 91 OTUs.

```{r}
data(QMP, package = "ANCOMBC")
set.seed(123)
n = 150
d = ncol(QMP)
diff_prop = 0.1
lfc_cont = 1
lfc_cat2_vs_1 = -2
lfc_cat3_vs_1 = 1

# Generate the true abundances
abn_data = sim_plnm(abn_table = QMP, taxa_are_rows = FALSE, prv_cut = 0.05, 
                    n = n, lib_mean = 1e8, disp = 0.5)
log_abn_data = log(abn_data + 1e-5)
rownames(log_abn_data) = paste0("T", seq_len(d))
colnames(log_abn_data) = paste0("S", seq_len(n))

# Generate the sample and feature meta data
# Sampling fractions are set to differ by batches
smd = data.frame(samp_frac = log(c(runif(n/3, min = 1e-4, max = 1e-3),
                                   runif(n/3, min = 1e-3, max = 1e-2),
                                   runif(n/3, min = 1e-2, max = 1e-1))),
                 cont_cov = rnorm(n),
                 cat_cov = as.factor(rep(seq_len(3), each = n/3)))
rownames(smd) = paste0("S", seq_len(n))
                      
fmd = data.frame(taxon = paste0("T", seq_len(d)),
                 seq_eff = log(runif(d, min = 0.1, max = 1)),
                 lfc_cont = sample(c(0, lfc_cont), 
                                   size = d,
                                   replace = TRUE,
                                   prob = c(1 - diff_prop, diff_prop)),
                 lfc_cat2_vs_1 = sample(c(0, lfc_cat2_vs_1), 
                                        size = d,
                                        replace = TRUE,
                                        prob = c(1 - diff_prop, diff_prop)),
                 lfc_cat3_vs_1 = sample(c(0, lfc_cat3_vs_1), 
                                        size = d,
                                        replace = TRUE,
                                        prob = c(1 - diff_prop, diff_prop))) %>%
    mutate(lfc_cat3_vs_2 = lfc_cat3_vs_1 - lfc_cat2_vs_1)

# Add effect sizes of covariates to the true abundances
smd_dmy = model.matrix(~ 0 + cont_cov + cat_cov, data = smd)
log_abn_data = log_abn_data + outer(fmd$lfc_cont, smd_dmy[, "cont_cov"] )
log_abn_data = log_abn_data + outer(fmd$lfc_cat2_vs_1, smd_dmy[, "cat_cov2"])
log_abn_data = log_abn_data + outer(fmd$lfc_cat3_vs_1, smd_dmy[, "cat_cov3"])

# Add sample- and taxon-specific biases
log_otu_data = t(t(log_abn_data) + smd$samp_frac)
log_otu_data = log_otu_data + fmd$seq_eff
otu_data = round(exp(log_otu_data))

# Create the tse object
assays = S4Vectors::SimpleList(counts = otu_data)
smd = S4Vectors::DataFrame(smd)
tse = TreeSummarizedExperiment::TreeSummarizedExperiment(assays = assays, colData = smd)
```

## 3.2 Run ancombc2 function

```{r}
set.seed(123)
output = ancombc2(data = tse, assay_name = "counts", tax_level = NULL,
                  fix_formula = "cont_cov + cat_cov", rand_formula = NULL,
                  p_adj_method = "holm", 
                  prv_cut = 0.10, lib_cut = 1000, s0_perc = 0.05,
                  group = "cat_cov", struc_zero = FALSE, neg_lb = FALSE,
                  alpha = 0.05, n_cl = 2, verbose = TRUE,
                  global = FALSE, pairwise = TRUE, 
                  dunnet = FALSE, trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, 
                                      verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = NULL, 
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

res_prim = output$res
res_pair = output$res_pair
```

## 3.3 Power and FDR

```{r}
res_merge1 = res_pair %>%
  dplyr::transmute(taxon, 
                   lfc_est1 = lfc_cat_cov2 * diff_cat_cov2,
                   lfc_est2 = lfc_cat_cov3 * diff_cat_cov3,
                   lfc_est3 = lfc_cat_cov3_cat_cov2 * diff_cat_cov3_cat_cov2) %>%
  dplyr::left_join(fmd %>%
                     dplyr::transmute(taxon, 
                                      lfc_true1 = lfc_cat2_vs_1,
                                      lfc_true2 = lfc_cat3_vs_1,
                                      lfc_true3 = lfc_cat3_vs_2),
                   by = "taxon") %>%
  dplyr::transmute(taxon, 
                   lfc_est1 = case_when(lfc_est1 > 0 ~ 1,
                                        lfc_est1 < 0 ~ -1,
                                        TRUE ~ 0),
                   lfc_est2 = case_when(lfc_est2 > 0 ~ 1,
                                        lfc_est2 < 0 ~ -1,
                                        TRUE ~ 0),
                   lfc_est3 = case_when(lfc_est3 > 0 ~ 1,
                                        lfc_est3 < 0 ~ -1,
                                        TRUE ~ 0),
                   lfc_true1 = case_when(lfc_true1 > 0 ~ 1,
                                         lfc_true1 < 0 ~ -1,
                                         TRUE ~ 0),
                   lfc_true2 = case_when(lfc_true2 > 0 ~ 1,
                                         lfc_true2 < 0 ~ -1,
                                         TRUE ~ 0),
                   lfc_true3 = case_when(lfc_true3 > 0 ~ 1,
                                         lfc_true3 < 0 ~ -1,
                                         TRUE ~ 0))
lfc_est1 = res_merge1$lfc_est1
lfc_true1 = res_merge1$lfc_true1
lfc_est2 = res_merge1$lfc_est2
lfc_true2 = res_merge1$lfc_true2
lfc_est3 = res_merge1$lfc_est3
lfc_true3 = res_merge1$lfc_true3

tp1 = sum(lfc_true1 == 1 & lfc_est1 == 1) +
  sum(lfc_true1 == -1 & lfc_est1 == -1)
fp1 = sum(lfc_true1 == 0 & lfc_est1 != 0) +
  sum(lfc_true1 == 1 & lfc_est1 == -1) +
  sum(lfc_true1 == -1 & lfc_est1 == 1)
fn1 = sum(lfc_true1 != 0 & lfc_est1 == 0)

tp2 = sum(lfc_true2 == 1 & lfc_est2 == 1) +
  sum(lfc_true2 == -1 & lfc_est2 == -1)
fp2 = sum(lfc_true2 == 0 & lfc_est2 != 0) +
  sum(lfc_true2 == 1 & lfc_est2 == -1) +
  sum(lfc_true2 == -1 & lfc_est2 == 1)
fn2 = sum(lfc_true2 != 0 & lfc_est2 == 0)

tp3 = sum(lfc_true3 == 1 & lfc_est3 == 1) +
  sum(lfc_true3 == -1 & lfc_est3 == -1)
fp3 = sum(lfc_true3 == 0 & lfc_est3 != 0) +
  sum(lfc_true3 == 1 & lfc_est3 == -1) +
  sum(lfc_true3 == -1 & lfc_est3 == 1)
fn3 = sum(lfc_true3 != 0 & lfc_est3 == 0)

tp = tp1 + tp2 + tp3
fp = fp1 + fp2 + fp3
fn = fn1 + fn2 + fn3

power1 = tp/(tp + fn)
fdr1 = fp/(tp + fp)

res_merge2 = res_pair %>%
  dplyr::transmute(taxon, 
                   lfc_est1 = lfc_cat_cov2 * diff_cat_cov2 * passed_ss_cat_cov2,
                   lfc_est2 = lfc_cat_cov3 * diff_cat_cov3 * passed_ss_cat_cov3,
                   lfc_est3 = lfc_cat_cov3_cat_cov2 * diff_cat_cov3_cat_cov2* passed_ss_cat_cov3_cat_cov2) %>%
  dplyr::left_join(fmd %>%
                     dplyr::transmute(taxon, 
                                      lfc_true1 = lfc_cat2_vs_1,
                                      lfc_true2 = lfc_cat3_vs_1,
                                      lfc_true3 = lfc_cat3_vs_2),
                   by = "taxon") %>%
  dplyr::transmute(taxon, 
                   lfc_est1 = case_when(lfc_est1 > 0 ~ 1,
                                        lfc_est1 < 0 ~ -1,
                                        TRUE ~ 0),
                   lfc_est2 = case_when(lfc_est2 > 0 ~ 1,
                                        lfc_est2 < 0 ~ -1,
                                        TRUE ~ 0),
                   lfc_est3 = case_when(lfc_est3 > 0 ~ 1,
                                        lfc_est3 < 0 ~ -1,
                                        TRUE ~ 0),
                   lfc_true1 = case_when(lfc_true1 > 0 ~ 1,
                                         lfc_true1 < 0 ~ -1,
                                         TRUE ~ 0),
                   lfc_true2 = case_when(lfc_true2 > 0 ~ 1,
                                         lfc_true2 < 0 ~ -1,
                                         TRUE ~ 0),
                   lfc_true3 = case_when(lfc_true3 > 0 ~ 1,
                                         lfc_true3 < 0 ~ -1,
                                         TRUE ~ 0))
lfc_est1 = res_merge2$lfc_est1
lfc_true1 = res_merge2$lfc_true1
lfc_est2 = res_merge2$lfc_est2
lfc_true2 = res_merge2$lfc_true2
lfc_est3 = res_merge2$lfc_est3
lfc_true3 = res_merge2$lfc_true3

tp1 = sum(lfc_true1 == 1 & lfc_est1 == 1) +
  sum(lfc_true1 == -1 & lfc_est1 == -1)
fp1 = sum(lfc_true1 == 0 & lfc_est1 != 0) +
  sum(lfc_true1 == 1 & lfc_est1 == -1) +
  sum(lfc_true1 == -1 & lfc_est1 == 1)
fn1 = sum(lfc_true1 != 0 & lfc_est1 == 0)

tp2 = sum(lfc_true2 == 1 & lfc_est2 == 1) +
  sum(lfc_true2 == -1 & lfc_est2 == -1)
fp2 = sum(lfc_true2 == 0 & lfc_est2 != 0) +
  sum(lfc_true2 == 1 & lfc_est2 == -1) +
  sum(lfc_true2 == -1 & lfc_est2 == 1)
fn2 = sum(lfc_true2 != 0 & lfc_est2 == 0)

tp3 = sum(lfc_true3 == 1 & lfc_est3 == 1) +
  sum(lfc_true3 == -1 & lfc_est3 == -1)
fp3 = sum(lfc_true3 == 0 & lfc_est3 != 0) +
  sum(lfc_true3 == 1 & lfc_est3 == -1) +
  sum(lfc_true3 == -1 & lfc_est3 == 1)
fn3 = sum(lfc_true3 != 0 & lfc_est3 == 0)

tp = tp1 + tp2 + tp3
fp = fp1 + fp2 + fp3
fn = fn1 + fn2 + fn3

power2 = tp/(tp + fn)
fdr2 = fp/(tp + fp)

tab_summ = data.frame(Comparison = c("Without sensitivity score filter", 
                                     "With sensitivity score filter"),
                      Power = round(c(power1, power2), 2),
                      FDR = round(c(fdr1, fdr2), 2))
tab_summ %>%
    datatable(caption = "Power/FDR Comparison")
```

# 4. Run ANCOM-BC2 on a real cross-sectional dataset
    
## 4.1 Import example data

The HITChip Atlas dataset contains genus-level microbiota profiling with 
HITChip for 1006 western adults with no reported health complications, 
reported in [@lahti2014tipping]. The dataset is available via the 
microbiome R package [@lahti2017tools] in phyloseq [@mcmurdie2013phyloseq] 
format. In this tutorial, we consider the following covariates:

* Continuous covariates: "age"

* Categorical covariates: "region", "bmi"

* The group variable of interest: "bmi"

    + Three groups: "lean", "overweight", "obese"
    
    + The reference group: "obese"

```{r}
data(atlas1006, package = "microbiome")
tse = mia::makeTreeSummarizedExperimentFromPhyloseq(atlas1006)

# subset to baseline
tse = tse[, tse$time == 0]

# Re-code the bmi group
tse$bmi = recode(tse$bmi_group,
                 obese = "obese",
                 severeobese = "obese",
                 morbidobese = "obese")
# Subset to lean, overweight, and obese subjects
tse = tse[, tse$bmi %in% c("lean", "overweight", "obese")]

# Note that by default, levels of a categorical variable in R are sorted 
# alphabetically. In this case, the reference level for `bmi` will be 
# `lean`. To manually change the reference level, for instance, setting `obese`
# as the reference level, use:
tse$bmi = factor(tse$bmi, levels = c("obese", "overweight", "lean"))
# You can verify the change by checking:
# levels(sample_data(tse)$bmi)

# Create the region variable
tse$region = recode(as.character(tse$nationality),
                    Scandinavia = "NE", UKIE = "NE", SouthEurope = "SE", 
                    CentralEurope = "CE", EasternEurope = "EE",
                    .missing = "unknown")

# Discard "EE" as it contains only 1 subject
# Discard subjects with missing values of region
tse = tse[, ! tse$region %in% c("EE", "unknown")]

print(tse)
```
    
## 4.2 Run ancombc2 function

```{r}
set.seed(123)
# It should be noted that we have set the number of bootstrap samples (B) equal 
# to 10 in the 'trend_control' function for computational expediency. 
# However, it is recommended that users utilize the default value of B, 
# which is 100, or larger values for optimal performance.
output = ancombc2(data = tse, assay_name = "counts", tax_level = "Family",
                  fix_formula = "age + region + bmi", rand_formula = NULL,
                  p_adj_method = "holm", 
                  prv_cut = 0.10, lib_cut = 1000, s0_perc = 0.05,
                  group = "bmi", struc_zero = TRUE, neg_lb = TRUE,
                  alpha = 0.05, n_cl = 2, verbose = TRUE,
                  global = TRUE, pairwise = TRUE, dunnet = TRUE, trend = TRUE,
                  iter_control = list(tol = 1e-2, max_iter = 20, 
                                      verbose = TRUE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
                  trend_control = list(contrast = list(matrix(c(1, 0, -1, 1),
                                                              nrow = 2, 
                                                              byrow = TRUE),
                                                       matrix(c(-1, 0, 1, -1),
                                                              nrow = 2, 
                                                              byrow = TRUE)),
                                       node = list(2, 2),
                                       solver = "ECOS",
                                       B = 10))
```

## 4.3 Structural zeros (taxon presence/absence)

```{r}
tab_zero = output$zero_ind
tab_zero %>%
    datatable(caption = "The detection of structural zeros")
```

## 4.4 ANCOM-BC2 primary analysis {.tabset}

Result from the ANCOM-BC2 methodology to determine taxa that are 
differentially abundant according to the covariate of interest. Results contain: 
1) log fold changes, 2) standard errors, 3) test statistics, 4) p-values, 
5) adjusted p-values, 6) indicators of whether the taxon is differentially 
abundant (TRUE) or not (FALSE).

```{r}
res_prim = output$res
```

### Results for age 

```{r}
df_age = res_prim %>%
    dplyr::select(taxon, ends_with("age")) 
df_fig_age = df_age %>%
    dplyr::filter(diff_age == 1, passed_ss_age == 1) %>% 
    dplyr::arrange(desc(lfc_age)) %>%
    dplyr::mutate(direct = ifelse(lfc_age > 0, "Positive LFC", "Negative LFC"))
df_fig_age$taxon = factor(df_fig_age$taxon, levels = df_fig_age$taxon)
df_fig_age$direct = factor(df_fig_age$direct, 
                           levels = c("Positive LFC", "Negative LFC"))
  
fig_age = df_fig_age %>%
    ggplot(aes(x = taxon, y = lfc_age, fill = direct)) + 
    geom_bar(stat = "identity", width = 0.7, color = "black", 
             position = position_dodge(width = 0.4)) +
    geom_errorbar(aes(ymin = lfc_age - se_age, ymax = lfc_age + se_age), 
                  width = 0.2, position = position_dodge(0.05), color = "black") + 
    labs(x = NULL, y = "Log fold change", 
         title = "Log fold changes as one unit increase of age") + 
    scale_fill_discrete(name = NULL) +
    scale_color_discrete(name = NULL) +
    theme_bw() + 
    theme(plot.title = element_text(hjust = 0.5),
          panel.grid.minor.y = element_blank(),
          axis.text.x = element_text(angle = 60, hjust = 1))
fig_age
```

### Results for bmi

```{r}
df_bmi = res_prim %>%
    dplyr::select(taxon, contains("bmi")) 
df_fig_bmi = df_bmi %>%
    dplyr::filter((diff_bmilean == 1 & passed_ss_bmilean == 1) | 
                    (diff_bmioverweight == 1 & passed_ss_bmioverweight == 1)) %>%
    dplyr::mutate(lfc_overweight = ifelse(diff_bmioverweight == 1, 
                                          lfc_bmioverweight, 0),
                  lfc_lean = ifelse(diff_bmilean == 1, 
                                    lfc_bmilean, 0)) %>%
    dplyr::transmute(taxon, 
                     `Overweight vs. Obese` = round(lfc_overweight, 2),
                     `Lean vs. Obese` = round(lfc_lean, 2)) %>%
    tidyr::pivot_longer(cols = `Overweight vs. Obese`:`Lean vs. Obese`, 
                        names_to = "group", values_to = "value") %>%
    dplyr::arrange(taxon)
  
lo = floor(min(df_fig_bmi$value))
up = ceiling(max(df_fig_bmi$value))
mid = (lo + up)/2
fig_bmi = df_fig_bmi %>%
  ggplot(aes(x = group, y = taxon, fill = value)) + 
  geom_tile(color = "black") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       na.value = "white", midpoint = mid, limit = c(lo, up),
                       name = NULL) +
  geom_text(aes(group, taxon, label = value), color = "black", size = 4) +
  labs(x = NULL, y = NULL, title = "Log fold changes as compared to obese subjects") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
fig_bmi
```

## 4.5 ANCOM-BC2 global test

ANCOM-BC2 global test aims to determine taxa that are differentially abundant 
between at least two groups across three or more experimental groups.

In this example, we want to identify taxa that are differentially abundant 
between at least two groups across "lean", "overweight", and "obese".
The result contains: 1) test statistics, 2) p-values, 3) adjusted p-values, 
4) indicators of whether the taxon is differentially abundant (TRUE) or not 
(FALSE).

```{r}
res_global = output$res_global
df_bmi = res_prim %>%
    dplyr::select(taxon, contains("bmi")) 
df_fig_global = df_bmi %>%
    dplyr::left_join(res_global %>%
                       dplyr::transmute(taxon, 
                                        diff_bmi = diff_abn, 
                                        passed_ss = passed_ss)) %>%
    dplyr::filter(diff_bmi == 1, passed_ss == 1) %>%
    dplyr::mutate(lfc_lean = lfc_bmilean,
                  lfc_overweight = lfc_bmioverweight) %>%
    dplyr::transmute(taxon, 
                     `Lean vs. Obese` = round(lfc_lean, 2), 
                     `Overweight vs. Obese` = round(lfc_overweight, 2)) %>%
    tidyr::pivot_longer(cols = `Lean vs. Obese`:`Overweight vs. Obese`, 
                        names_to = "group", values_to = "value") %>%
    dplyr::arrange(taxon)
  
lo = floor(min(df_fig_global$value))
up = ceiling(max(df_fig_global$value))
mid = (lo + up)/2
fig_global = df_fig_global %>%
  ggplot(aes(x = group, y = taxon, fill = value)) + 
  geom_tile(color = "black") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       na.value = "white", midpoint = mid, limit = c(lo, up),
                       name = NULL) +
  geom_text(aes(group, taxon, label = value), color = "black", size = 4) +
  labs(x = NULL, y = NULL, title = "Log fold changes for globally significant taxa") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
fig_global
```

## 4.6 ANCOM-BC2 multiple pairwise comparisons

ANCOM-BC2 multiple pairwise comparisons aim to determine taxa that are 
differentially abundant between any pair of two groups across three or more 
experimental groups, while controlling the mdFDR.

In this example, we want to identify taxa that are differentially abundant 
between any pair of two groups across "lean", "overweight", and "obese".
The result contains: 1) log fold changes, 2) standard errors, 3) test 
statistics, 4) p-values, 5) adjusted p-values, 6) indicators of whether the 
taxon is differentially abundant (TRUE) or not (FALSE).

```{r}
res_pair = output$res_pair

df_fig_pair = res_pair %>%
    dplyr::filter((diff_bmilean == 1 & passed_ss_bmilean == 1) | 
                    (diff_bmioverweight == 1 & passed_ss_bmioverweight == 1) |
                    (diff_bmilean_bmioverweight == 1 & passed_ss_bmilean_bmioverweight == 1)) %>%
    dplyr::mutate(lfc_lean = ifelse(diff_bmilean == 1, 
                                    lfc_bmilean, 0),
                  lfc_overweight = ifelse(diff_bmioverweight == 1, 
                                          lfc_bmioverweight, 0),
                  lfc_lean_overweight = ifelse(diff_bmilean_bmioverweight == 1, 
                                               lfc_bmilean_bmioverweight, 0)) %>%
    dplyr::transmute(taxon, 
                     `Lean vs. Obese` = round(lfc_lean, 2),
                     `Overweight vs. Obese` = round(lfc_overweight, 2),
                     `Lean vs. Overweight` = round(lfc_lean_overweight, 2)
                     ) %>%
    tidyr::pivot_longer(cols = `Lean vs. Obese`:`Lean vs. Overweight`, 
                        names_to = "group", values_to = "value") %>%
    dplyr::arrange(taxon)
df_fig_pair$group = factor(df_fig_pair$group, 
                           levels = c("Lean vs. Obese",
                                      "Overweight vs. Obese",
                                      "Lean vs. Overweight"))
  
lo = floor(min(df_fig_pair$value))
up = ceiling(max(df_fig_pair$value))
mid = (lo + up)/2
fig_pair = df_fig_pair %>%
  ggplot(aes(x = group, y = taxon, fill = value)) + 
  geom_tile(color = "black") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       na.value = "white", midpoint = mid, limit = c(lo, up),
                       name = NULL) +
  geom_text(aes(group, taxon, label = value), color = "black", size = 4) +
  labs(x = NULL, y = NULL, title = "Log fold change of pairwise comparisons") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
fig_pair
```

## 4.7 ANCOM-BC2 multiple pairwise comparisons against a pre-specified group (Dunnett's type of test)

The Dunnett’s test [@dunnett1955multiple; @dunnett1991step; @dunnett1992step] 
is designed for making comparisons of several experimental groups with the 
control or the reference group. ANCOM-BC2 Dunnett's type of test adopts the 
framework of Dunnett's test while controlling the mdFDR. Of note is that the
ANCOM-BC2 primary results do not control the mdFDR for the comparison of 
multiple groups.

In this example, we want to identify taxa that are differentially abundant 
between "lean", "overweight", and the reference group "obese".
The result contains: 1) log fold changes, 2) standard errors, 3) test 
statistics, 4) p-values, 5) adjusted p-values, 6) indicators of whether the 
taxon is differentially abundant (TRUE) or not (FALSE).

```{r}
res_dunn = output$res_dunn

df_fig_dunn = res_dunn %>%
    dplyr::filter((diff_bmilean == 1 & passed_ss_bmilean == 1) | 
                    (diff_bmioverweight == 1 & passed_ss_bmioverweight == 1)) %>%
    dplyr::mutate(lfc_lean = ifelse(diff_bmilean == 1, lfc_bmilean, 0),
                  lfc_overweight = ifelse(diff_bmioverweight == 1, 
                                          lfc_bmioverweight, 0)) %>%
    dplyr::transmute(taxon, 
                     `Lean vs. Obese` = round(lfc_lean, 2), 
                     `Overweight vs. Obese` = round(lfc_overweight, 2)) %>%
    tidyr::pivot_longer(cols = `Lean vs. Obese`:`Overweight vs. Obese`, 
                        names_to = "group", values_to = "value") %>%
    dplyr::arrange(taxon)
  
lo = floor(min(df_fig_dunn$value))
up = ceiling(max(df_fig_dunn$value))
mid = (lo + up)/2
fig_dunn = df_fig_dunn %>%
  ggplot(aes(x = group, y = taxon, fill = value)) + 
  geom_tile(color = "black") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       na.value = "white", midpoint = mid, limit = c(lo, up),
                       name = NULL) +
  geom_text(aes(group, taxon, label = value), color = "black", size = 4) +
  labs(x = NULL, y = NULL, title = "Log fold changes as compared to obese subjects") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
fig_dunn
```

## 4.8 ANCOM-BC2 pattern analysis {.tabset}

Sometimes the experimental groups are intrinsically ordered (e.g., in a 
dose-response study), and we would like to see if the microbial abundances
show some patterns accordingly. Examples of patterns include: monotonically
increasing, monotonically decreasing, and umbrella shape.

ANCOM-BC2 is able to identify potential patterns by testing the contrast:
$$Ax \ge 0$$
where $A$ is the contrast matrix and $x$ is the vector of parameters.

For instance, in this example, we want to identify taxa that are monotonically
increasing across "obese", "overweight", and "lean". Note that "obese" is the 
reference group, and the parameters we can estimate are the differences as 
compared to the reference group, i.e., 
$$x = (\text{overweight - obese}, \text{lean - obese})^T$$
To test the monotonically increasing trend:
$$H_0: \text{obese} = \text{overweight} = \text{lean} \\
H_1: \text{obese} \le \text{overweight} \le \text{lean} \quad \text{with at least one strict inequality}$$
We shall specify the contrast matrix $A$ as
$$A = \begin{bmatrix} 1 & 0 \\ -1 & 1 \end{bmatrix}$$
Similarly, to test for the monotonically decreasing trend:
$$H_0: \text{obese} = \text{overweight} = \text{lean} \\
H_1: \text{obese} \ge \text{overweight} \ge \text{lean} \quad \text{with at least one strict inequality}$$
We shall specify the contrast matrix $A$ as
$$A = \begin{bmatrix} -1 & 0 \\ 1 & -1 \end{bmatrix}$$
To test for monotonic trend (increasing or decreasing), one should specify the 
`node` parameter in `trend_control` as the last position of `x`.
In this example, the vector of parameters $x$ is of length 2, thus, the last 
position is 2. For testing umbrella shape, for instance, in this example:
$$H_0: \text{obese} = \text{overweight} = \text{lean} \\
H_1: \text{obese} \le \text{overweight} \ge \text{lean} \quad \text{with at least one strict inequality}$$
one should set `node` as the position of the turning point of `x`. 
In this example, the turning position is `overweight`, thus, `node = 1`.

We will test both the monotonically increasing and decreasing trends in this 
example. The result contains: 
1) log fold changes, 2) standard errors, 3) test  statistics, 
4) p-values, 5) adjusted p-values, 6) indicators of whether the taxon 
is differentially abundant (TRUE) or not (FALSE).

Note that the LFC and SE for the reference group ("obese" in this example) will
set to be 0s.

```{r, fig.width=10}
res_trend = output$res_trend

df_fig_trend = res_trend %>%
    dplyr::filter(diff_abn == 1, passed_ss == 1) %>%
    dplyr::transmute(taxon,
                     lfc = lfc_bmioverweight,
                     se = se_bmioverweight,
                     q_val,
                     group = "Overweight - Obese") %>%
    dplyr::bind_rows(
        res_trend %>%
            dplyr::filter(diff_abn, passed_ss == 1) %>%
            dplyr::transmute(taxon,
                             lfc = lfc_bmilean,
                             se = se_bmilean,
                             q_val,
                             group = "Lean - Obese")
    )
    
df_fig_trend$group = factor(df_fig_trend$group, 
                           levels = c("Overweight - Obese", "Lean - Obese"))
  
fig_trend = df_fig_trend %>%
  ggplot(aes(x = group, y = lfc, fill = group)) + 
  geom_bar(stat = "identity", position = position_dodge(), color = "black") +
  geom_errorbar(aes(ymin = lfc - se, ymax = lfc + se), width = .2,
                position = position_dodge(.9)) +
  facet_wrap(vars(taxon), nrow = 2, scales = "free") +
  labs(x = NULL, y = NULL, title = "Log fold change as compared to obese subjects") +
  scale_fill_brewer(palette = "Set2", name = NULL) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "bottom")
fig_trend
```

# 5. Run ANCOM-BC2 on a real longitudinal dataset {.tabset}

## 5.1 Import example data

A two-week diet swap study between western (USA) and traditional (rural Africa) 
diets [@lahti2014tipping]. The dataset is available via the 
microbiome R package [@lahti2017tools] in phyloseq [@mcmurdie2013phyloseq] 
format.

```{r}
data(dietswap, package = "microbiome")
tse = mia::makeTreeSummarizedExperimentFromPhyloseq(dietswap)
print(tse)
```

## 5.2 Run ancombc2 function

In this tutorial, we consider the following fixed effects:

* Continuous covariates: "timepoint"

* Categorical covariates: "nationality"

* The group variable of interest: "group"

    + Three groups: "DI", "ED", "HE"
    
    + The reference group: "DI"
    
and the following random effects:

* A random intercept

* A random slope: "timepoint"

Procedures of ANCOM-BC2 global test, pairwise directional test, Dunnett's type 
of test, and trend test are the same as those for the cross-sectional 
data shown above.

```{r}
set.seed(123)
# It should be noted that we have set the number of bootstrap samples (B) equal 
# to 10 in the 'trend_control' function for computational expediency. 
# However, it is recommended that users utilize the default value of B, 
# which is 100, or larger values for optimal performance.
output = ancombc2(data = tse, assay_name = "counts", tax_level = "Family",
                  fix_formula = "nationality + timepoint + group",
                  rand_formula = "(timepoint | subject)",
                  p_adj_method = "holm", 
                  prv_cut = 0.10, lib_cut = 1000, s0_perc = 0.05,
                  group = "group", struc_zero = TRUE, neg_lb = TRUE,
                  alpha = 0.05, n_cl = 2, verbose = TRUE,
                  global = TRUE, pairwise = TRUE, dunnet = TRUE, trend = TRUE,
                  iter_control = list(tol = 1e-2, max_iter = 20, 
                                      verbose = TRUE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
                  trend_control = list(contrast = list(matrix(c(1, 0, -1, 1),
                                                              nrow = 2, 
                                                              byrow = TRUE)),
                                       node = list(2),
                                       solver = "ECOS",
                                       B = 10))

res_prim = output$res %>%
    mutate_if(is.numeric, function(x) round(x, 2))
res_prim %>%
    datatable(caption = "ANCOM-BC2 Primary Results")
```

# 6. Bias-corrected log abundances

It is important to acknowledge that the estimation of sampling fractions in 
ANCOM-BC2 is limited to an additive constant. This means that only the 
difference between bias-corrected log abundances is meaningful, rather than the 
absolute values themselves. 

Furthermore, ANCOM-BC2 does not consider taxon-specific biases when calculating 
the bias-corrected log abundances. The method assumes that these biases vary 
across taxa but remain constant within a taxon across different samples. This 
assumption allows ANCOM-BC2 to account for variation between taxa while focusing 
on identifying differential abundance patterns across samples.

```{r}
bias_correct_log_table = output$bias_correct_log_table
# By default, ANCOM-BC2 does not add pseudo-counts to zero counts, which can 
# result in NAs in the bias-corrected log abundances. Users have the option to 
# either leave the NAs as they are or replace them with zeros. 
# This replacement is equivalent to adding pseudo-counts of ones to the zero counts. 
bias_correct_log_table[is.na(bias_correct_log_table)] = 0
# Show the first 6 samples
round(bias_correct_log_table[, 1:6], 2) %>% 
  datatable(caption = "Bias-corrected log abundances")
```

# Session information

```{r sessionInfo, message = FALSE, warning = FALSE, comment = NA}
sessionInfo()
```

# References