Year: 2021

Abstract:
Multi-omics data analysis has the potential to discover hidden molecular interactions, revealing potential regulatory and/or signal transduction pathways for cellular processes of interest when studying life and disease systems. One of critical challenges when dealing with real-world multi-omics data is that they may manifest heterogeneous structures and data quality as often existing data may be collected from different subjects under different conditions for each type of omics data. We propose a novel deep Bayesian generative model to efficiently infer a multi-partite graph encoding molecular interactions across such heterogeneous views, using a fused Gromov-Wasserstein (FGW) regularization between latent representations of corresponding views for integrative analysis. With such an optimal transport regularization in the deep Bayesian generative model, it not only allows incorporating view-specific side information, either with graph-structured or unstructured data in different views, but also increases the model flexibility with the distribution-based regularization. This allows efficient alignment of heterogeneous latent variable distributions to derive reliable interaction predictions compared to the existing point-based graph embedding methods. Our experiments on several real-world datasets demonstrate enhanced performance of MoReL in inferring meaningful interactions compared to existing baselines.

heading: 1 INTRODUCTION
text: Multi-view learning tries to fully leverage the information from multiple sources (i.e. different types of omics data in molecular biology) and represents them in a shared embedding space, which is beneficial for many downstream tasks with a limited number of training samples. In biomedical applications, the shared embedding space also enables better understanding of the underlying biological mechanisms by discovering interactions between different types of molecules, which is our focus in this paper.
Existing multi-omics data integration methods are limited in their applicability. First, most of them attempt to derive low-dimensional embeddings of the input samples and are not designed to infer a multi-partite graph that encodes the interactions across views. In unsupervised setting, matrix factorization based methods, such as Bayesian Canonical Correlation Analysis (BCCA) (Klami et al., 2013) and Multi-Omics Factor Analysis (MOFA) (Argelaguet et al., 2018), can achieve the similar goal of cross-view relational learning but often through two-step procedures, in which the factor loading parameters are used for downstream interaction analyses across views. Second, a very recent relational inference for multi-view data integration, BayRel (Hajiramezanali et al., 2020), is built on three strict assumptions, which may limit its practical application, including in multi-omics data integration: 1) A graph of dependency between features of each view is available; 2) Input dataset is complete on all views with no missing samples; 3) The samples in different views are well-paired. While the first limitation might be solved by learning a graph using an ad-hoc technique, the last two issues are common in many multi-omics data integration problems. Integrated samples commonly have one or more view with various missing patterns. This is mostly due to limitations of experimental designs or compositions from different data platforms. In addition, data might be collected in different laboratories or the sample IDs are not available due to patient identification or security concerns, leading to unpaired datasets. Apart from these, we might not have access to a priori graph structure data in some view(s) as the nature of data might not be structured, or we only have incomplete or very noisy prior knowledge. For such multi-omics data, leaving out such a view may lose some complementary information while enforcing graph structures may cause degraded performance.
In this work, we propose a new Multi-omics Relational Learning method, MoReL, based on the fused Gromov-Wasserstein (FGW) regularization, mitigating the dependency of multi-view learning on the aforementioned two assumptions. The proposed method contains four major contributions: 1) MoReL provides a new Bayesian multi-omics relational learning framework with efficient variational inference and is able to exploit non-linear transformations of data by leveraging deep learning models for either unstructured or graph-structured data; 2) MoReL learns a multi-partite graph across different features from multiple views using a FGW-based decoder, facilitating meaningful biological knowledge discovery from integrative multi-omics data analysis while accounting for arbitrarily permutation and/or transformation caused by processing features with different deep functions across the views; 3) MoReL can flexibly integrate both structured and unstructured heterogeneous views in one framework, in which only confident constraints need to be imposed to improve the model performance; 4) MoReL is able to integrate multiple views with unpaired samples and/or arbitrary sample-missing patterns.

