A unified mixed effects model for gene set analysis of time course microarray experiments

Abstract

Methods for gene set analysis test for coordinated changes of a group of genes involved in the same biological process or molecular pathway. Higher statistical power is gained for gene set analysis by combining weak signals from a number of individual genes in each group. Although many gene set analysis methods have been proposed for microarray experiments with two groups, few can be applied to time course experiments. We propose a unified statistical model for analyzing time course experiments at the gene set level using random coefficient models, which fall into the more general class of mixed effects models. These models include a systematic component that models the mean trajectory for the group of genes, and a random component (the random coefficients) that models how each gene's trajectory varies about the mean trajectory. We show that the proposed model (1) outperforms currently available methods at discriminating gene sets differentially changed over time from null gene sets; (2) provides more stable results that are less affected by sampling variations; (3) models dependency among genes adequately and preserves type I error rate; and (4) allows for gene ranking based on predicted values of the random effects. We describe simulation studies using gene expression data with "real life" correlations and we demonstrate the proposed random coefficient model using a mouse colon development time course dataset. The agreement between results of the proposed random coefficient model and the previous reports for this proof-of-concept trial further validates this methodology, which provides a unified statistical model for systems analysis of microarray experiments with complex experimental designs when re-sampling based methods are difficult to apply.

Figure 1

An illustration for the computation of the design matrix for the random effects {r_l; l = 1, …, p} in Model 1. This gene set has 3 genes (variables) and the dataset has 12 samples (observations). Covariance Matrix = estimated gene-gene covariance matrix ∑̂. Under “Eigenvectors”, Prin 1 = the estimated first eigenvector α̂₁ of ∑̂, and λ̂₁ = 0.09802458 is the estimated first eigenvalue of ∑̂. The column in design matrix corresponding to r₁ is then $\sqrt{0.098}{\widehat{\alpha }}_1$, note that they vary according to genes, so the random effects have sub-index i in Model 1.

Figure 2:

ROC Curves for Testing the Central Null Hypothesis H_0C: the average gene expression of a gene group is not differentially expressed over time. The receiver operating characteristic (ROC) curves show a trade-off between sensitivity and 1-specificity as the significance cutoff is varied. Among all models, the random coefficient model MMevct had the best sensitivities across all levels of specificity, the model rstudent performed comparably, Fisher’s exact test lacked sensitivity while globalANCOVA lacked specificity.