Statistical Methods for Microbiome Compositional Data Network Inference: A Survey

Abstract

Microbes can be found almost everywhere in the world. They are not isolated, but rather interact with each other and establish connections with their living environments. Studying these interactions is essential to an understanding of the organization and complex interplay of microbial communities, as well as the structure and dynamics of various ecosystems. A widely used approach toward this objective involves the inference of microbiome interaction networks. However, owing to the compositional, high-dimensional, sparse, and heterogeneous nature of observed microbial data, applying network inference methods to estimate their associations is challenging. In addition, external environmental interference and biological concerns also make it more difficult to deal with the network inference. In this article, we provide a comprehensive review of emerging microbiome interaction network inference methods. According to various research targets, estimated networks are divided into four main categories: correlation networks, conditional correlation networks, mixture networks, and differential networks. Their assumptions, high-level ideas, advantages, as well as limitations, are presented in this review. Since real microbial interactions can be complex and dynamic, no unifying method has, to date, captured all the aspects of interest. In addition, we discuss the challenges now confronting current microbial interaction study and future prospects. Finally, we point out several feasible directions of microbial network inference analysis and highlight that future research requires the joint promotion of statistical computation methods and experimental techniques.

1. Introduction

Microbes exist widely in the biosphere, from fertile soil to warm spring water and from hard rocks to thin atmosphere, and even in extreme environments, such as high-temperature deep-sea liquid vents, rocky crusts in the depths of the earth, and cold poles (Huber et al., 2007; Pikuta et al., 2007). Microbes play a vital role in the biochemical cycle of the earth. In the ecological cycle of the biosphere, microbes act not only as producers, but also as the main decomposers of organic matter, demonstrating their importance as a force driving the biochemical cycle (Falkowski et al., 2008; Madsen, 2011).

Massine numbers of microbes live in the human body, and their total number is estimated to be 10 times that of human cells (Micah et al., 2007). They are distributed in various parts of the human body, such as skin, mouth, and intestines. Among them, the intestine has the largest number and variety of microbes (Filyk and Osborne, 2016). Microbes often exist in the form of communities, and they affect their own environment through synergy, which, in turn, affects the composition and structure of the microbial community.

The genomes of these microbes constantly interact with the human genome to regulate the metabolism of the host, which is closely related to human health. In recent years, many studies have emerged to explore the relationship between microbes and human health. For example, intestinal microbes are closely related to obesity (Turnbaugh et al., 2009), diabetes (Qin et al., 2012), high blood pressure (Holmes et al., 2008), cardiovascular diseases (Wang et al., 2011), AIDS (Lozupone et al., 2013), and cancer (Schwabe and Jobin, 2013). Meanwhile, people's eating habits and lifestyles will also affect the structure of intestinal flora (Wu et al., 2011). The Human Microbiome Project initiated by the American National Institutes of Health uses metagenomics to study the relationship among human body surface, microbial flora, and human health (Turnbaugh et al., 2007; Huttenhower et al., 2012; Methé et al., 2012).

The microbiome refers to the collection of all microbial species and their genetic information and functions that exist in a specific environment. It includes not only interactions among microbes in the environment, but also the interactions between microbes and other species and even environments. In recent years, with the rapid development of sequencing technology, researchers can extract genomic information of all microbes from microbial communities in the natural environment. More and more novel bacterial species have been discovered, and microbial sequencing data have grown exponentially (Proal et al., 2017).

Such sequencing makes it easier to study the correlations among microbial communities and between microbes and hosts. For example, 16S rRNA sequencing is a widely used high-throughput sequencing technology for microbiome studies. The 16S rRNA gene is the most common genetic marker to identify microbes (Wooley et al., 2010), and its advantages are mainly reflected in three aspects. First, the 16S rRNA gene exists in all microbes and has the same function in most organisms, except for eukaryotic cells.

Second, the length of 16S rRNA gene can reach 1500 bp, which can provide enough variation to identify the types of microbes (Janda and Abbott, 2007). Third, 16S rRNA sequences have highly conserved regions among different species, which allow us to use universal primers to sequence multiple species at the same time and distinguish the type of microbes based on the given hypervariable regions (Chakravorty et al., 2007).

Compared with other sequencing technologies, 16S rRNA sequencing technology also provides fast sequencing speed and low sequencing cost, which are also important factors in its widespread application in microbiome research. In the specific sequencing process, the variable regions of the 16S rRNA gene sequence are first amplified and sequenced, and then the sequences with a similarity of 95%, 97%, or 99% are grouped into one operational taxonomic unit (OTU); the OTU count is then used as a proxy variable for potential microbial abundance (Schloss, 2010; Morgan and Huttenhower, 2012). Further, grouping 16S rRNA gene sequences to obtain OTUs can be achieved in three ways.

The first is to cluster similar sequences by using an unsupervised clustering algorithm (Schloss et al., 2009), and the second is to use a phylogenetic tree that integrates mutation and evolution information (Hamady et al., 2010). The third method is to use the supervised learning model obtained from the labeled training data to determine the OTU type corresponding to the sequences (Wang et al., 2007). This method is also applied to whole-genome shotgun sequencing. To control the experimental deviation caused, for example, by sequencing depth or biological sample size, the total count in the sample is typically used to normalize the count of each OTU (Ozsolak and Milos, 2011; Ni et al., 2013; Fernandes et al., 2014).

These normalized data only reflect the relative abundance information of each OTU. In this way, microbiome sequencing data typically consist of compositional data (Kurtz et al., 2015).

A wide range of interactions takes place among microbes coexisting in nature. These interactions are vital to the earth's biochemical cycle and have a profound impact on human health and disease. Based on their influence on members of the microbial community, the various interactions among microbes can be divided into three types: beneficial, neutral, and harmful (Lidicker, 1979; Faust and Raes, 2012; Zengler and Zaramela, 2018). Beneficial interactions include commensalism and mutualism. The former means that the interaction is beneficial to one microbe but has no effect on other microbes.

The latter indicates that microbes coexist and benefit from each other. Neutral interaction represents the absence of any effect on both microbes. Harmful interactions are mainly divided into four situations: amenalism, competition, parasitism, and predation. Amenalism means that the connection is harmful to one microbe but has no effect on other microbes. Competition implies that all microbes involved in the interaction are at a disadvantage owing to the existence of other microbes. Parasitism occurs when microbes live off of their host, benefiting themselves at the expense of their hosts.

Predation means that one microbe feeds on another microbe. Together, these interactions constitute a complex ecological network that determines the structure, composition, and function of the microbial community in the ecosystem.

In microbiology and systems biology research, probabilistic graph models are typically used to characterize the interactions between various microbes in a microbial community (Jordan, 2003; Friedman, 2004; Whittaker, 2009; Cao et al., 2019). We can regard each kind of microbe as a vertex. If an interaction occurs between two microbes, we can connect their corresponding vertices to form an edge. In this way, we can get a network model describing the interactive relationship between microbes, namely the microbial interaction network.

Further, we can determine the strength of the interaction between microbes by assigning weights to edges such that the greater the absolute value of the weight, the stronger the interaction. A positive weight indicates beneficial interaction, and a negative weight indicates harmful interaction. Zero weight means neutral interaction. In the probabilistic graphical model, the absolute abundance of each microbe is regarded as a random variable, making the absolute abundance of all microbes in the microbial community a random vector. Microbiome studies are concerned with how the interaction network corresponds to each component in this random vector.

The goal of microbiome compositional data network analysis is to make statistical inferences about the interaction network corresponding to their unobserved absolute abundance by utilizing the observed relative abundance data of various microbes. Such interactions are often inferred by estimating the covariance matrix. In detail, interactions can be divided into those that are direct and those that are indirect. Indirect interaction may be mediated between two microbes, or it may be transmitted through a third microbe.

Direct interactions require no mediator or third factor and are more intrinsic. Since the two correlative estimation methods noted earlier do not make such distinctions, conditional correlation network inference methods that can estimate direct interactions have gained much attention. In addition, the interaction network between microbes is not constant for a long time, but changes as the external environment changes, such as the health status of the host, environmental conditions, or genetic conditions.

Hence, research on differential networks and mixture networks has recently become an important topic in microbiome studies. Mixture networks account for the heterogeneity of microbiome data and generalize the single underlying biological network into multiple networks. Differential networks focus on the dynamics of interaction networks and explore interactions from the aspect of comparison.

