Skip to main content
  • Methodology Article
  • Open access
  • Published:

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

Abstract

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) [48], 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, 1218]. 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 HumanMethylation450 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.

Fig. 1
figure 1

msPLS identified a combination of 40 epigenomic markers (denoted as ξ1) and 52 transcriptomic markers (denoted as ξ2) that explain the most variance in the proteome variables. The color scale represents the strength of w regression weights

Fig. 2
figure 2

msPLS identified 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

Table 1 The weights of the epigenomic, transcriptomic and proteomic variables extracted by msPLS from the Marfan data

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 transcriptomic biomarkers and the proteomic variables with their corresponding weights can be found in Table 2.

Fig. 3
figure 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

Table 2 The second set of weights of the epigenomic, transcriptomic and proteomic variables extracted by msPLS from the Marfan data

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 p-values 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.

Table 3 Over representation analysis results of the msPLS analysis on Marfan data

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 associated with the classic type of Marfan syndrome. 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 (TGF-beta, associated pathway “TGF-beta receptor signaling in EMT”), is increased in aneurysmal vascular walls [31, 4446]. Finally, we examined the physical interaction and co-expression patterns of the list of all genes identified by the first set of LVs (see Table 1) with the online tool GeneMania (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.

Fig. 4
figure 4

The co-expression pattern of the resulting Marfan genes queried on their biological process based functions. The figure was produced with GeneMania (available at https://genemania.org)

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 state-of-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% 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.

Table 4 The percentage variation in the chronic lymphocytic leukemia (CLL) data sources explained by the subsequent LVs of msPLS and MOFA
Table 5 The weights of the genomic, epigenomic, and transcriptomic variables extracted by msPLS from CLL data sources
Table 6 The loadings of the three subsequent LVs extracted by msPLS from the genomic variables of the CLL data set

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 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” (p-value 1.93×10−4), “TP53 Regulates Transcription of Cell Death Genes” (p-value 8.1×10−4) [47], “HDACs deacetylate histones” (p-value 8.73×10−25) [48], “HS-GAG degradation” (p-value 1.22×10−4), “HS-GAG biosynthesis” (p-value 7.6×10−4), and “Heparan sulfate/heparin (HS-GAG) metabolism” (p-value 5.66×10−3) [49].

Fig. 5
figure 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]

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 be further appraised in the context of known biological pathways, for example via gene set enrichment analysis.

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 explain the most variance in response variables. True positive rates range from 0.61 to 0.99, indicating that the precision of finding truly associated variables improves with increasing sample size. Similarly, true negative rates are reported to quantify msPLS’s ability to exclude irrelevant variables from the final model. True negative rates are above 0.99 for each simulation studies, indicating that the final model excludes irrelevant variables with high precision, regardless of sample size.

The analysis of a genomewide omics dataset of 364,134 epigenomic, 18,424 transcriptomic and 47 proteomic 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 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.

Methods

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 (X1,...,Xq,...,XQ), with the help of latent variables (LVs). Each data source contains pq number of manifest variables (MVs), measured on the same n samples (i.e. \(\phantom {\dot {i}\!}\mathbf X_{q} \in \mathbb {R}^{n \times p_{q}}\)), each data source is assigned to its corresponding LV (ζ1,...,ζq,...,ζQ). The LVs are linear combinations of their MVs (\(\boldsymbol \zeta _{q} = \mathbf X_{q} \mathbf w_{q}\phantom {\dot {i}\!}\), where \(\boldsymbol \zeta _{q} \in \mathbb {R}^{n \times 1}\) and \(\phantom {\dot {i}\!}\mathbf w_{q} \in \mathbb {R}^{p_{q} \times 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;

$$ \boldsymbol \zeta_{q} = \sum_{m=1}^{M_{q}} \theta_{qm} \boldsymbol \zeta_{m \rightarrow q} + \mathbf v_{q}, $$
(1)

where \(\sum _{m=1}^{M_{q}} \boldsymbol \zeta _{m \rightarrow q}\) denotes the sum of Mq LVs that are explanatory for ζq, θqm is the coefficient capturing the effect of the mth ζmq on ζq, and vq 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 wq are estimated as

$$ \mathbf w_{q} = \left[\mathbf X_{q}^{\prime} \mathbf X_{q}\right]^{-1} \mathbf X_{q}^{\prime} \boldsymbol \zeta_{q}, $$
(2)

or as

$$ \mathbf w_{q} = (1/n) \mathbf X_{q}^{\prime} \boldsymbol \zeta_{q}, $$
(3)

depending on the mode of the regression. PLS-PM denotes Eq. (2) as Mode A and Eq. (3) as Mode B regression. For msPLS, Mode A (i.e. multiple univariate regression) is used for the weight vectors of MVs that do not have any response MVs, and Mode B (i.e. multivariate regression) is used for the weight vectors of MVs that do have response MVs. The descriptions of the objective functions of PLS-PM can be found in [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. pq>>n), the covariance matrix of Xq 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

$$ {}\mathbf w_{q} = arg \underset{\mathbf w_{q} }{\text{min}} \: \mathbf w_{q}^{\prime} \left(\frac{\mathbf X_{q}^{\prime} \mathbf X_{q} + \lambda_{2} \mathbf I}{1 + \lambda_{2}} \right) \mathbf w_{q} - 2 \tilde{\boldsymbol \zeta}_{q}^{\prime} \mathbf X_{q} \mathbf w_{q} + \lambda_{1} \mathbf w_{q}, $$
(4)

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 X1, X2, and X3 with p1, p2 and p3 number of variables, measured on n samples (i.e. \(\phantom {\dot {i}\!}\mathbf X_{1} \in \mathbb {R}^{n \times p_{1}}\), \(\phantom {\dot {i}\!}\mathbf X_{2} \in \mathbb {R}^{n \times p_{2}}\) and \(\phantom {\dot {i}\!}\mathbf X_{3} \in \mathbb {R}^{n \times p_{3}}\)), we consider the following relationships between the data sources: X1 and X2 have a symmetric relation (i.e. they are responses for each other). Furthermore, there are asymmetric relations between X1 and X3, and between X2 and X3, such that X3 is response for both X2 and X1 (Fig. 6). These relationships are encoded in a three dimensional connectivity matrix C (i.e. \(\phantom {\dot {i}\!}\mathbf {C} \in \{0,1\}^{3 \times 3}\)), where the entry \(\phantom {\dot {i}\!}c_{{qq}^{\prime }}\) is 1 if data source q is response for data source q, and 0 otherwise (where qq and \(c_{qq^{\prime }}\) indicates the element from qth row and qth column of matrix C). The objective of the analysis is then to simultaneously extract the MVs from X1 and X2 with the highest explanatory power for the variance in MVs of X3.

Fig. 6
figure 6

The proposed relationship between three data sources. X1 and X2 have a symmetric relation (i.e. they are responses for each other) and X3 have asymmetric relation with both X1 and X2 (i.e. X3 is response for both X2 and X1)

Three data sources msPLS algorithm Given data sources X1, X2, and X3, and \(\boldsymbol \Theta = \mathbf C = \left [\begin {array}{lll} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \end {array}\right ]\)

  1. 1.

    Preliminary steps

    1. (a)

      Center and scale X1, X2, and X3

    2. (b)

      Set \(\mathbf w_{1}^{(0)}\), \(\mathbf w_{2}^{(0)} \) and \(\mathbf w_{3}^{(0)}\) initial weight vectors to arbitrary vectors of [1,1,...,1] with length p1, p2 and p3, respectively

    3. (c)

      Define convergence criterion CRT=1 and a small positive tolerance γ=10−6

  2. 2.

    Iterative regression steps

    While CRTγ;

    • a. Estimate initial LVs\(\boldsymbol \zeta _{1} \propto \mathbf X_{1} \mathbf w_{1}^{(0)} \); where indicates that ζ1 is normalised to unit variance

    • \(\boldsymbol \zeta _{2} \propto \mathbf X_{2} \mathbf w_{2}^{(0)} \)

    • \(\boldsymbol \zeta _{3} \propto \mathbf X_{3} \mathbf w_{3}^{(0)} \)

    • b. Model the relationship between data sources

    • (i) Let vector cq be the q-th row of C that indicates the data sources that are explanatory for data source q, i.e. c1=[0,1,0],c2=[1,0,0],c3=[1,1,0]; indicating X1 has one explanatory, X2 has one explanatory and X3 has two explanatory data sources

    • If\({\sum ^{3}_{i=1} c_{qi}} { > 0}\), i.e. if data source q has any explanatory data sources:

    • \(\mathbf \Theta _{\mathbf c_{q} q} = \left [\mathbf Z_{\mathbf c_{q}}^{\prime } \mathbf Z_{\mathbf c_{q}}\right ]^{-1} \mathbf Z_{\mathbf c_{q}}^{\prime } \boldsymbol \zeta _{q},\)

      where Z is the matrix of column bind LVs, i.e. Z=[ζ1,ζ2,ζ3], and \(\mathbf Z_{\mathbf c_{q}}\) is the matrix of column bind explanatory LVs of data source q. Then \(\mathbf \theta _{\mathbf c_{q} q}\) is calculated as follows:

      For c1 we calculate

      \(\mathbf \Theta _{\mathbf {c}_{1} 1} = \theta _{2 1} = \left [\boldsymbol \zeta _{2}^{\prime } \boldsymbol \zeta _{2}\right ]^{-1} \boldsymbol \zeta _{2}^{\prime } \boldsymbol \zeta _{1}\) and the value of θ11 and θ31 remain 0.

      For c2 we calculate

      \(\mathbf \Theta _{\mathbf c_{2} 2} = \theta _{1 2} = [\boldsymbol \zeta _{1}^{\prime } \boldsymbol \zeta _{1}]^{-1} \boldsymbol \zeta _{1}^{\prime } \boldsymbol \zeta _{2}, \) and the value of θ22 and θ32 remain 0,

      and for c3 we calculate

      \(\mathbf \Theta _{\mathbf c_{3} 3} = \begin {array}{l} \theta _{1 3}\\ \theta _{2 3} \end {array}^{\prime } = \big [[\boldsymbol \zeta _{1}, \boldsymbol \zeta _{2}]^{\prime } [\boldsymbol \zeta _{1}, \boldsymbol \zeta _{2}]\big ]^{-1} [\boldsymbol \zeta _{1}, \boldsymbol \zeta _{2}]^{\prime } \boldsymbol \zeta _{3},\)

      where the entries θ13 and θ23 are obtained from the multiple regression step [[ζ1,ζ2][ζ1,ζ2]]−1[ζ1,ζ2]ζ3, and [ζ1,ζ2] is the matrix obtained by column binding ζ1 and ζ2. The value of θ33 remains 0.

    • (ii) Let vector \(\mathbf c_{q^{\prime }}\phantom {\dot {i}\!}\) be the q-th column of C that indicates the data sources that are response for data source q, i.e. c1=[0,1,1],c2=[1,0,1],c3=[0,0,0]; indicating X1 has two responses, X2 has two responses and X3 has no response data sources

    • If\({\sum ^{3}_{i=1} c_{iq^{\prime }}} { > 0}\), i.e. if data source q has any response data sources:

    • \(\boldsymbol \Theta _{\mathbf c_{q^{\prime }} q^{\prime }} = cor(\boldsymbol \zeta _{q^{\prime }},\boldsymbol \zeta _{\mathbf c_{q^{\prime }}}),\)

      i.e., for c1 we calculate

      \(\mathbf {\Theta }_{\mathbf c_1 1} = \begin {array}{l} \theta _{2 1}\\ \theta _{3 1} \end {array} = \begin {array} cor(\boldsymbol \zeta _{1},\boldsymbol \zeta _{2})\\ cor(\boldsymbol \zeta _{1},\boldsymbol \zeta _{3}) \end {array}\)

      and for c2 we calculate

      \(\mathbf \Theta _{\mathbf c_2 2} = \begin {array}{l} \theta _{1 2}\\ \theta _{3 2} \end {array} = \begin {array}{l} cor(\boldsymbol \zeta _{2},\boldsymbol \zeta _{1})\\ cor(\boldsymbol \zeta _{2},\boldsymbol \zeta _{3}) \end {array}\)

      After Steps (b-i) and (b-ii), the entries of Θ are;

      $$\mathbf \Theta =\left[ \begin{array}{ccc} 0 & cor(\boldsymbol{\zeta_{2}},\boldsymbol{\zeta_{1}}) & \theta_{13} \\ cor(\boldsymbol \zeta_{1},\boldsymbol \zeta_{\mathbf 2}) & 0 & \theta_{23} \\ cor(\boldsymbol \zeta_{1},\boldsymbol \zeta_{\mathbf 3}) & cor(\boldsymbol \zeta_{2},\boldsymbol \zeta_{\mathbf 3}) & 0 \end{array}\right],$$

      Notice that θ21 and θ12 in Step (b-i) are overwritten in Step (b-ii). This is because ζ1 and ζ2 are both responses to each other.

    • c. Re-estimate the the latent variables

    • \([\tilde {\boldsymbol \zeta }_{1},\tilde {\boldsymbol \zeta }_{2}, \tilde {\boldsymbol \zeta }_{3}] = [\boldsymbol \zeta _{1},\boldsymbol \zeta _{2}, \boldsymbol \zeta _{3}] \mathbf \Theta \)

    • d. Estimate the new w(1) weights

    • \( \mathbf w_{1}^{(1)} = arg \underset {\mathbf w_{1}^{(0)}}{\text {min}} \: \mathbf w_{1}^{'(0)} \bigg (\frac {\mathbf X_{1}^{\prime } \mathbf X_{1} + \lambda _{2} \mathbf I}{1 + \lambda _{2}}\bigg) \mathbf w_{1}^{(0)} - 2 \tilde {\boldsymbol \zeta }_{1}^{\prime } \mathbf X_{1} \mathbf w_{1}^{(0)} + \lambda _{1} \mathbf w_{1}^{(0)}\)

    • \( \mathbf w_{2}^{(1)} = arg \underset {\mathbf w_{2}^{(0)}}{\text {min}} \: \mathbf w_{2}^{'(0)} \bigg (\frac {\mathbf X_{2}^{\prime } \mathbf X_{2} + \lambda _{2} \mathbf I}{1 + \lambda _{2}}\bigg) \mathbf w_{2}^{(0)} - 2 \tilde {\boldsymbol \zeta }_{2}^{\prime } \mathbf X_{2} \mathbf w_{2}^{(0)} + \lambda _{1} \mathbf w_{2}^{(0)}\)

    • \( \mathbf w_{3}^{(1)} = \big [ [\tilde {\boldsymbol \zeta }_{3}^{\prime } \tilde {\boldsymbol \zeta }_{3}]^{-1} \tilde {\boldsymbol \zeta }_{3}^{\prime } \mathbf X_{3} \big ]^{\prime }\)

    • e. Evaluate the convergence criteria and discard the old w(0) weights

    • \(CRT = \sum ^{3}_{q=1} (\mathbf w_{q}^{(1)}-\mathbf w_{q}^{(0)})^{2}\)

    • \(\mathbf w_{1}^{(0)} = \mathbf w_{1}^{(1)}\), \(\mathbf w_{2}^{(0)} = \mathbf w_{2}^{(1)}\) and \(\mathbf w_{3}^{(0)} = \mathbf w_{3}^{(1)}\)

  3. 3.

    Upon convergence, return\(\mathbf w_{1}^{(0)},\mathbf w_{2}^{(0)}\), and\(\mathbf w_{3}^{(0)}\)

General case

The general case for msPLS can be described as follows. Given Q related data sources X1,...,Xq,...,XQ with p1,...,pq,...pQ corresponding MVs, measured on n samples (i.e. \(\phantom {\dot {i}\!}\mathbf {X}_1 \in \mathbb {R}^{n \times p_1}\),..., \(\phantom {\dot {i}\!}\mathbf {X}_q \in \mathbb {R}^{n \times p_q},...,\mathbf {X}_Q \in \mathbb {R}^{n \times p_Q}\)), and a Q dimensional connectivity matrix C (i.e. C{0,1}Q×Q), where the entry \(\phantom {\dot {i}\!}c_{{qq}^{\prime }}\) 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;

$$ OF \,=\, arg\text{max}\! \left\{\begin{array}{ll} \sum_{r=1}^{R_{q^{\prime}}} Cor (\boldsymbol \zeta_{q^{\prime} \rightarrow r},\boldsymbol \zeta_{q^{\prime}})^{2} & \text{if}\ {\sum^{Q}_{i=1} c_{iq^{\prime}}} { > 0} \\ \sum_{i=1}^{p_{q^{\prime}}} \sum_{m=1}^{M_{q^{\prime}}} Cor (\boldsymbol \zeta_{m \rightarrow q^{\prime}},\mathbf x_{q^{\prime}(i)})^{2} & \text{{otherwise}} \\ \end{array}\right. $$
(5)

where ζq is the LV of \(\phantom {\dot {i}\!}\mathbf X_{q'}\), \(\phantom {\dot {i}\!}\mathbf c_{q'}\) indicates the qth column of matrix C (i.e. \(\phantom {\dot {i}\!}|| \mathbf c_{q'}|| > 0\) indicates that data source q have at least one response data source), \(\mathbf x_{q'(i)}\phantom {\dot {i}\!}\) denotes the ith column of data source \(\phantom {\dot {i}\!}\mathbf X_{q'}\) (i.e. the ith MV of \(\phantom {\dot {i}\!}\mathbf X_{q'}\)), \(\phantom {\dot {i}\!}\sum _{r=1}^{R_{q'}} \boldsymbol \zeta _{q' \rightarrow r}\) denotes the sum of \(\phantom {\dot {i}\!}R_{q'}\) LVs that are response for \(\phantom {\dot {i}\!}\boldsymbol \zeta _{q'}\), and \(\sum _{m=1}^{M_{q'}} \boldsymbol \zeta _{m \rightarrow q'}\phantom {\dot {i}\!}\) denotes the sum of \(\phantom {\dot {i}\!}M_{q'}\) LVs that are explanatory for \(\phantom {\dot {i}\!}\boldsymbol \zeta _{q'}\). In other words, if data source q have at least one response data source, then the squared correlation between \(\phantom {\dot {i}\!}\boldsymbol \zeta _{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 \(\mathbf {X}_{q^{\prime }}\phantom {\dot {i}\!}\) and the combination of the explanatory LVs for \(\phantom {\dot {i}\!}\mathbf {X}_{q^{\prime }}\) is maximised. The symmetric relationship between Xq and \(\mathbf {X}_{q^{\prime }}\phantom {\dot {i}\!}\) is indicated as \(\phantom {\dot {i}\!}c_{{qq}^{\prime }} = c_{q^{\prime } q} = 1\), in which case the OF of their pairwise analysis is to maximise the correlation between their LVs ζq and \(\boldsymbol \zeta _{q^{\prime }}\phantom {\dot {i}\!}\), 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 \(\mathbf {X}_{q^{\prime }}\phantom {\dot {i}\!}\), 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 X1,..,Xq,...,XQ,

  • and Θ=C{0,1}Q×Q, where \(c_{q,q'} = \left \{\begin {array}{ll} 1 & \text { if } \mathbf X_{q} \text { response for } \mathbf X_{q^{\prime }}\\ 0 & {otherwise} \end {array}\right.\)

  1. 1.

    Preliminary steps

    1. (a)

      Center and scale X1,..,Xq,...,XQ

    2. (b)

      Set \(\mathbf w_{q}^{(0)}\) to arbitrary weight vectors [1,1,...,1] with length pq

    3. (c)

      Define convergence criterion CRT=1 and a small positive tolerance γ=10−6

  2. 2.

    Iterative regression steps

    While CRTγ;

    • a. Estimate initial LVs

    • \(\boldsymbol \zeta _{q} \propto \mathbf X_{q} \mathbf w_{q}^{(0)} \); 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 cq be the q-th row of C that indicates the data sources that are exploratory for data source q

    • If\({\sum ^{Q}_{i=1} c_{qi}} { > 0}\), i.e. if data source q has any explanatory data sources:

    • \( \mathbf \Theta _{\mathbf c_{q} q} = [\mathbf Z_{\mathbf c_{q}}^{\prime } \mathbf Z_{\mathbf c_{q}}]^{-1} \mathbf Z_{\mathbf c_{q}}^{\prime } \boldsymbol \zeta _{q},\)

      where Z is the matrix of column bind LVs, i.e. Z=[ζ1,ζ2,ζ3], and \(\mathbf Z_{\mathbf c_{q}}\) is the matrix of the column-bind explanatory LVs of data source q.

    • (ii) Let vector \(\mathbf c_{q'}\phantom {\dot {i}\!}\) be the q-th column of C that indicates the data sources that are response for data source q

    • If\(\phantom {\dot {i}\!}{\sum ^{Q}_{i=1} c_{iq^{\prime }}} { > 0}\), i.e. if data source q has any responses:

    • \(\phantom {\dot {i}\!}\boldsymbol \Theta _{\mathbf c_{q^{\prime }}q^{\prime }} = cor(\boldsymbol \zeta _{q^{\prime }},\boldsymbol \zeta _{\mathbf c_{q^{\prime }}})\)

    • c. Re-estimate the LVs

    • \(\tilde {\mathbf Z} = \mathbf Z \mathbf \Theta \)

    • d. Estimate the new \(\mathbf w_{q}^{(1)}\) weights

    • If Xq doesn’t have any response data sources:

    • \( \mathbf w_{q}^{(1)} = \big [ [\tilde {\boldsymbol \zeta }_{q}^{\prime } \tilde {\boldsymbol \zeta }_{q}]^{-1} \tilde {\boldsymbol \zeta }_{q}^{\prime } \mathbf X_{q} \big ]^{\prime } \)

    • otherwise:

    • \(\mathbf w_{q}^{(1)} = arg \underset {\mathbf w_{q}^{(0)} }{\text {min}} \: \mathbf w_{q}^{'(0)} \bigg (\frac {\mathbf X_{q}^{\prime } \mathbf X_{q} + \lambda _{2} \mathbf I}{1 + \lambda _{2}}\bigg)\mathbf w_{q}^{(0)} - 2 \tilde {\boldsymbol \zeta }_{q}^{\prime } \mathbf X_{q} \mathbf w_{q}^{(0)} + \lambda _{1} \mathbf w_{q}^{(0)}\)

    • e. Evaluate the convergence criteria and discard the old\(\mathbf w_{q}^{(0)}\)weights and calculate OF from Eq. (5) with respect to each data sources

    • \(CRT = \sum (\mathbf w_{q}^{(1)}- \mathbf w_{q}^{(0)})^{2}\)

    • \(\mathbf w_{q}^{(0)} = \mathbf w_{q}^{(1)}\)

  3. 3.

    Upon convergence, return \(\mathbf w_{q}^{(0)}\)

After the algorithm converges, the wq 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 R2 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

$$\begin{aligned} \mathbf{X}_{q(\alpha)}^{res} &= \mathbf{X}_{q(\alpha)} - \hat{\mathbf{X}}_{q(\alpha)}\\ &= \mathbf{X}_{q(\alpha)} - \boldsymbol \zeta_{q(\alpha)} [\boldsymbol \zeta_{q(\alpha)}^{\prime} \boldsymbol \zeta_{q(\alpha)}]^{-1} \boldsymbol \zeta_{q(\alpha)}^{\prime} \mathbf{X}_{q(\alpha)}. \end{aligned} $$

Selecting the optimal penalisation parameters and assessing the statistical significance of the resulting model

In order to obtain the wq 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, X1, X2 and X3, in such way that the relationship between data sources resembles the one we describe in “Chronic lymphocytic leukaemia data” section (Fig. 4).

All Xqs were assigned to pq number of MVs (i.e. p1 = p2 = 1000, p3 = 100) from which kq variables were associated with their LVs and response MVs (i.e. k1 = k2 = k3 = 10), and there were jq number of not associated MVs (i.e. j1 = 990, j2 = 990, j3 = 90). The number of samples are denoted by n samples (i.e. \(\mathbf X_q \in \mathbb {R}^{n \times p_{q}}\) with kq associated MVs and jq not associated MVs, pq = kq + jq), and in the first three simulation studies n varied from 1, 100, and 250.

X1 and X2 were generated from a multivariate normal distribution with mean 0 and covariance matrix Σ, and their response MVs in X3 was generated from LVs ζ1 and ζ2, as follows;

  1. (1)

    Σ=I2000

  2. (2)

    Replace \(\mathbf \Sigma _{1001 : 1010, 1 : 10} = \mathbf \Sigma _{1:10, 1001:1010}^{'} = \mathbf H\),

    where \(\mathbf H \in \mathbb {R}^{10 \times 10}\) distributed over \(\mathcal {N}(0.3,0.05)\)

  3. (3)

    \(\mathbf D \sim \mathcal {N}(0,\mathbf \Sigma)\) where \(\mathbf D \in \mathbb {R}^{n \times 2000}\)

  4. (4)

    X1=D1:n,1:1000 and X2=D1:n,1001:2000

Σ is a p1 + p2 dimensional identity matrix where elements \(\mathbf \Sigma _{1001 : 1010,1 : 10} = \mathbf \Sigma _{1:10, 1001:1010}^{'}\) were replaced with H, where \(\mathbf H \in \mathbb {R}^{10 \times 10}\) was distributed over \(\mathcal {N}(0.3,0.05)\). D was sampled from the multivariate normal distribution with mean 0 and covariance matrix Σ, and D was used to generate X1 and X2. Next, the weight vectors were generated;

  1. (5)

    \(\mathbf w_q = (w_{q(1)}, w_{q(2)},..., w_{q (k_1)}, w_{q (k_{1}+1)},... w_{q (p_1)}\phantom {\dot {i}\!}\)), \(\mathbf w_{q(1:k_1)} = w_{q}^{associated}\phantom {\dot {i}\!}\), \(\mathbf w_{q (k_{q}+1:p_q)} = 0\phantom {\dot {i}\!}\)

The associated k MVs had higher regression weights with their LVs (with weights \(\phantom {\dot {i}\!}w_{1}^{associated}\) = 0.7, \(\phantom {\dot {i}\!}w_{2}^{associated}\) = 0.6, \(w_{3}^{associated}\) = 0.3) than the not associated jq MVs (i.e. \(\mathbf w_q = (w_{q1}, w_{q2},..., w_{q k_q}, w_{q k_{q}+1},... w_{q p_q}\phantom {\dot {i}\!}\)), \(\phantom {\dot {i}\!}\mathbf w_{q (1:k_q)} = w_{q}^{associated}\), \(\mathbf w_{q (k_{q}+1:p_q)} = 0\phantom {\dot {i}\!}\)). The LVs were generated as a linear combination of the MVs and weights,

  1. (6)

    ζ1=X1w1 and ζ2=X2w2

X3 was generated with from ζ1 and ζ2. The k3 associated LVs were sampled from the normal distribution with mean θ1ζ1 + θ2ζ2 (where θq is the regression coefficient from Eq. (1), with θ1 = 0.8 and θ2 = 0.7) and standard deviation \(\sqrt []{1- (w_3)^{2}}\). The j3 not associated variables were sampled from the standard normal distribution;

  1. (7)

    \(\mathbf X_3 \in \mathbb {R}^{n \times 100}\)

  2. (8)

    For i=1,...,k3:

    • X3(i) distributed \(\mathcal {N}(\theta _1 \zeta _1 + \theta _2 \zeta _2, \sqrt []{1- (w_3)^{2}})\)

  3. (9)

    For i=k3+1,...,p3:

    • X3(i) distributed \(\mathcal {N}(0, 1)\)

In addition, we designed a fourth simulation study, where the size of the data resambled the size of the omics data sources, described in “Results” section (i.e. p1 = 360000, p2 = 18000, p3 = 47, k1 = k2 = 40, k3 = 10, and n = 37).

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 kq associated MVs from explanatory data sources X1 and X2, 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 = \sum _{i=1}^{k_q}I(w_{q(i)} \neq 0) / {\text {min}}(k_q, \sum _{i=1}^{p_q}I(w_{q(i)} \neq 0))\)). TPR ranges from 0.61 to 0.99 and increases with increasing sample size when the variable size held constant (Table 7).

Table 7 True-positive rate (TPR) and true-negative rate (TNR) results of the simulation study

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 = \sum _{i=k_{q}+1}^{p_q}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 X3 resulted in 114.21 (95% CI [85.02, 163.43], p-value <0.001), for n=100 the optimised OF resulted in 117.12 (95% CI [71.58, 142.69], p-value <0.001), and for n=250 the optimised OF resulted in 123.61 (95% CI [91.04, 130.94], p-value <0.001).

Fig. 7
figure 7

The null distributions of the optimisation criteria (with respect to X3) 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

Availability of data and materials

An R implementation of msPLS that contains instructions for data simulation and analysis is available at http://uva.csala.me/mspls. The Marfan data analysed during the current study are not publicly available due to privacy reasons, and the Chronic lymphocytic leukaemia data are available in the MOFA R package, https://github.com/bioFAM/MOFA.

Abbreviations

CCA:

Canonical correlation analysis

CCL:

Chronic lymphocytic leukaemia

CI:

Confidence interval

ENeT:

Elastic net

LVs:

Latent variables

LASSO:

Least absolute shrinkage and selection operator

MVs:

Manifest variables

MOFA:

Multi-omics factor analysis

msPLS:

Multiset sparse partial least squares path modeling

OF:

Objective function

OMIM:

Online mendelian inheritance in man

PLS-PM:

Partial least squares path modeling

RDA:

Redundancy analysis

TNR:

True-negative rate

TPR:

True-positive rate

UST:

Univariate soft thresholding

References

  1. Timpson NJ, Greenwood CMT, Soranzo N, Lawson DJ, Richards JB. Genetic architecture: the shape of the genetic contribution to human traits and disease. Nat Rev Genet. 2017; 19(2):110–24. https://0-doi-org.brum.beds.ac.uk/10.1038/nrg.2017.101.

    Article  PubMed  CAS  Google Scholar 

  2. Karczewski KJ, Snyder MP. Integrative omics for health and disease. Nat Rev Genet. 2018; 19(5):299–310. https://0-doi-org.brum.beds.ac.uk/10.1038/nrg.2018.4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Huang S, Chaudhary K, Garmire LX. More Is Better: Recent Progress in Multi-Omics Data Integration Methods. Front Genet. 2017; 8(JUN):1–12. https://0-doi-org.brum.beds.ac.uk/10.3389/fgene.2017.00084.

    Google Scholar 

  4. Tenenhaus A, Tenenhaus M. Regularized generalized canonical correlation analysis for multiblock or multigroup data analysis. Eur J Oper Res. 2014; 238(2):391–403. https://0-doi-org.brum.beds.ac.uk/10.1016/j.ejor.2014.01.008.

    Article  Google Scholar 

  5. Tenenhaus A, Philippe C, Guillemot V, Le Cao K-A, Grill J, Frouin V. Variable selection for generalized canonical correlation analysis. Biostatistics. 2014; 15(3):569–83. https://0-doi-org.brum.beds.ac.uk/10.1093/biostatistics/kxu001.

    Article  PubMed  Google Scholar 

  6. Li W, Zhang S, Liu C-C, Zhou XJ. Identifying multi-layer gene regulatory modules from multi-dimensional genomic data. Bioinformatics. 2012; 28(19):2458–66. https://0-doi-org.brum.beds.ac.uk/10.1093/bioinformatics/bts476.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Karaman I, Norskov NP, Yde CC, Hedemann MS, Bach Knudsen KE, Kohler A. Sparse multi-block PLSR for biomarker discovery when integrating data from LC–MS and NMR metabolomics. Metabolomics. 2015; 11(2):367–379. https://0-doi-org.brum.beds.ac.uk/10.1007/s11306-014-0698-y.

    Article  CAS  Google Scholar 

  8. Hotelling H. Relations Between Two Sets of Variates. Biometrika. 1936; 28(3/4):321. https://0-doi-org.brum.beds.ac.uk/10.2307/2333955.

    Article  Google Scholar 

  9. Csala A, Hof MH, Zwinderman AH. Multiset sparse redundancy analysis for high-dimensional omics data. Biom J. 2018; November 2017:1–18. https://0-doi-org.brum.beds.ac.uk/10.1002/bimj.201700248.

    Google Scholar 

  10. van den Wollenberg AL. Redundancy analysis an alternative for canonical correlation analysis. Psychometrika. 1977; 42(2):207–19. https://0-doi-org.brum.beds.ac.uk/10.1007/BF02294050.

    Article  Google Scholar 

  11. Argelaguet R, Velten B, Arnol D, Dietrich S, Zenz T, Marioni JC, Buettner F, Huber W, Stegle O. Multi-Omics Factor Analysis—a framework for unsupervisedintegration of multi-omics data sets. Mole Syst Biol. 2018; 14(6):8124. https://0-doi-org.brum.beds.ac.uk/10.15252/msb.20178124.

    Google Scholar 

  12. Kim M, Tagkopoulos I. Data integration and predictive modeling methods for multi-omics datasets. Mole omics. 2018; 14(1):8–25. https://0-doi-org.brum.beds.ac.uk/10.1039/c7mo00051k.

    Article  CAS  Google Scholar 

  13. Meng C, Zeleznik OA, Thallinger GG, Kuster B, Gholami AM, Culhane AC. Dimension reduction techniques for the integrative analysis of multi-omics data. Brief Bioinforma. 2016; October 2015:108. https://0-doi-org.brum.beds.ac.uk/10.1093/bib/bbv108.

    Google Scholar 

  14. Li Y, Wu FX, Ngom A. A review on machine learning principles for multi-view biological data integration. Brief Bioinforma. 2018; 19(2):325–40. https://0-doi-org.brum.beds.ac.uk/10.1093/bib/bbw113.

    Google Scholar 

  15. Camacho DM, Collins KM, Powers RK, Costello JC, Collins JJ. Next-Generation Machine Learning for Biological Networks. Cell. 2018; 173(7):1581–92. https://0-doi-org.brum.beds.ac.uk/10.1016/j.cell.2018.05.015. 0608246v3.

    Article  CAS  PubMed  Google Scholar 

  16. Yan J, Risacher SL, Shen L, Saykin AJ. Network approaches to systems biology analysis of complex disease: integrative methods for multi-omics data. Brief Bioinforma. 2017; 19(June 2017):1370–81. https://0-doi-org.brum.beds.ac.uk/10.1093/bib/bbx066.

    Google Scholar 

  17. Min S, Lee B, Yoon S. Deep learning in bioinformatics. Brief Bioinforma. 2016; 18(5):068. https://0-doi-org.brum.beds.ac.uk/10.1093/bib/bbw068.

    Google Scholar 

  18. Ching T, Himmelstein DS, Beaulieu-Jones BK, Kalinin AA, Do BT, Way GP, Ferrero E, Agapow P-M, Zietz M, Hoffman MM, Xie W, Rosen GL, Lengerich BJ, Israeli J, Lanchantin J, Woloszynek S, Carpenter AE, Shrikumar A, Xu J, Cofer EM, Lavender CA, Turaga SC, Alexandari AM, Lu Z, Harris DJ, DeCaprio D, Qi Y, Kundaje A, Peng Y, Wiley LK, Segler MHSS, Boca SM, Swamidass SJ, Huang A, Gitter A, Greene CS. Opportunities and obstacles for deep learning in biology and medicine. J Royal Soc Int. 2018; 15(141):142760. https://0-doi-org.brum.beds.ac.uk/10.1098/rsif.2017.0387. 142760.

    Google Scholar 

  19. Dihazi H, Asif AR, Bei βbarth T, Bohrer R, Feussner K, Feussner I, Jahn O, Lenz C, Majcherczyk A, Schmidt B, Schmitt K, Urlaub H, Valerius O. Integrative omics - from data to biology. Expert Rev Proteom. 2018; 15(6):463–6. https://0-doi-org.brum.beds.ac.uk/10.1080/14789450.2018.1476143.

    Article  CAS  PubMed  Google Scholar 

  20. Zhao Q, Shi X, Huang J, Liu J, Li Y, Ma S. Integrative analysis of ’-omics’ data using penalty functions. Wiley Interdiscip Rev: Comput Stat. 2015; 7(1):99–108. https://0-doi-org.brum.beds.ac.uk/10.1002/wics.1322. NIHMS150003.

    Article  Google Scholar 

  21. Ritchie MD, Holzinger ER, Li R, Pendergrass SA, Kim D. Methods of integrating data to uncover genotype–phenotype interactions. Nat Rev Genet. 2015; 16(2):85–97. https://0-doi-org.brum.beds.ac.uk/10.1038/nrg3868.

    Article  CAS  PubMed  Google Scholar 

  22. Vinzi VE, Trinchera L, Amato S. Pls path modeling: from foundations to recent developments and open issues for model assessment and improvement. In: Handbook of Partial Least Squares. Springer Berlin Heidelberg: 2010. p. 47–82. https://0-doi-org.brum.beds.ac.uk/10.1007/978-3-540-32827-8_3.

    Google Scholar 

  23. Sanchez G. Pls path modeling with r. Berkeley: Trowchez Editions. 2013.

  24. Vinzi VE, Russolillo G. Partial least squares algorithms and methods. Wiley Interdiscip Rev: Comput Stat. 2013; 5(1):1–19. https://0-doi-org.brum.beds.ac.uk/10.1002/wics.1239.

    Article  Google Scholar 

  25. Crick F. Central Dogma of Molecular Biology. Nature. 1970; 227(5258):561–3. https://0-doi-org.brum.beds.ac.uk/10.1038/227561a0.

    Article  CAS  PubMed  Google Scholar 

  26. Shapiro JA. Revisiting the central dogma in the 21st century. Ann NY Acad Sci. 2009; 1178(1):6–28.

    Article  CAS  PubMed  Google Scholar 

  27. Zou H, Hastie T. Regularization and variable selection via the elastic-net. J Royal Stat Soc. 2005; 67:301–20. https://0-doi-org.brum.beds.ac.uk/10.1111/j.1467-9868.2005.00503.x.

    Article  Google Scholar 

  28. Groenink M, Den Hartog AW, Franken R, Radonic T, De Waard V, Timmermans J, Scholte AJ, Van Den Berg MP, Spijkerboer AM, Marquering HA, Zwinderman AH, Mulder BJM. Losartan reduces aortic dilatation rate in adults with Marfan syndrome: A randomized controlled trial. Eur Heart J. 2013; 34(45):3491–500. https://0-doi-org.brum.beds.ac.uk/10.1093/eurheartj/eht334.

    Article  CAS  PubMed  Google Scholar 

  29. Dietrich S, Oleś M, Lu J, Sellner L, Anders S, Velten B, Wu B, Hüllein J, da Silva Liberio M, Walther T, et al. Drug-perturbation-based stratification of blood cancer. J Clin Investig. 2018; 128(1):427–45.

    Article  PubMed  Google Scholar 

  30. Radonic T, de Witte P, Groenink M, de Waard V, Lutter R, van Eijk M, Jansen M, Timmermans J, Kempers M, Scholte AJ, Hilhorst-Hofstee Y, van den Berg MP, van Tintelen JP, Pals G, Baars MJH, Mulder BJM, Zwinderman AH. Inflammation aggravates disease severity in marfan syndrome patients. PLoS ONE. 2012; 7(3):1–9. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pone.0032963.

    Article  CAS  Google Scholar 

  31. Jondeau G, Michel JB, Boileau C. The translational science of Marfan syndrome. Heart. 2011; 97(15):1206–14. https://0-doi-org.brum.beds.ac.uk/10.1136/hrt.2010.212100.

    Article  PubMed  Google Scholar 

  32. Yu E, Foote K, Bennett M. Mitochondrial function in thoracic aortic aneurysms. Cardiovasc Res. 2018; 114(13):1696–8. https://0-doi-org.brum.beds.ac.uk/10.1093/cvr/cvy180.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Ackermann MA, Petrosino JM, Manring HR, Wright P, Shettigar V, Kilic A, Janssen PML, Ziolo MT, Accornero F. TGF- β1 affects cell-cell adhesion in the heart in an NCAM1-dependent mechanism. J Mole Cell Cardiol. 2017; 112:49–57. https://0-doi-org.brum.beds.ac.uk/10.1016/j.yjmcc.2017.08.015.

    Article  CAS  Google Scholar 

  34. Balistreri CR, Ruvolo G, Lio D, Madonna R. Toll-like receptor-4 signaling pathway in aorta aging and diseases: “its double nature”. J Mole Cell Cardiol. 2017; 110:38–53. https://0-doi-org.brum.beds.ac.uk/10.1016/j.yjmcc.2017.06.011.

    Article  CAS  Google Scholar 

  35. Akdis M, Aab A, Altunbulakli C, Azkur K, Costa RA, Crameri R, Duan S, Eiwegger T, Eljaszewicz A, Ferstl R, Frei R, Garbani M, Globinska A, Hess L, Huitema C, Kubo T, Komlosi Z, Konieczna P, Kovacs N, Kucuksezer UC, Meyer N, Morita H, Olzhausen J, O’Mahony L, Pezer M, Prati M, Rebane A, Rhyner C, Rinaldi A, Sokolowska M, Stanic B, Sugita K, Treis A, van de Veen W, Wanke K, Wawrzyniak M, Wawrzyniak P, Wirz OF, Zakzuk JS, Akdis CA. Interleukins (from IL-1 to IL-38), interferons, transforming growth factor β, and TNF- α: Receptors, functions, and roles in diseases. J Allergy Clin Immun. 2016; 138(4):984–1010. https://0-doi-org.brum.beds.ac.uk/10.1016/j.jaci.2016.06.033.

    Article  CAS  PubMed  Google Scholar 

  36. Ju X, Ijaz T, Sun H, LeJeune W, Vargas G, Shilagard T, Recinos A, Milewicz DM, Brasier AR, Tilton RG. IL-6 Regulates Extracellular Matrix Remodeling Associated With Aortic Dilation in a Fibrillin-1 Hypomorphic mgR/mgR Mouse Model of Severe Marfan Syndrome. J Am Heart Assoc. 2014; 3(1):1–13. https://0-doi-org.brum.beds.ac.uk/10.1161/JAHA.113.000476.

    Article  CAS  Google Scholar 

  37. Lenk GM, Tromp G, Weinsheimer S, Gatalica Z, Berguer R, Kuivaniemi H. Whole genome expression profiling reveals a significant role for immune function in human abdominal aortic aneurysms. BMC Genomics. 2007; 8(1):237. https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2164-8-237.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Davis MR, Arner E, Duffy CRE, De Sousa PA, Dahlman I, Arner P, Summers KM. Expression of FBN1 during adipogenesis: Relevance to the lipodystrophy phenotype in Marfan syndrome and related conditions. Mol Genet Metab. 2016; 119(1-2):174–85. https://0-doi-org.brum.beds.ac.uk/10.1016/j.ymgme.2016.06.009.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Syyong H, Chung A, Yang H, van Breemen C. Dysfunction of endothelial and smooth muscle cells in small arteries of a mouse model of Marfan syndrome. British J Pharmacol. 2009; 158(6):1597–608. https://0-doi-org.brum.beds.ac.uk/10.1111/j.1476-5381.2009.00439.x.

    Article  CAS  Google Scholar 

  40. Rayner KJ. Cell Death in the Vessel Wall. Arterioscler Thromb Vasc Biol. 2017; 37(7):75–81. https://0-doi-org.brum.beds.ac.uk/10.1161/ATVBAHA.117.309229.

    Article  CAS  Google Scholar 

  41. Lukashev M. ECM signalling: orchestrating cell behaviour and misbehaviour. Trends Cell Biol. 1998; 8(11):437–41. https://0-doi-org.brum.beds.ac.uk/10.1016/S0962-8924(98)01362-2.

    Article  CAS  PubMed  Google Scholar 

  42. Soto ME, Guarner-Lans V, Herrera-Morales KY, Pérez-Torres I. Participation of Arachidonic Acid Metabolism in the Aortic Aneurysm Formation in Patients with Marfan Syndrome. Front Physiol. 2018; 9(FEB):1–13. https://0-doi-org.brum.beds.ac.uk/10.3389/fphys.2018.00077.

    Google Scholar 

  43. Chung AW, Au Yeung K, Sandor GG, Judge DP, Dietz HC, Van Breemen C. Loss of elastic fiber integrity and reduction of vascular smooth muscle contraction resulting from the upregulated activities of matrix metalloproteinase-2 and-9 in the thoracic aortic aneurysm in marfan syndrome. Circ Res. 2007; 101(5):512–22.

    Article  CAS  PubMed  Google Scholar 

  44. Neptune ER, Frischmeyer PA, Arking DE, Myers L, Bunton TE, Gayraud B, Ramirez F, Sakai LY, Dietz HC. Dysregulation of tgf- β activation contributes to pathogenesis in marfan syndrome. Nat Genet. 2003; 33(3):407.

    Article  CAS  PubMed  Google Scholar 

  45. Bolar N, Van Laer L, Loeys BL. Marfan syndrome: from gene to therapy. Curr Opin Pedia. 2012; 24(4):498–504.

    Article  CAS  Google Scholar 

  46. Judge DP, Dietz HC. Marfan’s syndrome. Lancet. 2005; 366(9501):1965–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Farooqui MZ, Valdez J, Martyr S, Aue G, Saba N, Niemann CU, Herman SE, Tian X, Marti G, Soto S, et al. Ibrutinib for previously untreated and relapsed or refractory chronic lymphocytic leukaemia with tp53 aberrations: a phase 2, single-arm trial. Lancet Oncol. 2015; 16(2):169–76.

    Article  CAS  PubMed  Google Scholar 

  48. Van Damme M, Crompot E, Meuleman N, Mineur P, Bron D, Lagneaux L, Stamatopoulos B. Hdac isoenzyme expression is deregulated in chronic lymphocytic leukemia b-cells and has a complex prognostic significance. Epigenetics. 2012; 7(12):1403–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Sebestyen A, Kovalszky I, Mihalik R, Gallai M, Bocsi J, Laszlo E, Benedek S, Sreter L, Kopper L. Expression of syndecan-1 in human b cell chronic lymphocytic leukaemia. Eur J Canc. 1997; 33(13):2273–7.

    Article  CAS  Google Scholar 

  50. Waaijenborg S, Zwinderman AH. Sparse canonical correlation analysis for identifying, connecting and completing gene-expression networks. BMC Bioinformatics. 2009; 10(1):315. https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2105-10-315.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Csala A, Voorbraak FPJM, Zwinderman AH, Hof MH. Sparse redundancy analysis of high-dimensional genetic and genomic data. Bioinformatics (Oxford, England). 2017; 33(20):3228–34. https://0-doi-org.brum.beds.ac.uk/10.1093/bioinformatics/btx374.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

This work was previously presented at the 40th Annual Conference of the International Society for Clinical Biostatistics. The authors would like to thank the anonymous reviewer and the editor for their valuable comments and suggestions.

Funding

Not applicable.

Author information

Authors and Affiliations

Authors

Contributions

AHZ designed the idea and supervised the study process. ACS analysed the data, implemented the results and wrote the manuscript. AHZ and MHH revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Attila Csala.

Ethics declarations

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1

Appendix: Supplementary materials for the proposed method and results in pdf file format.

Additional file 2

Supplementary materials for the gene set enrichment analysis results of the Marfan data, obtained from Reactome in coma separated file format.

Additional file 3

Supplementary OMIM Gene Map Retrieval results of the Marfan data, obtained from OMIM in xls file format.

Additional file 4

Supplementary co-expression pattern results of the Marfan data, obtained from Gene Mania in txt file format.

Additional file 5

Supplementary gene set enrichment analysis results of the Chronic lymphocytic leukaemia data analyzed by MOFA. Results are obtained by the MOFA R package and exported as txt file format.

Additional file 6

Supplementary gene set enrichment analysis results of the Chronic lymphocytic leukaemia data analyzed by msPLS. Results are obtained by the MOFA R package and exported as txt file format.

Additional file 7

Supplementary material on the overlapping gene set enrichment analysis results between msPLS and MOFA on the Chronic lymphocytic leukaemia data. Results are obtained by the MOFA R package and exported as txt file format.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Csala, A., Zwinderman, A.H. & Hof, M.H. Multiset sparse partial least squares path modeling for high dimensional omics data analysis. BMC Bioinformatics 21, 9 (2020). https://0-doi-org.brum.beds.ac.uk/10.1186/s12859-019-3286-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12859-019-3286-3

Keywords