Home | About | Journals | Submit | Contact Us | Français |

**|**Adv Bioinformatics**|**v.2013; 2013**|**PMC3600260

Formats

Article sections

- Abstract
- 1. Introduction
- 2. RIAA Algorithm
- 3. Methods
- 4. Results
- 5. Conclusions
- Supplementary Material
- References

Authors

Related links

Adv Bioinformatics. 2013; 2013: 171530.

Published online 2013 February 28. doi: 10.1155/2013/171530

PMCID: PMC3600260

Kwadwo S. Agyepong,^{
1
} Fang-Han Hsu,^{
1
} Edward R. Dougherty,^{
1
,}^{
2
} and Erchin Serpedin^{
1
,}^{*}

*Erchin Serpedin: Email: ude.umat.ece@nidepres

Academic Editor: Mohamed Nounou

Received 2012 October 26; Accepted 2013 January 23.

Copyright © 2013 Kwadwo S. Agyepong et al.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Time-course expression profiles and methods for spectrum analysis have been applied for detecting transcriptional periodicities, which are valuable patterns to unravel genes associated with cell cycle and circadian rhythm regulation. However, most of the proposed methods suffer from restrictions and large false positives to a certain extent. Additionally, in some experiments, arbitrarily irregular sampling times as well as the presence of high noise and small sample sizes make accurate detection a challenging task. A novel scheme for detecting periodicities in time-course expression data is proposed, in which a real-valued iterative adaptive approach (RIAA), originally proposed for signal processing, is applied for periodogram estimation. The inferred spectrum is then analyzed using Fisher's hypothesis test. With a proper *p*-value threshold, periodic genes can be detected. A periodic signal, two nonperiodic signals, and four sampling strategies were considered in the simulations, including both bursts and drops. In addition, two yeast real datasets were applied for validation. The simulations and real data analysis reveal that RIAA can perform competitively with the existing algorithms. The advantage of RIAA is manifested when the expression data are highly irregularly sampled, and when the number of cycles covered by the sampling time points is very reduced.

Patterns of periodic gene expression have been found to be associated with essential biological processes such as cell cycle and circadian rhythm [1], and the detection of periodic genes is crucial to advance our understanding of gene function, disease pathways, and, ultimately, therapeutic solutions. Using high-throughput technologies such as microarrays, gene expression profiles at discrete time points can be derived and hundreds of cell cycle regulated genes have been reported in a variety of species. For example, Spellman et al. applied cell synchronization methods and conducted time-course gene expression experiments on *Saccharomyces cerevisiae* [2]. The authors identified 800 cell cycle regulated genes using DNA microarrays. Also, Rustici et al. and Menges et al. identified 407 and about 500 cell cycle regulated genes in *Schizosaccharomyces pombe* and *Arabidopsis*, respectively [3, 4].

Signal processing in the frequency domain simplifies the analysis and an emerging number of studies have demonstrated the power of spectrum analysis in the detection of periodic genes. Considering the common issues of missing values and noise in microarray experiments, Ahdesmäki et al. proposed a robust detection method incorporating the fast Fourier transform (FFT) with a series of data preprocessing and hypothesis testing steps [5]. Two years later, the authors further proposed a modified version for expression data with unevenly spaced time intervals [6]. A Lomb-Scargle (LS) approach, originally used for finding periodicities in astrophysics, was developed for expression data with uneven sampling [7]. Yang et al. further improved the performance using a detrended fluctuation analysis [8]. It used harmonic regression in the time domain for significance evaluation. The method was termed “Lomb-Scargle periodogram and harmonic regression (LSPR).” Basically, these methods consists of two steps: transferring the signals into the frequency (spectral) domain and then applying a significance evaluation test for the resulting peak in the spectral density.

While numerous methods have been developed for detecting periodicities in gene expression, most of these methods suffer from false positive errors and working restrictions to a certain extent, particularly when the time-course data contain limited time points. In addition, no algorithm seems available to resolve all of these challenges. Microarray as well as other high-throughput experiments, due to high manufacturing and preparation costs, have common characteristics of small sample size [9], noisy measurements [10], and arbitrary sampling strategies [11], thereby making the detection of periodicities highly challenging. Since the number and functions of cell cycle regulated genes, or periodic genes, remain greatly uncertain, advances in detection algorithms are urgently needed.

Recently, Stoica et al. developed a novel nonparametric method, termed the “real-valued iterative adaptive approach (RIAA),” specifically for spectral analysis with nonuniformly sampled data [12]. As stated by the authors, RIAA, an iteratively weighted least-squares periodogram, can provide robust spectral estimates and is most suitable for sinusoidal signals. These characteristics of RIAA inspired us to apply it to time-course gene expression data and conduct an examination on its performance. Herein, we incorporate RIAA with a Fisher's statistic to detect transcriptional periodicities. A rigorous comparison of RIAA with several aforementioned algorithms in terms of sensitivities and specificities is conducted through simulations and simulation results dealing with real data analysis are also provided.

