Evaluation of O2PLS in Omics data integration
 Said el Bouhaddani^{1}Email author,
 Jeanine HouwingDuistermaat^{1},
 Perttu Salo^{2},
 Markus Perola^{2},
 Geurt Jongbloed^{3} and
 HaeWon Uh^{1}
https://doi.org/10.1186/s128590150854z
© Bouhaddani et al. 2016
Published: 20 January 2016
Abstract
Background
Rapid computational and technological developments made large amounts of omics data available in different biological levels. It is becoming clear that simultaneous data analysis methods are needed for better interpretation and understanding of the underlying systems biology. Different methods have been proposed for this task, among them Partial Least Squares (PLS) related methods. To also deal with orthogonal variation, systematic variation in the data unrelated to one another, we consider the Twoway Orthogonal PLS (O2PLS): an integrative data analysis method which is capable of modeling systematic variation, while providing more parsimonious models aiding interpretation.
Results
A simulation study to assess the performance of O2PLS showed positive results in both low and higher dimensions. More noise (50 % of the data) only affected the systematic part estimates. A data analysis was conducted using data on metabolomics and transcriptomics from a large Finnish cohort (DILGOM). A previous sequential study, using the same data, showed significant correlations between the LipoLeukocyte (LL) module and lipoprotein metabolites. The O2PLS results were in agreement with these findings, identifying almost the same set of covarying variables. Moreover, our integrative approach identified other associative genes and metabolites, while taking into account systematic variation in the data. Including orthogonal components enhanced overall fit, but the orthogonal variation was difficult to interpret.
Conclusions
Simulations showed that the O2PLS estimates were close to the true parameters in both low and higher dimensions. In the presence of more noise (50 %), the orthogonal part estimates could not distinguish well between joint and unique variation. The joint estimates were not systematically affected. Simultaneous analysis with O2PLS on metabolome and transcriptome data showed that the LL module, together with VLDL and HDL metabolites, were important for the metabolomic and transcriptomic relation. This is in agreement with an earlier study. In addition more gene expression and metabolites are identified being important for the joint covariation.
Keywords
Integration of Omics data Dimension reduction Latent variable regression O2PLSBackground
With rapid and continuous technological improvements large amounts of omics data from different levels (genome, transcriptome, proteome and metabolome) are now available. In an integrative systems biology approach, it is becoming increasingly clear that the integration of omics data will provide a better understanding of biological systems. Towards this end, the simultaneous analysis of two data sets is an important task to better understand the relationships between different biological functional levels.
Statistically, integrative approaches face theoretical and computational issues: the typical “large p, small n” problem as in high dimensional data. Some statistical methods require the inverse of matrices; often they are singular, this can be dealt with by penalization or dimension reduction. Interpretation of the results of the analysis is yet another major challenge. In terms of integrating two data sets the following questions need to be answered: (i) which variables in one data set are related to those in another data set, (ii) which variables are not related, but still important, in each of the data sets, and (iii) which variables are relevant, i.e. provide more insight into the biological systems?
A statistical solution is to perform variable selection while combining the two types of variables in the modeled integration process: for example, a regularized version of canonical correlation analysis (CCA) [1], and a variant of partial least squares (PLS) regression [2] called sparse PLS [3] to simultaneously integrate and select variables using lasso penalization [4].
The integration and the variable selection of two different types of omics data sets is nowadays an active research subject. For example, Inouye et al. [5] assessed metabonomic, transcriptomic, and genomic variation for a large populationbased cohort from the capital region of Finland. For an overview of the data integration and the different analyses in the study we refer to Figure 1 of their paper [5]. In this work we focus on the first part of data integration of the paper: ‘metabolite associations of gene modules’. First they identified the sets of highly correlated genes, such as the lipidleukocyte (LL) module, using network analysis of the transcriptomic data. Next a Spearman’s rank correlation was used to identify finescale detail of potentially causative/reactive effects between the LL module expression profile (defined by its first principal component) and the individual metabolites. The motivation of the present paper lies in this sequential analysis procedure. In other areas of biostatistics, simultaneous joint modeling of the variables is known to be more efficient than analyzing data sequentially: network construction, identifying the latent variable or module, and correlating this identified module with the individual metabolites.
Model estimates for integrative parts in the data are often not representing the true underlying biological relation when systematic variation unrelated to the outcome is present, the estimates are biased due to this variation. It has been demonstrated that PLS suffers from this [6]. To deal with this, extensions of PLS have been developed. The asymmetric Orthogonal PLS (OPLS) [7], tries to correct for systematic variation in the design matrix before presenting the data to PLS. The main advantage is an easier interpretation of the model: the model estimates focus more on the predictive variation in the design matrix. In order to integrate two data sets, we need a symmetric approach of OPLS. The Twoway Orthogonal PLS (O2PLS) model [6] is a symmetric method, modeling both predictive and systematic variation. The model decomposes the variation present in two data matrices, for example two omics data matrices X and Y, into three parts. In the first joint part, underlying latent variables in both data matrices are assumed to induce the relationship between X and Y. This joint part can be seen as a representation of the integration of the two data sets. The second part is called the orthogonal part. Underlying latent variables, independent from those in the joint part, are assumed to be responsible for the unique systematic variation in X (Y), which does not contribute to the prediction of Y (X). The third part indicates the noise part, and captures the unsystematic variation in the data.
The aim of this paper is twofold. Our first aim is to jointly model metabolomics and transcriptomics data, in the light of previous study by Inouye et al. [5], to gain a better insight in the interplay between the two omics by decomposing the data in three parts. We extract latent variables for the joint and orthogonal part, and summarize relevant information by looking at the amount of variation captured by these latent variables. Our second aim is to investigate the performance of the O2PLS estimates, in terms of accuracy, with a simulation study under different conditions. We will look at the accuracy in terms of bias, using settings similar to those present in real metabolomics and transcriptomics data.
Integrating metabolomics and transcriptomics using O2PLS is not new. A small scale integration, on 12 aspen grown in a controlled environment, of 453 metabolomic variables and 27,648 transcriptomic data has been performed in [8]. Our analysis is in a larger scale, namely human epidemiological study, consisting of 466 participants. In the metabolomics data set (containing 137 metabolites) we have a classical situation of more participants than variables; the transcriptomics data contains more variables (35419) than participants.
This paper is organized as follows: the “Methods” Section discusses the symmetric integration method O2PLS. A simulation study is set up to assess its performance. In the “Results” Section the simulation results are discussed, furthermore metabolomics and transcriptomics data are analyzed with O2PLS. The “Discussion” Section gives an interpretation of the results from the simulations and data analysis, as well as commenting on the O2PLS model and arguing for a probabilistic approach.
Methods
Previous methods
Here T contains the lower dimensional data. The matrix P contains the directions in the X space which optimizes the covariance T ^{T} Y (where Y has zero mean). The matrix T P ^{T} is to be seen as a ‘best’ approximation of X based on the covariance with Y. The proof for this is deferred to a separate paragraph later on in this section. The matrix E contains the residuals.
Again $\stackrel{~}{T}{\stackrel{~}{P}}^{\mathrm{T}}$ represents a best approximation based on the covariance with Y, but the direction vectors in $\stackrel{~}{P}$ are corrected for (i.e. do not contain directions of) orthogonal variation. The orthogonal variation in X is approximated with ${T}_{\perp}{P}_{\perp}^{\mathrm{T}}$.
Both PLS and OPLS deal with outcome vectors. While generalizations can be made to make them suitable for an outcome matrix, they focus on regressing Y on X, but not simultaneously the other way around. This symmetric approach is appropriate for integrating multiple omics data, while also prediction in both ways can be done.
The O2PLS model
They capture all ‘left over’ variation not captured by the scores.
To approximate Y with X (or X with Y), we need the corresponding inner relation defined via B _{ T } (or B _{ U } ) in (4). A description of the O2PLS algorithm can be found in Trygg’s paper [6]. The inner relation can be recognized as being an ordinary linear model.
 1.
