# ANCOM-BC2 Tutorial

#### Huang Lin\(^1\)

#### \(^1\)NICHD, 6710B Rockledge Dr, Bethesda, MD 20892

#### October 18, 2022

# 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 (Lin and Peddada 2020) 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) (Tusher, Tibshirani, and Chu 2001) 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 log transforms the observed counts. However, to deal with zero counts, a pseudo-count is added before the log transformation. Several studies have shown that differential abundance results could be sensitive to the choice of pseudo-count (Costea et al. 2014; Paulson, Bravo, and Pop 2014), resulting in an inflated false positive rate. To avoid such false positives, we conduct a sensitivity analysis and provide a sensitivity score for each taxon to determine if a particular taxon is sensitive to the choice of pseudo-count. The larger the score, the more likely the significant result is a false positive.

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

   - Multiple pairwise directional test: When performning pairwise directional test, 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 (Guo, Sarkar, and Peddada 2010; Grandhi, Guo, and Peddada 2016).

   - Multiple pairwise directional test against a pre-specified group (e.g., Dunnett’s type of test): We use the same set-up as in the multiple pairwise directional test but use the Dunnett-type modification described in (Grandhi, Guo, and Peddada 2016) to control the mdFDR which is more powerful.

   - Trend test 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 (Jelsema and Peddada 2016) to perform trend test 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.

    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("ANCOMBC")

Load the package.

    library(ANCOMBC)

# 3. Run ANCOM-BC2 on a simulated dataset

## 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 (Hu and Satten 2020). 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 (Vandeputte et al. 2017). This dataset contains quantitative microbiome count data of 106 samples and 91 OTUs.

    data(QMP)
    set.seed(12345)
    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, scale_mean = 1e8, scale_sd = 1e8/3)
    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(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)
    rownames(fmd) = paste0("T", seq_len(d))

    # Add effect sizes of covariates to the true abundances
    dmy = caret::dummyVars(" ~ cat_cov", data = smd)
    smd_dmy = data.frame(predict(dmy, newdata = smd))

    log_abn_data = log_abn_data + outer(fmd$lfc_cont, smd$cont_cov)
    log_abn_data = log_abn_data + outer(fmd$lfc_cat2_vs_1, smd_dmy$cat_cov.2)
    log_abn_data = log_abn_data + outer(fmd$lfc_cat3_vs_1, smd_dmy$cat_cov.3)

    # 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 = SimpleList(counts = otu_data)
    smd = DataFrame(smd)
    tse = TreeSummarizedExperiment(assays = assays, colData = smd)

## 3.2 Run ancombc2 function

    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", pseudo = 0, pseudo_sens = TRUE,
                      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 = FALSE,
                      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 = NULL,
                      trend_control = NULL)

    res_prim = output$res
    tab_sens = output$pseudo_sens_tab

## 3.3 Sensitivity scores for “cat_cov”

1. Extreme values (>5) of the sensitivity score are observed.

2. Significant taxa with extreme sensitivity scores are likely false positives.

3. We empirically set the threshold for the sensitivity score as 5. Taxa with sensitivity scores larger than the specified threshold will be considered nonsignificant.

    sens_cat = tab_sens %>%
        transmute(sens_cat = cat_cov2) %>%
        rownames_to_column("tax_id") %>%
        left_join(res_prim %>%
                      rownames_to_column("tax_id") %>%
                      transmute(tax_id, diff_cat = diff_cat_cov2),
                  by = "tax_id") %>%
        mutate(group = "Cat2 vs. Cat1") %>%
        bind_rows(
            tab_sens %>%
                transmute(sens_cat = cat_cov3) %>%
                rownames_to_column("tax_id") %>%
                left_join(res_prim %>%
                              rownames_to_column("tax_id") %>%
                              transmute(tax_id, diff_cat = diff_cat_cov3),
                          by = "tax_id") %>%
                mutate(group = "Cat3 vs. Cat1")
        )
    sens_cat$diff_cat = recode(sens_cat$diff_cat * 1,
                               `1` = "Significant",
                               `0` = "Nonsignificant")

    fig_sens_cat = sens_cat %>%
        ggplot(aes(x = tax_id, y = sens_cat, color = diff_cat)) +
        geom_point() +
        scale_color_brewer(palette = "Dark2", name = NULL) +
        facet_grid(rows = vars(group), scales = "free") +
        labs(x = NULL, y = "Sensitivity Score") +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 60, vjust = 0.5))
    fig_sens_cat

    sens_cut1 = 5
    sens_cut2 = 5

