Nucleotide sequences x=x[n] are sequences of elements (xi,i∈[n]) with values from the alphabet A:={A,C,G,T} . Here [n]:=[1,n]={1,…,n} is an interval of (positive) integers.

If x=(x1,…,xn) is the leading strand of a DNA sequence, then the complementary one (the lagging strand read in the opposite direction) is denoted by x∗=(x1∗,…,xn∗) , where xi∗ is the complementary nucleotide to xi in the ith base pair, and x∗=x∗[n]:=(x∗1,…,x∗n)=(xn∗,…,x1∗) . This determines the complementary transformation. For instance, (1) x=(x1,…,xn):…CGGATTTAGCTA…→, (2) x∗=(x1∗,…,xn∗):…GCCTAAATCGAT…←, (3) x∗=(xn∗,…,x1∗):…TAGCTAAATCCG…→. Chargaff and his colleagues (Rudner et al., 1968) have noticed that |{t∈[n]:xt=A}|≈|{t∈[n]:xt=T}|,|{t∈[n]:xt=C}|≈|{t∈[n]:xt=G}|, ( |A| is the number of elements in a set A) which actually means that |{t∈[n]:xt=ν}|≈|{t∈[n]:xt∗=ν}|,∀ν∈A. Thus, if x is treated as a random sequence, the last expression can be interpreted and generalized as follows: a probabilistic law generating x is invariant with respect to the complementary transformation x→x∗ .

Let us start with basic notation and notions. Set N=[n] , fix some positive integer m<n/2 and define the m-interior Nm∘ , the m-boundary ∂Nm of N , and a collection of neighbourhoods: (4) N∘=Nm∘:=[m+1,…,n−m],∂N=∂Nm:=N∖Nm∘,N(ℓ)=Nm(ℓ):=[ℓ−m,ℓ+m]∖{ℓ},ℓ∈N∘. Given x∈An and a set of indices J⊂N , let xJ:=(xi,i∈J) denote the corresponding subsequence of x treated as an element of A|J| .

Suppose that values of m-MRF x are fixed on the boundary ∂N : x∂N=b a.s. for some b∈A2m . Denote Xb:={w∈An:w∂N=b}. From Hammersley–Clifford theorem (Besag, 1974), we obtain the following statement.

The statement is well-known, it is just rewritten in the notation introduced above.

Let us recall that DNA strand symmetry means that probability distribution of oligonucleotides (sequences of adjacent nucleotides) of the both complementary strands of DNA, read in the respective direction, are similar in some sense. Having in mind the definition of m-MRF, the following formal definition of DNA strand symmetry can be given in terms of complementary transformation w→w∗,w∈A2m+1, defined in Section 2.1.

Thus, for locally symmetric sequence, the marginal conditional distributions given m nearest neighbours (from the each side) are invariant under the complementary transformation. Under the assumption that DNA sequence x is m-MRF, the local strand symmetry can be expressed in terms of the conditional distributions p(v|u) and/or the generalized logits Λv(u) .

For characterization of local symmetry in terms of the generalized logits, it is convenient to change the initial alphabet A={A,C,G,T} of nucleotides v to A12=A1×A1 , A1:={−1,+1} , via mapping v→z:=(s,y) by the rule indicated in Table 1. The components s=s(v)∈A1 and y=y(v)∈A1 of a nucleotide v∈A represent its bonding property strong versus weak and its hydrophobic property pyrimidine (large molecule, less hydrophobic) versus purine (small molecule, more hydrophobic), respectively.

