The self-similarity properties and multifractal analysis of DNA sequences

The discovery of deoxyribonucleic acid (DNA) structure lies as a fundamental topic of 20th century biology and continues to be the favourite question of some contemporary scientists. Deoxyribonucleic Acid or DNA is a long linear polymer, formed by a large number of nucleotides. Each nucleotide contains a phosphate group, a sugar group and a nitrogen base. The four types of nitrogen bases are adenine (A), thymine (T), guanine (G) and cytosine (C). The uracil base (U) is presented in the RNA molecule. The DNA adopts a double helix form, which is a helical structure constituted by two complementary strands of nucleic acids. For this reason, DNA is also called α-helix.

Shortly speaking, DNA is a macromolecule formed by a large number of nucleotides. Thus, the genetic information is stored in the sequence of the bases. Therefore, the order of the bases determines the information available for the RNA production and proteins.

Powerful methods for DNA sequencing have been developed. These methods serve to sequence complete genomes. For example, small genomes from viruses or fungi, and large genomes, such as the human genome, are made up of 3000 million base pairs [1].

In the National Center for Biotechnology Information (NCBI) directories are the databases that contain complete genomes, complete sequences of chromosomes, sequences of mRNAs, and proteins. The importance to analyse the large DNA databases in the Nonlinear Dynamics context is based on the work conducted earlier by Jeffrey [2], who proposed a graphic representation of these databases via an extended chaos game. Other contributions similarly based on a statistical description of DNA sequences take on a more structured form, such as Zu-Guo Yu et al [3], where the generalized dimensions D_q and its derivative, the ’analogous’ specific heat C_q , are calculated for the coding and noncoding length sequences of bacteria.

The Multifractal formalism was originally proposed to study various chaotic models derived from phenomena associated with turbulence [4]. This theory is used as an archetype of fractal measurements composed of interwoven fractal sets characterized by the Holder exponent, and it has become a crucial tool in the analysis of statistical data.

Section 2 is be dedicated to give the theoretical support of the DNA chaos game paraphrasing the ideas of Barnsley about definition of measures on fractals. In section 3 is showed a didactic approach using a pretty small genome applying the ideas of the last section. Sections 4 and 5 involves the use of these ideas to give a multifractal interpretation to the DNA sequences. In section 5, we discuss the coarse grained multifractal theory and using the curdling theorem. The alternative definition of the singularity spectra used in the Chabbra-Jensen algorithm are obtained. In section 6, we show the multifractal spectra for the six-mers of two bacterias, two archaea, the Homosapiens chromosome 21 and a fungus, and discuss the information contained in the singularity spectra for six DNA sequences.

We construct an unitary square Q with its corners labeled with a different basis of the genomic sequence V = {A,C,G,T}; where V indicates the vertices of Q on the cartesian plane which are A = (0,0), C = (0,1), G =(1,1) and T =(1,0). It is applied the Jeffrey’s "chaos game algorithm", and the Chaos Game Representation (CGR) [2] is obtained. It provides an objective meaning for studying DNA sequences in the context of fractal geometry and non-linear science.

The analytical process for the construction of CGR belongs to a DNA organism sequence that begins with the choice of an arbitrary point inside the square Q. Let P(x₀, y₀) be the starting point. For simplicity, we take the middle point of Q, and let V(x_V , y_V ) be any corner of Q. P(x₁, y₁) are the coordinates of a middle point of the segment PV. These coordinates are given by,

The matrix function of this process can be defined by

The next points are produced by successive applications of the same process.

It means that the DNA sequence can be expressed as a sequence of points produced by iterative applications of (3). Notice that the graph of the set of points P_n forms the CGR.

We plot the CGR of three different organisms in Fig. 1. This shows visually that the DNA sequence has some fractal patterns and a self-affine characteristics. Thus, it shows that the DNA sequence is not a random sequence.

Fig. 1