We choose a vector of values for the number of joint components a.
 2.For fixed a we choose the number of orthogonal components n _{ X } and n _{ Y } that maximize the sum of the two coefficients of determination (R ^{2} ) of the inner relation regression (4). Mathematically: we search in a two dimensional grid the integers n _{ X } and n _{ Y } that maximize$({n}_{X},{n}_{Y})\mapsto 1\frac{\sum {\left({H}_{\mathit{\text{UT}}}\right)}_{i,j}^{2}}{\sum \underset{i,j}{\overset{2}{U}}}+1\frac{\sum {\left({H}_{\mathit{\text{TU}}}\right)}_{i,j}^{2}}{\sum \underset{i,j}{\overset{2}{T}}}.$(8)
 3.
Two Mean Squared Errors (MSE) of Prediction concerning $\sum {(\u0176Y)}^{2}$ and $\sum {(\widehat{X}X)}^{2}$  are calculated with 10fold crossvalidation to determine a with the previously obtained n _{ X } and n _{ Y } fixed.
 4.
We go back to step 2 using for a the next element in the vector of values as chosen in step 1.
The quality of the O2PLS estimates depends on the accuracy of the estimated covariance matrix S=X ^{T} Y. Suppose X=E and Y=F, so X and Y are only noise. The covariance matrix S can be decomposed with SVD: S=W D C ^{T}, where W and C are unit norm. It may be that we will observe a ‘large’ positive loading value, since the norm of the loading vectors are forced to be one, and may mistakenly conclude that X and Y are related. However since X and Y are independent the projected data T and U are little correlated (due to noisy variation), thus the inner relation parameters B _{ T } and B _{ U } will have a small magnitude.
Orthogonal correction captures variation unrelated to the joint part. The residual data is hoped to correlate stronger, thus providing a better inner relation fit. Especially with a high number of variables, this may improve the fit (and thus interpretability of the obtained loadings) substantially. Estimation accuracy will not likely be improved by correcting for orthogonal variation, since we do not add information concerning the relation between X and Y. However the exact statistical implications of orthogonality correction on the joint part estimators is still an unclear matter.
Proof of optimality
To make clear why the singular value decomposition is important for O2PLS, some optimality properties are proven.
The maximum of the covariance (9) is attained only if α _{1}=β _{1}=±1. In that case all summands in (12) are zero except when i=1, yielding the maximum to be the first (and largest) singular value. The first singular vectors c=C _{ Y }_{;1} and w=W _{ X }_{;1} are the maximizers. Note that c=−C _{ Y }_{;1} and w=−W _{ X }_{;1} would also yield equivalently the maximum, this is a minor identifiability problem which does not alter the O2PLS model fit. To get the second direction vectors, we optimize the objective function (9) over the unit norm vectors c and w; we require also that c ^{T} C _{ Y }_{;1}=w ^{T} W _{ X }_{;1}=0. This last restriction, the orthogonality constraint, on c and w imply that α _{1}=β _{1}=0 in (12). The maximal covariance is then attained only if α _{2}=β _{2}=1, yielding c and w to be the second singular vectors C _{ Y }_{;2} and W _{ X }_{;2}. Continuing this argument we find the singular vectors in C _{ Y } and W _{ X } to be the maximizers of (9) satisfying the unit norm and orthogonality constraint. If we have a set of indices I for which d _{ i }_{,}i =d _{ j }_{,}j for all i,j∈I, we choose $c={C}_{Y;min\left(I\right)}$ and $w={W}_{X;min\left(I\right)}$ as maximizer. If we have more of those sets, we choose the maximizer in each set in the same fashion.
The maximum is obtained if w _{ Y }_{⊥} is the eigenvector of E ^{T} T T ^{T} E corresponding to the largest eigenvalue. This is the first leftsingular vector of E ^{T} T. Together with the constraint that W _{ Y }_{⊥} should have orthonormal columns, we find W _{ Y }_{⊥} to be the matrix with leftsingular vectors of E ^{T} T. The orthogonal scores can be constructed via T _{ Y }_{⊥}=E W _{ Y }_{⊥}. The same derivation can be used to find that the maximal covariance between U _{ X }_{⊥}:=F P _{ X }_{⊥} and U, where F=Y−U C ^{T}, is obtained if C _{ X }_{⊥} is the collection of leftsingular vectors of F ^{T} U.
Simulation study
A simulation study was performed to investigate the performance of the O2PLS loading estimates, W, C, P _{ Y }_{⊥} and P _{ X }_{⊥}. Although Trygg et al. included a simulation study in their paper [6], the exact simulation study design was not clearly presented. Therefore we could not reproduce their simulation results, and the parameters for our simulation study were arbitrarily chosen.
Simulation parameter choices. The loading value for variable i is the density value of a normal distribution with mean μ and standard deviation σ, denoted as N(i;μ,σ). The noise terms were drawn from a normal distribution with zero mean. The scores were drawn from a standard normal distribution. The variances of the noise terms are such that the expected sum of squares of the noise account for 100α % (equal to 5 or 50 %) of the total sum of squares
Parameter  ‘Low’dimensional case  ‘higher’dimensional case 

