article: Gene set analysis methods, which consider predefined groups of genes in the analysis of genomic data, have been successfully applied for analyzing gene expression data in cross-sectional studies. The time-course gene set analysis (TcGSA) introduced here is an extension of gene set analysis to longitudinal data. The proposed method relies on random effects modeling with maximum likelihood estimates. It allows to use all available repeated measurements while dealing with unbalanced data due to missing at random (MAR) measurements. TcGSA is a hypothesis driven method that identifies a priori defined gene sets with significant expression variations over time, taking into account the potential heterogeneity of expression within gene sets. When biological conditions are compared, the method indicates if the time patterns of gene sets significantly differ according to these conditions. The interest of the method is illustrated by its application to two real life datasets: an HIV therapeutic vaccine trial (DALIA-1 trial), and data from a recent study on influenza and pneumococcal vaccines. In the DALIA-1 trial TcGSA revealed a significant change in gene expression over time within 69 gene sets during vaccination, while a standard univariate individual gene analysis corrected for multiple testing as well as a standard a Gene Set Enrichment Analysis (GSEA) for time series both failed to detect any significant pattern change over time. When applied to the second illustrative data set, TcGSA allowed the identification of 4 gene sets finally found to be linked with the influenza vaccine too although they were found to be associated to the pneumococcal vaccine only in previous analyses. In our simulation study TcGSA exhibits good statistical properties, and an increased power compared to other approaches for analyzing time-course expression patterns of gene sets. The method is made available for the community through an R package.
  Microarray experiments are increasingly used for evaluating changes in gene expression over time. The analysis of the temporal change of gene expression should help in understanding the complex mechanisms of gene regulation. For instance, transcriptional profiles have been repeatedly measured to study the change in gene expression during the natural history of SIV/HIV infection [1,2] or to evaluate the effect of vaccines [3,4]. In the applications considered in this paper [5,6], the investigators wanted to detect the genes for which the abundance changed over time after a vaccination (against HIV, influenza or pneumococcus). In order to analyze such longitudinal high-dimensional data, several approaches have been suggested including a gene-by-gene statistical analysis [7,8], dimension reduction methods [9] or gene set analysis [10]. A gene set is a group of genes that are a priori co-regulated or functionally linked. Examples of such gene set relating to biological processes or pathways are those defined by KEGG [11], Gene Ontology [12] or Chaussabel’s functional modules [13]. The gene set analysis [14–16] is supposed to be more powerful than a gene-by-gene analysis because it can detect a change of expression of a group of genes although none of them show a very high absolute fold change. Furthermore, a change of all genes in a given pathway may be biologically more meaningful than a large increase of a single gene. Also, provided that the gene sets are well defined, the result should be more sound and comparable across studies than a gene-by-gene analysis [14]. Finally, gene set analysis avoids a second step for a global interpretation as described in the “bottom up” approach [10,17]. The analysis of longitudinal microarray experiments through a gene set approach is not trivial because the dynamics of gene expressions inside a gene set can be complex and heterogeneous. This has already been underlined in some of the approaches developed to analyze gene sets [15,18–20]. Fig 1 shows an example of a homogeneous gene set, whereas Fig 2 shows an example of a heterogeneous one. Actually, such a heterogeneity is frequently observed [20], and cannot be ignored, as genes inside a functional gene set are not expected to change their expression synchronously. Moreover this heterogeneity can be biologically meaningful by itself. Prieto et al. [21] provide an example from a cancer application, where deregulated pathways are of primary biological interest. They identified heterogeneous gene sets linked to acute promyelocytic leukemia. Another example is given by Hu et al. : pathways affected by the HER2, such as the KEGG pathways of ‘Ubiquitin mediated proteolysis’, ‘Glioma’, and ‘Prostate cancer’ were identified by studying heterogeneity [22]. The main advantage of detecting the heterogeneity inside a gene set is to detect any change over time whatever the specification of the model for the trends. In other words, the dynamics of gene expression inside a stable gene set will be summarized by a flat slope and no heterogeneity. Hence, in the spirit of [19], to find any significant change of the overall expression of genes inside a gene set over time, we suggest to look for any significant trend over time or any heterogeneity between gene trends inside the gene set. Several approaches have already been proposed to analyze longitudinal measurements of gene expression [7,23–28], but only a few include gene set analysis [10,18,29–31]. Among the latter, all but Nueda et al. [18] fail to account for possible heterogeneity inside a gene set. An extension of the popular Gene Set Enrichment Analysis (GSEA) method [14] is available for the analysis of time series data: GSEA for time series. Unfortunately, it does not account for the structure of longitudinal data, simply treating all observations as independent and calculating Pearson correlation of the genes with the time variable. Therefore this is quite a limited modeling strategy for longitudinal data. The globalANCOVA procedure developed by Hummel et al. [29] focuses on the comparison of groups, testing whether there is a group influence on change over time of any gene expression inside a gene set. In practice, the global null hypothesis tested is quite flexible relying on the ANOVA framework, but cannot accommodate missing values. Wang et al. [10] proposed to use a linear mixed effects model to explain gene expression inside a gene set. They considered a random effect for the array level rather than for the patient or the gene level. Zhang et al. [30] proposed a robust non-parametric approach to compare gene expression dynamics between different treatment-groups. Wu et al. [31] proposed the CAMERA procedure, a gene set test based on linear modeling that takes into account inter-gene correlations. However, CAMERA accounts for neither potential heterogeneity inside a gene set nor for repeated measurement correlation. Of note, it is not possible to look at the change of gene expression in only one group using Zhang et al. [30] approach. In contrast, the PCA-maSigFun procedure developed by Nueda et al. [18,32] can account for possible heterogeneity inside a gene set. It is based on a Principal Component Analysis (PCA) of each gene set followed by a linear regression of the significant principal components (i. e. components that have a variation above the mean gene variance) over time. However, they did not consider time-course experiments where repeated measures are available for multiple patient. Gene set analysis methods can also be distinguished by their choice of the null hypothesis. Those can be classified into two main types of hypothesis: i) the competitive null hypothesis, that tests the genes inside a given gene set against all the other genes outside the gene set; ii) the self-contained null hypothesis, that only uses the genes inside the gene set of interest [20,33,34]. In the present paper, interest is focused on self-contained null hypotheses because the question was “Which gene sets have a change of gene abundance over time? ”. According to Emmert-Streib et al. [35], self-contained gene set tests are biologically easier to interpret and can be more powerful compared to competitive tests. Self-contained gene set tests are especially appropriate in a hypothesis driven context where a priori defined gene sets are validated and relevant regarding the biological question. We propose the implementation of a hypothesis driven method that directly tests the time-course significance of predefined gene sets: the Time-course Gene Set Analysis (TcGSA). It relies on the use of linear mixed effect models, a very useful and well established statistical tool [36,37] especially suited for longitudinal settings. By using all available repeated measures, it provides increased statistical power. TcGSA can accommodate for heterogeneity of gene expression within the gene sets through random effects, and is robust to unbalanced designed due to missing (at random) values thanks to the maximum likelihood estimates. No previously proposed approach combines all of TcGSA features. A simulation study demonstrated the good statistical performance of the proposed method. It has been applied to two studies: one HIV vaccine trial, and one influenza and pneumoccocal vaccine study [6], using the same definition of gene sets [13] that is increasingly used in systems immunology research [38–42]. Compared to gene-by-gene analyses, TcGSA disclosed changes of additional gene sets that endorse previous conclusions [6], but also revealed common pathways across the three vaccines.
  The TcGSA method includes three steps: 1) modeling gene expression in a gene set with mixed models, 2) testing the significance of a gene set, and 3) estimating individual gene profiles. The TcGSA method has been implemented in R as a package TcGSA, whose latest release is available from the CRAN repository (http: //cran. r-project. org/web/packages/TcGSA/index. html).
  The DALIA-1 trial is a phase 1 therapeutic HIV vaccine trial whose details are described on http: //clinicaltrials. gov (ClinicalTrial. gov identifier: NCT00796770) and in [5]. The vaccine candidate was based on ex-vivo generated interferon-α dendritic cells loaded with HIV-1 lipopeptides and activated with lipopolysaccharide. The objectives of the trial were to evaluate the safety of the strategy and to evaluate the immune response to the vaccine. For the purpose of the present paper, we focus on the gene expression component of this study. Gene abundance in whole blood was measured through Illumina HumanHT-12 v4 Expression BeadChips. In a recent paper, Obermoser et al. [6] investigated the response to influenza and pneumococcal vaccines in healthy individuals at the gene expression level. In order to assess the behavior of the proposed method, a simulation study of TcGSA has been performed. The simulation scheme was chosen to be very close to the motivating exemple: the DALIA-1 trial. In each simulation run, gene expression data was simulated for 20 patients at 8 time points. 5,000 genes were simulated, divided into 100 non overlapping gene sets of 50 genes each. Each of the 100 gene sets was either simulated under (H0) or (H1). The proportion of genes under (H1) varied between 0%, 27% (which corresponds to results found in pre-ATI) and 85% (which corresponds to results found in post-ATI). When there are gene sets simulated under (H1), 75% of those were homogeneous (simulated with parameters close to those estimated for gene set M1. 2 in DALIA post-ATI—see Fig 4) while the remaining 25% were heterogeneous (simulated with parameters close to those estimated for gene set M7. 1 in DALIA post-ATI—see Fig 4). Statistical performances of the proposed method are presented in Table 2. Without correcting for the fact that 100 gene set were investigated by TcGSA at each simulation runs, the average Type I error (the probability of rejecting H0 given that H0 is actually true) over a hundred runs was between 0. 03 and 0. 07 depending on the situation. But as soon as a control of the FDR was used, the Type-I error rate dropped well below 1%, regardless of the flexibility of the time function estimated (linear or cubic polynomials). The average statistical power (the probability of rejecting H0 given that H0 is actually false) is very good, always above 0. 8 (dropping a little bit after multiple testing correction as expected). Three other methods were also evaluated on those simulations, namely globalANCOVA [29] (using either permutations or an approximation to compute p-values), GSEA for time series and CAMERA. Their statistical performances are also presented in Table 2. Type I error is always well controlled by those three methods. However, both CAMERA and GSEA for time series exhibit very low statistical power (as low as 10 times less than TcGSA after multiple testing correction when there is a high proportion of significant gene sets). globalANCOVA, whose global null hypothesis is not so different from the one tested in TcGSA, performs quite well in terms of statistical power. Nonetheless it is still about 10% below TcGSA performances. Those simulation results confirm that the higher number of selected gene sets by TcGSA in the two real-life examples presented in this paper are mainly due to the increased power of gene set analysis over gene-by-gene analysis (when repeated structure of the measurement is properly accounted for), and not to a large number of false positives.
  In this paper, we present a method to analyze repeated measurements of gene expression using a gene set approach. Provided that the definition of the gene sets is relevant, this method helps with detecting and interpreting subtle changes of gene expression over time. In our applications where the same definition of gene sets has been applied, we were able to compare the response to several vaccines (against HIV, Influenza and Pneumococcus). Interestingly, we found common pathways that were triggered by all three vaccines, mostly related to inflammation, as well as pathways specific to each vaccine. The capacity of the proposed approach to detect subtle changes of gene expression is due to two main factors: i) the use of a predefined gene sets that are functionally related ii) the use of all available information, taking advantage of repeated measurements using mixed models. Measurements of gene expression data in longitudinal studies may be missing because of missed visits or poor quality of the samples, leading to unbalanced data. Missing at random (MAR) processes (i. e. when the probability of missing data is associated to the previously measured information) may lead to biased estimates when using least squares or generalized estimating equations [37]. TcGSA can cope with such issues because of the use of Maximum Likelihood to estimate the parameters of the mixed models. This is an advantage of the TcGSA approach over those of Hummel et al. [29] or Nueda et al. [32]. An increasing number of gene sets databases are available, such as KEGG [11], Gene Ontology [12], Modules [13]. An immune related subset of Gene Ontology as well as an immune related subset of KEGG pathway have been used in additional analyses (see S3 Fig). The choice of the database used for the analysis impacts the interpretability but also the limitations of TcGSA. The efficiency of TcGSA will vary according to the number of genes represented in each gene set. The size of a given gene set has an impact on its significance, as the more genes it includes, the more likely a significant variation will be detected. The average size of the Chaussabel’s V2 modules is 55 genes. 17% of the 260 modules include more than 100 genes, and 31% less than 20 genes. For small gene sets, the normality assumptions of random effects of the Eqs (1), (1bis) and (2) are questionable (even though inference should still be robust [59]). Nevertheless, even though we expect that the models could be misspecified in many cases (if not all), the objective of such an analysis is to detect any significant variation over time (in the spirit of Shahbaba et al. [19], a significant variability of the trajectories between the genes inside a gene set indicates a change over time regardless of the fixed effects specification). The use of flexible time functions may help to get a better fit of dynamics although beyond cubic polynomials it did not have a substantial impact in our motivating example—see Table 1. These results vary according to the dataset and the number of time points available, and we recommend to try several models to check the robustness of the results. TcGSA focuses on self-contained hypothesis. While competitive gene sets tests are appropriate in the presence of a strong signal (e. g. when comparing pre-ATI versus post-ATI), self-contained approaches are able to detect more subtle and sparse signals (e. g in pre-ATI). The choice between a competitive or a self-contained approach is therefore highly dependent on the biological question. Several extensions of the TcGSA are possible for its use in other contexts. One can also model time trends with a random effect grouped on the patient level as of γp, ⋅, instead of on the gene level as in Eqs (1), (1bis) and (2). This identifies gene sets whose dynamic differs across the patients. This option is also implemented in the TcGSA R package. TcGSA could easily be adapted to mRNA counts data. In that case, generalized linear mixed effects models could be used, with a Poisson distribution for instance, instead of linear mixed effect models that rely on a Gaussian assumption. Or else, precision weights estimated with voom [60] could also be used to model the mean-variance relationship while fitting the proposed linear mixed effect directly on log-counts data from RNA-seq, thus taking into account heteroscedasticity. In conclusion, the method presented gives a solution for the full exploitation of any repeated measurements of gene expression data based on a gene set analysis where a great sensibility to detect subtle change, while controlling false discovery, is needed.
summary: Gene set analysis methods use prior biological knowledge to analyze gene expression data. This prior knowledge takes the form of predefined groups of genes, linked through their biological function. Gene set analysis methods have been successfully applied in transversal studies, their results being more sensitive and interpretable than those of methods investigating genomic data one gene at a time. The time-course gene set analysis (TcGSA) introduced here is an extension of such gene set analysis to longitudinal data. This method identifies a priori defined groups of genes whose expression is not stable over time, taking into account the potential heterogeneity between patients and between genes. When biological conditions are compared, it identifies the gene sets that have different expression dynamics according to these conditions. Data from 2 studies are analyzed: data from an HIV therapeutic vaccine trial, and data from a recent study on influenza and pneumococcal vaccines. In both cases, TcGSA provided new insights compared to standard approaches thanks to an increased sensitivity compared to other approaches. Those results highlight the benefits of the TcGSA method for analyzing gene expression dynamics.
title: Time-Course Gene Set Analysis for Longitudinal Gene Expression Data