This paper develops a Bayesian graphical model for fusing disparate types of count data. The motivating application is the study of bacterial  communities from diverse high dimensional features, in this case transcripts, collected from different treatments. In such datasets, there are no explicit correspondences between the communities and each correspond to different factors,
making data fusion challenging.
We introduce a flexible multinomial-Gaussian generative model for jointly modeling such count data. This latent variable model jointly characterizes the observed data through a common multivariate Gaussian latent space that parameterizes the set of multinomial probabilities of the transcriptome counts. The covariance matrix of the latent variables induces a covariance matrix of co-dependencies between all the transcripts, effectively fusing multiple data sources. We present a computationally scalable variational Expectation-Maximization (EM) algorithm for inferring the latent variables and the parameters of the model. The inferred latent variables provide a common dimensionality reduction for visualizing the data and the inferred parameters provide a predictive posterior distribution. In addition to simulation studies that demonstrate the variational EM procedure, we apply our model to a bacterial microbiome dataset.[This work was partially supported by grants from ARO W911NF-19-102 and DOE DE-NA0003921. ]

§ INTRODUCTION

We introduce a Bayesian graphical model for jointly modeling and fusing high dimensional count data collected from different sensors with no explicit correspondences between their feature sets. Our model is relevant to the many areas of multi-modality fusion where data is collected from diverse but incommensurate sensor modalities. Examples include multi-view learning in computer vision and automated language translation in natural language processing. However, this paper focuses on a particularly timely application: fusion of microbiome data from diverse microbial communities.

The analysis of meta-transcriptomic data has been of increasing interest to researchers in the biological and health sciences. Microbiomes exist in diverse environments and are critical to sustaining life, balancing ecosystems, and producing antibiotics, among many other functions.  Microbiomes consist of communities of microbes that interact with each other to maintain stability and resilience to environmental conditions and microbial intrusions from competitors. It has therefore been of great scientific interest to quantify changes in microbiome communities due to changing conditions using experimental data. For example, one area of study is the rhizosphere, which is a community of microbial species living around plant root systems, known to be sensitive to environmental factors <cit.>. Another area of study is the spectrum of responses of microbiomes to stressors, collectively called the microbial exposome <cit.>.

One of the principal sensing platforms used to study microbiome communities applies gene sequencing to a microbiome sample, e.g., collected from the gut, the soil, or other environment.
A common way to obtain a global profile of a microbial community is to perform gene sequencing on a biological sample.  For example, RNA-Seq
measures gene expression in a community by quantifying the number of times each gene transcript occurs in the pool of sequenced RNAs. Each microbial species in the community is represented by its own unique set of transcripts, i.e., its transcriptome, and fusing information from different transcriptomes yields the global profile of gene expression across all species in the community. This type of analysis is known as metatranscriptomics and it provides a functional profile of the community that can complement the gene taxonomic profiling provided by metagenomics <cit.>.  This paper introduces a Bayesian graphical model for the metatranscriptomics problem, where inference is performed using a scalable variational EM inference method.  Notably, our model can capture patterns of similarity between histograms of different species' gene expression without inter-species genome-to-genome mappings or knowledge of inter-species transcriptomic pathway correspondences.

A main feature of our model is that it estimates the global covariance structure of gene expression when the observations are in the form of count vectors produced by RNA-Seq. Correlations between transcript abundances are informative about the effect of environmental conditions on microbial communities

In particular, the global covariance matrix captures inter- and intra species interactions. For example, the expression of a single gene in a species can influence other gene expressions in that species or the gene expressions of other community members.

We propose a latent variable graphical model that can capture the hidden factors underlying such dependencies.

The main assumption underlying our proposed model is the existence of a hidden low-dimensional continuous latent space that can explain the observed data. We model the observations as conditionally multinomial distributed given the latent variables, which are assumed to be multivariate Gaussian with low rank covariance structure.

Due to the lack of conjugacy between the Gaussian and multinomial distributions, exact Bayes inference is not tractable. We therefore adopt a Bayes variational inference approach <cit.> to develop an algorithm for estimating the parameters of the proposed model and projecting the data to the latent space.

The proposed model can be contrasted with previously introduced latent variable models used in multi-view learning and dimensionality reduction. Factor analysis (FA) <cit.> is a classical method that is a  generalization of Principal Component Analysis (PCA) <cit.> and Probabilistic PCA <cit.>. FA decomposes the observed data matrix into a low dimensional set of factor loadings and factor scores, imposing a low-rank constraint on the covariance matrix. Like our proposed model, the FA model also assumes a low-dimensional Gaussian latent space but it does not account for the counting nature of the observed data.

Several latent variable models have been proposed for counting observations. These include Latent Semantic Analysis <cit.>, Multinomial PCA <cit.>, and Latent Dirichlet Allocation (LDA) <cit.>. LDA is the most closely related model to the model proposed here since

it is also a Bayesian graphical model for count data and uses multinomial distribution. The main difference is that LDA uses a Dirichlet distributed latent space instead of a Gaussian distributed latent space. Our Gaussian distributed latent space

makes it possible to recover a non-trivial covariance structure among the count variables, unlike LDA <cit.>.

Another way to capture the covariance structure of the observed variables is to ignore the counting nature of the data and use Gaussian Markov random fields (GMRF)  <cit.> to directly estimate the covariance, or Gaussian Graphical Models (GGM) <cit.> to enforce sparsity on the inverse of the covariance estimate.

Unlike our proposed multinomial-Gaussian model the GMRF and GGM do not incorporate a low rank latent structure on the covariance nor do they account for the counting nature of the data.

There have been extensions of the GGM to handle multinomial observations using copulas <cit.> that have been applied to microbiome analysis <cit.>.

Inference in latent variable models, like the one we propose here, can be challenging.

This is especially difficult when there is lack of conjugacy between the distributions of the latent variables and the observed variables. One approach is to perform point estimation for both the latent variables and the parameters in an alternating fashion <cit.>, but this is prone to over-fitting <cit.> and convergence issues. Another approach is to use Markov Chain Monte Carlo (MCMC) methods, which can be computationally expensive <cit.>, especially in high dimension. As an alternative, variational Bayes inference has shown much promise <cit.>. Note that
Variational Bayes

is not a general purpose method and must be tailored to the specific statistical model <cit.>.

When there is a lack of conjugacy, as is the case for the multinomial-Gaussian model in this paper, local variational bound approximations are often adopted <cit.>. Additionally, when there is a problematic expression in the joint density, such as LogSumExp or LogGamma function, which may prevent the inference of the latent variables, surrogate optimization transfer based on Taylor series expansion can be applied to approximate the non-linear function either with linear <cit.> or quadratic <cit.> functions.
We adopt such a local variational bound approach for deriving an inference algorithm for our proposed model.

We summarize our contributions as follows. First, we propose a novel multinomial-Gaussian graphical model to fuse and capture the low rank covariance structure in counting data of disparate types. Our low-dimensional continuous latent space formulation provides dimensionality reduction that can be used for visualization of the count vectors on a common space.