## 3.4 Power and FDR

    # Not considering the effect of pseudo-count addition
    lfc = res_prim %>%
        dplyr::select(starts_with("lfc"))
    diff = res_prim %>%
        dplyr::select(starts_with("diff"))
    res_merge = (lfc * diff) %>%
        rownames_to_column("tax_id") %>%
        left_join(fmd %>%
                      rownames_to_column("tax_id"),
                  by = "tax_id") %>%
        transmute(
            tax_id = tax_id,
            est_cov2 = case_when(
                lfc_cat_cov2 > 0 ~ 1,
                lfc_cat_cov2 < 0 ~ -1,
                TRUE ~ 0),
            est_cov3 = case_when(
                lfc_cat_cov3 > 0 ~ 1,
                lfc_cat_cov3 < 0 ~ -1,
                TRUE ~ 0),
            true_cat2 = case_when(
                lfc_cat2_vs_1 > 0 ~ 1,
                lfc_cat2_vs_1 < 0 ~ -1,
                TRUE ~ 0),
            true_cat3 = case_when(
                lfc_cat3_vs_1 > 0 ~ 1,
                lfc_cat3_vs_1 < 0 ~ -1,
                TRUE ~ 0)
            )

    # Cat 2 vs Cat 1
    lfc_est = res_merge$est_cov2
    lfc_true = res_merge$true_cat2
    tp = sum(lfc_true != 0 & lfc_est != 0)
    fp = sum(lfc_true == 0 & lfc_est != 0)
    fn = sum(lfc_true != 0 & lfc_est == 0)
    power1_nosens = tp/(tp + fn)
    fdr1_nosens = fp/(tp + fp)

    # Cat 3 vs Cat 1
    lfc_est = res_merge$est_cov3
    lfc_true = res_merge$true_cat3
    tp = sum(lfc_true != 0 & lfc_est != 0)
    fp = sum(lfc_true == 0 & lfc_est != 0)
    fn = sum(lfc_true != 0 & lfc_est == 0)
    power2_nosens = tp/(tp + fn)
    fdr2_nosens = fp/(tp + fp)

    # Considering the effect of pseudo-count addition
    res_merge2 = res_merge %>%
        left_join(tab_sens %>%
                      rownames_to_column("tax_id"),
                  by = "tax_id")

    # Cat 2 vs Cat 1
    lfc_est = res_merge2$est_cov2 * (res_merge2$cat_cov2 < sens_cut1)
    lfc_true = res_merge2$true_cat2
    tp = sum(lfc_true != 0 & lfc_est != 0)
    fp = sum(lfc_true == 0 & lfc_est != 0)
    fn = sum(lfc_true != 0 & lfc_est == 0)
    power1_sens = tp/(tp + fn)
    fdr1_sens = fp/(tp + fp)

    # Cat 3 vs Cat 1
    lfc_est = res_merge2$est_cov3 * (res_merge2$cat_cov3 < sens_cut2)
    lfc_true = res_merge2$true_cat3
    tp = sum(lfc_true != 0 & lfc_est != 0)
    fp = sum(lfc_true == 0 & lfc_est != 0)
    fn = sum(lfc_true != 0 & lfc_est == 0)
    power2_sens = tp/(tp + fn)
    fdr2_sens = fp/(tp + fp)

    tab_summ1 = data.frame(Comparison = c("Not considering the effect of pseudo-count addition",
                                          "Considering the effect of pseudo-count addition"),
                           Power = round(c(power1_nosens, power1_sens), 2),
                           FDR = round(c(fdr1_nosens, fdr1_sens), 2))

    tab_summ2 = data.frame(Comparison = c("Not considering the effect of pseudo-count addition",
                                          "Considering the effect of pseudo-count addition"),
                           Power = round(c(power2_nosens, power2_sens), 2),
                           FDR = round(c(fdr2_nosens, fdr2_sens), 2))
    tab_summ1 %>%
        datatable(caption = "Cat 2 vs Cat 1")

    tab_summ2 %>%
        datatable(caption = "Cat 3 vs Cat 1")

