Abstract
Background
Timecourse gene expression data such as yeast cell cycle data may be periodically expressed. To cluster such data, currently used Fourier series approximations of periodic gene expressions have been found not to be sufficiently adequate to model the complexity of the timecourse data, partly due to their ignoring the dependence between the expression measurements over time and the correlation among gene expression profiles. We further investigate the advantages and limitations of available models in the literature and propose a new mixture model with autoregressive random effects of the first order for the clustering of timecourse geneexpression profiles. Some simulations and real examples are given to demonstrate the usefulness of the proposed models.
Results
We illustrate the applicability of our new model using synthetic and real timecourse datasets. We show that our model outperforms existing models to provide more reliable and robust clustering of timecourse data. Our model provides superior results when genetic profiles are correlated. It also gives comparable results when the correlation between the gene profiles is weak. In the applications to real timecourse data, relevant clusters of coregulated genes are obtained, which are supported by genefunction annotation databases.
Conclusions
Our new model under our extension of the EMMIXWIRE procedure is more reliable and robust for clustering timecourse data because it adopts a random effects model that allows for the correlation among observations at different time points. It postulates genespecific random effects with an autocorrelation variance structure that models coregulation within the clusters. The developed R package is flexible in its specification of the random effects through userinput parameters that enables improved modelling and consequent clustering of timecourse data.
Keywords:
Timecourse data; Mixtures of linear mixed models; Autoregressive random effects; EMMIXWIRE procedureBackground
DNA microarray analysis has emerged as a leading technology to enhance our understanding of gene regulation and function in cellular mechanism controls on a genomic scale. This technology has advanced to unravel the genetic machinery of biological rhythms by collecting massive geneexpression data in a time course. Timecourse gene expression data such as yeast cell cycle data [1] appear to be periodically expressed. To associate the profile of gene expression with a physiological function of interest, it is crucial to cluster the types of gene expression on the basis of their periodic patterns. The identification of coexpressed genes also facilitates the prediction of response to treatment or toxic compounds [2]. Statistical modelling and algorithms play a central role in cataloguing dynamic geneexpression profiles.
Various computational models have been developed for gene clustering based on crosssectional microarray data [35]. Also, considerable attention has been paid to methodological derivations for detecting temporal patterns of gene expression in a time course based on functional principal component analysis or mixture model analysis [615], including the applications to identify differentially expressed genes over time [16,17].
Finite mixture models [18] have been widely used to model the distributions of a variety of random phenomena. Multivariate normality is generally assumed for multivariate data of a continuous nature. The multivariate normal mixture model is employed to detect different patterns in geneexpression profiles. However, when the two assumptions that are commonly adopted in practice, namely,
(1) there are no replications on any particular entity specifically identified as such and
(2) all the observations on the entities are independent of one another,
are violated, multivariate normal mixture models may not be adequate. For example, condition (2) will not hold for the clustering of gene profiles, since not all the genes are independently distributed, and condition (1) will generally not hold either as the gene profiles may be measured over time or on technical replicates. While this correlated structure can be incorporated into the normal mixture model by appropriate specification of the componentcovariance matrices, it is difficult to fit the model under such specifications. For example, the Mstep may not exist in closed form [19].
Accordingly, Ng et al. [13] have developed the procedure called EMMIXWIRE (^{2}EMbased MIXture analysis With Random Effects) to handle the clustering of correlated data that may be replicated. They adopted a mixture of linear mixed models to specify the correlation structure between the variables and to allow for correlations among the observations. It also enables covariate information to be incorporated into the clustering process [13]. Proceeding conditionally on the tissuespecific random effects as formulated in [13], the E and Msteps can be implemented in closed form. In particular, an approximation to the Estep by carrying out timeconsuming Monte Carlo methods is not required. A probabilistic or an outright clustering of the genes into g components can be obtained, based on the estimated posterior probabilities of component membership given the profile vectors and the estimated tissuespecific random effects; see [13].
Fourier series approximations have been used to model periodic gene expression, leading to the detection of periodic signals in various organisms including yeast and human cells [1,20,21]. If the genes studied are periodically regulated, their timedependent expression can be accurately approximated by a Fourier series approximation [20]. A general form of the kth order Fourier series expansion is given as
where a_{0} is the average value of g_{k}(t). The other coefficients a_{k}and b_{k} are the amplitude coefficients that determine the times at which the gene achieves peak and trough expression levels, respectively, and ω is the period of the signal of gene expression. While the timedependent expression value of a gene can be adequately modelled by a Fourier series approximation of the first three orders [14], recent results [13,14] demonstrate that the firstorder Fourier series approximation is sufficient to provide good results in terms of clustering the timecourse data into meaningful functional groups. Alternatively, the likelihood ratio test may be used to determine the order of the Fourier series approximation within the nested regression models.
The EMMIXWIRE procedure of Ng et al. [13] is developed primarily for clustering genes from general microarray experimental designs. On the other hand, Kim et al. [14] focus specifically on clustering periodic gene profiles and propose a special covariance structure to incorporate the correlation between observations at different time points. They also review current methods and compare their method with that of Ng et al. [13]. More recently, Scharl et al. [22] use integrated autoregressive (AR) models to create cluster centers in their simulation study of mixtures of regression models for timecourse gene expression data through the new version of software FlexMix in Leisch [23]. Wang and Fan [24] propose mixtures of multivariate linear mixed models with autoregressive errors to analyse longitudinal data. In this paper, we propose a new EMMIXWIRE normal mixture regression model with AR(1) random effects for the clustering of timecourse data. In particular, the model accounts for the correlation among gene profiles and models the dependence between expressions over time via AR(1) random effects.
The paper is organized as follow: we first present the development of the extension of the EMMIXWIRE model to incorporate AR(1) random effects which are fitted under the EM framework. Then in the following section, we conduct a simulation study and the data analysis with three real yeast cell datasets. In the last section some discussion is provided. The technical details of the derivations are provided in the Additional file 1.
Additional file 1. Supplementary file bmcbioinfsupp2012.pdf
Format: PDF Size: 73KB Download file
This file can be viewed with: Adobe Acrobat Reader
Methods
EMMIXWIRE Model with AR(1) Random Effects
We let X denote the design matrix and βthe associated vector of regression coefficients for the fixed effects. In the specification of the mixture of mixed linear components as adopted by Ng et al. [13], the vector y_{j} for the jth gene conditional on its membership of the hth component of the mixture is expressed as
where β_{h} is a (2k + 1) vector containing unknown parameters a_{0}a_{1},…,a_{k}b_{1},…,b_{k}; see (1), u_{jh}=(u_{jh1},…,u_{jhm})^{T}and v_{h}=(v_{h1},…,v_{hm})^{T}are the random effects, where m is the number of time points. In (2), Z_{1} and Z_{2} are m×m identity matrices. Without loss of generality, we assume ε_{jh}and v_{h} to be independent and normally distributed, N(0,Ω) and N(0,D), independent of u_{jh}. To further account for the time dependent random gene effects, a firstorder autoregressive correlation structure is adopted for the gene profiles, so that u_{jh}follows a N(0,θ^{2}A(ρ)) distribution, where
The inverse of A(ρ) can be expressed as
and
where I, J, and K are all m×mmatrices. Specifically, I is the identity matrix; J has its subdiagonal entries ones and zeros elsewhere, and K takes on the value 1 at the first and last element of its principal diagonal and zeros elsewhere. The expressions (4) and (5) are needed in the derivation of the maximum likelihood estimates of the parameters.
The assumptions (2) and (3) imply that our new model assumes an autocorrelation covariance structure under which measurements at each time point have a larger variance compared to the model of Kim et al. [14] under an AR(1) autocorrelation residual structure.
In the context of mixture models, we consider the gcomponent mixture with probability density function (pdf) as
where f_{h}is the componentpdf of the multivariate normal distribution with mean vector X_{h}β_{h} and covariance matrix
The vector of unknown parameters is denoted by Ψand can be estimated by maximum likelihood via the EM algorithm.
Maximum likelihood via the EM algorithm
In the EM framework adopted here, the observed data vector y=(y_{1},y_{2},…,y_{n})^{T} is augmented by the unobservable component labels, z_{1},z_{2},…,z_{n}of y_{1},y_{2},…,y_{n}, where z_{j} is the gdimensional vector with hth element z_{jh}, which is equal to 1 if y_{j}comes from the hth component of the mixture, and is zero otherwise. These unobservable values are considered to be missing data and are included in the socalled completedata vector. Finally, we take the random effect vectors u_{jh} and v_{h}(j=1,…,n; h=1,…,g), to be missing and include them too in the completedata vector. Now the socalled completedata loglikelihood l_{c} is the sum of four terms l_{c}=l_{1} + l_{2} + l_{3} + l_{4}, where
is the logarithm of the probability of the component labels z_{jh}, and where l_{2} is the logarithm of the density function of y conditional on u_{jh},v_{h}, and z_{jh}=1, and l_{3} and l_{4}is the logarithm of the density function of u and v, respectively, given z_{jh}=1,
where
To maximize the completedata log likelihood l_{c}, the above decomposition implies that each of l_{1},l_{2},l_{3}, and l_{4} can be maximized separately. The EM algorithm proceeds iteratively until the difference between successive values of the log likelihood is less than some specified threshold. All major derivations are given in the Additional file 1.
Results
Simulation study
To illustrate the performance of the proposed model, we present a simulation study based on synthetic timecourse data. In the following simulation, we consider an autocorrelation dependence for the periodic expressions and compare our model to that of Kim et al. [14]. Synthetic timecourse data from three different parametric models (the full model under our new extended EMMIXWIRE approach denoted by EMW in the tables, the extended model of Qin and Self [6], and the model of Kim et al. [14]), assuming a firstorder Fourier series of periodicity, are considered in the simulation study. Within each model, we consider two different settings of θ^{2}corresponding to low and high autocorrelation among the periodic gene expressions. We also assume that Ωand D are diagonal matrices, where the common diagonal elements are represented by σ^{2} and d^{2}, respectively.
There are three clusters of genes. The periods for each cluster are 6, 10, and 16, respectively. There are 24 measurements at time points 0, 1, …, 23, and the first order Fourier expansion is adopted in the simulation models. Parameters and simulation results are listed in Tables 1, 2, 3, 4, 5, 6. In each table, we summarize the results from 1000 simulated sets of data. The true values of the parameters and the biases of their estimates are given in these tables, along with the root mean square errors (RMSEs) in parentheses. We terminated the EM algorithm iterations when the absolute values of the relative changes in all estimates between consecutive iterations were smaller than 0.00001, with the maximum iteration of 1000. For our model, we started from the true partition; for the model of Kim et al. [14], we started from the true values of the parameters. Alternatively, initialization procedures have been considered for mixtures of regression models with and without random effects [22]. For the comparison, we consider the misclassified error rate, the Rand Index, and the adjusted Rand Index [25], where the latter two assess the degree of agreement between the partition and the true clusters of genes. A larger (adjusted) Rand Index indicates a higher level of agreement.
Table 1. Bias and RMSE in brackets from 1000 simulated datasets (generated from new EMMIXWIRE (EMW) model with equal to 0.5)
Table 2. Bias and RMSE in brackets from 1000 simulated datasets (generated from new EMMIXWIRE (EMW) model with equal to 1.3
Table 3. Bias and RMSE in brackets from 1000 simlated datasets (generated from new EMMIXWIRE (EMW) model with equal to 0.5 and d^{2 }equal to 0)
Specifically, we first investigate the performance of our new extended EMMIXWIRE model and that of Kim et al. [14] when the data are generated from the extended EMMIXWIRE model, in which gene expressions within a cluster are correlated. As listed in Tables 1 and 2, the estimates of the parameters pa_{0}a_{1}b_{1}θ^{2}ρ, and σ^{2} in the proposed model are approximately unbiased, except for d^{2}, which is slightly underestimated. In contrast, the model of Kim et al. [14] fails to capture the contributions from genespecific and tissuespecific effects on the autocorrelation among periodic gene expressions at each time point, and thus overestimates the correlation between different time points for each gene. Their method therefore leads to an inferior clustering performance in terms of higher error rates and smaller Rand Indices. From Tables 1 and 2, our proposed method performs better in more than 98% out of 1000 simulated datasets. It also has smaller RMSEs (relative to 0 for error rates and to 1 for Rand Indices) and the difference in performance is significant based on the standard deviation (SD) of the error rates and Rand Indices.
We now compare our model with that of Kim et al. [14], using the data from the extended model of Qin and Self [6], which is a special case of our EMMIXWIRE model (with d^{2} = 0), where gene expressions are independent. The results are presented in Tables 3 and 4, where it can be seen that our method provides essentially unbiased estimates of all the parameters. On the other hand, the model of Kim et al. [14] still overestimates the residual variance at different time points and underestimates the correlation between different time points for each gene, as it fails to capture the contribution from genespecific effects to the autocorrelation among periodic gene expressions at each time point. Their method again produces slightly larger error rates and smaller Rand Indices, though the difference is not significant. From Tables 3 and 4, the proposed method indeed performs better with slightly larger Rand Indices in more than 80% of 1000 simulated datasets.
Table 4. Bias and RMSE in brackets from 100 simulated datasets (generated from new EMMIXWIRE (EMW) model with equal to 1.3 and d^{2 }equal to 0)
Lastly, we generate the data from the model of Kim et al. [14] and provide comparative results in Tables 5 and 6. It is observed from Tables 5 and 6 that the clustering performances are comparable between the two models, as the difference is not significant. Out of 1000 simulated datasets, the proposed method is not worse in more than 35%.
Our model again provides unbiased estimates for all parameters. In contrast to the model of Kim et al. [14], our model accounts for the correlation among gene profiles via the linear effects modelling. As presented in Tables 1 to 6, our model outperforms the model of Kim et al. [14] when the genetic profiles are correlated. When the genetic profiles are generated independently, our model has slightly better performance in cases where the variability in gene expressions at each time point is large. In cases where the residual covariance structure follows an AR(1) model as in Kim et al. [14], our model still provides comparative results and unbiased estimates as with Kim et al. the model of [14]. The advantage of our model to provide more reliable and robust clustering of timecourse data is apparent. With microarray experiments including those timecourse studies, gene expression levels measured from the same tissue sample (or time point) are correlated [19], clustering methods which assume independently distributed gene profiles, such as the model of Kim et al. [14], may overlook important sources of variability in the experiments, resulting in the consequent possibility of misleading inferences being made [13].
Applications: Yeast cell cycle datasets
Yeast cell cycle dataset 1
The first example considers the yeast cell cycle data analysed recently by Wong et al. [26]. This dataset (extracted from Cho et al. [27]) is available from Yeung et al. [28]. It contains 237 genes and 17 samples. These genes are categorized with respect to the four categories in the MIPS database (DNA synthesis and replication, organization of centrosome, nitrogen, and sulphur metabolism, and ribosomal proteins). These categories are assumed to represent the true clusters. In this illustration, we fit our new extended EMMIXWIRE model and the model of Kim et al. [14] to the yeast cell cycle data, with the period of 85 in the Fourier extension [9].
In Table 2 of Wong et al. [26], it shows that the Rand and adjusted Rand Indices for their twostage method are 0.7087 and 0.3697, respectively, and these indices are higher than other methods considered in their paper. Using the model of Kim et al. [14], the Rand indices are 0.7330 and 0.4721, respectively. With the EMMIXWIRE model of Ng et al. [13], we have the Rand and adjusted Rand Indices 0.7799 and 0.5568, respectively. Using the proposed new model, the Rand and adjusted Rand Indices are 0.8123 and 0.6189, respectively, and are the best matches (the largest index) compared with the aforementioned models. The four clusters of genes timecourse profiles are presented in Figure 1. It can be seen that the genes have very similar expression patterns within each cluster, except in cluster 2, where there is greater individual variation by some of the genes. The estimation using the proposed model is listed in Table 7. It can be seen that the correlations in the first three components are from 0.27 to 0.72, indicating a significant correlation among gene expressions at different time points. Ignoring this correlation may therefore lead to a lower Rand Index, that is, a worse clustering. We can see the estimates of d^{2}in clusters 1 and 4 are large and are greater than the corresponding estimates of θ^{2}, indicating coregulation in these two clusters. If we ignore such withincluster coregulation, we will have Rand Indices similar to those for the model of Kim et al. [14]. Our model considers both autocorrelation and coregulation, and thus obtains the best clustering performance.
Figure 1. Clustering of gene expression profiles into four groups for the yeast dataset 1.
Table 7. Estimation of parameters for the yeast cell cycle dataset 1 (237 genes)
Yeast cell cycle dataset 2
The second example is the subset of 384 genes from the yeast cell cycle data in Cho et al. [27], corresponding to five functional groups [28].
Each of gene is assigned a “phase”. We call each “phase” a “Main Group”. There are five “Main Groups” in this dataset, namely, early G1, late G1, S, G2, and M. We now compare and assess the cluster quality with the external criterion (the 5 phases). The raw data are log transformed and normalized by columns and rows. Figure 2 presents the five clusters of genes profiles obtained using the proposed model. It can be seen that the genes have very similar expression patterns within each cluster. The estimations are listed in Table 8. The Rand and adjusted Rand Indices are 0.8102 and 0.4484, respectively. They are 0.8108 and 0.4592 for the model of Kim et al. [14]. The error rates are the same (0.2813) for the two models. The performances of the two models are very similar because the correlation among gene profiles is weak in this dataset. As indicated in Table 8, the estimates of d^{2} are all very small compared to the estimates of θ^{2}.
Figure 2. Clustering of gene expression profiles into five groups for the yeast dataset 2.
Table 8. Estimation of parameters for the yeast cell cycle dataset 2 (384 genes)
A complete Yeast dataset
With this third example, we demonstrate how the proposed method can be adopted to cluster a large amount of yeast genes of which only a small proportion shows periodicity. The original dataset consists of more than 6000 genes, where the yeast cells were sampled at 7 min intervals for 119 min with a total of 18 time points after synchronization [20]. By comparing the ‘aggregate’ numerical score (on the basis of a Fourier algorithm for testing periodicity) for each gene relative to a threshold score, 800 genes were identified as periodically regulated [20]. The threshold score was determined empirically, where 91% of previously known cell cycleregulated genes have a higher score than the threshold. The number of false positives among these 800 cell cycle regulated genes was estimated to be between 3% and 10% [20]. In this study, we worked with 4489 genes that have no missing expression levels across any of the 18 time points. Of these 4489 genes, 612 are periodically regulated and 3877 are not periodic.
The new mixture model with AR(1) random effects and Fourier series approximations was fitted to the periodic gene expression data with the number of clusters g=1 to g=20, where the cell cycle period ω=63 was determined using a global grid search method described in the Discussion section. The optimal number of clusters was determined using the Bayesian Information Criterion (BIC) of Schwarz [29]. The use of BIC for model selection has been considered in the analysis of gene expression data including [8,13,28]. Based on the BIC, there are eight clusters of periodically regulated genes. With the 3877 nonperiodic genes, we adopted the same mixture model with AR(1) random effects, but replacing the Fourier series approximations (1) by a timeseries regression form with Bsplines [8]. Model selection via BIC indicated that there are thirteen clusters of nonperiodic genes. Figure 3 presents the expression profiles of genes in each of the twentyone clusters. From Figure 3(a), it can be seen that the genes have very similar expression patterns within each cluster, except in clusters 3 and 8, where there is greater individual variation by some of the genes. In Table 9, we give the composition of the eight clusters with respect to the five phases of peak expression. Our clusters 2 and 6 consist mainly of those genes with typical G1 peak expression, while most genes in cluster 1 show G2/M or S/G2 phases of peak expression. On the other hand, a majority of genes in cluster 5 has a typical M/G1 phase and those in cluster 7 have a G2/M phase.
Figure 3. Clustering of gene expression profiles into twentyone groups for the complete yeast dataset: (a) eight clusters of periodic genes; (b) thirteen clusters of nonperiodic genes.
Table 9. Distribution of five phases of peak expression over eight clusters obtained (complete yeast data)
With reference to the findings by Spellman et al. [20], a majority of genes in our identified clusters are coregulated. For example, cluster 1 contains genes previously classified to the “CLB2” cluster of Spellman et al. These genes, including ACE2, BUD4, CDC5, and CLB1, are regulated by the MCM1 and SFF transcription factors that induce genes during mitosis [20]. Cluster 8 contain genes described by Spellman et al. as the “MET” cluster. These genes, such as six MET genes and ECM17, are likely to be involved in methionine metabolism [20]. In addition, genes in our clusters 2 and 5 are the major members of the “CLN2” and “SIC1” clusters described in Spellman et al., respectively. The former cluster includes CLN2, CDC9, CDC45, RNR1, POL12, POL30, and are involved in DNA replication and repair. Cluster 5 contains genes that are strongly cell cycle regulated, such as SIC1, TEC1, ASH1, PIR1, and PIR3 [20].
Discussion
We have presented a new mixture model with AR(1) random effects for the clustering of timecourse gene expression profiles. Our new model involves three elements taking important role in modelling timecourse periodic expression data, namely, (a) Fourier expansion which models the periodic patterns; (b) autocorrelation variance structure that accounts for the autocorrelation among the observations at different time points; and (c) the clusterspecific aandom effects which incorporate the coregulation within the clusters. In particular, the latter two elements corresponding to the correlations between timepoints and between genes are crucial for reliable and accurate clustering of timecourse data. We have demonstrated in the simulation and real examples that the accuracy of clustering is improved if the autocorrelation among the time dependent gene expression profiles has been accounted for along the time points; this is also demonstrated in Kim et al. [14]. Furthermore, better results are obtained if the coregulation within the clusters is modelled appropriately. To justify whether an autoregressive correlation structure for the gene expression profiles is appropriate, besides reporting the sample correlation as suggested by one of the reviewers, one can also compare the estimated random components and residual variance with each other. When the correlation between genetic profiles is not small, which is the case for typical timecourse data, ignorance of this dependency may lead to less accurate clustering results. To further illustrate this, we generated 50 gene expression profiles for each of the three clusters from the three parametric models presented in Tables 1, 3, and 5, respectively. The simulated profiles displayed in Figure 4 highlight the differences between the three models we considered in this paper. The results indicate that the proposed time dependent random gene effects are able to capture different degrees of correlation between gene expression profiles within a cluster. Our method performs better when the degree of correlation is large (Figures 4(a) and 4(b)) and provides comparable results when the degree of correlation is small (Figure 4(c)). This finding is also pronounced in the analyses of the real datasets, where yeast cell cycle dataset 1 (Figure 1) shows a high degree of correlation between profiles and yeast cell cycle dataset 2 (Figure 2) shows a relatively small degree of correlation between profiles.
Figure 4. Simulated gene expression profiles for the three models.
As an additional empirical comparison, we applied a simple kmeans clustering procedure to all the simulated and real datasets considered in the paper. We found that the kmeans procedure gave higher error rate and smaller adjusted Rand index (and both with higher variability), especially when the correlation between genetic profiles is not small. For example, the mean (SE) of the error rate and adjusted Rand index obtained for the kmeans procedure for the model in Table 1 are 0.062 (0.057) and 0.849 (0.103), rsepectively. For the model in Table 2, they are 0.357 (0.074) and 0.361 (0.103), respectively. With the model in Table 3, they are 0.187 (0.048) and 0.551 (0.077), respectively, while they are 0.464 (0.038) and 0.153 (0.038), respectively, for the model in Table 4. With the models in Tables 5 and 6, where the degree of correlation is small, the mean (SE) of the error rate and adjusted Rand index obtained for the kmeans procedure are 0.045 (0.018) and 0.876 (0.041), respectively, in Table 5, while they are 0.443 (0.039) and 0.187 (0.041), respectively, in Table 6. For the Yeast 1 real dataset, the adjusted Rand index is 0.509 for the kmeans procedure, which is smaller than the two methods that are based on the EMMIXWIRE model. Using the kmeans procedure, the error rate and adjusted Rand index are 0.404 and 0.442, respectively, for the Yeast 2 dataset. This error rate is the highest among the methods considered in the paper, while the adjusted Rand index is comparable to the other methods. With the complete yeast dataset, the results obtained using the kmeans procedure are somewhat different from those for our method. For example, we have identified a majority of yeast genes (81%) in cluster 5 which show a typical M/G1 phase, while the cluster obtained by the kmeans procedure contains only 69% of genes with a M/G1 phase. Moreover, the clusters obtained by the kmeans procedure for the nonperiodic genes are very different from those presented in Figure 3(b) using our method. These findings indicate that a more complex method is generally required for the clustering of timecourse data, especially when the correlation between the expression levels is not weak.
For the purpose of comparison, the periods of the signal of gene expression are assumed to be known in the simulation study and applications to real data. In practice, there are several ways to estimate the periods for each cluster [9,13,14,20]. For example, in Kim et al. [14], the periods are estimated using simplex algorithm at the Mstep during the EM algorithm. However, when the periods are estimated during the EM iterations, we find that the periods depend also on other parameters. In addition, when we start from an initial period and get the design matrix X, then with higher possibility the best period will be the initial periods. So we change the strategy to a slow one, and we call it global grid search method, which guarantees the highest maximum log likelihood at the best periods. It is implemented as follows. Let S be the space with typical element (a vector) (ω_{1}ω_{2}, …, ω_{g})^{T}, representing the component periods, where ω_{h}can take all possible values (grid points). For example, for the yeast cell cycle data, the possible periods are 60,61, …, 90. Then for each fixed (ω_{1}ω_{2}, …, ω_{g})^{T}, we estimate the parameters as if the periods for each component were known. Finally, we compare the log likelihood and choose the one with the highest log likelihood as the final result. Since it is very slow if there are too many elements in S when we have no prior information about the periods, we recommend using other methods to obtain the periods in such cases, such as the weighted leastsquares estimation approach considered in [15]. In all the calculations in this paper, we assume the periods are fixed.
The proposed model is very flexible through the different specifications of design matrices or model options as originally available in Ng et al. [13]. For example, besides the full model, it enables us to incorporate the model of Qin and Self [6] as a special case. Specifically, we can obtain their model by assuming zero cluster effects (v = 0) and that random effects u be autocorrelated for each gene. Furthermore, when both random effects u and v are assumed to be zero, then we have normal mixture of regression models. In the program we have developed, there are many options and parameters for users to specify the models they want to use in addition to the models we list in our paper. For example, the developed program is applicable to cluster timecourse gene expression profiles that are not periodic (see Figure 3(b)). When periodicity is not obvious, Fourier seris approximations (Equation (1)) can be replaced by a timeseries regression form (such as cubic or spline function) to model the conditional mean expression profiles for each component. With reference to Equation (2), the proposed mixture model framework with time dependent AR(1) random gene effects is again desriable to capture the dependence between the expression measurements over time. The program is written in an R package and is available in Additional file 2, which also contains the data.
Additional file 2. Supplementary file for code and data supp2.zip
Format: ZIP Size: 1.6MB Download file
Conclusions
Our new extended EMMIXWIRE model is more reliable and robust for clustering timecourse data because it postulates genespecific random effects with an autocorrelation variance structure that models coregulation within the clusters. The developed R package is flexible in its specification of the random effects through userinput parameters that enables improved modelling and consequent clustering of timecourse data.
Availability
An Rprogram is available on request from the corresponding author.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
All authors contributed to the production of the manuscript. SKN and GJM directed the research. KW wrote the Rprogram and analysed the simulated and real data. All authors read and approved the final manuscript.
Acknowledgements
This research was supported by a grant from the Australian Research Council.
References

Wichert S, Fokianos K, Strimmer K: Identifying periodically expressed transcripts in microarray time series data.
Bioinformatics 2004, 20:520. PubMed Abstract  Publisher Full Text

Hafemeister C, Costa IG, Schönhuth A, Schliep A: Classifying short gene expression timecourses with Bayesian estimation of.
Bioinformatics 2011, 27:946952. PubMed Abstract  Publisher Full Text

McLachlan GJ, Bean RW, Peel D: A mixture model based approach to the clustering of microarray expression data.

Ramoni MF, Sebastiani P, Kohane IS: Cluster analysis of gene expression dynamics.
Proc Natl Acad Sci USA 2002, 99:91219126. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fan J, Ren Y: Statistical analysis of DNA microarray data in cancer research.
Clin Cancer Res 2006, 12:44694473. PubMed Abstract  Publisher Full Text

Qin LX, Self SG: The clustering of regression models method with applications in gene expression data.
Biometrics 2006, 62:526533. PubMed Abstract  Publisher Full Text

Xu XL, Olson JM, Zhao LP: A regressionbased method to identify differentially expressed genes in microarray time course studies and its application in an inducible Huntingtons disease transgenic model.
Human Mol Genet 2002, 11:19771985. Publisher Full Text

Luan Y, Li H: Clustering of timecourse gene expression data using a mixedeffects model with Bsplines.
Bioinformatics 2003, 19:474482. PubMed Abstract  Publisher Full Text

Luan Y, Li H: Modelbased methods for identifying periodically expressed genes based on time course microarray gene expression data.
Bioinformatics 2004, 20:332339. PubMed Abstract  Publisher Full Text

Storey JD, Xiao W, Leek JT, Tompkins RG, Davis RW: Significant analysis of time course microarray experiments.
Proc Natl Acad Sci USA 2005, 102:1283712842. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hong F, Li H: Functional hierarchical models for identifying genes.
Biometrics 2006, 62:534544. PubMed Abstract  Publisher Full Text

Ma P, CastilloDavis CI, Zhong W, Liu JS: A data driven clustering method for time course gene expression data.
Nucleic Acids Res 2006, 34:12611269. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ng SK, McLachlan GJ, Wang K, B T Jones L, Ng SW: A mixture model with randomeffects components for clustering correlated geneexpression profiles.
Bioinformatics 2006, 22:17451752. PubMed Abstract  Publisher Full Text

Kim BR, Zhang L, Berg A, Fan J, Wu R: A computational approach to the functional clustering of periodic geneexpression profiles.
Genetics 2008, 180:821834. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Booth JG, Casella G, Hobert JP: Clustering using objective functions and stochastic search.
J R Statist Soc 2008, 70:119139. Publisher Full Text

Park T, Yi SG, Lee S, Lee SY, Yoo DH, et al.: Statistical tests for identifying differentially expressed genes in timecourse microarray experiments.
Bioinformatics 2003, 19:694703. PubMed Abstract  Publisher Full Text

Sun W, Wei Z: Multiple testing for pattern identification, with applications to microarray timecourse experiments.
J Am Stat Assoc 2011, 106:7388. Publisher Full Text

McLachlan GJ, Peel D: Finite Mixture Models. New York: John Wiley & Sons; 2000.

McLachlan GJ, Do KA: Analyzing Microarray Gene Expression Data. New Jersey: Wiley; 2004.

Spellman PT, Sherlock S, Zhang MO, Iyer VR, Aners K, et al.: Comprehensive identification of cell cycleregulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization.
Mol biol cell 1998, 9:32733297. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kim BR, Littell RC, Wu RL: Clustering the periodic pattern of gene expression using Fourier series approximations.
Curr Genomics 2006, 7:197203. Publisher Full Text

Scharl T, Grun B, Leisch F: Mixtures of regression models for timecourse gene expression data: evaluation of initialization and random effects.
Bioinformatics 2010, 26:370377. PubMed Abstract  Publisher Full Text

Leisch F: FlexMix: A general framework for finite mixture models and latent class regression in R.

Wang W, Fan T: ECMbased maximum likelihood inference for multivariate linear mixed models with autoregressive errors.
Comput Stat Data Anal 2010, 54:13281341. Publisher Full Text

Hubert L, Arabie P: Comparing partitions.
J Classif 1985, 2:193218. Publisher Full Text

Wong DSV, Wong FK, Wood GR: A multistage approach to clustering and imputation of gene expression profiles.
Bioinformatics 2007, 23:9981005. PubMed Abstract  Publisher Full Text

Cho RJ, Huang M, Campbell MJ, Dong H, Steinmetz L, Sapinoso L, Hampton G, Elledge SJ, Davis RW, Lockhart DJ: Transcriptional regulation and function during the human cell cycle.
Nat Genet 2001, 27:4854. PubMed Abstract  Publisher Full Text

Yeung KY, Fraley C, Murua A, Raftery AE, Ruzzo WL: Modelbased clustering and data transformations for gene expression data.
Bioinformatics 2001, 17:977987. PubMed Abstract  Publisher Full Text

Schwarz G: Estimating the dimension of a model.
Ann Stat 1978, 6:461464. Publisher Full Text