In this study, we found that the RIAA algorithm can provide robust spectral estimates for the detection of periodic genes regardless of the sampling strategies adopted in the experiments or the nonperiodic nature of noise present in the measurement process. We show through simulations that the RIAA can outperform the existing algorithms particularly when the data are highly irregularly sampled, and when the number of cycles covered by the sampling time points is very few. These characteristics of RIAA fit perfectly the needs of time-course gene expression data analysis. This paper is organized as follows. In Section 2, we begin with an overview of RIAA. In Section 3, a scheme for detecting periodicities is proposed, and simulation models for performance evaluation and a real data analysis for validation purposes are presented. A complete investigation of the performance of RIAA and a rigorous comparison with other algorithms are provided in Section 4.

RIAA is an iterative algorithm developed for finding the least-squares periodogram with the utilization of a weighted function. The essential mathematics involved in RIAA is introduced in this section with the algorithm input being time-course expression data; for more details regarding RIAA, the readers are encouraged to check the original paper by Stoica et al. [12].

Suppose that the signals associated with the periodic gene expressions are composed of noise and sinusoidal components. Let *y*
_{h}(*t*
_{i}), *i* = 1,…, *n*, denote the time-course expression ratios of gene *h* at instances *t*
_{1},…, *t*
_{n}, respectively; *y*
_{h}(*t*
_{i}) are real numbers; ∑_{i=1}
^{n}
*y*
_{h}(*t*
_{i}) = 0. The least-squares periodogram Φ_{lsp} is given by

$$\begin{array}{c}{\mathrm{\Phi}}_{lsp}={\left|\widehat{\alpha}\left(\omega \right)\right|}^{\mathrm{2}},\end{array}$$

(1)

where $\widehat{\alpha}(\omega )$ is the solution to the following fitting problem:

$$\begin{array}{c}\widehat{\alpha}\left(\omega \right)=\mathrm{arg}\underset{\alpha \left(\omega \right)}{\mathrm{min}}\sum _{i=\mathrm{1}}^{n}{\left[{y}_{h}\left({t}_{i}\right)-\alpha \left(\omega \right){e}^{j\omega {t}_{i}}\right]}^{\mathrm{2}}.\end{array}$$

(2)

Let *α*(*ω*) = |*α*(*ω*) | *e*
^{jϕ(ω)} = *βe*
^{jθ}, where *β* = |*α*(*ω*)|≥0 and *θ* = *ϕ*(*ω*) [0,2*π*] refer to the amplitude and phase of *α*(*ω*), respectively. The criterion in (2) can then be rewritten as

$$\begin{array}{c}\sum _{i=\mathrm{1}}^{n}{\left[{y}_{h}\left({t}_{i}\right)-\beta \mathrm{cos}\left(\omega {t}_{i}+\theta \right)\right]}^{\mathrm{2}}+{\beta}^{\mathrm{2}}\sum _{i=\mathrm{1}}^{n}{\mathrm{sin}}^{\mathrm{2}}\left(\omega {t}_{i}+\theta \right).\end{array}$$

(3)

The second term in the above equation is data independent and can be omitted from the minimization operation. Hence, the criterion (2) is simplified to

$$\begin{array}{c}\left(\widehat{\beta},\widehat{\theta}\right)=\mathrm{arg}\underset{\beta ,\theta}{\mathrm{min}}\sum _{i=\mathrm{1}}^{n}{\left[{y}_{h}\left({t}_{i}\right)-\beta \mathrm{cos}\left(\omega {t}_{i}+\theta \right)\right]}^{\mathrm{2}}.\end{array}$$

(4)

We further apply *a* = *β*cos (*θ*) and *b* = −*β*sin(*θ*) and derive an equivalent of (4) as follows:

$$\begin{array}{c}\left(\widehat{a},\widehat{b}\right)=\mathrm{arg}\underset{a,b}{\mathrm{min}}\sum _{i=\mathrm{1}}^{n}{\left[{y}_{h}\left({t}_{i}\right)-a\mathrm{cos}\left(\omega {t}_{i}\right)-b\mathrm{sin}\left(\omega {t}_{i}\right)\right]}^{\mathrm{2}}.\end{array}$$

(5)

The target of interest to the fitting problem now becomes $\widehat{a}$ and $\widehat{b}$ (instead of *α*(*ω*)), and the solution is well known to be

$$\begin{array}{c}\left[\begin{array}{c}\widehat{a}\\ \widehat{b}\end{array}\right]={\mathbf{R}}^{-\mathbf{1}}\mathbf{r},\end{array}$$