Now, let x=(x1,…,xk)∈Ak be a nucleotide sequence in the leading strand of DNA and let x∗ be its complement read from the left to the right but taken in the common direction. Set (10) z=z(x):=(s→,y→)∈A12k,(s→)i:=s(xi),(y→)i:=y(xi),i=1,…,k, (11) z∗=z∗(x):=(s→∗,y→∗)=(s←,−y←), (12) (s←)i=(s→)k−i+1,(y←)i=(y→)k−i+1,i=1,…,k. Then z(x∗)=z∗(x) . To illustrate the notation we apply them to the nucleotide sequence from (1) (to save space here and below we will omit the numeral 1): x=(x1,…,xn):…CGGATTTAGCTA…,s=(s1,…,sn):…+++—–++–…,y=(y1,…,yn):…+—+++–++-…,x∗=(xn∗,…,x1∗):…TAGCTAAATCCG…,s→∗=(sn,…,s1):…–++—–+++…,y→∗=−(yn,…,y1):…+–++—+++-…. In what follows we identify p(v∣u) with p(z(v)∣z(u)) and Λv|r(u),r=A, with Λz(v)(z(u)) .

Let us introduce functions that are symmetric (antisymmetric) with respect to the complementary transformation z→z∗,z∈A12k, defined in (10)–(12).

When estimating the generalized logits Λτ(w) one needs some parametrization. Below convenient parametric representations for symmetric and antisymmetric functions are presented.

According to the recoding rule defined in Table 1 and (10)–(12), z=(s,y),s,y∈A12m, and hence in the sequel we deal with functions ψ(s,y) , ψ:A1k×A1k→R , k:=2m .

Let J⊂K:={1,…,k} . Define the conjugate set J∗ of the set J by (26) J∗:=k+1−J={k+1−j:j∈J}. For a given sequence s=(s1,…,sk)∈A1k , denote (27) sJ:=∏i∈Jsi,s∅:=1. Any function ψ:A1k×A1k→R has the unique representation (28) ψ(s,y)=∑J′,J⊂KaJ′JsJ′yJ,s,y∈A1k, where summation is over all subsets of K (including the empty set ∅), aJ′J=aJ′J(ψ),J′,J⊂K, are free parameters determining the function ψ. In general, there are 4k free parameters.

For a symmetric (antisymmetric) function ψ, we have (29) ψ(s←,−y←)=ψ(s,y)(respectively,ψ(s←,−y←)=−ψ(s,y)). Consequently, in the case of the symmetric ψ, for all s,y∈A1k , (30) ∑J′,J⊂KaJ′JsJ′yJ=∑J′,J⊂K(−1)|J|aJ′JsJ∗′yJ∗=∑J′,J⊂K(−1)|J|aJ∗′J∗sJ′yJ, and hence (31) aJ′J=(−1)|J|aJ∗′J∗,J′,J⊂K. If J∗=J and J∗′=J′ (i.e. the both subsets J′ and J are self-conjugate), the set J has an even number of elements and the equations (31) become the identities. Thus, there are no restrictions on the parameter aJ′J values in this case. Let k∗=k∗(k) denote the total number of the self-conjugate subsets of K.

Let τ be some total order (enumeration of elements) in the class of pairs (J′,J) of the set K. Equations (31) imply that, for not self-conjugate pairs (J′,J) , (J′,J)≠(J∗′,J∗) , values of the coefficients aJ′,J,τ(J′,J)<τ(J∗′,J∗), uniquely determine values of the remaining coefficients aJ′,J,τ(J′,J)>τ(J∗′,J∗) \tau ({J^{\prime }_{\ast }},{J_{\ast }})$]]>. Define (32) K2:={(J′,J):τ(J′,J)<τ(J∗′,J∗)}, (33) K20:={(J′,J):τ(J′,J)=τ(J∗′,J∗)}. From (28), (31), (32) and (33) we derive a general parametric form of a symmetric function ψ: (34) ψS(s,y)=∑(J′,J)∈K20aJ′JsJ′yJ+∑(J′,J)∈K2aJ′J(sJ′yJ+(−1)|J|sJ∗′yJ∗). It has (35) kS=kS(k):=k∗2+(4k−k∗2)/2 free parameters.

The case of antisymmetric function differs from that of symmetric function only in additional minus sign in equations (31). For self-conjugate pairs (J′,J) , these equations hold if and only if aJ′,J=0 . Thus, the first summand in (34) and in (35) disappears giving the function (36) ψA(s,y)=∑(J,J′)∈K2aJ,J′(sJyJ′−(−1)|J|sJ∗′yJ∗) with (37) kA=kA(k):=(4k−k∗2)/2 free parameters.

