Nonlinear Model-Based Method for Clustering Periodically Expressed Genes

ABSTRACT

Clustering periodically expressed genes from their time-course expression data could help understand the molecular mechanism of those biological processes. In this paper, we propose a nonlinear model-based clustering method for periodically expressed gene profiles. As periodically expressed genes are associated with periodic biological processes, the proposed method naturally assumes that a periodically expressed gene dataset is generated by a number of periodical processes. Each periodical process is modelled by a linear combination of trigonometric sine and cosine functions in time plus a Gaussian noise term. A two stage method is proposed to estimate the model parameter, and a relocation-iteration algorithm is employed to assign each gene to an appropriate cluster. A bootstrapping method and an average adjusted Rand index (AARI) are employed to measure the quality of clustering. One synthetic dataset and two biological datasets were employed to evaluate the performance of the proposed method. The results show that our method allows the better quality clustering than other clustering methods (e.g., k-means) for periodically expressed gene data, and thus it is an effective cluster analysis method for periodically expressed gene data.

1. BACKGROUND

Many biological processes such as cell-cycle division exhibit periodic behaviors. To understand the mechanisms of these biological processes, DNA microarray experiments have been employed to produce gene expression profiles at a series of time points, for example, the cell division cycle processes of yeast Saccharomyces cerevisiae [1, 2], bacterium Caulobacter crescentus [3], and human being [4]. Such time-course gene expression data provides a dynamic snapshot of most (if not all) of the genes related to the biological development process. It is believed that clustering periodically expressed gene from their time-course expression data could help understand the molecular mechanisms of those biological processes.

In past decade, a number of methods have been proposed for identifying and clustering periodically expressed genes. The discrete Fourier transform method is the earliest method for identifying and clustering periodically expressed genes [1–4]. In these papers, the discrete Fourier transform is applied to gene expression data to get a two-dimensional vector. One component of the vector is the sum of all coefficients of sine functions while another component is the sum of all coefficients of cosine functions. Then the magnitude of the two-dimensional vector is used to measure periodicity of time-course gene expression profile. The rather subjective cut-off value is taken to determine if a gene is periodically expressed. By this way, Spellman et al. determine that 800 genes are periodically expressed out of more 6000 gene expression profiles from yeast Saccharomyces cerevisiae. After performing cluster analysis, these 800 genes are divided into five groups [2]. However, microarray experiments typically generate short time-course data. As pointed in [5, 6], the frequency resolution obtained on such short time-course data by the discrete Fourier transform is often not adequate for resolving periodicities of interest.

Authors in [7] propose a method called CORRCOS to find periodically expressed genes. CORRCOS generates totally 101000 periodic synthetic models. Each gene expression profile is compared to each of these 101000 models. Although it can identify periodically expressed gene, CORRCOS is too time consuming and the cross-correlation is not real metric. In [6], authors develop another algorithm named RAGE for detecting periodically expressed genes. Like CORRCOS, RAGE is a synthetic model-based method. Compared with CORRCOS, RAGE is less time consuming [6]. Wichert et al. [8] propose a statistical method to identify periodically expressed genes from their time-course gene expression profiles. The method models gene expression profiles also as sine functions use the Fisher 𝑔-test for statistical analysis. Given a time-course gene expression profile 𝑦_𝑡(𝑡=1,2,…,𝑚), the 𝑔-static is defined as 𝑔=max^𝑘𝐼𝜔_𝑘²∑^{[𝑚/2]𝑘=1}𝐼𝜔_𝑘,(1.1) where 𝐼(𝜔)=1𝑚|||||^𝑚_𝑡=1𝑦_𝑡exp(−𝑖𝑤𝑡)|||||²,𝜔∈[0,𝜋](1.2) is called the periodogram. It is assumed that if a time-course gene expression profile has a significant sinusoidal component with frequency 𝜔₀∈[0,𝜋], the periodogram exhibits a peak at that frequency with a high probability. On the other hand, if a time-course gene expression profile is purely random, the periodogram reduces to a straight line. Based on Fisher 𝑔-test [9], Chen [10] proposes a C&G procedure to identify periodically expressed genes from their time-course expression profiles. The 𝑔-statistic is effective only for evenly spaced gene expression profiles. For unevenly spaced gene expression profiles, Chen et al. propose to use Lomb-Scargle periodograms to discover statistically significant periodic gene expression [11, 12]. However, a recent research [13] has concluded that the Fisher 𝑔-test is poor if the time-course data is short and/or that data length is not an integer number of periods. Therefore, one can not expect to get a good clustering based on periodically expressed genes identified from these methods.

