# A method to identify differential expression profiles of time-course gene data with Fourier transformation

- Jaehee Kim
^{1}Email author, - Robert Todd Ogden
^{2}and - Haseong Kim
^{3}

**14**:310

**DOI: **10.1186/1471-2105-14-310

© Kim et al.; licensee BioMed Central Ltd. 2013

**Received: **15 February 2013

**Accepted: **10 October 2013

**Published: **18 October 2013

## Abstract

### Background

Time course gene expression experiments are an increasingly popular method for exploring biological processes. Temporal gene expression profiles provide an important characterization of gene function, as biological systems are both developmental and dynamic. With such data it is possible to study gene expression changes over time and thereby to detect differential genes. Much of the early work on analyzing time series expression data relied on methods developed originally for static data and thus there is a need for improved methodology. Since time series expression is a temporal process, its unique features such as autocorrelation between successive points should be incorporated into the analysis.

### Results

This work aims to identify genes that show different gene expression profiles across time. We propose a statistical procedure to discover gene groups with similar profiles using a nonparametric representation that accounts for the autocorrelation in the data. In particular, we first represent each profile in terms of a Fourier basis, and then we screen out genes that are not differentially expressed based on the Fourier coefficients. Finally, we cluster the remaining gene profiles using a model-based approach in the Fourier domain. We evaluate the screening results in terms of sensitivity, specificity, FDR and FNR, compare with the Gaussian process regression screening in a simulation study and illustrate the results by application to yeast cell-cycle microarray expression data with alpha-factor synchronization.

The key elements of the proposed methodology: (i) representation of gene profiles in the Fourier domain; (ii) automatic screening of genes based on the Fourier coefficients and taking into account autocorrelation in the data, while controlling the false discovery rate (FDR); (iii) model-based clustering of the remaining gene profiles.

### Conclusions

Using this method, we identified a set of cell-cycle-regulated time-course yeast genes. The proposed method is general and can be potentially used to identify genes which have the same patterns or biological processes, and help facing the present and forthcoming challenges of data analysis in functional genomics.

## Background

Time-course gene expression data are often measured to study dynamic biological systems and gene regulatory networks. Array technologies have made it straightforward to monitor the expression pattern of thousands of genes simultaneously. The challenge now is to interpret such massive data sets. The first step is to extract the fundamental patterns of gene expression inherent in the data. Gene-expression levels can be monitored with cDNA or oligonucleotide chips over a time-course for a temporal process. Following a microarray time series experiment, a key challenge is to extract the continuous representation of all genes throughout the time-course. Identifying significant or differentially expressed genes is challenging because different genes may have different profiles, and because of the noise present in time series expression data. A comprehensive review about time series expression data analysis and the related computational challenges may be found in [1].

Microarrays have recently been used for the purpose of monitoring expression levels of thousands of genes simultaneously and for identifying genes that are differentially expressed. With the number of inferences made in the analysis of microarray data, it is natural to be concerned about multiple testing. This problem of multiplicity can be dealt with by controlling the false discovery rate (FDR) [2].

In the past decade, many approaches to gene selection have been considered including a two sample t-test [3], a regression approach [4], and a mixture model approach [5]. Other approaches to this problem include the Empirical Bayesian (EB) method [6] and the Significance Analysis of Microarray (SAM) method [7]. The multiplicity problem is addressed in adopting a resampling-based approach to controlling FDR [8]. Also an ANOVA formulation and an empirical Bayes adjustment to the t-statistics [9] and an empirical Bayes screening procedure have been proposed [10].

There has been considerable research about discovering patterns using clustering and testing including clustering after transformation and smoothing as a technique for nonparametrically estimating and clustering a large number of curves [11] and clustering short time series gene expression data by selecting a set of potential expression profiles [12].

Smoothing away noise-induced wiggles of gene expression data with Fourier series for microarray data has been considered including an improved Fourier method with irregular or monotonic components of cell-cycle expression [13], a two-step procedure for clustering periodic patterns of gene expression profiles using a Fourier series approximation with frequency and amplitude of order one [14], a multivariate modeling approach using partial least squares (PLS) regression to identify genes with periodic fluctuations in the budding yeast cell cycle data [15], a Hidden Markov Models (HMMs) approach to account for the horizontal dependencies along the time axis [16], and a model-based clustering of the Fourier coefficients calculated on the first difference of the time-course data [17].

Model-based hierarchical clustering was proposed in character recognition problems using a multivariate normal model [18] and it may be used to guide the choice of the model based on computing an approximate maximum for the classification likelihood [19].

There has been much work done on clustering microarray data, mostly on grouping common expression patterns. However, less attention has been paid to time-course gene studies. Currently the analysis of GETS (gene expression time-series) is commonly performed using a GP (Gaussian process) [20-24]. Also a Bayesian analysis of microarray time series has been developed with the software package BATS [25].