For m=1 , k∗=2 , thus kA=(42−22)/2=6 and kS=k∗2+kA=10 . Then symmetric (34) and antisymmetric (36) functions in a general form are given by ψS(s,y)=a∅∅+a∅{12}y1y2+a{12}∅s1s2+a{12}{12}s1s2y1y2+a∅{1}(y1+y2)+a{1}∅(s1−s2)+a{1}{1}(s1y1−s2y2)+a{1}{2}(s1y2−s2y1)+a{1}{12}(s1y1y2−s2y1y2)+a{12}{1}(s1s2y1+s1s2y2),ψA(s,y)=a∅{1}(y1−y2)+a{1}∅(s1+s2)+a{1}{1}(s1y1+s2y2)+a{1}{2}(s1y2+s2y1)+a{1}{12}(s1y1y2+s2y1y2)+a{12}{1}(s1s2y1−s1s2y2), respectively.

Let us introduce the following data structure of the observed sequence x∈An with n=(nm+1)·(m+1)−1 , the quantity nm being an integer: (38) D:={(vℓ,zℓ),ℓ∈S},S=Sn,m={1,2,…,nm}, where vℓ:=x(m+1)ℓ is a response variable and zℓ=xUm((m+1)ℓ)∈A2m is a vector of explanatory variables, ℓ∈S .

Assumption (Am):

Assumption (Am) ensures that usual conditions of the generalized logit model with the response variable v∈A and the vector z∈A2m of explanatory variables are satisfied, see Agresti (1990), Stokes et al. (2001). Remark 2

Suppose that a DNA sequence x is an outcome of a “long” homogeneous Markov evolution and hence has a stationary distribution. Assumption (Am) imposed on x is compatible with some common DNA evolutionary models. In particular, assumption (Am) with m=2 hold for the independent codon evolution (Goldman and Yang, 1994). Assumption (Am) is also fulfilled if x is generated by m-MRF. Thus, it is valid in case of time-reversible, site-homogeneous and context-dependent Markov evolution model with m-order nearest neighbour interactions (see, e.g. Jensen, 2005). However, it is satisfied for some non-homogeneous, say (m+1) -periodic, MRF of order m as well.

In general, the introduced regression-type data structure supplemented with a saturated generalized logit model for m-order conditional probabilities does not determine the distribution of x. However, if assumption (Am) holds for S=N∘ (to be precise, for all shifts ((m+1)S+ℓ)∩N∘ of the set of central nucleotides (m+1)S by ℓ, ℓ=1,…,m−1 , simultaneously), then, due to Hammersley–Clifford theorem (Proposition 1), x is m-MRF, and m-order generalized logits take the form of (8) and determine the distribution of x.

