Map-invariant spectral analysis for the identification of DNA periodicities

Many signal processing based methods for finding hidden periodicities in DNA sequences have primarily focused on assigning numerical values to the symbolic DNA sequence and then applying spectral analysis tools such as the short-time discrete Fourier transform (ST-DFT) to locate these repeats. The key results pertaining to this approach are however obtained using a very specific symbolic to numerical map, namely the so-called Voss representation. An important research problem is to therefore quantify the sensitivity of these results to the choice of the symbolic to numerical map. In this article, a novel algebraic approach to the periodicity detection problem is presented and provides a natural framework for studying the role of the symbolic to numerical map in finding these repeats. More specifically, we derive a new matrix-based expression of the DNA spectrum that comprises most of the widely used mappings in the literature as special cases, shows that the DNA spectrum is in fact invariable under all these mappings, and generates a necessary and sufficient condition for the invariance of the DNA spectrum to the symbolic to numerical map. Furthermore, the new algebraic framework decomposes the periodicity detection problem into several fundamental building blocks that are totally independent of each other. Sophisticated digital filters and/or alternate fast data transforms such as the discrete cosine and sine transforms can therefore be always incorporated in the periodicity detection scheme regardless of the choice of the symbolic to numerical map. Although the newly proposed framework is matrix based, identification of these periodicities can be achieved at a low computational cost.

Many researchers have noted that the occurrence of repetitive structures in a DNA sequence is symptomatic of a biological phenomena. Specific applications of this observation include identification of diseases [1], DNA forensics [2], and detection of pathogen exposure [3]. Some of these structures are simple repetition of short DNA segments such as exons [4], tandem repeats [5], dispersed repeats [6], and unstable triplet repeats in the noncoding regions [7] while other forms more elaborate patterns such as palindromes [8] and the period-3 component [9-13], a strong periodic characteristic found primarily in genes and pseudogenes [14]. Methods that detect these DNA periodicities are either probabilistic or deterministic. Most of the deterministic techniques rely on spectral analysis of the DNA sequence using the short-time discrete Fourier transform (ST-DFT) [15-17]. The main idea is as follows: given a DNA sequence of length N, numerical values are first assigned to every element in ={ACGT}, where these letters denote the four nucleotides in the DNA, namely the two purines: adenine (A) and guanine (G) and the two pyrimidines: thymine (T) and cytosine (C). A typical DNA double helix is shown in Figure 1.

The symbolic to numerical map is clearly not unique, typically has a biological interpretation, and needs to preserve the specific structure of the DNA sequence under study. One such popular map is the Voss representation ↦={0,1}, where four binary indicator sequences x_l(n), l∈, are generated with 1 indicating the presence of a nucleotide and 0 its absence [18]. An example of the mapping of a single DNA strand to is shown in Figure 2.

Once the DNA symbolic sequence is mapped into numerical version(s), a set of discrete time sequences are generated and are the numerical equivalence of the DNA sequence. These numerical sequences can then by processed using standard signal processing techniques. In particular, the ST-DFT for each elementary sequences can be computed as

, where n is the window starting point, R is the amount of window shift, and h(m)=1 for −M + 1≤m≤0 and zero otherwise. If R=1, then, the window slides one nucleotide at a time whereas if R=3, the displacement of the window is on a 3-nucleotide basis. Note that the all-ones function h(m) does not affect the value of X_l(Rn,k). However, it serves as a place holder for other filters that can be used to replace it, as will be shown in the following section. One popular application of the ST-DFT based technique that has received considerable attention in the past is the identification of the period-3 component using the DNA spectrum, defined for R=3 as follows

A number of researchers have advocated the use of the period-3 component to discriminate between coding and non coding regions (see for example [11,13,16,19-23] to name a few) but the subject remains highly controversial as it is successful for certain genes but does not work for others. To better comprehend the underlying reasons behind this disparity in performance, a new multirate DSP model that provides a full understanding of the inner workings of the DNA periodicity has been first proposed in [24], and studied in details in [25]. This model is shown in Figure 3.

This model provides closed form expressions for the DNA spectrum that generalize and unify some of the already existing results in the literature were obtained. One of these expressions in particular clearly shows that the identification of the period-3 component in the DNA spectrum, a signal processing problem, is equivalent to the detection of the nucleotide distribution disparity in the codon structure of a DNA sequence, a genomic problem. The disparity in the nucleotide distribution within the codon structure of a DNA sequence is termed the codon bias. Using this model, the DNA spectrum is completely characterized by a set of digital sequences, termed the filtered polyphase sequences. By processing these sequences, signal processing techniques can potentially have an impact on understanding and detecting biological structures of this nature. From a computational cost perspective, the computation of the DNA spectrum using this model does not require any complex valued operations [26]. This finding is rather surprising given the existence of complex multipliers in the proposed DSP model as clearly illustrated in Figure 3. It is shown that the direct computation of the DNA spectrum using (2) requires essentially double the amount of arithmetic operations compared to the DSP model approach.