# 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 (Lahti et al. 2014). The dataset is also available via the microbiome R package (Lahti et al. 2017) in phyloseq (McMurdie and Holmes 2013) 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”

    data(atlas1006)

    # Subset to baseline
    tse = atlas1006[, atlas1006$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)

    class: TreeSummarizedExperiment
    dim: 130 873
    metadata(0):
    assays(1): counts
    rownames(130): Actinomycetaceae Aerococcus ... Xanthomonadaceae
      Yersinia et rel.
    rowData names(3): Phylum Family Genus
    colnames(873): Sample-1 Sample-2 ... Sample-1005 Sample-1006
    colData names(12): age sex ... bmi region
    reducedDimNames(0):
    mainExpName: NULL
    altExpNames(0):
    rowLinks: NULL
    rowTree: NULL
    colLinks: NULL
    colTree: NULL

## 4.2 Run ancombc2 function

    set.seed(123)
    output = ancombc2(data = tse, assay_name = "counts", tax_level = "Family",
                      fix_formula = "age + region + bmi", rand_formula = NULL,
                      p_adj_method = "holm", pseudo = 0, pseudo_sens = TRUE,
                      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 = 100))

## 4.3 Structural zeros

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

## 4.4 Sensitivity scores

    tab_sens = output$pseudo_sens_tab
    tab_sens %>%
        datatable(caption = "Sensitivity Scores") %>%
        formatRound(colnames(tab_sens), digits = 2)

## 4.5 ANCOM-BC2 primary analysis

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).

    res_prim = output$res

### Results for age

#### Original results

    df_age = res_prim %>%
        rownames_to_column("tax_id") %>%
        dplyr::select(tax_id, ends_with("age"))
    df_fig_age = df_age %>%
        filter(diff_age == 1) %>%
        arrange(desc(lfc_age)) %>%
        mutate(direct = ifelse(lfc_age > 0, "Positive LFC", "Negative LFC"))
    df_fig_age$tax_id = factor(df_fig_age$tax_id, levels = df_fig_age$tax_id)
    df_fig_age$direct = factor(df_fig_age$direct,
                               levels = c("Positive LFC", "Negative LFC"))

    fig_age = df_fig_age %>%
        ggplot(aes(x = tax_id, 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

#### Pseudo-count sensitivity analysis for age

For the covariate of age, no outlying sensitivity scores are observed. All significant taxa have low sensitivity scores.

    sens_age = tab_sens %>%
        transmute(sens_age = age) %>%
        rownames_to_column("tax_id") %>%
        left_join(df_age, by = "tax_id")
    sens_age$diff_age = recode(sens_age$diff_age * 1,
                               `1` = "Significant",
                               `0` = "Nonsignificant")

    fig_sens_age = sens_age %>%
        ggplot(aes(x = tax_id, y = sens_age, color = diff_age)) +
        geom_point() +
        scale_color_brewer(palette = "Dark2", name = NULL) +
        labs(x = NULL, y = "Sensitivity Score") +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 60, vjust = 0.5))
    fig_sens_age

### Results for bmi

#### Original results

    df_bmi = res_prim %>%
        rownames_to_column("tax_id") %>%
        dplyr::select(tax_id, contains("bmi"))
    df_fig_bmi = df_bmi %>%
        filter(diff_bmilean == 1 | diff_bmioverweight == 1) %>%
        mutate(lfc_overweight = ifelse(diff_bmioverweight == 1,
                                       lfc_bmioverweight, 0),
               lfc_lean = ifelse(diff_bmilean == 1,
                                 lfc_bmilean, 0)) %>%
        transmute(tax_id,
                  `Overweight vs. Obese` = round(lfc_overweight, 2),
                  `Lean vs. Obese` = round(lfc_lean, 2)) %>%
        pivot_longer(cols = `Overweight vs. Obese`:`Lean vs. Obese`,
                     names_to = "group", values_to = "value") %>%
        arrange(tax_id)

    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 = tax_id, 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, tax_id, 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

#### Pseudo-count sensitivity analysis for bmi

