Abstract
Background
Over the past decade, many investigators have used sophisticated time series tools for the analysis of genomic sequences. Specifically, the correlation of the nucleotide chain has been studied by examining the properties of the power spectrum. The main limitation of the power spectrum is that it is restricted to stationary time series. However, it has been observed over the past decade that genomic sequences exhibit nonstationary statistical behavior. Standard statistical tests have been used to verify that the genomic sequences are indeed not stationary. More recent analysis of genomic data has relied on timevarying power spectral methods to capture the statistical characteristics of genomic sequences. Techniques such as the evolutionary spectrum and evolutionary periodogram have been successful in extracting the timevarying correlation structure. The main difficulty in using timevarying spectral methods is that they are extremely unstable. Large deviations in the correlation structure results from very minor perturbations in the genomic data and experimental procedure. A fundamental new approach is needed in order to provide a stable platform for the nonstationary statistical analysis of genomic sequences.
Results
In this paper, we propose to model nonstationary genomic sequences by a timedependent autoregressive moving average (TDARMA) process. The model is based on a classical ARMA process whose coefficients are allowed to vary with time. A series expansion of the timevarying coefficients is used to form a generalized YuleWalkertype system of equations. A recursive leastsquares algorithm is subsequently used to estimate the timedependent coefficients of the model. The nonstationary parameters estimated are used as a basis for statistical inference and biophysical interpretation of genomic data. In particular, we rely on the TDARMA model of genomic sequences to investigate the statistical properties and differentiate between coding and noncoding regions in the nucleotide chain. Specifically, we define a quantitative measure of randomness to assess how far a process deviates from white noise. Our simulation results on various gene sequences show that both the coding and noncoding regions are nonrandom. However, coding sequences are "whiter" than noncoding sequences as attested by a higher index of randomness.
Conclusion
We demonstrate that the proposed TDARMA model can be used to provide a stable time series tool for the analysis of nonstationary genomic sequences. The estimated timevarying coefficients are used to define an index of randomness, in order to assess the statistical correlations in coding and noncoding DNA sequences. It turns out that the statistical differences between coding and noncoding sequences are more subtle than previously thought using stationary analysis tools: Both coding and noncoding sequences exhibit statistical correlations, with the coding regions being "whiter" than the noncoding regions. These results corroborate the evolutionary periodogram analysis of genomic sequences and revoke the stationary analysis' conclusion that coding DNA behaves like random sequences.
Background
Understanding the statistical properties of genomic sequences helps recreate the dynamical processes that led to the current DNA structure, and determine generelated diseases like cancer and Alzheimer disease. For instance, based on the view that noncoding DNA exhibits longrange correlations [16], Li [7] proposed an expansionmodification model of gene evolution. The model incorporates the two basic features of DNA evolution: (i) sequence elongation due to gene duplication and (ii) mutations. It can be shown that the limiting sequence created by this dynamical process exhibits a longrange correlation structure, as attested by a 1/f^{α }spectrum, where the exponent α is a function of the probability of mutation. To understand the relationship between the DNA correlation structure and possible gene abberations, Dodin et al. [8] designed a simple correlation function intended to visualize the regular patterns encountered in DNA sequences. This function is used to revisit the intriguing question of triplet repeats with the aim of providing a visual estimate of the propensity of genes to be highly expressed and/or to lead to possible aberrant structures formed upon strand slippage.
Statistical analysis of genomic sequences has, however, relied, for a long time, on signal processing techniques for stationary signals (correlation and power spectrum) [2,4,9,10], and statistical tools for slowlyvarying trends within stationary signals (Detrended Fluctuation Analysis or DFA) [1,5,6]. Stationarity can be argued as a valid assumption for timeseries of short duration. However, such an assumption rapidly loses its credibility in the enormous databases maintained by various genome projects. Standard statistical tests (e.g., Priestley's test for stationarity) have been used to verify that genomic sequences are not stationary and the nature of their nonstationarity varies and is often more complex than a simple trend [11,12]. Subsequently, more recent analysis of genomic data [1] has relied on timevarying power spectral methods (the evolutionary spectrum and periodogram) to capture the statistical characteristics of genomic sequences [11,12]. The main difficulty in using timevarying spectral methods is that they are extremely unstable and very noisy. Typically, the power spectrum and the evolutionary spectrum are averaged over time in order to obtain smooth and less jittery curves. Moreover, as was pointed out in [13], the evolutionary spectrum is restricted to the class of oscillatory processes. A stochastic process, X(t), is oscillatory if it has a representation of the form
Where Z(λ) is an orthogonal increment process, and the evolutionary power spectrum of the process is defined by P (t, λ) = A(t, λ)^{2}. This definition of the evolutionary power spectrum has the following disadvantages [13]:
(i) It is not uniquely defined for a given nonstationary process.
(ii) The estimation procedure for the evolutionary spectrum depends greatly on the nature of theamplitude function A(t, λ), which is not known a priori.
(iii) An increase in the number of observations does not provide added information on the local behavior of the evolutionary spectrum, and thus does not improve estimation accuracy.
We propose to model nonstationary genomic sequences by a timedependent autoregressive moving average (TDARMA) process. Cramer [14] showed that a nonstationary process still possesses a Wold decomposition in terms of its innovation and its generating system. However, the linear system generating the nonstationary signal, x(t), when driven by the innovation, w(t), is no longer shiftinvariant; the parameters of the impulse response, h_{u}, of this system are timedependent so that
The existence of a timevarying ARMA representation of this process is ensured by two theorems due, independently, to Grenier [15] and Huang and Aggarwal [16]. The uniqueness of the TDARMA representation is obtained by constraining the ARMA model to be invertible, but this leads to conditions on the timevarying impulse response {h_{u}(t)} and its inverse (namely to be absolutely summable at any time t), which cannot be easily expressed in terms of the timedependent coefficients of the ARMA model. In this paper, we estimate the timedependent coefficients of the general TDARMA model using meansquares, leastsquares and recursive leastsquares algorithms. The meansquares estimation leads to generalized YuleWalker type equations [15]. Once the nonstationary parameters are estimated (as time series), we use them to provide a basis for statistical inference by defining an index of randomness, which quantitatively assesses how close the nonstationary signal is to white noise. Our simulation results on various gene sequences show that (i) both the coding and noncoding segments of a gene are not random, and (ii) the coding segments are "closer" to random sequences than noncoding segments. Our results support the view that both coding and noncoding sequences are not random [11,12,9,1720], and revokes the stationary study that maintains that noncoding DNA sustains longrange correlations whereas coding DNA behaves like random sequences [13,5,6,10].
Methods
Numerical representation of genomic sequences
Converting the DNA sequence into a digital signal offers the opportunity to apply powerful signal processing methods for the handling and analysis of genomic information. This is, however, not an easy task as the analysis results might depend on the chosen map. Various numerical mappings have been adopted in the literature. To cite few, Peng et al. [1] construct a onedimensional map of nucleotide sequences onto a walk, u(i), which they termed "DNA walk". The DNA walk is defined by the rule that the walker steps up (u(i) = +1) if a pyrimidine resides at position i, and steps down (u(i) = 1) otherwise. Voss [9] represents a DNA sequence by four binary indicator sequences, which indicate the locations of the four nucleotides A, T, C and G. Berthelsen et al. [21] proposed a twodimensional representation of DNA sequences, characterized by a Hausdorff dimension (also called fractal dimension) that is invariant under (i) complementarity, (ii) reflection symmetry, (iii) compatibility and (iv) substitution symmetry of A͘T and C↔G. The corresponding embedding assignment is given by A = (1; 0), T = (1; 0), C = (0; 1) and G = (0; 1). In this paper, since we are interested in timedependent ARMA modeling of onedimensional nonstationary genomic sequences, we adopt the widely used "DNA walk" map proposed by Peng et al [1]. The DNA walk provides a nice graphical representation for each gene. For instance, Figure 1 shows the structure of the Human gene 276 located in chromosome 1, and its DNA walk is displayed in Fig. 2.
Figure 1. Gene Structure. Gene structure of the Human gene 276 located in chromosome 1: The boxes correspond to the exons (coding regions), and the lines between them represent the introns (noncoding regions). The total length of the gene is N = 8208 bases, including 1536 coding bases and 6672 noncoding bases.
Figure 2. DNA Walk. DNA walk of the Human gene 276.
Timedependent ARMA model
Grenier [22] showed that a discrete nonstationary signal, {x [n]}, can be represented by finiteorder timevarying ARMA processes of the form
where N is the length of the sequence x [n], a_{i }[n] and b_{i }[n] are the timedependent model parameters, p and q are the model orders and w [n] is a white noise process. Observe that the coefficients a_{i }[n] and b_{i }[n] appear with an argument n  i depending on i. This is purely arbitrary since any time origin can be chosen, without restraining the generality of the model. We assume that the timedependent coefficients a_{i }[n] and b_{i }[n] can be expressed as linear combinations of some basis functions ,
The advantage of the basis parametrization is clear from the fact that the identification of the timedependent coefficients a_{i }[n] and b_{i }[n] reduces to the identification of the constant coefficients and , and therefore the linear nonstationary problem reduces to a linear timeinvariant problem. The basis functions do not have to be limited to the standard choices of Legendre, Fourier, or the prolate spheroidal basis but can also take advantage of any prior information, such as the presence of a jump in the coefficients at a known instant [22].
Estimation of the timedependent ARMA coefficients
From Eqs. (4) and (5), the unknown parameters of the TDARMA model are the weights of the linear combinations defining each timevarying parameter. The linearity is the key to the algorithms proposed in this paper. We will derive meansquares, leastsquares and recursive leastsquares solutions to estimate the timedependent coefficients and .
Meansquares estimation
Define the process
and define the vector
where the symbol ^{t }stands for the transpose of a vector or a matrix. It is possible to derive for this process orthogonality conditions similar to the stationary ARMA model conditions [23]. Observe that the process v [n], defined in Eq. (6), is orthogonal to [w[n  q  1], w [n  q  2], ⋯]; hence, it is orthogonal to x [n  q  i] for all i > 0, and orthogonal to X [n  q  i] for all i > 0 [22]. This orthogonality condition leads to a generalized YuleWalker equation [22]
Although the process x [n] is nonstationary, the stationarity and ergodicity of the process w [n], together with the linearity of the model, allow us to replace in Eq. (8) the expectation by a summation. However, although consistent, the above estimator is often considered a poor one [22].
Leastsquares estimation
Equations (4) and (5) can be written in vector format as follows
where
Define
Then, we have
Using this vector notation, Eq. (3) can be written as
Or equivalently
where φ^{t }[n] is the row vector
and
Observe that the vector θ contains all the unknowns of the TDARMA model. Writing Eq. (10) for n = 0, 1, ⋯, N  1 leads to
where
The leastsquares solution of Eq. (11) is given by
To overcome the computational complexity associated with the leastsquares estimation (involving in particular the inversion of a square (m + 1)(p + q) × (m + 1)(p + q) matrix), we opted for a recursive leastsquares estimation as follows.
Recursive leastsquares estimation
The recursive least squares algorithm is summarized as [24]
The initial conditions can be chosen arbitrarily.
Index of randomness
Over the past decade, there has been a flow of conflicting papers about the longrange powerlaw correlations detected in eukaryotic DNA [13,5,6,912,1720]. The controversy is generated by conflicting views that either advocate that noncoding DNA sustains longrange correlations whereas coding DNA behaves like random sequences [1,10,3,5,6], or maintains that both coding and noncoding DNA exhibit longrange powerlaw correlations [11,12,9,1720]. Based on the analysis of the timedependent power spectrum of genomic sequences, Bouaynaya and Schonfeld [11,12] showed that the statistical differences between coding and noncoding sequences are more subtle than previously concluded using stationary analysis tools. In fact they found that both coding and noncoding sequences are nonrandom. However, coding sequences are "whiter" than noncoding sequences.
We propose to qualitatively assess the degree of randomness of both coding and noncoding sequences using the timedependent ARMA coefficients a_{i }[n] and b_{i }[n]. Consider the system function, H (z), of a stationary ARMA model (whose coefficients a_{i }and b_{i }are constant, i.e., independent of time). We know that
where (resp. ) are the zeros (resp. poles) of the system function. From the fact that a stationary white noise process has a at spectrum, we observe that the closer (in L_{2 }distance) the zeros and poles are, the flatter is the spectrum of the process. Following the same reasoning locally for nonstationary processes, we define the curve of randomness, CR [n], of the nonstationary process x [n] by
where the minimum is taken over all pairs (r_{k }[n], p_{k }[n]). Observe that the case p <q is obtained from the p >q case by interchanging the roles of r_{k }and p_{k}, and the indices p and q. The curve of randomness defined in Eq. (17) is a measure of how close the zeros and the poles of the system function are, and therefore, is a measure of how flat the system function is, and how close is the process from a white noise. The index of randomness, IR(p, q), of a TDARMA(p, q), is then defined as the average of the curve of randomness, i.e.,
In particular, the index of randomness of a TDARMA(1,1) (x [n] + a[n  1]x[n  1] = w[n] + b[n]w[n  1]) is given by
Observe that the index of randomness of a white noise process is equal to zero. We say that the sequence x_{1 }[n] with index of randomness IR_{1 }is more random than the sequence x_{2 }[n] with index of randomness IR_{2 }if the index of randomness of the former is lower than the index of randomness of the latter, i.e., IR_{1 }<IR_{2}.
Results
All genome sequences considered in this paper have been extracted from the NIH website http://www.ncbi.nlm.nih.gov webcite. The algorithms were implemented in MATLAB. The DNA sequences were mapped to numerical sequences using the purinepyrimidine DNA walk proposed in [1]. In our simulations, the recursive least squares algorithm was found to best estimate the timedependent coefficients of the TDARMA model. We used the MATLAB function rarmax, which implements the recursive leastsquares algorithm for TDARMA models. The choice of the orders p and q of the model were determined experimentally as follows: For each genomic sequence, we computed 100 TDARMA models corresponding to the orders (1, 1) up to (10, 10). The best model was chosen to be the one that minimizes the average squared error between the actual and the fitted sequences. Our extensive simulations on various DNA sequences from different organisms show that most of the sequences are best fitted with loworder TDARMA models, e.g., TDARMA(1,1), TDARMA(1,2) and TDARMA(2,1). Figure 3 shows the DNA walk of the Human gene 276 and its TDARMA(1,1) fitted sequence. Observe that the TDARMA(1,1) model accurately fits this gene sequence. The estimated timevarying coefficients a [n] and b [n] are displayed in Fig. 4 for both the coding and noncoding regions of this gene. Their statistical differences are not clear from the plot of the timeseries coefficients. The curves of randomness of the coding and noncoding regions are displayed in Fig. 5. Table 1 shows the index of randomness of various gene sequences. For concise representation, the column titles have been abbreviated as follows: "C. length" (resp."N.C. length") denotes the length (in base pairs) of the coding (resp. noncoding) segment of the gene. The total length of the gene is the sum of the lengths of its coding and noncoding regions. "C. (p, q)" (resp. "N.C. (p, q)") denotes the optimal TDARMA parameters (p, q) for the coding (resp. noncoding) region of the gene. Finally, "C. IR" (resp. "N.C. IR") is the index of randomness of the coding (resp. noncoding) segment of the gene. Observe that, in all considered genes, the index of randomness of both the coding and noncoding segments are strictly positive, and the index of randomness of the coding region is consistently lower than the index of randomness of the noncoding region (recall that the index of randomness of a white noise is zero). These observations bring to bear two important conclusion: (i) Both the coding and noncoding sequences are not random, as attested by an index of randomness greater than zero. (ii) The coding sequences are "whiter" than the noncoding sequences. This conclusion revokes previous work on statistical correlation of DNA sequences, which, based on stationary timeseries analysis, presumed that coding DNA behaves like random sequences [13,5,6,10]; and supports the conflicting view that both coding and noncoding sequences are not random [11,12,9,1720]. In particular, our conclusion is in accordance with the evolutionary periodogram analysis conducted in [11,12].
Figure 3. TDARMA modeling. TDARMA modeling of the Human gene 276: The blue signal is the DNA walk, and the red signal is the TDARMA(1,1) fitted signal. The TDARMA(1,1) model accurately fits the genomic signal.
Figure 4. TDARMA coefficients estimation. Estimation of the TDARMA(1,1) coefficients of the Human gene 276. The TDARMA(1,1) model is given by x [n] + a [n  1] x [n  1] = w [n] + b [n  1] w [n  1]. The blue and black (resp. red and green) curves depict the time series a[n] (resp. b[n]) for the coding and noncoding regions of the gene, respectively.
Figure 5. Curve of randomness. The curves of randomness of the coding and noncoding regions of the Human gene 276 are shown in blue and red, respectively. The index of randomness of the coding sequence is equal to 1.0603, whereas its corresponding value for the noncoding sequence is equal to 1.0627.
Table 1. Index of Randomness of the Coding and NonCoding segments of Various Gene Sequences
Conclusion
In this paper, we modelled the nonstationary genomic sequences by a timedependent autoregressive moving average (TDARMA) model. By expressing the timedependent coefficients as linear combinations of parametric basis functions, we were able to transform a linear nonstationary problem into a linear timeinvariant problem. Subsequently, we proposed three methods to estimate the timedependent coefficients: Mean square, leastsquares, and recursive leastsquares algorithms. Based on the estimated TDARMA coefficients, we defined an index of randomness to quantify the degree of randomness of both coding and noncoding sequences. We found that both coding and noncoding sequences are not random. However, a higher index of randomness attests that coding sequences are "whiter" than noncoding sequences. These results corroborate the evolutionary periodogram analysis of genomic sequences performed in [11] and [12], and revoke the stationary analysis' conclusion that coding DNA behaves like random sequences.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
JSZ derived the different estimation algorithms of the TDARMA parameters and performed the simulations. NB proposed the use of the nonstationary analysis and the index of randomness as a basis for statistical inference and biophysical interpretation of genomic data, derived the different estimation algorithms of the TDARMA parameters, and drafted the manuscript. DS proposed the use of the nonstationary analysis and the index of randomness as a basis for statistical inference and biophysical interpretation of genomic data and derived the different estimation algorithms of the TDARMA parameters. WO proposed the use of TDARMA modeling as a nonstationary model of genomic sequences. All authors read and approved the final manuscript.
Acknowledgements
The publication of this paper was partially supported by the Arkansas IDeA Networks of Biomedical Research Excellence (INBRE).
This article has been published as part of BMC Bioinformatics Volume 9 Supplement 9, 2008: Proceedings of the Fifth Annual MCBIOS Conference. Systems Biology: Bridging the Omics. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/9?issue=S9
References

Peng CK, Buldyrev SV, Goldberger AL, Havlin S, Sciortino F, Simons M, Stanley HE: Longrange correlations in nucleotide sequences.
Nature 1992, 356(6365):168170. PubMed Abstract  Publisher Full Text

Buldyrev SV, Goldberger AL, Havlin S, Mantegna RN, Matsa ME, Peng CK, Simons M, Stanley HE: Longrange correlation properties of coding and noncoding DNA sequences: GenBank analysis.

Stanley HE, Buldyrev SV, Goldberger AL, Havlin S, Peng CK, Simons M: Scaling features of noncoding DNA.
Physica A 1999, 273:118. PubMed Abstract

Li W, Holste D: Universal 1/f noise, crossovers of scaling exponents, and chromosomespecific patterns of guaninecytosine content in DNA sequences of the human genome.

Podobnik B, Shao J, Dokholyan NV, Zlatic V, Stanley HE, Grosse I: Similarity and dissimilarity in correlations of genomic DNA.

Carpena P, BernaolaGalvan P, Coronado AV, Hackenberg M, Oliver JL: Identifying chracteristic scales in the human genome.

Li W: Expansionmodification systems: a model for spatial 1/f spectra.
Physical Review A 1991, 43(10):52405260. PubMed Abstract  Publisher Full Text

Dodin G, Levoir P, Cordier C: Triplet Correlation in DNA Sequences and Stability of Heteroduplexes.
Journal of Theoretical Biology 1996, 183:341343. PubMed Abstract  Publisher Full Text

Voss RF: Evolution of longrange fractal correlations and 1/f noise in DNA base sequences.
Physical Review Letters 1992, 68:38053808. PubMed Abstract  Publisher Full Text

Li W, Kaneko K: Longrange correlation and partial 1/f spectrum in a noncoding DNA sequence.

Bouaynaya N, Schonfeld D: NonStationary Analysis of Genomic Sequences. In IEEE Statistical Signal Processing Workshop. Madison, WI; 2007:200204.

Bouaynaya N, Schonfeld D: Nonstationary Analysis of Coding and Noncoding Regions in Nucleotide Sequences.

ADAK S: Timedependent spectral analysis of nonstationary time series.
Journal of the American Statistical Association 1998, 93(444):14881501.

Cramer H: On some classes of nonstationary stochastic processes. In Proceedings of the Berkeley Symppsium on Math, Statistics, and Probability. Los Angeles, CA; 1961.

Grenier Y: Rational nonstationary spectra and their estimation.

Huang NC, Aggarwal JK: On linear Shiftvariant digital filters.
IEEE Transactions on Circuits and Systems 1980, 27(8):672679.

ChatzidimitriouDreismann CA, Larhammar D: Longrange correlations in DNA.
Nature 1993, 361:212. PubMed Abstract  Publisher Full Text

Pande VS, Grosberg AY, Tanaka T: Nonrandomness in protein sequences – evidence for a physically driven stage of evolution.
Proceedings of the National Academy of Sciences 1994, 91(26):1297212975.

Guharay S, Hunt BR, York JA, White OR: Correlations in DNA sequences across the three domains of life.

Berthelsen CL, Glazier JA, Skolnick MH: Global fractal dimension of human DNA sequences treated as pseudorandom walks.
Physical Review A 1992, 45(12):89028913. PubMed Abstract  Publisher Full Text

Grenier Y: TimeDependent ARMA Modeling of Nonstationary Signals.
IEEE Transactions on Acoustics, Speech, and Signal Processing 1983, 31(4):899911.

Statistical digital signal processing and modeling. Wiley. 1996.

Ljung L: System Identification – Theory for the User. second edition. Prentice Hall; 2006.