On the other hand, a number of clustering methods have been proposed for cluster analysis on gene expression data. These include distance/correlation-based clustering methods (e.g., hierarchical clustering [14], 𝑘-means clustering [15], and self-organizing maps [16]) and static model-based clustering methods [17, 18]. In these methods, gene expression profiles are viewed as multidimensional vectors. Distance/correlation-based clustering methods cluster genes based on the distance/correlation among their expression profiles. Static model-based clustering methods assign genes to one of clusters if their expression profiles may be generated by a multivariate normal distribution. These methods do not take into account the dynamic of time-course gene expression data and thus are not efficient for periodically expressed gene data.

Recently, some dynamic model-based clustering methods have been proposed to analyze time-course gene expression data [19, 20]. These methods employ autoregressive models to describe the dynamics of time-course gene expression data. As periodically expressed genes are associated with periodic biological processes, it is natural to model a periodically expressed gene data by periodic (nonlinear) function. This paper proposes a nonlinear model based method for clustering periodically expressed genes from their time-course expression profiles. The proposed method assumes that a periodically expressed gene dataset is generated by a number of periodical processes which are modelled by a linear combination of trigonometric sine and cosine functions in time plus a Gaussian noise term. A two-stage method is proposed to estimate the model parameters, and a relocation-iteration algorithm is employed to assign each gene to an appropriate cluster. A bootstrapping method and an average adjusted Rand index (AARI) are employed to measure the quality of clustering. One synthetic dataset and two biological datasets were employed to evaluate the performance of the proposed method.

2. METHODS

2.1. Model for Periodically Expressed Gene Profiles

Let 𝑥(𝑡) (𝑡=1,2,…,𝑚) be a time-course gene expression profile generated from a periodical biological process, where m is the number of time points at which gene expression is measured. After shifting the mean of gene expression profiles to 0, the periodicity of this time-course gene expression profile can be modeled by a linear combination of trigonometric sine and cosine functions in time plus a Gaussian noise term as follows [21] 𝑥(𝑡)=𝑎cos(𝜔𝑡)+𝑏sin(𝜔𝑡)+𝜀(𝑡),(2.1) where 𝑎 and 𝑏 are the coefficients of sine and cosine function, respectively; 𝜔 is the frequency of periodic expression data, and 𝜀(𝑡) represent random errors. This study assumes that the errors have a normal distribution independent of time with the mean of 0 and the variance of 𝜎². This model is equivalent to sinusoidal function model [7, 8, 10–13] 𝑥(𝑡)=𝐴sin(𝜔𝑡+Φ)+𝜀(𝑡)(2.2) which are widely used to generate the synthetic periodic gene expression profiles [7] and to detect the periodically expressed genes [2, 8, 10–12]. In model (2.2), 𝐴=√𝑎²+𝑏² is called magnitude and Φ=arctan(𝑎/𝑏)is called the phase.

Given a time-course gene expression profile 𝑥(𝑡)(𝑡=1,2,…,𝑚), estimating parameters 𝑎, 𝑏, and 𝜔 in model (2.1) is a nonlinear estimation problem as 𝜔 is nonlinear in the model. In general, all nonlinear optimization programs can be used to estimate parameters in model (2.1), for example, Gauss-Newton iteration method and its variants such as Box-Kanemasu interpolation method, Levenberg damped least squares methods, and Marquardt’s method [22]. However, these iteration methods are sensitive to initial values. Another main shortcoming is that these methods may converge to the local minimum of the least squares cost function and thus cannot find the true values of the parameters.