For the covariate of bmi, no outlying sensitivity scores are observed. All significant taxa have low sensitivity scores.

    sens_bmi = tab_sens %>%
        transmute(sens_bmi = bmilean) %>%
        rownames_to_column("tax_id") %>%
        left_join(df_bmi %>%
                      transmute(tax_id, diff_bmi = diff_bmilean),
                  by = "tax_id") %>%
        mutate(group = "Lean vs. Obese") %>%
        bind_rows(
            tab_sens %>%
                transmute(sens_bmi = bmioverweight) %>%
                rownames_to_column("tax_id") %>%
                left_join(df_bmi %>%
                              transmute(tax_id, diff_bmi = diff_bmioverweight),
                          by = "tax_id") %>%
                mutate(group = "Overweight vs. Obese")
        )
    sens_bmi$diff_bmi = recode(sens_bmi$diff_bmi * 1,
                               `1` = "Significant",
                               `0` = "Nonsignificant")

    fig_sens_bmi = sens_bmi %>%
        ggplot(aes(x = tax_id, y = sens_bmi, color = diff_bmi)) +
        geom_point() +
        scale_color_brewer(palette = "Dark2", name = NULL) +
        facet_grid(rows = vars(group), scales = "free") +
        labs(x = NULL, y = "Sensitivity Score") +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 60, vjust = 0.5))
    fig_sens_bmi

## 4.6 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).

    res_global = output$res_global

### Results for BMI

    df_bmi = res_prim %>%
        rownames_to_column("tax_id") %>%
        dplyr::select(tax_id, contains("bmi"))
    df_fig_global = df_bmi %>%
        left_join(res_global %>%
                      rownames_to_column("tax_id") %>%
                      transmute(tax_id, diff_bmi = diff_abn)) %>%
        dplyr::filter(diff_bmi == 1) %>%
        mutate(lfc_lean = ifelse(diff_bmilean == 1,
                                 lfc_bmilean, 0),
               lfc_overweight = ifelse(diff_bmioverweight == 1,
                                       lfc_bmioverweight, 0)) %>%
        transmute(tax_id,
                  `Lean vs. Obese` = round(lfc_lean, 2),
                  `Overweight vs. Obese` = round(lfc_overweight, 2)) %>%
        pivot_longer(cols = `Lean vs. Obese`:`Overweight vs. Obese`,
                     names_to = "group", values_to = "value") %>%
        arrange(tax_id)

    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 = tax_id, 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, tax_id, 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

### Pseudo-count sensitivity analysis for the global test

For the global test of bmi, no outlying sensitivity scores are observed. All significant taxa have low sensitivity scores.

    sens_global = tab_sens %>%
        transmute(sens_global = global) %>%
        rownames_to_column("tax_id") %>%
        left_join(res_global %>%
                      rownames_to_column("tax_id") %>%
                      transmute(tax_id, diff_global = diff_abn * 1),
                  by = "tax_id")
    sens_global$diff_global = recode(sens_global$diff_global,
                                     `1` = "Significant",
                                     `0` = "Nonsignificant")

    fig_sens_global = sens_global %>%
        ggplot(aes(x = tax_id, y = sens_global, color = diff_global)) +
        geom_point() +
        scale_color_brewer(palette = "Dark2", name = NULL) +
        labs(x = NULL, y = "Sensitivity Score") +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 60, vjust = 0.5))
    fig_sens_global

## 4.7 ANCOM-BC2 pairwise directional test

ANCOM-BC2 pairwise directional test aims 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).

    res_pair = output$res_pair

### Results for BMI

    df_pair = res_pair %>%
        rownames_to_column("tax_id")
    df_fig_pair = df_pair %>%
        filter(diff_bmilean == 1 | diff_bmioverweight == 1 |
                   diff_bmilean_bmioverweight == 1) %>%
        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)) %>%
        transmute(tax_id,
                  `Lean vs. Obese` = round(lfc_lean, 2),
                  `Overweight vs. Obese` = round(lfc_overweight, 2),
                  `Lean vs. Overweight` = round(lfc_lean_overweight, 2)
                  ) %>%
        pivot_longer(cols = `Lean vs. Obese`:`Lean vs. Overweight`,
                     names_to = "group", values_to = "value") %>%
        arrange(tax_id)
    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 = tax_id, 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, tax_id, 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

### Pseudo-count sensitivity analysis for bmi