In this research, we propose a new method for gene screening using Fourier coefficients to cluster time-courses of genes that exhibit similar patterns.

This paper introduces a methodology for gene selection based on time-course data. The first step is screening, in which we seek to isolate the inactive genes from the active ones, while properly taking into account the serial dependence in the time course data and controlling the FDR, all in the Fourier domain. The second step involves a model-based clustering of the “active” genes, also in the Fourier domain. We evaluate the performance of the methodology using both simulated data and yeast cell cycle data.

## Results and discussion

### Simulated data

Since real expression data sets are generally noisy and their clusters may not be fully reflective of the class information, we first evaluate the performance of our method with simulated data, for which the “true” classes are known.

*n*= 800 genes,

*m*= 20 time points, and

*E*(

*ϵ*

_{ iu }) = 0 and

*ϵ*

_{ iu }’s from an autoregressive AR(1) process with a variety of values of the AR parameter. The regression functions for

*f*are:

Each simulated dataset consists of 800 curves originating from the 6 functions: 400 *f*_{1}, s and 80 curves of each *f*_{2}, ⋯, *f*_{6}, to reflect typical gene expression data. Thus, there are 5 sets of differential genes and 1 set of non-differential genes. The standard deviation of the innovation process was set to *σ* = 0.5 and *σ* = 1.5 to represent low and high noise situations, respectively.

The cosine system $\left\{\sqrt{2}cos(\pi \mathit{\text{jt}})\right\}$ is used as the set of basis functions. Though the optimal choice for *J* could vary from function to function, we choose to use a single smoothing parameter that performs reasonably well for all of the curves. In the simulation ten Fourier coefficients are used for the spectrum estimation. Several numbers of Fourier coefficients are considered for test statistics in the proposed screening method.

The number of clusters is determined according to the Bayesian Information Criterion (BIC).

*T*be a clustering map defined as

*η*(

*K*) depending on

*K*clusters is defined as

where {*f*_{1}, ⋯, *f*_{
n
}} denote the true curves and $\left\{{\widehat{f}}_{1},\cdots ,{\widehat{f}}_{n}\right\}$ denote the estimated curves. Let *T* and $\widehat{T}$ represent the corresponding cluster maps, and *K* denote the number of clusters. *n(K)* is the fraction of all pairs that are incorrectly placed in different clusters.

We compare our screening method with the recently proposed GPR (Gaussian process regression) screening [24]. GPR seeks to quantify the true signal and noise in a gene expression time series, allows us to rank the differential expression of the gene profile. A Gaussian process using a squared-exponential covariance function is based on the assumption that the underlying true signal in a profile is a smooth function. GPR is applied to each gene curve with the assumption that each gene can be categorized as either quiet or differentially expressed. The genes may then be ranked in decreasing order according to their values of likelihood ratio (LR) test statistics for testing between these two characterizations. The optimal critical value for the LR test statistics for GPR screening was empirically determined to be log(1.5) based on the true positive rate (sensitivity) since in the simulation, this cutoff allowed more than 80% of differentially expressed genes to be detected. There is no clear-cut critical value for this scheme. The genes greater than the critical value are detected as differential. After the differential genes are obtained with GPR screening, clustering procedure can be done with those genes.

- i.
Sensitivity: proportion of differentially expressed genes that are declared significant

- ii.
Specificity: proportion of non-differentially expressed genes that are not declared significant

- iii.
False discovery rate (FDR): proportion of genes declared significant that are not differentially expressed

- iv.
False non-discovery rate (FNR): proportion of genes not declared significant that are differentially expressed

**Comparison of screening and clustering results (low noise)**

Without screening | With screening | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

AR(1) parameter | Method |
| Error | Sil | ARI | Error | Sil | ARI | Sensitivity | Specificity | FDR | FNR |

| FC* | 2 | .037 | .509 | .909 | .015 | .560 | .918 | .878 | .723 | .121 | .276 |

3 | .020 | .471 | .921 | .016 | .484 | .932 | .860 | .783 | .139 | .216 | ||

4 | .015 | .430 | .963 | .017 | .438 | .937 | .863 | .842 | .136 | .157 | ||

5 | .015 | .388 | .964 | .014 | .403 | .944 | .854 | .798 | .145 | .201 | ||

8 | .015 | .305 | .964 | .017 | .317 | .940 | .851 | .836 | .148 | .163 | ||

GPR** | .855 | .779 | .220 | .144 | ||||||||

| FC | 2 | .052 | .471 | .871 | .021 | .523 | .875 | .871 | .722 | .128 | .277 |

3 | .036 | .423 | .912 | .026 | .443 | .888 | .846 | .783 | .153 | .217 | ||