Previously, several researchers have reviewed the inference methods of microbiome networks, but they mainly focused on correlation networks (Cardona et al., 2016; Weiss et al., 2016; Dohlman and Shen, 2019). In contrast, in this article, we review some of the statistical calculation methods developed in recent years based on all four main network situations: correlation networks, conditional correlation networks, mixture networks, and differential networks. We summarize the high-level idea of these methods and discuss their model assumptions and scope of applications. We also provide some comments and prospects for future development of algorithms. Our article is organized as follows.

In the second section, we illustrate some background knowledge. From the third to the sixth section, we introduce typical estimation methods of correlation networks, conditional correlation networks, mixture networks, and differential networks accordingly. After that, we discuss the existing challenges to microbiome interaction network inference and aspects that can be further studied in the future. The detailed algorithm procedure of each method is summarized in Supplementary Appendix SA1.

2. Background

Sequencing depth is defined as the total sequencing count summed over all taxa of a sample and can vary across samples. Such variation of sequencing depth may be mainly caused by experimental techniques. Considering this fact, the microbial count data only carry information of relative abundances instead of absolute abundances of species, that is, the compositionality of microbiome data. First denote the observed count matrix as and the unobserved absolute abundance matrix as . Let $x_{i j} = \frac{z_{i j}}{\sum_{j = 1}^{p} z_{i j}},$ (1)

then the relative abundance matrix is . The sample version of $x_{i j}$ can be written as $\frac{y_{i j}}{\sum_{j = 1}^{p} y_{i j}}$ . When the count value of a taxon increases in a sample, the relative abundances of all other taxa in this sample decrease, even if their absolute abundances do not change. This phenomenon is also regarded as compositional bias. Owing to the sum-constant constraint, traditional correlation statistics may give spurious results. Take Pearson correlation as an example:

Thus, the correlation is characterized by a negative trend, even though the true abundances may be independent. Such problems can be managed by ratio transformation such as the centered log-ratio (clr) transformations and the additive log-ratio (alr) transformations.

Moreover, the number of features in the microbiome data is usually far larger than the sample size. In other words, microbiome data are high dimensional. This high dimensionality can lead to a heavy computational burden. Algorithms need to be designed carefully for efficient computation. Since the sequencing technique keeps improving, the scalability onto high-dimensional data having a larger scale is necessary. However, it is interesting to notice that from another perspective, high dimensionality plays a positive role in the analysis of compositional data. Although the number of features increases, the influence of false negative correlation trends may alleviate, which is called as “the blessing of dimensionality.”

Another feature of microbiome data is excessive zeros. In fact, zeros in the count data can be generally divided into two types. One type results from taxa that are present, but in proportions too low to be detected, and such zeros are often called technical zeros. Another type caused by the true absence of certain taxa, leading to essential zeros, is also called structural zeros (Peschel et al., 2020). Some methods choose to model these zeros, but the performance depends on the rationality of distribution assumptions.

In addition, many methods add a pseudo-count to avoid invalid computation. However, the choice of pseudo-count is subjective and may affect the performance of the method. Moreover, some methods choose to filter some nonsense taxa to preserve more informative features to handle problems caused by excessive zeros.

The composition of microbial communities can drastically vary from hosts to environments (Jiang et al., 2019). Thus, when collecting the data of many samples, a taxon that is relatively abundant in one sample may not be observed in other samples. This is the sparsity mentioned earlier. Another phenomenon is the tremendous variation among the dominant taxa collected in just one sample. This is termed as taxonomic heterogeneity (Jiang et al., 2019). In addition to taxonomic heterogeneity, multiple underlying microbial networks may determine the sample-taxa matrix in reality, rather than one single network. For example, in the human gut, the interaction between two taxa may depend on the enterotypic context of the individual's gut microbiome that can vary from individual to individual (Jiang et al., 2019). This is another type of heterogeneity that is functional.

3. Correlation Network

To identify the interactions among members of natural microbial communities, correlation analysis is proposed to measure the dependencies among microbes. However, the observed data only provide relative fractions of species. For compositional data, standard correlation analytic methods could result in spurious correlations. Therefore, various specialized methods based on compositional data have been developed to analyze the microbiome correlation network. See Table 1 for summary of correlation network methods.

Table 1.

Summary of Correlation Network Methods

Methods	Models	Advantages	Drawbacks	Availability
SparCC (2012)	Pearson correlations via Bayesian framework	1. Free of compositional bias;	1. Covariance matrix not positive semi-definite;	https://bitbucket.org/yonatanf/sparcc
SparCC (2012)		2. Performs well under low diversity and high sparsity.	2. High computational complexity.
CCLasso (2015)	Latent variable model	1. Free of compositional bias	1. Nonlinear interaction undetected;	https://github.com/huayingfang/CCLasso
		2. Guaranteeing the correlation matrix positive.	2. Lack of theoretical performance guarantee;
			3. Only study pairwise correlations.
REBACCA (2015)	Linear system using pairwise log ratios	1. Fitting with large sample size;	1. High computational requirement;	http://faculty.wcas.northwestern.edu/~hji403/REBACCA.htm
		2. Free of compositional bias.	2. Lack of theoretical performance guarantee;
			3. Strongly affected by noise.
COAT (2019)	Thresholding estimator for compositional data	1. Free of optimization and efficient;	1. Nonlinear interaction undetected;	https://github.com/yuanpeicao/COAT
		2. Theoretical guarantees provided;	2. No direct way to handle excessive zeros.
		3. Estimate large covariance on constrained data.

3.1. SparCC

3.1.1. High-level idea

SparCC (Friedman and Alm, 2012) is just about the earliest method that takes compositionality into consideration. It introduces correlation between log-transformed absolute abundances based on compositional data. Then based on assumptions that the number of different OTUs is large and the actual correlation network is sparse, SparCC estimates the linear Pearson correlations between the log-transformed components via approximation. It proposes both a basic procedure and an iterative version to infer such linear Pearson correlations. The iterative one identifies component pairs with high correlations at each step to guarantee sparsity until the identification is stable.

Comment. SparCC estimates the correlations of log-transformed absolute abundances without the involvement of relative information, but there are several limitations of this method. First, although absolute abundances cannot be observed, SparCC can be free of compositional bias using this concept, implying that it relies heavily on the quality of component counts data. Second, the inferred covariance matrix by SparCC cannot be guaranteed positive definite, and correlation coefficients may not locate in $[- 1, 1]$ . Lastly, SparCC may not be computationally efficient with high species dimension, and as pointed out (Friedman and Alm, 2012), it only estimates linear relationships between log-transformed abundances.

3.2. CCLasso

3.2.1. High-level idea

Similar to SparCC, CCLasso (Fang et al., 2015) considers the compositional property of the observed data in correlation analysis explicitly, and it has an extra advantage in that the estimated correlation matrix of the latent variables for compositional data is positive definite. CCLasso estimates correlations of log-transformed absolute abundance via weighted least-square regression with $ℓ_{1}$ penalty based on log-transformed relative abundance. To solve the optimization problem, it proposes an alternating direction algorithm from the augmented Lagrangian method.

Comment. However, the limitations of CCLasso are similar to those of SparCC. It also needs reliable component fraction estimation, only explains linear relationships, and does not offer theoretical performance guarantees. In addition, in terms of managing excessive zeros, CCLasso also adopts a simple pseudo count to avoid zero components, which may affect its performance. Lastly, the performance of CCLasso may not be satisfying if the number of taxa is much larger than that of samples.

3.3. REBACCA

3.3.1. High-level idea

REBACCA (Ban et al., 2015) identifies important co-occurrence patterns by finding sparse solutions. Different from SparCC, REBACCA constructs a linear system that is equivalent to log-ratio transformations and solves the system by using the $ℓ_{1}$ -norm shrinkage method under the sparsity assumption.

Comment. Similar to SparCC, REBACCA estimates pairwise relationships of log-transformed absolute abundances, so it is free of compositional bias. However, the iterative process may not scale to p, resulting in a high computational requirement. Besides, no theoretical performance guarantees are provided for REBACCA to work effectively. In addition, REBACCA may also fail when the sample size is small and it can be strongly affected by noise.

3.4. COAT

3.4.1. High-level idea

To solve the problems existing in the methods cited earlier and to provide statistical insights into the issue of identifiability and the impacts of dimensionality, COAT (Cao et al., 2019) was developed to estimate the sparse covariance matrix in a principled way and does provide theoretical guarantees on support recovery. It shows that the basis covariance is approximately identifiable if it belongs to the large sparse covariance family. Under the assumption that the basis covariance matrix is sparse, a composition-adjusted thresholding estimator is defined.

