The common origin of symmetry and structure in genetic sequences.

Biologists have long sought a way to explain how statistical properties of genetic sequences emerged and are maintained through evolution. On the one hand, non-random structures at different scales indicate a complex genome organisation. On the other hand, single-strand symmetry has been scrutinised using neutral models in which correlations are not considered or irrelevant, contrary to empirical evidence. Different studies investigated these two statistical features separately, reaching minimal consensus despite sustained efforts. Here we unravel previously unknown symmetries in genetic sequences, which are organized hierarchically through scales in which non-random structures are known to be present. These observations are confirmed through the statistical analysis of the human genome and explained through a simple domain model. These results suggest that domain models which account for the cumulative action of mobile elements can explain simultaneously non-random structures and symmetries in genetic sequences.

Introduction

Compositional inhomogeneity at different scales has been observed in DNA since the early discoveries of long-range spatial correlations, pointing to a complex organisation of genome sequences[1-3]. While the mechanisms responsible for these observations have been intensively debated[4-9], several investigations indicate the patchiness and mosaic-type domains of DNA as playing a key role in the existence of large-scale structures[4,10,11]. Another well-established statistical observation is the symmetry known as “Second Chargaff Parity Rule”[12], which appears universally over almost all extant genomes[13-15]. In its simplest form, it states that on a single strand the frequency of a nucleotide is approximately equal to the frequency of its complement[16-20]. This original formulation has been later extended to the frequency of short (n ≃ 10) oligonucleotides and their reverse-complement[20-22]. While the first Chargaff parity rule[23] (valid in the double strand) was instrumental for the discovery of the double-helix structure of the DNA, of which it is now a trivial consequence, the second Chargaff parity rule remains of mysterious origin and of uncertain functional role. Different mechanisms that attempt to explain its origin have been proposed during the last decades[19,24-27]. Among them, an elegant explanation[27,28] proposes that strand symmetry arises from the repetitive action of transposable elements.

Structure and symmetry are in essence two independent observations: Chargaff symmetry in the frequency of short oligonucleotides (n ≃ 10) does not rely on the actual positions of the oligonucleotides in the DNA, while correlations depend on the ordering and are reported to be statistically significant even at large distances (thousands of bases). Therefore, the mechanism shaping the complex organization of genome sequences could be, in principle, different and independent from the mechanism enforcing symmetry. However, the proposal of transposable elements[29,30] as being a key biological processes in both cases suggests that these elements could be the vector of a deeper connection.

In this paper we start with a review of known results on statistical symmetries of genetic sequences and proceed to a detailed analysis of the set of chromosomes of Homo Sapiens. Our main empirical findings are: (i) Chargaff parity rule extends beyond the frequencies of short oligonucleotides (remaining valid on scales where non-trivial structure is present); and (ii) Chargaff is not the only symmetry present in genetic sequences as a whole and there exists a hierarchy of symmetries nested at different structural scales. We then propose a model to explain these observations. The key ingredient of our model is the reverse-complement symmetry for domain types, a property that can be related to the action of transposable elements indiscriminately on both DNA strands. Domain models have been used to explain structures (e.g., the patchiness and long-range correlations in DNA), the significance of our results is that it indicates that the same biological processes leading to domains can explain also the origin of symmetries observed in the DNA sequence.

Results

Statistical Analysis of Genetic Sequences

We explore statistical properties of genetic sequences s = α1α2 … α, with α ∈ {A, C, T, G}, by quantifying the frequency of appearance in s of a given pattern of symbols (an observable X). For instance, we may be interest in the frequency of the codon ACT in a given chromosome. More generally, we count the number of times a given symbol α0 is separated from another symbol α1 by a distance τ1, and this from α2 by a distance τ2, and so on. The case of ACT corresponds to α0 = A, α1 = C, α2 = T, τ1 = 1, τ2 = 1. We denote a selected finite sequence of symbols, and by a sequence of gaps. For shortness, we denote this couple by and the size of the observable X by . The frequency of occurrence of an observable X in the sequences s is obtained counting how often it appears varying the starting point i in the sequence:where . As a simple example, for the choice of X = ((A, C, G), (1, 2)) in the sequence s = GGACCGGCCACAGGAA we have N = 16, N′ = 13, and P(X) = 2/13. All major statistical quantities numerically investigated in literature can be expressed in this form, as we will recall momentarily.

The main advantage of the more general formulation presented above is that it allows to inspect both the role of symmetry (varying ) and structure (varying scale separations ) and it thus permits a systematic exploration of their interplay. We say that a sequence has the symmetry S at the scale if for any observable X with length we have, in the limit of infinitely long s,where S(X) is the observable symmetric to X.

We start our exploration of different symmetries S with a natural extension to observables X of the reverse-complement symmetry considered by Chargaff. The reverse complement of an oligonucleotide α1α2, …, α of size n is , where (e.g., the reverse-complement of CGT is ACG). For our more general case it is thus natural to consider that the observable symmetric toisThis motivates us to conjecture the validity of an extended Chargaff symmetryThis is an extension of Chargaff’s second parity rule because X may in principle be an observable involving (a large number of) distant nucleotides and thus equation (4) symmetrically connects structures even at large scales. One of the goals of our manuscript is to investigate the validity of Eq. (4) at different scales, which will be done by choosing observables X of size of up to millions of base pairs.

By combining P(X) of different observables X this extended Chargaff symmetry applies to the main statistical analyses already investigated in literature, unifying numerous previously unrelated observations of symmetries. As paradigmatic examples we have:

the frequency of a given oligonucleotide can be computed as P(X) with the choice and τ = (1, 1, … 1). Equation (4) implies that the frequency of an oligonucleotide is equal to the frequency of its reverse-complement symmetric and thus implies the second Chargaff parity rule, a feature that has been extensively confirmed[20-22] to be valid for short oligonucleotides . We report few examples of frequencies of dinucleotides () in human chromosome 1: P(AG) = 7.14% ≈ P(CT) = 7.13% ≠ P(GA) = 6.01% ≈ P(TC) = 6.01%, in agreement with symmetry (4). Note that, the validity of Second Chargaff Parity rule at small scales ( for dinucleotides) is not enough to enforce equation (4) for generic observables X (e.g., of size );

the autocorrelation function C(t) of nucleotide ω at delay t is the central quantity in the study of long-range correlations in the DNA. It corresponds to the choice and . Equation (4) predicts the symmetry . In the specific case of dinucleotides, such relation has been remarked in ref.[31]. Our result holds for any oligonucleotide ω;

the recurrence-time distribution R(t) of the first return-time between two consecutive appearances of the oligonucleotide ω is studied in refs[32,33]. By using elementary arithmetic and common combinatorial techniques, it is easy to see that R(t) can be in fact written as a sum of different P(X). Equation (4) hence predicts . This symmetry was observed for oligonucleotides in ref.[34].

This brief review of previous results shows the benefits of our more general view of Chargaff’s second parity rule and motivates a more careful investigation of the validity of different symmetries at different scales .

Symmetry and structure in Homo Sapiens

We now investigate the existence of new symmetries in the human genome. In order to disentangle the role of different symmetries at different scales we construct a family of observables for which we can scan different length scales by varying the gaps vector . Particularly useful is to fix all gaps in but a chosen one τ, and let it vary through different scales. To be more specific consider the following construction: given two patterns and we look for their appearance in a sequence, separated by a distance . This is equivalent to look for composite observable or, for simplicity, . We consider two patterns X, X of small (fixed) size and we vary their separation to investigate the change in the role of different symmetries.

To keep the analysis feasible, we scrutinize the case where X and X are dinucleotides separated by a distance from each other. This goes much beyond the analysis of the frequencies of short oligonucleotides mentioned above because with ranging from 1 to 107 we span ranges of interests for structure and long-range correlations. This choice has two advantages: by keeping the number of nucleotides in each X small we improve statistics, but at the same time we still differentiate from the simpler complement transformation.

In order to compare results for pairs X, X with different abundance, we normalise our observable by the expectation of independence appearance of X, X obtainingDeviations from z = 1 are signatures of structure (correlations). Table 1 shows the results for chromosome 1 of Homo Sapiens, using a representative set of eight symmetrically related pairs of dinucleotides at a small scale . The results show that z is significantly different from one and that Chargaff symmetric observables appear with similar frequency, in agreement with conjecture (4). Figure 1 shows the same results of the Table varying logarithmically the scale from up to (more precisely we use with i ∈ {0, 1, 2, .., 24}). At different scales we see that a number of lines (observables X, X) coincide with each other, reflecting the existence of different types of symmetries.

