StackbarExtended: a user-friendly stacked bar-plot representation incorporating phylogenetic information and microbiota differential abundance analysis

Operation

StackbarExtended is implemented in native R, required version (≥ 4.0), and is platform independent, with its source code available on GitHub and archived on Zenodo (Cuisiniere, 2024a).

Implementation

StackbarExtended is designed for microbiota data analysis and visualization. It requires three main inputs: a taxonomic table providing phylogenetic information, a count table detailing the abundance of the taxa, and metadata containing biological or experimental group data (Navgire et al., 2022). These tables can be obtained through established pipelines such as DADA2 (Callahan et al., 2016), Mothur (Schloss et al., 2009), and QIIME2 (Bolyen et al., 2019) in conjunction with reference databases such as Greengenes (DeSantis et al., 2006) or SILVA (Quast et al., 2013). The package supports a phyloseq S4 class object as an input, promoting a data management approach that regroups the count, taxonomy, and metadata tables into a single entity, thus simplifying data handling and enhancing reproducibility. The Phyloseq S4 object can be obtained through the phyloseq() function of the phyloseq R package (McMurdie and Holmes, 2013).

One of the StackbarExtended functionalities is the ability to handle multiple taxonomic levels simultaneously, typically focusing on user-defined levels set by default to Family (X) and Phylum (X+1). This is achieved through the use of the tax_glom() function from the phyloseq package. Most importantly, StackbarExtended uses a color-coding mechanism where taxa at a specified level (X) are visually differentiated using user-defined color palettes that match those of their Phylum level (X+1), making visual identification of the taxonomic hierarchies straightforward. This is achieved by applying different shades of the same color to represent individual taxa (X) within the same Phylum (X+1), creating a gradient effect that maps the phylogenetic taxonomy.

Another challenge in analyzing microbiota data is the large number of taxa present in the dataset. In microbiota communities, a small number of taxa typically dominate, while numerous others are present at significantly lower abundance levels (Neu et al., 2021). In classic stacked bar-plot representations such as the ones included in Phyloseq (McMurdie and Holmes, 2013), MicrobiomeAnalyst (Chong et al., 2020) or MicroEco (Liu et al., 2021) pipelines, these low-abundance taxa are simply filtered out before plotting. StackbarExtended allows users to classify taxa based on their relative abundance, and, by default, the package plots only taxa representing more than 1% of the total abundance. In addition, low-abundance taxa (X) can be grouped into their respective “Others” taxonomic level and the information about their taxonomy (i.e. Phylum) is kept while taxa belonging to low abundance Phyla (X+1) are grouped into a general “Others” category. Thus, by still including the low-abundance taxa, StackbarExtended provides a more accurate representation of the taxa in a microbiota community.

Finally, StackbarExtended includes DESeq2 differential abundance analysis functionality (Love et al., 2014) which allows users to statistically assess the difference in taxa abundance between experimental groups and apply fdr corrections. Significant features (i.e. taxa with significant differences in abundance between 2 groups) are highlighted in the legend using bold font, and the significance levels (p-value) after fdr corrections are automatically shown using stars in the legend. When more than two treatment groups are compared, multiple pairwise comparisons can be computed and data-frames are produced with results providing information about the taxa identified as significant at each phylogenetic level analyzed through the DESeq2 analysis. This information includes the taxa names and phylogeny, their corresponding abundance levels, statistical metrics such as log2 fold-change, p-values and fdr-corrected p-values.

Use cases

To demonstrate the use of this package, we have provided ready-to-use example data “ps” from our previously published study (Cuisiniere, 2024a; Cuisiniere et al., 2021) which assessed shifts in mouse gut microbiota composition after antibiotic treatment (Figure 1).

Figure 1. StackbarExtended graphical output.

Mice (n = 9) received oral antibiotic treatment with neomycin and metronidazole for one week. Fecal samples were collected before (Day 0) and after antibiotic treatment (Day 7). 16S rRNA of DNA extracted from fecal samples was sequenced using the Illumina MiSeq platform. The 4 most abundant phyla and families are represented. Families with a mean abundance lower than 1% across the samples are regrouped into “Others”. Taxa represented in bold within the legends are statistically significant after fdr correction (*P < 0.05, **P < 0.01, ***P < 0.001). Data used are from the “ps” dataset and are included in the StackbarExtended R package.

Animal experiments were approved by the Institutional Animal Care committee of the Centre de recherche du Centre hospitalier de l’Université de Montréal (CRCHUM) in agreement with the guidelines of the Canadian Council of Animal Care. The study was carried out in compliance with the ARRIVE guidelines (Cuisiniere, 2024b). No criteria were set for including and excluding animals during the experiment or data points during the analysis. No exclusions of animals, experimental units, or data points were applied for the analysis.