It is important, however, to keep in mind that the above conclusions and results were obtained using the Voss symbolic to numerical transformation. A fundamental research issue is to therefore determine the sensitivity of the signal processing based method to the choice of the symbolic to numerical map. In particular, the core question here is: how dependent are the above results on the Voss representation? Are these results invariant with respect to the other popular maps in the literature? Can we derive necessary and/or sufficient conditions for the invariance of the DNA spectrum to the symbolic to numerical transformation? Is there a general mathematical framework that can help us generate new symbolic to numerical maps for which the DNA spectrum remains essentially the same? These are the type of questions we address in this article and provide answers to. One approach to answer this question was presented in [27], where a novel framework for the analysis of the equivalence of the mappings used for numerical representation of symbolic data based on signal correlation was presented, along with strong and weak equivalence properties. In [28], we attempted to answer the same question starting at the aforementioned DSP model for a limited set of mappings. Our main goal in this study is to de-embed the symbolic to numerical mapping process from the DNA spectrum computation process. We answer a set of other relevant questions along the way.

A key remark is in order at this point: while the DSP model approach proposed in Figure 3 has many advantages, it is not well suited for investigating the role of the symbolic to numerical map in the identification of DNA harmonics. It follows that a completely new paradigm for detecting DNA harmonics is required. The main contribution of this article is therefore the derivation of a novel matrix-based framework for the computation of the DNA spectrum that is extremely well fitted to the study of the symbolic to numerical transformation. Specifically, we first derive a new matrix-based expression of the DNA spectrum that:

1. comprises most of the existing mappings in the literature as special cases,

2. shows that the DNA spectrum is in fact invariable under all these mappings,

3. generates a necessary condition for the invariance of the DNA spectrum to the symbolic to numerical mapping used to compute it.

Furthermore, the new algebraic framework presented here decomposes the frequency identification problem into several fundamental components that are totally independent of each other. It follows that sophisticated digital filters and/or alternative transformations to the DFT such as the discrete cosine, sine, and Hartley transforms can always be easily incorporated in the harmonics detection scheme irrespective of the choice of the symbolic to numerical map. Finally, although the newly proposed framework is matrix based, we show that similar to the DSP model approach, the computation of the DNA spectrum using this new framework is very efficient.

The article is organized as follows. In Section 2, we derive a new matrix based framework to efficiently compute the ST-DFT-based spectrum. New expressions for the ST-DFT and its magnitude squared are obtained and indicate that these quantities are completely parameterized by some pre-defined matrices. The numerical values of these matrices simply depend on our choice of filtering (e.g., rectangular window versus non-rectangular one versus general FIR filters) as well as our choice of data transform (e.g., the DFT versus the DCT versus the DST).

Using these results, in Section 3, a new expression of the DNA power spectrum is derived and is also completely defined by these matrices. The elegance of this matrix based approach is that it allows the incorporation of general symbolic to numerical maps into the newly derived DNA spectrum expression provided these generic maps can be expressed as affine transformations of the Voss representation. This last assumption is motivated by the fact that all the popular maps that are available in the literature satisfy the affine condition. Furthermore, the maps are now completely characterized by the affine transformation (two matrices A and b) and can be therefore changed without affecting the remaining matrices in the DNA spectrum expression. In conclusion, the newly derived DNA spectrum expression is stated as a function of a number of matrices. Each of these matrices captures an essential component of the process (filtering, data transform, symbolic to numerical map) and the elements of each matrix can be changed without affecting the other matrices.

In Section 4 and using the above results, we show that the Voss-based DNA spectrum is essentially invariant under some of the most popular maps in the literature. A necessary and sufficient condition for the invariance of the DNA spectrum under any map is also derived.

In Section 5, we show how the special structure of the filtering matrix allows the efficient use of sophisticated digital filters to improve the detection performance of DNA harmonics through the computation of the DNA spectrum. We also show how to replace the DFT by other fast transforms such as the discrete cosine transform (DCT), the discrete sine transform (DST), and the discrete Hartley transform (DHT). Finally, some concluding remarks are mentioned in Section 6. A list of the different notation used in the article is summarized in Table 1.

Given a sequence x(n) of length N, the ST-DFT is typically implemented using a sliding window approach as shown in Figure 4. Windows of length M that overlap with a factor R are first generated to form x_r(n),r=1,2,…,N_w, where N_w=⌈(N−M + 1)/R⌉ is the number of resulting windows. Once we map the DNA sequence into an integer number of numeric sequences γ, given by x_l(n), l=1,…,γ(), the ST-DFT’s X_l(n), l=1,…,γcan be found and their squared magnitudes are added to result in the DNA Spectrum S(n) as summarized in Figure 5.

It was shown in [26] that the ST-DFT of x(n) can be written as

where the quantities X_r(n),∀ r∈{0,1,…,R−1} are the so-called filtered polyphase sequences given by

∀ r∈{0,1,…,R−1}. The impulse response h_r(m) is the inverse -transform of H_r(z) in Figure 3. Equations (3) and (4) can be used to compute the ST-DFT of a discrete time sequence, and subsequently its magnitude squared. In this section, we re-express these equations in matrix form, and then use the new formula to derive an expression for . Throughout the article, vectors and matrices (arrays) are always expressed in bold letters. The notation for the various matrix operations is given in Table 2.