An IFS is associated with a DNA sequence. We construct four affine transformations [5] to explore regions in Q.

It is important to note that there is a relationship between the chaos game algorithm and the Bernoulli mapping. This is possible, thanks to the chaos game properties. The sequence of points P = {P₁P₂ . . .P_n} and the sequence of DNA have a one-to-one relation. The P is unique for each genome. Knowing the coordinates of P_n, the coordinates of P_n−1 can be determined by using a binary representation for the x and y coordinates. It is possible to show that the inverse relationship of (1) is:

Fig. 2

Any point of P contains the information of which are all the bases of DNA sequence up to that point. Since, knowing P_n, all the points prior to it can be calculated, and using each point corresponds to a base of the genome which are the bases before the n −th.

The Hutchinson’s operator of (4) is W(B)=ωA(B)∪ωC(B)∪ωG(B)∪ωT(B)whereB⊂Q. $W\left( B \right)={{\omega }_{A}}\left( B \right)\bigcup {{\omega }_{C}}\left( B \right)\bigcup {{\omega }_{G}}\left( B \right)\bigcup {{\omega }_{T}}\left( B \right)\,\,\text{where}\,\,B\subset Q.$The application of each function to Q provides an square Q_V = ω_V (Q) with side 1/2. The action of Hutchinson’s operator on Q forms the unitary square, because the union of four squares of side 1/2 is again Q. In other words, Q is invariant under W [6].

When we apply W twice to Q, we have

in this way, a grid with 4² squares is obtained„ each one with side (1/2)². For m applications of W, we obtain 4^m sub-squares of side (1/2)^m, which conform Q.

We have considered subregions of Q instead of chaos game points, due to the fact that P₀ is the midpoint of Q, P₁ = ω_V (P₀) is the midpoint of QV1 , P₂ = ωV2 ° ωV1(P₀) is the midpoint of QV2V1 , and so on. Thus, PN = ωVN ◦· · ·◦ωV1(P₀) is the midpoint of QVN···V1.

Then, each DNA sequence generates a large collection of squares of sides [ (1/2),(1/2)2,…,(1/2)N ]. $\left[ \left( {1}/{2}\; \right),{{\left( {1}/{2}\; \right)}^{2}},\ldots ,{{\left( {1}/{2}\; \right)}^{N}} \right].$The squares satisfy the following property:

Consequently, we can observe that the number of monomers V = (A,C,G,T) in a DNA sequence matches the number of points inside Q_V and the number of dimers V₁V₂ matches with the number of points inside QV2V1 . On the whole, the number of R-mers is given by the number of points inside of QVR···V2V1.

The assertion that the DNA sequence is not a random sequence involves the appearance of statistical techniques for the analysis of the different distributions of the DNA sequences.

This section exemplifies how chaos game works. Through a playful model, we propose to illustrate the theoretical support of DNA chaos game. In this context, a small DNA sequence GAATC it is a toy genome.

Lets us take our toy GAATC. As we have pointed out in (2), the sequence points of GAATC are: {P₁ = ω_G(P₀),P₂ = ω_A(P₁),P₃ = ω_A(P₂),P₄ = ω_T (P₃),P₅ = ω_C(P₄)}, we can see its graph or CGR in figure 3(a).

Fig. 3

The application of Hutchinson’s operator to Q and the SFI (4) may be illustrated by the figure 3(b). Therefore, we can show the subregion’s descriptions for this small sequence, where we have obtained sub-squares of side [ (1/2),(1/2)2,(1/2)3,(1/2)4,(1/2)5 ] $\left[ \left( {1}/{2}\; \right),{{\left( {1}/{2}\; \right)}^{2}},{{\left( {1}/{2}\; \right)}^{3}},{{\left( {1}/{2}\; \right)}^{4}},{{\left( {1}/{2}\; \right)}^{5}} \right]$in figure 4. A tiny sequence provides a didactic point of view to motivate the reader to make his or her own games on the computer.