Our observation is that noise-free model  (2.1) 𝑥(𝑡)=𝑎cos(𝜔𝑡)+𝑏sin(𝜔𝑡)(2.3) can be viewed as the general solution of a following second-order ordinary differential equation 𝑥̈(𝑡)+𝜔²𝑥(𝑡)=0(2.4) and that 𝜔² is linear in equation (2.4) which is independent of 𝑎 and 𝑏. Therefore, we propose the following two-step parameter estimation methods to estimate parameters 𝑎, 𝑏, and 𝜔 in model  (2.2).

Step 1. Numerically calculate the second derivative of 𝑥(𝑡). Then, based on equation (2.4), use linear least squares method to estimate parameter 𝜔². In details, let 𝑋2=[𝑥̈(1),…,𝑥̈(𝑙)]^𝑇,𝑋1=[𝑥(1),…,𝑥(𝑙)]^𝑇,(2.5) then, by the least squares method, 𝜔² is estimated as 𝜔²=−𝑋1^𝑇𝑋2𝑋1^𝑇𝑋1,𝜔=𝜔²,(2.6) as time-course gene expression data are discrete, the second derivative 𝑥̈(𝑡) is estimated by the central finite difference formula as follows: 𝑥̈(𝑡)=𝑥(𝑡+1)+𝑥(𝑡−1)−2𝑥(𝑡)Δ²for𝑡=2,…,𝑚−1,(2.7) where Δ is time difference between two consecutive gene expression data points. From (2.7), the length of vectors 𝑋2 and 𝑋1 is 𝑚−2. Note that if the value of 𝜔² calculated by (2.6) for a gene is negative, this gene will be judged not to be periodically expressed.

Step 2. Substitute the estimated value of 𝜔 into (2.2). Apply the maximum likelihood method to model (2.1) to estimate parameters 𝑎 and 𝑏. In detail, let 𝑋=[𝑥(1),…,𝑥(𝑚)],𝐴=⎡⎢⎢⎣csoinsΔΔ𝜔𝜔,,……,,csions𝑚𝑚ΔΔ𝜔𝜔⎤⎥⎥⎦(2.8) by the least squares method, 𝑎 and 𝑏 are estimated as ⎡⎢⎢⎣𝑎𝑏̂⎤⎥⎥⎦̂=𝐴𝐴_𝑇⁻¹𝐴𝑋^𝑇.(2.9)

2.2. Nonlinear Model-Based Clustering

2.2.1. The Mixture Model

In this study, it is assumed that a time-course gene expression dataset is a collection of periodically expressed gene profiles which belong to several clusters, and profiles in each cluster can be described by model (2.1) or (2.2) with different parameters. Let 𝜃_𝑘=[𝑎_𝑘,𝑏_𝑘,𝜔_𝑘,𝜎^2𝑘] be parameters of model (2.1) for the 𝑘th cluster. Then the task of nonlinear model-based clustering is as follows: for a given number of cluster 𝐾, divide a time-course gene expression dataset into a partition 𝐶={𝐶₁,…,𝐶_𝑘,…,𝐶_𝐾} using model (2.1) with parameters 𝜃_𝑘=[𝑎_𝑘,𝑏_𝑘,𝜔_𝑘,𝜎^2𝑘] (𝑘=1,…,𝐾) which minimize 𝑓(Θ)=^𝐾_𝑘=1_{𝑥∈𝐶𝑘}^𝑚_𝑖=1𝑥(𝑖)−𝑎_𝑘cos𝑖Δ𝜔_𝑘−𝑏_𝑘sin𝑖Δ𝜔_𝑘²,(2.10) where the parameters Θ consist of {𝜃_𝑘,𝑘=1,…,𝐾}.

2.2.2. Estimation of Model Parameters