Matrix formulation of the ST-DFT

Using the defined matrix notation, we can restate Equation (3) as

The real valued array

is the vector whose elements are the R filtered polyphase components. Similarly, the complex valued R-element array

is the vector whose elements are the R equispaced phasors located on the unit circle with phase deviations as shown in Figure 6 for R=3 and R=8. Note that

∀R≠1, which implies that the sum of elements in C is equal to 0. This is a key feature of the complex array C that will be used in later sections to simplify important expressions.

On the other hand, we observe that (4) can be written in the following matrix format

∀ r∈{0,1,…,R−1}, where h is an all-one vector of length M/R, and of length M/R is the r^thpolyphase component of the window x(n) of length M. Using (9), the R filtered polyphase components X_r(n) can be arranged in the following array format

Using the identity

it follows that

which can be restated in matrix format as

where H≐I_R⊗h^T is an R×R matrix of blocks, given by

The window of length M is a block interleaved version of the sliding window x(n) of length M starting at index n. Generating can be accomplished by blocking the window x(n) into an array of R elements per row (hence M/R rows), and then reading the array out column by column. The ST-DFT can therefore be completely identified as a function of C, h, and as follows

The complex row vector C^TH is an array of R blocks, each of length as given by

which represents M/R repetitions of the elements in C. Similar to C, the sum of elements in C^TH is equal to 0.

A matrix based expression for the magnitude squared of the ST-DFT

Using (5), the magnitude squared of the ST-DFT can be expressed as

where matrix D≐C^⋆C^T is an R×Rmatrix given by

D is obviously a right circulant (hence Toeplitz) matrix whose rows and columns are rotated versions of C. Obviously, the sum of any row or column elements in D is equal to 0. Substituting (12) in (14), or equivalently using (13), implies that the spectrum S(n) can be stated as

where

is an R×R matrix of blocks, given by

Matrix W can be represented as a Kronecker product of D and an all-one matrix. Note that any row or column in W is a rotated version of C^TH, therefore, the sum of the elements of any row or column in Wis equal to 0.

A first step towards finding the DNA spectrum S(n) is the symbolic to numeric mapping as was shown in Figure 5. Once the symbolic DNA sequence is mapped into γnumeric sequence(s), the short-time discrete Fourier transform is applied to each of them and the sum of the squared magnitudes of the ST-DFTs will result in the DNA spectrum at the frequency point as given by

For simplicity, we denote as S(n) in the following sections. Several mappings were introduced in the literature using both real and complex numerical values with typical number of sequences γ=1 up to 4 to maintain reasonable computation complexity. In this section, we use the results of Section 2 to derive general expressions for the M/R ST-DFT and spectrum for any symbolic to numeric mapping.

The Voss-based DNA spectrum

The simplest and most commonly used map of a DNA sequence is the Voss representation : that is to form γ=4 binary indicator sequences x_A(n),x_C(n), x_G(n), and x_T(n) where a 1 would indicate the presence of a base and 0 indicates its absence [18]. This approach has been extensively used in relevant genomic research. Note that the four sequences are not linearly independent since for any index n, the four sequences will add up to one. That is

This redundancy plays an important role in the derivations of this section. Moreover, it follows that for any length-M window starting at n, the four mapped Voss windows will add up to an all-one length-M sequence and the same fact holds for the interleaved windows

For illustration, Figure 7a shows a sample DNA window that is mapped into the corresponding numeric windows in Figure 7b,d,f,h. With an example interleaving factor R=3, the interleaved windows are shown in Figure 7c,e,g,i. Each of the four sequences is a discrete time sequence that can be processed using the analysis of Section 2.

Therefore, the ST-DFT of each sequence can be found using (13) to be

, and the power spectrum of each sequence can hence be derived as in (16) to be

. It follows that the Voss-based DNA spectrum S_v(n) is

An obvious step at this point is to simplify (19) to avoid the summation over different bases. To do this, we use Equation (18) to arrange the ST-DFT’s of x_l(n), in the following format

Using (11), it follows that

We define ϒ_v(n): the array of the four Voss-based ST-DFTs. It can now be written as

where I₄ is the 4×4 identity matrix, and the vector of length 4M is an array of the four Voss interleaved windows starting at index n: . Using the identity

the Voss-based DNA power spectrum can be manipulated into

In (26), I₄ and W are constant matrices ∀n. Hence the computation of the spectrum S_v(n) for different windows of a DNA sequence needs only the evaluation of the Voss interleaved array .

Computing the DNA spectrum under general symbolic to numerical maps

Similar to the Voss representation case, any map of γsequences can be processed using the analysis of Section 2. It directly follows that the ST-DFT and spectrum of a single sequence are given by

where l=1,2,…,γ. The array of -mapped ST-DFTs ϒ_d(n) is therefore given by

The -based DNA spectrum can easily be shown to be

