Perhaps the most important thing in microbiome research is the ability to look into the (relative) abundance of taxa within and across samples. It provides a first impression of which taxa (functions) are most prominent within a sample, group or population, and can guide us into certain directions of interests. Due to the compositional structure of the data in microbiome research, one tends to prefer looking into relative abundances rather than absolute abundances. An exploration of the (relative) abundances involves no complex statistical modelling, and can be easily done by means of some descriptive statistics and a visual representation of the data.

Visualization—Stacked or regular bar-charts seem to be the most prevalent visual encodings to do so, although they are limited in the number of species (functions) they can visualize for the chart to still be readable (Knaflic, 2015). Heatmaps are a frequently used alternative that allow us to visualize all species (functions) at once. The use of color intensity as a channel in heatmaps on the other hand makes the comparison in terms of relative abundance a bit harder than using length (bars) (Munzner, 2014). Nonetheless, does the use of color allows us to easily include (relative) abundance visualization in other microbiome aspect oriented visualizations [e.g., alongside taxonomic classification (Ondov et al., 2011)]. Other alternative visual encodings found in literature include the use of angle [e.g., sunburst chart (Ondov et al., 2011)] and area [e.g., bubble plot (Dussud et al., 2018)] to display (relative) abundance. An overview of how visualization is been used to represent (relative) abundance in literature is shown in Figure 3.

Microbiome analysis can be done up to different levels depending on the interest of the study, and the sequencing process used to sample the data. In general, sequencing up to a deeper level provides more detailed information. On the other hand, does it bring more problems into the analysis due to sparseness. Most statistical models are not suited to handle many zero counts in the data (Knight et al., 2018).

Visualization—In the analysis of microbiome samples, it can be interesting to visually represent the hierarchical level of the taxonomies (domain, kingdom, phylum, class, order, family, genus, species), hierarchical level of the functions (category e.g., metabolism, superpathway e.g., carbohydrate metabolism, subpathway e.g., glycolysis), or even the phylogenetic relationship of the species. Tree structures (including radial trees, cladograms, etc.) are the typical visual encodings used, and are basically the only visual encoding found in literature (Figures 3B,E–G).

Alpha diversity provides an idea of the diversity of species within a particular sample. This metric is often used as a biomarker (Prehn-Kristensen et al., 2018) in disease association studies, but also as a check of sample quality (Schloss et al., 2009).

Analysis—Looking into alpha diversity calculations and visual representations, no clear evolution could be found. Many different options exist and are used, but no uniform standard has emerged yet. Typically, alpha diversity metrics can be distinguished into two types: richness- and evenness-measures; Chao1 being the most used richness metric, and Shannon the most used evenness metric. A full list of alpha diversity measures is provided by Hagerty et al. (2020). The authors advocate for the use of a composite metric based on exploratory factor analysis (EFA), taking into account both richness and evenness metrics unified in one.

Visualization—box-plots are widely used to display alpha diversity if the objective is to make a comparison between sample cohorts. Line-charts (rarefaction curves) and scatter-plots tend to be used more frequently when visualizing the metrics across samples; the rarefaction curve presenting the (predicted) sample richness by sequence size, often used for re-sampling. Venn diagrams are used to display which part of the microbial taxa are present in multiple samples in relation to the total diversity within those samples. An overview of the visuals used to represent the within sample diversity is given in Figure 4A–C.

Beta diversity represents the diversity of species across samples, commonly used to find clusters of similar samples. Typically, this feature is calculated in the exploratory analysis, as it provides a first impression on which taxa are important to distinguish samples, but also on how microbial compositions are related to environmental and personal meta data. With regard to the research objectives listed above, social contact networks could for instance be revealed based on similar microbiome compositions of the skin.

Analysis—Beta diversity is expressed as a distance matrix calculation on relative OTU abundance, which serves as an input for visual exploration of sample divergence and similarity. Often occurring distance metrics are: (weighted) UniFrac, Jaccard, Bray-Curtis and Jenson-Shannon (Oliveira et al., 2018; Chong et al., 2020; Shamsaddini et al., 2020). An important note however is that none of these measures account for the compositionality of the data. Compositional replacements for these distance metrics have been developed; philr (Silverman et al., 2017) as a replacement for (weighted) UniFrac, and Aitchison distance (Aitchison et al., 2000) for Jensen-Shannon divergence and the Bray-Curtis dissimilarity metrics. Nevertheless, implementation is lacking in the microbiome visual analysis tools.

From 2019 onward, a new trend seemed to develop, which is to test for statistical significance of the between-sample differences (ordination measures). Statistical tests used for this include AMOVA, HOMOVA, ANOSIM, PERMANOVA, PERMDISP, and LIBSSHUFF (Buza et al., 2019; Chong et al., 2020; Shamsaddini et al., 2020). One important recent development is that ordination analysis techniques can be performed on sample functional potentials rather than their taxonomic proportions (Nagpal et al., 2019).