For the pairwise analysis of bmi, no outlying sensitivity scores are observed. All significant taxa have low sensitivity scores.

    sens_pair = tab_sens %>%
        transmute(sens_pair = `obese - lean`) %>%
        rownames_to_column("tax_id") %>%
        left_join(df_pair %>%
                      transmute(tax_id,
                                diff_pair = diff_bmilean),
                  by = "tax_id") %>%
        mutate(group = "Lean vs. Obese") %>%
        bind_rows(
            tab_sens %>%
                transmute(sens_pair = `obese - overweight`) %>%
                rownames_to_column("tax_id") %>%
                left_join(df_pair %>%
                              transmute(tax_id,
                                        diff_pair = diff_bmioverweight),
                          by = "tax_id") %>%
                mutate(group = "Overweight vs. Obese")
        ) %>%
        bind_rows(
            tab_sens %>%
                transmute(sens_pair = `overweight - lean`) %>%
                rownames_to_column("tax_id") %>%
                left_join(df_pair %>%
                              transmute(tax_id,
                                        diff_pair = diff_bmilean_bmioverweight),
                          by = "tax_id") %>%
                mutate(group = "Lean vs. Overweight")
        )
    sens_pair$diff_pair = recode(sens_pair$diff_pair * 1,
                                 `1` = "Significant",
                                 `0` = "Nonsignificant")
    sens_pair$group = factor(sens_pair$group,
                             levels = c("Lean vs. Obese",
                                        "Overweight vs. Obese",
                                        "Lean vs. Overweight"))

    fig_sens_pair = sens_pair %>%
        ggplot(aes(x = tax_id, y = sens_pair, color = diff_pair)) +
        geom_point() +
        scale_color_brewer(palette = "Dark2", name = NULL) +
        facet_grid(rows = vars(group), scales = "free") +
        labs(x = NULL, y = "Sensitivity Score") +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 60, vjust = 0.5))
    fig_sens_pair

## 4.8 ANCOM-BC2 Dunnett’s type of test

The Dunnett’s test (Dunnett 1955; Dunnett and Tamhane 1991, 1992) 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).

    res_dunn = output$res_dunn

### Results for BMI

    df_dunn = res_dunn %>%
        rownames_to_column("tax_id")
    df_fig_dunn = df_dunn %>%
        dplyr::filter(diff_bmilean == 1 | diff_bmioverweight == 1) %>%
        mutate(lfc_lean = ifelse(diff_bmilean == 1,
                                 lfc_bmilean, 0),
               lfc_overweight = ifelse(diff_bmioverweight == 1,
                                       lfc_bmioverweight, 0)) %>%
        transmute(tax_id,
                  `Lean vs. Obese` = round(lfc_lean, 2),
                  `Overweight vs. Obese` = round(lfc_overweight, 2)) %>%
        pivot_longer(cols = `Lean vs. Obese`:`Overweight vs. Obese`,
                     names_to = "group", values_to = "value") %>%
        arrange(tax_id)

    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 = tax_id, 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, tax_id, 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

### Pseudo-count sensitivity analysis for bmi

For the covariate of bmi, no outlying sensitivity scores are observed. All significant taxa have low sensitivity scores.

    sens_dunn = tab_sens %>%
        transmute(sens_bmi = bmilean) %>%
        rownames_to_column("tax_id") %>%
        left_join(df_dunn %>%
                      transmute(tax_id, diff_bmi = diff_bmilean),
                  by = "tax_id") %>%
        mutate(group = "Lean vs. Obese") %>%
        bind_rows(
            tab_sens %>%
                transmute(sens_bmi = bmioverweight) %>%
                rownames_to_column("tax_id") %>%
                left_join(df_dunn %>%
                              transmute(tax_id, diff_bmi = diff_bmioverweight),
                          by = "tax_id") %>%
                mutate(group = "Overweight vs. Obese")
        )
    sens_dunn$diff_bmi = recode(sens_dunn$diff_bmi * 1,
                                `1` = "Significant",
                                `0` = "Nonsignificant")

    fig_sens_dunn = sens_dunn %>%
        ggplot(aes(x = tax_id, y = sens_bmi, color = diff_bmi)) +
        geom_point() +
        scale_color_brewer(palette = "Dark2", name = NULL) +
        facet_grid(rows = vars(group), scales = "free") +
        labs(x = NULL, y = "Sensitivity Score") +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 60, vjust = 0.5))
    fig_sens_dunn

## 4.9 ANCOM-BC2 trend test

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 trends accordingly. Examples of trends include: monotonically increasing, monotonically decreasing, and umbrella shape.

ANCOM-BC2 is able to identify potential trends 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).

    res_trend = output$res_trend