N  500  500 
p,q  [100,50]  [500,250] 
W  [N(i;60,10)]_{ i }_{=1,…,100}  [N(i;300,50)]_{ i }_{=1,…,500} 
C  [N(i;70,5)]_{ i }_{=1,…,50}  [N(i;175,25)]_{ i }_{=1,…,250} 
P _{ Y } _{⊥}  [N(i;20,20)]_{ i }_{=1,…,100}  [N(i;100,100)]_{ i }_{=1,…,500} 
P _{ X } _{⊥}  [N(i;15,10)]_{ i }_{=1,…,50}  [N(i;75,50)]_{ i }_{=1,…,250} 
B _{ T }  2  2 
${\sigma}_{T}^{2},{\sigma}_{{T}_{\mathit{\text{Y}}\perp}}^{2},{\sigma}_{{U}_{\mathit{\text{X}}\perp}}^{2}$  [1,1,1]  [1,1,1] 
${\sigma}_{E}^{2},{\sigma}_{F}^{2},{\sigma}_{H}^{2}$  $\frac{\alpha}{(1\alpha )}[0.02,0.104,4]$  $\frac{\alpha}{(1\alpha )}[0.004,0.021,4]$ 
Implementation of the O2PLS algorithm, calculations and analyses were conducted in R [9].
Availability of supporting data
The metabonomic measures are available as Supplementary Table 4 in [5]. The raw and normalized gene expression intensities have been deposited in ArrayExpress which can be found at: http://www.ebi.ac.uk/arrayexpress/under the accession number ETABM1036. ArrayExpress is hosted by the European Bioinformatics Institute.
Results
Results of simulation study
Firstly in both “low”(p=100, q=50) and “higher”(p=500, q=250) dimensions, the accuracy of the estimates were very similar, as can be seen from the location and range of the boxplots. Secondly at the variables with a high joint loading value but low orthogonal loading value, the orthogonal part estimates followed the true orthogonal loading profiles. The joint part estimates also followed the true joint loading profiles regardless of the value of the orthogonal loadings at those variables. Thirdly, the difference between the estimates for the Xand Y components was minor. There was slightly more variation present in the X data at variables with a low loading value.
Application to DILGOM data
Samples on metabolome (137 variables) and transcriptome (35,419 variables) were collected as part of the ‘Dietary, Lifestyle, and Genetic determinants of Obesity and Metabolic syndrome’ (DILGOM) study [5]. Study participants were aged 25–74 years, median age was 53, and were sampled from the region of Helsinki, Finland. A total of 506 participants were present in both studies, of which 232 male and 274 female. In this analysis, we excluded participants whenever they had a missing value for one or more measurements in either the metabolomics or transcriptomics data. This resulted in 40 omitted participants, the used data thus finally consisted of N=466 participants.
The metabolomics data were derived from nuclear magnetic resonance (^{1}H NMR), providing absolute quantitative measurements on the serum metabolome. The transcriptomics data were derived from averaged gene expression counts on technical replicates. The raw counts were quantile normalized at strip level. For more detailed info, see [5], [10]. In transcriptomics filters are proposed to reduce the amount of uninformative (low variance and expression level) variables, which are often interpreted as containing noise. The original study [5] used a filter retaining only the 10 % highest expression levels, and considered 3520 gene expression variables for analysis. To model the orthogonal noise components we were less stringent and extracted the top 25 % of the absolute values of the gene expressions, and we intersected this set of expressions with the set containing the 25 % expressions with the largest interquantile range conform [11]. The reduced transcriptomics data contained 6272 variables. Results of the analysis with all 35,419 variables were very similar (not shown).
A BoxCox transformation [12] with parameter $\frac{1}{4}$ was performed for the metabolomics data, to reduce skewness. The ‘best’ choice for the BoxCox parameter has been investigated by many, we observed from the first four central moments that $\frac{1}{4}$ was sufficient to continue the data analysis. Inouye et al. [5] also applied a BoxCox transformation per variable, but the powers of the transformations were not stated. A scaling here would amplify the effect of noise on the estimates, so the data were only mean centered.
We continued our data analysis with the integration of metabolomics (X) and transcriptomics (Y), using O2PLS. To determine the optimal number of components, we utilized the proposed alternative crossvalidation procedure as discussed in Section “Methods”, initializing with a=1,2,…,10. The optimal number of model components were found a=1, n _{ X } =1, n _{ Y } =8. The modeled variations per component is shown in Table 2. In terms of explained variances (R ^{2}) we observed the following:

The variation in X and Y explained by the model was 58 and 51 % respectively. The rest of the variation was estimated as noise.Table 2
Absolute and relative variations in O2PLS
a
${R}_{X}^{2}$
${R}_{Y}^{2}$
${R}_{\text{Xcorr}}^{2}$
${R}_{\text{Xcorr}}^{2}$
${R}_{\text{Xhat}}^{2}$
${R}_{\text{Yhat}}^{2}$
${R}_{\text{Xhat}}^{2}$/${R}_{\text{Xcorr}}^{2}$
${R}_{\text{Yhat}}^{2}$/${R}_{\text{Ycorr}}^{2}$
1
57.97
50.81
46.31
1.37
26.74
0.80
57.74
58.55
2
67.94
53.40
60.80
4.24
29.52
1.45
48.55
34.25
3
74.08
54.79
68.99
7.35
26.70
2.00
38.69
27.23
4
78.06
55.62
72.94
9.63
29.23
2.40
40.07
24.87
5
80.93
56.69
76.51
11.30
29.81
3.32
38.97
29.43
$1\frac{\sum \underset{i,j}{\overset{2}{E}}}{\sum \underset{i,j}{\overset{2}{X}}}$
$1\frac{\sum \underset{i,j}{\overset{2}{F}}}{\sum \underset{i,j}{\overset{2}{Y}}}$
$\frac{\sum \underset{i,j}{\overset{2}{\left(T{W}^{\mathrm{T}}\right)}}}{\sum \underset{i,j}{\overset{2}{X}}}$
$\frac{\sum \underset{i,j}{\overset{2}{\left(U{C}^{\mathrm{T}}\right)}}}{\sum \underset{i,j}{\overset{2}{Y}}}$
$\frac{\sum \underset{i,j}{\overset{2}{\left(U{B}_{U}{W}^{\mathrm{T}}\right)}}}{\sum \underset{i,j}{\overset{2}{X}}}$
$\frac{\sum \underset{i,j}{\overset{2}{\left(T{B}_{T}{C}^{\mathrm{T}}\right)}}}{\sum \underset{i,j}{\overset{2}{Y}}}$

The joint correlated part in X explained 46 % of the variation in X. Further 1 % of the total variation in Y was explained by the joint correlated part in Y. This means that 46 % of X and 1 % of Y could be explained with one another.

Of the 46 %, Y explained 27 % of X. This could be seen relatively as 57 % of the joint variation in X. Furthermore 0.8 % of Y was explained by X, which was 58 % of the explainable variation in Y.
Absolute and relative variations of the scores and noise in O2PLS
T  T _{ Y } _{⊥}  E  U  U _{ X } _{⊥}  F  H  

Absolute  2551  642  2316  3852  138502  137837  2061 
Relative  46.3 %  11.7 %  42.0 %  1.4 %  49.4 %  49.2 %  53.5 % 
Gene composition of the LL module identified by Inouye et al.
Gene annotation  Ilumina ID 

C1ORF186  ILMN_1690209 
CPA3  ILMN_1766551 
ENPP3  ILMN_1749131 
FCER1A  ILMN_1688423 
GATA2  ILMN_2102670 
HDC  ILMN_1792323 
HS.132563  ILMN_1899034 
MS4A2  ILMN_1806721 
SLC45A3  ILMN_1726114 
SPRYD5  ILMN_1753648 
CACNG6  ILMN_1779043 
LL module and top 10 gene expressions. Identified gene expressions in the top 10 most important variables for the joint variation in the transcriptome. The corresponding genes are shown. Four gene expressions fall into the earlier identifies LipidLeukocyte module
Gene annotation  Ilumina ID  Module 

CPA3  ILMN_1766551  LL and top 10 
FCER1A  ILMN_1688423  LL and top 10 
GATA2  ILMN_2102670  LL and top 10 
HDC  ILMN_1792323  LL and top 10 
DEFA1B  ILMN_1725661  top 10 
DEFA1B  ILMN_1679357  top 10 
DEFA1B  ILMN_2102721  top 10 
SNORD13  ILMN_1892403  top 10 
DEFA3  ILMN_2165289  top 10 
IFIT1  ILMN_1707695  top 10 
Discussion
The integrative systems biology approach is becoming increasingly popular and integration of omics data will provide more insight into the biological systems. The PLS method is widely known in chemometrics and provides data integration and simultaneous modeling, but as shown in [6] the estimates are sensitive to structural noise. While OPLS [7] provides correction for such orthogonal variation, it is oriented towards predicting an outcome and thus lacks symmetry. We considered here the O2PLS method [6]; it is a symmetric data integration method, accounting for structural noise in both matrices. We particularly aimed to integrate two omics data sets for embedding a high dimensional data set in terms low dimensional ‘latent’ variables. To extract relevant information in the data sets, we decompose the two data sets into three parts: joint part in which variables in one data set are related to those in another data set; orthogonal part in which variables are not related, but still important, in each of the data sets; and noise. Simultaneously we searched for the relevant variables in each part.
Several approaches similar to O2PLS are available. To handle more than two data sets, a generalization of O2PLS has been proposed in [13], called OnPLS. Methods to deal with the general idea of decomposing data sets in a joint and systematic part have been proposed. They differ in methodology and estimation. For example, DISCOSCA [14] can handle multiple data sets and may perform better when prior information about the configuration of the joint and orthogonal components is available. An essential assumption in this model is that the components scores or loadings in each data set are exactly the same. Another method providing data decomposition in a joint and orthogonal part is JIVE [15], which can also handle more than two data sets. JIVE may be used if the common source underlying all data sets are similar/homogeneous. One should note that that JIVE restricts the joint part to be orthogonal to the systematic parts. Though it may be argued that the joint and systematic loadings in the population are orthogonal, when obtaining a sample from this population the joint and systematic loadings will typically not be orthogonal. This orthogonality of the joint and systematic loadings is not essential in O2PLS. More research is needed to assess the impact of these methods.
A simulation study is conducted to assess the accuracy of the O2PLS estimates, see Figs. 1, 2, 3 and 4. The estimates were accurate if “little” noise was present (proportion of noise in the data is α=0.05). The model can distinguish well between joint and orthogonal variation. This is the case in both “low”(p=100, q=50) and “higher”(p=500, q=250) dimensional simulated data. The presence of “much” noise (α=0.5) did not cause a substantial decrease in accuracy of the joint part estimates. They followed the true underlying loading profile well. The orthogonal part estimates were affected by more noise in a negative way. Especially in the “higher” dimensional case, the orthogonal part estimates concerning Y (q=250) are biased upwards. The model cannot distinguish well joint and orthogonal variation, it mixes up both loading profiles. It may be argued that the estimation method of the joint loadings is borrowing accuracy from both two data sets, while the orthogonal loadings estimation method is less precise since it uses noisy remaining (total minus joint) variation. Similar to any method, under noisy circumstances it will be difficult to estimate the true orthogonal loadings. This effect was less in the orthogonal part in X (p=500), which has higher dimensions. It is not clear why the orthogonal part estimates with less parameters (the orthogonal part in Y) degrade more than those with more parameters (the orthogonal part in X) in the presence of noise.
We integrate data on the metabolome and transcriptome, extracting both the joint and the orthogonal part, provided in the O2PLS fit, in both data sets. Finding the optimal number of components is a computationally expensive task. A balance between computation time and accuracy is sought by maximizing the explained variance in the inner relation to determine the number of orthogonal parts, and then minimizing the prediction error for determining the number of joint parts. Investing more time in this particular subject will aid in choosing a more accurate method, without compromising computational efficiency. We find four of the eleven LL module gene expressions among the top ten, in terms of importance for the joint variation (Fig. 8). Moreover, the two gene expressions with the highest absolute loading values are in the LL module. Furthermore in the metabolomics data we find the VLDL subgroup together with the HDL subgroup to be important for the joint variation in the metabolomics data (Fig. 7). This shows a contribution of the LL module to the joint variation, partially induced by the VLDL and HDL subgroups. This result can be found back in [5]. The simultaneous data analysis approach identifies more expressed genes important for the joint variation, the ID’s are in Table 5. All genes except SNORD13 are involved in immune/defence system pathways, but information for SNORD13 is at the time of writing unavailable. Also there is large contribution from the mobile lipids MOBCH2 and MOBCH3 to the joint metabolite variation. The orthogonal variation in this data is difficult to interpret, no noticeable trends or clusters were found in the loading values (Figs. 9 and 10). Including orthogonal components in the model does improve the crossvalidated prediction error (which depends on the joint components), which makes it still useful to include in the model. As we saw from the simulation results in the“higher” noise (50 %) case (the estimated amount of noise in the metabolomics and transcriptomics data is also around 50 %), the joint loading estimates still follow the profile of the true loadings. The orthogonal loading estimates are performing worse, indicating a loss of accuracy and thus interpretation in the orthogonal components.
To meet the challenge of interpretation of the results and to infer the relative importance of the variables a structured and tractable probabilistic framework is required. It is beyond the scope of this paper to propose a new method; nevertheless, we argue for the necessity and the feasibility of such a framework. Due to a lack of an explicit probabilistic model in O2PLS, it is not straightforward how to perform statistical tests on the loadings. For PLS, a bootstrap approach is proposed in [16]. In the O2PLS model we must take into account the orthogonal loadings, which are correlated with the joint loadings due to the nature of the estimation algorithm. This may invalidate the bootstrap results. Furthermore a potential problem of multiple testing may exist, which needs to be correctly addressed. The assumptions made in the model imply that the orthogonal scores T _{ Y }_{⊥} and U _{ X }_{⊥} cannot be seen as realisations of random variables, which is a fundamental property in statistical inference. Furthermore without additional assumptions on the orthogonal part loadings P _{ Y }_{⊥} and P _{ X }_{⊥} the model is unidentifiable. Also, the probabilistic approach gives insight in hidden flaws of the estimators, which are very difficult to discover with the current O2PLS algorithm. These potential problems may invalidate statistical inference on the whole population.
Providing a probabilistic framework to nonprobabilistic methods was done earlier. Probabilistic PCA has been developed in [17], and for the factor analysis model there is a well written probabilistic approach in [18]. A novel probabilistic approach for the O2PLS method, which puts the O2PLS method in a statistical framework, is currently being developed. The optimization criterion will be maximum likelihood. The use of a parametric model and a likelihood are indeed restricting the researcher, as one needs to assume a distribution on the data. However we expect that the probabilistic O2PLS model, just as the ordinary linear model, will be robust against small violations of the assumptions. The resulting likelihood can be easily optimized, using a factorization of the probability density which allows for seperately optimizing the likelihood.
A new derivation in multiplatform data analysis we intend to do is the use of a likelihood information score, which will rely on PO2PLS, indicating how much or little two data sets are related. Combining the data integration approach with a probabilistic framework will aid interpretability and inference in more general epidemiological studies.
Declaration
Publication costs for this article were funded by the European Union’s Seventh Framework Programme (FP7HealthF52012) under grant agreement number 305280 (MIMOmics).
This article has been published as part of BMC Bioinformatics Volume 17 Supplement 1, 2016: Selected articles from the Fourteenth Asia Pacific Bioinformatics Conference (APBC 2016). The full contents of the supplements are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/17/S1.
Declarations
Acknowledgements
The research leading to these results has received funding from the European Union’s Seventh Framework Programme (FP7HealthF52012) under grant agreement number 305280 (MIMOmics).
Authors’ Affiliations
References
 González I, Déjean S, Martin PGP, Gonçalves O, Besse P, Baccini A: Highlighting relationships between heterogeneous biological data through graphical displays based on regularized canonical correlation analysis. J Biol Syst. 2009, 17 (02): 17399.View ArticleGoogle Scholar
 Wold H: Estimation of principal components and related models by iterative least squares. Multivariate Analysis (Proc. Internat. Sympos., Dayton, Ohio, 1965). 1966, Academic Press, New YorkGoogle Scholar
 Lê Cao K, Le Gall C: Integration and variable selection of ‘omics’ data sets with pls: a survey. J de la Société Française de Stat. 2011, 152 (2): 7796.Google Scholar
 Tibshirani R: Regression shrinkage and selection via the lasso. J Roy Statist Soc Ser B. 1996, 58 (1): 26788.Google Scholar
 Inouye M, Kettunen J, Soininen P, Silander K, Ripatti S, Kumpula LS, et al.Metabonomic, transcriptomic, and genomic variation of a population cohort. Mol Syst Biol. 2010; 6(1). doi:10.1038/msb.2010.93.Google Scholar
 Trygg J, Wold S: O2pls, a twoblock (x–y) latent variable regression (lvr) method with an integral osc filter. J Chemometr. 2003, 17 (1): 5364.View ArticleGoogle Scholar
 Trygg J, Wold S: Orthogonal projections to latent structures (opls). J Chemometr. 2002, 16 (3): 11928.View ArticleGoogle Scholar
 Bylesjö M, Eriksson D, Kusano M, Moritz T, Trygg J: Data integration in plant biology: the o2pls method for combined modeling of transcript and metabolite data. The Plant J. 2007, 52 (6): 118191.View ArticlePubMedGoogle Scholar
 R Core Team: R: A Language and Environment for Statistical Computing. 2014, R Foundation for Statistical Computing, Vienna, Austria, http://www.Rproject.org/Google Scholar
 Inouye M, Silander K, Hamalainen E, Salomaa V, Harald K, Jousilahti P, et al: An immune response network associated with blood lipid levels. PLoS Genet. 2010, 6 (9): 1001113View ArticleGoogle Scholar
 Liu H, D’Andrade P, FulmerSmentek S, Lorenzi P, Kohn KW, Weinstein JN, et al: mrna and microrna expression profiles of the nci60 integrated with drug activities. Mol Cancer Ther. 2010, 9 (5): 108091.View ArticlePubMedPubMed CentralGoogle Scholar
 Box G, Cox D: An analysis of transformations. J Roy Statist Soc Ser B. 1964, 26: 21152.Google Scholar
 Löfstedt T, Trygg J: Onpls–a novel multiblock method for the modelling of predictive and orthogonal variation. J Chemometr. 2011, 25 (8): 44155.Google Scholar
 Schouteden M, Van Deun K, Wilderjans T, Van Mechelen I: Performing discosca to search for distinctive and common information in linked data. Behav Res Methods. 2014, 46 (2): 57687.PubMedGoogle Scholar
 Lock EF, Hoadley KA, Marron J, Nobel AB: Joint and individual variation explained (jive) for integrated analysis of multiple data types. Ann Appl Stat. 2013, 7 (1): 52310.1214/12AOAS597.View ArticlePubMedPubMed CentralGoogle Scholar
 Wehrens R, van der Linden WE: Bootstrapping principal component regression models. J Chemometr. 1997, 11 (2): 15771.View ArticleGoogle Scholar
 Tipping ME, Bishop CM: Probabilistic principal component analysis. J R Stat Soc Ser B. 1999, 61: 61122. 10.1111/14679868.00196.View ArticleGoogle Scholar
 Rubin D, Thayer D: EM algorithms for ML factor analysis. Psychometrika. 1982, 47 (1): 6976.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. 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.