(6)

where

$$\begin{array}{c}\mathbf{R}=\sum _{i=\mathrm{1}}^{n}\left[\begin{array}{cc}\mathrm{cos}{\left(\omega {t}_{i}\right)}^{\mathrm{2}}& \mathrm{cos}\left(\omega {t}_{i}\right)\mathrm{sin}\left(\omega {t}_{i}\right)\\ \mathrm{sin}\left(\omega {t}_{i}\right)\mathrm{cos}\left(\omega {t}_{i}\right)& \mathrm{sin}{\left(\omega {t}_{i}\right)}^{\mathrm{2}}\end{array}\right],\\ \mathbf{r}=\sum _{i=\mathrm{1}}^{n}\left[\begin{array}{c}\mathrm{cos}\left(\omega {t}_{i}\right)\\ \mathrm{sin}\left(\omega {t}_{i}\right)\end{array}\right]{y}_{h}\left({t}_{i}\right).\end{array}$$

(7)

After $\widehat{a}$ and $\widehat{b}$ are estimated, the least-squares periodogram can be derived.

Prior to implementation of RIAA for periodogram estimation, the observation interval [0, *ω*
_{max }] and the resolution in terms of grid size have to be selected. To this end, the maximum frequency *ω*
_{max } in the observation interval without aliasing errors for sampling instances *t*
_{1},…, *t*
_{n}, can be evaluated by

$$\begin{array}{c}{\omega}_{\mathrm{max}}=\frac{{\omega}_{\mathrm{0}}}{\mathrm{2}},\end{array}$$

(8)

where *ω*
_{0} is given by

$$\begin{array}{c}{\omega}_{\mathrm{0}}=\frac{\mathrm{2}\left(n-\mathrm{1}\right)\pi}{{\sum}_{i=\mathrm{1}}^{n-\mathrm{1}}\left({t}_{i+\mathrm{1}}-{t}_{i}\right)}.\end{array}$$

(9)

The observation interval [0, *ω*
_{max }] is hence chosen after *ω*
_{max } is obtained.

To ensure that the smallest frequency separation in time-course expression data with regular or irregular sampling can be adequately detected, the grid size Δ*ω* is chosen to be

$$\begin{array}{c}\mathrm{\Delta}\omega =\frac{\mathrm{2}\pi}{{t}_{n}-{t}_{\mathrm{1}}},\end{array}$$

(10)

which, in fact, is the resolution limit of the least-squares periodogram. As a result, the frequency grids *ω*
_{g} considered in periodogram are

$$\begin{array}{c}{\omega}_{g}=g\mathrm{\Delta}\omega ,\hspace{1em}g=\mathrm{1},\dots ,G,\end{array}$$

(11)

where the number of grids *G* is given by

$$\begin{array}{c}G=\lfloor \frac{{\omega}_{\mathrm{max}}}{\mathrm{\Delta}\omega}\rfloor .\end{array}$$

(12)

The following notations are introduced for the implementation of RIAA at a specific frequency *ω*
_{g}:

$$\begin{array}{c}\mathbf{Y}={\left[\begin{array}{c}{y}_{h}\left({t}_{\mathrm{1}}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\cdots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{y}_{h}\left({t}_{n}\right)\end{array}\right]}^{T},\\ {\rho}_{g}={\left[\begin{array}{c}a\left({\omega}_{g}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}b\left({\omega}_{g}\right)\end{array}\right]}^{T},\\ {\mathbf{A}}_{g}=\left[\begin{array}{c}{\mathbf{c}}_{g}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathbf{s}}_{g}\end{array}\right],\end{array}$$

(13)

where