### Results for BMI

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

    df_trend = res_trend %>%
        rownames_to_column("tax_id")
    df_fig_trend = df_trend %>%
        dplyr::filter(diff_abn) %>%
        transmute(tax_id,
                  lfc = lfc_bmioverweight,
                  se = se_bmioverweight,
                  q_val,
                  group = "Overweight - Obese") %>%
        bind_rows(
            df_trend %>%
                dplyr::filter(diff_abn) %>%
                transmute(tax_id,
                          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(tax_id), 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 = c(0.9, 0.2))
    fig_trend

### Pseudo-count sensitivity analysis for bmi

For the covariate of bmi, no outlying sensitivity scores are observed. All significant taxa have low sensitivity scores.

    sens_trend = tab_sens %>%
        transmute(sens_bmi = bmilean) %>%
        rownames_to_column("tax_id") %>%
        left_join(df_trend %>%
                      transmute(tax_id, diff_bmi = diff_abn),
                  by = "tax_id") %>%
        mutate(group = "Lean vs. Obese") %>%
        bind_rows(
            tab_sens %>%
                transmute(sens_bmi = bmioverweight) %>%
                rownames_to_column("tax_id") %>%
                left_join(df_trend %>%
                              transmute(tax_id, diff_bmi = diff_abn),
                          by = "tax_id") %>%
                mutate(group = "Overweight vs. Obese")
        )
    sens_trend$diff_bmi = recode(sens_trend$diff_bmi * 1,
                                `1` = "Significant",
                                `0` = "Nonsignificant")

    fig_sens_trend = sens_trend %>%
        ggplot(aes(x = tax_id, y = sens_bmi, color = diff_bmi)) +
        geom_point() +
        scale_color_brewer(palette = "Dark2", name = NULL) +
        facet_grid(rows = vars(group), scales = "free") +
        labs(x = NULL, y = "Sensitivity Score") +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 60, vjust = 0.5))
    fig_sens_trend

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

## 5. 1 Import example data

A two-week diet swap study between western (USA) and traditional (rural Africa) diets (Lahti et al. 2014). The dataset is also available via the microbiome R package (Lahti et al. 2017) in phyloseq (McMurdie and Holmes 2013) format.

    data(dietswap)
    tse = dietswap
    print(tse)

    class: TreeSummarizedExperiment
    dim: 130 222
    metadata(0):
    assays(1): counts
    rownames(130): Actinomycetaceae Aerococcus ... Xanthomonadaceae
      Yersinia et rel.
    rowData names(3): Phylum Family Genus
    colnames(222): Sample-1 Sample-2 ... Sample-221 Sample-222
    colData names(8): subject sex ... timepoint.within.group bmi_group
    reducedDimNames(0):
    mainExpName: NULL
    altExpNames(0):
    rowLinks: NULL
    rowTree: NULL
    colLinks: NULL
    colTree: NULL

## 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.

    set.seed(123)
    output = ancombc2(data = tse, assay_name = "counts", tax_level = "Family",
                      fix_formula = "nationality + timepoint + group",
                      rand_formula = "(timepoint | subject)",
                      p_adj_method = "holm", pseudo = 0, pseudo_sens = TRUE,
                      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 = 100))