Second, we develop a novel and computationally scalable optimization algorithm based on variational inference to fit the proposed model, which exploits variational local bound approximations. Third, we validate and illustrate the model and its inference algorithm on a synthetic dataset and a real-world bacterial microbiome dataset [The code and the dataset are available at <https://github.com/maktukmak/microbiome-thor>.].

§ PROPOSED MODEL

In this section, we formally define our proposed model and its corresponding variational inference algorithm. Lastly, we discuss computational complexity.

roundnode/.style=circle, draw=green!60, fill=green!5, very thick, minimum size=7mm,
squarednode/.style=rectangle, draw=black!60, very thick, minimum size=5mm,

[latent]       (z) $\bm{z}_{k,i}$;
[obs, above=2cm of z, xshift=0cm]        (x) $\bm{x}_{kl,i}$;

[squarednode, left=2cm of z]        (S) $\bm{\Sigma}_{k}$;
[squarednode, left=2cm of x]    (Th) $\Theta_{k,l}$;

[squarednode, below=2cm of z]        (mu) $\bm{\mu}_{k}$;

z x;

S z;
Th x;

mu z;

plate caption/.append style=above=5pt of #1.north
[inner sep=0.35cm, xshift=-0.0cm, yshift=0cm] zxy (z)(x) $i = 1:I_k $ ;
plate caption/.append style=below=5pt of #1.south west
xTh (Th)(x) $l = 1:L $ ;

plate caption/.append style=below=5pt of #1.south
k (z)(mu)(x)(Th)(S)(zxy)(xTh) $k = 1:K $ ;

Graphical model representation of the proposed latent variable model. $\bm{x}_{kl,i}$ corresponds to the $i$th sample of community $l$ collected from environment $k$. The variables $\{\bm{x}_{kl,i}\}_{l=1:L}$ share a common low-dimensional latent variable $\bm{z}_{k,i}$ that captures the hidden causes of the observations.

§.§ Notation

We denote the $i$th data sample of the $l$th species as $\bm{x}_{kl,i} \in \mathbb{Z}_+^{d_{l}}$, where $k$ indexes the experimental condition, and $d_{l}$ denotes the total number of transcripts for species $l$. The total number of experimental condition from which the samples are collected is denoted as $K$, and the total number of species in the model community is denoted as $L$, hence $l=1:L$ and $k=1:K$. For each experimental condition, different numbers of identically distributed samples are collected. Hence, we denote the total number of samples for the experimental condition $k$ as $I_k$.
The count vector for experimental condition $k$, species $l$, and replicate $i$ is denoted  $\bm{x}_{kl,i}$. The dataset for experimental condition $k$ is $D_{k} = \{ \{\bm{x}_{kl,i}\}_{l=1}^L \}_{i=1}^{I_k}$.

§.§ Latent Variable Model

We model the observed multi-species model community data as generated from a low dimensional latent variable generative model. Under this model the data are conditionally multinomial distributed given the latent variables, which are themselves Gaussian distributed with mean $\bm{\mu}_{k}$ and covariance matrix $\bm{\Sigma_k}$. As will be shown below, the model fuses the observed data across species and it induces a low rank decomposition of the population transcriptome covariance.

Let $\bm{z}_{k,i} \in \mathbb{R}^{d_z}$ be the latent variable assigned for the data sample $D_{k,i}$. $\bm{z}_{k,i}$ thus has the following multivariate normal prior distribution:

\begin{equation} \label{eq:z_pr}
p(\bm{z}_{k,i} | \bm{\mu}_{k} , \Sigma_k) = \mathcal{N} (\bm{z}_{k,i}  | \bm{\mu}_{k} , \Sigma_k),
\end{equation}

where $\bm{\mu}_{k} \in \mathbb{R}^{d_z}$ is the prior mean vector and $\Sigma_k \in \mathbb{S}_{++}^{d_z}$ is the positive definite prior covariance matrix. The observed data consists of count vectors of the transcriptomes, which are modeled as multinomial distributed <cit.>. We model the conditional distributions of the observed count vectors of species $l$ as follows:

\begin{equation} \label{eq:x_pr}
p(\bm{x}_{kl,i} | \bm{z}_{k,i}, \Theta_{kl}) = Mu(\bm{x}_{kl,i} | N_{kl,i}, \mathcal{S}(\Theta_{kl}^T \bm{z}_{k,i}),
\end{equation}

where $Mu$ denotes the multinomial distribution, and $N_{kl,i}$ is the total number of counts of the $i$th data sample of species $l$. Note that we have introduced one more model parameter for each species, specifically $\Theta_{kl} \in \mathbb{R}^{d_{l} \times d_{z}}$, that maps lower dimensional latent space to the higher dimensional observation space of species $l$. Also note that both the latent variable $\bm{z}_{k,i}$ and the parameter $\Theta_{kl}$ are real-valued. Therefore, to provide a proper simplex support set for the multinomial distribution, we use the soft-max function, $\mathcal{S}(\bm{\eta})_d = \exp{\eta_d} / \sum_{d^{'}=1}^D \exp{\eta_{d^{'}}} $. The output of this function is a proper probability vector, i.e., $\sum_{d=1}^D \mathcal{S}(\bm{\eta})_d = 1$ and $\mathcal{S}(\bm{\eta})_d \ge 0$ for all $d = 1:D$. See Fig. <ref> for a graphical representation of the proposed model.

Although it is natural to model the observed counts as multinomial distributed, it may not be obvious why we use  Gaussian latent variables for the latent space. A conjugate distribution such as Dirichlet may seem more natural than the Gaussian distribution, which is not conjugate to Multinomial. However, the components of the Dirichlet distribution are nearly independent <cit.>, hence it is non-trivial to capture the correlations between the hidden components. On the other hand, the Multivariate normal distribution has a covariance parameter which specifically captures the correlation between the hidden components. This is useful for modeling the correlation between multiple datasets. Similar model assumptions are also adopted in topic models <cit.>, categorical PCA <cit.>, and Gaussian process classification <cit.>.

Note that, although the communities are dependent through the latent variables, the experimental conditions are modeled as independent. Hence, there is no coupling between the experimental conditions and thus we fit independent  models for each condition.

The joint log likelihood of the proposed model is of the form $\sum_{k=1}^{K}\sum_{i=1}^{I_k} \log p(\bm{z}_{k,i}, D_{k,i})$, where:

\begin{align}
\begin{split}
\log p(&\bm{z}_{k,i}, D_{k,i}) = \log p(\bm{z}_{k,i}) + \sum_{l=1}^L \log p(\bm{x}_{kl,i} | \bm{z}_{k,i}) \\
= & -\frac{1}{2} \big[(\bm{z}_{k,i} - \bm{\mu}_{k})^T \Sigma_k^{-1} (\bm{z}_{k,i} - \bm{\mu}_{k}) + \log|\Sigma_k|\big]\\
& + \sum_{l=1}^L\sum_{d=1}^D \bm{x}_{kl,id} ( \Theta_{kl,d} \bm{z}_{k,i} - \text{lse}(\Theta_{kl}  \bm{z}_{k,i}) ) + const,
\end{split}
\end{align}

in which lse denotes the log-sum-exp function, i.e., log of the denominator of the soft-max function, and we suppress the deterministic parameters to avoid clutter.
Taking the expectation with respect to $\bm{z}_{k,i}$ is tractable for the linear and quadratic terms, but intractable for the lse term.  We describe an asymptotic approximation in the next section.

§.§ Optimization

Next, we develop a variational EM maximum likelihood algorithm <cit.> to infer the deterministic parameters $\bm{\mu}_{k}$, $\bm{\Sigma}_{k}$, and $\Theta_{kl}$. The main objective is to maximize the likelihood of the observations under the model. The algorithm comprises two alternating steps: i) Expectation step (E-step), where we integrate out the latent variables, ii) Maximization step (M-step), where we optimize the model parameters to maximize the marginal likelihood.

§.§.§ Objective

