Abstract
Background
High throughput ’omics’ experiments are usually designed to compare changes observed between different conditions (or interventions) and to identify biomarkers capable of characterizing each condition. We consider the complex structure of repeated measurements from different assays where different conditions are applied on the same subjects.
Results
We propose a twostep analysis combining a multilevel approach and a multivariate approach to reveal separately the effects of conditions within subjects from the biological variation between subjects. The approach is extended to twofactor designs and to the integration of two matched data sets. It allows internal variable selection to highlight genes able to discriminate the net condition effect within subjects. A simulation study was performed to demonstrate the good performance of the multilevel multivariate approach compared to a classical multivariate method. The multilevel multivariate approach outperformed the classical multivariate approach with respect to the classification error rate and the selection of relevant genes. The approach was applied to an HIVvaccine trial evaluating the response with gene expression and cytokine secretion. The discriminant multilevel analysis selected a relevant subset of genes while the integrative multilevel analysis highlighted clusters of genes and cytokines that were highly correlated across the samples.
Conclusions
Our combined multilevel multivariate approach may help in finding signatures of vaccine effect and allows for a better understanding of immunological mechanisms activated by the intervention. The integrative analysis revealed clusters of genes, that were associated with cytokine secretion. These clusters can be seen as gene signatures to predict future cytokine response. The approach is implemented in the R package mixOmics (http://cran.rproject.org/ webcite) with associated tutorials to perform the analysis^{a}.
Background
Recent advances in high throughput ‘omics’ technologies enable quantitative measurements of expression or abundance of biological molecules in a whole biological system. Various popular omics platforms in systems biology include transcriptomics, proteomics, cytomics and metabolomics. These experiments are usually designed to compare changes observed between different conditions or groups and are often used to identify biomarkers capable of characterising pathological states or response to treatment.
The decreasing costs of these highthroughput platforms now enable repeated measures experiments on the same individuals or biological samples. Such experiments allow a substantial gain in information. For instance, longitudinal designs are more powerful as they reduce the noise due to inter individual variability, as long as the correlation between repeated observations is taken into account. There exists an abundant literature on the analysis of repeated measurements of omics data [1,2]. In this context, a common approach is to apply a univariate mixed model on each gene followed by multiple testing correction [3]. However, this approach disregards the dependency between genes, and due to the high dimensionality of the data, numerous hypotheses tests must be performed.
The mixed model approach has been used for the analysis of one single data type (e.g. gene expression). However, a growing number of highthroughput data are generated in standard clinical trials. For example, the evaluation of HIV vaccine in phase I/II trials incorporates measurements of counts of numerous types of cell, of the production of intra and extracellular cytokines and of gene expression [4]. The integration of such multilayer information can help unravel the complexities of a biological system, as each functional level is hypothesized to be related to each other [5]. However the integration of omics data is a challenging task. Firstly, the large number of measured biological entities makes it very difficult to obtain a good overview or understanding of the system under study. Secondly, the small number of samples or patients makes statistical inference difficult and argue for using the maximum amount of available information. Thirdly, the integration of heterogeneous data represents an analytical and numerical challenge when trying to find common patterns in data from different origins.
In recent years, several multivariate approaches have been proposed to combine two omics data, often in an unsupervised framework. In contrast to univariate repeated measures analysis, these linear multivariate approaches take into account the dependency between genes, are able to handle large and noisy data sets and do not face computational issues in the high dimensional case as matrix inversions are avoided. Most importantly in the context of this study, they enable the integration of data coming from different platforms and provide interpretable visualisation tools. These approaches aim at selecting correlated biological entities from two [611] or more data sets [12]. In particular, with sparse Partial Least Squares (sPLS) we have shown that the integrative analysis of large scale omics datasets could generate new knowledge not accessible by the analysis of a single data type alone [7,8]. The biological relevance of this approach has been illustrated recently in some studies [13,14].
The flexibility and versatility of PLS also enable a supervised framework through PLSDiscriminant Analysis (PLSDA [15]). A variant of which has recently been proposed to select discriminative features that best separate the different conditions (sPLSDA, [16]). sPLSDA was shown to give similar performances to classical classification methods such as Machine Learning approaches and variants of Linear Discriminant Analysis and was recently applied in a biological study [17].
In this paper, we consider a twostep approach to model the correlation between repeated measurements while taking advantage of the multivariate approaches. We first propose to extract the withinsample variation [1820] before analysing this transformed data set using sPLSDA for a discriminant analysis or sPLS for an integrative analysis.
Starting from the classical mixedmodel, we present the principle of a multilevel analysis to extract the withinsample deviation of the data and we extend the approach to a twofactor analysis. The within data set is then analysed with either sPLSDA to select discriminative genes between the groups of subjects on a single data set, or with sPLS to select subsets of correlated variables from two data sets. A simulation study is performed which demonstrates the good performance of multilevel sPLSDA compared to a classical sPLSDA. The approach is then illustrated on an HIV vaccination study, where the effect of a lipopeptide based vaccine was explored by measuring before and after vaccination various components of the immune response, including gene expression and cytokine secretion. These repeated measurement were made in several in vitro conditions on Peripheral Blood Mononuclear Cells: ‘NS’ (no stimulation); HIV Gag peptides ‘GAG+’ (peptides included in the vaccine), HIV Gag peptides ‘GAG’ (peptides not included in the vaccine) and ‘LIPO5’ (all five peptides included in the vaccine).
Methods
Notations
Let X(N × p) and Z(N × q) represent two data matrices (e.g. gene expression and cytokine secretion). We denote by N the total number of samples (or rows) in the data, and by n the number of experimental units (or unique subjects), p (resp. q) is the total number of genes (resp. cytokines), also called variables or predictors. The dummy matrix Y(N × G) indicates the group/treatment of the samples, with G the total number of groups.
Multilevel approach
We first present the mixedeffect model as a pedagogical tool and then introduce the concept of the multilevel approach based on the “splitup” variation. Despite the fact that some similarities exist between the mixedeffect model and the multilevel “splitup” variation approach, we emphasize that the latter is performed completely independently from the estimation of the mixedeffect model. Moreover, the mixed model relies on certain assumptions (such as Gaussian distribution of random effects) that the splitup variation approach does not require.
Mixedeffect model
Let be the gene expression of a given gene k for subject s with stimulation j. In this context, the mixed model is defined by:
where for a given gene k, measures the fixed effect of stimulation j, which can be further decomposed into , the overall mean stimulation effect, plus which is the differential effect for stimulation j. The are independent random variables following a normal distribution , which take into account the dependency between the repeated measures made on the same subject s, the residuals are independent random variables following a distribution. Note that the number of stimulation for each subject (G_{s}) may differ. However, for each subject s we observe no more than one observation for each stimulation j, thus the subject effect interactions with the stimulation factor are confounded with the residuals. We also assume that and are independent.
This model is known also as the oneway unbalanced randomeffects ANOVA. A simple approach for identifying differentially expressed (DE) genes in this model is to test the stimulation effect for each gene and apply a multiple testing correction (FDR from [21] set to 5%). Following this global test, pairwise comparison can then be applied between two stimulations, followed by multiple correction (e.g. FDR, 5%). Some limitations of this standard approach are discussed in the Results and discussion Section. The main advantage of the mixed model is the introduction of the random element , which is specific to the subject s and represents the betweensubject deviation. In the same spirit, the multilevel approach based on the splitup variation focusses on separating the different sources of variation: the withinsubject deviation (“variation”) and the betweensubject deviation.
Splitup variation
As suggested by Westerhuis et al. [19] in the mixed model framework, the observation can be decomposed into:
where and . The offset term is an estimation of , the betweensubject deviation is an estimation of and the withinsubject deviation is an estimation of and can be further decomposed as:
where , with n_{j} the number of subject undergoing stimulation j. Therefore, a part of the withinsubject deviation is explained by the stimulation effect.
Let X be the (N× p) gene expression matrix on s=1,…,n subjects with G_{s} stimulations (in the balanced case N=n× G, otherwise ). According to equation (1):
The matrix X_{··} represents the offset term defined as , where 1_{N} is the (N× 1) matrix containing ones and ; X_{b} is the betweensubject matrix of size (N× p) defined by concatenating for each subject into X_{b} with ; X_{w} =X−X_{s·} is the withinsubject matrix of size (N× p), with X_{s·} the matrix defined by concatenating the matrices for each subject into X_{s·}, with .
Similarly to the Analysis of Variance, it is easy to show that the sum of squares can be separated into three parts:
where X^{2}=trace(X^{T}X). Equation (2) can be used to evaluate the magnitude of the different sources of variation.
The mixedmodel described earlier can provide an analysis for repeated measurements data in an unbalanced design. It can be viewed as an extension of a paired ttest to test the differences between paired observations. However, to tackle some of the previously mentioned limitations of the approach, we propose to combine a multilevel approach and a multivariate approach as an interesting alternative. Indeed, the multilevel step splits the different parts of the variation while taking into account the repeated measurements on each subject. Since the stimulation effect from each subject can be separated from the between subject deviation (variation), it is possible to examine the differences in stimulation effect within the subjects in a much easier way than without the separation of the difference sources of variation [19]. Westerhuis et. al (2010) provided the rationale and showed the benefit of the multilevel approach in the analysis of multivariate paired (crossover) data. In this paper where we aim at identifying genes discriminating the different stimulations, we propose to apply a multivariate approach on the within matrix X_{w} which includes the stimulation effect, in the same spirit as in [1820]. This approach is more powerful, as it takes into account not only the dependency between genes via the multivariate approach, but also the repeated measures between individuals and the stimulation effects via X_{w}.
Extended method for two factors
We propose to extend this approach for data with two factors: the time (‘before’ and ‘after’ vaccination), in addition to the stimulation factor. Let be the expression of a given gene k for subject s with stimulation j at time t=1,2. In this context, the mixed model is defined as:
where for a given gene k, is the gene population mean (offset term); measures the fixed effect of stimulation j; measures the fixed effect of time t; is the interaction effect between the stimulation j and the time t; is the random subject effect; measures the random interaction effect between the subject s and the stimulation j; measures the random interaction effect between the subject s and the time t; the residuals and the variables , , and are assumed to be independent. In the context of our application, the potential subject interactions effect with the stimulation and the time effect are confounded with the residuals terms since only one observation is available per subject for each level of both time and stimulation factors.
According to the mixed model, we have:
where the withinsubject deviation can be further decomposed as:
The matrix representation gives:
Similar to the onefactor decomposition, the multivariate approach will be applied on the within matrix , which includes stimulation, time and interaction effects.
Discriminant analysis of one data set
Once the multilevel approach has been applied to split up the variation in the data, a variant of the multivariate approach PLS Discriminant Analysis (called sparse PLSDA) is applied on the within matrix X_{w} or in order to select discriminative genes between the groups of subjects on a single data set.
Sparse PLSDA
Linear Discriminant Analysis (LDA) and Partial Least Squares Discriminant Analysis (PLSDA, [15]) are exploratory approaches seeking the optimal linear combinations of variables (genes) which best separate the sample groups. PLSDA has been found to be a promising alternative to LDA since the latter faces numerical limitations when dealing with too many correlated predictors. Let X(N × p) be the within predictor matrix (to improve readability, the subscript _{w} is removed) and Y(N × G) the response dummy matrix indicating the group of each sample. In PLSDA, X is column standardized. The PLSDA objective function to solve can be written as [15]:
where we denote by ξ=Xuthe discriminant direction vector, which is a linear combination of the original variables. The vector u is the associated loading vector indicating the weights of each variable in the linear combination ξ. Once step (3) has been performed and the first weight vector u_{1} has been extracted, both matrices Y and X are deflated such that the following loading vector u_{2} is orthogonal to the previous one. PLSDA therefore outputs a set of loading weight vectors u_{1}u_{2},…,u_{H} and associated discriminant direction vectors ξ_{1}ξ_{2},…,ξ_{H}, where H is the number of PLSDA dimensions (or deflations).
The sparse version proposed by Lê Cao et al. [16] uses the Lagrange form of PLSDA to include a L_{1} constraint on u in order to ensure that some u_{j} will be estimated as exactly zero (j=1,…,p). Thus, these corresponding variables will not contribute to the discriminant direction. sPLSDA therefore allows variable selection for choosing the variables that best discriminate/separate the sample groups.
Parameters tuning
Two parameters need to be tuned in sPLSDA: the number of discriminant vectors H and the number of variables to select on each dimension (PLS component). Lê Cao et al. [16] showed that for most cases, the user could effectively set H=G−1. The number of variables to select is, however, a challenging issue as tuning criteria are often limited by the very small number of samples. In this study, we considered two criteria to guide the choice of this parameter, both of them are applied sequentially, dimension per dimension.
Tuning criterion 1. One option is to use crossvalidation to choose the optimal number of selected variables in order to avoid selection bias [22]. After this step, the full data are analyzed given this tuned parameter. We propose to estimate the generalization error rate using k fold crossvalidation. In the case of a very small sample size (<15 subjects), the leaveoneout crossvalidation (denoted “loo”) can be used instead. In the specific context of repeated measures and in order to respect the data structure, the training set is composed of the measurements on all experimental units except the measurements on one subject s which defines the test set denoted . The test set prediction is defined by , where β is the regression coefficient matrix from sPLSDA (see [16] for more details). The process is repeated for each subject and the classification error rate is averaged across all subjects. This process is tested for each number of variables to select (see Additional file 1: Figure S3 and Figure S6), and the “optimal” number of variables is then determined when the lowest error rate is obtained.
Additional file 1. Supplementaries results regarding the VAC18 study experiments from two assays.
Format: PDF Size: 3.1MB Download file
This file can be viewed with: Adobe Acrobat Reader
Tuning criterion 2. In the case where the number of subjects is too small, an adhoc alternative approach was used on the whole data set by computing cor(YXu)var(Xu) for each deflated matrix and with respect to the number of selected variables. This is similar to that proposed by Waaijenborg et al. [10] and Parkhomenko et al. [9]. The number of variables selected is chosen to maximize the criterion value.
Integrative analysis of two data sets
Similarly to the PLSDA analysis, a more general PLS multivariate approach can be applied on the matching within matrices X_{w} (or ) and Z_{w} (or ). For this analysis however, the aim is to integrate two data sets in a nonsupervised manner and select correlated variables from both data sets across the subjects.
Sparse PLS
Partial Least Square regression (PLS, [23]) is in fact the ancestor of PLSDA and is applied in a non supervised context, where X(N × p) and Z(N × q) are two continuous within matrices of two different types of predictors (e.g. gene expression and cytokine secretion). In PLS, both X and Z are column standardized. To improve readability, the subscript _{w} is removed from both these matrices. PLS relates X and Z by a linear multivariate model, while also modelling the structure of X and Z. PLS is particularly useful for analysing noisy, collinear, even incomplete, high dimensional data, see [24] for a review. PLS performs successive decompositions of Xand Z into new variables (component scores) denoted by (ξ_{1},…,ξ_{H}) for the Xscores and (ω_{1},…,ω_{H}) for the Zscores. These scores should be few in number (H small), orthogonal to each other within each data set, and estimated as linear combinations of the original variables from X and Z with their weights coefficients indicated in the associated loading vectors u_{h} and v_{h}(h=1,…,H) respectively. In matrix representation, we have X=ΞC^{T} + E, Z=ΩD^{T} + F, where E and F are the residual matrices, and the column matrices in C and D are the coefficients from the local regressions of the score vectors ξ_{h}(ω_{h}) onto the current deflated matrices defined as and , where and .
PLS relates both matrices by maximising the covariance between each pair of scores (ξ_{h},ω_{h}). The PLS objective function is:
This PLS form is often referred to as “PLS2 mode A” in the literature [25] where, similar to Canonical Correlation Analysis, the aim is to model a ‘bidirectional’ relationship between the two data sets (to maximise the common information between the two data sets), as opposed to a ‘unidirectional’ relationship when using a regression model. The sparse version, sPLS, enables variable selection from both sets by including L_{1} penalizations on both u_{h} and v_{h} simultaneously in (4), which is solved with a Lagrangian form (see [7,8] for more details about the methodology and the algorithm). The result is a subset of correlated variables from both X and Z indicated in the loading vectors (u_{h}v_{h}) for each PLS dimension h, and a set of score vectors (ξ_{h}ω_{h}) that are useful for graphical representations.
Parameter tuning
As an extension to the tuning criterion 2 from the previous section, and similar to what was proposed by Waaijenborg et al. [10] and Parkhomenko et al. [9], the number of PLS components and number of variables to select in each step can be tuned by computing cov(X_{h}u_{h}Z_{h}v_{h}), which is the criterion maximized in sPLS2 mode A for each PLS dimension h (see equation (4)). For an optimal number of selected variables from both datasets, one would expect this criterion to achieve also a maximum.
Results and discussion
We first present the results of a short simulation study to show the importance of using a multilevel approach in comparison to a standard sparse partial least square analysis on the original data. We then apply the proposed multilevel approach on an HIVvaccination study.
Simulation study
Simulated model
A simulation study based on the following mixed effects model was performed:
with , , where and are independent. From this model, 10 clusters of 100 genes each were generated (k=1,…,1000). For any given pair of genes k_{1} and k_{2} in the same cluster, a pairwise correlation for and is specified by assuming and , while genes belonging to different clusters are taken to be uncorrelated. The random variables and from same cluster are generated from the multivariate normal distribution where the variancecovariance matrix Σ_{π} is a (100× 100) matrix with along the diagonal and for the others terms; and from the multivariate normal distribution where the variancecovariance matrix Σ_{ε} is a (100×100) matrix with along the diagonal and for the others terms.
To mimic the application, clusters of genes discriminating 4 conditions were generated (the 4 stimulations denoted LIPO5, GAG+, GAG and NS) , where the mean effect of each stimulation is specified by , according to the following:
• 2 gene clusters discriminate (LIPO5, GAG+) versus (GAG, NS) with μ^{k}=(4,4,0,0)^{T}and μ^{k}=(3,3,0,0)^{T}.
• 2 gene clusters discriminate LIPO5 versus GAG+, while GAG+ and NS have the same effect: μ^{k}=(5,2,0.2,0.2)^{T}and μ^{k}=(5,2,0,0)^{T}.
• 2 gene clusters discriminate GAG versus NS, while LIPO5 and GAG+ have the same effect: μ^{k}=c(1,1,5,2)^{T}and μ^{k}=c(0,0,5,2)^{T}.
• the 4 remaining clusters represent noisy signal (no stimulation effect): μ^{k}=c(0,0,0,0)^{T}and μ^{k}=(0.5,0.5,0.5,0.5)^{T}.
The intra cluster correlation was either set to ρ=0.7 or 0.8. Different values for and were studied, but for the sake of conciseness the results are only presented for and .
Numerical results
From the simulated data, the within matrix was computed and applied to multilevel sPLSDA. Figure 1 displays the sample representation for the first 3 axes or dimensions for one simulation run.
Figure 1. Simulation study. Sample representation from multilevel sPLSDA. Samples were projected onto a subspace spanned by the first 3 sPLSDA components, based on the 200 genes selected on each of the 3 components.
Firstly, in order to highlight the benefit of the multilevel approach in comparison to the multivariate approach without the splitup variation step, a prespecified number of genes was selected on each dimension in order to assess the ability of each approach to select the true relevant genes. As expected, 3 components (linear combinations of 200 genes) were sufficient to discriminate the effect of the 4 stimulations. Multilevel sPLSDA (applied on the within matrix) selected 92% of the true simulated discriminative genes as compared to 75% of the true discriminative genes for classical sPLSDA (applied on the original matrix), see Table 1. The hierarchical clustering of the genes selected by sPLSDA on the within matrix (Figure 2) confirmed the discriminatory ability of these genes to separate the 4 groups of samples. As expected, a group of 6 gene clusters can be observed. On the contrary, we did not observe such clusters when applying sPLSDA on the original matrix (not shown).
Figure 2. Stimulation study. Hierarchical clustering (Euclidian distance and Ward method aggregation) of the genes selected with multilevel sPLSDA. Samples are represented in columns and genes in rows.
Table 1. Simulation study
Secondly, leaveoneout crossvalidation was performed on each simulation run to evaluate the error rate of classification of classical sPLSDA or multilevel sPLSDA (Table 2). The classification error rate was evaluated for different number of genes selected on each component. For example, the average classification error rate for multilevel sPLSDA was of 0.009 for 200 genes selected on each of the 3 axes compared to an error rate of 0.268 for the same parameters with classical sPLSDA.
Table 2. Simulation study
Application to HIV vaccine evaluation
Description of the study
The data come from a trial evaluating a vaccine based on HIV1 lipopeptides in HIVnegative volunteers [26]. The vaccine (HIV1 LIPO5 ANRS vaccine) contains five HIV1 amino acid sequences coding for Gag, Pol and Nef proteins. A subsample of 12 vaccinated participants was randomly selected and experiments were performed before and after vaccination. The data consist of monitored cytokine secretion and gene expression measurements from purified in vitro stimulated Peripheral Blood Mononuclear Cells (PBMC). Cytokine secretion was analysed by cytokine multiplex (millipore) in the supernatant of PBMC after 11dayculture, and the data set consists of 10 cytokines measurements (IFNγ, IL1β, IL2, IL5, IL6, IL10, IL13, IL17, IL21, and TNFα). Gene expression was analysed using the Illumina HumanHT12 v4 Expression BeadChip on PBMC before (W0) and 14 weeks after vaccination (W14), 6 hours after in vitro stimulation by either (1) all the peptides included in the vaccine (LIPO5), or (2) the Gag peptides included in the vaccine (GAG+) or (3) the Gag peptides not included in the vaccine (GAG) or (4) without any stimulation (NS).
Preprocessing
Background correction, log_{2} transformation and quantile normalisation were applied on the gene expression data using the Rlimma package. Probes were further prefiltered for each time point (before and after vaccination) using a Pvalue detection (<1% in all samples). The preprocessed data set contained the expression of 25,109 probes for 12 subjects for 4 types of stimulation before vaccination (W0) and the expression of 24,687 probes after vaccination (W14). Some samples were not available due to DNA quality issues, resulting in 44 samples at W0 and 42 samples at W14. For the multilevel approach with two factors, the analysis was performed on the common prefiltered probes before and after vaccination (21,350 probes in total).
The statistical analysis was performed on the probe expression, but the results were biologically interpreted at the gene level.
Discriminant analysis on the transcriptomics data
First we present results obtained using a mixed model and discuss some potential limitations of this method in the context of small sample size. Then we present the results obtained using multilevel sPLSDA for one and twofactor analyses. To shorten the length of the paper, some results have been moved in Additional file 1. The R code used for the analysis of this study is provided in Additional file 2.
Additional file 2. R code used for the analysis of the VAC18 study.
Format: R Size: 15KB Download file
Mixed model
The onelevel mixed model was applied to the W14 transcriptomics data. We used the mle function from the R package nlme with the maximum likelihood method for the estimation of the different models. A global test (likelihood ratio test) followed by an FDR multiple correction (5%) identified 2308 DE genes in at least one of the stimulation. Pairwise comparisons based on Wald test (FDR = 5%) were then performed to compare LIPO5 vs. NS (2108 DE genes), GAG+ vs. NS (1087 DE genes) and GAGvs. NS (209 DE genes). The summary of the results is available in Section 1 of the Additional file 1. In our case study, the clustering analysis of the 100 most significant differentially expressed genes selected by the mixed model failed to discriminate the four stimulations (Additional file 1: Figure S2).
The univariate mixed model approach is commonly used to analyse data with repeated measurement with an unbalanced design. However, several reasons favor the use of a multilevel approach in this high dimensional setting. Apart from the already mentioned problem of numerous independent tests and the requirement to apply multiple correction [27], another limitation is the sensitivity of the FDR threshold (and therefore the number of declared DE genes) to the total number of test performed. The latter depends on the preprocessing method used to filter the probes. Another issue encountered was problems of convergence with both the maximum likelihood (ML) and the restricted ML methods due to the small number of samples in this data set. The asymptotic likelihood ratio test used for fixed effects has been reported to be anticonservative in [28]. The authors recommended to use the Ftest which still poses the issue of the choice of the number of degrees of freedom with a small number of samples [29,30].
Multilevel approach with one factor
A multilevel sPLSDA analysis was performed on the W14 transcriptomics data, with H=3. Respectively 30, 137 and 123 genes were selected with the approach on each dimension according to the tuning criterion 1 for the most parcimonious model. Although k fold crossvalidation would have been preferable to use, loo was used in this study given the small number of subjects. The following ‘loo’ classification error rates (0.48, 0.26, 0.24) were obtained on the first three sPLSDA dimensions compared to (0.48, 0.36, 0.38) when applying sPLSDA on the original matrices (see Additional file 1: Figure S3).
Given the expression of these 290 selected genes, Figures 3(b) and 3(c) highlight a good separation between the four stimulations. These sample representations obtained from sPLSDA reveal that the first component discriminates the stimulation LIPO5 versus the other stimulations, while the second component discriminates the stimulation GAG+ versus the other stimulations and the third component discriminates the stimulation GAG versus the others. Therefore, the first two components separated the stimulations according to the peptides included in the vaccine. As expected there was a clear separation between LIPO5 and other stimulation conditions. Especially, a part of the differential effect of LIPO5 compared to GAG+ could be due to the lipid tail or perhaps to the effect of the other peptides of LIPO5. GAG and NS were not distinguishable on the first two components of the sPLSDA (Figure 3(b)). This last result is not surprising as no specific response is expected from peptides not included in the vaccine after vaccination.
Figure 3. Multilevel sPLSDA analysis on the transcriptomics data with one factor (W14).(a) Unsupervised clustering analysis with Euclidian distance and Ward method of the 290 genes selected by sPLSDA. Samples are represented in columns and genes in rows. (b) and (c) sPLSDA sample representation for dimensions 12 (b) or 13 (c).
Figure 3(a) displays the unsupervised clustering of the 290 selected probes. A cluster of 11 genes MT1M, C20ORF127, MT2A, MT1A, MT1G, MT1F, LOC441019, MT1X, MT1H, MTE, MT1E, was removed to improve visualisation. These genes selected on dimension 1 were all overexpressed in LIPO5 stimulation (see Additional file 1: Figure S5).
Several clusters of genes which expression seemed related to each type of stimulation could be identified. Cluster 1 included a subset of genes downregulated in GAG, in cluster 2 the genes were overexpressed in GAG, while cluster 3 included a subset of genes overexpressed in LIPO5 and GAG+, and cluster 4 was composed of a subset of genes mainly overexpressed in GAG+. The advantage of sPLSDA is its ability to select genes related to a specific stimulation group on each component. For instance, clusters 1 and 2 included 126 out of the 137 probes selected on the third dimension which separated GAG from the other stimulation groups (Figure 3(c)). Cluster 3 included 19 out of the 30 probes selected on the first component, and 12 probes from the second component in order to discriminate stimulations LIPO5 and GAG+, while the fourth cluster included 72 out of the 123 probes selected on the second component which separated GAG+ vs. the other stimulations (Figure 3(b)). Interestingly, some of the genes in this cluster belong to the TNF family (TNFSF13B) or interferon family (ISG20L2) demonstrating a specific effect of the GAG peptides on gene expression related to the immune response.
Note that the same analysis was also performed on W0 but identified much fewer discriminative genes (30 genes in total), indicating that there was a change in expression level after vaccination (see Additional file 1: Figure S8).
Multilevel approach with two factors
A multilevel sPLSDA analysis was performed on the within matrix including the time factor W0 and W14 in addition to the stimulation factor for the transcriptomics data. The complexity of this crossover design implied more conditions (4× 2=8) to be compared for 12 unique subjects. Therefore, the tuning criterion 2 gave for each sPLSDA dimension a maximum correlation of (0.94, 0.96, 0.95) for variable selection sizes of 30, 40, 150 genes on each dimension. A sudden drop in the correlation value in the fourth dimension (0.62) indicated that 3 sPLSDA dimensions should be chosen for this analysis.
The hierarchical clustering of the 220 selected genes indicated a very satisfying separation of both time and stimulation factors (Figure 4(a)). The first component discriminated the stimulations GAG+/LIPO5 vs. GAG/NS irrespective of the time, whereas the second component reflected the time effect W0 vs. W14 (Figure 4(b)). This suggests that the stimulation groups are easier to separate than the time points by the approach. On this second dimension, relevant genes related to the immune response were selected (CD8a, CD79a, CD19, SLAMF6). The third component (Figure 4(c)) separated GAG+ vs. LIPO5 irrespective of the time and several of the genes selected on this third dimension were found to be metallothionein genes (MT1M, MT2A, MT1A, MT1G, MT1F, MT1X, MT1H, MTE, MT1E) that may be stimulated by the lipid tail of LIPO5. From a biological point of view, the significance of genes from metallothionein family in the context of HIV is not clear although some results have been recently reported [31]. These authors showed an increased resistance to apoptosis of immuneactivated monocyte linked to the increase in Metallothionein (MT) gene expression and intracellular zinc levels.
Figure 4. Multilevel sPLSDA analysis on the transcriptomics data with two factors stimulation and time.(a) Unsupervised clustering analysis with Euclidian distance and Ward method of the 220 genes selected by sPLSDA. sPLSDA sample representations for dimensions 12 (b) or 13 (c).
Integrative analysis
Multilevel sPLS enables the integration of data measured using different assays. This approach differs from multilevel sPLSDA as the aim is to select subsets of genes and cytokines which are highly correlated (positively or negatively) across the samples. While the paired structure of the data is still taken into account in the analysis via the decomposition of the within matrices and , the analysis is completely unsupervised: no prior knowledge about the samples groups is included.
Multilevel approach
Multilevel sPLS was applied on the within matrices of the gene and cytokine data sets after vaccination. Given the very small number of cytokines, all cytokines were selected in the model, and the tuning of the number of variables to select was only performed on the gene expression data set. Respectively, a selection of 50, 1 and 60 genes was performed each of the sPLS dimension, corresponding to a correlation of (0.86, 0.62 and 0.84). A drop of the subsequent correlations for the other dimensions guided the choice of 3 components in the model.
Although unexpected and indicated by the tuned correlation value of 0.62, the selection of one single gene on the second dimension was not surprising given the sample representation that was obtained (see Additional file 1: Figure S11): while the first and third dimensions separated LIPO5, GAG+ and GAG/NS, the second dimension did not seem to highlight any interesting pattern in the data. The approach might reveal some unknown phenomenon in the data for this component that would need to be further investigated. Nonetheless, sPLS multilevel was able to identify very relevant information from both data sets. Graphical tools help to unravel the correlation structure between the two data set such as Clustered Image Maps (CIM). Figure 5 reveals clusters of selected genes associated with cytokines secretion. These genes were not known to participate in the cytokines pathways but can be seen as gene signatures to predict future cytokine response. For example Figure 5 highlights relevant clusters of cytokines, such as the proximity of the two THelper type 2 (Th2) cytokines IL5 and IL13. Also, IL17 and IL21 have often been associated in the type 17 response. The correlations between genes and cytokines were similar for the pairs (IL5, IL13), (IL21,IL1b) and (TNF,IL6) underlying potential similar pathways related to the production of these cytokines.
Figure 5. Integrative analysis of gene expression and cytokine secretion for W14. Clustered Image Maps (CIM) obtained from multilevel sPLS. Selected genes are represented in columns and cytokines in rows.
Conclusion
In this paper, we have proposed a twostep analysis combining a multilevel approach and a multivariate approach to analyze repeated measures of gene expression. The multilevel approach first extracts the withinsample variation while the multivariate approach applied on the within matrix takes into account the dependency between the variables. The multilevel approach was extended for one and two factors analyses.
Two multilevel variants were proposed with either sPLSDA or sPLS. The multilevel sPLSDA approach selects genes separating the groups of subjects on a single data set. The simulation study comparing multilevel sPLSDA and the sPLSDA applied on the original data demonstrated the good performance of the model. The multilevel sPLS approach integrates two experiments made on different platforms but on the same subjects, and selects subsets of correlated variables from both sets.
The application of both types of approaches on the HIV1 vaccine trial showed their ability to highlight the stimulation groups and to select biologically relevant genes related to immune response. Hence, our combined multilevel approach may help in finding signatures of vaccine effect and allows for a better understanding of immunological mechanisms activated by the intervention. Future work will include a thorough analysis on the gene/probe annotations to fully understand the mechanistic link between gene differential expression, cytokine secretion according to the various stimulations.
Endnote
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
BL and KL developed the methodology, the R code, performed the simulation and the analysis on the dataset as well as wrote the manuscript. RT developed the methodology, interpreted the dataset as well wrote the Manuscript. HH collected the data and wrote the Application section. All authors read and approved the final manuscript.
Acknowledgements
This work was supported, in part, by the Wound Management Innovation CRC (established and supported under the Australian Government’s Cooperative Research Centres Program) for KA.LC. Financial support of the VAC18 trial was provided by the French Natioanl Agency for Research on AIDS and Hepatitis (ANRS); Sanofi Pasteur provided HIVLIPO5 vaccine. The authors would like to thank the VAC18 study group to provide the data.
References

Storey JD, Xiao W, Leek JT, Tompkins RG, Davis RW: Significance analysis of time course microarray experiments.

Li H, Wood C, Getchell T, Getchell M, Stromberg A: Analysis of oligonucleotide array experiments with repeated measures using mixed models.
BMC Bioinformatics 2004, 5:209.
[http://www.biomedcentral.com/14712105/5/209 webcite]

Karlovich C, DuchateauNguyen G, Johnson A, McLoughlin P, Navarro M, Fleurbaey C, Steiner L, Tessier M, Nguyen T, WilhelmSeiler M, Caulfield J: A longitudinal study of gene expression in healthy individuals.

Palermo RE, Patterson LJ, Aicher LD, Korth MJ, RobertGuroff M, Katze MG: Genomic Analysis Reveals Pre and Postchallenge Differences in a Rhesus Macaque AIDS Vaccine Trial: Insights into Mechanisms of Vaccine Efficacy.

Chun H, Keleş S: Sparse partial least squares regression for simultaneous dimension reduction and variable selection.
J R Stat Soc: Series B (Statistical Methodology) 2010, 72:325.

Lê Cao KA, Rossouw D, RobertGranié C, Besse P: Sparse PLS Variable Selection when Integrating Omics data.

Lê Cao KA, Martin P, RobertGranié C, Besse P: Sparse canonical methods for biological data integration: application to a crossplatform study.

Parkhomenko E, Tritchler D, Beyene J: Sparse canonical correlation analysis with application to genomic data integration.

Waaijenborg S, de Witt Hamer V, Philip C, Zwinderman A: Quantifying the Association between Gene Expressions and DNAMarkers by Penalized Canonical Correlation Analysis.

Witten DM, Tibshirani R, Hastie T: A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis.

Tenenhaus A, Tenenhaus M: Regularized Generalized Canonical Correlation Analysis.

Morine M, Tierney A, van Ommen B, Daniel H, Toomey S, Gjelstad I, Gormley I, PérezMartinez P, Drevon C, LópezMiranda J, et al.: Transcriptomic Coordination in the Human Metabolic Network Reveals Links between n3 Fat Intake, Adipose Tissue Gene Expression and Metabolic Health.

Rose M, Rose T, PariascaTanaka J, Wissuwa M, et al.: Revisiting the role of organic acids in the bicarbonate tolerance of zincefficient rice genotypes.

Barker M, Rayens W: Partial least squares for discrimination.

Lê Cao K, Boitard S, Besse P: Sparse PLS Discriminant Analysis: biologically relevant feature selection and graphical displays for multiclass problems.

Dauwe R, Holliday J, Aitken S, Mansfield S: Metabolic dynamics during autumn cold acclimation within and among populations of Sitka spruce (Picea sitchensis).

de Noord OE, Theobald EH: Multilevel component analysis and multilevel PLS of chemical process data.
J Chemom 2005, 19(57):301307.
[http://dx.doi.org/10.1002/cem.933 webcite]

Westerhuis JA, van Velzen EJ, Hoefsloot HC, Smilde AK: Multivariate paired data analysis: multilevel PLSDA versus OPLSDA.
Metabolomics 2010, 6:119128.
[http://dx.doi.org/10.1007/s113060090185z webcite]

van Velzen EJJ, Westerhuis JA, van Duynhoven JPM, van Dorsten FA, Hoefsloot HCJ, Jacobs DM, Smit S, Draijer R, Kroner CI, Smilde AK: Multilevel Data Analysis of a Crossover Designed Human Nutritional Intervention Study.

Benjamini Y, Hochberg Y: Controlling the False Discovery Rate : A Practical and Powerful Approach to Multiple Testing.

Ambroise C, McLachlan G: Selection bias in gene extraction on the basis of microarray geneexpression data.

Wold H: Multivariate Analysis. Academic Press. Wiley, krishnaiah, p.r. (ed.) edition, New York; 1966.

Boulesteix A, Strimmer K: Partial least squares: a versatile tool for the analysis of highdimensional genomic data.

Vinzi V, Trinchera L, Amato S: PLS path modeling: from foundations to recent developments and open issues for model assessment and improvement.

SalmonCéron D, Durier C, Desaint C, Cuzin L, Surenaud M, Hamouda N, Lelièvre J, Bonnet B, Pialoux G, PoizotMartin I, Aboulker J, Lévy Y, Launay O, trial group AV: Immunogenicity and safety of an HIV1 lipopeptide vaccine in healthy adults: a phase 2 placebocontrolled ANRS trial.

Dudoit S, Shafer JP, Boldrick JC: Multiple Hypothesis Testing in Microarray Experiments.

Pinheiro J, Bates D: MixedEffects Models in S and SPLUS. New York. Springer, USA; 2000.

Kenward MG, Roger JH: Small Sample Inference for Fixed Effects from Restricted Maximum Likelihood.

Schaalje GB, McBride JB, Fellingham GW: Adequacy of Approximations to Distributions of Test Statistics in Complex Mixed Linear Models.

Raymond AD, Gekonge B, Giri MS, Hancock A, Papasavvas E, Chehimi J, Kossevkov AV, Nicols C, Yousef M, Mounzer K, Shull J, Kostman J, Showe L, Montaner LJ: Increased metallothionein gene expression, zinc, and zincdependent resistance to apoptosis in circulating monocytes during HIV viremia.