# Session information

    sessionInfo()

    R version 4.2.1 (2022-06-23 ucrt)
    Platform: x86_64-w64-mingw32/x64 (64-bit)
    Running under: Windows Server x64 (build 20348)

    Matrix products: default

    locale:
    [1] LC_COLLATE=C
    [2] LC_CTYPE=English_United States.utf8
    [3] LC_MONETARY=English_United States.utf8
    [4] LC_NUMERIC=C
    [5] LC_TIME=English_United States.utf8

    attached base packages:
    [1] stats4    stats     graphics  grDevices utils     datasets  methods
    [8] base

    other attached packages:
     [1] caret_6.0-93                   lattice_0.20-45
     [3] DT_0.25                        mia_1.4.0
     [5] MultiAssayExperiment_1.22.0    TreeSummarizedExperiment_2.4.0
     [7] Biostrings_2.64.1              XVector_0.36.0
     [9] SingleCellExperiment_1.18.1    SummarizedExperiment_1.26.1
    [11] Biobase_2.56.0                 GenomicRanges_1.48.0
    [13] GenomeInfoDb_1.32.4            IRanges_2.30.1
    [15] S4Vectors_0.34.0               BiocGenerics_0.42.0
    [17] MatrixGenerics_1.8.1           matrixStats_0.62.0
    [19] forcats_0.5.2                  stringr_1.4.1
    [21] dplyr_1.0.10                   purrr_0.3.5
    [23] readr_2.1.3                    tidyr_1.2.1
    [25] tibble_3.1.8                   ggplot2_3.3.6
    [27] tidyverse_1.3.2                ANCOMBC_1.6.4

    loaded via a namespace (and not attached):
      [1] estimability_1.4.1          ModelMetrics_1.2.2.2
      [3] coda_0.19-4                 bit64_4.0.5
      [5] knitr_1.40                  irlba_2.3.5.1
      [7] multcomp_1.4-20             DelayedArray_0.22.0
      [9] data.table_1.14.4           rpart_4.1.16
     [11] hardhat_1.2.0               RCurl_1.98-1.9
     [13] doParallel_1.0.17           generics_0.1.3
     [15] ScaledMatrix_1.4.1          TH.data_1.1-1
     [17] RSQLite_2.2.18              proxy_0.4-27
     [19] future_1.28.0               bit_4.0.4
     [21] tzdb_0.3.0                  xml2_1.3.3
     [23] lubridate_1.8.0             assertthat_0.2.1
     [25] DirichletMultinomial_1.38.0 viridis_0.6.2
     [27] gargle_1.2.1                gower_1.0.0
     [29] xfun_0.33                   hms_1.1.2
     [31] jquerylib_0.1.4             evaluate_0.17
     [33] fansi_1.0.3                 dbplyr_2.2.1
     [35] readxl_1.4.1                DBI_1.1.3
     [37] htmlwidgets_1.5.4           googledrive_2.0.0
     [39] Rmpfr_0.8-9                 CVXR_1.0-10
     [41] ellipsis_0.3.2              crosstalk_1.2.0
     [43] backports_1.4.1             energy_1.7-10
     [45] permute_0.9-7               deldir_1.0-6
     [47] sparseMatrixStats_1.8.0     vctrs_0.4.2
     [49] cachem_1.0.6                withr_2.5.0
     [51] checkmate_2.1.0             emmeans_1.8.1-1
     [53] vegan_2.6-4                 treeio_1.20.2
     [55] cluster_2.1.4               gsl_2.1-7.1
     [57] ape_5.6-2                   lazyeval_0.2.2
     [59] crayon_1.5.2                recipes_1.0.2
     [61] pkgconfig_2.0.3             labeling_0.4.2
     [63] nlme_3.1-160                vipor_0.4.5
     [65] nnet_7.3-18                 rlang_1.0.6
     [67] globals_0.16.1              lifecycle_1.0.3
     [69] sandwich_3.0-2              modelr_0.1.9
     [71] rsvd_1.0.5                  cellranger_1.1.0
     [73] rngtools_1.5.2              Matrix_1.5-1
     [75] boot_1.3-28                 zoo_1.8-11
     [77] reprex_2.0.2                base64enc_0.1-3
     [79] beeswarm_0.4.0              googlesheets4_1.0.1
     [81] png_0.1-7                   viridisLite_0.4.1
     [83] rootSolve_1.8.2.3           bitops_1.0-7
     [85] pROC_1.18.0                 blob_1.2.3
     [87] DelayedMatrixStats_1.18.2   doRNG_1.8.2
     [89] decontam_1.16.0             parallelly_1.32.1
     [91] jpeg_0.1-9                  DECIPHER_2.24.0
     [93] beachmat_2.12.0             scales_1.2.1
     [95] memoise_2.0.1               magrittr_2.0.3
     [97] plyr_1.8.7                  zlibbioc_1.42.0
     [99] compiler_4.2.1              RColorBrewer_1.1-3
    [101] lme4_1.1-30                 cli_3.4.1
    [103] lmerTest_3.1-3              listenv_0.8.0
    [105] htmlTable_2.4.1             Formula_1.2-4
    [107] MASS_7.3-58.1               mgcv_1.8-40
    [109] tidyselect_1.2.0            stringi_1.7.8
    [111] highr_0.9                   yaml_2.3.6
    [113] BiocSingular_1.12.0         latticeExtra_0.6-30
    [115] ggrepel_0.9.1               grid_4.2.1
    [117] sass_0.4.2                  tools_4.2.1
    [119] lmom_2.9                    future.apply_1.9.1
    [121] parallel_4.2.1              rstudioapi_0.14
    [123] foreach_1.5.2               foreign_0.8-83
    [125] gridExtra_2.3               gld_2.6.5
    [127] prodlim_2019.11.13          farver_2.1.1
    [129] digest_0.6.30               lava_1.6.10
    [131] Rcpp_1.0.9                  broom_1.0.1
    [133] scuttle_1.6.3               httr_1.4.4
    [135] Rdpack_2.4                  colorspace_2.0-3
    [137] rvest_1.0.3                 fs_1.5.2
    [139] splines_4.2.1               yulab.utils_0.0.5
    [141] tidytree_0.4.1              expm_0.999-6
    [143] scater_1.24.0               Exact_3.2
    [145] xtable_1.8-4                gmp_0.6-6
    [147] jsonlite_1.8.2              nloptr_2.0.3
    [149] timeDate_4021.106           ipred_0.9-13
    [151] R6_2.5.1                    Hmisc_4.7-1
    [153] pillar_1.8.1                htmltools_0.5.3
    [155] glue_1.6.2                  fastmap_1.1.0
    [157] minqa_1.2.4                 BiocParallel_1.30.4
    [159] BiocNeighbors_1.14.0        class_7.3-20
    [161] codetools_0.2-18            mvtnorm_1.1-3
    [163] utf8_1.2.2                  bslib_0.4.0
    [165] numDeriv_2016.8-1.1         ggbeeswarm_0.6.0
    [167] DescTools_0.99.46           interp_1.1-3
    [169] survival_3.4-0              rmarkdown_2.17
    [171] munsell_0.5.0               e1071_1.7-11
    [173] GenomeInfoDbData_1.2.8      iterators_1.0.14
    [175] haven_2.5.1                 reshape2_1.4.4
    [177] gtable_0.3.1                rbibutils_2.2.9