The proposed model uses Gaussian latent variables for the multinomial observations. Due to lack of conjugacy between Gaussian and Multinomial distributions, the likelihood function is not closed form. Specifically, integrating out the latent variables becomes intractable (See Section <ref> for the details). Hence, we resort to variational inference, in which a lower bound on the likelihood function is derived and maximized. This lower bound is obtained by approximating the posterior distributions of the latent variables. In variational inference, the objective is to minimize the distance (KL-divergence) between the approximate and exact posterior distributions. This objective can be expressed for a single latent variable $\bm{z}_{k,i}$ as follows:

\begin{align}
\begin{split}
\mathbb{KL} (q_{\lambda_{k,i}} | p  ) &= \mathbb{E}_{q_{\lambda_{k,i}}} \log \big[ \frac{q(\bm{z}_{k,i} | \bm{\lambda}_{k,i})}{p(\bm{z}_{k,i} | D_{k,i})} \big] \\
&= \mathbb{E}_{q_{\lambda_{k,i}}} \log \big[ \frac{q(\bm{z}_{k,i} | \bm{\lambda}_{k,i})}{p(\bm{z}_{k,i}, D_{k,i})} p(D_{k,i}) \big] \\
&= \mathbb{E}_{q_{\lambda_{k,i}}} \big[ \log q(\bm{z}_{k,i} | \bm{\lambda}_{k,i}) - \log p(\bm{z}_{k,i}, D_{k,i}) \big] \\
& \;\;\;\; + \log p(D_{k,i}),
\end{split}
\end{align}

where $\bm{\lambda}_{k,i}$ corresponds to the set of parameters of the approximate posterior distribution $q(\bm{z}_{k,i} | \bm{\lambda}_{k,i})$. The expectation operator is defined as $\mathbb{E}_{q_{\lambda}}f(z) = \int f(z) q(\bm{z} | \bm{\lambda}) dz$. Note that the evidence (marginal likelihood) $p(D_{k,i})$ does not depend on $\bm{z}_{k,i}$. Hence, the negative of the expectation term forms a lower bound on the log evidence since the KL distance is always non-negative. This function is known as evidence lower bound (ELBO) and it is the objective function that is maximized in variational EM. The ELBO has the following form:

\begin{equation} \label{eq:obj}
\mathcal{L} = \sum_{i=1}^I \sum_{k=1}^K \mathbb{E}_{q_{\lambda_{k,i}}} \big[ \log p(\bm{z}_{k,i}, D_{k,i}) - \log q(\bm{z}_{k,i} | \bm{\lambda}_{k,i}) \big],
\end{equation}

where the first term in the expectation corresponds to the joint distribution of the latent variable $\bm{z}_{k,i}$ and the associated observed data $D_{k,i}$. The second term corresponds to log of the approximate posterior distribution. The joint distribution has the following form:

\begin{equation} \label{eq:joint}
\log p(\bm{z}_{k,i}, D_{k,i}) = \sum_{l=1}^L \log p(\bm{x}_{kl,i} | \bm{z}_{k,i}) + %\log p(\bm{z}_{k,i} | \bm{u}_{k,i}),
\log p(\bm{x}_{kl,i} | \bm{z}_{k,i}).
\end{equation}

The expressions for $ p(\bm{x}_{kl,i} | \bm{z}_{k,i})$ and $p(\bm{x}_{kl,i} | \bm{z}_{k,i})$ are given in Eq. <ref> and Eq. <ref>, respectively.
We approximate the posterior distribution of $\bm{z}_{k,i}$ as Gaussian with the following form:

\begin{equation}
q(\bm{z}_{k,i} | \bm{\lambda}_{k,i}) =  \mathcal{N}(\bm{z}_{k,i} | \bm{m}_{k,i}, \bm{S}_{k,i}),
\end{equation}

where $\bm{\lambda}_{k,i} = \{ \bm{m}_{k,i}, \bm{S}_{k,i} \}$ is the set of free parameters. Specifically, $\bm{m}_{k,i}$ is the posterior mean and $\bm{S}_{k,i}$ is the posterior covariance. The expectation of the approximate posterior distribution in Eq. <ref> corresponds to the Gaussian entropy function, which has a closed form expression. However, the expectation of the joint distribution is intractable to compute. Next, we present an approximation to resolve the issue.

§.§.§ An upper bound on the LSE

To see why the conditional expectation is intractable, note that the explicit form of the log likelihood of $\bm{x}_{kl,i}$ is a multinomial distribution:

\begin{equation}  \label{eq:mult_ll}
\log p(\bm{x}_{kl,i} | \bm{z}_{k,i}) = \sum_{d=1}^D \bm{x}_{kl,id} ( \Theta_{kl,d} \bm{z}_{k,i} - \text{lse}(\Theta_{kl} \bm{z}_{k,i}) ).
\end{equation}

Taking expectation corresponds to integrating out Gaussian distributed $\bm{z}_{k,i}$. The conditional expectation of the first term is easily determined since it linearly depends on $\bm{z}_{k,i}$. However, the expectation of the second term, which requires integrating $\bm{z}_{k,i}$ over the $\text{lse}$ function, is intractable to compute in a closed form. To overcome this issue, we perform quadratic surrogate optimization transfer, in which a quadratic approximation to the $\text{lse}$ function <cit.> is applied. This results in an upper bound on the multinomial log likelihood. This approximation uses the second order Taylor series expansion with a fixed Hessian matrix.
Particularly, the quadratic upper bound takes the following form:

\begin{equation} \label{eq:lse_app}
\text{lse}(\Theta_{kl} \bm{z}_{k,i}) \le \frac{1}{2} \bm{z}_{k,i}^T \Theta_{kl}^T \bm{A}_{l} \Theta_{kl} \bm{z}_{k,i} - \bm{b}_{kl,i}^T \Theta_{kl} \bm{z}_{k,i} + \bm{c}_{kl,i},
\end{equation}

\begin{equation} \label{eq:a}
\bm{A}_{l} = 0.5 [ \bm{I}_{D_{xl}} - (1/(D_{xl}+1)) \bm{1}_{D_{xl}} \bm{1}_{D_{xl}}^T ]
\end{equation}

is a constant Hessian matrix, whose entries depend only on the dimension of the observation space. The other intermediate parameters $\bm{b}_{kl,i}$ and $\bm{c}_{kl,i}$ are given as follows:

\begin{equation}
\bm{b}_{kl,i}= \bm{A}_{l} \Phi_{kl,i} - \mathcal{S}(\Phi_{kl,i}),
\end{equation}

\begin{equation}
\bm{c}_{kl,i} = \frac{1}{2} \Phi_{kl,i}^T \bm{A}_{l} \Phi_{kl,i} - \mathcal{S}(\Phi_{kl,i})^T \Phi_{kl,i} + \text{lse}(\Phi_{kl,i}),
\end{equation}

where $\Phi_{kl,i}$ is the Taylor series expansion point, which is optimized as a free variational parameter. Note that intermediate parameters are deterministic function of $\Phi_{kl,i}$. Plugging the approximation in Eq. <ref> to Eq. <ref> results in a convex lower bound on ELBO, denoted as $\mathcal{L'}$, which is $\le \mathcal{L}$ and tight at $\Phi_{kl,i}$. Using $\mathcal{L'}$ resolves the intractable integration in Eq. <ref>, resulting in closed form posterior parameter estimates, as described in the next section.

§.§.§ Posterior Distributions - E-step

The E-step in the variational EM algorithm computes approximate posterior distributions of the latent variables, which are subsequently used to compute the expectations in Eq. <ref>. Particularly, there are two parameters to be estimated for each latent variable $\bm{z}_{k,i}$, which are the mean vector $\bm{m}_{k,i}$ and the covariance matrix $\bm{S}_{k,i}$. It is straightforward to maximize over these parameters by using the completing-the-square approach <cit.> (See Appendix <ref>). The terms that quadratically depend on $\bm{z}_{k,i}$ in the joint log-likelihood yield the posterior covariance update:

\begin{equation} \label{eq:post_cov}
\bm{S}_{k,i} = \big[ \Sigma_k^{-1}  +  \sum_{l=1}^{L} N_{kl,i} \Theta_{kl}^T \bm{A}_{l} \Theta_{kl} \big]^{-1},
\end{equation}

where $N_{kl,i}$ is the total number of counts of the $i$th data sample. Similarly, the terms that linearly depend on $\bm{z}_{k,i}$ yield the posterior mean update:

\begin{equation} \label{eq:post_mean}
\bm{m}_{k,i} = \bm{S}_{k,i} \big[ \Sigma_k^{-1} \bm{\mu}_{k}  +
\sum_{l=1}^{L} (\bm{x}_{kl,i} + N_{kl,i} \bm{b}_{kl,i}) \Theta_{kl}
\big].
\end{equation}

Lastly, we update the Taylor series expansion point as:

\begin{equation} \label{eq:post_taylor}
\Phi_{kl,i} = \Theta_{kl} \bm{m}_{k,i}.
\end{equation}

Note that the update of $\Phi_{kl,i}$ depends on the posterior mean. Hence, the algorithm repeats the updates in Eq. <ref>, Eq. <ref>, and Eq. <ref>, respectively, until convergence of the expansion point $\Phi_{kl,i}$.

§.§.§ Point Estimates - M-step