$$\begin{array}{c}{\mathbf{c}}_{g}={\left[\begin{array}{c}\mathrm{cos}\left({\omega}_{g}{t}_{\mathrm{1}}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\cdots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{cos}\left({\omega}_{g}{t}_{n}\right)\end{array}\right]}^{T},\\ {\mathbf{s}}_{g}={\left[\begin{array}{c}\mathrm{sin}\left({\omega}_{g}{t}_{\mathrm{1}}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\cdots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{sin}\left({\omega}_{g}{t}_{n}\right)\end{array}\right]}^{T},\end{array}$$

(14)

and *a*(*ω*
_{g}) and *b*(*ω*
_{g}) denote variables *a* and *b* at frequency *ω*
_{g}, respectively.

RIAA's salient feature is the addition of a weighted matrix **Q**
_{g} to the least-squares fitting criterion. The weighted matrix **Q**
_{g} can be viewed as a covariance matrix encapsulating the contributions of noise and other sinusoidal components in **Y** other than *ω*
_{g} to the spectrum; it is defined as

$$\begin{array}{c}{\mathbf{Q}}_{g}=\mathrm{\Sigma}+\sum _{m=\mathrm{1},m\ne g}^{G}{\mathbf{A}}_{m}{\mathbf{D}}_{m}{\mathbf{A}}_{m}^{T},\end{array}$$

(15)

where

$$\begin{array}{c}{\mathbf{D}}_{m}=\frac{{a}^{\mathrm{2}}\left({\omega}_{g}\right)+{b}^{\mathrm{2}}\left({\omega}_{g}\right)}{\mathrm{2}}\left[\begin{array}{c}\mathrm{1}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{0}\\ \mathrm{0}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{1}\end{array}\right],\end{array}$$

(16)

and Σ denotes the covariance matrix of noise in expression data **Y**, given by

$$\begin{array}{c}\mathrm{\Sigma}=\left[\begin{array}{c}{\sigma}^{\mathrm{2}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\dots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{0}\\ \vdots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\ddots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\vdots \\ \mathrm{0}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\dots \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\sigma}^{\mathrm{2}}\end{array}\right].\end{array}$$

(17)

Assuming that **Q**
_{g} is invertible, in RIAA, a weighted least-squares fitting problem is formulated and considered for finding $\widehat{a}$ and $\widehat{b}$ (instead of using (5)), and it is written in the form of matrices using (13) as follows:

$$\begin{array}{c}{\widehat{\rho}}_{g}=\mathrm{arg}\underset{{\rho}_{g}}{\mathrm{min}}{\left[\mathbf{Y}-{\mathbf{A}}_{g}{\rho}_{g}\right]}^{T}{\mathbf{Q}}_{g}^{-\mathrm{1}}\left[\mathbf{Y}-{\mathbf{A}}_{g}{\rho}_{g}\right].\end{array}$$

(18)

In Stoica et al. [12], the solution to (18) has been shown to be

$$\begin{array}{c}{\widehat{\rho}}_{g}=\frac{{\mathbf{A}}_{g}^{T}{\mathbf{Q}}_{g}^{-\mathrm{1}}\mathbf{Y}}{{\mathbf{A}}_{g}^{T}{\mathbf{Q}}_{g}^{-\mathrm{1}}{\mathbf{A}}_{g}},\end{array}$$

(19)

and the RIAA periodogram at *ω* = *ω*
_{g} can be derived by

$$\begin{array}{c}{\mathrm{\Phi}}_{\text{riaa}}\left({\omega}_{g}\right)=\frac{\mathrm{1}}{n}{\widehat{\rho}}_{g}^{T}\left({\mathbf{A}}_{g}^{T}{\mathbf{A}}_{g}\right){\widehat{\rho}}_{g}.\end{array}$$

(20)

From (15) and (19), it is obvious that **Q**
_{g} and ${\widehat{\rho}}_{g}$ are dependent on each other. An iterative approach (i.e., RIAA) is hence a feasible solution to get the estimate ${\widehat{\rho}}_{g}$ and the weighted matrix **Q**
_{g}.

The iteration for estimating spectrum starts with initial estimates ${\widehat{\rho}}_{g}^{\mathrm{0}}$, in which the elements $\widehat{a}$ and $\widehat{b}$ are given by (6) with *ω* = *ω*
_{g}, *g* = 1,…, *G*. After initialization, the first iteration begins. First, the elements $\widehat{a}$ and $\widehat{b}$ of ${\widehat{\rho}}_{g}^{\mathrm{0}}$ are applied to obtain ${\widehat{\mathbf{D}}}_{m}^{\mathrm{1}}$ using (16). Secondly, to get a good estimate of ${\widehat{\sigma}}^{\mathrm{1}}$, the frequency *ω*
_{p} at which the largest value-*p* is located in the temporary periodogram Φ^{0}(*ω*
_{g}), *g* = 1,…, *G*, derived using (20) with ${\widehat{\rho}}_{g}={\widehat{\rho}}_{g}^{\mathrm{0}}$, is applied for obtaining a reversed engineered signal ${\widehat{\mathbf{Y}}}^{\mathrm{0}}$. The elements ${\widehat{y}}_{h}({t}_{i}),\hspace{0.17em}i=\mathrm{1},\dots ,n$, in ${\widehat{\mathbf{Y}}}^{\mathrm{0}}$ are given by

$$\begin{array}{c}{\widehat{y}}_{h}\left({t}_{i}\right)=\sqrt{\mathrm{2}P}\mathrm{cos}\left({\omega}_{p}{t}_{i}+s\right),\hspace{1em}i=\mathrm{1},\dots ,n.\end{array}$$

(21)

The phase of the cosine function *s* is unknown; however, ${\widehat{\sigma}}^{\mathrm{1}}$ is estimable using

$$\begin{array}{c}{\widehat{\sigma}}^{\mathrm{1}}=\underset{s\in \left[\mathrm{0,2}\pi \right]}{\mathrm{min}}\frac{{||\mathbf{Y}-{\widehat{\mathbf{Y}}}^{\mathrm{0}}||}^{\mathrm{2}}}{n},\end{array}$$

(22)

where ||·|| is the Euclidean norm. With estimates ${\widehat{\mathbf{D}}}_{m}^{\mathrm{1}}$ and ${\widehat{\sigma}}^{\mathrm{1}}$, the estimates ${\widehat{\mathbf{Q}}}_{g}^{\mathrm{1}}$, *g* = 1,…, *G*, in the first iteration are hence given by (15). After this, ${\widehat{\mathbf{Q}}}_{g}^{\mathrm{1}}$ are inserted into the right-hand side of (19) and updated estimates ${\widehat{\rho}}_{g}^{\mathrm{1}}$, *g* = 1,…, *G*, are derived. The algorithm consists of repeating these steps and updating ${\widehat{\mathbf{Q}}}_{g}^{k}$ and ${\widehat{\rho}}_{g}^{k}$ iteratively, where *k* denotes the number of iterations, until a termination criterion is reached. If the process stops at the *K*th iteration, then the final RIAA periodogram is given by (20) using ${\widehat{\rho}}_{g}^{K}$. The pseudocode in Algorithm 1 represents a concise description of the iterative RIAA process.

Figure 1 demonstrates our scheme for periodicity detection and algorithm comparison. The first step involves a periodogram estimation, which converts the time-course gene expression ratios into the frequency domain. Three methods are considered for comparison: RIAA, LS, and a detrend LS (termed DLS), which uses an additional detrend function (developed in LSPR) before regular LS periodogram estimation is applied. The derived spectra are then analyzed using hypothesis testing. This study is conducted using a Fisher's test, with the null hypothesis that there are no periodic signals in the time domain and hence no significantly large peak in the derived spectra. The algorithm performance is evaluated and compared via simulations and receiver operating characteristic (ROC) curves. In real microarray data analysis, three published benchmark sets are utilized as standards of cell cycle genes for performance comparison.

After the spectrum of time-course expression data is obtained via periodogram estimation, a Fisher's statistic *f* for gene *h* with the null hypothesis *H*
_{0} that the peak of the spectral density is insignificant against the alternative hypothesis *H*
_{1} that the peak of the spectral density is significant is applied as

$$\begin{array}{c}{f}_{h}=\frac{{\mathrm{max}}_{\mathrm{1}\le g\le G}\left(\mathrm{\Phi}\left({\omega}_{g}\right)\right)}{{G}^{-\mathrm{1}}{\sum}_{g=\mathrm{1}}^{G}\mathrm{\Phi}\left({\omega}_{g}\right)},\end{array}$$

(23)

where Φ refers to the periodogram derived using RIAA, LS, or DLS. The null hypothesis *H*
_{0} is rejected, and the gene *h* is claimed as a periodic gene if its *p*-value, denoted as *p*
_{h}, is less than or equal to a specific significance threshold. For simplicity, *p*
_{h} is approximated from the asymptotic null distribution of *f* assuming Gaussian noise [13] as follows:

$$\begin{array}{c}{p}_{h}=\mathrm{1}-{e}^{-n{e}^{-{f}_{h}}}.\end{array}$$

(24)

In real data analysis, deviation might be invoked for the estimation of *p*
_{h} when the time-course data is short. This issue was carefully addressed by Liew et al. [14], and, as suggested, alternative methods such as random permutation may provide less deviation and better performance. However, permutation also has limitations such as tending to be conservative [15]. While finding the most robust method for the *p*-value evaluation remains an open question, it gets beyond the scope of this study since the algorithm comparison via ROC curves is threshold independent [16], and the results are unaffected by the deviation.

Simulations are applied to evaluate the performance of RIAA. The simulation models and sampling strategies used for simulations are described in the following paragraphs.

Three models, one for periodic signals and two for nonperiodic signals, are considered as transcriptional signals. Since periodic genes are transcribed in an oscillatory manner, the expression levels *y*
_{s} embedded with periodicities are assumed to be

$$\begin{array}{c}{y}_{s}\left({t}_{i}\right)=M\mathrm{cos}\left({\omega}_{s}{t}_{i}\right)+{\u03f5}_{{t}_{i}},\hspace{1em}i=\mathrm{1},\dots ,n,\end{array}$$

(25)

where *M* denotes the sinusoidal amplitude; *ω*
_{s} refers to the signal frequency; *ϵ*
_{ti} are Gaussian noise independent and identically distributed (i.i.d.) with parameters *μ* and *σ*. For nonperiodic signals, the first model *y*
_{n} is simply composed of Gaussian noise, given by

$$\begin{array}{c}{y}_{n}\left({t}_{i}\right)={\u03f5}_{{t}_{i}},\hspace{1em}i=\mathrm{1},\dots ,n.\end{array}$$

(26)

Additionally, as visualized by Chubb et al., gene transcription can be nonperiodically activated with irregular intervals in a living eukaryotic cell, like pulses turning on and off rapidly and discontinuously [17]. Based on this, the second nonperiodic model *y*
_{n}′ incorporates one additional transcriptional burst and one additional sudden drop into the Gaussian noise, which can be written as

$$\begin{array}{c}{y}_{n}^{\prime}\left({t}_{i}\right)={I}_{b}\left({t}_{i}\right)-{I}_{d}\left({t}_{i}\right)+{\u03f5}_{{t}_{i}},\hspace{1em}i=\mathrm{1},\dots ,n,\end{array}$$

(27)

where *I*
_{b} and *I*
_{d} are indicator functions, equal to 1 at the location of the burst and the drop, respectively, and 0 otherwise. The transcriptional burst assumes a positive pulse while the transcriptional drop assumes a negative pulse. Both of them may be located randomly among all time points and are assumed to last for two time points. In other words, the indicator functions are equal to 1 at two consecutive time points, say, *I*
_{b} = 1 at *t*
_{i} and *t*
_{i+1}. The burst and the drop have no overlap.

As for the choices of sampling time points *t*
_{i}, *i* = 1,…, *n*, four different sampling strategies, one with regular sampling and three with irregular sampling, are considered. First, regular sampling is applied in which all time intervals are set to be 1/*c*, where *c* is a constant. Secondly, a bio-like sampling strategy is invoked. This strategy tends to have more time points at the beginning of time-course experiments and less time points after we set the first 2/3 time intervals as 1/*c* and set the next 1/3 time intervals as 2/*c*. Third, time intervals are randomly chosen between 1/*c* and 2/*c*. The last sampling strategy, in which all time intervals are exponentially distributed with parameter *c*, is less realistic than the others but it is helpful for us to evaluate the performance of RIAA under pathological conditions.

ROC curves are applied for performance comparison. To this end, 10,000 periodic signals were generated using (25) and 10,000 nonperiodic signals were generated using either (26) or (27). Sensitivity measures the proportion of successful detection among the 10,000 periodic signals and specificity measures the proportion of correct claims on the 10,000 nonperiodic simulation datasets. Sampling time points are decided by one of the four sampling strategies and the number of time points *n* is chosen arbitrarily. For all ROC curves in Section 4, *c* = 2 and *n* = 24.

Two yeast cell cycle experiments synchronized using an alpha-factor, one conducted by Spellman et al. [2] and one conducted by Pramila et al. [18], are considered for a real data analysis. The first time-course microarray data, termed dataset alpha and downloaded from the Yeast Cell Cycle Analysis Project website (http://genome-www.stanford.edu/cellcycle/), harbors 6,178 gene expression levels and 18 sampling time points with a 7-minute interval. The second time-course data, termed dataset alpha 38, is downloaded from the online portal for Fred Hutchinson Cancer Research Center's scientific laboratories (http://labs.fhcrc.org/breeden/cellcycle/). This dataset contains 4,774 gene expression levels and 25 sampling time points with a 5-minute interval. Three benchmark sets of genes that have been utilized in Lichtenberg et al. [19] and Liew et al. [20] as standards of cell cycle genes are also applied herein for performance comparison. These benchmark sets, involving 113, 352, and 518 genes, respectively, include candidates of cycle cell regulated genes in yeast proposed by Spellman et al. [2], Johansson et al. [21], Simon et al. [22], Lee et al. [23], and Mewes et al. [24] and are accessible in a laboratory website (http://www.cbs.dtu.dk/cellcycle/).

RIAA performed well in the conducted simulations. As shown in Figure 2(a), a periodic signal (solid line) with amplitude *M* = 1 and frequency *ω*
_{s} = 0.4*π*is sampled using the bio-like sampling strategy, which applies 16 time points in (0,8] and 8 more time points in (8,16]. Gaussian noise with parameters *μ* = 0 and *σ* = 0.5 is assumed during microarray experiments. The resulting time-course expression levels (dots), at a total of 24 time points and the sampling time information were treated as inputs to the RIAA algorithm. Figure 2(b) demonstrates the result of periodogram estimation. In this example, the grid size Δ*ω* was chosen to be 0.065 and a total of 11 amplitudes corresponding to different frequencies were obtained and shown in the spectrum. Using Fisher's test, the peak at the third grid (frequency = 0.195) was found to be significantly large (*p*-value = 2.4 × 10^{−3}), and hence a periodic gene was claimed.

(a) A time-course periodic signal with frequency = 0.2 sampled by the bio-like sampling strategy; 16 time points are assigned to the interval (0,8], and 8 time points are assigned to the interval (8,16]. (b) The periodogram derived using RIAA. The maximum **...**

ROC curves strongly illustrate the performance of RIAA. In Figures Figures33 and and4,4, subplots (a)-(b), (c)-(d), (e)-(f), and (g)-(h) refer to the simulations with regular, bio-like, binomially random, and exponentially random sampling strategies, respectively. Additionally, in the left-hand side subplots (a), (c), (e), and (g), nonperiodic signals were simply Gaussian noise with parameters *μ* = 0 and *σ* = 0.5, while in the right-hand side subplots (b), (d), (f), and (h), nonperiodic signals involve not only the Gaussian noise but also a transcriptional burst and a sudden drop (27). Periodic signals were generated using (25) with amplitude *M* = 1, *c* = 2, and *n* = 24. The only difference in simulation settings between Figures Figures33 and and44 is the frequency of periodic signals; they are *ω*
_{s} = 0.4*π* and 0.1*π*, respectively. As shown in these figures, LS and DLS can perform well as RIAA when the time-course data are regularly sampled, or mildly irregularly sampled; however, when data are highly irregularly sampled, RIAA outperforms the others. The superiority of RIAA over DLS is particularly clear when the signal frequency is small.

Figure 5 illustrates the results of the real data analysis when these three algorithms, namely, the RIAA, LS, and DLS, were applied. On the *x*-axis, the numbers indicate the thresholds *η* that we preserved and classified as periodicities among all yeast genes; on the y-axis, the numbers refer to the intersection of *η* preserved genes and the proposed periodic candidates listed in the benchmark sets. Figures 5(a)–5(c) demonstrate the results derived from dataset alpha when the 113-gene benchmark set, 352-gene benchmark set, and 518-gene benchmark set were applied, respectively. Similarly, Figures 5(d)–5(f) demonstrate the results derived from dataset alpha 38. The RIAA does not result in significant differences in the numbers of intersections when compared to those corresponding to LS and DLS in most of these cases. However, RIAA shows slightly better coverage when the dataset alpha 38 and the 113-gene benchmark set was utilized (Figure 5(d)).

In this study, the rigorous simulations specifically designed to comfort with real experiments reveal that the RIAA can outperform the classical LS and modified DLS algorithms when the sampling time points are highly irregular, and when the number of cycles covered by sampling times is very limited. These characteristics, as also claimed in the original study by Stoica et al. [12], suggest that the RIAA can be generally applied to detect periodicities in time-course gene expression data with good potential to yield better results. A supplementary simulation further shows the superiority of RIAA over LS and DLS when multiple periodic signals are considered (see Supplementary Figures1 available online at http://dx.doi.org/10.1155/2013/171530). From the simulations, we also learned that the addition of a transcriptional burst and a sudden drop to nonperiodic signals (the negatives) does not affect the power of RIAA in terms of periodicity detection. Moreover, the detrend function in DLS, designed to improve LS by removing the linearity in time-course data, may fail to provide improved accuracy and makes the algorithm unable to detect periodicities when transcription oscillates with a very low frequency.

The intersection of detected candidates and proposed periodic genes in the real data analysis (Figure 5) does not reveal much differences among RIAA, LS, and DLS. One possible reason is that the sampling time points conducted in the yeast experiment are not highly irregular (not many missing values are included), since, as demonstrated in Figures 3(a)–3(d), the RIAA just performs equally well as the LS and DLS algorithms when the time-course data are regularly or mildly irregularly sampled. Also, the very limited time points contained in the dataset may deviate the estimation of *p*-values [14] and thus hinder the RIAA from exhibiting its excellence. Besides, the number of true cell cycle genes included in the benchmark sets remains uncertain. We expect that the superiority of RIAA in real data analysis would be clearer in the future when more studies and more datasets become available.

Besides the comparison of these algorithms, it is interesting to note that the bio-like sampling strategy could lead to better detection of periodicities than the regular sampling strategy (as shown in Figures 3(c) and 3(d)). It might be beneficial to apply loose sampling time intervals at posterior periods to prolong the experimental time coverage when the number of time points is limited.

Supplementary Figure s1: demonstrates the simulation results considering multiple periodic signals. Two sinusoidal waves with different frequency settings are superimposed as the periodic signal, and RIAA, LS, and DLS are applied for comparison.

Click here for additional data file.^{(83K, pdf)}

The authors would like to thank the members in the Genomic Signal Processing Laboratory, Texas A&M University, for the helpful discussions and valuable feedback. This work was supported by the National Science Foundation under Grant no. 0915444. The RIAA MATLAB code is available at http://gsp.tamu.edu/Publications/supplementary/agyepong 12a/.

1. Zhao W, Agyepong K, Serpedin E, Dougherty ER. Detecting periodic genes from irregularly sampled gene expressions: a comparison study. *EURASIP Journal on Bioinformatics and Systems Biology*. 2008;2008769293 [PMC free article] [PubMed]

2. Spellman PT, Sherlock G, Zhang MQ, et al. Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. *Molecular Biology of the Cell*. 1998;9(12):3273–3297. [PMC free article] [PubMed]

3. Rustici G, Mata J, Kivinen K, et al. Periodic gene expression program of the fission yeast cell cycle. *Nature Genetics*. 2004;36(8):809–817. [PubMed]

4. Menges M, Hennig L, Gruissem W, Murray JAH. Cell cycle-regulated gene expression in Arabidopsis. *Journal of Biological Chemistry*. 2002;277(44):41987–42002. [PubMed]

5. Ahdesmäki M, Lähdesmäki H, Pearson R, Huttunen H, Yli-Harja O. Robust detection of periodic time series measured from biological systems. *BMC Bioinformatics*. 2005;6(article 117) [PMC free article] [PubMed]

6. Ahdesmäki M, Lähdesmäki H, Gracey A, et al. Robust regression for periodicity detection in non-uniformly sampled time-course gene expression data. *BMC Bioinformatics*. 2007;8(article 233) [PMC free article] [PubMed]

7. Glynn EF, Chen J, Mushegian AR. Detecting periodic patterns in unevenly spaced gene expression time series using Lomb-Scargle periodograms. *Bioinformatics*. 2006;22(3):310–316. [PubMed]

8. Yang R, Zhang C, Su Z. LSPR: an integrated periodicity detection algorithm for unevenly sampled temporal microarray data. *Bioinformatics*. 2011;27(7):1023–1025. [PubMed]

9. Dougherty ER. Small sample issues for microarray-based classification. *Comparative and Functional Genomics*. 2001;2(1):28–34. [PMC free article] [PubMed]

10. Tu Y, Stolovitzky G, Klein U. Quantitative noise analysis for gene expression microarray experiments. *Proceedings of the National Academy of Sciences of the United States of America*. 2002;99(22):14031–14036. [PubMed]

11. Bar-Joseph Z. Analyzing time series gene expression data. *Bioinformatics*. 2004;20(16):2493–2503. [PubMed]

12. Stoica P, Li J, He H. Spectral analysis of nonuniformly sampled data: a new approach versus the periodogram. *IEEE Transactions on Signal Processing*. 2009;57(3):843–858.

13. Fan J, Yao Q. *Nonlinear Time Series: Nonparametric and Parametric Methods*. New York, NY, USA: Springer; 2003.

14. Liew AWC, Law NF, Cao XQ, Yan H. Statistical power of Fisher test for the detection of short periodic gene expression profiles. *Pattern Recognition*. 2009;42(4):549–556.

15. Berger V. Pros and cons of permutation tests in clinical trials. *Statistics in Medicine*. 2000;19(10):1319–1328. [PubMed]

16. Bradley AP. The use of the area under the ROC curve in the evaluation of machine learning algorithms. *Pattern Recognition*. 1997;30(7):1145–1159.

17. Chubb JR, Trcek T, Shenoy SM, Singer RH. Transcriptional pulsing of a developmental gene. *Current Biology*. 2006;16(10):1018–1025. [PMC free article] [PubMed]

18. Pramila T, Wu W, Noble W, Breeden L. Periodic genes of the yeast Saccharomyces cerevisiae: a combined analysis of five cell cycle data sets. 2007.

19. Lichtenberg U, Jensen LJ, Fausbøll A, Jensen TS, Bork P, Brunak S. Comparison of computational methods for the identification of cell cycle-regulated genes. *Bioinformatics*. 2005;21(7):1164–1171. [PubMed]

20. Liew AWC, Xian J, Wu S, Smith D, Yan H. Spectral estimation in unevenly sampled space of periodically expressed microarray time series data. *BMC Bioinformatics*. 2007;8(article 137) [PMC free article] [PubMed]

21. 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(4):467–473. [PubMed]

22. Simon I, Barnett J, Hannett N, et al. Serial regulation of transcriptional regulators in the yeast cell cycle. *Cell*. 2001;106(6):697–708. [PubMed]

23. Lee TI, Rinaldi NJ, Robert F, et al. Transcriptional regulatory networks in *Saccharomyces cerevisiae*
. *Science*. 2002;298(5594):799–804. [PubMed]

24. Mewes HW, Frishman D, Güldener U, et al. MIPS: a database for genomes and protein sequences. *Nucleic Acids Research*. 2002;30(1):31–34. [PMC free article] [PubMed]

Articles from Advances in Bioinformatics are provided here courtesy of **Hindawi**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |