A novel approach for biomarker selection and the integration of repeated measures experiments from two assays
 Benoit Liquet^{1, 2, 5}Email author,
 KimAnh Lê Cao^{3},
 Hakim Hocini^{4, 5} and
 Rodolphe Thiébaut^{1, 2, 5}
DOI: 10.1186/1471210513325
© Liquet et al.; licensee BioMed Central Ltd. 2012
Received: 21 May 2012
Accepted: 26 November 2012
Published: 6 December 2012
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/) 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 [6–11] 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 [18–20] 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
where for a given gene k, $\left(\right)close="">{\mu}_{j}^{k}$ measures the fixed effect of stimulation j, which can be further decomposed into $\left(\right)close="">{\mu}_{\xb7\xb7}^{k}$, the overall mean stimulation effect, plus $\left(\right)close="">{\alpha}_{j}^{k}$ which is the differential effect for stimulation j. The $\left(\right)close="">{\pi}_{s}^{k}$ are independent random variables following a normal distribution $\left(\right)close="">\mathcal{N}(0,{\sigma}_{\pi ,k}^{2})$, which take into account the dependency between the repeated measures made on the same subject s, the residuals $\left(\right)close="">{\epsilon}_{\mathrm{sj}}^{k}$ are independent random variables following a $\left(\right)close="">\mathcal{N}(0,{\sigma}_{\epsilon ,k}^{2})$ 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 $\left(\right)close="">{\pi}_{s}^{k}$ and $\left(\right)close="">{\epsilon}_{\mathrm{sj}}^{k}$ 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 $\left(\right)close="">{\pi}_{s}^{k}$, 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
where $\left(\right)close="">{x}_{\xb7j}^{k}=\frac{1}{{n}_{j}}\sum _{s=1}^{{n}_{j}}{x}_{\mathrm{sj}}^{k}$, with n_{ j } the number of subject undergoing stimulation j. Therefore, a part of the withinsubject deviation is explained by the stimulation effect.
The matrix X_{··} represents the offset term defined as $\left(\right)close="">{\mathbf{1}}_{\text{N}}{\mathit{x}}_{\xb7\xb7}^{T}$, where 1_{N} is the (N× 1) matrix containing ones and $\left(\right)close="">{\mathit{x}}_{\xb7\xb7}^{T}=({x}_{\xb7\xb7}^{1},\dots ,{x}_{\xb7\xb7}^{p})$; X_{ b } is the betweensubject matrix of size (N× p) defined by concatenating $\left(\right)close="">{\mathbf{1}}_{{\mathrm{G}}_{s}}{\mathit{x}}_{\mathit{b}s}^{T}$ for each subject into X_{ b } with $\left(\right)close="">{\mathit{x}}_{\mathit{bs}}^{T}=({x}_{s\xb7}^{1}{x}_{\xb7\xb7}^{1},\dots ,{x}_{s\xb7}^{p}{x}_{\xb7\xb7}^{p})$; X_{ w } =X−X_{s·} is the withinsubject matrix of size (N× p), with X_{s·} the matrix defined by concatenating the matrices $\left(\right)close="">{\mathbf{1}}_{{\text{G}}_{s}}{\mathit{x}}_{\mathit{s}\xb7}^{T}$ for each subject into X_{s·}, with $\left(\right)close="">{\mathit{x}}_{\mathit{s}\xb7}^{T}=({x}_{s\xb7}^{1},\dots ,{x}_{s\xb7}^{p})$.
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 [18–20]. 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
where for a given gene k, $\left(\right)close="">{\mu}_{\xb7\xb7}^{k}$ is the gene population mean (offset term); $\left(\right)close="">{\alpha}_{j}^{k}$ measures the fixed effect of stimulation j; $\left(\right)close="">{\beta}_{t}^{k}$ measures the fixed effect of time t; $\left(\right)close="">{\left(\mathrm{\alpha \beta}\right)}_{\mathrm{jt}}^{k}$ is the interaction effect between the stimulation j and the time t; $\left(\right)close="">{\pi}_{s}^{k}\sim \phantom{\rule{1em}{0ex}}N(0,{\sigma}_{{\pi}_{k}}^{2})$ is the random subject effect; $\left(\right)close="">{\left(\mathrm{\alpha \pi}\right)}_{\mathrm{sj}}^{k}\sim \phantom{\rule{1em}{0ex}}N(0,{\sigma}_{\alpha {\pi}_{k}}^{2})$ measures the random interaction effect between the subject s and the stimulation j; $\left(\right)close="">{\left(\mathrm{\beta \pi}\right)}_{\mathrm{st}}^{k}\sim \phantom{\rule{1em}{0ex}}N(0,{\sigma}_{\beta {\pi}_{k}}^{2})$ measures the random interaction effect between the subject s and the time t; the residuals $\left(\right)close="">{\epsilon}_{\mathrm{sjt}}\sim \phantom{\rule{1em}{0ex}}N(0,{\sigma}_{{\epsilon}_{k}}^{2})$ and the variables $\left(\right)close="">{\pi}_{s}^{k}$, $\left(\right)close="">{\epsilon}_{\mathrm{sjt}}^{k}$, $\left(\right)close="">{\left(\mathrm{\alpha \pi}\right)}_{\mathrm{sj}}^{k}$ and $\left(\right)close="">{\left(\mathrm{\beta \pi}\right)}_{\mathrm{st}}^{k}$ 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.
Similar to the onefactor decomposition, the multivariate approach will be applied on the within matrix $\left(\right)close="">{\mathit{X}}_{{w}^{\ast}}$, 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 $\left(\right)close="">{\mathit{X}}_{{w}^{\ast}}$ in order to select discriminative genes between the groups of subjects on a single data set.
Sparse PLSDA
where we denote by ξ=X u the 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 $\left(\right)close="">{\mathit{X}}_{w,s}^{\mathrm{test}}$. The test set prediction is defined by $\left(\right)close="">{\mathit{Y}}^{\mathrm{test}}={\mathit{X}}_{w,s}^{\mathrm{test}}\mathit{\beta}$, 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.
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(Y X u)var(X u) 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 $\left(\right)close="">{\mathit{X}}_{{w}^{\ast}}$) and Z_{ w } (or $\left(\right)close="">{\mathit{Z}}_{{w}^{\ast}}$). 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 X and 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 $\left(\right)close="">{\mathit{X}}_{h}={\mathit{X}}_{h1}{\mathit{\xi}}_{h}{\mathit{c}}_{h}^{\prime}$ and $\left(\right)close="">{\mathit{Z}}_{h}={\mathit{Z}}_{h1}\mathit{\omega}{\mathit{d}}_{h}^{\prime}$, where $\left(\right)close="">{\mathit{c}}_{h}={\mathit{X}}_{h1}^{T}{\mathit{\xi}}_{h}/{\mathit{\xi}}_{h}^{\prime}{\mathit{\xi}}_{h}$ and $\left(\right)close="">{\mathit{d}}_{h}={\mathit{Y}}_{h1}^{T}{\mathit{\omega}}_{h}/{\mathit{\omega}}_{h}^{\prime}{\mathit{\omega}}_{h}$.
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
with $\left(\right)close="">{\pi}_{s}^{k}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\sim \phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}N(0,{\sigma}_{{\pi}_{k}}^{2})$, $\left(\right)close="">{\epsilon}_{\mathrm{sj}}^{k}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\sim \phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}N(0,{\sigma}_{{\epsilon}_{k}}^{2})$, where $\left(\right)close="">{\pi}_{s}^{k}$ and $\left(\right)close="">{\epsilon}_{\mathrm{sj}}^{k}$ 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 $\left(\right)close="">{X}_{\mathrm{sj}}^{{k}_{1}}$ and $\left(\right)close="">{X}_{\mathrm{sj}}^{{k}_{2}}$ is specified by assuming $\left(\right)close="">\mathrm{cor}({\pi}_{s}^{{k}_{1}},{\pi}_{s}^{{k}_{2}})=\rho $ and $\left(\right)close="">\mathrm{cor}({\epsilon}_{\mathrm{sj}}^{{k}_{1}},{\epsilon}_{\mathrm{sj}}^{{k}_{2}})=\rho $, while genes belonging to different clusters are taken to be uncorrelated. The random variables $\left(\right)close="">{\pi}_{s}^{k}$ and $\left(\right)close="">{\epsilon}_{\mathrm{sj}}^{k}$ from same cluster are generated from the multivariate normal distribution $\left(\right)close="">({\pi}_{s}^{{k}_{1}},\dots ,{\pi}_{s}^{{k}_{100}})\sim {N}_{100}({0}_{100},{\Sigma}_{\pi})$ where the variancecovariance matrix Σ_{ π } is a (100× 100) matrix with $\left(\right)close="">{\sigma}_{{\pi}_{k}}^{2}$ along the diagonal and $\left(\right)close="">\rho {\sigma}_{{\pi}_{k}}^{2}$ for the others terms; and from the multivariate normal distribution $\left(\right)close="">({\epsilon}_{\mathrm{sj}}^{{k}_{1}},\dots ,{\epsilon}_{\mathrm{sj}}^{{k}_{100}})\sim {N}_{100}({0}_{100},{\Sigma}_{\epsilon})$ where the variancecovariance matrix Σ_{ ε } is a (100×100) matrix with $\left(\right)close="">{\sigma}_{\epsilon}^{2}$ along the diagonal and $\left(\right)close="">\rho {\sigma}_{{\epsilon}_{k}}^{2}$ 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 $\left(\right)close="">{\mathit{\mu}}^{k}={({\mu}_{1}^{k},{\mu}_{2}^{k},{\mu}_{3}^{k},{\mu}_{4}^{k})}^{T}$, 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 $\left(\right)close="">{\sigma}_{{\pi}_{k}}^{2}$ and $\left(\right)close="">{\sigma}_{{\epsilon}_{k}}^{2}$ were studied, but for the sake of conciseness the results are only presented for $\left(\right)close="">{\sigma}_{{\pi}_{k}}=2$ and $\left(\right)close="">{\sigma}_{{\epsilon}_{k}}=0.5$.
Numerical results
Simulation study
Component 1  Component 2  Component 3  All  