The M-step in the variational EM algorithm maximizes the ELBO with respect to the model parameters. Using the posterior distributions computed in the E-step, we compute the lower bound $\mathcal{L}^{'}$ by taking the expectations with respect to the posterior distributions. Afterwards, taking the derivatives with respect to the model parameters yields closed form update equations for the model parameters. Specifically, the updates for each $\Theta_{kl}$ are given as follows:

\begin{align} \label{eq:Theta}
\begin{split}
\Theta_{kl} = \Big[\sum_{i = 1}^{I_k} (\bm{x}_{kl,i} &+ N_{kl,i} \bm{b}_{kl,i}) \bm{A}_{l}^{-1} \bm{m}_{k,i} \Big] \\
&\Big[ \sum_{i = 1}^I N_{kl,i} ( \bm{m}_{k,i} \bm{m}_{k,i}^T + \bm{S}_{k,i}) \Big]^{-1}.
\end{split}
\end{align}

The update equations for the mean parameter and covariance of the prior distribution of $\bm{z}_{k,i}$ then follow as:

\begin{equation} \label{eq:prior_mean}
\bm{\mu}_k = \frac{1}{I_k} \sum_{i = 1}^I \bm{m}_{k,i},
\end{equation}

\begin{equation}\label{eq:prior_cov}
\Sigma_k = \frac{1}{I} \sum_{i = 1}^I (\bm{m}_{k,i} - \bm{\mu}_k) (\bm{m}_{k,i} - \bm{\mu}_k)^T + \bm{S}_{k,i},
\end{equation}

respectively. The variational EM algorithm is summarized in Algorithm <ref>.

Note that the model parameters of the proposed model are unidentifiable. Due to isotropic Gaussian prior on the latent variables, arbitrary rotation on $\Theta_{kl}$ results in same likelihood <cit.>. This makes direct interpretation of the inferred latent variables ambiguous. Fortunately, this does not affect the predictive performance nor the predictor of covariance, which are the main focus of this paper.

Proposed Variational EM algorithm

Input: $\{D_k\}_{k=1:K}$
Initialize $\{\bm{\mu}_k, \bm{\Sigma}_k, \{\bm{\Theta}_{kl}, \{\bm{\Phi}_{kl,i}\}_{i=1:I_k}\}_{l=1:L}\}_{k=1:K}$
not $\mathcal{L}^{'}$ converged
$k=1$ to $K$
$i=1$ to $I_k$
Infer posterior covariance $\bm{S}_{k,i}$ by Eq.
Infer posterior mean $\bm{m}_{k,i}$ by Eq.<ref>
$l=1$ to $L$
Update variational parameter $\bm{\Phi}_{kl,i}$ by Eq. <ref>
$l=1$ to $L$
Estimate $\Theta_{kl}$ by Eq. <ref>
Estimate $\bm{\mu}_k$  by Eq. <ref>
Estimate $\bm{\Sigma}_k$ by Eq. <ref>
Compute $\mathcal{L}^{'}$ by Eq. <ref>

§.§ Model-predicted Density

A predictor of the population covariance matrix can be extracted from the inferred model. There are two fundamental choices in the proposed model that pave the way to a predictor of the covariance matrix of the transcriptomes from the abundance data. First, we select Gaussian latent space that is common for all species, which models the observation covariance matrix with a low-rank decomposition due to inherent low-dimensional latent space. Second, we adopt a quadratic lower bound on the multinomial likelihood, hence the marginal likelihood of the observations can be approximated with a Gaussian distribution whose covariance matrix reveals the correlations between genomes. Particularly, we define a transformed version of the sample $\bm{x}_{kl,i}$ as $\tilde{\bm{x}}_{kl,i}$ with the following function:

\begin{equation}
\tilde{\bm{x}}_{kl,i} = \bm{A}_{l}^{-1} ( \bm{b}_{kl,i} + \bm{x}_{kl,i} ),
\end{equation}

where $\bm{A}_{l}$ is the matrix defined in Eq. <ref>. Then, it is straightforward to show that the likelihood of the transformed data $\tilde{\bm{x}}_{kl,i}$ is given as follow:

\begin{align} \label{eq:incuded_dist}
\begin{split}
p(\tilde{\bm{x}}_{kl,i}&|\Theta_{kl}, \bm{\mu}_k , \Sigma_k) \\ &= \int \mathcal{N} ( \tilde{\bm{x}}_{kl,i} | \Theta_{kl} \bm{z}_{k,i}, \bm{A}_{l}^{-1} ) \mathcal{N}(\bm{z}_{k,i} | \mathcal{N}(\bm{\mu}_k , \Sigma_k) d \bm{z}_{k,i} \\
&= \mathcal{N} ( \tilde{\bm{x}}_{kl,i} | \Theta_{kl} \bm{\mu}_k, \bm{A}_{l}^{-1} + \Theta_{kl} \Sigma_k \Theta_{kl}^T ),
\end{split}
\end{align}

where the covariance matrix $\bm{C}_{kl, \text{intra}} = \bm{A}_{l}^{-1} + \Theta_{kl} \Sigma_k \Theta_{kl}^T $ and the mean vector $\bm{\phi}_{kl} = \Theta_{kl} \bm{\mu_k}$ are of interest to us, in which $\bm{C}_{kl, \text{intra}}$ captures intra-species correlations of species $l$ in condition $k$. To obtain inter-species correlations, define $\tilde{\bm{A}}^{-1} = \text{diag}(\bm{A}_{1}^{-1}, \dots, \bm{A}_{L}^{-1})$ and $\tilde{\Theta}_{k} =  [ \Theta_{k1}, \dots, \Omega_{kL}]$, then  $\bm{C}_{k, \text{inter}} = \tilde{\bm{A}}^{-1} + \tilde{\Theta}_{k} \Sigma_k \tilde{\Theta}_{k}^T $ gives a covariance matrix for both inter-species and intra-species. To convert any covariance matrix to a proper correlation matrix, which is useful for visualization and analysis,  one can use the transformation $\text{Corr} = \text{diag}(\bm{C})^{-1/2} \bm{C} \text{diag}(\bm{C})^{-1/2} $.

§.§ Computational Complexity

The computational complexity of the variational EM algorithm determines the algorithm's scalability to large datasets. For notational simplicity, we assume that there is only one discrete condition, hence we use $I$ instead of $I_k$. In the E-step, Eq. <ref> computes posterior covariance, which requires multiplication of a $d_z \times d_{l}$ matrix with its transpose resulting $O(d_z^2d_{l})$ complexity. This process is repeated for each species resulting in $O(L d_z^2 d_{l})$. Inverting the matrix for each sample costs $ O(Id_z^3)$. Hence, overall asymptotic complexity for the posterior covariance computation is $O(I ( d_z^3 + L d_z^2 d_{l}))$.
The posterior mean computation in Eq. <ref> involves matrix-vector multiplications that require $O(L d_z d_{l})$, and $O(d_z^2)$ due to covariance posterior covariance multiplication. Hence, the total cost per sample is $O(L d_z d_{l} + d_z^2 )$ and the overall cost is $O(I (L d_z d_{l} + d_z^2))$. Consequently, the complexity of the E-step is $O(I (L d_z d_{l} + d_z^2 + d_z^3 + L d_z^2 d_{l}))$. Removing non-dominant terms results in $O(I (d_z^3 + L d_z^2 d_{l}))$. One can see that this scales linearly in terms of $L$, $d_{l}$, and $I$.

On the other hand, the dominant computation in the M-step is for $\Theta_{kl}$. Eq. <ref> comprises two terms. The first term requires $O(I d_{l} d_z)$ due to $I$ times vector-vector outer products. The second term requires $O(I d_z^2 + d_z^3)$ due to vector-vector outer products and subsequently matrix inversion. Multiplying these terms costs $O(d_{l} d_z^2)$, hence resulting total complexity of $O(L(I d_{l} d_z + I d_z^2 + d_z^3 + d_{l} d_z^2))$ for all $l=1:L$. It is also clear that this computation scales linearly in terms of $L$, $d_{l}$, and $I$. Modeling the conditions independently also induces linear complexity in terms of K. In summary, both E and M steps scale linearly in terms of $K$, $L$, $d_{l}$, and $I$, which suggests that the proposed optimization algorithm is scalable for large datasets as long as the latent space dimension $d_z$ is relatively small.

BIC approximation to the evidence and RMSE of the predicted covariance matrix with respect to the latent space dimension $d_z$. True dimensions are 4, 8, and 12. Blue, orange, and green curves show RMSE and the BIC penalized log likelihood (BIC), respectively. Note that the BIC exhibits a clear maximum over latent space dimension $d_z$. BIC values are scaled by factor $10^{-3}$.

RSME of the covariance estimation with respect to average total number of counts observed in the metatranscriptomic data. As the counts increase the errors decrease until the counts reach a saturation limit.

RMSE of the predicted covariance matrix with respect to the latent space dimension for three different observation space dimensions.

True embeddings

Proposed model embeddings

2D Latent space visualization of 100D count vectors

§ EXPERIMENTS

In this section, we perform numerical experiments to illustrate the proposed model. We start with simulation studies, then conclude with experiments on a bacterial microbiome dataset.

§.§ Simulations

We generate synthetic datasets i) to explain the model selection strategy, ii) to demonstrate the accuracy of the latent embeddings, and iii) to show the ability to capture the covariance structure from observed data.

True Covariance

Proposed Model Covariance

Graphical Lasso Covariance

Emprical Covariance

Ledoit-Wolf Covariance

FA Covariance

Estimated covariance matrices produced by the considered algorithms for a three species community with simulated transcript data. For more details on the model see Section <ref>. The proposed model provides a much more accurate estimated covariance than than do the other methods.

§.§.§ Model Selection

The proposed algorithm estimates the covariance matrix with a low-rank decomposition. The rank of the matrix is equal to the number of components $d_z$ in the latent space. which is a model hyper-parameter to be determined. We use the Bayesian Information Criterion (BIC) to estimate this parameter. The BIC arises from the Laplace approximation to the model posterior $p(M | D_k)$ <cit.>, where $M$ is the complete model including the latent dimension $d_z$. This results in a Bayesian estimate of $d_z$: $d_z = \text{argmax}_{d_z} ( \log p(D_{k}) - \text{BIC}/2)$, where BIC is a function of the total number of unknown parameters, which penalizes the log likelihood with a model complexity penalty term. In the proposed model, we use ELBO lower bound to the the likelihood by following <cit.>. The unknown parameters of the model are  $\{\Theta_{kl}\}_{l=1}^L$, $\bm{\mu}_k$, and $\bm{\Sigma_k}$. Hence, the total number of parameters is $\text{dof} = K \times (d_{l} +  K)$, which is used in the BIC expression as $0.5 \times \text{dof} \times \log(I_k) $. To illustrate the BIC model selection for the proposed model, we simulate three datasets with true latent space dimensions 4, 8, and 12, respectively, and then train multiple models while varying the dimensions $d_z$ over $\{2,3,\dots,12\}$ as the search range. We repeat the experiment 10 times to report the performance. The panel on top of Fig. <ref> shows the average BIC values obtained after convergence of the variational EM algorithm. We see that maximum BIC is obtained on the vicinity of the true ranks for all the datasets. On the other hand, the panel on the bottom shows the RMSE values of the estimated covariance matrices. One can see that the lowest error is achieved by maximizing $d_z$, which is likely to result in over-fitting the model. This over-fitting is mitigated by using BIC estimate of $d_z$.

§.§.§ Embedding Characteristics

We generate a synthetic dataset with a 2 dimensional latent space having 3 different classes, i.e., experimental conditions, according to the model specification in Section <ref>. The latent variables are sampled for each class from different Gaussian distributions. The associated means are predefined as $[0,0], [1.5,1.5], [-1,-1]$ and the variances of the isotropic covariances are selected as $0.5, 0.5, 0.1$, respectively. Three class conditional densities are generated with different affine transformation parameters. The observation space is 25 dimensional. The observations are sampled from the conditional multinomial distributions with soft-max link function as in Eq. <ref>. We generate 200 observations for each class with fixed total number of counts, which is 100, per observation, then stack all the observations. Fig. <ref>.a shows the true embeddings of the resulting dataset. We trained the proposed algorithm with the true latent space dimension. Fig. <ref>.b shows the embeddings of the model, which are obtained through the posterior distributions. Due to non-identifiability of the model, the latent variables can only be recovered up to a rotation. The distorted shape of the latent clusters in  Fig. <ref>.b is due to the use of the soft-max link function. If there is a large component in the affine transformed latent vector, the other components are washed out, hence such points would map to very close points in the observation space. Notwithstanding the differences between Fig. <ref>.a and Fig. <ref>.b, the model preserves the clustering structure accurately.

§.§.§ Influence of the Total Counts and Dimensions on Performance

The number of counts of the observed vector $\bm{x}_{kl,i}$ is an observation-specific parameter, which affects the accuracy of the proposed algorithm. Figure <ref> shows the effect of the number of counts $N_{kl,i}$ on the RMSE values of the covariance estimator under three different latent dimension settings. We sample the total counts of a simulated vector from the Poisson distribution with fixed mean.
We also fix the observation dimension to 128. RMSE is reported based on averaging 20 experiments. Figure <ref> shows that increasing the mean number of counts improves performance. In particular, we see that the total counts $N_{kl,i}$ and the mean error are inversely proportional. This is expected since the number of counts directly affects the posterior uncertainty (Eq. <ref>) and mean (Eq. <ref>). The contribution to the ELBO of the observations increases as the total count increases. Furthermore, for low number of counts, the covariance matrix becomes harder to predict due to higher vulnerability to over-fitting. On the other hand, Fig. <ref> demonstrates the opposite trend when the dimension of  the observation space dimension is increased. Here the total mean counts is fixed to 100. In higher dimensional datasets, the model struggles to estimate the covariance structure when the rank is low. However, this phenomenon diminishes when we observe more counts as can be seen in Fig. <ref>.

§.§.§ Baseline Algorithms

Here we show performance comparisons of the proposed method relative to four baseline methods for estimating the underlying covariance and inverse covariance matrices. i) Empirical covariance, which is the sample covariance. ii) The Ledoit-Wolf estimator <cit.>, which uses shrinkage regularization to perform MAP estimation for the covariance matrix by assigning an inverse Wishart prior on the covariance matrix. iii) Gaussian Copula GraphicalLasso <cit.>, which penalizes the covariance matrix in a different way, particularly, with L1-norm constraints on the precision matrix after transforming the data by using Gaussian copulas. Regularization forces the entries of the precision matrix to be sparse. iv) Factor Analysis <cit.> uses another form of regularization of the covariance matrix by imposing low-rank structure. Each of these baseline methods is expected to perform best when the data, or its transformed version, is normally distributed.