We analyse data of bacterial genomes (1221 genomes) taken from the database GenBank (https://www.ncbi.nlm.nih.gov/genbank/). In order to bypass the data sparsity problem the longest non-coding (for the both strands) DNA sequences are extracted from each genome. Assuming that the extracted sequences satisfy assumption (Am) with m=2 we test their first order local symmetry.

In Fig. 1, the length distribution density of the extracted sequences is plotted in a logarithmic scale. The sequence lengths range from 1891 to 42901 with median 6605 and mean 7721. About a half of the sequences have length between 6000 and 8000. Since we assume (Am) with m=2 , the logit analysis is based on three-dimensional contingency tables (64 cells) of nonintersecting triplets in the DNA sequences. The average and median of cell counts in the tables are 40 and 34, respectively. The percentage of cells with less than 6 counts does not exceed 1%. Thus we can ignore p-value approximation problems incident to statistical analysis of sparse contingency tables (Agresti, 1990).

Generalized logit model is fitted to the data and the Wald criterion is applied to test if the coefficients of the generalized logit model satisfy conditions implied by antisymmetric (36) and symmetric (34) components of generalized logits specified in Proposition 2.

In Figs. 2–4, values of the logarithmized Student statistic (the Student statistic S transformed by S→sgn(S)log2(1+|S|) ) for testing the significance of the coefficients of response functions ψ−,ψ+ and ψ defined in (13)–(15), respectively, are presented. For better visibility of the logarithmized Student statistic distributions, we use the violin plot (Hintze and Nelson, 1998; Wickham, 2016), which combines a box plot and a kernel density plot that is rotated and placed on each side, to show the distribution shape of the data. The first 6 coefficients represent the antisymmetric part of the response functions and the last 10 represent the symmetric part. According to Proposition 2, in case of the local symmetry, the first 2 response functions should be antisymmetric while the last one should be symmetric. Hence the last 10 and, respectively, the first 6 coefficients should be insignificant. In the figures, the approximate critical value obtained by 3σ rule (i.e. for the significance level ≈0.0054) corresponds to y-coordinates ±2.

First response function (expected to be antisymmetric). The distributions of its coefficient estimates are represented in Fig. 2. The coefficient estimates of the antisymmetric part (white violins) have skewed distributions, especially the second, which is left-skewed and has large positive bias, and the third, which is right-skewed and has large negative bias. The distributions in the symmetric part (grey violins) are quite symmetric about zero. A large proportion of the non-coding DNA sequences (>40%) has significant (at the approximate significance level of 0.005) 7th coefficient (7th parameter) expected to be zero in the case of local symmetry.

In what follows only violations of local symmetry (grey violins) are discussed.

Second response function (expected to be antisymmetric). A major part (>70%) of the non-coding sequences has significant 7th coefficient (23rd parameter) expected to be zero in the case of local symmetry. A large proportion of the sequences also exhibits significant deviations from 0 of the 8th coefficient (24th parameter).

Third response function (expected to be symmetric). The second coefficient (34th parameter) expected to be zero in case of local symmetry shows a clear tendency to deviate significantly from 0.

In Fig. 5, centres of 8 clusters obtained using the standard R function for k-means clustering (, 2018) of 48-dimensional vectors of the estimated model parameters (i.e. estimated coefficients of the all three response functions) are drawn. The coordinates of each centre are joint thus representing 8 different patterns of their interrelationships. The centre of the 8th cluster represents DNA sequences which approximately satisfy the local symmetry hypothesis. The sequences of the third cluster are also rather close to symmetry. Clusters 8 and 3, however, apparently differ in the regions [17,19] and [19,41] . All the clusters are similar in [1,6] . In the grey zones (regions [7,16] and [23,38] ), we have two triplets of similar clusters: (1,2,6) and (4,5,7) . The 39th parameter for cluster 1 clearly differs from that of clusters 2 and 6 having the opposite sign. The same applies to clusters 4, 5 and 7, respectively. Clusters 2 and 6, as well as 5 and 7, exhibit some discrepancy in values of parameter 41. Cluster 5 also has specific values in the region [18,19] .

Note that the deviations of the parameter estimates in the grey region, i.e. their deviations from the DNA local symmetry hypothesis, are quite symmetric, see also Figs. 2–4. This observation is consistent with the ISFDP property noticed in Powdel et al. (2009).

Elements of DNA sequences x are treated as random variables taking values from the alphabet A:={A,C,G,T} . A definition of the local symmetry of x of order m is given and is characterized in terms of generalized logits (Židanavičiūtė, 2010). To test the first order local symmetry of non-coding sequences of bacteria genoms a special regression-type structure is imposed on probability distribution of x (assumption (Am) with m=2 ). It defines a generalized logit model with 48 scalar parameters. In the case of the first order local symmetry, 22 of them should vanish.

The generalized logit model was fitted to the longest non-coding sequences of 1221 bacteria genomes taken from GenBank and Wald test was applied to check the null hypothesis of the first order local symmetry.

Conclusions: 1.

Most of the non-coding sequences of bacteria genomes do not possess the first order local symmetry.