heading: 2 RELATED WORKS
text: Optimal transport. There have been extensive efforts to utilize Gromov-Wasserstein (GW) discrepancy to solve the alignment problems in shape and object matching (Mémoli, 2009; 2011). A similar attempt has been made recently to investigate its potential for more diverse applications, such as aligning vocabulary sets between different languages (Alvarez-Melis & Jaakkola, 2018), and graph matching (Chowdhury & Mémoli, 2019; Vayer et al., 2018b; Xu et al., 2019b). Peyré et al. (2016) have proposed a fast Sinkhorn projection-based algorithm (Cuturi, 2013) to compute the entropy-regularized GW distance. Following this direction, Xu et al. (2019b) have replaced the entropy regularizer with a Bregman proximal term. To further reduce the computational complexity, the recursive GW distance (Xu et al., 2019a) and the sliced GW distance (Vayer et al., 2019) have been proposed. In Bunne et al. (2019), a pair of generative models are learned for incomparable spaces by defining an adversarial objective function based on the GW discrepancy. It imposes an orthogonal assumption on the transformation between the sample and its latent space. However, it can not incorporate the graph structured data. Similar to our model in this paper, Vayer et al. (2018a) and Xu et al. (2020) have proposed to impose the fused GW regularization in their objective functions by combining GW and Wasserstein discrepancies.
Graph CCA (gCCA). In order to utilize a priori known information about geometry of the samples, gCCA methods (Chen et al., 2019; 2018) have been proposed to construct a dependency graph between samples and directly impose it into a regularizer. Similar to classical CCA, gCCA learns an unstructured shared latent representation. Unlike our MoReL, though, they can neither take advantage of the dependency graph between features, nor explicitly model relational dependency between features across views. Therefore, they rely on ad-hoc post-processing procedures as a second step to infer inter-relations.
Graph representation learning. Graph neural network architectures have been shown to be effective for link prediction (Hamilton et al., 2017; Kipf & Welling, 2016; Hasanzadeh et al., 2019) as well as matrix completion for recommender systems (Berg et al., 2017; Monti et al., 2017; Kalofolias et al., 2014; Ma et al., 2011). The first group of models is dealing with a single graph and is not able to deal with heterogeneous graphs, with multiple types of nodes and edges, and node attributes (Zhang et al., 2019). The second group utilizes the known item-item and user-user relationships and their attributes to complete the user-item rating matrix. However, they rely on two strict assumptions: 1) the inter-relation matrix is partially observed; and 2) both views have structured information. The proposed MoReL achieves robust multi-view learning without these assumptions, making it more practical in multi-omics data integration.

heading: 3 PRELIMINARIES

heading: 3.1 WASSERSTEIN DISTANCE
text: Wasserstein distance (WD) is a measure for comparing probability distributions on the same or different spaces where a meaningful distance across domains can be computed (Vayer et al., 2019). Given two distributions Λ and ∆ defined on the corresponding space X and Y, and a transportation
cost c : X× Y→ R+, WD is the solution to the following optimization problem:
inf π∈Π E(x,y)∼π[c(x, y)] = inf π∈Π
∫ c(x, y) dπ(x, y),
where π ∈ Π(X×Y) is the transport map such that its marginals are Λ and ∆, respectively. Assuming that the probability distributions are discrete, i.e. Λ = ∑n i=1 aiδxi and ∆ = ∑m j=1 bjδyj with δ as the Dirac delta function, WD optimization could be simplified as follows:
DW(Λ,∆) = min T∈Π(a, b) n∑ i=1 m∑ j=1 Ti,j c(xi, yj),
where Ti,j is an element of the transport matrix T whose row-wise and column-wise sums equal to [ai] n i=1 and [bj ] m j=1, respectively.

heading: 3.2 GROMOV-WASSERSTEIN DISTANCE
text: Gromov-Wasserstein distance (GWD) has been proposed as a natural extension of WD when a meaningful transportation cost between the distributions cannot be defined for WD. For example when two distributions are defined in Euclidean spaces with different dimensions (Vayer et al., 2019). Instead of measuring inter-domain distances, GWD measures the distance between pairs of samples in one distribution and compares it to those in the other domain. More specifically, given two distributions Λ and ∆ defined on spaces X and Y, as well as two domain-specific transportation costs c(X) and c(Y), GWD is the solution to the following optimization problem:
inf π∈Π
E(x,y)∼π,(x′,y′)∼π[L(x,x ′,y,y′)] = inf
π∈Π
∫ ∫ L(x,x′,y,y′) dπ(x, y) dπ(x′, y′),
where L(x, x′, y, y′) =‖ c(X)(x, x′) − c(Y)(y, y′) ‖ and π ∈ Π(X × Y) is the transport map. Likewise, this can be derived for discrete distributions Λ = ∑n i=1 aiδxi and ∆ = ∑m j=1 bjδyj , as follows:
DGW(Λ,∆) = min T∈Π(a, b) n∑ i,i′=1 m∑ j,j′=1 Ti,j Ti′,j′ L(xi, xi′ , yj , yj′), (1)
where Ti,j is an element of transport matrix T whose row-wise and column-wise sums equal to [ai] n i=1 and [bj ] m j=1, respectively.

heading: 4 METHOD

heading: 4.1 PROBLEM FORMULATION AND NOTATIONS
text: We propose a novel hierarchical generative model for multi-omics data integration that incorporates view-specific structure information when it is available. Given observations from structured and unstructured views, our model, Multi-omics Relational Learning (MoReL), aims to infer the interrelations among entities, i.e. features, across all of the views. More specifically, assume that multiple
views, V , of data are given. Without loss of generality, we assume that the structure information, provided as a graph, is available for some of the views Vs ⊂ V , and the remaining views Vu = V \Vs are unstructured. We note that every structure could be represented as a graph. For example, image and sequential data could be represented over grid and directed path graphs, respectively.
We represent the set of graphs for structured views by Gs = {G(v)}v∈Vs and their adjacency matrices by As = {A(v)}v∈Vs . We also define Xs = {X(v)}v∈Vs as the set of node attributes for structured views, and Xu = {X(v)}v∈Vu as the set of data for unstructured views. Moreover, Nv denotes the number of nodes in structured views and number of features for unstructured views. MoReL infers the interactions among the nodes in Gs and features in Xu. We represent these inter-relations by a multi-partite graph with ∑ v∈V Nv nodes and a multi-adjacency tensor A = {A(vv
′)}v,v′∈V,v 6=v′ , where A(vv ′) is the Nv ×Nv′ bi-adjacency matrix between views v and v′.

heading: 4.2 MOREL GENERATIVE MODEL
text: We define a hierarchical Bayesian model for MoReL with three sets of latent variables: 1) H = Hs ∪ Hu = {H(v)}v∈Vs∪Vu , which captures the (hidden) structural information; 2) A, which encodes the interaction among features across views; and 3) Z = Zs ∪Zu = {Z(v)}v∈Vs∪Vu , which summarizes the feature/attribute specific information. The joint probability of observations and latent variables factorizes as follows:
pθ(Xu, Xs, As, H, A, Z) = pθx(Xu | Zu) pθx(Xs | Zs) pθg (As |Hs) pθz (Z |H, A) pθa(A |H) p(H).
(2)
Figure 1 depicts the generative model of MoReL with structured and unstructured views. In the following subsections, we define different parts of the generative and inference model.

heading: 4.2.1 OPTIMAL TRANSPORT FOR MULTI-PARTITE GRAPH DECODER
text: In this subsection, we define the generative distribution of the multi-adjacency tensor, A. We note that inferring A is the main goal of our model. Given the structural latent variablesH, we introduce a fused Gromov-Wasserstein (FGW) distance based mapping to generate A. FGW refers to distance metrics defined by combining WD and GWD, which has been proposed to compare structured distributions (Vayer et al., 2018b; Chen et al., 2020). Considering graphs with node attributes as structured distributions, WD compares node distributions in two graphs (i.e, node similarity), GWD measures the distance between pairs of nodes in one graph and compares it to those in the other (i.e., edge/path similarity).
FGW distance. Given two sets of structured distributions, Λ and ∆, FGW is defined as follows: DFGW(Λ, ∆) =αDW(Λ, ∆) + βDGW(Λ, ∆)
=α inf πw∈Π
E(x,y)∼πw [c (XY)(x, y)]
+ β inf πgw∈Π
E(x,y),(x′,y′)∼πgw [‖ c (X)(x, x′)− c(Y)(y, y′) ‖],
(3)
where α, β ∈ [0, 1] are scalar hyper-parameters, Π is the set of all admissible couplings between Λ and ∆, and c(XY), c(X), and c(Y) are corresponding transportation cost functions. DFGW can be further simplified by choosing πw to be equal to πgw (Chen et al., 2020).
Relational learning via FGW. We are interested in aligning the nodes/features in every pair of views, i.e. (v, v′). Hence, we will have a FGW distance based decoder for every pair of views, in which each view independently belongs to either structured or unstructured views, i.e. (v, v′) ∈ V . To that end, we first define the transportation cost functions c(vv
′) and c(v), and then approximate DFGW. We define the (inter-)cost function for the first term of FGW, i.e. DW, as follows:
c(vv ′)(H (v) i,: , H (v′) j,: ) = 1− σ ( H (v) i,: (H (v′) j,: ) T ) ; v, v′ ∈ V, (4)
where σ denotes the sigmoid function, and H(v)i,: represents the structural latent variable of node/feature i in view v. To calculate the DGW, we define two different transportation costs based
on the nature of the inputs. For the structured views, we define the cost function as a combination of the shortest path distance from graph and the distance between structural latent variables. More specifically, given the normalized shortest path distance matrix between every pair of nodes in the input graph D(v):
c(v)(H (v) i,: , H (v) j,: ) = D
(v) ( 1− σ ( H
(v) i,: (H (v) j,: )
T ))
; for v ∈ Vs,
where denotes the Hadamard product. This construction ensures both graph and attributes information are incorporated in the distance function. For unstructured views, we define the cost function between two features as follows:
c(v)(H (v) i,: , H (v) j,: ) = 1− σ
( H
(v) i,: (H (v) j,: )
T )
; for v ∈ Vu.
Noting the definitions of WD and GWD in Section 3, we rewrite DFGW between two views of data with shared transport matrixas follows:
DFGW ( p(H(v)), p(H(v ′)) ) =
Nv∑ i=1 Nv′∑ j=1 min T (vv′) gw ∈Π ∑ H
(v) i,: ,H (v′) j,: ,H (v) i,: ′,H (v′) j,: ′
[ α c(vv
′)(H (v) i,: , H (v′) j,: ) +
β ‖ c(v)(H(v)i,: , H (v) i,: ′)− c(v ′)(H (v′) j,: ,H (v′) j,: ′) ‖ ] .
(5)
To approximate the FGW distance, we first deploy GW algorithm in equation (1) to obtain T(vv ′)
gw
and DGW, and then utilize T(vv ′) gw along with the defined transportation cost c(vv ′) to calculate Wasserstein distance term in DFGW (Chen et al., 2020). The pseudo-code in Algorithm 1 in the Appendix provides the details of the FGW distance calculation procedure. Please note that we use the same Sinkhorn solver as in Chen et al. (2020) and Alvarez-Melis & Jaakkola (2018).
We further can generate A for every pair of views based on T(vv ′) gw as follows:
p(A |H) = ∏
v,v′∈V v 6=v′
p(A(vv ′) |H(v),H(v ′)) = ∏
v,v′∈V v 6=v′
Ber ( A(vv ′) | γT(vv ′) gw /max(T (vv′) gw ) ) , (6)
where γ ∈ [0, 1] is a normalizing hyper-parameter, and Ber is short for Bernoulli. We note that the sum of the elements in each of the transport matrices T(vv ′) gw equals to one. Hence each of its elements has a small value. Therefore, we normalize the transport matrices (as γT(vv ′) gw /max(T (vv′) gw )) to avoid very sparse and trivial solutions. To use the reparametrization trick during training, we sample from concrete relaxation of Bernoulli (Gal et al., 2017). We emphasize that our proposed FGW-based decoder is the key in aligning features/nodes across structured and unstructured views via accurate and efficient distribution matching scheme.