where the vector of length γM is an array of the γ-mapped and interleaved windows starting at index n: . It is clear that for every different map , a new interleaved windows array has to be evaluated in order to compute a spectrum point S_d(n). In this following, we introduce a different new approach to recompute (28) without updating for every map. Basically, we derive a new expression for S_d(n) in terms of and a new constant matrix so that we incorporate the map dependance in the matrix part rather than the interleaved array part. In other words, since the map is already well-defined, we use the map to complete the chain and hence find the spectrum S_d(n). Consider the following affine transformation from Voss sequences to a general array of -mapped sequences

where A_γ×4and b_γ×1=[b₁b₂ … b_γ]^Tare constant possibly complex valued arrays. It follows that the array of the -mapped interleaved windows can be written in terms of the array the Voss-mapped interleaved windows in the following form

where defined as

is an array of γM-element blocks, each block is M repetitions of one element of b. Substituting for in (24) results in a new formula for the array of -mapped ST-DFTs ϒ_d(n) into

An important result at this point is that the second term in ϒ_d(n) is actually equal to 0. This can be verified by reducing it into the following form

Recall that the sum of elements in C^THis equal to 0. Therefore, since is a constant vector, the product is equal to 0, ∀ l=1,2,…,γand hence

The ST-DFTs array ϒ_d(n) can therefore be simplified using the Kronecker product identity (22) into

It follows that the -based DNA spectrum S_d(n) is

where B≐A^HA. Equation (35) indicates that when a certain symbolic to numeric mapping is used, the DNA power spectrum S_d(n) is completely defined in terms of the Voss-based interleaved array along with constant matrices W and B which is a function of the transformation matrix A (). Note that if A=I₄then B=I₄ at which (35) reduces to (26) which is the Voss-based spectrum case.

The results found in Section 3 can be applied to some mappings that are widely used in the literature. In specific, by defining the corresponding transformation matrices A and B (), closed form expressions for S_d(n) are obtained. Furthermore, for a number of mappings, we show that the -mapped spectrum S_d(n) is in fact a scaled version of the Voss-based spectrum S_v(n).

Four-to-four (γ=4) representations

In this scheme, each Voss sequence is scaled by a possibly complex coefficient according to the following transformations matrices

where a, a, g, and t are real or complex coefficients used to scale x_A(n),x_C(n),x_G(n), and x_T(n), respectively. The corresponding array of ST-DFT’s ϒ_d(n) is subsequently given by

and the DNA spectrum S_d(n) is

Now, we extend this result to certain transformations where numeric values of the scale factors a, a, g, and t are specified.

§Tetrahedral mapping.

The so-called tetrahedral representation has been proposed in [13,29]. In this mapping scheme, the four nucleotides are represented by four equal length vectors oriented towards the corners of a tetrahedron. Projecting the basic tetrahedron on a plane will reduce the dimensionality of the representation to two. This mapping can be defined by the mapping matrix

It can be easily seen that in this case: which implies that B=2I₄. The corresponding DNA spectrum is

Since B=αI₄(α=2), the tetrahedral-based DNA spectrum is a scaled version of the Voss-based spectrum.

§Quaternion mapping.

A more involved step is to replace the complex number set of the tetrahedral mapping with its algebraic generalization, the set of quaternions. Quaternions have been used to map DNA sequences [30] and are simply defined as hypercomplex numbers given by , where i,j,k are complex coefficients such that i²=j²=k²=ijk=−1 and . The transformation matrix is given by

In this case, . The corresponding DNA spectrum is

§Higher order mappings.

An alternative Quaternion transformation is given by A=diag(1 + i + j + k,1 + i−j−k,1−i−j + k,1−i + j−k), which results in B=4I₄ and consequently S_d(n)=4S_v(n). In general, for a complex representation system with η dimensions and equal amplitude coefficients: B=ηI₄and hence the spectrum S_d(n)=ηS_v(n).

Four-to-three (γ=3) mappings

In order to reduce the DNA spectrum computational cost, several mappings have been proposed with smaller numbers of sequences.

§-curve mapping.

One such important symbolic-to-numeric map is the -curve mapping [24], which is a unique 3-dimensional curve representation whose sequences have values 1 and −1. One advantage of the -curve mapping is that each of its three sequences has a biological interpretation. This scheme is given by

Therefore, the transformation matrices are

Matrix B in this case can be written as

Note that the term involving B₁in S_d(n) can be manipulated into

Recall from (19) that . Take also into consideration that the sum of elements of any row or column in W is equal to 0. This implies that , at which it is easy to see that . Similarly, . Therefore, only the first term in B contributed to S_d(n) at which the -curve mapped DNA spectrum is a scaled version of the Voss-based DNA spectrum

This ratio is consistent with the result we first derived in [24] for R=3, but is now shown to be general for any value of R. We are now ready to state an important result.

Theorem.Necessary and Sufficient condition for the invariance of the DNA spectrum. Consider the following affine transformation from Voss sequences to a general array of -mapped sequences