Fig. 4

We must emphasize that, there exist a subtle distinction between two processes. One of them is the functions’ application, it means that the application of (4) is given according to the characters apparition order in the sequence and it is described graphically in the figure 4(a). The other corresponds to the location of sub-squares, where given a sub-square, we can find it, using the property (8). It means that in the localization processes, the sequence must be read in a reverse order as shown in the figure 4(b).

Figure 1 shows that a fractal F could be generated by a certain DNA sequence and F is contained in Q. From the graphical point of view, 21 Human chromosome’s CGR is the practical proof of self-similarity. If we look closely 1(d), we will notice that the upper right quadrant Q_G, is repeated at every subregion on Q but in reduced scales and different tonalities. We see that the tonalities vary in some nontrivial manner; it can be related to a complex measure. Therefore, this is our motivation to use the multifractal theory.

The description of subregions provides a systematic way for the construction of a cover C_R, this is called the optimal cover of F. The C_R is composed of squares Q_R of side λ_R, each Q_R has a statistical measure μ(Q_R) ≠ 0 given by the frequency with which the DNA sequence visits Q_R. The aim of this description is that each Q_R of the cover C_R is associated with a R-mer of the DNA sequence.

Summarizing, we can say that a R-mer is a DNA subsequence of size R. It means, that this subsequence hasR basis, if it is within the DNA sequence, then the probabilistic measure of a R-mer depends on the total number with which this appears in the DNA sequence. The cover C_R gives the number of all the different possible R-mers: 4^R.

In order to get the frequency of each different R-mer of a DNA sequence with N basis, we denote NV1V₂...VR as the number of times R-mer appears in the total sequence. It can be easily noticed that the total number of R-mers in the entire sequence is N −R +1. Therefore the frequency of a given subsequence is

Now, we may say by the fact that the R-mer is inside of the square QV1V2...VR , that is its probabilistic measure μ and its side λ are:

To understand R-mers statistical behaviour more thoroughly, it is necessary to make an intermediate step where we establish a one-to-one relationship between sub-squares of Q and sub-intervals of the unit interval [6].

Thus, we use a quaternary basis associating each nitrogenous base with a number: A = 0, C = 1, G = 2 and T = 3, then the DNA sequence becomes a sequence of digits. This makes it easy the calculation of the frequencies in the environment of a programming language (python 3.6.0). We have plotted the histograms of 21 Human chromosome in figure 5. We see the same behavioral patterns as in a multiplicative process’ measure.

Fig. 5

In this section we use the coarse grained multifractal theory to characterize the different R-mers that are contained in the DNA sequence, as each R-mer has associated a sub-square Q_R of size (1/2)^R , these sub-squares conform the cover C_R = {Q_R} of the fractal F generated by the DNA sequence. Each Q_R has a coarse grained Holder’s exponent, which is given by:

The fractal F is covered by interwoven subsets J_α, which are defined as:

The pre-fractal dimension of each J_α is given by

where N_R(α) is the cardinality of J_α. The multifractal spectra of the cover means the knowledge of the prefractal dimension of all subsets J_α. We determine these spectra using the Chabbra-Jensen algorithm [7].

More details of the algorithm than we present in this paper can be found in [8]. In this algorithm, is fundamental the definition of the escort probability measure of μ(QV1V2...VR) of order q is fundamental, which is given by [9]:

where

The sum over j is taken over all R-mers with a measure different from zero. Taking into account that:

and that Z_R(q) satisfies the power law:

The q-measure of a set J_α is given by:

As the q-measure only takes values between zero and one, then

For all values of λ_R, the q-measure obtains its maximal value for a special value α^∗(q), which satisfies:

We note that (19) and (20) imply that

where the infimum is taken over all possible values of α. This expression establishes that τ(q) and f (q) are a couple of Legendre transforms.

The behaviour of the q-measure around α^∗(q) is found using a logarithmic expansion of (18):