heading: 4.2.2 PRIOR CONSTRUCTION AND LIKELIHOODS
text: Prior. We impose independent zero-mean unit-variance Gaussian priors on elements of H. The prior for Z is a multivariate Gaussian distribution whose mean and diagonal covariance matrix are constructed from the inferred multi-partite graph and the structural latent variable H. We use two graph neural networks (GNNs) g(µ)pz and g (σ) pz to map H and A to the parameters of pθz (Z). Specifically,
pθz (Z |H,A) = ∏ v∈Vs Nv∏ i=1 pθz (Z (v) i,: |H,A); pθz (Z (v) i,: |H,A) = N (µ (v,i) pz ,σ (v,i) pz ),
with [µ(v,i)pz ]v,i = g (µ) pz (H,A), [σ(v,i)pz ]v,i = g(σ)pz (H,A).
We note that in this setting,H is considered as node attributes of the multi-partite interaction graph.
Likelihood of observations. To reconstruct the input graphs in the structured views, we assume that the views and edges are conditionally independent. More specifically, we employ an inner-product decoder as follows:
pθg (As |Hs) = ∏ v∈Vs Nv∏ i,j=1 pθg ( A (v) i,j |H (v) i,: , H (v) j,: ) ;
pθg ( A(v)i, j |H(v)i,: ,H (v) j,: ) = Ber ( σ(H (v) i,: H (v) j,: T ) ) .
To generate the features in unstructured views and node attributes in structured views, we assume that views are conditionally independent. Hence we can expand the feature reconstruction terms in the the equation (2) as follows:
pθx(Xu | Zu) = ∏ v∈Vu pθx(X (v) |Z(v)), pθx(Xs | Zs) = ∏ v∈Vs pθx(X (v) |Z(v)).
We note that pθx could also be view specific depending on whether the node attributes/features in a view are discrete or continuous. In our experiments, we have deployed the Gaussian likelihood with the unit variance. The mapping from Z to the parameters of pθx(X ), in our case, the mean of the Gaussian distribution, can be any highly expressive function such as neural networks. We denote these functions by f (v,s)px and f (v,u) px .