where A_γ×4 and b_γ×1=[b₁b₂ … b_γ]^T are constant possibly complex valued arrays. Define the 4×4 matrix B=A^HA. The DNA spectrum is invariant under this map, i.e., S_d(n)=αS_v(n) if the transformation matrix B can be written as , where B_i holds constant rows and/or constant columns ∀ i.

The proof follows by simply observing that if B_ihas constant rows and/or constant columns, then . We remind the reader at this point that the vector b_γ×1has no bearing on the invariance of the DNA spectrum.

§Simplex mapping.

The simplex mapping is essentially another tetrahedron structured mapping that aims to eliminate the computational redundancy. Its transformations matrices are

Matrix B in this case can be written as

Similar to the -curve case, . It follows that the simplex-based DNA spectrum is also a scaled version of the Voss-based spectrum, and is given by

This ratio is consistent with the result in [31] which was limited to direct DFT and is now shown to be extended to M/R ST-DFT with any value of R.

Four-to-two (γ=2) mappings

Pairing couples of nucleotides together was proposed in the literature in order to exploit certain biological features in addition to complexity reduction. For example, it was suggested that exons are rich in nucleotides C and G, while introns have more A and T[29]. This claim inspired the transformation

It is obvious that the DNA spectrum in this case can be simplified to

which obviously is not a scaled version of S_v(n) since B in this case can not be written as , where B_i holds constant rows and/or constant columns ∀ i.

Four-to-one (γ=1) mappings

Single sequence representations can be generated by assigning each nucleotide a certain coefficient [4,13] in order to keep the single sequence structure using the transformation array and matrix

Note that the coefficients chosen for the tetrahedral, quaternion, and paired coupled mappings can be reused along with the single sequence formulation. For example, the paired couples case can be reformulated in a single sequence of 1’s and −1’s using and

at which the DNA spectrum is

Similar to the previous case, S_d(n) is not a scaled version of S_v(n).

Experimental verification

To briefly verify the results of this section experimentally, we apply Equation (35) to real DNA sequences, when the Voss, tetrahedral, quaternion, Z-curve, and simplex maps are employed. For comparison with previous study, we consider first the DNA sequence F56F11.4 in the C. elegans chromosome III. This sequence is 8060 nucleotides and has been used as a benchmark by many researchers [13] to extract the periodicity component at R=3. The DNA spectra at R=3 are shown in Figure 8 for the five former mappings, and are obviously related by the constant scale factors derived earlier in the section which clearly verifies our results. Although we lack the space for more general simulations, it is important to state that all the spectra relations are maintained experimentally at other values of R associated with higher order periodicities.

For generality purposes, we test two more sequences extracted from the well known Burset-Guigo database [32]. In specific, DNA spectra at R=3 of the zeta globin gene (ECZGL2) of length 1563, and the Alouatta seniculus epsilon-globin gene (ALOEGLOBIM) of length 1691 are shown in Figures 9 and 10, respectively, for the five former mappings. It can be seen that the relations are still preserved.

Alternative DNA periodicity measures using fast data transforms [33-35], wavelets, and finite impulse response (FIR) digital filters [25,36] were recently proposed to improve the detection performance of these periodicities. However, each method was obtained separately from the other using seemingly a different approach. In this section, we show that our proposed framework can systematically generate all these results by simply changing a number of matrices. It therefore provides a generic unified framework for generating alternative measures of DNA periodicities. For example, we can re-express the matrices D and W in terms of general digital filters and use these filters to modify (35) in order to generate new spectrum formulas. Furthermore, using symmetry based decompositions of D and W, we simplify (35) into a formula with low computational complexity.

Modified periodicity measures

Recall from Section ‘2’ that matrix W is given by

Obviously, W is completely defined by the real array h and the generally complex array C. Note that h and C can be viewed as the impulse responses of two FIR filters defined by the z-transforms H(z) and C(z).

Updating the real filter h

The FIR filter H(z) is the standard rectangular window filter and has a low pass frequency response with a −13 dB attenuation. To improve its filtering performance, we can use a more general FIR filter, denoted by and expressed as

which is the -transform of the general array given by

From a signal processing perspective, achieving better performance can be obtained by replacing the rectangular window with another one, , that has slightly wider main lobes but much more attenuated side lobes, as shown in Table 3. The impulse responses of such windows are depicted in Figure 11a for R=8 and M=96. Better harmonics characterization can be achieved by giving each nucleotide position within the window a relative weight in contrast to the rectangular where equal weighting is given to all nucleotides. It turns out that the Blackman window has the best main-to-first side lobe attenuation behavior as shown in Figure 11b compared to the rectangular window case and therefore provides the best smoothing of the DNA spectrum.

By replacing h with , the matrix H can be in turn expressed as

and the complex row vector is now given by

It can be easily seen that the sum of elements in is still equal to zero as was the case for C^TH. Consequently, it follows that the sum of any row or column in is still equal to zero. This is a fundamental result which, in turn, implies that all the derivations of Section 3 are still the same even when replaces h. In particular, the -based DNA spectrum and the -based one can be stated as