classical sPLSDA  58.0  75.0  87.2  78.2 
multilevel sPLSDA  82.8  95.6  93.1  92.0 
Simulation study
Number of genes  Original matrix  Within matrix  

1 component  2 components  3 components  1 component  2 components  3 components  
25  0.535  0.369  0.312  0.500  0.271  0.024 
50  0.530  0.364  0.311  0.500  0.265  0.016 
75  0.527  0.360  0.306  0.500  0.261  0.013 
100  0.524  0.354  0.300  0.500  0.258  0.011 
125  0.522  0.351  0.296  0.500  0.257  0.009 
150  0.520  0.343  0.285  0.500  0.250  0.008 
175  0.518  0.335  0.281  0.500  0.243  0.009 
200  0.516  0.327  0.268  0.500  0.234  0.009 
225  0.514  0.323  0.269  0.500  0.227  0.009 
250  0.512  0.316  0.267  0.500  0.220  0.008 
275  0.510  0.314  0.266  0.500  0.207  0.007 
300  0.510  0.306  0.262  0.500  0.196  0.007 
325  0.509  0.299  0.260  0.500  0.182  0.007 
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 R limma 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.
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).
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 $\left(\right)close="">{\mathit{X}}_{{w}^{\ast}}$ 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.
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 $\left(\right)close="">{\mathit{X}}_{\mathit{w}}^{\mathbf{\ast}}$ and $\left(\right)close="">{\mathit{Z}}_{\mathit{w}}^{\mathbf{\ast}}$, 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.
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
Declarations
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.
Authors’ Affiliations
References
 Storey JD, Xiao W, Leek JT, Tompkins RG, Davis RW: Significance analysis of time course microarray experiments. Proc National Acad Sci USA 2005, 102(36):12837–12842.View ArticleGoogle Scholar
 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/1471–2105/5/209] []PubMed CentralView ArticlePubMedGoogle Scholar
 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. BMC Med Genomics 2009, 2: 33.PubMed CentralView ArticlePubMedGoogle Scholar
 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. J Virol January 15 2011, 85(2):1099–1116.PubMed CentralView ArticlePubMedGoogle Scholar
 Kitano H: Computational Systems Biology. Nature 6912, 420: 206–210.View ArticleGoogle Scholar
 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: 3–25.View ArticleGoogle Scholar
 Lê Cao KA, Rossouw D, RobertGranié C, Besse P: Sparse PLS Variable Selection when Integrating Omics data. Stat App Mol Biol 2008, 7(1):37.Google Scholar
 Lê Cao KA, Martin P, RobertGranié C, Besse P: Sparse canonical methods for biological data integration: application to a crossplatform study. BMC Bioinformatics 2009., 10(34):Google Scholar
 Parkhomenko E, Tritchler D, Beyene J: Sparse canonical correlation analysis with application to genomic data integration. Stat App Genet Mol Biol 2009, 8(1):1–34.Google Scholar
 Waaijenborg S, de Witt Hamer V, Philip C, Zwinderman A: Quantifying the Association between Gene Expressions and DNAMarkers by Penalized Canonical Correlation Analysis. Stat App Genet Mol Biol 2008., 7(3):Google Scholar
 Witten DM, Tibshirani R, Hastie T: A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis. Biostatistics 2009, 10(3):515–534.PubMed CentralView ArticlePubMedGoogle Scholar
 Tenenhaus A, Tenenhaus M: Regularized Generalized Canonical Correlation Analysis. Psychometrika 2011, 76(2):257–284.View ArticleGoogle Scholar
 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. PLoS comput biol 2011, 7(11):e1002223.PubMed CentralView ArticlePubMedGoogle Scholar
 Rose M, Rose T, PariascaTanaka J, Wissuwa M, et al.: Revisiting the role of organic acids in the bicarbonate tolerance of zincefficient rice genotypes. Funct Plant Biol 2011, 38(6):493–504.Google Scholar
 Barker M, Rayens W: Partial least squares for discrimination. J Chemom 2003, 17(3):166–173.View ArticleGoogle Scholar
 Lê Cao K, Boitard S, Besse P: Sparse PLS Discriminant Analysis: biologically relevant feature selection and graphical displays for multiclass problems. BMC bioinformatics 2011, 12: 253.PubMed CentralView ArticlePubMedGoogle Scholar
 Dauwe R, Holliday J, Aitken S, Mansfield S: Metabolic dynamics during autumn cold acclimation within and among populations of Sitka spruce (Picea sitchensis). New Phytologist 2012, 194(1):192–205.View ArticlePubMedGoogle Scholar
 de Noord OE, Theobald EH: Multilevel component analysis and multilevel PLS of chemical process data. J Chemom 2005, 19(5–7):301–307. [http://dx.doi.org/10.1002/cem.933] []View ArticleGoogle Scholar
 Westerhuis JA, van Velzen EJ, Hoefsloot HC, Smilde AK: Multivariate paired data analysis: multilevel PLSDA versus OPLSDA. Metabolomics 2010, 6: 119–128. [http://dx.doi.org/10.1007/s11306–009–0185z] []PubMed CentralView ArticlePubMedGoogle Scholar
 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. J Proteome Res 2008, 7(10):4483–4491.View ArticlePubMedGoogle Scholar
 Benjamini Y, Hochberg Y: Controlling the False Discovery Rate : A Practical and Powerful Approach to Multiple Testing. J R Stat Soc. Series B 2009, 57: 289–300.Google Scholar
 Ambroise C, McLachlan G: Selection bias in gene extraction on the basis of microarray geneexpression data. Proc National Acad Sci 2002, 99(10):6562.View ArticleGoogle Scholar
 Wold H: Multivariate Analysis. Academic Press. Wiley, krishnaiah, p.r. (ed.) edition, New York; 1966.Google Scholar
 Boulesteix A, Strimmer K: Partial least squares: a versatile tool for the analysis of highdimensional genomic data. Briefings in Bioinf 2007, 8: 32.View ArticleGoogle Scholar
 Vinzi V, Trinchera L, Amato S: PLS path modeling: from foundations to recent developments and open issues for model assessment and improvement. Handbook of Partial Least Squares 2010, chapter 2: 47–82.View ArticleGoogle Scholar
 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. AIDS 2010, 24(14):2211–2223.View ArticlePubMedGoogle Scholar
 Dudoit S, Shafer JP, Boldrick JC: Multiple Hypothesis Testing in Microarray Experiments. Stat Sci 2003, 18: 71–103.View ArticleGoogle Scholar
 Pinheiro J, Bates D: MixedEffects Models in S and SPLUS. New York. Springer, USA; 2000.View ArticleGoogle Scholar
 Kenward MG, Roger JH: Small Sample Inference for Fixed Effects from Restricted Maximum Likelihood. Biometrics 1997, 53(3):983–997.View ArticlePubMedGoogle Scholar
 Schaalje GB, McBride JB, Fellingham GW: Adequacy of Approximations to Distributions of Test Statistics in Complex Mixed Linear Models. J Agric, Biol, and Environ Stat 2002, 7(4):512–524.View ArticleGoogle Scholar
 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. J Leukocyte Biol 2010, 88(3):589–596.PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.