Table 1

Chargaff symmetric observables appear with similar frequency in the human chromosome 1.

X = XA **** XB	P(X)	z_[X_A,X_B](ℓ=4)
CC **** TCGA **** GG	0.004010.00404	1.2361.242
GG **** GATC **** CC	0.003660.00366	1.1271.128
GG **** TCGA **** CC	0.003020.00299	0.9290.922
CC **** GATC **** GG	0.002650.00264	0.8180.813

Each line contains an observable X constructed combining oligonucleotides α1α2 … α8 where α1α2 equal to X, α3α4α5α6 are arbitrary (any in {A, C, T, G}), and α7α8 equal to X. Observables related by the extended Chargaff symmetry (3) appear on top of each other (separated by an horizontal line). The frequency of each observable P(X) was computed using Eq. (1) and the normalized version (cross correlation) using Eq. (5).

Figure 1

Symmetrically related cross-correlations in Homo Sapiens - Chromosome 1. The normalized cross-correlations as a function of the scale , together with that of its symmetrical related companions. Symmetries are significant also at scales where non trivial correlations are present z ≠ 1.

In order to understand the observations reported above it is necessary to formalize the symmetries that arise as composition of basic transformations. Starting from a reference observable , these symmetries are defined as compositions of the following two transformations:

(R) reverses the order in the pair: ;

(C) applies our extended symmetry equation (3) to the first of the two observable in the pair: .

Note that RC ≠ CR (i.e. R and C do not commute), RR = CC = Id (i.e. R, C are involutions), and CRC is the symmetry equivalent to equation (3). A symmetry S is defined by a set of different compositions of C and R. We denote by the set of observables obtained applying to observable Y all combinations of transformations in S. For example if S1 = {CRC} then ; if S2 = {CRC, R} then in addition to the set obtained from CRC we should add the observables obtained by applying R to every element of , that are and thus obtaining . The four symmetries we consider here are shown in Fig. 2 and correspond to: S1 (blue, obtained from {CRC} and corresponding to the extended Chargaff (4)), S2 (green, obtained from {CRC, R}), S3 (red, obtained from {R, C}), and S4 (black, obtained from {RCR, C}). We can now come back to Fig. 1 and interpret the observations as follows: at scales curves are significantly different from z = 1 and appear in pairs (same symbol, symmetry S1) which almost coincide even in the seemingly random fluctuations; around two pairs merge forming two groups of four curves each (symmetry S2). At larger scales all curves coincide (symmetry S3) at z ≈ 1 (no structure). At very large scales two groups of four observables separate (symmetry S4). Similar results are obtained for all choice of dinucleotides and for all chromosomes (see SI: Supplementary data)[35]. These results suggest that: (i) the extended Chargaff symmetry we conjectured in Eq. (4) is valid up to a critical scale ; (ii) there are other characteristic scales connected to the other symmetries.

Figure 2

Sets of symmetrically related observables. Starting from a reference observable , the different colors illustrate the nested sets related by symmetries S1–S4. The symmetries shown in the figure correspond to: S1 (blue, obtained from {CRC}), S2 (green, obtained from {CRC, R}), S3 (red, obtained from {R, C}), and S4 (black, obtained from {RCR, C}).

The scale-dependent results discussed above motivate us to quantify the strength of validity of symmetries at different scale . This is done computing for each symmetry S an indicator that measures the distance between the curves of symmetry-related pairs (X, X) and compares this distance to the ones that are not related by S. More precisely, for a given reference pair Y  = (X, X) and symmetry S ∈ {S1, S2, S3, S4}, we consider the following distance of Y to the set of observables obtained from symmetry S:where denotes the standard deviation of over all Y. We then average over the set of all Yref (all possible pairs X, X) to obtain a measure of the strength of symmetry S at scale given byNote that indicates full validity of the symmetry S at the scale (z is the same for all Y in ) and indicate that S is not valid at scale (z varies in as much as it varies in the full set).