heading: 4.3 INFERENCE NETWORK AND LEARNING
text: Posterior. We model the posterior of the structural latent variables as a Gaussian distribution and infer its parameters independently for each view. More specifically,
qφh(Hu | Xu) = ∏ v∈Vu qφh(H (v) |X(v)), qφh(Hs | Xu, As) = ∏ v∈Vs qφh(H (v) |X(v), A(v)).
We use two GNNs for each structured view, {g(µ,v)qh (X(v),A(v)), g (σ,v) qh (X (v), Ab(v))}v∈Vs , and two fully connected neural networks per unstructured view, {f (µ,v)qh (X(v)), f (σ,v) qh (X
(v))}v∈Vu , to map inputs to the mean and variance of the posteriors. We consider the variational distribution of Z to be a multivariate Gaussian distribution, and it is factorized as follows:
qφz (Zu | Xu) = ∏ v∈Vu qφz (Z (v) |X(v)), qφz (Zs | Xu, As) = ∏ v∈Vs qφz (Z (v) |X(v), A(v)).
We use two GNNs per structured view, {g(µ,v)qz (X(v), A(v)), g(σ,v)qz (X(v),A(v))}v∈Vs , and two fully connected neural networks for each unstructured view, {f (µ,v)qz (X(v)), f (σ,v)qz (X(v))}v∈Vu , in the same fashion as qφh to infer parameters of qφz .
Objective function. Having defined the prior and posterior distributions as well as the likelihood, we write the overall loss function as the sum of the negative variational ELBO and FGW regularization terms. Specifically,
L =− ELBO + LFGW = Eqφz (Zu | Xu)log pθx(Xu | Zu) + Eqφz (Zs | Xs)log pθx(Xs | Zs)
+ Eqφz (Zu,Hu | Xu) log pθ(Zu | A,H) + Eqφz (Zs,Hs | Xs,As) log pθ(Zs | A,H) − Eqφz (Zu | Xu) log qφz (Zu | Xu)− Eqφz (Zs | Xs,As) log qφz (Zs | Xs, As) − Eqφh (Hu | Xu) log qφh(Hu | Xu)− Eqφh (Hs | Xs,As) log qφh(Hs | Xu, As)) + Eqφh (Hu | Xu) log p(Hu) + Eqφh (Hs | Xs,As) log p(Hs)
+ ∑ v∈V ∑ v′∈V v′ 6=v DFGW ( p(H(v)), p(H(v ′)) ) .
(7)
While, as mentioned previously, we use the Sinkhorn algorithm to calculate the DFGW, the overall loss is optimized using stochastic gradient descent based optimization algorithms such as Adam (Kingma & Ba, 2014).

heading: 5 EXPERIMENTS

heading: 5.1 DATASETS AND EVALUATION METRICS
text: Datasets. We use the same datasets, i.e. microbiome-metabolite interactions in cystic fibrosis (CF) and gene-drug interactions in precision medicine, as BayReL. Dataset description and graph construction procedure are detailed in the Appendix. We want to emphasize that although these datasets have structured views, have no missing samples and their samples are completely paired, in many real-world cases these assumptions are not satisfied. These datasets were chosen merely to get a better understanding of the advantages of MoReL specially compared to BayReL. We evaluate MoReL in different settings. More specifically, we demonstrate the performance of MoReL when: 1) one or both views are unstructured, 2) there are missing samples, and 3) samples are not paired. Furthermore, we have a comprehensive comparison with BayReL when both views are structured.
Evaluation metrics. To quantify the performance of the methods, we use the same evaluation metrics as the ones introduced in BayReL (Hajiramezanali et al., 2020). Since in these datasets, the true negatives, i.e. non-interactions, are not known; and there are only a small subset of true positives, i.e. true interactions, well-known classification metrics cannot be used for evaluation. Therefore, positive accuracy and negative accuracy have been defined to evaluate microbiome-metabolite experiments. Positive accuracy refers to the accuracy of identifying validated interactions with P. aeruginosa. Negative accuracy exploits the fact that there should not be any common metabolite targets between known anaerobic microbes (Veillonella, Fusobacterium, Prevotella, and Streptococcus) and notable pathogen P. aeruginosa. Let B denote the set of all microbes and A1 and A2 represent two disjoint sets of metabolites. Negative accuracy is defined as 1 − ∑ i∈A1 ∑ j∈A2 ∑ l∈B 1(i and j are connected to l)
|A1|×|A2|×|B| , where 1(·) is the indicator function. Having both higher positive and negative accuracy is desired. For precision medicine, we compare the prediction sensitivity of identifying known interactions in the test sets while tracking the average density of the overall constructed graphs. We note that inferring very dense graphs would lead to high prediction sensitivity as it will includes most of the possible interactions. Therefore, tracking the sparsity of the inferred graphs is the key to properly evaluate the models’ capability in predicting meaningful interactions.