Moreover, all the mathematical relations derived in Section 3 between the -based spectrum and the Voss-based one are all still valid even when h is replaced by .

Experimental verification

To experimentally verify this result, we consider finding the DNA spectrum of the three DNA sequences used in the previous section when is set to a Blackman window. The relations between the spectra when using the Voss, tetrahedral, quaternion, Z-curve, and simplex mappings are still the same as shown in Figures 12, 13, and 14.

Updating the complex filter C

Similar to H(z), the FIR filter C(z) can be replaced by a more sophisticated filter expressed as

which is the -transform of the general array given by

Note that, in this case, the elements in array do not necessarily add to zero anymore. Consequently, the sum of elements in any row or any column in or is not necessarily zero. We also note that unlike the case of , using instead of C keeps the spectrum formulas in (58) correct but does not preserve the mathematical relations between the different -mapped spectra and the Voss-based spectrum.

Joint optimization of and

It should be clear at this point that better DNA harmonics detection performance can be potentially achieved through a joint “optimization” of and . For example, a learning paradigm can be used with a least-mean-square (LMS) criterion to find the optimal set, and . Alternatively, a biologically induced criterion can yield a substantial boost in performance but it is not clear which criterion to use. This interesting but challenging research topic is however outside the scope of this article and will not be further pursued here.

Example

Standard discrete time transforms have been proposed to replace the ST-DFT in the periodicity detection problem. In particular, the short time discrete cosine transform (ST-DCT), sine transform (ST-DST), and Hartley transform (ST-DHT) were introduced and analyzed for this purpose [33]. In this example, we show that these three transforms fit naturally within our proposed analysis when the two arrays and are adjusted correctly for each case. Although these standard transforms are not optimized for certain data sets, they can serve as preliminary tests for better periodicity detection. In [33], the short time DFT, DCT, DST, and DHT at k=M/R where shown to be given by

where t∈{fcsh} indicates Fourier, cosine, sine, and Hartley transforms, respectively, are possibly complex coefficients, and h^(t)(m)=(α^t)^m. Values of the parameters α, a, b, and θ_rfor every transform are adjusted according to Table 4. For illustration, setting α=1, a=1, b=0, and θ_r=−2Πr/R in (59) results in the ST-DFT case. An efficient implementation to calculate Equation (59) is shown in Figure 15 which generalizes Figure 3.

This model provides a general framework that encapsulates the computation of the short-time Fourier, cosine, sine, and Hartley transforms at frequency point k=M/R. Therefore, the same matrix-based analysis of Sections 2 and 3 can be used. Matrix W will be updated into

and therefore the -based DNA spectrum when one of the ST- DFT, DCT, DST, or DHT is employed can be stated as

where the values of and are adjusted according to Table 5.

Note that similar to the Fourier case, the sum of elements in for the cosine and Hartley transforms cases is equal to zero. Therefore, under these two cases, the relations between different -based DNA spectra and the -based DNA spectrum are still the same as given in Section 3.

At this point, it can be concluded that the -based DNA spectrum is completely defined in terms of the Voss-based array of interleaved windows , the mapping matrix A, the real array , and the generally complex array . This conclusion is summarized in Figure 16.

A real approach for the spectrum computation

A real computationally-efficient alternative for the evaluation of S_d(n) can be found by observing the special properties of the circulant/toeplitz matrix D or equivalently the block matrix W. We use the fact that for a generally-complex matrix Q: y^HQy=0, , if Q is an antisymmetric matrix. We start by splitting Dinto its symmetric and antisymmetric parts

where D_s is a circulant and Toeplitz real R×Rmatrix given by

and D_as is a circulant and Toeplitz complex R×Rmatrix given by

Substituting for D in (15), we get a simple form of the spectrum S(n)

where , . The block matrix

is an R×R matrix of blocks. Using (63) to update the DNA spectrum (21), S_v(n) simplifies into

Following the same analysis of Section 3, (64) can be easily manipulated into a more elegant completely real form given by

or more generally, (35) can be updated into

which provides a completely real approach for the computation of the -mapped spectrum S_d(n). Note that all results and different spectra relations in Section 3 still hold when W_s replaces W as in (65).

Computational complexity comparison

To quantify the computational credit of this real approach, we compare the complexity of (63) to that of (16) of a single discrete time sequence. Since can be complex as well according to the mapping used, we find the number of real multiplications and additions needed to evaluate (63) when each of and W is either real or complex, as given in Table 6. Recall that the multiplication of the complex numbers x and y, where x=a + jb and y=c + jd requires the computation of ac−bdand ad + bc, which requires four real multiplications and two real additions.

Example

For illustration, we evaluate the spectrum S_v(n) using W_swhen R=3, and compare the result to the formula derived in [37]. In specific, we use (64) to find the spectrum S(n) as follows

Expanding and completing the square, it follows that

where q=(r + 1) mod 3. The matrix-based DNA spectrum formula in (67) is consistent with the result derived using a different approach in [37].