4 | .029 | .386 | .931 | .028 | .398 | .895 | .847 | .839 | .152 | .160 | ||

5 | .027 | .348 | .935 | .022 | .366 | .906 | .837 | .798 | .162 | .205 | ||

8 | .028 | .274 | .936 | .029 | .287 | .895 | .830 | .836 | .169 | .163 | ||

GPR | .826 | 678 | .321 | .173 | ||||||||

| FC | 2 | .073 | .430 | .822 | .030 | .487 | .815 | .863 | .723 | .136 | .276 |

3 | .056 | .380 | .865 | .042 | .402 | .814 | .828 | .783 | .171 | .217 | ||

4 | .052 | .339 | .875 | .045 | .356 | .825 | .827 | .834 | .172 | .165 | ||

5 | .049 | .306 | .883 | .036 | .326 | .845 | .817 | .790 | .182 | .209 | ||

8 | .047 | .244 | .888 | .049 | .257 | .823 | .803 | .832 | .196 | .167 | ||

GPR | .798 | .571 | .428 | .201 | ||||||||

| FC | 2 | .159 | .340 | .610 | .056 | .414 | .633 | .835 | .717 | .165 | .201 |

3 | .139 | .287 | .663 | .093 | .329 | .591 | .775 | .768 | .224 | .231 | ||

4 | .124 | .255 | .702 | .113 | .280 | .578 | .766 | .811 | .233 | .188 | ||

5 | .132 | .226 | .682 | .093 | .259 | .615 | .762 | .773 | .237 | .226 | ||

8 | .143 | .181 | .649 | .134 | .205 | .562 | .730 | .815 | .269 | .184 | ||

GPR | .756 | .410 | .589 | .244 | ||||||||

| FC | 2 | .266 | .287 | .345 | .088 | .357 | .351 | .755 | .704 | .244 | .295 |

3 | .264 | .224 | .347 | .153 | .272 | .314 | .682 | .738 | .317 | .261 | ||

4 | .258 | .190 | .370 | .186 | .230 | .303 | .668 | .771 | .331 | .228 | ||

5 | .258 | .171 | .375 | .161 | .211 | .317 | .676 | .745 | .324 | .255 | ||

8 | .267 | .137 | .339 | .220 | .172 | .287 | .641 | .769 | .358 | .230 | ||

GPR | .731 | .335 | .664 | .268 |

**Comparison of screening and clustering results (high noise)**

Without screening | With screening | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

AR(1) parameter | Method |
| Error | Sil | ARI | Error | Sil | ARI | Sensitivity | Specificity | FDR | FNR |

| FC* | 2 | .326 | .259 | .200 | .235 | .292 | .149 | .589 | .708 | .411 | .291 |

3 | .321 | .192 | .199 | .298 | .214 | .151 | .561 | .716 | .439 | .283 | ||

4 | .321 | .151 | .194 | .320 | .166 | .156 | .553 | .722 | .447 | .278 | ||

5 | .323 | .125 | .185 | .269 | .142 | .156 | .571 | .714 | .428 | .285 | ||

8 | .324 | .084 | .175 | .324 | .094 | .148 | .552 | .722 | .448 | .277 | ||

GPR** | .483 | .779 | .221 | .517 | ||||||||

| FC | 2 | .343 | .253 | .164 | .234 | .291 | .117 | .567 | .697 | .432 | .302 |

3 | .339 | .185 | .155 | .287 | .208 | .117 | .545 | .703 | .454 | .296 | ||

4 | .338 | .149 | .151 | .306 | .160 | .120 | .539 | .708 | .461 | .291 | ||

5 | .337 | .125 | .146 | .261 | .138 | .120 | .555 | .702 | .445 | .297 | ||

8 | .338 | .086 | .132 | .307 | .094 | .113 | .538 | .706 | .462 | .293 | ||

GPR | .536 | .677 | .323 | .463 | ||||||||

| FC | 2 | .359 | .248 | .128 | .284 | .290 | .090 | .546 | .681 | .453 | .318 |

3 | .351 | .185 | .119 | .329 | .208 | .089 | .531 | .683 | .468 | .316 | ||

4 | .350 | .148 | .115 | .347 | .159 | .092 | .526 | .685 | .473 | .314 | ||

5 | .350 | .127 | .108 | .304 | .137 | .088 | .537 | .681 | .462 | .318 | ||

8 | .357 | .083 | .089 | .351 | .091 | .079 | .526 | .685 | .474 | .314 | ||

GPR | .584 | .572 | .427 | .415 | ||||||||

| FC | 2 | .383 | .246 | .073 | .330 | .284 | .053 | .517 | .632 | .482 | .367 |

3 | .375 | .183 | .066 | .356 | .198 | .051 | .512 | .634 | .488 | .365 | ||