heading: 5.2 BASELINES AND EXPERIMENTAL SETUPS
text: Baselines. We compare MoReL with three baselines including Spearman’s Rank Correlation Analysis (SRCA), BCCA (Klami et al., 2013), and BayReL (Hajiramezanali et al., 2020). While SRCA applies to raw data, BCCA first finds low-dimensional latent representations of views via matrix factorization and then the interactions are discovered based on the correlation between representations. While BCCA and SRCA could not incorporate the structure of data and need a two-step procedure to infer the interaction between features across the views, BayReL is able to use the structured of data and infer the relations without any ad-hoc post-processing procedure. However, BayReL suffers from three strict assumptions: 1) all views of data are structured; 2) there are no missing samples in any views; and 3) samples are paired, i.e. the id of samples are known. We re-emphasize that MoReL is the very first model that not only can infer interactions across structured and unstructured views but also is able to handle missing and unpaired samples in different domains, making it more applicable in real-world multi-omics data integration A widely used method for multi-omics data integration is MOFA. The mathematical modeling of MOFA is the same as BCCA except for the data likelihood part. While BCCA only supports continuous data, MOFA can have discrete likelihoods such as Poisson. Since our datasets do not have discrete features, we are only reporting BCCA results.
Hyper-parameters. In all of our experiments, to have a fair comparison, architectural hyperparameters (i.e. number of layers and number of neurons) were set to be the same as in BayReL. Other hyper-parameters that are unique to MoReL were tuned using the validation set. More specifically, the number of hidden layers as well as their dimensions are the same for the corresponding functions in both structured and unstructured views. We use graph convolutional layers (Kipf & Welling, 2017) for structured views and fully connected layers for unstructured views except for reconstructing X from Z , for which we use fully connected layers in all of the views. The mapping from inputs to the mean and variance parameters ofH are two 2-layer neural networks (16 and 8 dimensional layers) with a shared first layer for each view. We use two 2-layer neural networks (16 and 8 dimensional
layers) with a shared first layer for each view for the mapping fromH to the mean and variance of Z . We use a 3-layer fully connected neural network (8 and 16 dimensional hidden layers) for each view as the reconstruction function mapping Z to X . The temperature for relaxed Bernoulli distribution is set to 0.3. The normalizing parameter γ in equation 6 is 0.9 while α and β in DFGW are set to 1 and 0.5, respectively. We used the exponential decaying learning rate with the decay rate of 0.01 and initial learning rate of 0.01. All of our results are averaged over multiple runs with different random seeds. We have implemented MoReL and all the competing methods in Tensorflow (Abadi et al., 2015). All the experiments are performed on a workstation with a single NVIDIA P100 GPU.

heading: 5.3 DISCUSSION, DATASETS WITH UNSTRUCTURED VIEWS
text: Table 1 shows the performance of three variants of MoReL and competing methods for microbiomemetabolite data integration with the CF data. In these experiments, we assume that the samples are paired and all are available in both views. In MoReL uu, we report the results when both views are unstructured. In MoReL us, we have the graph of interactions between microbiomes while the metabolite view is assumed to be unstructured. Comparing MoReL uu and baselines that do not incorporate any graph-structured data as input, we observe an almost 30% improvement in positive accuracy while maintaining higher than 97% negative accuracy. This demonstrates that our proposed MoReL, even without any structural information, is effective in inferring meaningful interactions. Further incorporating the network between microbiomes (i.e. MoReL us) leads to a 37% and 7% improvement compared to the baselines and MoReL uu, respectively. This shows not only the importance of incorporating view-specific side information, but also the effectiveness of FGW-based decoder in aligning structured and unstructured views. Further results on interpretability and robustness of MoReL on CF dataset is provided in the Appendix.
The results for prediction sensitivity of two variants of MoReL and competing methods in the precision medicine experiments are shown in Table 2. We observe that both MoReL uu, where both views are unstructured, and MoReL us, where the graph structure between genes is given, consistently outperform the baselines by a significant margin using graphs of different densities. This, once again, proves that MoReL is able to learn meaningful relations both in sparse and dense graphs. Comparing the results for MoReL us and BCCA, the difference between their performance increases as the the density of the bipartite graph increases. This confirms that MoReL us can identify potential gene-drug interactions more robustly.

