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 settings, 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) The 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 views 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 privacy/security concerns, leading to unpaired datasets. Apart from these, we might not have access to a priori graph structured 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 lead to degraded performances.
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.
Optimal transport. There have been extensive efforts to utilize Gromov-Wasserstein (GW) discrepancy to solve the alignment problems in shape and object matching (Mé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 & Mémoli, 2019;Vayer et al., 2018b;Xu et al., 2019b). Peyré 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). 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 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;Hajiramezanali et al., 2019;Hasanzadeh et al., 2020) 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.
Wasserstein distance (WD) quantifies the geometric discrepancy between two probability distributions by measuring the minimal amount of "work" needed to move all the mass contained in one distribution onto the other (Solomon et al., 2015). More specifically, given two probability measures Λ ∈ P(X) and ∆ ∈ P(Y), and a transportation cost c : X × Y → R + , WD is the solution to the following
where T i,j is an element of the transport matrix T whose row-wise and column-wise sums equal to [a i ] n i=1 and [b j ] m j=1 , respectively.
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 example, when two distributions are defined in Euclidean spaces with different dimensions or more generally when X and Y are unaligned, i.e. when their features are not in correspondence (Vayer et al., 2019). Instead of measuring inter-domain distances, GWD measures the distance between pairs of samples in one domain and compares it to those in the other domain. More specifically, given two probability measures Λ ∈ P(X) and ∆ ∈ P(Y), as well as two domain-specific transportation costs c (X) : X × X → R + and c (Y) : Y × Y → R + , GWD is the solution to the following optimization problem:
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 ) , π is the transport map, and Π(X×Y) := {π ∈ P(X × Y) | π(x, y) dy = Λ(x), π(x, y) dx = ∆(y)} is the set of all admissible couplings. Likewise, this can be derived for discrete distributions with probability mass functions n i=1 a i δ xi and m j=1 b j δ y j , as follows:
where T i,j is an element of transport matrix T whose row-wise and column-wise sums equal to
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 V s ⊂ V, and the remaining views V u = V \ V s 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 G s = {G (v) } v∈Vs and their adjacency matrices by A s = {A (v) } v∈Vs . We also define X s = {X (v) } v∈Vs as the set of node attributes for structured views, and X u = {X (v) } v∈Vu as the set of data for unstructured views. Moreover, N v denotes the number of nodes in structured views and number of features for unstructured views. MoReL infers the interactions among the nodes in G s and features in X u . We represent these inter-relations by a multi-partite graph with v∈V N v nodes and a multi-adjacency tensor vv ) is the N v × N v bi-adjacency matrix between views v and v .
We define a hierarchical Bayesian model for MoReL with three sets of latent variables: 1)
} v∈Vs∪Vu , which captures the (hidden) structural information; 2) A, which encodes the interaction among features across views; and 3) Z = Z s ∪ Z u = {Z (v) } v∈Vs∪Vu , which summarizes the feature/attribute specific information. The joint probability of observations and latent variables factorizes as follows:
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.
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 variables H, 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 structured probability distributions, Λ ∈ P(X) and ∆ ∈ P(Y), FGW is defined as follows:
where α, β ∈ [0, 1] are scalar hyper-parameters, Π(X × Y) is the set of all admissible couplings between Λ and ∆, and c (XY) , c (X) , and c (Y) are corresponding transportation cost functions. D FGW 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 D FGW . We define the (inter-)cost function for the first term of FGW, i.e. D W , as follows:
where σ denotes the sigmoid function, and
i,: represents the structural latent variable of node/feature i in view v. To calculate the D GW , 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) :
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:
Noting the definitions of WD and GWD in Section 3, we rewrite D FGW between two views of data with shared transport matrixas follows:
j,: ) .
(5)
To approximate the FGW distance, we first deploy GW algorithm in equation ( 1) to obtain T along with the defined transportation cost c (vv ) to calculate Wasserstein distance term in D FGW (Chen et al., 2020). The pseudo-code in Algorithm 1 (Appendix A.1) 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:
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.
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 to map H and A to the parameters of p θz (Z). Specifically,
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:
;
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:
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
Posterior. We model the posterior of the structural latent variables as a Gaussian distribution and infer its parameters independently for each view. More specifically,
We use two GNNs for each structured view, {g
qh (X (v) , A (v) )} v∈Vs , and two fully connected neural networks per unstructured view, {f
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:
We use two GNNs per structured view, {g
qz (X (v) , A (v) )} v∈Vs , and two fully connected neural networks for each unstructured view, {f 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,
While, as mentioned previously, we use the Sinkhorn algorithm to calculate the D FGW , the overall loss is optimized using stochastic gradient descent based optimization algorithms such as Adam (Kingma & Ba, 2014).
Datasets. We use the same datasets as BayReL (Hajiramezanali et al., 2020), i.e. microbiomemetabolite interactions in cystic fibrosis (CF) and gene-drug interactions in precision medicine. Dataset description and graph construction procedure are detailed in Appendix A.2. 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. 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 A 1 and A 2 represent two disjoint sets of metabolites. Negative accuracy is defined as 1 -i∈A 1 j∈A 2 l∈B 1(i and j are connected to l)
, 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.
Baselines. We compare MoReL with three baselines including Spearman's Rank Correlation Analysis (SRCA), BCCA (Klami et al., 2013), and BayReL. 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. 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. In contrast, BayReL is able to use the structure 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 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 (Argelaguet et al., 2018). 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. 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 of H 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 from H 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 D FGW 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.
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 Appendix A.3.
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 in graphs with different densities. This 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, showing that MoReL us can identify gene-drug interactions more robustly.
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 and4. 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 and2) 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. the Cancer Cell Line Encyclopedia1 (CCLE) using R package GENIE3 (Vâ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.
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 MoReL ss and MoReL us , 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 MoReL us identifies 70% of the validated edges of P. aeruginosa, MoReL ss 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.