For the Ledoit-Wolf, GraphicalLasso, and FA, as is customary, we first normalize the data by subtracting the mean and dividing by the variance, before running these methods. On the other hand, as the non-Gaussian counting nature of the data is explicitly modeled in our proposed model, our algorithm is run on the raw observations. For model selection in both FA and proposed model, we use the exact rank of the simulated dataset. The regularization coefficient of Gaussian Copula GraphicalLasso algorithm is estimated by using 5-fold cross-validation. For the Ledoit-Wolf algorithm, we used the expression for the shrinkage coefficient given in <cit.>. We use scikit-learn implementations of the baseline methods.

§.§.§ Simulating Model Communities

Next, we generate a synthetic dataset which contains the transcript abundance data (an estimate of gene expression) of two different species existing in the same community, hence $L=2$.  The latent variables $\bm{z}_{i}$ with dimension $d_z = 5$ are generated for each measurement site by sampling from $\bm{z}_{i} \sim \mathcal{N} (\bm{0}_{d_z} , \bm{I}_{d_z})$, where $i$ indexes the replicate for $i=1:I$. These latent variables have elements that correspond to the hidden factors generating the data, such as environmental variables, mediator species effects, and direct associations. We transform the latent variables to the probabilities in the observation space, whose dimensions (abundance of transcripts) are chosen as $d_1 = 20, d_2 = 10$ by using affine and subsequently soft-max transformations as described in Section <ref>. The parameters $\Theta_l \in \mathbb{R}^{d_l\times d_z}$, are chosen randomly by sampling from a zero mean multi-variate normal distribution. Then, we sample the observed data $\bm{x}_{l,i}$ from the multinomial distribution. The total counts $N_{l,i}$ of a sample is chosen randomly by sampling from a Poisson distribution with rate parameter $1000$. We simulate a total of $I = 200$ replicates for each dataset. The true covariance matrix is then given as $\tilde{\Theta}_{l}\tilde{\Theta}_{l}^T$, where $\tilde{\Theta}_{l} =  [ \Theta_{1}, \Theta_{2}]$.

§.§.§ Correlation Results

Fig. <ref> show the estimated covariance matrices of the baseline algorithms, the proposed algorithm, alongside with the ground truth matrix, when simulated dataset is realized. The proposed model can recover the covariance structure accurately. The relatively poorer accuracy of the other methods can be attributed  to several factors. First, these models do not exploit the counting nature of the data. The second reason is that the covariance matrix is simulated with low-rank structure, which is not taken into account by the Gaussian Copula Graphical-Lasso, Ledoit-Wolf, or standard sample covariance estimation methods.
As the data was simulated from the proposed model, the proposed algorithm naturally performs better.

Note also that, for multiple species, the proposed model can discover both inter-species and intra-species correlations. Table <ref> shows the resulting RMSE values between the estimated and the ground truth covariance matrices for the aforementioned simulation setting. The proposed model achieves lower error in overall. This is expected since the model uses an ELBO approximation to the true marginal likelihood function.

Algorithm     Covariance     Precision

Empirical                .231 $\pm$ .013     .510 $\pm$ .021

Ledoit-Wolf              .225 $\pm$ .011     .150 $\pm$ .004

Cov    FA                   .237 $\pm$ .014     .086 $\pm$ .003

GLasso                   .173 $\pm$ .013     .073 $\pm$ .002

Proposed                 .096 $\pm$ .022     .003 $\pm$ .001

Mean and standard deviation of RMSE between the estimated covariance matrices and ground truth over 10 different realizations of the simulated abundance dataset.

§.§ Bacterial Community Experiment

In this section, we demonstrate a real world use-case of the proposed model: transcript analysis of a bacterial model community called THOR <cit.>.

Microbial model communities are useful to understand principles that govern community behaviours <cit.>. The Hitchhikers Of the Rhizosphere (THOR) is a model community consisting of three microbial species, Bacillus cereus, Flavobacterium johnsoniae, and P. koreensis that co-isolate from field-grown soybean roots. The organisms in THOR represent three dominant rhizosphere taxa (at the phylum level), and are common in soil and the mammalian gut.

B. cereus is a Firmicute that carries F. johnsoniae, a member of the Bacteriodetes, and  P. koreensis, a member of the Proteobacteria, as hitchhikers <cit.>. Due to their abundance in several environments, their may demonstrated interactions in the lab and field, and their genetic tractability, these species make a useful model community with relevance to the natural world. The model community provides a simple system in which to study and model community level interactions, which are poorly understood. Developing governing principles of community behavior may lead to strategies to manipulate microbiomes for human or environmental health.

B. cereus

F. johnsoniae

P. koreensis

Effect of koreenciene removal on the the centrality (vertex degree) of vertices in the transcriptional orthology correlation networks inferred from our model for the experimental THOR dataset.  For each species, the ortholog IDs are sorted in decreasing order of the wildtype vertex degree.  The upper row shows plots of the degree of each vertex (transcriptional orthology ID), in descending order of magnitude, for the wildtype condition. The bottom row shows  corresponding plots of the vertex degree when the koreenceine pathway is removed (mutant condition), under the same ordering of vertices as in the top row. P. koreensis preserves its network connectivity better than the other two species. The network connectivity of  F. johnsoniae is the most affected by koreenciene removal.

The dataset is collected under two conditions associated with the treatments applied to P. koreensis. In the first condition the THOR community contains the wild type P. koreensis strain and in the second condition the wildtype is replaced with a mutant of P. koreensis that does not produce koreenceine antibiotics. Production of koreenceines is an important factor in community interactions because they inhibit growth of F. johnsoniae <cit.> and B. cereus protects F. johnsoniae by modulating koreenceine levels.

By using our proposed model, in particular the associated estimated joint probability density of the data, we will be able to reveal effects of the treatment. Since the joint probability density model is parameterized by the mean and covariance of a multivariate Gaussian latent variable (See Section <ref>), the mean and covariance parameters play the principal role in our metatranscriptomic analysis. For brevity we focus our discussion on the inferred covariance parameters here (See Supplementary for discussion of the mean parameters inferred by the model).