Figure 3 shows the results for chromosome 1 and confirms the existence of a hierarchy of symmetries at different structural scales. The estimated relevant scales in chromosome 1 (of total length N ≈ 2 × 108) are L ≈ 102, L ≈ 103, and L ≈ 106. Note that L and L are compatible with the known average-size of transposable elements and isochores respectively[36,37]. Moreover, the results for all Homo-Sapiens chromosomes, summarised in Fig. 4, show that not only the hierarchy is present, but also that the scales L, L, and L are comparable across chromosomes. This remarkable similarity (see also[38-40]) suggests that some of the mechanisms shaping simultaneously structure and symmetry work similarly in every chromosomes and/or act across them (e.g. chromosome rearrangements mediated by transposable elements). This, and the scales of L, L, and L, provide a hint on the origin of our observations, which we explore below through the proposal of a minimal model.

Figure 3

Hierarchy of symmetries in Homo Sapiens - Chromosome 1. [Upper panel] The symmetry index as a function of the scale , the smaller the value the larger the importance of the symmetry. [Bottom panel] The color bars helps visualise the onset of the different symmetries: symmetry is considered present if 0 ≤ I ≤ 0.025 and bar is (linearly interpolated) from full color to white, correspondingly.

Figure 4

Hierarchy of symmetries in Homo Sapiens - All chromosomes. The symmetry index as a function of the scale for the full set of chromosomes in Homo Sapiens. The brighter the color, the larger is the relevance of the symmetry S1, S2, S3, or S4 (more precisely, if I is the minimum I in each chromosome, I ≤ 1.05I is set to full color, I ≥ 6.5I is set to white, with intermediate values interpolated between these extremes).

A minimal model

We construct a minimal domain model for DNA sequences s that aims to explain the observations reported above. The key ingredient of our model is the reverse-complement symmetry of domain-types, suggested by the fact that transposable elements act on both strands. Mobile elements are recognised to play a central role in shaping domains and other structures up to the scale of a full chromosome, as well as being considered responsible for the appearance of Chargaff symmetry[27]. Our model accounts for structures (e.g., the patchiness and long-range correlations in DNA) in a similar way as other domain models do, the novelty is that it shows the consequences to the symmetries of the full DNA sequence.

Motivated by our finding of the three scales L, L, and L, our model contains three key ingredients at different length scales (see Fig. 5 for an illustration):

Figure 5

Structure and symmetry at different scales: domain model. Structure and symmetry at different scales of genetic sequences can be explained using a simple domain model. Our model considers that the full sequence is composed of macro-structures (of size L) made by the concatenation of domains (of average size L < L), which are themselves correlated with neighbouring domains (up to a scale L < L < L). The biological processes that shapes domains imposes that, in each macrostructure, the types of domains comes in symmetric pairs. As a consequence, we show that four different symmetries S1–S4 are relevant at different scales (see text for details).

at small scales, a genetic sequence (of average size 〈n〉 ≈ L) is generated as a realization of a given process p. We do not impose a priori restrictions or symmetries on this process. We consider that one realization of this process builds a domain of type p. For a given domain type, the symmetrically related type is defined by the process as follows: take a realization (α1α2 … α) of the process p, revert its order (αα … α1), and complement each base , where , , , ;

at intermediate scales, a macro-structure is composed as a concatenation of domains d1d2 … d (of average size 〈m〉 ≈ L), each domain belonging to one of a few types. We assume that symmetrical related domains (generated by p and ) appear with the same relative abundance and size-distribution in a given macro-structure. The concatenation process is done so that domains of the same type tend to form clusters of average size L such that ;

at large scales , the full genetic sequence is composed by concatenations of macro-structures, each of them governed by different processes and statistics (e.g. different CG content)[10,11].

Statistical properties of the model and predictions

We now show how the model proposed above accounts for our empirical observation of a nested hierarchy of four symmetries S1-S4 at different scales. We start generating a synthetic sequence for a particular choice of parameters of the model described above (see section Methods for details). Figure 6 shows that such synthetic sequence reproduces the same hierarchy of symmetries we detected in Homo Sapiens.

Figure 6

Hierarchy of symmetries in a synthetic sequence generated by the domain model. The analysis of a synthetic genetic sequence generated by our model reproduces the hierarchy of symmetries observed in the human genome (compare the two panels to Figs 1 and 3). The synthetic sequence is obtained following steps (1)–(3) of the main text. As main stochastic processes p we use Markov chains with invariant probabilities μ such that μ(A) ≠ μ(T) and μ(C) ≠ μ(G) (no symmetries) (see Materials for details).

We now argue analytically why these results are expected. The key idea is to note that for different separations (between the two observables X and X) different scales of the model above dominate the counts used to compute through Eq. (1) (see Supplementary Information for a more rigorous derivation):Note that .

: is dominated by X and X in the same domain. As domain-types appear symmetrically in each macro-structure, . This is compatible with the conjecture (4).

: is dominated by X and X in different domains. As domains are independent realizations, the order of X and X becomes irrelevant and therefore R becomes a relevant symmetry (in addition to CRC). If domains of the same type tend to cluster, then for the main contribution to comes from X and X in different domains of the same type (i.e., on different realizations of the same process p).Note that .

: is dominated by X and X in different domains inside the same macro-structure. For the domains of X and X of different types can be considered independent form each other. Therefore, in addition to the previous symmetries, C is valid.Note that .

: is dominated by X and X in different macro-structures. Note that the frequency of X in one macro-structure and in a different macro-structure are, in general, different. Therefore, for generic X, X we have , meaning that S1 (and thus S2 and S3) is no longer valid. On the other hand, our conjectured Chargaff symmetry, Eq. (4), is valid for both X and X separately (because they are small scale observables). Therefore X and X can be interchanged in the composite observable Y.

Discussion

The complement symmetry in double-strand genetic sequences, known as the First Chargaff Parity Rule, is nowadays a trivial consequence of the double-helix assembly of DNA. However, from a historical point of view, the symmetry was one of the key ingredients leading to the double-helix solution of the complicated genetic structure puzzle, demonstrating the fruitfulness of a unified study of symmetry and structure in genetic sequences. In a similar fashion, here we show empirical evidence for the existence of new symmetries in the DNA (Figs 1–4) and we explain these observations using a simple domain model whose key features are dictated by the role of transposable elements in shaping DNA. In view of our model, our empirical results can be interpreted as a consequence of the action of transposable elements that generate a skeleton of symmetric domains in DNA sequences. Since domain models are known to explain also much of the structure observed in genetic sequences, our results show that structural complex organisation of single-strand genetic sequences and their nested hierarchy of symmetries are manifestations of the same biological processes. We expect that future unified investigations of these two features will shed light into their (up to now not completely clarified) evolutionary and functional role. For this aim, it is crucial to extend the analyses presented here to organisms of different complexity[21]. In parallel, we speculate that the unraveled hierarchy of symmetry at different scales could play a role in understanding how chromatin is spatially organised, related to the puzzling functional role of long-range correlations[41,42].

Methods

Algorithm used to generate the synthetic sequence

We create synthetic genetic sequences through the following implementation of the three steps of the model we proposed above:where columns (and rows) corresponds to the following order: [A, C, G, T].

The processes p we use to generate genetic sequences are Markov processes of order one such that the nucleotide s at position i is drawn from a probability , where M is a 4 by 4 stochastic matrix. The matrices M are chosen such that the processes’ invariant measures μ do not satisfy the Chargaff property: μ(A) ≠ μ(T) and μ(C) ≠ μ(G). The exponential decay of correlations of the Markov chains determines the domain sizes L (in our case ).

We use the processes p to generate chunks of average size 150 (the length of each chunck was drawn uniformly in the range [130, 170]. With probability 1/2, we applied the reverse-complement (CRC) operation to the chunck before concatenating it to the previous chunck (process ). This choice implies that the typical cluster size is . The process of concatenating chunks together is repeated to form a macrostructure of length .

We concatenate two different macrostructures, obtained from steps (1) and (2) with two different matrices M and M:

Data handling

Genetic sequences of Homo Sapiens were downloaded from the National Center for Biotechnology Information (ftp://ftp.ncbi.nih.gov/genomes/H_sapiens). We used reference assembly build 38.2. The sequences were processed to remove all letters different from A, C, G, T (they account for ≈1.66% of the full genome and thus their removal has no significant impact on our results).

Codes

Reference[35] contains data and codes that reproduce the figures of the manuscript for different choices of observables and chromosomes.