Then, the results of this expansion are:

Therefore, the q-measure for values of α near α^∗ takes the form:

It shows that q-measure concentrates around Jα*(q), $J_{\alpha }^{*}\left( q \right),$and when λR→0,Jα* ${{\lambda }_{R}}\to 0,J_{\alpha }^{*}$is the support of the q-measure, i.e.

then, J_α∗ is the curdling set of the q-measure, and the result (21) is called the curdling theorem [8].

The value of α^∗(q) is found taking the derivative of (20) and using (22).

where in the second term of the r.h.s. of (26), (17) was used. Finally, using (14) and (15), we obtain that:

The fractal dimension of J_α∗ is determined using (20). i.e.

introducing (25) in the last expression, we obtained:

The relations (27) and (29) are the alternative definition of the singularity spectrum used in the Chhabbra and Jensen algorithm. The curve f_R(α) versus α is obtained when the q-parameter is eliminated, but these expressions contain more information about the singularity spectra, as we will discuss in the next section.

We begin discussing the information contained in the Chhabbra-Jensen algorithm, considering different values of the q-parameter. The simplest case is q = 0, when the q-measure is given by:

It corresponds to consider that all the R-mers have the same probability. The value of the fractal dimension of the curdling set, is evaluated using (29), then

This is the fractal dimension of Q which is the support of F , and it is the maximum value of the fractal dimension. As follows from (27), the value of the α^∗(q = 0) is given by:

This is the average of the Holder’s exponent of all the members of the cover C_R. When q = 1, the q-measure reduces to the measure of the Q_R, i.e. μ(q = 1,QV1V2...VR) = μ(QV1V2...VR). From (27) and (29), we have

Thus, for the measure μ(Q_R), the fractal dimension of the curdling set J_α∗ is identical with the value of the Holder’s exponent which characterize this set, but this value is given by the average entropy of the cover C_R of it. We analyse the behavior of μ(q,Q_R) when q goes to infinity. For large values of q, the measure that gives the maximum contribution to Z_R(q) is

Then, for large values of q, we have

Thus, using this result into (14), we obtain

Therefore,

Thus, when it is taking the limit of (27) when q →∞, we obtain

Similarly, using the fact that μmin(QR)=λRαmax, ${{\mu }_{\min }}\left( {{Q}_{R}} \right)=\lambda _{R}^{{{\alpha }_{\max }}},$the behavior of the Holder’s exponent can be obtained when q →−∞, i.e.

The important characteristics of the fractal spectra f_R(α) are the following:

We conclude that there are at least four important values of the Holder’s exponent, α₀,α₁,α_min and α_max, which characterize globally the fractal spectra. The value α₀ where f_R(α) reaches its maxima; the value of α₁ which determines the fractal dimension of the curdling set, the other two values are used for determining the width of the spectrum W = α_max−α_min, and the skew shape of the spectrum r = (α_max−α₀)/(α₀−α_min). The skew parameter r determines which fractal exponents are dominant. Right skew shape of the spectrum with r > 1 is more complex than the left skew shape with r < 1 [10, 11, 12].

The geometrical shape of the multifractal spectra illustrates the level of multifractality of the DNA sequences [10]. In figure 6, we show the multifractal spectra of the sixth-mers of several DNA sequences, all of them were constructed with the range of q’s between q_min = −30 and q_max = 30; with steps Δq = 0.5.

Fig. 6

The values of the quantities that we select for evaluating the "randomness" are reported in Table 1. The maximum of all curves is 2; which is due to that all take the square Q as the support of F. The parameter r shows that the Halobacterium salinarum’ spectrum has a greater symmetry than the others, and the rest have a right-skewed shape. The bacteria have a higher value of W and r, while the archaea have a lower value of them. This suggests that the multifractal spectrum of bacteria presents a greater complexity than the spectrum of archaea. This suggests that the use of the multifractal technique could be used for a quantitative classification of archaea and bacteria.