heading: 5.4 COMPARISON WITH BAYREL
text: While the primary goal of experiments so far was showing the effectiveness of MoReL in integrating unstructured and structured views, here we investigate the advantages of MoReL over BayReL.
All structured. As mentioned earlier in the manuscript, BayReL assumes that all of the views are structured. To show the expressive power of MoReL, we train it in the same setting as BayReL where all of the views are structured. Particularly, we assume that for CF dataset both metabolites network and microbiome network are observed in the microbiome-metabolite experiment. Also, in precision medicine experiment both drug network and gene regulatory network are known a priori. For fair
comparison, we set the number of layers as well as hidden dimensions to be the same in both models. We train MoReL with the exponential decaying learning rate with the initial rate of 0.01 and decay rate of 0.001 for 120 training epochs. For BayReL, we use the setting reported in Hajiramezanali et al. (2020). The results for CF and precision medicine are summarized in Tables 3 and 4. We see that MoReL outperforms BayReL on CF dataset by a margin of 7% which indicates that knowing the metabolic pathways can greatly improve interaction learning. In precision medicine experiment, we observe a consistent 2% improvement by MoReL compared to BayReL.
We emphasize that the declined performance of MoReL uu and MoReL us (shown in Tables 1 and 2) compared to BayReL is expected, as they uses less information than BayReL. Incorporating this extra information in MoReL enhances its performance substantially. Note that BayReL is bound to use the same set of functions for all views to account for arbitrarily rotations and transformations, which limits its expressive power. However, the FGW based decoder in MoReL allows to have different processing functions for each view. We argue that this increases the expressive power and plays the key role in enhancing the performance.
Paired vs. unpaired. To show that MoReL can handle unpaired input samples, we perform an ablation study on CF dataset. We reverse the order of samples in metabolite view while keeping the order of samples in microbiome. We report the performance of BayReL and MoReL us where we don’t use the structure of metabolite view. The results are shown in Table 5. While MoReL performed virtually the same as a completely paired scenario (shown in Table 1), BayReL’s performance drastically declined. We note that the reported negative accuracy is the best one achieved by BayReL.
Missing samples. We should again point out that in a setting where all views are structured but the number of node attributes are not the same in different views, BayReL cannot be deployed (as it uses the same processing functions for all views). To see how MoReL us performs in such a scenario, we randomly remove 10% of samples in metabolite view of CF dataset. MoReL achieves positive accuracy
(in %) of 61.36±3.74 with negative accuracy of 97%. This again shows the robustness of FGW-based decoder in aligning nodes with different number of samples.
Computational complexity. We have also benchmarked computational complexity of MoReL and BayReL by tracking their runtime on CF dataset on the same hardware. While BayReL takes 0.6 seconds per training epoch, MoReL takes 2.7 seconds per training epoch. This is due to the computational overhead caused by deploying FGW regularization. Considering the model flexibility and significant prediction performance improvement, such computational overhead is acceptable.

heading: 6 CONCLUSIONS
text: We have proposed MoReL, a novel Bayesian deep generative model that efficiently infers hidden molecular relations across heterogeneous views of data. By using a fused Gromov-Wasserstein (FGW) based decoder, MoReL addresses several main shortcomings of the state-of-the-art omics data integration model. Specifically, MoReL can 1) integrate both structured and unstructured omics datasets while accounting for arbitrarily permutation and/or transformation caused by processing features with different deep functions across the views; 2) handle unpaired samples across the views of data; 3) combine multiple views from different data sources with any number of missing samples. Our experiments on two real-world datasets have demonstrated substantial improvement in inferring meaningful relations as well as improving prediction sensitivity compared to the competing methods. MoReL has shown the promising potential for multi-view learning, in particular multi-omics data integration for biological knowledge discovery, when facing heterogeneous data from different views.