Visualization—The visual representation of beta diversity can be either directly through heatmaps of the distance matrix (Lei et al., 2017), through ordination based methods (e.g., PCoA, NMDS) which present the samples in a 2 or 3 dimensional space using dimensionality reduction techniques (Vázquez-Baeza et al., 2013; Wang et al., 2016; Bolyen et al., 2019), or by means of network visualizations based on topological data analysis (TDA) (Liao et al., 2019) or cut-off based edges (McMurdie and Holmes, 2013). Note that because of the compositional ignorance in the commonly used distance metrics, samples will be almost exclusively discriminated based on the features that are most abundant realtive to the others features and not on the most variable ones between samples. Therefore, sample location could vary a lot in ordination plots when different features are included or excluded (Gloor et al., 2017). An example of the visual encodings listed above is shown in Figure 4D–F.

With differential abundance analysis, OTUs that differ significantly between samples, cohorts or populations are identified using statistical hypothesis testing. In doing so, taxa can be related to a certain response (e.g., disease state, growth process).

Analysis—The search for the ideal analysis method for differential abundance is still ongoing (Hawinkel et al., 2019). To date, it has been proven that distributional assumptions do not hold for the majority of the taxa, leading to poor performance of parametric models (Hawinkel et al., 2020). The problem with non parametric rank alternatives such as Wilcoxon is that they are typically less powerful in comparison to parametric tests due to their vulnerability to ties in the data (Jonsson et al., 2016). Custom methods have been developed to test on significant differences between microbiome data, taking the compositionality of the data into account (e.g., ANCOM, ALDEx2) (Gloor et al., 2017). In comparison to the complete lack of awareness in Beta diversity analyses, differential relative abundance analysis methods relying on these compositional assumptions are present in some visual analysis tools (Zakrzewski et al., 2017). Yet, another possible solution lies in semiparametric models, such as Probabilistic Index Models (PIM) (Thas et al., 2012). These are based on rank tests (non parametric), but allow for estimates of effect sizes and inclusion of continuous covariates. So far, they haven’t been introduced in microbiome visual analysis tools in a significant way. An important note that came up during one of the workshops, is that the methods used in visual analysis tools are all limited to cross sectional analysis. To the awareness of the expert panel, methods that do allow differential abundance testing in longitudinal studies are sparse, and mostly parametric. Besides, with the currently offered methods, conclusions can only be drawn about associations between taxa and meta data identifying sample cohorts, whereas inference on causality would be of major interest. In recent years, several methods have been proposed relying on structural equation models to reveal the direct and mediation effect of the microbiome on a certain response (Sohn and Li, 2019; Wang et al., 2020). These however cannot be found in the current visual analysis tools. Nonetheless, these methods suffer from validity issues (Vanderweele and Vansteelandt, 2009).

Visualization—To visualize statistical significance, several visual encodings have been used; ranging from simple heatmaps and box-plots, to more complex visuals like the Manhattan plot (Harris et al., 2015), rocky mountain plot (Carpenter et al., 2021), volcano plot (Shamsaddini et al., 2020) or heat tree (Foster et al., 2017). An overview of some of the visualizations found in literature is given in Figures 5A,C,E,F.

Biomarker discovery focuses on finding specific parameters or indicators, called biomarkers, that can be related (assigned) to a particular condition (disease).

Analysis—When it comes to biomarker discovery, two schools of thought can be distinguished: one using predictive models such as machine learning classifiers, and the other based on hypothesis testing. Among the predictive models, LEfSE (Swenson and Swenson, 2014) is by far the most offered method in the visual analysis tools, followed by some other machine learning algorithms. Methods based on hypothesis testing include methods for statistical difference testing between groups (both parametric and non-parametric). Similar to differential abundance testing, models for clinical studies that take into account the effect of an intervention on both the response (immune response) and biomarkers can be of interest as well. The primary difference however is that their focus is merely on association rather than causal relationships. To the best of our knowledge, there are only two tools that test for association between biomarkers (microbiome taxa compositions) and clinical response variables: NetShift using an algorithmic visualization (Kuntal et al., 2019a), and PhyloSeq using supervised methods (i.e., canonical correspondence analysis, discriminant correspondence analysis, sparse linear discriminant analysis, etc.) (McMurdie and Holmes, 2013). The authors of IVikodak listed the quantification of association between specific sets of bacteria with disease state as a planned future enhancement (Nagpal et al., 2019). None of them however allow for longitudinal analysis, taking into account the effect of an intervention on both the biomarkers and disease response.