In order to obtain the data, four-week-old female C57Bl/6 mice were purchased from Charles River Laboratories (Saint-Constant, QC, Canada). Constant efforts were made to minimize the suffering of the animals. Nine mice were kept under controlled specific pathogen free (SPF) conditions in the CRCHUM animal facility at a temperature of 22°C, 45-60% humidity with a light-dark cycle of 12-12. They were housed at three mice per cage with ad libitum access to chow and water. Cages were enriched with nesting material. Mice were allowed one week of acclimation following arrival to the CRCHUM animal facility, after which oral antibiotics consisting of metronidazole (1 mg.ml⁻¹, Hospira, St-Laurent, QC, Canada) and neomycin (1 mg.ml⁻¹, Sigma, St-Louis, MO, USA) were added to the drinking water for one week. Fecal samples were collected before (Day 0) and after antibiotic treatment (Day 7), snap-frozen and stored at -80°C. Mice were then euthanized using CO₂ followed by cervical dislocation. DNA was extracted using the Qiagen DNeasy PowerSoil^® kit (Qiagen, Toronto, ON) and quantified using a spectrophotometer (DeNovix DS-11 FX, Wilmington, DE). The 16S ribosomal RNA (rRNA) library preparation and sequencing was performed using the Illumina MiSeq platform at Genome Québec targeting the V3-V4 (Primers: 341F, 805R) region of the 16S rRNA gene. Forward and reverse, raw, demultiplexed 16S rRNA reads were denoised, chimera filtered, and clustered into sequence variants using the Dada2 package (version 1.16) (Callahan et al., 2016) in R (version 4.0.1). Reads were trimmed at the first instance of a quality score less than or equal to 2 or removed if they contained ambiguous nucleotides (N) or if two or more errors were expected based on the quality of the trimmed read. After taxonomic assignment using Silva training set v132 (Quast et al., 2013), ASV (Table 1), taxonomy (Table 2) and metadata (Table 3) were combined into a phyloseq object (McMurdie and Holmes, 2013).

Users can then use the following code to create graphical representation of the gut microbiota composition at the phylum and family levels comparing mice before and after antibiotic treatment and performing differential abundance with fdr correction.

# The plot and the output tables are stored into a single list object

my_plot <- plot_microbiota(
  ps_object = ps,
  exp_group = 'timepoint',
  sample_name = 'SampleID',
  hues = c("Purples", "Blues", "Greens", "Oranges"),
  differential_analysis = T,
  sig_lab = T,
  fdr_threshold = 0.05
)

print (my_plot$plot)

In addition to the graphical representation, two output tables (one for each level) are created containing statistical information concerning the differentially abundant taxa (Table 4). The tables can be accessed as follows:

#Display statistically significant taxa at the phylum level
print(my_plot$significant_table_main)

#Display statistically significant taxa at the family level
print(my_plot$significant_table_sub)

The primary objective of our example dataset is to identify which taxa at the Family and Phylum levels were impacted by antibiotic treatment in mice gut microbiota. This analysis is crucial for understanding how antibiotic interventions alter microbial communities. The graphical output (Figure 1) provided a clear assessment, revealing that among the four most abundant Phyla, three were significantly affected by the antibiotic treatment. Specifically, significant alterations were observed in Verrucomicrobia, Bacteroidetes, and Proteobacteria (fdr < 0.05), indicating a substantial shift in the microbial composition due to antibiotic exposure. Among the abundant Families (>1% of the total abundance) within these Phyla, 10 out of 11 exhibited significant alterations (fdr < 0.05). On average, the six remaining low-abundance Phyla, grouped into the “Others” category, accounted for 2.5% of the total relative abundance. Similarly, within the abundant Phyla, the 14 low-abundance Families accounted for 2.4% of the total relative abundance.

The two generated output tables provide detailed statistical information about the differentially abundant taxa at each level, including log2 fold-change, exact p-values, and fdr-corrected p-values. This comprehensive data allows for identification of taxa significantly impacted by the treatment. Notably, four additional low-abundance Phyla - Cyanobacteria, Deferribacteres, Tenericutes, and Actinobacteria - were found to be significantly affected (fdr < 0.05). Furthermore, a total of 24 families were identified as statistically different (fdr < 0.05) (Cuisiniere, 2024a).

These findings are consistent with previous studies that have reported significant shifts in gut microbiota composition following antibiotic treatment (Fishbein et al., 2023). Hence, this use case demonstrated the advantage of using StackbarExtended to present a clear and interpretable graphical representation while retaining the capacity to perform detailed statistical analysis.