4 | .369 | .151 | .062 | .365 | .158 | .052 | .510 | .633 | .490 | .366 | ||

5 | .369 | .126 | .056 | .338 | .137 | .046 | .514 | .633 | .485 | .367 | ||

8 | .370 | .086 | .046 | .370 | .092 | .042 | .509 | .634 | .490 | .365 | ||

GPR | .646 | .409 | .590 | .353 | ||||||||

| FC | 2 | .395 | .248 | .035 | .356 | .275 | .030 | .505 | .504 | .495 | .422 |

3 | .384 | .186 | .034 | .368 | .193 | .028 | .503 | .503 | .496 | .424 | ||

4 | .383 | .148 | .031 | .373 | .155 | .026 | .502 | .502 | .497 | .421 | ||

5 | .381 | .125 | .028 | .358 | .134 | .024 | .504 | .504 | .496 | .419 | ||

8 | .377 | .092 | .027 | .370 | .097 | .023 | .503 | .502 | .497 | .415 | ||

GPR | .679 | .337 | .662 | .320 |

As shown in Table 2 the error rates of the high noise case are bigger and the sensitivity is smaller than those of the lower noise case in Table 1. However the high noise level does not diminish the performance of the screening procedure. From Tables 1 and 2 it can be seen that the clustering estimation error is smaller after screening than it is without screening. Also, the clustering estimation error becomes smaller as the number of Fourier coefficients *J* becomes larger and for smaller values of the AR parameter. Table 2 (higher noise level) shows that the screening percentage of our method is a little higher than GPR for some choices of *J*. Overall, the proposed method demonstrates improved sensitivity according to the number of Fourier coefficients *J* and much improved specificity and FDR as compared to GPR. In addition, the proposed method has an advantage that it does not require estimation of the covariance structure

### Yeast cell cycle data analysis

We have used alpha synchronized yeast cell expression data [26] available at http://genome-www.stanford.edu/cellcycle/ to test our algorithm. After removing genes with missing values, there were 4,489 genes remaining out of 6178 genes. This data contained 18 time points sampled uniformly every 7 min between 0 and 119 min.

*J*. Table 3 shows the median and the average silhouette values with Euclidean distance between samples by model-based clusterings for various

*J*values both with and without the screening step using a significance level of FDR 5%. The number of genes with differentially expression is 2,227 out of 4,489 at the significance level

*α = 0.05*.

**Silhouette values for model-based clustering with Fourier coefficients of yeast data**

Number. of Fourier coeff | Number of clusters without-screening | Number of clusters with screening | No. of sig. genes | Without-screening | Screening | ||
---|---|---|---|---|---|---|---|

Med sil | Avg. sil | Med sil | Avg. sil | ||||

J = 2 | 5 | 4 | 1715 | .160 | .112 | .451 | .388 |

J = 3 | 5 | 4 | 1735 | .204 | .146 | .505 | .426 |

J = 4 | 5 | 4 | 2227 | .174 | .119 | .552 | .485 |

J = 5 | 5 | 4 | 2792 | .029 | .015 | .048 | .028 |

J = 6 | 5 | 4 | 3071 | .196 | .119 | .041 | .043 |

J = 8 | 5 | 4 | 3050 | .136 | .055 | .032 | .024 |

J = 10 | 5 | 4 | 3142 | .153 | .092 | .031 | .010 |

*J = 4*, the model-based clustering approach results in 4 clusters consisting of 51, 29, 2077, and 70 genes, respectively. Figure 1 shows plots of the sample Fourier coefficients. Figure 2 shows pointwise means of Fourier estimated gene scores in each cluster with

*J = 4*sample Fourier coefficients. The graph in the bottom left-hand corner of Figure 2 shows the estimated mean of gene scores of screened out. Gaussian mixture model clustering allows clusters to have different orientation or sizes while preserving some common features, such as an ellipsoidal shape. Clusters 1, 2 and 4 in particular show a cycle while cluster 3 consists of less active genes.

We also applied GPR screening for yeast data with the critical value log (1.5) and identified 1,620 non-quiet differentially expressed genes. A Gaussian mixture model clustering algorithm was applied to these genes indicating 4 clusters. The clusters have 209, 477, 598, and 336 genes, respectively. For this clustering the median silhouette value is 0.508 and the average silhouette value is 0.370, each less than the corresponding silhouette summaries for the proposed procedure.

Owing to noise and the high dimensionality of data, careful consideration of statistical and biological validity is needed when analyzing microarray data. From our review we have found that without plausible interpretation and biological validation, the number of partitions produced by numerical analysis is highly unreliable, and sometimes even misleading.

**Number of genes in each cluster with J = 4**

Screened out | Cluster 1 | Cluster 2 | Cluster 3 | Cluster 4 | |
---|---|---|---|---|---|

Number of genes | 2262 | 51 | 29 | 2077 | 70 |

Filtered genes | 2041 | 46 | 28 | 1881 | 64 |

Significant GO terms | 1 | 36 | 17 | 6 | 37 |

Even though the number of genes affects to the hypergeometric test, cluster 3 and the screened-out group have only one and six significantly overrepresented GO terms, respectively. This means that all the genes (or some of them) in those clusters rarely share their biological processes and these clusters can be considered as random collections of genes. So, we focused on the three clusters 1, 2, and 4 to find their biological meanings.

**GO terms connected sequentially in their GO relationship graph (C: cluster number, S: subset number)**

C | S | GO ID | Terms | S | GO ID | Terms |
---|---|---|---|---|---|---|

1 | 1 | GO:0006334 | Nucleosome assembly | 2 | GO:0000280 | Nuclear division |

GO:0031497 | Chromatin assembly | GO:0000087 | M phase of mitotic cell cycle | |||

GO:0006323 | DNA packaging | GO:0048285 | Organelle fission | |||

GO:0034728 | Nucleosome organization | GO:0022402 | Cell cycle process | |||

GO:0006333 | Chromatin assembly or disassembly | GO:0000278 | Mitotic cell cycle | |||

GO:0016043 | Cellular component organization | GO:0007049 | Cell cycle | |||

GO:0051276 | Chromosome organization | GO:0007067 | Mitosis | |||

GO:0006325 | Chromatin organization | GO:0000279 | M phase | |||

GO:0006996 | Crganelle organization | GO:0022403 | Cell cycle phase | |||

2 | 1 | GO:0007109 | Cytokinesis, completion of separation | 2 | GO:0071554 | Cell wall organization or biogenesis |

GO:0000920 | Cell separation during cytokinesis | GO:0007047 | Cellular cell wall organization | |||

GO:0032506 | Cytokinetic process | GO:0071555 | Cell wall organization | |||

GO:0000910 | Cytokinesis | GO:0070882 | Cellular cell wall organization or biogenesis | |||

GO:0031505 | Cungal-type cell wall organization | |||||

4 | 1 | GO:0051716 | Cellular response to stimulus | 2 | GO:0006302 | Couble-strand break repair |

GO:0050896 | Response to stimulus | GO:0006281 | DNA repair | |||

GO:0033554 | Cellular response to stress | GO:0000724 | Couble-strand break repair via homologous recombination | |||

GO:0006950 | Response to stress | GO:0006974 | Cesponse to DNA damage stimulus | |||

GO:0034605 | Cellular response to heat | 3 | GO:0006260 | DNA replication | ||

GO:0009408 | Response to heat | GO:0006273 | Cagging strand elongation | |||

GO:0009628 | Response to abiotic stimulus | GO:0006261 | DNA-dependent DNA replication | |||

GO:0009266 | Response to temperature stimulus | GO:0006271 | DNA strand elongation during DNA replication | |||

GO:0022616 | DNA strand elongation | |||||

GO:0006259 | DNA metabolic process |

## Conclusions

A number of recent studies in this field have focused on the analysis of time series of gene expression data, collected by performing DNA microarray experiments at several or more points in time. We have proposed a significance method to identify differentially expressed genes in time course microarray gene expression data using time series screening based on Fourier coefficients controlling FDR and model based clustering with the sample genewise Fourier coefficients, and have compared our screening method with GP screening. Recently spectral mixture kernels [32] have been introduced with a Gaussian mixture as a Fourier transform of kernels and these kernels are able to discover patterns and extrapolate and model negative covariances, illustrating the relationship between the GP and the Fourier approach.

We demonstrated the effectiveness of our approach using model-based clustering of gene profiles. Although we assumed that the residuals follow an AR process, we found that it is not necessary to assume a specific correlation structure since the sample Fourier coefficient estimates themselves do not depend heavily on the underlying covariance structure. The most commonly used techniques are clustering (unsupervised) techniques, which are particularly well suited for an exploratory investigation of this kind of data. The main advantage of the model-based methods is their reliance on a highly structured theoretical basis. Model-based clustering methods are based on the assumption that the data were generated by some underlying model and attempt to infer these models from data. Data generated by the same model is then considered to be “similar” and clustered together. Also, the choice of the optimal number of clusters and the selection of the best model can be performed using sound statistical criteria.

The proposed method is able to identify a set of cell-cycle-regulated genes in yeast and time-course genes. The proposed method is general and can be potentially used to identify genes which have the same patterns or biological processes, and help facing the present and forthcoming challenges of data analysis in functional genomics.

## Methods

### The Fourier representation model

*Y*

_{ iu },

*u*th observation on the

*i*th curve, of the form