In this article, we have introduced a matrix based framework for locating hidden DNA periodicities using spectral analysis techniques that are invariant to the choice of the symbolic to numerical map. The primary advantage of the presented approach over some of the previous study is the decomposition of the spectrum expression into key matrices whose values can be set independently from each other. Each matrix represents one of the essential components involved in the computation of the spectrum such as the symbolic to numerical map, the data transform, and the filtering scheme. The above framework is derived under the assumption that the symbolic to numerical map can be obtained from the Voss representation using an affine transformation. This assumption is however quite loose given that most (if not all) of the proposed maps in the literature satisfy this requisite. Using the new framework, we have then shown that the DNA spectrum expression is invariant under these maps. We have also derived a necessary and sufficient condition for the invariance of the DNA spectrum in terms of the affine transformation matrix A (the b vector in the affine transformation does not affect the DNA spectrum).

This condition can serve as the basis for generating novel symbolic to numerical map that preserve the DNA spectrum expression. Finally, in the latter sections of the article, we have shown the potential of using different filtering schemes, e.g., windows other than the rectangular one as well as alternate fast data transforms, e.g., the DCT, DST, and the Hartley transform. A number of simulation results that verify the findings of this article and a brief quantitative analysis of the computational complexity of the new approach were given in the same sections. Future research study would consider the optimization of the different building blocks, namely the symbolic to numerical map, the data transform, and the filtering scheme. This, in turn, requires a deep understanding of the biological significance of different DNA periodicities in order to set up a meaningful objective function and appropriate constraints. Ultimately, the framework proposed here can be incorporated in a more sophisticated system to study the complex structure of genomic sequences and understand the functionality of its various components. Finally, this efficient framework can be extended to the analysis of other types of symbolic sequences of various limited alphabets, either biological sequences (such as protein sequences) or even non-biological ones.

The authors declare that they have no competing interests.

TS acknowledges partial support from the NSF via grants DMS 0811169 and DMS-1042939.

1. G Benson, Tandem repeat finder: a program to analyze DNA sequences. Nucleic Acids Res 27(2), 573–580 (1999). PubMed Abstract | Publisher Full Text | PubMed Central Full Text

2. J Butler, Forensic DNA Typing: Biology and Technology behind STR Markers (Academic Press, MA, Burlington, 2003)

3. CA Cummings, DA Relman, Microbial forensics: cross-examining pathogens. Science 296, 1976–1979 (2002). PubMed Abstract | Publisher Full Text

4. P Ramachandran, W Lu, A Antoniou, Filter-based methodology for the location of hot spots in proteins and exons in DNA. IEEE Trans. Biomed. Eng 59(6), 1598–1609 (2012). PubMed Abstract | Publisher Full Text