Comment. COAT is equivalent to thresholding the sample clr covariance matrix. So, it is free of optimization and efficient even when p is large. This work also fills the gap of large covariance estimation on constrained data by adopting covariance thresholding methods to compositional data. The work of COAT includes explicit statistical assumptions that are needed for it to be efficient, and it guarantees its performance theoretically. However, the estimated correlation relationship by COAT is linear. And there is no direct way for COAT to handle excessive zeros.

4. Conditional Correlation Network

Taking the compositionality of microbiome data into consideration, a microbial interaction network describes the pairwise correlations between two taxa. To restate, such correlations can be divided into direct interactions and indirect interactions. Indirect interactions involve two taxa correlated by some confounding factor, such as a certain environment condition, or a third species. Direct interaction means the nonexistence of confounding factors, and one taxon influences another directly. Since direct interactions are usually more fundamental and intrinsic (Friedman, 2004), methods modeling the conditional dependencies among taxa given all other taxa have gained much attention.

To infer such conditional dependencies, many methods turn to estimating a probabilistic graphical model in which each node represents a taxon, and each edge is undirected, indicating the conditional dependency of two taxa (Chen et al., 2020). If data follow a distribution in the exponential family, a straightforward statistical interpretation of such models can be made (Wainwright and Jordan, 2008). Taking the multivariate Gaussian distribution as an example, the element of the inverse covariance matrix is zero if and only if its two corresponding variables are conditionally independent given all other variables.

Therefore, to infer the conditional correlation network where an edge between two nodes represents conditional dependency given all other nodes, many researchers turn to estimating the nonzero entries of the inverse covariance matrix. The inverse covariance matrix is also termed as precision matrix. For real microbiome data, which are high-dimensional and sparse, the calculated covariance matrix is usually degenerate. In addition to this, we can only utilize the observed compositional data, which also makes estimating the precision matrix under absolute abundances challenging. See Table 2 for summary of conditional correlation network methods.

Table 2.

Summary of Conditional Correlation Network Methods

Methods	Approaches/models	Advantages	Drawbacks	Availability
SPIEC-EASI (2015)	Graphical models using selecting neighborhoods or inverse covariance	1. Compositional bias tackled;	1. Hard to recover networks having large hubs;	https://github.com/zdk123/SpiecEasi
		2. Prior knowledge incorporated.	2. Hard to handle excessive zeros.
gCoda (2017)	Logistic normal distribution	1. More accurate than SPIEC-EASI;	1. Nonconvex optimization function;	https://github.com/huayingfang/gCoda
		2. More computational, efficient, and stable than SPIEC-EASI.	2. Cannot guarantee the consistency of estimators.
MPLasso (2017)	Multivariate normal distribution incorporating prior information	1. Prior knowledge incorporated;	1. Untested on a dynamic model;	https://github.com/ChiehLo/MPLasso_RPackage
		2. The need to perform exhaustive searches eliminated.	2. Cannot directly extended to estimate associations among communities.
CD-Trace (2019)	A loss function for compositional data with D-trace loss	1. Convex loss making computation easier than gCoda;	1. Heavy computation complexity;	https://github.com/coamo2/CD-trace
		2. Unique solution under the exchangeable condition.	2. Exchangeable condition needs further checking.
CDTr (2019)	A loss function for compositional data with D-trace loss	1. Simpler form of loss function than that of CD-trace;	1. No theoretical guarantee of consistency;	https://github.com/coamo2/CDTr
		2. Unique solution under the exchangeable condition.	2. Exchangeable condition needs further checking.
MInt (2016)	Poisson-multivariate normal hierarchical model	1. Flexible mean-variance relationship;	1. Compositional bias unrectified;	https://cran.r-project.org/web/packages/MInt/
		2. Over-dispersion handled.	2. Indirect interactions may be regarded as direct interactions.
mLDM (2017)	Lognormal-Dirichlet Multinomial hierarchical model	1. Environmental factors incorporated;	1. Heavy computation burden;	https://github.com/tinglab/mLDM
		2. Compositional bias considered.	2. No theoretical guarantee.
HARMONIES (2020)	Zero-inflated negative binomial distribution with Dirichlet process prior	1. Able to capture zeros and over-dispersion.	1. Performance affected by small sample size.	https://github.com/shuangj00/HARMONIES
BC-GLASSO (2020)	Logistic normal multinomial distribution	1. Bias-corrected estimator proposed;	1. Potential interactions between the reference taxon and all other taxa not modeled;
		2. Precision matrix well defined.	2. Covariates not considered.
CompGlasso (2020)	Logistic normal multinomial distribution	1. Reference-invariant property.	1. Probably depending on the algorithm that is used.

4.1. SPIEC-EASI

4.1.1. High-level idea

SPIEC-EASI (Kurtz et al., 2015) estimates the covariance of the log-ratio transformed absolute abundance approximately by calculating that of clr transformed data. It focuses on direct interactions and estimates’ conditional independence by inferring an underlying graphical model. The overall inference procedure of SPIEC-EASI includes estimating graphical models, as well as selecting neighborhoods and sparse inverse covariance. It is worth noting that SPIEC-EASI scales as the dimension of data grows and can recover the graph even when $p ≫ n$ . Although not specifically introduced here, it also proposes a pipeline of generating synthetic data.

Comment. The inference engine of SPIEC-EASI has the advantage of incorporating prior knowledge about the underlying data or network structure from independent scientific studies in a principled manner. Nevertheless, this method has several weaknesses. First, SPIEC-EASI relies heavily on the conditional number of the inverse covariance matrix. It also ignores the uncertainty of estimates, which assumes that the proportion of a taxon in a sample approximates the true value. Second, it is difficult for clr transformation to handle excessive zeros, and the transformed data are not even close to obeying Gaussian or sub-Gaussian graphical models. Third, the consistency of estimated covariance matrix cannot be guaranteed.

4.2. gCoda

4.2.1. High-level idea

gCoda (Fang et al., 2017) assumes that the absolute abundances follow a multivariate log-normal distribution and that the underlying network is sparse. Therefore, it can model the generation of relative abundance by a logistic normal distribution, and then estimate conditional dependence structures, which indicate direct interactions in microbial communities by maximum likelihood, using $ℓ_{1}$ -penalty to ensure sparsity. Compared with SPIEC-EASI, gCoda is more computationally efficient and stable.

Comment. gCoda is just about the first method to derive the likelihood of compositional data for Gaussian graphical models, However, one shortcoming of gCoda is that the likelihood function is nonconvex. Besides, the consistency of the estimator cannot be guaranteed.

4.3. MPLasso

4.3.1. High-level idea

MPLasso (Lo and Marculescu, 2017) uses biological knowledge mined from the literature as prior information in the graphical Lasso. The prior information including microbial co-occurrence patterns and associations is searched directly from a large amount of scientific literature by using automated text mining. It is utilized to restrict the search space to get more possible existing associations and make interaction inference more reliable.

Comment. MPLasso proposes to use experimentally verified biological knowledge as prior knowledge via text mining to enhance accuracy in inferring associations. Prior information plays a role in penalizing the elements of precision matrix, which may not dominate the results but restrict the search space. If the prior information mined is rare, then the method may be considered a graphical lasso problem. However, it has not been tested on a dynamic model of microbial communities. Extending MPLasso to estimate associations among communities may, therefore, require a new algorithm.

4.4. CD-trace

4.4.1. High-level idea

Considering that graphical lasso is not computationally efficient, D-trace loss (Zhang and Zou, 2014) is, hence, proposed. It is a convex and smooth loss function with a unique minimizer at the inverse covariance matrix. The sparse precision matrix can then be estimated by minimizing the lasso penalized D-trace loss under a positive-definite constraint. D-trace loss can be viewed as an analogue of the least-squares loss used in regression to estimate the precision matrix. Based on D-trace loss, CD-trace (Yuan et al., 2019) introduces an empirical loss function for compositional data called compositional D-trace loss. Under the positive-definite constraint, minimizing lasso-penalized CD-trace loss can get a sparse matrix estimator for the direct microbial interaction network.

Comment. The loss function proposed in CD-trace is convex, which makes the solution easier than that in gCoda. However, it has several limitations. First, during application, CD-trace adds a pseudo-count to avoid zero counts, which may result in inflated composition and the Gaussian distribution assumption may be oversimplified. This method still needs further study in contending with excessive zeros and information from these zeros. Second, theoretical support for the consistency of the estimator is absent. Third, the computational complexity of CD-trace is $O (p^{3})$ , the same as that for gCoda and SPIEC-EASI, which may be heavy when p is too large.