where *E*(*ϵ*_{
iu
}) = 0 and the *ϵ*_{
iu
} values arise from a covariance-stationary process with mean zero and covariance function *γ*_{
i
}, *γ*_{
i
}(*k*) = *E*(*ϵ*_{
iu
}*ϵ*_{i,u + k}) for all *u* and *k*. In a microarray experiment *Y*_{
iu
} is the log gene expression of gene *i* at time *u*. We assume that the data from one curve are independent from those of other curves.

*f*

_{ i }belongs to a class of smooth functions as defined below:

*b*

_{ j }} is an orthonormal basis system and

*f*

_{ i }using its empirical Fourier coefficients:

which is the projection onto the first *J* basis functions where *J* , 1 ≤ *J* ≤ *m*, is a smoothing parameter to be chosen based on the data.

with *t*_{
r
} *= r/m*. To estimate the true Fourier coefficients, the covariance structure is not considered since the covariance matrix of a finite set of estimated Fourier coefficients is asymptotically proportional to the identity matrix.

### Screening out flat genes

Many microarray experiments are aimed at finding ‘active’ genes that vary significantly in expression. Differential expression indicates the changing of transcription levels across different time points, and it is thought that these transcription changes might be responsible for the change in phenotype. For example, the genes responsible for the presence of a certain disease will be transcribed at a different rate than when the disease is absent. Cluster analysis often fails to detect differentially expressed genes that belong to clusters for which most genes do not change because most of the other genes in their clusters do not change significantly.

The problem can be formulated as hypothesis-testing for individual genes as follows:

*H*_{0} : *f*_{
i
}(·) = *C* versus *H*_{1} : *f*_{
i
}(·) ≠ *C* for *i* = 1, 2, ⋯ *n*,i.e., for *n* genes, we are considering *n* pairs of mutually exclusive hypotheses:

*H*_{
0
}: Gene *i* is not differentially expressed.

*H*_{
1
}: Gene *i* is differentially expressed.

In a microarray setting, it is typical to consider thousands of tests simultaneously. In this situation the familywise error rate (FWER) or FDR (false discovery rate), the average proportion of inactive genes among those that were declared active, should be controlled. The FDR procedure [2] is as follows:

*k*be the largest

*g*, 0 ≤

*g*≤

*n*, for which

*H*

_{(g)0}, for

*g*= 1, 2, ⋯,

*k*, where

*H*

_{(g)0}is the associated null hypothesis and

*P*

_{(g)}is the

*g*th smallest p-value among all the p-values calculated for each of the hypotheses. For all genes, we apply a first-order auto-regressive (AR(1)) process to model the time dependency of the data. For testing change in the mean function of time series data, the test [32] rejects the null hypothesis of no change for large values of

*k*is

where $\widehat{\mathit{\u03f5}}$ is the residual from the Fourier estimation in (4).

*Ts*can be calculated from the asymptotic distribution. Since each $m{\widehat{\mathit{\phi}}}_{j}^{2}/\widehat{S}\left(0\right)$ has an approximate chi-squared distribution ([33, 34]),

where ${\mathit{\chi}}_{j}^{2}$ is a chi-squared random variable with *j* degrees of freedom.

### Clustering differentially expressed genes

By the Central Limit Theorem for dependent data [35], the sample Fourier coefficient ${\widehat{\mathit{\phi}}}_{j}$ is asymptotically normally distributed [36] as *m* → *∞*. With this asymptotic property, we can use the Gaussian mixture model for clustering.

With the large number of genes monitored in these studies, clustering is a key task for microarray data analysis. It seeks to identify groups of genes with similar expression profiles across samples. Clustering can reduce the effort of studying individual genes and more importantly it can unmask the functional groups among genes.

Considering that the empirical Fourier coefficients of the gene profiles have an asymptotic multivariate normal distribution enables the use of an efficient algorithm to compute the posterior probability that a gene belongs to a certain cluster.

The geometric features (shape, volume, orientation) of each group *k* are determined by the covariance matrix *∑*_{
k
} of the Fourier coefficients. A general framework for exploiting the representation of the covariance matrix is done in terms of its eigenvalue composition [37]. Each elliptical model for the covariance matrix is implemented in Mclust [19]. Model-based hierarchical agglomerative clustering is an approach to compute an approximate maximum of the classification likelihood. Each component is weighted by the probability that a sample Fourier coefficient belongs to that component. Our clustering strategy involves model-based agglomerative hierarchical clustering and selection of the model and the number of clusters using approximate Bayes factors with the BIC approximation.

### Performance metrics

For evaluating the performance of clustering algorithms, the adjusted Rand Index [38] and the Silhouette index [39] are used for the simulated data and for the yeast data.

*T*is the true clustering of a gene expression data set based on domain knowledge and

*C*a clustering result given by some clustering algorithm applied to the observed data. Let

*a*,

*b*,

*c*and

*d*respectively denote the number of gene pairs belonging to the same cluster in both

*T*and

*C*, the number of pairs belonging to the same cluster in

*T*but to different clusters in

*C*, the number of pairs belonging to different clusters in

*T*but to the same cluster in

*C*and the number of pairs belonging to different clusters in both

*T*and

*C*. The adjusted Rand index

*ARI(T,C)*is defined as.

*i*th sample in the

*j*th cluster is defined as.

where *a* (*i*) is the average distance between the *i*th sample and all other samples included in the *j*th cluster, *b*(*i*) is the minimum average distance between the *i*th sample and all of the samples clustered in *k*th cluster for *k ≠ j*

## Declarations

### Acknowledgements

JK’s research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF 2012-0001879). HK’s research was supported by the KRIBB Research Initiative Program (2012-0111).

## Authors’ Affiliations

## References

- Bar-Joseph Z: Analyzing time series gene expression data. Bioinformatics. 2004, 20 (16): 2493-2503. 10.1093/bioinformatics/bth283.View ArticlePubMedGoogle Scholar
- Benjamin Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Statist Soc B. 1995, 57: 289-300.Google Scholar
- Devore JL, Peck R: Statistics: the exploration and analysis of data. 1997, Belmont: Duxbury PressGoogle Scholar
- Thomas JG, Olson JM, Tapscott SJ, Zhao LP: An efficient and robust statistical modeling approach to discover differentially expressed genes using genomic expression profiles. Genome Res. 2001, 11: 1227-1236. 10.1101/gr.165101.PubMed CentralView ArticlePubMedGoogle Scholar
- Pan W, Lin J, Le C: A mixture model approach to detecting differentially expressed genes with microarray data. 2001, Technical report: Division of Biostatistics, University of MinnesotaGoogle Scholar
- Efron B, Tibshirani R, Gross V, Tusher VG: Empirical Bayes analysis of a microarray experiment. J Am Stat Assoc. 2001, 96: 1151-1160. 10.1198/016214501753382129.View ArticleGoogle Scholar
- Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci USA. 2001, 98: 5116-5121. 10.1073/pnas.091062498.PubMed CentralView ArticlePubMedGoogle Scholar
- Reiner A, Yekutieli D, Benjamini Y: Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinformatics. 2003, 19 (3): 368-375. 10.1093/bioinformatics/btf877.View ArticlePubMedGoogle Scholar
- Datta S, Satten GA, Benos DJ, Xia J, Heslin MJ, Datta S: An empirical bayes adjustment to increase the sensitivity of detecting differentially expressed genes in microarray experiments. Bioinformatics. 2004, 20: 235-242. 10.1093/bioinformatics/btg396.View ArticlePubMedGoogle Scholar
- Datta S, Datta S: Empirical Bayes screening of many p-values with applications to microarray studies. Bioinformatics. 2005, 21 (9): 1987-1994. 10.1093/bioinformatics/bti301.View ArticlePubMedGoogle Scholar
- Serban N, Wasserman L: CATS: clustering after transformation and smoothing. J Am Stat Assoc. 2005, 471: 990-999.View ArticleGoogle Scholar
- Ernst J, Nau G, Bar-Joseph Z: Clustering short time series gene expression data. Bioinformatics. 2005, 21: 159-168. 10.1093/bioinformatics/bti1022.View ArticleGoogle Scholar
- Murthy KRK, Hua LJ: Improved Fourier transform method for unsupervised cell-cycle regulated gene prediction. Proc IEEE Comp Sys Bioinformatics Conf Proceedings. 2004, 194-203.Google Scholar
- Kim B, Littell RC, Wu R: Clustering periodic patterns of gene expression based on Fourier approximations. Curr Genomics. 2006, 7: 197-203. 10.2174/138920206777780229.View ArticleGoogle Scholar
- Johansson D, Lindgren P, Berglund A: A multivariate approach applied to microarray data for identification of genes with cell cycle-coupled transcription. Bioinformatics. 2003, 19: 467-473. 10.1093/bioinformatics/btg017.View ArticlePubMedGoogle Scholar
- Schliep A, Schonhuth A, Steinhoff C: Using hidden Markov models to analyze gene expression time course data. Bioinformatics. 2003, 19 (Suppl): i255-i263.View ArticlePubMedGoogle Scholar
- Kim J, Kim H: Clustering of change patterns using Fourier coefficients. Bioinformatics. 2008, 24: 184-191. 10.1093/bioinformatics/btm568.View ArticlePubMedGoogle Scholar
- Murtage C, Raftery AE: Fitting straight lines to point patterns. Pattern Recogn. 1984, 17: 479-483. 10.1016/0031-3203(84)90045-1.View ArticleGoogle Scholar
- Fraley C, Raftery AE: MCLUST: software for model-based cluster analysis. J of Classification. 1999, 16: 297-306. 10.1007/s003579900058.View ArticleGoogle Scholar
- Yuan M: Flexible temporal expression profile modelling using the Gaussian process. Comput Stat Data Anal. 2006, 51 (3): 1754-1764. 10.1016/j.csda.2005.11.017.View ArticleGoogle Scholar
- Kirk PDW, Stumpf MPH: Gaussian process regression bootstrapping: exploring the effects of uncertainty in time course data. Bioinformatics. 2009, 25 (10): 1300-1306. 10.1093/bioinformatics/btp139.PubMed CentralView ArticlePubMedGoogle Scholar
- Lawrence ND, Sanguinetti G, Rattray M: Modelling transcriptional regulation using Gaussian processes. Adv Neural Inf Process Syst. 2007, 19: 785-792.Google Scholar
- Gao P, Honkela A, Rattray M, Lawrence ND: Gaussian process modelling of latent chemical species: applications to inferring transcription factor activities. Bioinformatics. 2008, 24 (16): i70-i75. 10.1093/bioinformatics/btn278.View ArticlePubMedGoogle Scholar
- Kalaitzis AA, Lawrence ND: A simple approach to ranking differentially expressed gene expression time course through Gaussian process regression. BMC Bioinforma. 2011, 12: 180-10.1186/1471-2105-12-180.View ArticleGoogle Scholar
- Angelini C, Cutillo L, Canditiis D, Mutarelli M, Pensky M: BATS: a Bayesian user-friendly software for analyzing time series microarray experiments. BMC Bioinforma. 2008, 9: 415-10.1186/1471-2105-9-415.View ArticleGoogle Scholar
- Spellman PT, Sherlock G, Zhang MQ, Lyer VR, Anders K, Eisen MB, Brown PO, Botstein D, Futcher B: Comprehensive identification of cell cycle-regulated genes of the yeast saccaromyces cerevisiae by microarray hybridization. Mol Biol Cell. 1998, 9: 3273-3297. 10.1091/mbc.9.12.3273.PubMed CentralView ArticlePubMedGoogle Scholar
- Luan Y, Li H: Clustering of time-course gene expression data using mixed-effects models with B-splines. Bioinformatics. 2003, 19: 474-482. 10.1093/bioinformatics/btg014.View ArticlePubMedGoogle Scholar
- Beissbarth T, Speed TP: GOstat: Find Statistically Overrepresented Gene Ontologies Within a Group of Genes. 2004, London: Oxford Univ Press, 1464-1465.Google Scholar
- Brown PO, Botstein D, Futcher B: Comprehensive identification of cell cycle-regulated genes of the yeast saccaromyces cerevisiae by microarray hybridization. Mol Bio Cell. 1998, 9: 3273-3297. 10.1091/mbc.9.12.3273.View ArticleGoogle Scholar
- Lindquist S, Craig EA: The heat-shock proteins. Annu Rev Genet. 1998, 22: 631-677.View ArticleGoogle Scholar
- Smith B, Yaffe M: A mutation in the yeast heat-shock factor gene causes temperature-sensitive defects in both mitochondrial protein import and the cell cycle. Mol Cell Biol. 1991, 11: 2647-2655.PubMed CentralView ArticlePubMedGoogle Scholar
- Wilson AG, Adams RP: Proceedings of the 30th International Conference on Machine Learning. Gaussian Process Kernels for Pattern Discovery and Extrapolation. 2013, Atlanta: International Machine Learning SocietyGoogle Scholar
- Kim J, Hart JD: Test for change when the data are dependent. J Time Series Anal. 1998, 19 (4): 399-424. 10.1111/1467-9892.00100.View ArticleGoogle Scholar
- Eubank R, Hart JD: Testing goodness-of-fit via order selection criteria. Ann Stat. 1992, 20 (3): 1412-1425. 10.1214/aos/1176348775.View ArticleGoogle Scholar
- Berk KN: A central limit theorem for
*m*-dependent random variables with unbounded*m*. Anal Prob. 1973, 1973 (1): 352-354.View ArticleGoogle Scholar - Freedman D, Lane D: The empirical distribution of Fourier coefficients. Anal stat. 1980, 8 (6): 1244-1251. 10.1214/aos/1176345197.View ArticleGoogle Scholar
- Banfield JD, Raftery AE: Model-based Gaussian and non-Gaussian clustering. Biometrics. 1993, 49: 803-821. 10.2307/2532201.View ArticleGoogle Scholar
- Hubert L: Comparing partitions. J. of Classification. 1985, 2: 193-218. 10.1007/BF01908075.View ArticleGoogle Scholar
- Rousseeuw PJ: Silhouettes: graphical aid to the interpretation and validation of cluster analysis. J Comp App Math. 1987, 20: 53-65.View ArticleGoogle 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.