Multiset sparse partial least squares path modeling for high dimensional omics data analysis

Background Recent technological developments have enabled the measurement of a plethora of biomolecular data from various omics domains, and research is ongoing on statistical methods to leverage these omics data to better model and understand biological pathways and genetic architectures of complex phenotypes. Current reviews report that the simultaneous analysis of multiple (i.e. three or more) high dimensional omics data sources is still challenging and suitable statistical methods are unavailable. Often mentioned challenges are the lack of accounting for the hierarchical structure between omics domains and the difficulty of interpretation of genomewide results. This study is motivated to address these challenges. We propose multiset sparse Partial Least Squares path modeling (msPLS), a generalized penalized form of Partial Least Squares path modeling, for the simultaneous modeling of biological pathways across multiple omics domains. msPLS simultaneously models the effect of multiple molecular markers, from multiple omics domains, on the variation of multiple phenotypic variables, while accounting for the relationships between data sources, and provides sparse results. The sparsity in the model helps to provide interpretable results from analyses of hundreds of thousands of biomolecular variables. Results With simulation studies, we quantified the ability of msPLS to discover associated variables among high dimensional data sources. Furthermore, we analysed high dimensional omics datasets to explore biological pathways associated with Marfan syndrome and with Chronic Lymphocytic Leukaemia. Additionally, we compared the results of msPLS to the results of Multi-Omics Factor Analysis (MOFA), which is an alternative method to analyse this type of data. Conclusions msPLS is an multiset multivariate method for the integrative analysis of multiple high dimensional omics data sources. It accounts for the relationship between multiple high dimensional data sources while it provides interpretable results through its sparse solutions. The biomarkers found by msPLS in the omics datasets can be interpreted in terms of biological pathways associated with the pathophysiology of Marfan syndrome and of Chronic Lymphocytic Leukaemia. Additionally, msPLS outperforms MOFA in terms of variation explained in the chronic lymphocytic leukaemia dataset while it identifies the two most important clinical markers for Chronic Lymphocytic Leukaemia Availability http://uva.csala.me/mspls.https://github.com/acsala/2018_msPLS