According to the parameter estimation method proposed in previous section for a single time-course expression profile, for the 𝑘th cluster parameters, 𝜃_𝑘=[𝑎_𝑘,𝑏_𝑘,𝜔_𝑘,𝜎^2𝑘]can be estimated as 𝜔^2𝑘=−∑^{𝑥∈𝐶𝑘}𝑋1^𝑇𝑋2∑_{𝑥∈𝐶𝑘}𝑋1^𝑇𝑋1,𝜔^𝑘⎡⎢⎢⎣𝑎_𝑘̂=𝜔^2𝑘,𝑏̂_𝑘⎤⎥⎥⎦=1𝜎_2𝑘||𝐶_𝑘||_{𝑥∈𝐶𝑘}𝐴𝐴^𝑇⁻¹𝐴𝑋^𝑇,=1𝑚||𝐶_𝑘||_{𝑥∈𝐶𝑘}^𝑚_𝑖=1𝑥(𝑖)−𝑎_𝑘̂cos(𝑖Δ𝜔_𝑘)−𝑏_𝑘sin(𝑖Δ𝜔_𝑘)²,(2.11) where |𝐶_𝑘| represents the number of time series in cluster 𝐶_𝑘, ∑^𝐾𝑘=1|𝐶_𝑘|=𝑁.

2.2.3. Algorithm

This study employs a relocation-iteration algorithm as shown in Algorithm 1 to estimate the parameters such that the cost function (2.10) is minimized. In 2(a) of Algorithm 1, Θ^𝑡 represents the estimated parameters in cost function (2.10) at iteration 𝑡 while, in 2(b), parameters𝑎^𝑡𝑘,𝑏̂^𝑡𝑘, and 𝜔^𝑡𝑘represent the parameters of model 𝑘 at iteration 𝑡.

(1) Select an initial partition for given the number of clusters, 𝐾 (2) Iteration (𝑡=1,2,…):  (a) estimate the parameter Θ𝑡 based on the present partition by using (2.11);  (b) generate a new partition by assigning each sequence 𝑥 to cluster 𝑘 for which the value of 𝑠2=∑𝑚𝑖=1[𝑥(𝑖)−𝑎𝑡𝑘̂cos(𝑖Δ𝜔𝑡𝑘)−𝑏̂𝑡𝑘sin(𝑖Δ𝜔𝑡𝑘)]2 is minimum. (3) Stop if the improvement of the cost function (2.10) is below a given threshold, the cluster memberships of time series do not change.

Algorithm 1 

Algorithm for nonlinear model-based clustering.

2.3. Evaluation

In this study, we use the adjusted Rand index (ARI) [23] to evaluate the quality of the clustering. Consider two partitions of 𝑁 objects: the 𝑟-cluster partition 𝑈={𝑢₁,…𝑢_𝑟} and the 𝑠-cluster partition 𝑉={𝑣₁,…,𝑣_𝑠}. One may construct a contingency table (matrix) as in Table 1.

Table 1 

Contingency table for two partitions of 𝑛 objects.

In Table 1, entry 𝑛_𝑖𝑗 is the number of objects that are both in clusters 𝑢_𝑖 and 𝑣_𝑗, 𝑖=1,…,𝑟, 𝑗=1,…,𝑠. Let 𝑛_𝑖.=∑^𝑠𝑗=1𝑛_𝑖𝑗 and 𝑛_.𝑗=∑^𝑟𝑖=1𝑛_𝑖𝑗 denote the sum of row 𝑖(𝑖=1,…,𝑟) and the sum of column 𝑗 (𝑗=1,…,𝑠) in the contingency matrix, respectively, and let𝑉=𝑁(𝑁−1)/2 (the number of pairs of 𝑁 objects). Based on the contingency matrix of two partitions, the ARI is defined as [23] ARI=∑^𝑟𝑖=1∑^𝑠𝑖=1^𝑛𝑖𝑗₂−(1/𝑉)∑^𝑟𝑖=1^𝑛𝑖.₂∑^𝑠𝑖=1^𝑛.𝑗₂(1/2)∑^𝑟𝑖=1^𝑛𝑖.₂+∑^𝑠𝑖=1^𝑛.𝑗₂−(1/𝑉)∑^𝑟𝑖=1^𝑛𝑖.₂∑^𝑠𝑖=1^𝑛.𝑗₂.(2.12) The expected value of ARI is 1 when they matched perfect and 0 when the two partitions are selected at random.

