matrix that is obtained by sequencing multiple biological samples and identifying the taxa counts in each sample. However, it is known that microbial associations are impacted by environmental and/or host factors. Thus, a sample-taxa matrix generated in a microbiome study involving a wide range of values for the environmental and/or clinical metadata variables may in fact be associated with more than one microbial network. Here we consider the problem of inferring multiple microbial networks from a given sample-taxa count matrix. Each sample is a count vector assumed to be generated by a mixture model consisting of component distributions that are Multivariate Poisson Log-Normal. We present a variational Expectation Maximization algorithm for the model selection problem to infer the correct number of components of this mixture model. Our approach involves reframing the mixture model as a latent variable model, treating only the mixing coefficients as parameters, and subsequently approximating the marginal likelihood using an evidence lower bound framework. Our algorithm is evaluated on a large simulated dataset generated using a collection of different graph structures (band, hub, cluster, random, and scale-free).
The structure of a microbial community and the organization of its constituent members (taxa) are determined by a combination of their mutual interactions and other factors such as the availability of carbon sources, energy, and nutrients, and the characteristics of the surrounding environment (Falkowski et al., 2008;Hibbing et al., 2010;Williamson and Yooseph, 2012). These variables determine whether the community contains only a handful of taxa with very little strain variation (such as in limiting environments like acid mine drainage (Tyson et al., 2004)) or whether they contain a moderate to high number of taxa with a large number of variants (such as in the human environment (Methe et al., 2012), oligotrophic open oceans (Rusch et al., 2007), and nutrient rich soils (Vogel et al., 2009)).
Microbial community composition can be obtained by sequencing the DNA extracted from a biological sample collected from the environment of interest.
Microbiome sequence data are generated either using a targeted approach, involving the sequencing of a taxonomic marker gene (for instance, the 16S ribosomal RNA gene, which is found in all bacteria (Woese and Fox, 1977)) or using a whole-genome shotgun sequencing approach (Venter et al., 2004); the later approach can be used to deduce both the taxonomic composition and the functional potential of the community.
Here, we study the computational problem of inferring microbial associations from microbiome data. We use the term microbial association to capture both influences and interactions between microbial taxa. In a microbial community, the presence and abundance of one taxonomic group may either directly or indirectly influence the abundance of another taxonomic group (Hibbing et al., 2010). For instance, two microbial taxa may directly influence each other through interactions involving exchange of metabolites or other products, or by competing for the same resources. Alternately, two microbial taxa may not directly communicate or compete for the same resources, but instead one taxon could interact with other members of the community, and these interactions could indirectly influence resource availability for the other taxon. Information about associations between taxa can provide important insights into the ecology of the microbial community.
Microbial associations can be represented using a weighted graph (micro-
bial network) whose nodes represent taxa and undirected edges between nodes represent associations. Edge weights capture the strength of the associations, and the edge weight sign reflects whether the association is positive or negative (Layeghifard et al., 2017). This graph representation can be used to model a variety of microbial interactions, including competition and co-operation (Loftus et al., 2021). Microbial associations can be inferred from the underlying covariance structure of the community which can be calculated using taxa abundances.
The covariance matrix is estimated from a sample-taxa count matrix; this count matrix is generated by sequencing biological samples collected from the environment of interest and identifying the counts of taxa in each sample.
Typically, the study of a microbial community in a particular environment assumes a single covariance structure, and computational methods have been developed to address this estimation problem (Layeghifard et al., 2017), including approaches based on probabilistic graphical models (Biswas et al., 2016;Kurtz et al., 2015;Loftus et al., 2021) and on latent variable models (Fang et al., 2015;Friedman and Alm, 2012). However, with the use of high-throughput next-generation DNA sequencing technologies (Quail et al., 2012) that allow for cost-effectively obtaining data from biological samples, microbiome studies now routinely collect, and generate data from, a large number of samples. In these situations, a microbiome study involving a large cohort or including a wide range of metadata variables (environmental and/or clinical) may in fact be sampling from a community where the microbial associations between taxa are not the same across all intervals of the metadata values. In other words, the microbiome samples in the study may be associated with more than one underlying covariance structure (and thus, more than one microbial network).
Motivated by this scenario, we recently developed an extension to the single network inference problem. In this extension, we treat the inference problem in a mixture model framework based on generative models (Tavakoli and Yooseph, 2019) and solve the following computational problem: given a sample-taxa count matrix generated by a mixture model with K component distributions, estimate the mixing coefficients and the parameters of the K component distributions.
The component distributions model taxa count data, and each component is associated with one precision matrix (and thus, one microbial network). In our framework, referred to as MixMPLN (Tavakoli and Yooseph, 2019), we assume that the taxa counts are generated by Multivariate Poisson Log-Normal (MPLN) distributions (Aitchison and Ho, 1989;Inouye et al., 2017). We estimate the parameters of the MixMPLN model in a maximum likelihood setting using an optimization technique based on the minorization-maximization principle (Lange, 2016). We note that the MPLN distribution has been used previously for the single network inference problem (Biswas et al., 2016;Chiquet et al., 2019). While distributions like the multinomial or the Dirichlet-Multinomial have been popular choices for modeling microbial count data in certain situations (Holmes et al., 2012;La Rosa et al., 2012), these distributions cannot capture both positive and negative associations between taxa. On the other hand, the MPLN distribution can be used to model multivariate count data and its covariance matrix can capture both types of microbial associations.
While a mixture model framework can be used to study the multiple network scenario, in practice however, we do not have a priori knowledge of the value of K, the number of component distributions. In this paper, we propose a variational approximation algorithm to determine the correct value of K for the MixMPLN framework, and solve the model selection problem in a statistically principled manner. As part of our approach, we reformulate the mixture model as a latent variable model (Corduneanu and Bishop, 2001), and treat only the mixing coefficients as parameters, while treating all other variables,
including the means and precision matrices of the component distributions, as latent variables. We use suitable factor distributions involving the latent variables and provide a variational Expectation Maximization (EM) algorithm to compute the parameters of these factor distributions in order to approximate the true marginal likelihood. We evaluate our approach using simulated sampletaxa count matrices generated using different classes of microbial network graph structures.
Notation: Given a matrix X, we use X :i to denote its i th column, X j: to denote its j th row, and x ji to denote its entry in row j and column i. We use n to denote the number of samples, d to denote the number of taxa, and K to denote the number of mixture components. Unless otherwise specified, all vectors are assumed to be column vectors. In the equations below, we associate the variables i, j, and l with samples, taxa, and mixture components respectively.
The MPLN distribution can be used to model count data (Aitchison and Ho, 1989). It has parameters µ and Ω, where µ is the d-dimensional mean vector and Ω d×d is the precision matrix of the distribution. A sample A = (a 1 , ..., a d ) T generated by this distribution is a d-dimensional count vector with the following property: The probability density function p(A|µ, Ω) of the MPLN distribution with parameters µ and Ω can be written as (Aitchison and Ho, 1989):
where,
e -e λ j e λj aj a j ! (2π
(3)
and |Ω| denotes the determinant of Ω. No simplification of the integral in Equation ( 2) is known.
The probability of a sample A = (a 1 , ..., a d ) T generated by a mixture model with K MPLN component distributions and a mixing coefficient vector φ = (φ 1 , φ 2 , .., φ K ) can be written as
where
In Equation ( 4), the variables µ (l) and Ω (l) denote the mean vector and precision matrix respectively, of the l th component distribution, and λ (l) denotes the latent variable vector of sample A that is associated with the l th component. Now, let X d×n denote a sample-taxa count matrix with d taxa and n samples that is generated by mixture of K MPLN distributions. Also, let Λ (l) denote the d × n matrix of latent variable vectors of the n samples that is associated with the l th component. That is, column vector Λ (l):i is associated with sample X :i . We also use λ lji to denote the j th entry in column vector Λ (l):i . Then, the probability of the observed sample-taxa count matrix X, given the parameters of the mixture model, can be written as
This mixture model is associated with K microbial networks, where the l th network (1 ≤ l ≤ K) has adjacency matrix equal to the precision matrix Ω (l) .
We reformulate the mixture model given in Equation ( 5) as a latent variable model (Bishop, 2006), in which we treat only the mixing coefficients φ l 's as parameters while all other variables, including Λ (l) , µ (l) , and Ω (l) , where 1 ≤ l ≤ K, are treated as latent variables.
Let Θ = L ∪ M ∪ T ∪ S denote the set of all latent variables in our model,
where
and S = {S i |1 ≤ i ≤ n}. The set S denotes component membership information for samples, where S i = (s i1 , ..., s iK ) T is a K-dimensional binary vector, also called a 1-of-K binary vector (Bishop, 2006), that is associated with sample X :i .
This vector has the property that if X :i was generated by component r then s ir = 1, and that s il = 0, for all l = r.
We now describe the different parts of the generative model. Each S i is drawn from a multinomial distribution; that is, S i ∼ Multinomial(1, φ). We have that
Conditional on S, each sample is assumed to be independently drawn from an MPLN distribution with parameters µ (l) and Ω (l) . Upon selection of component l, both sample X :i and its associated latent variable vector Λ (l):i are generated.
Thus,
Marginalizing the function p(X, L|M, T , S) × p(S|φ) over S and L results in Equation ( 5). We also introduce conjugate priors over each µ (l) and Ω (l) . We
assume that µ (l) ∼ N d (0, βI) where I is the d × d identity matrix, and β is a fixed parameter. We also assume that Ω l ∼ W(ν, V ), where W(ν, V ) is the Wishart distribution (Wishart, 1928) with fixed degrees of distribution ν and fixed scale matrix V . The density function for the Wishart distribution is given as
where Tr(.) and Γ(.) denote the matrix trace function and Gamma function respectively.
Finally, we set p(M) = K l=1 p(µ (l) ) and p(T ) = K l=1 p(Ω (l) ). Taken together, these specifications allow us to describe the joint distribution of X and all latent variables, conditioned on the mixing coefficients, as
The marginal likelihood function p(X|φ) can be obtained by integrating Equation (7) over all latent variables in Θ. However, the function defined in this manner is not analytically tractable. Instead, we employ a variational approximation method involving a lower bound on the marginal log-likelihood function.
This lower bound, called the Evidence Lower Bound (ELBO) function (Bishop, 2006;Tzikas et al., 2008), will then be maximized with respect to the mixing coefficients. The ELBO function Q(Θ) is defined as
where . Θ denotes the expectation over the distribution Q(Θ). For any function Q(Θ), the following identity holds (Bishop, 2006):
Our goal is to maximize ELBO(Q) using some choice of Q(Θ). The difference between ELBO(Q) and log p(X|φ) can be shown to be the Kullback-Leibler distance (KL) between Q(Θ) and the posterior distribution p(Θ|X, φ). Thus,
is equal to the posterior (Bishop, 2006). Let
is the product of independent factor distributions q(θ t ). With this assumption (Parisi, 1988), the form of the optimal factor distributions that minimize the KL distance can be computed (Bishop, 2006). For each t, the optimal distribution q(θ t ) can be shown to have the form q(θ t ) = exp( log p(X, Θ|φ) θv =θt ) exp( log p(X, Θ|φ) θv =θt ) dθ t (9)
As will be shown, the optimal distributions q(.) for the latent variables S i , µ (l) , and Ω (l) have the same functional forms as their respective priors p(S i |φ), p(µ (l) ), and p(Ω (l) ). Specifically, q(S i ) = Multinomial(1, α i ), with parameter α i = (α i1 , α i2 , .., α iK ) where
), a multivariate Gaussian with d-dimensional mean vector m (l) and d × d precision matrix T (l) , q(Ω (l) ) = W(η (l) , C (l) ), a Wishart distribution with degrees of freedom η (l) and
For the Λ (l):i latent variable vectors, the optimal form of q(Λ (l):i ) is quite unwieldy to work with; thus, instead, we define each q(Λ (l):i ) to be a multivariate
Gaussian distribution with a diagonal precision matrix. Specifically,
Since D (l)i is a diagonal matrix, the multivariate distribution q(Λ (l):i ) can be written as a product of d independent univariate Gaussian distributions q(λ lji ), 1 ≤ j ≤ d, where δ (l)i = (a l1i , a l2i , .., a ldi ) T and diag(D (l)i ) =
That is, a lji and b lji denote the mean and variance respectively, of the random variable λ lji .
We use Equations ( 8) and ( 9) to derive update equations for the parameters of the factor distributions. These update equations are linked, in the sense that the update equation for a variational parameter is a function of other variational parameters. Given these equations, our proposed variational EM algorithm involves maximizing the ELBO function using an iterative procedure. In each iteration, we cycle through the set of update equations to update the values of the variational parameters. While the variational parameters α i , m (l) , T (l) , η (l) , and C (l) have closed forms for their update equations (in terms of other parameters), the update values for variational parameters a lji and b lji are not closed form but are obtained using the Newton-Raphson method (Press et al., 1992). Finally, the mixing coefficients are re-estimated to improve the approximation to the marginal log-likelihood. Convergence of this iterative procedure is guaranteed since the ELBO function increases with each update, unless it is already at a (local) maximum value (Bishop, 2006;Boyd and Vandenberghe, 2004). The algorithm is run on the input sample-taxa matrix using a reasonably large value of K, and after convergence, the number of mixing coefficients that are above a preset threshold value denotes the optimal number of components in the model.
The update equations for the factor distribution parameters α i , m (l) , T (l) , η (l) , and C (l) are given in Table 1. The derivations of the update equations are provided in the Appendix. The estimates of the parameters for each q(Λ (l):i ) are obtained by maximizing the ELBO function restricted to each Λ (l):i .
The ELBO function given in Equation ( 8) when restricted to Λ (l):i has the form
The above integral can be rewritten as
Let ω lrt denote the entry in row r and column t of matrix Ω (l) , and h lt denote the t th entry in vector µ (l) . Since λ lji has a univariate Gaussian distribution, it follows that λ lji = a lji , λ 2 lji = a 2 lji + b lji , and e λ lji = e a lji + 1 2 b lji .
We use these observations to simplify the above expression to produce a function
We identify the values of the variables a l1i , b l1i , a l2i , b l2i , .., a ldi , b ldi that maximize F (.). This is accomplished by cycling through each of the 2d variables and maximizing the corresponding univariate function on that variable.
Restricted to variables a lji and b lji (and excluding constants), the function
As part of maximizing F (.), we first compute the derivatives of G(.) with respect to a lji and b lji separately, and set the two resulting derivatives to 0.
The roots of these two equations are then computed using the Newton-Raphson method. The corresponding equations for a lji and b lji are respectively
Once we have estimates for the variational parameters, we can compute the ELBO function using an expansion of Equation ( 8) as
-log q(S) -log q(M) -log q(T ) -log q(L)
The formulas for the expected values in Equation ( 10) are given in Table 2.
Since the ELBO function approximates the true marginal log-likelihood function log p(X|φ), after we have cycled through and estimated the variational parameters, we can then maximize the resulting ELBO with respect to the mixing
coefficients. This can be done by taking the derivatives of Equation ( 10) with respect to the φ l 's and using a Lagrange multiplier to enforce the constraint that , 2006).
Input: Sample-taxa matrix X d×n , number of components K, the prior parameters β, ν, and V .
Output: Values of the mixing coefficients and the variational parameters that maximize the ELBO function, and the maximum ELBO function value.
Initialization: Initialize the mixing coefficient vector φ = (φ 1 , φ 2 , .., φ K ) and the variational parameters α i , m (l) , T (l) , η (l) , C (l) , δ (l)i , and
Repeat until convergence (that is, the ELBO function does not increase any further):
Cycle through the variational parameters and update their estimates.
We implemented MS MixMPLN in the R programming language (R Development Core Team, 2013). The program is available at https://github.com/syooseph/YoosephLab/tree/master/ MixtureMicrobialNetworks/MS MixMPLN.
The performance of the model selection algorithm MS MixMPLN was evaluated using a collection of synthetic sample-taxa count matrices with d taxa and n
samples. The samples (count vectors) in each sample-taxa matrix were generated from a mixture model consisting of K MPLN component distributions and with mixing coefficient vector φ. The precision matrices of the component distributions were generated from an underlying graph structure. Five different types of graph structures were considered and these were band, cluster, hub, random, and scale-free. The R huge package (Zhao et al., 2012) was used to generate the precision matrices associated with each graph structure. Sample-taxa matrices were generated with number of taxa d = 50, number of components K = 2, 3, 4, and number of samples n = sK, where s is the number of samples per component (s = 200, 1000). The mixing coefficient vectors for K = 2, 3 and 4 were φ
respectively. For each graph type and combination of parameter values, twenty replicates were generated. Thus, 600 synthetic sample-taxa matrices were generated in total.
For each input sample-taxa matrix, MS MixMPLN was run with a larger value for the number of components (5, 6, and 7 respectively for ground-truth K = 2, 3, and 4). The values for the priors were as follows: β = 10 -6 , ν = 50, and V set to a diagonal matrix with all entries equal to 51. MS MixMPLN was run using 26 different starting points on each input, where 25 starting points were generated using random partitions of the samples in the sampletaxa matrix and one starting point was generated using the K-means algorithm (MacQueen, 1967) to partition the samples. Each partition was used to initialize the estimates for the mixing coefficients and the variational parameters. For each starting point, the algorithm was allowed to run for 30 iterations (or less, if it reached convergence). The output with the maximum ELBO value was selected.
The predicted number of components was determined by applying a threshold τ to the mixing coefficient values; that is, any component with mixing coefficient value < τ was not counted towards the predicted number of components. Table 3 shows the accuracy of MS MixMPLN on the simulated datasets at different values for the threshold τ . Here, accuracy for a particular τ value is measured as the proportion of times the predicted number of components at that threshold is equal to the ground-truth K.
We see from Table 3 that, for each combination of n, K, and τ values, MS MixMPLN shows fairly similar accuracy levels for all graph structures.
While the accuracy is highest for K = 2, in the region of 0.75 to 0.95 (at τ = 0.01), it decreases for K = 3 and 4, to smaller values in the region of 0.05 to 0.25 (at τ = 0.01). The accuracy also generally increases with an increase in the number of samples per component (from 200 to 1000). At a higher threshold value (τ = 0.06), even for K = 4, the accuracy estimates are moderately high (in the region of 0.5 to 0.75) for the sample sizes explored.
As an alternate strategy, we evaluated the AIC, BIC, and EBIC criteria (Epskamp and Fried, 2018;Zhu and Cribben, 2018) for model selection. As part of this assessment, we ran the MixMPLN algorithm (Tavakoli and Yooseph, 2019) with different values for K (from 1 through 6) on each sample-taxa matrix, and used the minimum values of AIC, BIC, and EBIC (see Appendix) to predict the number of components. We observed that none of the three criteria returned the correct answer on any of the input datasets.
Here we presented a variational EM algorithm for selecting the number of component distributions for the MixMPLN framework. The proposed algorithm was evaluated using a large simulated dataset. For the sample sizes evaluated, the prediction accuracy decreased as the number of components increased. Future work will explore further the relationship between the number of components and the number of samples in the context of improving approximation of the marginal log-likelihood estimate and the accuracy of the algorithm. It will also include a more comprehensive examination of the parameter space and their effect on model selection accuracy. This will include exploring additional values for the prior parameters and the mixing coefficients. Finally, we will explore the ELBO function landscape further, including the use of local maxima solutions to inform the model selection process. We will also evaluate the proximity of the variational approximation to the true posterior distribution (Yao et al., 2018;Huggins et al., 2020).
Parameter α i for q(S i ):
, where
Parameters m (l) and T (l) for q(µ (l) ):
Parameters η (l) and C (l) for q(Ω (l) ):
Expected values used in the above updates:
Table 1:
The update formulas for the parameters α i , m (l) , T (l) , η (l) , and C (l) . Notation: Tr(.) and ψ(.) denote the matrix trace function and the di-gamma function respectively;
F (θ t ) denotes the expectation of function F (θ t ) over the factor distribution q(θ t ). where the constant is independent of Ω (l) .
From the above equation, we can deduce that q(Ω (l) ) is a Wishart distribution with degrees of freedom η (l) and scale matrix C (l) defined as where log L is the log-likelihood score, K is the number of components, k is the total number of non-zero elements in the precision matrices of the K components, d is the number of taxa, and γ is a constant (set to 0.5).