4.5. CDTr

4.5.1. High-level idea

CDTr (He and Deng, 2019) is another loss function based on D-trace loss. Similar to CD-trace, CDTr also assumes that the log-transformed absolute abundance follows multivariate normal distribution. The sparse matrix estimators are defined as the minimizer of the corresponding lasso penalized loss, which has a unique solution under the exchangeable condition. For numerical solution, CDTr is based on the Alternating Direction Method of Multipliers (ADMM) algorithm to obtain a straightforward application.

Comment. We should consider the following limitations of CDTr. First, CDTr has a more concise form than CD-trace loss, but it still cannot guarantee the consistency of results. Second, the exchangeable condition used in CD-trace and CDTr needs theoretical verification and its biological reasonableness needs further checking.

4.6. MInt

4.6.1. High-level idea

MInt (Biswas et al., 2016) proposes a Poisson-multivariate normal hierarchical model to capture the conditional correlation structure. Confounding factors are controlled at the Poisson layer, and the multivariate normal layer with an $ℓ_{1}$ penalized precision matrix depicts taxon–taxon interactions. The optimization is obtained by an iterative conditional mode algorithm, which means we update parameters sequentially given the current estimates of other parameters.

Comment. MInt is based on a Poisson-multivariate normal model, which models the Poisson mean as a multivariate normal random variable, leading to the flexible mean-variance structure. It can also readily handle over-dispersion. When inferring direct interactions, it takes controlling for covariates into consideration. However, these potential confounding variables must be diligently measured; otherwise, indirect interactions may be regarded as direct interactions. Besides, MInt directly handles count data without considering the compositionality of microbiome data explicitly.

4.7. mLDM

4.7.1. High-level idea

mLDM (Yang et al., 2017) is a metagenomic Lognormal-Dirichlet-Multinomial hierarchical model. It first models the count matrix with a multinomial distribution and then uses a Dirichlet distribution to model the multinomial probabilities. Lastly, a multivariate log-normal distribution is utilized to fit the absolute abundances of samples. mLDM estimates both the conditional associations among microbes and the associations between microbes and environmental factors and considers compositional bias and variance of the metagenomic data. The sparse association network is estimated by sparse learning techniques.

Comment. We should take the following limitations of mLDM into consideration. First, mLDM assumes that microbes have linear relationships with environmental factors, which may not actually be the case. Second, many parameters are introduced in mLDM, which can limit its efficiency and scalability. The complexity of mLDM also reduces its interpretability. Third, mLDM lacks theoretical justification of the estimators.

4.8. HARMONIES

4.8.1. High-level idea

Considering that microbiome sequencing data usually have different sequence depths across samples, many zeros, and over-dispersion, normalization is an important step before analysis. Unlike traditional normalization, which gets the compositional data first and then applies log-ratio transformation (Kurtz et al., 2015), HARMONIES (Jiang et al., 2020b) is a model-based normalization method designed for microbiome count data. It uses the zero-inflated negative binomial (ZINB) model with the Dirichlet process as a prior (DPP) to normalize data. After normalization, HARMONIES estimates the partial correlations to infer direct associations among taxa.

Comment. HARMONIES takes into account excessive zeros and over-dispersion by modeling the count data with a ZINB distribution. It also takes sample heterogeneity into account by further implementing DPPs. However, in some cases, the sample size obtained is always small, and this limitation has a direct effect on the estimated normalized matrix from the ZINB model.

4.9. BC-GLASSO

4.9.1. High-level idea

To reduce bias in the estimator coming from the proportion-based approach, BC-GLASSO (Jiang et al., 2020) directly models compositional count data by a multinomial distribution. Then, it assumes that multinomial probabilities follow a logistic normal distribution. The procedure of BC-GLASSO is similar to that of SPIEC-EASI, but it defines the true covariance matrix based on multinomial probabilities after the alr transformation rather than abundances after clr transformation. Afterward, BC-GLASSO conducts a bias-correction procedure for the estimator of true covariance matrix, during which the unevenness of sequencing depths can be handled. Based on this bias-corrected estimator, the graphical lasso is applied to infer a sparse inverse covariance matrix.

Comment. The alr transformation used in BC-GLASSO has the advantage of leading to a positive definite covariance matrix and thus makes the precision matrix well-defined. BC-GLASSO is, in fact, a hierarchical model. The multinomial distribution can capture technical variation of the count data, and the logistic normal distribution can detect biological variation of microbiome composition. BC-GLASSO considers the bias induced by the estimator and corrects such bias. However, we should notice that it also has several drawbacks, First, BC-GLASSO does not now take covariates, such as sex and age, into consideration. Second, since alr transformation needs a reference taxon, the potential interactions between the reference taxon and other taxa may not be estimated and analyzed by BC-GLASSO.

4.10. CompGlasso and its reference invariant version

4.10.1. High-level idea

CompGlasso (Tian et al., 2020a) also adopts the logistic normal multinomial distribution to model the count data, which is two-level hierarchical. However, the objective function used in it is derived from log-likelihood, which is different from that of BC-GLASSO. Alr transformation involves choosing a reference taxon, which may affect network estimation. Therefore, further analyses are needed to determine whether using a different reference taxon could change the nature of the estimated network in compositional graphical lasso and further affect its robustness. The reference-invariance property established (Tian et al., 2020) shows that the specific submatrix of the precision matrix is invariant to the choice of reference taxon. According to this property, the reference-invariant version of compositional graphical lasso is proposed, which only penalizes the invariant submatrix (Tian et al., 2020).

Comment. CompGlasso establishes the theoretical convergence of the whole algorithm, but the theoretical property of the estimator remains to be elucidated. The reference-invariant property clearly shows how different reference taxon choices affect the estimated network. Based on this property, a reference-invariant optimization problem is designed to make the alr transformation free from the choice of the reference taxon. However, the empirical reference-invariant inverse covariance matrix and network may depend on the algorithm used to optimize the objective function.

5. Mixture Network

The methods described earlier all assume that the sample-taxa matrix is only related to a single network topology and set of edge weights. A sample-taxa matrix, however, is constructed by various samples in which taxa may be involved in more than one biological network. Moreover, different environmental conditions can lead to different association networks. In particular, for human gut microbiome data, other variables, such as age, sex, and diet, can also influence microbial network interactions, changing them according to these variables. See Table 3 for summary of mixture network methods.

Table 3.

Summary of Mixture Network Methods

Methods	Approaches/models	Advantages	Drawbacks	Availability
MixMPLN (2019)	Mixture of multivariate Poisson lognormal distribution	1. Different microbial clusters considered;	1. Heavy computation complexity;	https://github.com/sahatava/MixMPLN
		2. Compositionality is handled.	2. Biological significance of identified interactions needs to be further investigated.
MixMCMC (2019)	Mixture of multivariate Poisson lognormal distribution with MCMC	1. Different microbial clusters considered;	1. Significantly slower than MixMPLN;	https://github.com/sahatava/MixMCMC
		2. Able to recover the topologies for graphs more accurately.	2. Biological significance of identified interactions needs to be further investigated.
MixGGM (2019)	Mixture of multivariate Gaussian distribution	1. Different microbial clusters considered;	1. Depending strongly on the conditional number;	https://github.com/sahatava/MixGGM
		2. Compositionality considered.	2. clr transformed data may not fit distribution well;
			3. Heavy computation complexity.
kLDM (2020)	k-Lognormal Dirichlet-Multinomial hierarchical model	1. Compositional bias considered;	1. Heavy computation burden;	https://github.com/tinglab/kLDM
		2. Large variance of count data considered.	2. Gaussian distribution may not hold.

clr, centered log-ratio; MCMC, Markov Chain Monte Carlo; MPLN, Multivariate Poisson Log-Normal.

5.1. MixMPLN and MixMCMC

5.1.2. High-level idea

MixMPLN is a mixture model that comprises K Multivariate Poisson Log-Normal (MPLN) distributions (Tavakoli and Yooseph, 2019a). The mixture model directly tackles the count data, and the precision matrix of each component represents a microbial network. The MixMPLN framework is also extended by adding $ℓ_{1}$ - penalty for sparsity onto K precision matrices and during iteration, graphical lasso (Friedman et al., 2008) is applied to each component of the mixture model to estimate the precision matrix. Further, another algorithm, MixMCMC (Tavakoli and Yooseph, 2019), also models the sample-taxa count data by a mixture of MPLN distributions. It differs from MixMPLN in estimating the latent variables. MixMPLN estimates all parameters via the gradient descent method, whereas MixMCMC utilizes gradient descent along with Markov Chain Monte Carlo (MCMC) sampling, to estimate the latent variables.

Comment. MixMPLN tackles the count data directly, which considers different microbial clusters. Technical covariates, or environmental factors, which may influence a microbe network, can be included in the MixMPLN structure to control for confounding factors. However, there are several drawbacks of these two methods. First, as a mixture model, MixMPLN contains many parameters that can increase the time of optimization. Second, although MixMPLN could identify both positive and negative interactions between pairs of taxa for selected threshold values, the biological significance of these interactions should be investigated further. Lastly, the MCMC process in MixMCMC may be time-consuming.

5.2. MixGGM

5.2.1. High-level idea

MixGGM applies a log-ratio transformation to the sample count matrix and models the data as a mixture of Multivariate Gaussian distributions (Tavakoli and Yooseph, 2019). The iterative process of estimating parameters is similar to that in MixMPLN. The authors also extend MixGGM by adding sparse assumption onto the precision matrix for each multivariate Gaussian distribution. This is realized by applying graphical lasso to all K components.

Comment. MixGGM assumes that the covariance matrix of clr transformed data approximates that of log-transformed absolute abundances. Therefore, it may share the same problem as SPIEC-EASI, since approximation strongly depends on the conditional number of the inverse covariance matrix. Besides, the clr transformed data may not fit sub-Gaussian or Gaussian distribution well. Lastly, the computation complexity of MixGGM has not been well analyzed.

5.3. k-Lognormal-Dirichlet-multinomial model

5.3.1. High-level idea

As mentioned earlier, mLDM estimates both direct associations between two taxa and associations between taxa and environmental factors (Yang et al., 2017). However, this model only assumes one underlying biological network. Thus, the authors extended mLDM into k-Lognormal-Dirichlet-Multinomial model (kLDM) (Yang et al., 2020), which is a hierarchical Bayesian model and estimates multiple association networks with respect to specific environmental conditions. kLDM assumes that environmental factors regulate associations among microbes. These environmental factors tend to be stable given identical environmental conditions, but they vary with environmental conditions.

It also approximates the association pattern with multiple clusters, where each cluster represents an environmental condition. kLDM further assumes that the association networks are sparse. Then, it proposes a split-and-merge algorithm to estimate the number of environmental conditions and association networks.

Comment. For modeling, kLDM also takes compositional bias and the large variance of count data into consideration. Similar to mLDM, kLDM estimates both associations among microbes and associations between microbes and environmental factors. It can determine the number of environmental conditions automatically. However, kLDM has several limitations. First, the mathematical score used to merge clusters may not be sensitive enough to separate two clusters that represent two environmental conditions. Second, kLDM assumes that environmental factors follow a Gaussian distribution, which may not hold for categorical data. Lastly, kLDM shares a problem in common with mLDM in that too many parameters may lead to inefficient computation. This, in turn, results in unsatisfactory scalability of kLDM.

6. Differential Network

Mixture network inference methods estimate multiple underlying microbial networks simultaneously. In addition to considering the general case that multiple networks are associated with the sample-taxa matrix, discovering how the interaction network varies with environmental conditions or other confounding factors is another study that is of much significance to researchers. The differential network is often defined as the difference between two precision matrices. See Table 4 for summary of differential network methods.

Table 4.

Summary of Differential Network Methods

Methods	Approaches/models	Advantages	Drawbacks	Availability
DCDTr (2017)	A loss function for differential network estimation based on D-trace loss	1. Estimating the precision matrix difference directly;	1. Reasonability of exchangeable condition remaining doubtful;	https://github.com/coamo2/DCDTr
		2. Convex loss-making computation easier.	2. No theoretical guarantee.
MDiNE (2020)	Hierarchical Bayesian model	1. Estimation of the differential network using a binary covariate considered;	1. Long running time;	https://github.com/kevinmcgregor/mdine
		2. Able to estimate the co-occurrence network of two groups simultaneously.	2. Unable to control data with excessive zeros;
			3. Unable to explore how networks change with a continuous covariate.
Differential Markov random field analysis (2019)	Markov random field model	1. Free of distribution assumptions;	1. Strongly affected by the chosen discretization threshold;	https://rdrr.io/github/drjingma/TestBMN/man/TestBMN-package.html
		2. Statistical significance can be assigned to each element.	2. Each element can only take two values;
			3. Consistency of results needs examination.

However, using the difference between these two estimators of a precision matrix as the estimated differential network can be problematic. One reason is that the estimator of a single precision matrix is usually obtained under the sparsity assumption. Directly minusing such two estimators may lead to false positive and negative interactions (Yuan et al., 2017). Even if the sparsity assumption were to be held, these methods would have most power in exploring strong interactions in a single condition. Thus, when detecting interactions that are not strong in one condition, but change strongly in other conditions, the power of these methods is limited (Yuan et al., 2017).

6.1. DCDTr

6.1.1. High-level idea

As mentioned earlier, CDTr incorporates the clr transformation into D-trace (Zhang and Zou, 2014) to derive the conditional correlation network. Further, to infer microbiome differential networks, while avoiding the use of unobserved absolute abundance data, DCDTr loss (He and Deng, 2019) is proposed by combining clr transformation with new D-trace (Yuan et al., 2017). DCDTr assumes that the difference between the precision matrices is sparse, rather than both matrices are sparse. The DCDTr loss proposed for differential network estimation is smooth and convex under the sparse differential network assumption. It can estimate the precision matrix difference directly with compositional data instead of estimating the precision matrices individually.

Comment. Similar to CDTr, the loss function proposed in DCDTr has a unique minimizer as the precision matrix if the exchangeable condition is satisfied. However, the reasonableness of the exchangeable condition still needs to be further studied. Most importantly, no theoretical guarantee on the consistency of the estimator is provided.

6.2. MDiNE

6.2.1. High-level idea

MDiNE (Mcgregor et al., 2020) is a Bayesian hierarchical model, which can estimate the co-occurrence network of two groups simultaneously with respect to a binary covariate and thus estimate the difference between them. MDiNE models the counts of individual taxa through a multinomial distribution, with its probabilities depending on a latent Gaussian random variable. In addition, the co-occurrence network among taxa is determined by a sparse precision matrix over all the latent items. The model is estimated through the Hamiltonian Monte Carlo algorithm.

Comment. MDiNE makes interval estimations for both parameters of the network and the difference between networks available. However, it has several limitations. First, an underestimation of true associations may occur when the sample size is relatively small. Second, MDiNE has no zero-inflation parameter, so it may encounter problems in controlling data with excessive zeros. Third, MCMC sampling may take too much time and increase computational burden.

6.3. Differential Markov random field analysis

6.3.1. High-level idea

Since previous graphical models cannot assign statistical significance for estimated associations, a Markov random field model and a hypothesis testing framework are proposed (Cai et al., 2019) to detect microbial network structure and differences between networks. They first discretize each composition vector into a binary one. Then, conditional dependencies between components in the vector are modeled by a binary Markov random field. To test whether two precision matrices are identical, they conduct a global test first. If the null hypothesis is rejected, then the multiple tests of each entry are performed. The binary Markov random field accounts for the excessive zeros and unit-norm constraint in the compositional data.

Comment. This model does not require the assumption of data distribution, which, in turn, increases its robustness. It can be used to interpret the co-occurrence and co-exclusiveness of the microbial communities naturally. And it can assign statistical significance to each element of the differential network. Nevertheless, there are several drawbacks that may limit its usage. First, this method is strongly affected by the discretization threshold chosen. Taking different cut-off values to transform raw data into binary data may cause different results. Second, each element of the discretized vector can only take two values. It should be extended to study the case in which each node can take more than two discrete values. Lastly, the consistency of results needs to be further examined, such as results when dichotomizing the data in different ways.

7. Discussion

7.1. Ways to handle environmental confounding factors

Environmental factors, such as pH, humidity, temperature, and soil type, can strongly affect the composition of the microbial community. Understanding interactions among microbes and between microbes and environmental factors is of great importance. However, whether the edge in the microbial network is caused by an environmental factor or a direct interaction between two taxa is hard to distinguish, and many detected associations between microbes may be spurious without considering environmental factors (Yang et al., 2017).

Microbial network inference can handle environmental factors by regarding them as additional species and calculating their association with microbes, or treating them as extra variables in regression and inferring the association of the remaining abundances unaffected by the environment (Faust and Raes, 2012), such as mLDM (Yang et al., 2017) and kLDM (Yang et al., 2020). In addition, as mentioned in Faust (2021), we can cluster samples according to key environmental factors and construct networks for each sample group separately to decrease environmental variation within a group.

Multiple ways are available to predict the associations between species and environmental factors, but which one to use relies on concrete goals during the research. For instance, treating environmental factors as additional nodes helps us realize how they influence microbial community composition, whereas stratifying samples or regressing out environmental confounders is beneficial for inferring microbial interactions. However, the problem lies in establishing a reliable evaluation system of various treatments to properly manage environmental factors (Faust, 2021).

7.2. Ways to detect complex relationships

As mentioned earlier, methods that estimate interaction by means of Pearson correlations (Friedman and Alm, 2012), covariance matrix, or precision matrix all assume that the microbial correlation relationship is pairwise linear, which is not always the case. Sometimes, nonlinear complex interactions may occur in the microbial community. An interaction among some species could be influenced by an additional one (Faust, 2021), and one species could be affected by other species (Faust and Raes, 2012).

These cases could not be handled correctly with direct pairwise interaction detection methods mentioned in this article. In fact, high-order interactions have been shown to affect the microbial community both in simulations (Bairey et al., 2016) and in experiments (Gould et al., 2018). Hence, it is necessary to find ways to manage complex relationships, which could provide evidence for nonrandom co-association patterns (Cardona et al., 2016).

We can detect complex relationships by only considering simple relationships first. After the model is constructed, the biases between model estimations and actual microbial community behaviors may indicate the presence of complex relationships (Faust, 2021). However, improper model construction also can lead to the skewing of model predictions (Billick and Case, 1994), which makes it hard to distinguish the existence of complex relationships.

We can also handle complex relationship inferences by association rule mining, which involves enumerating all logical rules in the presence or absence of data to find valuable rules (Faust and Raes, 2012). For instance, if species A exists only in the presence of both species B and C, then neither the interaction between A and B nor the interaction between A and C exists unless a third species appears. This can be regarded as an extreme example of relationship correction (Faust, 2021).

People have reported several association rules containing more than two taxa in the past. However, possible reasons behind relationships detected may include true complex relationships and confounding results such as the assembly of pairwise interactions and overfitting (Faust, 2021). Thus, it is still a problem to uncover the underlying reasons that can lead to the detection of complex relationships, while still avoiding meaningless results.

Another approach worth trying is multiple regression, which predicts the abundance of a species through incorporating the abundance of other species (Faust and Raes, 2012). The principle of multiple regression is simple and clear and is, therefore, widely applied. However, this method is limited by the difficulty of discovering the implication of results when we sort factors through automatic feature selection methods such as sparse restriction. In this way, although we may be able to provide a roughly accurate prediction of the abundance of the species mathematically (Faust and Raes, 2012), the existence of causality in biology is doubtful.

In addition, people also attempt to capture complex forms of ecological interactions through the mutual information criterion, which measures the additional information known about one variable when the other variable is given (Butte and Kohane, 2000). Maximum information coefficient (MIC) is widely used to capture novel correlations on large datasets (Cardona et al., 2016), which puts mutual information into a continuous distribution and identifies diverse interactions between variables, both linear and nonlinear (Reshef et al., 2011).

However, sample size must be large enough to calculate the value of MIC accurately, whereas this can increase the difficulty of implementation (Li et al., 2016). Besides, recent research has questioned the equitability and statistical performance of MIC (Gorfine et al., 2012).

Finally, a microbial network involves complex relationships that likely consist of edges jointing more than two nodes. Consequently, it is not an easy task to visualize them (Faust and Raes, 2012). Hypergraphs should be taken into consideration, and how to investigate and explain them needs further exploration (Faust, 2021).

7.3. Tradeoff between model accuracy and computational complexity

We realize the fact that simple models may fail to fit the microbial community well, whereas complex models may lead to computational burden. For instance, compared with Pearson correlation, SparCC is proven to be more suitable for avoiding false interactions, while inducing higher computational complexity (Matchado et al., 2021). Further, correlation-based methods usually cannot distinguish direct from indirect interactions. To solve this problem, a large number of methods based on conditional dependence have been developed.

Compared with correlation-based methods, some conditional dependence methods have a sharp increase in computational complexity, which results in unsustainable running times (Matchado et al., 2021), such as hierarchical Bayesian models, for instance, mLDM (Yang et al., 2017) and kLDM (Yang et al., 2020). Therefore, the tradeoff between model accuracy and computational complexity is another point we should take into account.

To solve this dilemma, we should consider the computational burden brought by concrete datasets and available computation resources when choosing a proper analytic method to use. We should also consider whether it is feasible for the microbiome to be analyzed on a higher taxonomic level.

Since different methods can detect different edges in the same data (Weiss et al., 2016), it is assumed that model accuracy can be improved by using a combination of various analytical tools. As a matter of fact, some studies have reported experiments on linearly modeled data with foregone engineered interactions (Weiss et al., 2016). Compared with most individual methods, it has been shown that the accuracy of the ensemble method is significantly improved at the cost of sensitivity. Hence, when the given research requires us to obtain results of high precision, for instance, utilizing hypotheses of microbial interactions to identify co-culturing, integrated methods should be applied.

7.4. Lack of universal standardized simulation and evaluation

Simulated datasets are often needed to test the performance of methods. Various ways are currently available for simulation. One famous approach is generating data from statistical distributions, such as the logistic normal distribution adopted by CCLasso (Fang et al., 2015), COAT Cao et al. (2019), and gCoda (Fang et al., 2017). Another widely used method involves adopting specific models, especially differential equations, such as the famous generalized Lotka–Volterra model (Buffie et al., 2015; Bucci et al., 2016; Venturelli et al., 2018).

However, we are still confounded by the lack of a universal standard for simulating microbial community with ground-truth correlations. This makes it difficult to evaluate microbial network inference methods in a unified way since the performance of the proposed methods is largely affected by the data generation approach. If the artificial data are generated in a way that approximates assumptions of the microbial correlation network inference model, then the model would gain a relatively good performance. That is, one way of simulating data would favor methods similar to its generation approach (Faust, 2021). However, the true generation process of microbial data in reality is unknown, which may be more complex than current simulation approaches.

To avoid the evaluation bias caused by a single simulation method and judge the performance of various methods more objectively, it is necessary to adopt different simulation methods when testing. It is also important to separate developing new microbiome correlation network inference methods from comparison studies to ensure neutrality (Faust, 2021). Weiss et al. (2016) compare the performance of eight correlation techniques, as well as their ability of handling the issues in compositional data. Hirano and Takemoto (2019) investigate more recent methods, such as REBACCA (Ban et al., 2015) and CCLasso (Fang et al., 2015).

Comparative studies show that correlation estimation methods based on compositional data do not always outperform classical methods, such as Spearman correlation-based and MIC-based methods. This may result from the fact that newly proposed methods are developed under certain conditions or for particular problems.

7.5. Lack of benchmark microbial network dataset

In addition to simulated datasets, real datasets with verified biological microbial associations are also vital to evaluate those interaction inference methods. Biological gold standard network datasets for interaction inference are considered those datasets that have available sequencing data and validated known interactions. Lima-Mendez et al. (2015) built a literature-curated collection of 574 interactions in marine eukaryotic plankton, in which both mutualism and parasitism interactions are included. Besides, Durán et al. (2018) examined root-associated microbial communities and detected some microbial interactions that have been verified experimentally. However, benchmark datasets for microbiome interaction network inference are still scarce.

Several obstacles interfere with constructing benchmark real datasets. First, some interactions that happen in the real world may not be observed at all. Thus, the verified interactions are inadequate and incomplete. Using these datasets as benchmarks to test microbiome network inference methods would improve our evaluation of how these methods discover known associations. Even if a method could catch sight of underlying interactions among various taxa, such interactions may be viewed as false positives since they were not observed experimentally.

Second, as previously pointed out, the interactions can be either direct or indirect. In reality, they are complex, and it is difficult for some methods to tell which kind of interactions are inferred. Especially, correlation network inference methods can obtain direct interactions or indirect ones, such as high-order interactions where associations among some species are altered by an additional species. The gap between inferred interactions and expected ones brings a challenge to the development of benchmark microbial datasets.

8. Conclusion

Microbial interaction estimation provides an insight into the complex interplay that exists in microbial communities and the connections between the microbiome and its living environment. Statistical computational methods can infer such association networks and guide further biological studies.

In this review, we follow the development of microbiome network inference and analyze typical methods of correlation networks, conditional correlation networks, mixture networks, and differential networks. These methods tackle inference problems from various aspects, and they are based on different assumptions, even while facing some common challenges coming from the compositional data.

In practice, we recommend the following principle to choose a proper method in a certain scenario. First, taking factors such as type of interactions, heterogeneity of data, and the dynamics of networks into consideration, we should decide which network structure we are concerned about. Second, we should consider the characteristics of data, such as data distribution, the existence of noise, as well as the availability of prior information. Finally, sample size and computation complexity are also determining factors for users to make a decision.

After estimating correlations among taxa, their biological interpretation is further required, which is not an easy task. Most current methods fail to recognize the nature of relationships between microbes and cannot credibly distinguish authentic biological correlations from artifact ones. More specific explanations require the support of biological experiments and the literature.

In addition, we have mentioned that excessive zeros as well as the curse of dimensionality bring in challenges in microbiome compositional data network inference. To mitigate their negative effects, people often choose to group relative abundances on a higher taxonomic level. However, this may lead to hiding insights into smaller communities.

In addition to types of networks analyzed in the main paper, microbiome network inference methods can be improved by accounting for more factors affecting the living microbiome, estimating more complex relationships, and considering time-variable networks. A more advanced network theory is required to offer a further comprehensive understanding into the mechanism for the microbiome to influence the host and vice versa.

For instance, we should not ignore the fact that microbes from various life domains may have mutual effects and usually live in the same habitat. Compared with a single domain, joint cross-domain network inference can result in marked differences of the features and integral structure of the network. Thus, it is a vital developing direction to extend existing methods to cross-domain network analysis.

Nowadays, microbiome multi-omics datasets can be built by collecting multiple types of biological data, such as metagenomics, metatranscriptomics, 16S, mass spectrometry-based metaproteomics, and metabolomics (Helbling et al., 2012). Different types of biological data provide information from different views, so consolidating them together can bring a more comprehensive view. However, in reality, it is rare that all samples have all types of features measured, owing to the difficulty of gathering multi-omics data (Jiang et al., 2019). Therefore, integration methods should be able to manage data with missing features and not rely on paired multi-omics data.

How to make methods more scalable with less computational complexity, while, at the same time, guaranteeing efficiency, is also a task for future studies. As discussed earlier, some challenges in benchmarking real microbiome datasets need resolution to simulate reliable data for evaluation and propose universal valuation standards. All these problems require improvements in biological techniques and statistical methods together.

Footnotes

Acknowledgments

The authors sincerely thank all the reviewers and editors for their efforts and suggestions. They wish Professor Michael Waterman “A happy birthday!”

Author Disclosure Statement

The authors declare they have no conflicting financial interests.

Funding Information

This work was supported by the National Key Research and Development Program of China (2021YFF1200902) and the National Natural Science Foundation of China (12126305, 31871342).

Supplementary Material

References

Ban

, An

, and Jiang

2015. Investigating microbial co-occurrence patterns based on metagenomic compositional data. Bioinformatics, 31, 3322–3329.

Bairey

, Kelsic

, and Kishony

2016. High-order species interactions shape ecosystem diversity. Nat. Commun. 7, 1–7.

Billick

, and Case

T.J.

1994. Higher order interactions in ecological communities: What are they and how can they be detected?. Ecology, 75, 1529–1543.

Biswas

, McDonald

, Lundberg

D.S.

, et al. 2016. Learning microbial interaction networks from metagenomic count data. J. Comput. Biol. 23, 526–535.

Bucci

, Tzen

, Li

, et al. 2016. MDSINE: Microbial Dynamical Systems INference Engine for microbiome time-series analyses. Genome Biol. 17, 1–17.

Buffie

, Bucci

, Stein

, et al. 2015. Precision microbiome reconstitution restores bile acid mediated resistance to Clostridium difficile. Nature, 517, 205–208.

Butte

A.J.

, and Kohane

I.S.

2000. Mutual information relevance networks: Functional genomic clustering using pairwise entropy measurements. Pac. Symp. Biocomput. 5, 418–429.

Cai

T.T.

, Li

, Ma

, et al. 2019. Differential Markov random field analysis with an application to detecting differential microbial community networks. Biometrika, 106, 401–416.

Cao

, Lin

, and Li

2019. Large covariance estimation for compositional data via composition-adjusted thresholding. J. Am. Stat. Assoc. 114, 759–772.

10.

Cardona

, Weisenhorn

, Henry

, et al. 2016. Network-based metabolic analysis and microbial community modeling. Curr. Opin. Microbiol. 31, 124–131.

11.

Chakravorty

, Helb

, Burday

, et al. 2007. A detailed analysis of 16S ribosomal RNA gene segments for the diagnosis of pathogenic bacteria. J. Microbiol. Methods, 69, 330–339.

12.

Chen

, He

, Zhai

, et al. 2020. Direct interaction network inference for compositional data via codaloss. J. Bioinform. Comput. Biol. 18, 2050037.

13.

Dohlman

A.B.

, and Shen

2019. Mapping the microbial interactome: Statistical and experimental approaches for microbiome network inference. Exp. Biol. Med. 244, 445–458.

14.

Durán

, Thiergart

, Garrido-Oter

, et al. 2018. Microbial interkingdom interactions in roots promote Arabidopsis survival. Cell, 175, 973–983.

15.

Falkowski

P.G.

, Fenchel

, Delong

E.F.

2008. The microbial engines that drive Earth's biogeochemical cycles. Science, 320, 1034–1039.

16.

Fang

, Huang

, Zhao

, et al. 2015. CCLasso: Correlation inference for compositional data through lasso. Bioinformatics, 31, 3172–3180.

17.

Fang

, Huang

, Zhao

, et al. 2017. gCoda: Conditional dependence network inference for compositional data. J. Comput. Biol. 24, 699–708.

18.

Faust

2021. Open challenges for microbial network construction and analysis. ISME. J. 15, 3111–3118.

19.

Faust

, and Raes

2012. Microbial interactions: From networks to models. Nat. Rev. Microbiol. 10, 538–550.

20.

Fernandes

A.D.

, Reid

J.N.

, Macklaim

J.M.

, et al. 2014. Unifying the analysis of high-throughput sequencing datasets: Characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis. Microbiome, 2, 15.

21.

Filyk

H.A.

, and Osborne

L.C.

2016. The multibiome: The intestinal ecosystem's influence on immune homeostasis, health, and disease. EBioMedicine, 13, 46–54.

22.

Friedman

, and Alm

E.J.

2012. Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8, e1002687.

23.

Friedman

J.H.

, Hastie

, and Tibshirani

2008. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9, 432–441.

24.

Friedman

2004. Inferring cellular networks using probabilistic graphical models. Science, 303, 799–805.

25.

Gorfine

, Heller

, and Heller

2015. Comment on “Detecting Novel Associations in Large Data Sets.” Science 80, 1–6.

26.

Gould

A. L.

, Zhang

, Lamberti

, et al. 2018. Microbiome interactions shape host fitness. Proc. Natl. Acad. Sci. USA. 115, E11951–E11960.

27.

Hamady

, Lozupone

, and Knight

2010. Fast UniFrac: Facilitating high-throughput phylogenetic analyses of microbial communities including analysis of pyrosequencing and PhyloChip data. ISME. J. 4, 17–27.

28.

, and Deng

2019. Direct interaction network and differential network inference from compositional data via lasso penalized D-trace loss. PLoS One, 14, e0207731.

29.

Helbling

D.E.

, Ackermann

, Fenner

, et al. 2012. The activity level of a microbial community function can be predicted from its metatranscriptome. ISME. J. 6, 902–904.

30.

Hirano

, and Takemoto

2019. Difficulty in inferring microbial community structure based on co-occurrence network approaches. BMC Bioinf. 20, 1–14.

31.

Holmes

, Loo

R.L.

, Stamler

, et al. 2008. Human metabolic phenotype diversity and its association with diet and blood pressure. Nature, 453, 396–400.

32.

Huber

J.A.

, Welch

D.B.M.

, Morrison

H.G.

, et al. 2007. Microbial population structures in the deep marine biosphere. Science, 318, 97–100.

33.

Huttenhower

, Gevers

, Knight

, et al. 2012. Structure, function and diversity of the healthy human microbiome. Nature, 486, 207.