heading: A APPENDIX
text: A.1 FUSED GROMOV-WASSERSTEIN (FGW)
The algorithm to calculate the fused Gromov-Wasserstein distance between two views in our decoder is provided as pseudo-code in Algorithm 1. The algorithm takes C(v), C(v ′), and C(vv ′), which are intra/inter-costs between nodes in a matrix form, as well as β, which is a hyper-parameter. It returns the Wasserstein distance, Gromov-Wasserstein distance, as well as the transport matrix.
Algorithm 1: Computing fused Gromov-Wasserstein distance.
1 Input: C(v)n×n, C (v′) m×m, C (vv′) n×m, β 2 Definitions: = Hadamard product, 〈·, ·〉 = Frobenius dot-product
3 // Cross-view similarity: 4 Ĉ(vv ′) = (C(v))21n1m > + 1n1m >((C(v ′))2)>
5 // Initializing variables:
6 T = 1n1m >, σ = 1
m 1m, Bi,j = exp(Ĉ (vv′) i,j )/β
7 for t1 = 1, 2, . . . do 8 L = Ĉ(vv ′) − 2C(v)T (C(v ′))>
9 for t2 = 1, 2, . . . do 10 M = B T 11 for t3 = 1, 2, . . . do 12 δ = 1
nMσ , σ = 1 nM>δ
13 T = diag(δ)M diag(σ)
14 DW = 〈(C(vv ′))>,T 〉 15 DGW = 〈L>,T 〉
16 Return T , DW , DGW
A.2 DATA DESCRIPTION
Microbiome-metabolome interactions. The goal studying this dataset is to detect the microbemetabolite interactions in patients with Cystic Fibrosis (CF). This dataset includes the 16S ribosomal RNA (rRNA) sequencing and metabolomics for 172 patients diagnosed with CF. We follow the same preprocessing steps as in Morton et al. (2019); Hajiramezanali et al. (2020), and filter out microbes that appear in less than ten samples, which results in 138 unique microbial taxa and 462 metabolite features. To construct the microbiome network, we perform a taxonomic enrichment analysis using Fisher’s test and calculating p-values for each pairs of microbes as in Hajiramezanali et al. (2020). More specifically, the Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995) is adopted for multiple test correction and an edge is added between two microbes if the adjusted p-value is lower than 0.01, The microbiome graph has 984 edges with the graph density of 0.102. For the metabolomics network, there are 1185 edges in total, with each edge representing a connection between metabolites via a same chemical construction (Morton et al., 2019). The graph density of the metabolite network is 0.011. We use 80% of the reported target molecules of P. aeruginosain studies in Quinn et al. (2015) and Morton et al. (2019) as a test set to evaluate the predicted microbiomemetabolome interactions. The remaining 20% of the reported molecules are considered as a validation set and are only used for the early stopping purpose.
Precision medicine. Here we aim to identify genetic markers of cancer drug responses. This is a very challenging task due to the very limited number of observations with respect to the system complexity and huge number of biological and experimental confounders, which often leads to significant false positive associations (Barretina et al., 2012). We consider a dataset from 30 acute myeloid leukemia (AML) patients that contains gene expression and drug sensitivity data of 160 chemotherapy drugs and targeted inhibitors (Lee et al., 2018). For gene expression, we preprocessed the RNA-Seq data resulting in 9073 genes (Lee et al., 2018). Following Hajiramezanali et al. (2020), we construct the gene regulatory network based on the publicly available expression data of the 14 AML cell lines from
the Cancer Cell Line Encyclopedia1 (CCLE) using R package GENIE3 (Vân Anh Huynh-Thu et al., 2010). Moreover, We construct drug-drug interaction networks based on their action mechanisms. Specifically, the selected 53 drugs are categorized into 20 broad pharmacodynamics classes (Lee et al., 2018); 14 classes contain more than one drugs. Only 16 out of the 53 drugs are shared across two classes. We consider that two drugs interact if they belong to the same class. We use the area under the drug response curve reported in the CCLE dataset to indicate drug sensitivity across a range of drug concentrations (Barretina et al., 2012; Lee et al., 2018). Following Lee et al. (2018), we only consider the drugs that have less than 50% cell viability in at least half of the samples, resulting in 53 drugs. We use 797 reported drug-gene interactions in The Drug–Gene Interaction Database (DGIdb) (Wagner et al., 2016) in order to evaluate different models. We note that our test and validation sets only include the interactions for 43 of the 53 drugs in the dataset. We use 20% of the evaluation set as the validation set. Please note that the validation set has been only used for early stopping.
A.3 ADDITIONAL RESULTS FOR CF DATASET
In this section, we provide additional results for CF dataset demonstrating interpretability and robustness of MoReL.
Figure 2 shows two sub-networks of the inferred bipartite relational graphs by MoReLss and MoReLus, consisting P. aeruginosa, anaerobic microbes, and validated target nodes of P. aeruginosa and all of the inferred interactions between them. Based on the biology knowledge, the expected interactions in these sub-networks should be that the four highlighted nodes in the bottom row are connected to all of the nodes in the top row, and any other nodes in the bottom row are not connected to any of the top nodes. At the sub-network negative accuracy of 100% (i.e. any nodes in the bottom row other than the four highlighted ones are not connected to any of the top nodes), while MoReLus identifies 70% of the validated edges of P. aeruginosa, MoReLss identifies 86.8% of the edges. We note that BayReL identifies 78% of the validated interactions (Hajiramezanali et al., 2020). This clearly shows the effectiveness of our proposed FGW-based decoder and interpretability of MoReL to identify inter-relations.
In the main manuscript, we only reported the results for one specific threshold value of negative accuracy (97%). Here we provide additional results with other threshold values, which show similar improvements over competing methods and similar trends by MoReL ss and MoReL us, as clearly observed in Figure 3. We note that there is a trade-off between positive and negative accuracy, and the optimal point can be chosen depending on the application.

Review summaries:
The paper proposed a deep Bayesian model for heterogeneous multi-omics data integration. The Gromov-Wasserstein (FGW) regularization between latent representations of heterogeneous views is used to align nodes/features in every pair of views. The experimental results have demonstrated improvement in inferring meaningful relations.

The authors propose a Bayesian framework to learn relations among multi-omic datasets. The main unique advantage over existing methods is the proposed method is able to learn without a priori dependency structure and it allows a certain degree of missingness and mismatching. Experiments on two biomedical multi-omics data sets partially demonstrate the effectiveness of the proposed method.