5. T Strachan, AP Read, Human Molecular Genetics (John Wiley and Sons, New York, 1999) (http://www, 1999), . ncbi.nlm.nih.gov/books/NBK7580/ webcite

6. AF Smit, The origin of interspersed repeats in the human genome. Curr. Opin. Genet. Dev 6, 743–748 (1996). PubMed Abstract | Publisher Full Text

7. DC Rubinsztein, MR Hayden, Analysis of TRIPLET REPEAT Disorders (Bios Scientific Pub, Oxford, England, 1999)

8. R Gupta, A Mittal, S Gupta, An efficient algorithm to detect palindromes in DNA sequences using periodicity transform. Signal Process 18(4), 8–20 (2001). Publisher Full Text

9. VR Chechetkin, AY Turygin, Size-dependence of three-periodicity and long-range correlations in DNA sequences. Phys. Lett A 199, 75–80 (1995). Publisher Full Text

10. VR Chechetkin, AY Turygin, Search of hidden periodicities in DNA sequences. J. Theor. Biol 175, 477–494 (1995). PubMed Abstract | Publisher Full Text

11. BD Silverman, R Linsker, A measure of DNA periodicity. J. Theor. Biol 118(3), 295–300 (1986). PubMed Abstract | Publisher Full Text

12. D Holste, I Grosse, S Beirer, P Schieg, H Herzel, Repeats and correlations in human DNA sequences. Physic. Rev. E 67(06913) (2003)

13. D Anastassiou, Genomicsignalprocessing Process, IEEE Signal Mag 18(4), 8–20 (2001). Publisher Full Text

14. VR Chechetkin, VV Lobzin, Anticodons, frameshifts, and hidden periodicities in tRNA sequences. J. Biomol. Struct. Dyn 24(2), 189–202 (2006). PubMed Abstract | Publisher Full Text

15. D Anastassiou, Frequency domain analysis of biomolecular sequences. Bioinformatics 16(12), 1073–1082 (2000). PubMed Abstract | Publisher Full Text

16. S Tiwari, S Ramachandran, A Bhattacharya, S Bhattacharya, R Ramaswamy, Prediction of probable genes by Fourier analysis of genomic sequences. Comput. Appl. Biosci 13(3), 263–270 (1997). PubMed Abstract

17. M Akhtar, E Ambikairajah, Time and frequency domain methods for gene and exon prediction in Eukaryotes. Proceedings of ICASSP (International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 2007, Honolulu, Hawaii, USA, 2007), pp. 573–576

18. RF Voss, Evolution of long-range fractal correlations and 1/f noise in DNA base sequences. Phys. Rev. Lett 68(25), 3805–3808 (1992). PubMed Abstract | Publisher Full Text

19. JW Fickett, The gene identification problem: an overview for developers. Comput. Chem 20, 103–118 (1996). PubMed Abstract | Publisher Full Text

20. N Bouaynaya, D Schonfeld, Non-stationary analysis of coding and non-coding regions in nucleotide sequences. IEEE J. Sel. Top. Signal Process 2(3), 357–364 (2008)

21. PP Vaidyanathan, BJ Yoon, Gene and exon prediction using all pass-based filters. Gensips Proc ((Workshop on Genomic Signal Processing and Statistics (GENSIPS), Raleigh, North Carolina, USA, 2003)), pp. 1–4

22. PP Vaidyanathan, B Yoon, Digital filter for gene prediction applications. in Proc, ed. by . Asilomar conference (Asilomar Conference on Signals, Systems and Computers (ACSSC), Pacific Grove, CA, USA, 2003), pp. 306–310

23. M Akhtar, J Epps, E Ambikairajah, On DNA numerical representations for period-3 based exon prediction. Proceedings of the workshop on Genomic Signal Processing and Statistics (International Workshop on Genomic Signal Processing and Statistics (GENSIPS) 2007, Tuusula, Finland, 2007), pp. 1–4

24. A Rushdi, J Tuqan, Gene identification using the Z-curve representation. Proceedings of the 31st IEEE ICASSP conference (International Conference on Acoustics, Speech, and Signal Processing (ICASSP) 2006, Toulouse, France, 2006), pp. 1024–1027

25. J Tuqan, A Rushdi, A DSP approach for finding the codon bias in DNA sequences. IEEE J. Sel. Top. Signal Process 2(3), 343–356 (2008)

26. A Rushdi, J Tuqan, An efficient algorithm for DNA discrete Fourier analysis. Proceedings of the 3rd IEEE Cairo International Biomedical Engineering Conference (CIBEC) (Cairo International Biomedical Engineering Conference, Cairo, Egypt, 2006), pp. 1–4

27. L Wang, D Schonfeld, Mapping equivalence for symbolic sequences: Theory and applications. IEEE Trans. Signal Process 57(12), 4895–4905 (2009)

28. A Rushdi, J Tuqan, The role of the Symbolic-to-Numerical Mapping in the detection of DNA Periodicities. Proceedings of the workshop on Genomic Signal Processing and Statistics (International Workshop on Genomic Signal Processing and Statistics (GENSIPS) 2008, Phoenix, AZ, USA, 2008), pp. 1–4

29. PD Cristea, Conversion of nucleotides sequences into genomic signals. J. Cellul. Mol. Med 6(2), 279–303 (2002). Publisher Full Text

30. AK Brodzik, O Peters, Symbol-balanced quaternionic periodicity transform for latent pattern detection in DNA sequences. in Proceedings of the IEEE Int, ed. by . Conf. on Acoustics, Speech, and Signal Processing (International Conference on Acoustics, Speech, and Signal Processing (ICASSP) 2005, Philadelphia, PA, USA, 2005), pp. 373–376

31. E Coward, Equivalence of two Fourier methods for biological sequences. J. Math. Biol 36, 64–70 (1997). Publisher Full Text

32. M Burset, R Guigo, Evaluation of gene structure prediction programs. Genomics 34, 353–357 (1996). PubMed Abstract | Publisher Full Text

33. A Rushdi, J Tuqan, Trigonometric Transforms for Finding Repeats in DNA sequences. Proceedings of the workshop on Genomic Signal Processing and Statistics (International Workshop on Genomic Signal Processing and Statistics (GENSIPS) 2008, Phoenix, AZ, USA, 2008), pp. 1–4

34. JA Berger, SK Mitra, J Astola, Power spectrum analysis for DNA sequences. in Proc, ed. by . of the Int. Sym. on Signal Processing and its App (Conference: International Symposium on Signal Processing and Its Applications (ISSPA) 2003, Paris, France, 2003), pp. 29–32

35. D Kotlar, Y Lavner, Gene prediction by spectral rotation measure: a new method for identifying protein-coding regions. Genome Res 13(8), 1930–1937 (2003). PubMed Abstract | Publisher Full Text | PubMed Central Full Text

36. A Rushdi, J Tuqan, The filtered spectral rotation measure. Proceedings of the 40th IEEE Asilomar Conference on Signals, Systems, and Computers (Asilomar Conference on Signals, Systems and Computers (ACSSC) 2006, CA, USA, 2006), pp. 1875–1879

37. S Datta, A Asif, A fast DFT based gene prediction algorithm for identification of protein coding regions. in Proc, ed. by . of the ICASSP (International Conference on Acoustics, Speech, and Signal Processing (ICASSP) 2005, Philadelphia, PA, USA, 2005), pp. 113–116