If the true cluster labels for some dataset are known, the proposed clustering methods can be applied these datasets to obtain new cluster labels. Then ARI can be calculated for these two partitions. If ARI is close to 1, one can say that the proposed clustering method is in agreement with the true clusters. However, for real-life gene expression datasets, the true cluster labels are typically unknown. For this case, this study adopts a bootstrapping approach as shown in Algorithm 2 [20] to evaluate the proposed clustering methods. For the given number of clusters,𝐾, the average ARI (AARI) reports the quality of the clustering result obtained from the evaluated clustering methods. Accordingly, the larger AARI, the better the quality of the clustering is, that is, the better the performance of the clustering method is.

(1) Repeat the following B times (where B is a preset integer number).  (a) Randomly divide the original dataset into two nonoverlapping sets, a learning set 𝐿, and a test set 𝑇.  (b) Apply the evaluated method to the learning set 𝐿 to obtain a partition𝑃(∘,𝐿).  (c) Construct a predictor (classifier) 𝐶(∘,𝐿) using the cluster labels from the partition𝑃(∘,𝐿).  (d) Apply the predictor C(∘,𝐿) to the test set T to get the predicted partition 𝑃(∘,𝑇).  (e) Apply the evaluated method to the test set 𝑇 to obtain a partition 𝑃(∘,𝑇).  (f) Calculate the ARI of partitions 𝑃(∘,𝑇) and 𝑃(∘,𝑇). (2) Calculate the average ARI (AARI) over the B times as the measure index of the proposed clustering method. (3) For the various number of clusters, 𝐾, repeat the procedure described in steps (1) and (2) above to get AARI(𝐾), and then plot AARI(𝐾) with respect to 𝐾.

Algorithm 2 

The procedure for evaluating proposed clustering method.

3. EXPERIMENTAL RESULTS AND DISCUSSION

This study employs a synthetic dataset and two biological datasets to investigate the performance of the proposed method in different aspects.

3.1. Synthetic Dataset (SYN)

The synthetic dataset is generated by model (2.1). Let 𝑥_𝑖𝑡 be the simulated expression (log-ratio) values of gene 𝑖 at time point 𝑡 in the dataset, that is, 𝑥_𝑖𝑡=𝑎_𝑘cos𝑤_𝑘𝑡+𝑏_𝑘sin𝑤_𝑘𝑡+𝜀_𝑖𝑡,𝑖=1,2,…,𝑛,𝑡=1,2,…,𝑚;𝑘=1,2,…,𝐾,(3.1) where 𝑛 is the number of genes, 𝑚 is the number of time points, and 𝐾 is the number of clusters.

In this study, parameters for synthetic data 𝑎_𝑘 , 𝑏_𝑘 , and 𝑤_𝑘 are randomly chosen as follows: 𝑎_𝑘∼𝑁3𝑘,𝑑²,𝑏_𝑘∼𝑁3𝑘,𝑑²,𝑤_𝑘∼𝑈2𝜋𝐾=5,𝑛=20,𝑛_𝑘𝑘+1,2𝑘𝜋,𝜀∼𝑁0,𝜎²,∼𝑈(100,300),(3.2) where 𝑛_𝑘 is the number genes in the 𝑘th cluster. The resulted parameters for synthetic data are shown in Table 2.

Table 2 

Parameters for synthetic data.

For various numbers of clusters, we run the proposed method described in Algorithm 1 with randomly chosen initial partitions, with the initial partitions from 𝑘-means results as and to the 𝑘-means methods. The ARI between clustering results and the known true cluster labels is calculated. The values of AARI are calculated over 20 runs and shown in the Table 3 and Figure 1.

Table 3 

The values of AARI for different clustering methods on synthetic data.

Figure 1 

Plot of ARRI over different number of clusters for synthetic dataset.

From Figure 1, the proposed method with both initial partitions randomly chosen and those from 𝑘-means results has greater value of AARI than 𝑘-means when the number of clusters is greater than 3. Furthermore, when the number of clusters is the true value of 5, the AARI of the proposed method with both initial partitions reaches its maximum, which makes sense. However, the AARI of 𝑘-means method did not reach its maximum when the number of clusters is 5. Therefore, we can conclude that the proposed method outperforms the 𝑘-means in terms of AARI.

3.2. Real-Life Datasets

In this study, two real-life datasets are employed to illustrate the proposed method: ELU and BAC. ELU consist of expression profiles of 4304 genes without missing data. Expression profiles are obtained from yeast cell cycle division process through Eluration-synchronized experiments conducted by Spellman et al. [2]. Each expression profile has 14 equally spacing time points. BAC consists of expression profiles of 1590 genes without missing data. Expression profiles are measured during the cell cycle division process of the bacterium Caulobacter crescentus [3]. The measurements were taken at 11 equally spaced time points over 150 minutes. Both datasets are preprocessed in the following two steps.

Step 1. Shift the mean of each gene expression profile to 0.

Step 2. Filter the dataset with 𝐹-test at the significance level 𝛼, that is, 𝐻₀∶𝑥_𝑖(𝑡)=𝐹𝜀=_𝑖𝑡,𝑅_𝐻0𝐻₁∶𝑥_𝑖𝑡=𝑎cos(𝑤𝑡)+𝑏sin(𝑤𝑡)+𝜀_𝑡,−𝑅_𝐻1/2𝑅_𝐻1/(𝑚−2)∼𝐹(2,𝑚−2),(3.3) where 𝑅_𝐻𝑖 is the sum of squared errors under the specific hypothesis and 𝑚 is the number of time points. Keep the genes which reject the null hypothesis (show periodical behaviours) [21].

After these two steps, the number of genes remains for different significant level as in Table 4. Then we run the evaluation procedure proposed in Algorithm 2 on these selected gene expression profiles. The AARIs of the proposed method and 𝑘-means over various numbers of clusters are plotted in Figures 2 and 3 for dataset ELU and BAC, respectively. From Figures 2 and 3 the results from both real-life datasets show that the proposed method outperforms the 𝑘-means in terms of AARI.

Table 4 

The number of genes after filtering.

(a)

(b)

Figure 2 

Plot of ARRI over different number of clusters for ELU after filtering.

(a)

(b)

Figure 3 

Plot of ARRI over different number of clusters for BAC after filtering.

4. CONCLUSIONS

This paper has presented a nonlinear model-based method for clustering periodically expressed genes from their time-course expression profiles. In this method, profiles of periodically expressed genes and thus the cluster of profiles are modelled by a linear combination of trigonometric sine and cosine functions in time plus a Gaussian noise term which is equivalent to a sinusoidal function model [1–4, 6–13, 17–19]. Although this model is not new, the existing methods are not based on parameter estimation technique, especially not estimating the frequency in the model as it is nonlinear in parameter. In the presented method, a two step linear least squares method is proposed to estimate all model parameters including the frequency for each clusters. Computational experiments on one synthetic dataset and two biological datasets show that the proposed method outperforms the traditional clustering methods such as 𝑘-means in terms of AARI, which indicate that the proposed method can effectively cluster periodically expressed genes from their time-course expression profiles.

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

ACKNOWLEDGMENTS

This research is supported by Science and Technology Funds of Beijing Ministry of Education (SQKM201210037001) through the first author and Natural Sciences and Engineering Research Council of Canada (NSERC) through other authors.