Visualization—A wide variety of visual encodings have been used to represent the result of biomarker discovery analysis; ranging from simple heatmaps and bar charts, to more complex visuals like the volcano plot (Shamsaddini et al., 2020) and heat trees (Foster et al., 2017). An ongoing search noted by one of the experts in the focus group discussions is on how to visually represent the results of clinical longitudinal intervention studies: how do microbial composition and clinical response variables change over time given a particular intervention. In Figures 5B,D,F, some of the visualizations used in the visual analysis tools are shown.

Classification is used to classify samples in predefined groups based on their microbial composition. It provides information on the most important features (taxa) within sample cohorts, and is therefore often returning as a method for biomarker identification as well.

Analysis—Classification methods are fairly new in microbiome research, as only the more recently developed visual analysis tools cover these methods (Chong et al., 2020; Shamsaddini et al., 2020). Machine learning algorithms such as random forest classifiers or support vector machines are typically used for this type of analysis.

Visualization—Line charts (expressed as ROC curves) are typically used to represent model performance, whereas bar charts are used to display the most important features. An example of how this is been shown in literature is given in Figure 6C.

The analysis of microbial interaction is focused on identifying the relationship between species. Different types of relations can exist between microbes: mutualistic, commensal, parasitic and competitive (Faust et al., 2012). The goal is to find a method that reveals all of them at once. Identifying these relationships is important for all research objectives listed above. It provides more context on why certain taxa abundances differ in certain situations, and guides us towards possible causal relationships (e.g., is the drug altering the relative OTU abundance or is it altering its relative abundance through another taxa that contains a specific relationship with the OTU of interest).

Analysis and Visualization—Looking at the development of microbial interaction analysis within the microbiome visual analysis tools, new methods have been introduced during recent years, which gives an indication that the use of different methods is still further explored. At the moment, three schools of thought can be distinguished: 1) correlation based methods. Problem however with correlation is that it doesn’t correct for the compositionality of the data, and thus leads to spurious correlations (Gloor et al., 2017). Therefore, methods like SparCC, SPIEC-EASI and FastSpar were developed which result in network visualizations based on cut-off values (Chong et al., 2020). 2) Predator-Prey based methods using (generalized) Lotka Volterra equations to model relationships (Shaw et al., 2016; Kuntal et al., 2019b). 3) Topology based methods using topological data analysis (TDA) to construct the networks (Liao et al., 2019). All of these methods result in a graph visualized as a node-link diagram. Figures 6A,B provides an overview of how networks are used to represent microbial interactions.

As mentioned above on (relative) abundance, one could also look into the metabolic functions of microbial populations.

Analysis—Depending on the type of sequencing, different programs and methods can be used for functional profiling. Galloway-Peña and Hanson (2020) provide a nice overview including use cases and shortcomings. Using 16S rRNA sequencing, methods such as PICRUSt (Langille et al., 2013) and Tax4Fun (Aßhauer et al., 2015) allow to predict the gene content potential functionality based on a comparison between relative abundances and the reference genome of the taxa present. An important note of the authors that came up in the expert panel discussions as well is that these however are rough approximations, as they don’t take into account actual protein expressions. Using shotgun and metatranscriptome sequencing approaches, tools such as MetaGeneMark (Zhu et al., 2010) and Glimmer-MG (Kelley et al., 2012) carry out protein sequence homology based searches against databases of orthologues, enzymes, or protein domains and families for gene identification and annotation. The results could then be used for pathway enrichment analysis.

Visualization—The link between taxa and functions can be visualized using bipartite graphs (Figure 7C) or interactive stacked bar charts using highlighting, as was done in Burrito (McNally et al., 2018). The result of functional profiling are typically represented in a metabolic pathway network (Figure 7D) (Zhang et al., 2019).

As mentioned before in the section on differential abundance and repeated in the section on biomarker discovery, to gain a deeper understanding of causal relationships between the microbiome and various sample cohorts (e.g., grouped by disease state), longitudinal studies are required (Secrier and Schneider, 2013). Given the literature reviewed in this study, two tools were found to allow for longitudinal microbiome time series analysis; TIME (Baksi et al., 2018), and q2-longitudinal (Bokulich et al., 2018), which is an extension on QIIME2.

Analysis—In q2-longitudinal, linear mixed effect models are used to test for differential abundance. Changes of microbial sample compositions are captured across time using unweighted UniFrac, whereas in TIME dynamic time warping distance is used to capture groups of taxa showing similar trends over time. TIME identifies causal relationships among taxa using Granger Lasso causality. Stationary taxonomic groups (meaning no inter-microbial competition) are identified using an augmented dickey fuller test.

Visualization—Both tools allow for exploration of feature volatility using volatility plots (line charts) (Figure 7A). causal relationships between taxa are displayed using node-link diagrams (Figure 7B); clustering of taxa showing similar trends over time is visualized using a radial tree structure (Figure 7E).

Still, to the best of our knowledge no methods for longitudinal mediation analysis allowing for the identification of causal relationships between intervention, microbiome and response are incorporated yet.