The microbial community dataset consists of a total of 17244 gene transcripts associated with three species. There were respectively 38 and 36 replicates for the community with wildtype and mutant strains of P. koreensis. 343 transcripts were removed from the analysis as they had zero counts over all experimental replicates.  After removing these transcripts, B. cereus, F. johnsoniae, and P. koreensis express 5903, 5146, 5852 transcripts, respectively.
We reduced the dimension of the feature space using orthological groupings of gene transcripts into metabolic  pathways[The transcriptional orthology mappings of the THOR gene transcripts to metabolic pathways were obtained using Kegg <https://www.genome.jp/kegg/>. See supplementary for an example.].

Specifically, after pathway mapping each feature corresponds to a transcritpional orthology ID, and the associated data is the summation of the counts of the transcripts tagged with that ID. We aggregated all the transcripts that were not mapped to any Kegg ortholog into a single non-assigned orthology ID, denoted KXXXXX, and we only considered those ortholog IDs that are present in all 3 species. This filtering resulted in a set of 613 ortholog IDs, which corresponds to the dimension of the feature space used in our model.

The rank of the proposed model was determined by successively fitting the model to  latent spaces of dimensions ranging between 5 and 50 with increments of 5.

Then, the optimal model rank was determined as the latent dimension that yields the highest value of the BIC as described in Section <ref>. The optimal model rank was found to be 40. The parameters (mean and covariance) of the models were subsequently refitted with the optimal dimension. The probability distribution of the data is computed under the wildtype and mutant conditions, whose explicit form is given in Eq. <ref> as a marginalization over the latent variables.

Network centrality changes: We evaluate the effect of removal of koreenciene (mutant) on the centrality of the inferred $613 \times 613$ correlation network of metabolic pathways. Here the centrality of a vertex of the network is measured by vertex degree, i.e., the number of edges connecting the vertex.  To ensure that the networks contain only the most biologically significant edges in the networks, we applied a very high correlation threshold ($0.95$) to the respective inferred wild-type and mutant correlation matrices produced by fitting our proposed graphical model to the data. Using such a high threshold is in line with established RNA-Seq network inference practices <cit.>.

Figure <ref> illustrates the effect of removal of koreenceine on the degrees of the nodes (transcriptional orthology IDs) in these networks. Comparison of the upper panels with the lower panels of the figure indicates that the vertex degree distribution of F. johnsoniae is most affected, followed by B. cereus, with P. koreensis the least affected.
This relative ordering of sensitivity of the three species to koreenceine removal shown for vertex degree in Fig. <ref> mirrors the relative ordering of sensitivity shown for the mean changes (See Fig. 2
and associated discussion

in the Supplementary).

B. cereus

F. johnsoniae

P. koreensis

Effect of koreenceine removal on vertex centrality and vertex mean counts for the transcriptional orthology correlation network. For each species, the ortholog IDs are sorted in decreasing order of the vertex degree difference between mutant and wildtype. It is notable that, with few exceptions,  all  orthology IDs with significant changes in vertex mean also have changes vertex degree, but not conversely.  Furthermore, the asymmetry of the blue curve suggests that the removal of koreenceine is associated with an increase in network connectivity (many more vertices whose degrees increase than decrease), especially in F. johnsoniae.

Fig. <ref> illustrates the relative effect of koreenceine removal on increases vs decreases in vertex degree of the transcriptional orthology correlation network for each species.  In the figure the transcriptional orthology IDs are sorted according to the difference between mutant vs wildtype vertex degree. The blue curve shows the resultant vertex degree difference and the orange curve shows the vertex mean difference. Observe that the order of decreasing differences of vertex degree does not correspond to the order of decreasing differences in vertex mean. However, a change in the vertex mean almost always accompanies a change in vertex degree, although the converse is not true.

Also note from the asymmetry of the blue curves in Fig. <ref> that the mutant's networks have many more vertices that increase than decrease in vertex degree as compared to the wildtype.  Thus koreenceine removal seems to increase network centrality of a large number of transcriptional orthologs, especially for F. johnsoniae.
We point out that the large spikes that appear in the orange curves (vertex mean difference) for F. johnsoniae and P. koreensis, correspond to the ID KXXXXX, which are genes that were not mapped to any Kegg  transcriptional ortholog. Further discussion can be found in the supplementary.

In summary, the  proposed model can provide two important data analysis components for microbiome model community analysis. First, we can assess transcriptional orthology composition changes under the treatment by observing the means of the marginal distributions provided by the proposed model. Second, we can assess the second order interaction changes by using the correlation networks that are obtained from the covariance matrices of the marginal distributions. These two components along with the abundance ratio analysis in <cit.> provide a complementary analysis of microbial model communities, which can further be interpreted by microbiologists.

§ CONCLUSION

A hierarchical Bayesian latent variable model was proposed for the joint analysis of multiple discrete datasets. We explained the associations between the features of the datasets with a common lower dimensional latent space, represented by a set of independent identically distributed Gaussian random variables. To overcome the lack of conjugacy between the multinomial observation distribution and the Gaussiam latent space distribution, we developed a variational EM algorithm based on quadratic bound approximations for estimating the parameters in the model. The computation of the algorithm scales linearly with the number of features, samples, and datasets. Simulation studies show that the proposed model can recover low-rank covariance structures accurately. Furthermore, our real-world microbiome experiment demonstrates  the potential real-world utility of the model for exploration of correlation and associated networks for dichotomous microbiome data.

There are several promising directions for future work.
One possible area of future work is to incorporate system dynamics into the latent space so as to explicitly capture temporal correlations. In particular, there is increasing interest in collecting longitudinal microbiome data for studying adaptation, resilience, and dynamics over time.

Incorporation of a state-space dynamical model into our framework can reveal temporal evolution of the interactions between the genomes.

Another future direction is to improve the parsimony of the model by incorporating sparsity into the latent representation by using sparsity inducing priors for the covariance or inverse covariance (precision) matrices.

§ ESTIMATION OF POSTERIOR PARAMETERS

Log-likelihood of Multivariate Normal distribution $\log \mathcal{N}(\bm{x}|\bm{\mu}, \bm{\Sigma})$ can be written as:

\begin{equation*}
-\frac{1}{2}\bm{x}^T\bm{\Sigma}^{-1}\bm{x} + \bm{x}^T\bm{\Sigma}^{-1}\bm{\mu} + const
\end{equation*}

in which the second order term in $\bm{x}$ corresponds to the inverse of covariance matrix $\bm{\Sigma}$, and the linear term corresponds to the mean when multiplied with $\bm{\Sigma}$. Inferring the mean and covariance from linear and quadratic terms is called completing the square approach. We make use of this method to infer the posterior distributions of $\bm{z}_{k,i}$, which is denoted as $q(\bm{z}_{k,i} | \bm{m}_{k,i}, \bm{S}_{k,i})$. Given the joint likelihood in Eq. <ref> and quadratic approximation in Eq. <ref>, one can collect the quadratic terms in $\bm{z}_{k,i}$ as follows:

\begin{equation*}
-\frac{1}{2}\bm{z}_{k,i}^T\bm{\Sigma}_k^{-1}\bm{z}_{k,i} - \sum_{l=1}^L \frac{N_{kl,i}}{2} \bm{z}_{k,i}^T \Theta_{kl}^T \bm{A}_{l} \Theta_{kl} \bm{z}_{k,i},
\end{equation*}

which follows Eq. <ref> for posterior covariance $\bm{S}_{k,i}$ estimate. Similarly, the linear terms are collected as:

\begin{equation*}
\bm{\Sigma}_k^{-1} \bm{\mu}_k \bm{z}_{k,i} + \sum_{l=1}^L \bm{x}_{kl,i}\bm{\Theta}_{kl}\bm{z}_{k,i} + N_{kl,i}\bm{b}_{kl,i}^T\bm{\Theta}_{kl}\bm{z}_{k,i}.
\end{equation*}

Collecting the terms and multiplying with the posterior covariance estimate yields posterior mean $\bm{m}_{k,i}$ estimate as given in Eq. <ref>.

§ SUPPLEMENTARY

Functional orthologies of THOR genomes

Kegg Orthology K00001, an alcohol dehydrogenase

Examples of classification of THOR transcripts based on metabolic functions inferred from gene sequences. Multiple species may have the same transcriptional orthology. For instance, the Kegg functional orthology
K00001 is a group of transcipts that encode alcohol dehydrogenase (AD), an enzyme that plays a role in several cellular processes. AD homologues are found in genomes of

both Bacillus cereus and Flavobacterium johnsoniae.

B. cereus

F. johnsoniae

P. koreensis

Comparison of the effect of koreenceine removal on the mean transcriptional orthology compositions for each of the three species. The ortholog IDs are sorted in decreasing order of the wildtype's Hellinger distance to the mutant composition distribution. Top panel corresponds to wild-type compositions whereas bottom panel shows the mutant compositions, in which the koreenceine pathway of P. koreensis has been removed, with aligned indexes. Hellinger distances between the estimated probability vectors of the treatments are 0.12, 0.18, 0.05 for B. cereus, F. johnsoniae, and P. koreensis, respectively, which suggests that the transcriptional orthology composition of F. johnsoniae is most affected by koreenceine removal.

Mean changes over mutant/wildtype conditions.
For each species, and under either of the wildtype or mutant conditions, we estimate the $613$-dimensional mean vector of Kegg functional orthology counts, as determined by the iterative estimate Eq. (17).
We then transform this estimated mean vector to a probability distribution over IDs using the multivariate softmax transformation, yielding what we call the transcriptional orthology probability distribution that quantifies the composition of orthologs for each species under each condition.
The change to the species-specific probability distribution caused by the removal of koreenceine is quantified by comparing the wildtype and mutant distributions, shown in Fig. <ref>.

To quantify the effect of koreenceine removal, we evaluate the Hellinger distances between the respective wildtype and mutant transcriptional orthology probability distributions. The distances where found to be 0.12, 0.18, and 0.05 for B. cereus, F. johnsoniae, and P. koreensis, respectively. This suggests that the transcriptional orthology composition of F. johnsoniae is most affected under the mutant condition.  In contrast the transcritpional composition of P. koreensis changes the least of the three species.

We also observed that the removal of koreenceine had a substantial effect on the relative abundance of different microbial species. For wildtype and mutant conditions the relative abundances of different species in the model community were computed by taking the total number of RNA-Seq counts for each of the 3 species and normalizing it by its sum. The change in relative abundances was quantified by relative change abundance ratio, which are found to be -0.01 (-1%) for B. cereus, 0.10 (10%) for F. johnsoniae, and -0.08(-8%) for P. koreensis, respectively. These quantities are compatible with <cit.>, which suggests that deleting the koreenceine biosynthetic pathway in  P. koreensis enhances growth of F. johnsoniae and does not affect B. cereus.

It is interesting that while the  removal of koreenceine does not affect the abundance of B. cereus it has substantial effect on the transcriptional orthology probability distribution, as predicted by our model,  discussed in the previous paragraph.

Orthology ID     DegDiff     DegW     DegM     MeanDiff     MeanW     MeanM

Emerging Connections

K15777 (4,5-DOPA dioxygenase extradiol)    	275    	1      	276    	-0.08        	0.10    	0.01

K07646 (OmpR family, sensor histidine kinase KdpD)    	229    	94     	323       0.00       	0.00    	0.00

K03040 (DNA-directed RNA polymerase subunit alpha)    	227    	32     	259       6.16       	1.68    	7.84

Extinguished Connections

K01610 (phosphoenolpyruvate carboxykinase (ATP))    	-132    	135    	3    	-3.29    	4.64    	1.35

K01358 (ATP-dependent Clp protease, protease subunit)    	-136    	137    	1    	-5.88    	6.89    	1.01

K00274 (monoamine oxidase)    	-184    	186    	2    	-0.00    	0.00    	0.00

Up-regulated Functions

K02358 (elongation factor Tu)    	104    	225    	329    	11.57     	23.02    	34.59

K02355 (Not found in KEGG)    	104    	189    	293    	10.22    	10.59    	20.81

K03076 (preprotein translocase subunit SecY)    	165    	78     	243    	7.12    	2.35     	9.48

Down-regulated Functions

KXXXXX (Unannotated)    	139    	207    	346    	-57.08    	536.04    	478.95

K04043 (molecular chaperone DnaK)    	30     	212    	242    	-10.73    	16.28    	5.55

K01358 (ATP-dependent Clp protease, protease subunit)    	-136    	137    	1      	-5.88    	6.89    	1.01

Emerging Connections w/o regulations

K07646 (OmpR family, sensor histidine kinase KdpD)    	229    	94     	323    	0.00    	0.00    	0.00

K02238 (Not found in KEGG)    	161    	151    	312    	0.00    	0.00    	0.00

K00788 (thiamine-phosphate pyrophosphorylase)    	129    	177    	306    	0.00    	0.00    	0.00

Extinguished Connections w/o regulations

K00274 (monoamine oxidase)    	-184    	186    	2     	-0.00    	0.00    	0.00

K03781 (catalase)    	-109    	175    	66    	0.00    	0.00    	0.00

K02575 (MFS transporter, nitrate/nitrite transporter)    	-50     	150    	100    	0.00    	0.00    	0.01

Up-regulated Functions w/o degree changes

K02495 (oxygen-independent coproporphyrinogen III oxidase)    	0    	1    	1    	0.81    	0.16   	0.98

K03530 (Not found in KEGG)    	0    	1    	1    	0.67    	0.22   	0.89

K02899 (large subunit ribosomal protein) L27    	1    	1    	2    	0.64    	0.91   	1.56

Down-regulated Functions w/o degree changes

K03695 (ATP-dep Clp protease ATP-binding subunit ClpB)    	-2    	3    	1    	-1.37    	1.86    	0.48

K00265 (glutamate synthase (NADPH) large chain)    	0     	1    	1    	-0.73    	0.96    	0.22

K03733 (Not found in KEGG)    	0     	1    	1    	-0.65    	0.91    	0.25

Top 3 transcripts of Flavobacterium johnsoniae that change the most according to the eight different treatment effect definitions (emerging connection, extinguished connections, up-regulated functions, down-regulated functions, emerging connections w/o regulations, extinguished connections w/o regulations) . Note that the numbers in the 3 last columns are scaled by $10^3$ to reduce clutter in rendering the table. This table provides a snapshot of the transcriptional orthologs in Flavobacterium johnsoniae that are the most sensitive to removal of koreenceine under 8 different sensitivity criteria.  We highlight Flavobacterium johnsoniae here because it is the species that is most affected although analogous CSV tables for other species can be found in the supplementary files. Emerging and extinguished connections correspond to, respectively, positive and negative centrality changes, as measured by vertex degree, in the transcriptional ortholog correlation graph. Up-regulated and down regulated functions correspond to the orthologies whose compositions change positively and negatively, respectively. Emerging and extinguished connections w/o mean change denote transcriptional orthologs for which there are changes in vertex degree but no change in vertex mean. Up-regulated and down regulated w/o degree changes denote transcriptional orthologs that change have changes in mean but not in vertex degree. The  proposed model can provide two important data analysis components for microbiome model community analysis. First, we can assess transcriptional orthology composition changes under the treatment by observing the means of the marginal distributions provided by the proposed model. Second, we can assess the second order interaction changes by using the correlation networks that are obtained from the covariance matrices of the marginal distributions. These two components  provide a complementary analysis of microbial model communities, which can further be interpreted by microbiologists.