# References

Costea, Paul I, Georg Zeller, Shinichi Sunagawa, and Peer Bork. 2014. “A Fair Comparison.” Nature Methods 11 (4): 359–59.

Dunnett, Charles W. 1955. “A Multiple Comparison Procedure for Comparing Several Treatments with a Control.” Journal of the American Statistical Association 50 (272): 1096–1121.

Dunnett, Charles W, and Ajit C Tamhane. 1991. “Step-down Multiple Tests for Comparing Treatments with a Control in Unbalanced One-Way Layouts.” Statistics in Medicine 10 (6): 939–47.

———. 1992. “A Step-up Multiple Test Procedure.” Journal of the American Statistical Association 87 (417): 162–70.

Grandhi, Anjana, Wenge Guo, and Shyamal D Peddada. 2016. “A Multiple Testing Procedure for Multi-Dimensional Pairwise Comparisons with Application to Gene Expression Studies.” BMC Bioinformatics 17 (1): 1–12.

Guo, Wenge, Sanat K Sarkar, and Shyamal D Peddada. 2010. “Controlling False Discoveries in Multidimensional Directional Decisions, with Applications to Gene Expression Data on Ordered Categories.” Biometrics 66 (2): 485–92.

Hu, Yi-Juan, and Glen A Satten. 2020. “Testing Hypotheses About the Microbiome Using the Linear Decomposition Model (Ldm).” Bioinformatics 36 (14): 4106–15.

Jelsema, Casey M, and Shyamal D Peddada. 2016. “CLME: An R Package for Linear Mixed Effects Models Under Inequality Constraints.” Journal of Statistical Software 75.

Lahti, Leo, Jarkko Salojärvi, Anne Salonen, Marten Scheffer, and Willem M De Vos. 2014. “Tipping Elements in the Human Intestinal Ecosystem.” Nature Communications 5 (1): 1–10.

Lahti, Leo, Sudarshan Shetty, T Blake, J Salojarvi, and others. 2017. “Tools for Microbiome Analysis in R.” Version 1: 10013.

Lin, Huang, and Shyamal Das Peddada. 2020. “Analysis of Compositions of Microbiomes with Bias Correction.” Nature Communications 11 (1): 1–11.

McMurdie, Paul J, and Susan Holmes. 2013. “Phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data.” PloS One 8 (4): e61217.

Paulson, Joseph N, Héctor Corrada Bravo, and Mihai Pop. 2014. “Reply to:" A Fair Comparison".” Nature Methods 11 (4): 359–60.

Tusher, Virginia Goss, Robert Tibshirani, and Gilbert Chu. 2001. “Significance Analysis of Microarrays Applied to the Ionizing Radiation Response.” Proceedings of the National Academy of Sciences 98 (9): 5116–21.

Vandeputte, Doris, Gunter Kathagen, Kevin D’hoe, Sara Vieira-Silva, Mireia Valles-Colomer, João Sabino, Jun Wang, et al. 2017. “Quantitative Microbiome Profiling Links Gut Community Variation to Microbial Load.” Nature 551 (7681): 507–11.