34.

Janda

J.M.

, and Abbott

S.L.

2007. 16S rRNA gene sequencing for bacterial identification in the diagnostic laboratory: Pluses, perils, and pitfalls. J. Clin. Microbiol. 45, 2761–2764.

35.

Jiang

, Armour

C.R.

, Hu

, et al. 2019. Microbiome multi-omics network analysis: Statistical considerations, limitations, and opportunities. Front. Genet. 10, 995.

36.

Jiang

, Sharpton

, and Jiang

2020a. Microbial interaction network estimation via bias-corrected graphical lasso. Stat. Biosci. 13, 329–350.

37.

Jiang

, Xiao

, Koh

A.Y.

, et al. 2020b. Harmonies: A hybrid approach for microbiome networks inference via exploiting sparsity. Front. Genet. 11, 445.

38.

Jordan

M.I.

2003. An introduction to probabilistic graphical models.

39.

Kurtz

Z.D.

, Müller

C.L.

, Miraldi

E.R.

, et al. 2015. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput. Biol. 11, e1004226.

40.

, Lim

, Chng

K.R.

, et al. 2016. Predicting microbial interactions through computational approaches. Methods, 102, 12–19.

41.

Lidicker Jr

W.Z.

1979. A clarification of interactions in ecological systems. Bioscience, 29, 475–477.

42.

Lima-Mendez

, Faust

, Henry

, et al. 2015. Determinants of community structure in the global plankton interactome. Science, 348, 1262073.

43.

, and Marculescu

2017. MPLasso: Inferring microbial association networks using prior microbial knowledge. PLoS Comput. Biol. 13, e1005915.

44.

Lozupone

C.A.

, Li

, Campbell

T.B.

, et al. 2013. Alterations in the gut microbiota associated with HIV-1 infection. Cell Host Microbe. 14, 329–339.

45.

Madsen

E.L.

2011. Microorganisms and their roles in fundamental biogeochemical cycles. Curr. Opin. Biotechnol. 22, 456–464.

46.

Matchado

M.S.

, Lauber

, Reitmeier

, et al. 2021. Network analysis methods for studying microbial communities: A mini review. Comput. Struct. Biotechnol. J. 19, 2687–2698.

47.

Mcgregor

, Labbe

, and Greenwood

C.M.

2020. MDiNE: A model to estimate differential co-occurrence networks in microbiome studies. Bioinformatics, 36, 1840–1847.

48.

Methé

B.A.

, Nelson

K.E.

, Pop

, et al. 2012. A framework for human microbiome research. Nature, 486, 215.

49.

Micah

, Claire

F.L.

, and Rob

2007. The human microbiome project: Exploring the microbial part of ourselves in a changing world. Nature, 449, 804–810.

50.

Morgan

X.C.

, and Huttenhower

2012. Human microbiome analysis. PLoS Comput. Biol. 8, e1002808.

51.

, Yan

, and Yu

2013. How much metagenomic sequencing is enough to achieve a given goal. Sci. Rep. 3, 1–7.

52.

Ozsolak

, and Milos

P.M.

2011. RNA sequencing: Advances, challenges and opportunities. Nat. Rev. Genet. 12, 87–98.

53.

Peschel

, Müller CL von Mutius

, et al. 2020. NetCoMi: Network construction and comparison for microbiome data in R. Brief. Bioinform. 22, bbaa290.

54.

Pikuta

E.V.

, Hoover

R.B.

, and Tang

2007. Microbial extremophiles at the limits of life. Crit. Rev. Microbiol. 33, 183–209.

55.

Proal

A.D.

, Lindseth

I.A.

, and Marshall

T.G.

2017. Microbe-microbe and host-microbe interactions drive microbiome dysbiosis and inflammatory processes. Discov. Med. 23, 51–60.

56.

Qin

, Li

, Cai

, et al. 2012. A metagenome-wide association study of gut microbiota in type 2 diabetes. Nature, 490, 55–60.

57.

Reshef

D.N.

, Reshef

Y.A.

, Finucane

H.K.

, et al. 2011. Detecting novel associations in large data sets. Science, 6062, 1518–1524.

58.

Schloss

P.D.

2010. The effects of alignment quality, distance calculation method, sequence filtering, and region on the analysis of 16S rRNA gene-based studies. PLoS Comput. Biol. 6, e1000844.

59.

Schloss

P.D.

, Westcott

S.L.

, Ryabin

, et al. 2009. Introducing mothur: Open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541.

60.

Schwabe

R.F.

, and Jobin

2013. The microbiome and cancer. Nat. Rev. Cancer, 13, 800–812.

61.

Tavakoli

, and Yooseph

2019a. Learning a mixture of microbial networks using minorization–maximization. Bioinformatics 35, i23–i30.

62.

Tavakoli

, and Yooseph

2019b. Algorithms for inferring multiple microbial networks, 223–227. In IEEE International Conference on Bioinformatics and Biomedicine (BIBM). San Diego, California, USA. DOI: 10.1109/BIBM47256.2019.8983194.

63.

Tian

, Jiang

, and Jiang

2020a. Microbial Network Recovery by Compositional Graphical Lasso. (In preparation). Available at: http://sites.science.oregonstate.edu/~jiangd/papers/CompGLasso.pdf Last viewed August 2021.

64.

Tian

, Jiang

, and Jiang

2020b. Reference-invariance property of inverse covariance matrix estimation under additive log-ratio transformation and its application to microbial network recovery. (In preparation). Available at: http://sites.science.oregonstate.edu/~jiangd/papers/InverseCovariance.pdf Last viewed August 2021.

65.

Turnbaugh

P.J.

, Hamady

, Yatsunenko

, et al. 2009. A core gut microbiome in obese and lean twins. Nature, 457, 480–484.

66.

Turnbaugh

P.J.

, Ley

R.E.

, Hamady

, et al. 2007. The human microbiome project. Nature, 449, 804–810.

67.

Venturelli

O.S.

, Carr

A.C.

, Garth

, et al. 2018. Deciphering microbial interactions in synthetic human gut microbiome communities. Mol. Syst. Biol. 14, e8157.

68.

Wainwright

M.J.

, and Jordan

M.I.

2008. Graphical models, exponential families, and variational inference. Found Trends Mach. Learn. 1, 1–305.

69.

Wang

, Garrity

G.M.

, Tiedje

J.M.

, et al. 2007. Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267.

70.

Wang

, Klipfell

, Bennett

B.J.

, et al. 2011. Gut flora metabolism of phosphatidylcholine promotes cardiovascular disease. Nature, 472, 57–63.

71.

Weiss

, Van Treuren

, Lozupone

, et al. 2016. Correlation detection strategies in microbial data sets vary widely in sensitivity and precision. ISME. J. 10, 1669–1681.

72.

Whittaker

2009. Graphical models in applied multivariate statistics. Wiley Publishing, Hoboken, New Jersey, USA.

73.

Wooley

J.C.

, Godzik

, and Friedberg

2010. A primer on metagenomics. PLoS Comput. Biol. 6, e1000667.

74.

G.D.

, Chen

, Hoffmann

, et al. 2011. Linking long-term dietary patterns with gut microbial enterotypes. Science, 334, 105–108.

75.

Yang

, Tan

, Han

, et al. 2020a. Correlation-Centric Network (CCN) representation for microbial co-occurrence patterns: New insights for microbial ecology. NAR Genomics Bioinf. 2, lqaa042.

76.

Yang

, Chen

, and Chen

2017. Inference of environmental factor-microbe and microbe-microbe associations from metagenomic data using a hierarchical Bayesian statistical model. Cell Syst. 4, 129–137.

77.

Yang

, Wang

, Xie

, et al. 2020b. Inferring multiple metagenomic association networks based on variation of environmental factors. bioRxiv. 2020.03.04.976423. DOI: https://doi.org/10.1101/2020.03.04.976423.

78.

Yuan

, He

, and Deng

2019. Compositional data network analysis via lasso penalized D-trace loss. Bioinformatics, 35, 3404–3411.

79.

Yuan

, Xi

, Chen

, et al. 2017. Differential network analysis via lasso penalized D-trace loss. Biometrika, 104, 755–770.

80.

Zengler

, and Zaramela

L.S.

2018. The social network of microorganisms—How auxotrophies shape complex communities. Nat. Rev. Microbiol. 16, 383–390.

81.

Zhang

, and Zou

2014. Sparse precision matrix estimation via lasso penalized D-trace loss. Biometrika, 101, 103–120.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.34 MB