Background
Technological developments have enabled the measurement and storage of a plethora of biomolecular data extracted from various omics domains, such as data from the genome, epigenome, proteome or metabolome. It has become common to measure hundreds of thousands of biomolecular variables. To explore biological pathways across multiple omics domains, which might be associated with phenotypic (e.g. disease) outcomes, a natural research direction is to simultaneously analyse these omics domains. Complex diseases, such as obesity, diabetes, and schizophrenia have genetic architectures that involve many biological pathways, since they are a result of interactions between genomic, epigenomic and environmental variables [1,2]. Therefore, modeling biological pathways across multiple omics domains might help to better understand the underlying genetic architecture and biological processes of complex phenotypes, which in turn leads to improved diagnosis, prognosis and therapy [1].
There is ongoing research for suitable statistical methods that could help leverage the available omics data to better model and understand biological pathways and genetic architectures of complex phenotypes on the biomolecular level [3].
Some of the first statistical methods developed for the integrated (i.e. simultaneous) analysis of multiple high dimensional omics datasets are generalizations of well known multivariate methods; e.g. sparse Canonical Correlation Analysis (CCA) [4][5][6][7][8], sparse Redundancy Analysis (RDA) [9,10], and Multi-Omics Factor Analysis (MOFA) [11]. Detailed reviews and discussions on multivariate methods for omics data analysis can be found in [3,[12][13][14][15][16][17][18]. Although there are various statistical methods available to analyse omics data, recent reports argue that the simultaneous analysis of multiple (i.e. three or more) omics data sources is still challenging and current statistical methods are suboptimal. Among the challenges are the lack of accounting for the hierarchical structure between omics domains (i.e. relationship between data sources) and the difficulty of interpretation of genomewide results [2,3,19,20].
To address those challenges, we propose a multiset multivariate statistical method, called multiset sparse Partial Least Squares path modeling (msPLS). msPLS is the penalised extension of multi-block Partial Least Squares path modeling (PLS-PM). Given the situation where biomolecular variables from multiple omics domains are measured on the same patients with shared phenotypes of interest, msPLS models biological pathways by identifying biomarkers (i.e. biomolecular variables that are associated with the phenotypes of interest) in each omics domain. The omics domains are assumed to have a hierarchical structure between each other, and their relationship is modelled in terms of dependencies through explanatory and response domain pairs. The explanatory and response omics data source pairs can be determined through the hypothesised information transfer between data sources as follows [21]. In an asymmetric relationship, a response data source is dependent on a explanatory data source if the prevalent way of information transfer is from the explanatory to the response data source. In a symmetric relationship, there is a recursive information transfer between data sources, and both data sources are dependent on each other. In PLS-PM, latent variables (LVs) are used to model the relationships between explanatory and response manifest variables (MVs) [22,23]. Similarly to PLS-PM, the LVs in msPLS are linear combinations of the MVs, and are estimated in an iterative regression framework [24]. The LVs are constructed so that the combination of the explanatory MVs account for the most variance either directly in the response MVs (in an asymmetric relationship), or in the LVs of the response MVs (in a symmetric relationship). In general, Partial Least Squares path modeling distinguishes between these two types of relationships between data sources (i.e. symmetric or asymmetric relationships) the same way as the two well known multivariate statistical methods Canonical Correlation Analysis [8] and Redundancy Analysis [10] do. In the "Methods" section, we describe msPLS's direct correspondence with those two well known multivariate methods. We give a detailed description of msPLS in the "Methods" section.
To illustrate such an explanatory and response dependency structure, consider that we have biomolecular variables (i.e. genomewide epigenomic, transcriptomic and proteomic variables) measured in patients with Marfan syndrome. The goal of this analysis is to use msPLS to explore biological pathways associated with Marfan syndrome, through the simultaneous analysis of the data sources. For this setting, we assume that the proteomic variables are responses for both the epigenomic and transcriptomic variables. Thus the proteome data source has an asymmetric relationship with both the epigenome and the transcriptome data sources. Additionally, there is a symmetric relationship between the epigenome and the transcriptome data sources, assuming a recursive information transfer between the epigenome and transcriptome. These assumptions are based on the special biological sequential information transfers of the central dogma of molecular biology and its elaborated versions [25,26]. Given the above relationship between omics domains, msPLS identifies the combination of epigenomic and transcriptomic biomarkers that explain the most variance in the proteomic variables, while the combination of the epigenomic and transcriptomic biomarkers have maximum possible correlation with each other. This example is elaborated in more detail in the "Results" section of this paper.
To provide interpretable results from analyses of hundreds of thousands of MVs is addressed through sparse variable selection. msPLS enforces sparse variable selection through penalization methods, such as through the Least Absolute Shrinkage and Selection Operator (LASSO), Ridge, and Elastic Net (ENet) penalization methods [27]. These penalization methods are introduced to PLS-PM by regularising the multivariate regression steps in the iterative regression framework. Introducing regularisation allows msPLS to deal with the characteristic high dimensionality of omics datasets, where the number of variables are much higher than the number of samples. In addition, regularisation improves the interpretability of the final model in the form of sparse variable selection. Once the final model is obtained, the identified biomarkers can be interpreted in terms of biological pathways that are associated with the interest of phenotypes. In the "Methods" section, we quantify msPLS's ability to identify a handful of associated variables from multiple data sources among thousands of irrelevant variables.
The rest of the paper is structured as follows. In the next section, the results of the real data analyses are described, where msPLS was applied to geneomewide biomolecular variables measured in Marfan patients in order to explain the variance in the phenotypic proteomic variables with the combination of biomarkers from the epigenome and transcriptome, while accounting for a hypothesised relationship in omics domains. Additionally, msPLS was applied to a second omics dataset containing data from patients with Chronic lymphocytic leukaemia, and its results were compared to the results of MOFA. We discuss these findings in the "Discussion" section. In the "Methods" section, we describe msPLS and its implementation in an iterative regression framework, along with a working example of the analysis of three related data sources. In addition, we describe how msPLS, and PLS in general, relate to two well known multivariate methods, CCA and RDA. Finally, we show the results from a simulation study that was performed to assess the ability of msPLS to deal with high dimensional data and its ability to extract explanatory MVs that explain the most variance in the response MVs and LVs.

Results
We applied msPLS to genomewide epigenomic, transcriptomic and proteomic data sources measured in Marfan patients [28]. In addition, we applied msPLS to genomic, epigenomic, transcriptomic, and drug response data sources measured in Chronic Lymphocytic Leukaemia (CLL) patients [29].

Marfan data
The goal of this analysis was to explore biological pathways associated with Marfan disease based on epigenomic, transcriptomic and proteomic data measured in 37 Marfan patients [30]. The 364,134 epigenomic methylation variables were obtained by Illumina Infinium Human-Methylation450 BeadChip from blood leukocytes, the 18,424 transcriptomic gene expression variables were obtained by Affymetrix Human Exon 1.0ST Arrays from skin biopsy, and the 47 proteomic cytokine variables were measured in blood plasma.
The model was constructed by extracting the combination of LVs from the epigenome and transcriptome that explain the most variance in the phenotypic proteome MVs (Fig. 1). We hypothesised a symmetric relationship between the epigenome and transcriptome and asymmetric relationships from the proteome to both the epigenome and the transcriptome, so that the proteomic variables were set as response MVs for both the epigenomic and transcriptomic MVs. We used Univariate Soft Thresholding (UST) penalisation with 10-fold cross validation (see "Methods" section) to find the penalisation parameter that optimised the sum of squared correlations between the combination of LVs from the epigenome and transcriptome with respect to the proteome variables (see Eq. (5) in Methods). The final model extracted 40 methylation markers and 52 gene expression markers that optimised the sum of squared correlation of the explanatory LVs of the epigenome and transcriptome with the MVs from the proteome (Fig. 2). The sum of squared correlations was 9.32. Through bootstrapping, we obtained a 95% confidence interval of [9.03, 9.56] and a p-value <0.01 after permutation (see "Methods" section). The best fitting model resulted in a set of LVs that captured 49% of variance in both the epigenome and transcriptome variables and 65% of variance in the proteome variables. The extracted biomarkers with their corresponding individual contribution towards the overall explained variance in the proteomic variables (i.e. illustrated by the methylation and gene expression weights) and the proteomic variables with their corresponding individual correlation strength with the combination of the explanatory LVs (i.e. illustrated by the cytokine weights) are listed in Table 1.
Subsequent set of LVs can be extracted by applying msPLS to the residual data of the epigenome and transcriptome data sources and to the original proteome data source (see "Methods" section). Doing so, we obtained a second set of LVs that explain a different portion of variance in the MVs than the first set of LVs (Fig. 3). After optimizing the model on the residual data, we obtained the second set of LVs that captured 67% of the remaining variance in both the epigenome and transcriptome variables and 91% of the remaining variance in the proteome variables. Thus the first two sets of LVs captured a total of 83% variance in the epigenome and transcriptome variables and a total of 97% variance in the proteome variables. The list of the second set of epigenomic and  Table 2.
A gene set enrichment analysis (available at https:// reactome.org) was used to test the association of the resulting gene expression markers (see Table 1) with already known biological pathways. The gene set enrichment analysis identified 208 pathways (see Additional file 2). We ordered the pathways on their respective pvalues from an over-representation analysis (see https:// reactome.org). For the sake of interpretability, we assessed the pathways with p-values only lower than 5×10 −2 . From the 208 pathways, 58 (28%) had a p-value < 5 × 10 −2 (see Table 3). From these pathways, 44 (76%) can be associated with Marfan disease. From the 58 pathways there are 14 (24%) not known to be associated with Marfan disease, and from these 14 there are 12 pathways that can be associated with the Influenza Virus. This might suggest that Influenza as co-morbidity was present in the patients during data gathering.
Among the pathways that were identified, already known pathophysiological pathways associated with Marfan disease [31] were found, such as the "Extracellular matrix organization" (p-value 4.8 × 10 −3 ), the "Crosslinking of collagen fibrils" (p-value 1.2 × 10 −3 ), the "TGF-beta receptor signaling in EMT (epithelial to mesenchymal transition)" (p-value 3.92 × 10 −2 ), and the "Loss of Function of TGFBR2" (p-value 8.39 × 10 −3 ) pathway. The identified pathways can be further appraised in the context of known interactions of genes and genetic phenotypes. We queried the curated database of Online Mendelian Inheritance in Man (OMIM, available at https://www.omim.org). The OMIM query yielded 372 results (the full list can be found in Additional file 3). Among others, OMIM identified the TGF-beta, Collagen IV, Interleukin-6 loci. The identified pathways from these analysis suggest that some patients suffered from Marfan syndrome type 2, which is based on mutations in the TGFBR2 gene (associated pathway "Loss of Function of TGFBR2"). The mutation in the FBN1 Although MFS2 is phenotypically not separable from classic Marfan syndrome, both disease types include thoracic aortic aneurysm, and more generally aortic risk as the main common feature of the disease [31,43]. This aortic risk is reportedly caused by the loss of function of extracellular matrix proteins (associated pathway "Extracellular matrix organization"), such as collagens and elastin of the vascular wall (associated pathway "Crosslinking of collagen fibrils"), that leads to the loss of solidity and elasticity of the blood vessels, including the aorta, ultimately causing thoracic aortic aneurysm. In addition, it has been reported that the activity of transforming growth factor beta (TGFbeta, associated pathway "TGF-beta receptor signaling in EMT"), is increased in aneurysmal vascular walls [31,[44][45][46]. Finally, we examined the physical interaction and coexpression patterns of the list of all genes identified by the first set of LVs (see Table 1) with the online tool Gen-eMania (available at https://genemania.org). We queried the list of genes based on their biological functions. The analysis resulted in a rich interaction and co-expression pattern (see Fig. 4) with 403 reference studies describing these relationships. The full results of the GeneMania query is available in Additional file 4.

Chronic lymphocytic leukaemia data
We used msPLS for the simultaneous analysis of 69 genomic, 4248 epigenomic, 5000 transcriptomic and 310 drug response variables measured in 200 chronic lymphocytic leukaemia (CLL) patients. This data is publicly available through the Multi-Omics Factor Analysis (MOFA) R package [11]. We used MOFA to impute the missing variables as described in [11]. A detailed description of this dataset can be found in [29]. The goal of this analysis was to compare msPLS performance in terms of explained variance to the performance of MOFA, a stateof-art unsupervised statistical method for the integrative analysis of multiple omics data sources.
To construct the hierarchical structure between the data sources for the msPLS analysis, we hypothesised the following relationship structure between the data source pairs. We assumed symmetric relationships between the genomic, epigenomic and transcriptomic MVs, and the drug response variables were set as response to both the epigenomic and transcriptomic MVs. We used UST penalisation and we compared our results to the results of MOFA. MOFA's model selected 5 non-zero biomolecular variables in each LVs. To compare the results to the msPLS, we enforced the model to extract 5 genomic, 30 epigenomic and 30 transcriptomic MVs from the omics sources. Also, we extracted multiple set of LVs per data source, and compared the total captured variation of msPLS's LVs to MOFA's LVs. The final model of msPLS resulted in 3 set of LVs that together explained 92% of variance in the genomic variables, 97% of variance in the epigenomic variables, 98% of variance in the transcriptomic variables and 85% of variance in the drug response variables. In comparison, MOFA's first 10 LVs (i.e. referred to as factors in the MOFA model) together explained 23% of variance in the genomic variables, 24%

Fig. 3
The resulting model from Section 3.2 extended to two LVs per dataset. The first set of LVs ξ 1(1) and ξ 2(1) partition out a different portion of variance in the proteome MVs than the second set of LVs ξ 1(2) and ξ 2(2) . The colour scale represents the strength of w regression weights of variance in the epigenomic variables, 38% of variance in the transcriptomic variables and 40% of variance in the drug response variables (Table 4). We compared the correlations of Table 5 the selected MVs with their corresponding LVs (i.e. these correlations are referred to as loadings in the MOFA model) from msPLS's and MOFA's models. The biomarkers extracted with msPLS are listed with their corresponding loadings in Table 6.
We also compared the results of MOFA and msPLS in terms of clinical assessment of the outputs of both models (the full clinical assessment of MOFA's results can be found in [11]). For this, we used the gene set enrichment analysis in MOFA's environment. This query resulted in total more than 10,000 pathways, from which 241 pathways with p-values < 0.05 were identified with the gene sets obtained on the CLL data with MOFA, and 298 pathways with p-values < 0.05 were identified with the gene sets obtained on the CLL data with msPLS.
The first 1000 pathways (ordered by their corresponding p-values) for the gene sets from MOFA and msPLS can be found in Additional file 5 and 6. Out of these 1000 pathways, 811 (81%) were identified by both methods, and there are 158 (66% and 53%) overlapping pathways with p-values < 0.05 (see Additional file 7). Similarly to MOFA, msPLS extracted biomarkers from the genomic variables that can be associated with the pathphysiological pathways of CLL. After querying the gene sets from msPLS, the gene set enrichment analysis identified associations with biological pathways such as the "Transcriptional regulation of white adipocyte differentiation" (p-value 3.72 × 10 −4 ), the "Glycerophospholipid biosynthesis" (p-value 5.92 × 10 −5 ), and the "TP53 Regulates Metabolic Genes" (p-value 4.39 × 10 −4 ) pathway in the first LV. In the second LV, the pathways "Keratan sulfate/keratin metabolism" (p-value 5.16 × 10 −5 ), "Post NMDA receptor activation events" (p-value 1.15 × 10 −4 ), and "Activation of NMDA receptor upon glutamate   Arachidonic acid metabolism 4.76 ×10 −2 Arachidonic acid metabolism [42] The pathway names and p-values are obtained from https://reactome.org. Not known associations marked with asterisk (*) are all biomolecular pathways associated with reactions to Influenza virus binding and postsynaptic events" (p-value 2.03 × 10 −4 ) are among the identified ones. Finally, some of the pathways identified in the thirds LV are "Downstream TCR signaling" (p-value 7.27 × 10 −71 ), "Translocation of ZAP-70 to Immunological synapse" (p-value 1.52 × 10 −59 ), "TCR signaling" (p-value 3.14 × 10 −41 ), and "Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell" (p-value 8.68 × 10 −14 ). The two most important clinical markers for CLL, namely the immunoglobulin heavy chain gene (IGHV) and the trisomy of chromosome 12 (trisomy12) were extracted as the first and second LV, respectively (Table 6) [11]. Thus similarly to MOFA, the first two set of LVs from msPLS are aligned among IGHV and trisomy12 (the absolute loading of IGHV is 0.66 in the first LV and the absolute loading of trisomy12 is 0.65 in the second LV), and these can be seen as axis of disease heterogeneity. The samples can be clearly clustered based on their IGHV and trisomy 12 status (Fig. 5). Also, there were 140 pathways with p-values < 0.05 discovered by the gene sets from msPLS that are not overlapping with the pathways discovered by the gene sets from MOFA. Notable pathways that might signal new knowledge discovery are "Regulation of TP53 Activity through Phosphorylation"

Discussion
In this paper, we propose a penalised extension of multiset Partial Least Squares path modeling in response to recent reports pointing out the lack of appropriate statistical methods for the simultaneous analysis of multiple high dimensional omics data sources. msPLS addresses two challenges of integrated high dimensional omics data analysis; namely, it accounts for the relationships between multiple data sources and it provides interpretable results from analyses of hundreds of thousands of biomolecular variables.
Firstly, msPLS accounts for the hierarchical relationship between multiple high dimensional data sources in terms of a explanatory-response dependency structure. It can model dependencies between data sources, such as a hypothesised sequential information transfer in biomolecular domains, through explanatory-response data source pairs. This relationship structure can be easily redefined prior to the analysis, based on most recent biological knowledge. When the relationship is set according to biological knowledge, the biologically relevant biomarkers are identified instead of the variables that explain the most variance in the (combination of ) phenotypic variables. Secondly, msPLS provides interpretable results in the form of combinations of biomarkers that have the highest explanatory power for the variance in the phenotypic variables. The biomarkers are extracted along with their weights that indicate their strength of contribution to the overall explained variance. These biomarkers can Through simulation studies and analyses of omics datasets, we show that msPLS is able to find the combination of biomarkers with the highest explanatory power for the variance in the phenotypic variables, and it can capture a higher proportions of variance in data sources than MOFA, a state-of-art LV based method for multiset omics data analysis. True positive rates of msPLS are reported from the simulation studies (see "Methods" section) to quantify the ability of finding the combination of explanatory variables from the data sources that   variables resulted in biological relevant pathways. msPLS identified a combination of 40 epigenomic biomarkers and 52 transcriptomic biomarkers that has the highest explanatory power for the variance in the phenotypic proteome variables. Despite the low sample size of 37, msPLS identified biomarkers that can be found in known biological pathways associated with the pathophysiology of Marfan disease. Similarly to other LV based multivariate methods, it is possible to extract subsequent LVs with msPLS in a way that they explain a different portion of variance in the data sources. These subsequent LVs are orthogonal Fig. 5 The samples of the CLL data clustered around on their IGHV and trisomy 12 status, extracted by the first and second LV of the msPLS model. The figure was produced by the MOFA R package [11] to each other, thus the newly obtained biomarkers can be interpreted as biological pathways independent from the ones that were discovered in the previous set of LVs. Comparing the results of msPLS and MOFA on the analyses of the CLL dataset, we found that the three set of LVs from the msPLS model captured 92%, 97%, 98% and 85% of the variation in the genomic, epigenomic, transcriptomic and drug response data sources, respectively, while the first ten LVs of MOFA captured a total of 24%, 24%, 38% and 41% of variation in those same data sources, respectively. msPLS, similarly to MOFA, identified the two most important clinical markers for CLL in its first two LVs, and in the "Results" section we additionally report many highly associated and possible novel pathways found through gene enrichment analysis using the MOFA R package.
Note that the present framework of msPLS assumes linear relationships between data sources and that the omics data is measured on a single homogeneous population. As an interesting future direction to extend msPLS is to incorporate non-linear relations in the model or to extend the model such that it can identify different subgroups in the samples.

Conclusions
In summary, msPLS is an appropriate multiset multivariate method that can account for the relationships between high dimensional data sources while it provides interpretable results through its sparse solutions. In the "Methods" section we also describe the algorithm for msPLS and we provide an implementation of the algorithm in the open source R software, which is uploaded with the manuscript and available upon request from the authors. We provide open source code that facilitates the use of our msPLS method on new data with the aim to leverage more and more biomolecular data to model and better understand the genetic architectures and biological processes of complex phenotypes, and ultimately to transition the information synthesised from omics data analyses into medical knowledge to improve diagnosis, prognosis and therapy.

Multiset sparse partial least squares path modeling
Multiset sparse Partial Least Squares path modeling (msPLS) is a multivariate approach for modeling the relationship between Q related data sources (X 1 , ..., X q , ..., X Q ), with the help of latent variables (LVs). Each data source contains p q number of manifest variables (MVs), measured on the same n samples (i.e. X q ∈ R n×p q ), each data source is assigned to its corresponding LV (ζ 1 , ..., ζ q , ..., ζ Q ). The LVs are linear combinations of their MVs (ζ q = X q w q , where ζ q ∈ R n×1 and w q ∈ R p q ×1 ). The relationship between the data sources is encoded in a connectivity matrix, like in Partial Least Squares path modeling (PLS-PM), and modelled through a multiple regression model between the LVs; where M q m=1 ζ m→q denotes the sum of M q LVs that are explanatory for ζ q , θ qm is the coefficient capturing the effect of the mth ζ m→q on ζ q , and v q is white noise, following the notation of [22,24] for PLS-PM. A full description for the PLS-PM algorithm can be found in [24] (Algorithm 6). The weight vectors w q are estimated as or as  [22,24] and the objective function for msPLS is given by Eq. (5) in the "General case" section. In a high dimensional setting (i.e. p q >> n), the covariance matrix of X q in Eq. (2) is non-invertible. To solve this problem, we propose to replace Eq. (2) with Elastic Net (ENet) penalization. Replacing the ordinary least square estimator in Eq. (2) with ENet penalisation has two advantages; not only we overcome the multicollinearity issue encountered in a high dimensional setting, but ENet also enforces sparse variable selection, which ease the interpretability of the final model. Equation (2) then becomes where λ 1 denotes the LASSO penalty and λ 2 denotes the Ridge penalty parameters [27].

An example of msPLS with three data sources
Let us first examine an application of msPLS to three data sources. Given data sources X 1 , X 2 , and X 3 with p 1 , p 2 and p 3 number of variables, measured on n samples (i.e. X 1 ∈ R n×p 1 , X 2 ∈ R n×p 2 and X 3 ∈ R n×p 3 ), we consider the following relationships between the data sources: X 1 and X 2 have a symmetric relation (i.e. they are responses for each other). Furthermore, there are asymmetric relations between X 1 and X 3 , and between X 2 and X 3 , such that X 3 is response for both X 2 and X 1 (Fig. 6). These relationships are encoded in a three dimensional connectivity matrix C (i.e. C ∈ {0, 1} 3×3 ), where the entry c qq is 1 if data source q is response for data source q , and 0 otherwise (where q = q and c qq indicates the element from qth row and q th column of matrix C). The objective of the analysis is then to simultaneously extract the MVs from X 1 and X 2 with the highest explanatory power for the variance in MVs of X 3 .
Three data sources msPLS algorithm Given data sources X 1 , X 2 , and X 3 , and = C = a. Estimate initial LVs Model the relationship between data sources (i) Let vector c q be the q-th row of C that indicates the data sources that are explanatory for data source q, i.e. indicating X 1 has one explanatory, X 2 has one explanatory and X 3 has two explanatory data sources If 3 i=1 c qi > 0, i.e. if data source q has any explanatory data sources: where Z is the matrix of column bind LVs, i.e. Z = ζ 1 , ζ 2 , ζ 3 , and Z c q is the matrix of column bind explanatory LVs of data source q. Then θ c q q is calculated as follows: For c 1 we calculate where the entries θ 13 and θ 23 are obtained from the multiple regression step is the matrix obtained by column binding ζ 1 and ζ 2 . The value of θ 33 remains 0. (ii) Let vector c q be the q -th column of C that indicates the data sources that are response for data source q , i.e. c 1 =[ 0, 1, 1] , c 2 =[ 1, 0, 1] , c 3 =[ 0, 0, 0] ; indicating X 1 has two responses, X 2 has two responses and X 3 has no response data sources If 3 i=1 c iq > 0, i.e. if data source q has any response data sources: Fig. 6 The proposed relationship between three data sources. X 1 and X 2 have a symmetric relation (i.e. they are responses for each other) and X 3 have asymmetric relation with both X 1 and X 2 (i.e. X 3 is response for both X 2 and X 1 ) c q q = cor(ζ q , ζ c q ), i.e., for c 1 we calculate cor(ζ 1 , ζ 3 ) and for c 2 we calculate and (b-ii), the entries of are; Notice that θ 21 and θ 12 in Step (b-i) are overwritten in Step (b-ii). This is because ζ 1 Evaluate the convergence criteria and discard the old w (0) weights CRT = 3 q=1 (w (1) 2 , and w (0) 3

General case
The general case for msPLS can be described as follows. Given Q related data sources X 1 , ..., X q , ..., X Q with p 1 , ..., p q , ...p Q corresponding MVs, measured on n samples (i.e. X 1 ∈ R n×p 1 ,..., X q ∈ R n×p q , ..., X Q ∈ R n×p Q ), and a Q dimensional connectivity matrix C (i.e. C ∈ {0, 1} Q×Q ), where the entry c qq is 1 if data source q is a response data source for data source q and 0 otherwise. The goal of the analysis then is to optimise the following objective function () in respect to data source q ; where ζ q is the LV of X q , c q indicates the q th column of matrix C (i.e. ||c q || > 0 indicates that data source q have at least one response data source), x q (i) denotes the ith column of data source X q (i.e. the ith MV of X q ), R q r=1 ζ q →r denotes the sum of R q LVs that are response for ζ q , and M q m=1 ζ m→q denotes the sum of M q LVs that are explanatory for ζ q . In other words, if data source q have at least one response data source, then the squared correlation between ζ q and the combination of its response LVs is maximised, and if data source q does not have any response data sources, the correlation between the MVs of X q and the combination of the explanatory LVs for X q is maximised. The symmetric relationship between X q and X q is indicated as c qq = c q q = 1, in which case the OF of their pairwise analysis is to maximise the correlation between their LVs ζ q and ζ q , corresponding to the characteristic objective function of Canonical Correlation Analysis (CCA) [8,22,50]. In an asymmetric relationship, the OF of a pairwise analysis is to maximise the sum of squared correlation between the explanatory LV ζ q and the response MVs in X q , corresponding to the characteristic objective function of Redundancy Analysis (RDA) [10,22,51]. This direct correspondence with CCA and RDA is described in Additional file 1 under the Modes of relationships between data sources section.
Next we describe the general algorithm for Q data sources.

General msPLS algorithm
Given Q data sources X 1 , .., X q , ..., X Q , and = C ∈ {0, 1} Q×Q , where c q,q = 1 if X q response for X q 0 otherwise a. Estimate initial LVs ζ q ∝ X q w (0) q ; where q is the index from 1 to Q and ∝ indicates that ζ q is normalised to unit variance b. Model the relationship between data sources (i) Let vector c q be the q-th row of C that indicates the data sources that are exploratory for data source q If Q i=1 c qi > 0, i.e. if data source q has any explanatory data sources: where Z is the matrix of column bind LVs, i.e. Z =[ ζ 1 , ζ 2 , ζ 3 ], and Z c q is the matrix of the column-bind explanatory LVs of data source q.
(ii) Let vector c q be the q -th column of C that indicates the data sources that are response for data source q If Q i=1 c iq > 0, i.e. if data source q has any responses: c. Re-estimate the LVs Z = Z d. Estimate the new w (1) q weights If X q doesn't have any response data sources: Evaluate the convergence criteria and discard the old w (0) q weights and calculate OF from Eq. ( 5) with respect to each data sources After the algorithm converges, the w q weights indicate the contribution of explanatory MVs from the qth data source towards the overall explained variance in the response MVs or LVs (see Additional file 1 under Modes of relationships between data sources section). Through the penalisation of the multivariate regression in Step (2-d), a small subset of explanatory MVs are extracted, namely those with the highest explanatory power for the variance in their response MVs or LVs. The extracted set of MVs can be further explored in terms of known biological pathways, for example through gene enrichment analysis.

Multiple latent variables per dataset
It is possible to extract multiple LVs per data source in a way that they explain a different portion of variance in the MVs. The explained variance is based on the R 2 statistic obtained from the regression model from Step (2-d) in the general msPLS algorithm. The subsequent latent variables can be obtained by applying msPLS to the residual data sources, where the residuals data sources are calculated as

Selecting the optimal penalisation parameters and assessing the statistical significance of the resulting model
In order to obtain the w q weights that optimise OF in Eq. (5), the optimal LASSO and Ridge penalisation parameters can be selected through k-fold cross validation. Given the usual size of omics data and the multiset approach of the analysis, searching for the optimal penalisation parameters is often too computationally expensive. As a solution, we propose to use Univariate Soft Thresholding (UST), by setting λ 2 → ∞ in Eq. (4) [27].
To assess the statistical significance of a resulting model in respect to the OF in Eq. (5), we use a standard permutation approach. The null distribution of the optimisation criterion is estimated by applying msPLS to permuted dataset, where we permute the rows of each dataset. Permuting the samples removes the correlation between data sources while the internal correlation structure of each data source is preserved. The weights obtained from the permutation are used to calculate OF, and the null distribution of the optimisation criterion can be approximated by repeating the permutation a large number of times. In addition, we use bootstrapping to approximate the confidence intervals for the optimised OF. During bootstrapping, the observations are sampled with replacement and the penalisation parameters from the original model are used for the bootstrap samples. In contrast to permutation, with bootstrapping the correlation between data sources is also preserved. After repeating the bootstrapping many times, the selected quantiles of the resulting distribution are reported.

Assessing msPLS's ability to identify associated variables among multiple high dimensional data sources
Before we applied msPLS to omics data sources, we analysed simulated data to assess msPLS's ability to extract the associated MVs from multiple high dimensional data sources that optimise the OF in Eq. (5). Then we applied msPLS to omics data sources to see whether the resulting model can be interpreted in terms of known biological pathways. Below, we describe the simulation studies, and the real data analysis can be found in the "Results" section.

Simulation studies
We conducted simulation studies to assess msPLS's ability to identify associated MVs (i.e. explanatory MVs that are highly correlated with their response MVs and thus have the highest explanatory power for the variance in the response MVs) when those MVs are spread over multiple data sources. We repeated the simulations 1000 times and used UST penalisation for which the optimal penalty parameter (λ 1 ) was selected through 10-fold cross validation. Additionally, we assessed the statistical significance of the resulting models through permutations and the confidence interval of the optimisation criterion was approximated through bootstrapping.

Data generation for simulation studies
For all simulation studies, we generated three data sources, X 1 , X 2 and X 3 , in such way that the relationship between data sources resembles the one we describe in "Chronic lymphocytic leukaemia data" section (Fig. 4).

Simulation study results
We generated data as described in above with three different sample sizes, i.e. n = 50, n = 100, n = 250. To assess msPLS's ability to identify the k q associated MVs from explanatory data sources X 1 and X 2 , we used the true-positive rate (TPR) and true-negative rate (TNR) measures over 1000 simulations.
TPR measures the proportion of associated MVs included in the final model (i.e. those that are assigned to non-zero w weights) to either the number of associated MVs that were generated, or to the total number of non-zero w weights, whichever is smaller (i.e. TPR q = k q i=1 I(w q(i) = 0)/min(k q , p q i=1 I(w q(i) = 0))). TPR ranges from 0.61 to 0.99 and increases with increasing sample size when the variable size held constant (Table 7).
TNR measures the proportion of not associated MVs excluded from the model to the number of not associated MVs that were generated (i.e. TNR q = p q i=k q +1 I(w q(i) = 0)/j q ). TNR rates resulted in 0.99 and were not affected by the sample size (Table 7).
We assessed the statistical significance of the resulting models in respect to the optimised OFs through permutation, and the confidence intervals of the optimised OFs were constructed through bootstrapping (see the "Selecting the optimal penalisation parameters and assessing the statistical significance of the resulting model" section). All the three models obtained on the three different sample sizes with constant variable size were statistically significant, and the confidence interval of the optimised OFs shrank with increased sample size (Fig. 7). For n = 50, the optimised OF with respect to X 3 . 7 The null distributions of the optimisation criteria (with respect to X 3 ) for the simulated data with different sample sizes (n = 50, 100, 250), obtained after 1000 permutations. The red bars indicate the optimisation criteria obtained applying msPLS to the original data with the optimal λ 1 parameters for UST. The red dots are the bootstrapped values, and the dashed red bars are the 95% confidence intervals

Availability and Requirements
Project name: msPLS implementation Project home page: http://uva.csala.me/mspls and https:// github.com/acsala/2018_msPLS Operating system(s): Platform independent Programming language: R Other requirements: additional R packages listed in the source code, freely available from the Comprehensive R Archive Network online License: MIT Any restrictions to use by non-academics: MIT license applies