Due to symmetry, the nature is perfect.
Spices of asymmetry make it beautiful.
1 Introduction
Genetically (or biologically) informative sequences can be defined as those which are either close to a known genetically important sequence or are far from sequences known to be noninformative. The first criterion seems to be more practical, however it is limited since it tries to reproduce what is already known. The second principle is more fundamental and more convenient for mathematical formalization and statistical inference. When employing this principle, the problem is how to define the noninformative genetic sequence (we call it the genetic noise), i.e. the sequence which has no genetically or biologically important information.
A model of the genetic noise is also crucial for statistical hypotheses testing, the phylogenetic tree reconstruction, simulations of the (neutral) evolutions, and in assessing the variability and uncertainty.
Genome regions whose evolution is not subjected to natural selection pressure and hence evolve with a neutral mutation rate can be viewed as the genetic noise. Those regions could be parts of non-coding regions of genoms of primitive species.
A generic formulation of empirical findings is sometimes called a stylized fact. The definition of the genetic noise should be consistent with the stylized facts about non-coding DNA sequences as well as with a probabilistic model of their evolution. Thus, the general aim of our investigation is to specify and to test statistically the basic properties of non-coding DNA sequences implied by a model of DNA evolution (Markov property, homogeneity, long-range dependence, reverse-complement symmetry, CpG content, etc.). In this work we focus on symmetry/asymmetry properties of two complementary DNA strands.
Chargaff’s second parity rule. The simplest hypothesis of DNA strand symmetry (sometimes referred to as Chargaff’s second parity rule) states that proportions of nucleotides of the same base pair are approximately equal within single DNA strands (Rudner et al., 1968), i.e. %A ≈ %T and %C ≈ %G. Since the lagging strand is read in the reverse order, an extension of this first-order symmetry to higher-orders is called reverse-complement symmetry, or intra-strand parity (ISP) (Powdel et al., 2009), or simply strand symmetry (Baisnée et al., 2002; Zhang and Huang, 2008). Although rather natural, this universal phenomenon of strand symmetry in the chromosomes needs explicit description and explanation. Actually, it may be the effect of a wide range of mechanisms operating at multiple orders and length scales (Baisnée et al., 2002).
Thus far the issue about strand symmetry, its origins and biological significance is controversial. On the one hand, results of empirical studies using various asymmetry measures and visualization tools show that for long DNA sequences (approximate) strand symmetry generally holds with rather rare exceptions. The fact that the strand symmetry should hold at the equilibrium state is also derived theoretically (Sueoka, 1995; Lobry, 1995). Baisnée et al. (2002) defined strand symmetry indices through relative ${L_{1}}$ distance between the observed frequencies of respective reverse-complementary oligonucleotides and compare them with critical values calculated for completely random sequences. In Kong et al. (2009), various symmetry indices (reverse, complement and inverse symmetry indices, global as well as segmental) based on ${L_{2}}$ distance have been calculated for 786 complete chromosomes. The authors have found that reverse-complement symmetry (inverse-complement plus reverse-symmetry in terms of the authors) is prevalent in complex patterns in most chromosomes. Rosandić et al. (2016) considered 20 symbolic quadruplets of trinucleotides obtained via interstrand mirror symmetry mappings (direct, reverse complement, complement, and reverse) and demonstrated quadruplet’s symmetries in chromosomes of wide range of organisms, from Escherichia coli to human genomes. Powdel et al. (2009) have noticed another strand symmetry manifestation, intra-strand frequency distribution parity (ISFDP), which represents closeness of frequency distributions between the complementary mono/oligonucleotides. This general feature (with rare exceptions) was observed in chromosomes of bacteria, archaea and eukaryotes. It has been also noticed that the frequency of an genomic word is more similar to the frequency of its reversed complement than to the frequencies of other words of equivalent composition. This phenomenon is called exceptional symmetry. Afreixo et al. (2017) proposed a new measure to evaluate the exceptional symmetry effect based on discrepancy between frequency of symmetric word pair and frequencies of word pairs of equivalent composition. They identified words that show high symmetry effect across the 31 species, and across the 9 animal species studied. Fractal-like symmetry structures are considered in Petoukhov et al. (2018). Sobottka and Hart (2011) proposed a model based on a hidden Markov process for approximating the distributions of primitive DNA sequences. The model provides an alternative interpretation of strand symmetry and describes new symmetries in bacterial genomes. Cristadoro et al. (2018) introduced flexible statistical measures of symmetry and used them to define an extended Chargaff symmetry. The definition actually coincides with global strand symmetry of genoms defined and studied in Simons et al. (2005). Domain models introduced in Cristadoro et al. (2018) alow to explain simultaneously symmetries as well as non-random structures in genetic sequences and unravel previously unknown symmetries, which are organized hierarchically through different scales.
On the other hand, statistical analyzes of the genomic sequences (Shporer et al., 2016; Tavares et al., 2018), especially those based on Markov-type models (Hart and Martínez, 2011; Hart et al., 2012), have demonstrated significant deviations from the second Chargaff’s parity rule and its extensions. A statistical IS-Poisson model introduced in Shporer et al. (2016) assumes that frequencies of oligonucleotides (DNA k-mers) follow the Poisson distribution. The model allows to conclude that for k-mers with low k (even for nucleotides, $k=1$) violations of symmetry, although extremely small, are significant. In Tavares et al. (2018), both the distance distributions and the frequencies of symmetric words in the human DNA have been compared. The results obtained suggest that some asymmetries in the human genome go far beyond Chargaff’s rules.
One of the explanations of strand asymmetry (skew), i.e. violation of symmetry, is mutation bias. When investigating asymmetries in mutation patterns, phylogenetic estimation based on maximum likelihood can be applied. Usually independent evolution models completely determined by nucleotide substitution rates are employed, see, e.g. Faith and Pollock (2003), Marin and Xia (2008). Note that mathematical models for evolutionary inference considered in Parks (2015) also assume independent evolution. However, Siepel and Haussler (2004) presented extensions of standard phylogenetic models with context-dependent substitution and showed that the new models improve goodness of fit substantially for both coding and non-coding data. Moreover, considering context dependence leads to much larger improvements than does using a richer substitution model or allowing for rate variation across sites, under the assumption of site independence. We refer to Bérard and Guéguen (2012) for a more recent application of context-dependent substitution models in a phylogenetic context.
In this paper, DNA strand local symmetry introduced in Židanavičiūtė (2010) is tested for the longest non-coding (in the both leading and lagging strands) sequences of bacterial genomes taken from GenBank (https://www.ncbi.nlm.nih.gov/genbank/). Validity of a special regression-type probabilistic structure of the data is supposed. This structure is compatible with probability distribution of random nucleotide sequences at a steady state of a context-dependent reversible Markov evolutionary process (Jensen, 2005), see also Arndt et al. (2003), Lunter and Hein (2004). The null hypothesis of strand local symmetry is rejected in majority of bacterial genomes suggesting that even neutral mutations are skewed with respect to leading and lagging strands.
The rest of the paper is organized as follows. In the next section the definition of strand local symmetry is presented and characterization of this property in terms of generalized logits is given. Results of statistical analysis are discussed in Section 3. We end with some concluding remarks.
2 Local Symmetry
In this section we present the formal definition of local symmetry (Židanavičiūtė, 2010) and recall necessary notions and facts about discrete Markov random fields and loglinear modelling.
2.1 Complementary Transformation
Nucleotide sequences $x={x_{[n]}}$ are sequences of elements $({x_{i}},i\in [n])$ with values from the alphabet $\mathcal{A}:=\{\texttt{A},\texttt{C},\texttt{G},\texttt{T}\}$. Here $[n]:=[1,n]=\{1,\dots ,n\}$ is an interval of (positive) integers.
If $x=({x_{1}},\dots ,{x_{n}})$ is the leading strand of a DNA sequence, then the complementary one (the lagging strand read in the opposite direction) is denoted by ${x^{\ast }}=({x_{1}^{\ast }},\dots ,{x_{n}^{\ast }})$, where ${x_{i}^{\ast }}$ is the complementary nucleotide to ${x_{i}}$ in the ith base pair, and ${x_{\ast }}={x_{\ast [n]}}:=({x_{\ast 1}},\dots ,{x_{\ast n}})=({x_{n}^{\ast }},\dots ,{x_{1}^{\ast }})$. This determines the complementary transformation. For instance,
Chargaff and his colleagues (Rudner et al., 1968) have noticed that
(1)
\[\begin{aligned}{}x& =({x_{1}},\dots ,{x_{n}}):\hspace{1em}\overrightarrow{\dots \texttt{CGGATTTAGCTA}\dots }\hspace{0.1667em},\end{aligned}\]
\[\begin{aligned}{}& \big|\big\{t\in [n]:{x_{t}}=\texttt{A}\big\}\big|\approx \big|\big\{t\in [n]:{x_{t}}=\texttt{T}\big\}\big|,\\ {} & \big|\big\{t\in [n]:{x_{t}}=\texttt{C}\big\}\big|\approx \big|\big\{t\in [n]:{x_{t}}=\texttt{G}\big\}\big|,\end{aligned}\]
($|A|$ is the number of elements in a set A) which actually means that
\[ \big|\big\{t\in [n]:{x_{t}}=\nu \big\}\big|\approx \big|\big\{t\in [n]:{x_{t}^{\ast }}=\nu \big\}\big|,\hspace{1em}\forall \nu \in \mathcal{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\to {x_{\ast }}$.2.2 Basics of Markov Random Fields
Let us start with basic notation and notions. Set $\mathcal{N}=[n]$, fix some positive integer $m<n/2$ and define the m-interior ${\mathcal{N}_{m}^{\circ }}$, the m-boundary $\partial {\mathcal{N}_{m}}$ of $\mathcal{N}$, and a collection of neighbourhoods:
Given $x\in {\mathcal{A}^{n}}$ and a set of indices $J\subset \mathcal{N}$, let ${x_{J}}:=({x_{i}},i\in J)$ denote the corresponding subsequence of x treated as an element of ${\mathcal{A}^{|J|}}$.
(4)
\[\begin{aligned}{}& {\mathcal{N}^{\circ }}={\mathcal{N}_{m}^{\circ }}:=[m+1,\dots ,n-m],\hspace{2em}\partial \mathcal{N}=\partial {\mathcal{N}_{m}}:=\mathcal{N}\setminus {\mathcal{N}_{m}^{\circ }},\\ {} & \mathcal{N}(\ell )={\mathcal{N}_{m}}(\ell ):=[\ell -m,\ell +m]\setminus \{\ell \},\hspace{1em}\ell \in {\mathcal{N}^{\circ }}.\end{aligned}\]Definition 1.
A random sequence $x\in {\mathcal{A}^{n}}$ is called an m-order Markov random field (MRF) with the state space $\mathcal{A}$ and the collection of neighbourhoods (4) iff $\forall a\in {\mathcal{A}^{n}}$ and $\forall \ell \in {\mathcal{N}_{m}^{\circ }}$
A MRF x is called an m-order homogeneous MRF (m-MRF) if its m-order marginal conditional probabilities given in the right-hand side of (5) are independent of the site $\ell \in {N^{\circ }}$.
(5)
\[ \mathbf{P}\big\{{x_{\ell }}={a_{\ell }}\mid {x_{i}}={a_{i}},\hspace{0.1667em}i\in [n],\hspace{0.1667em}j\ne \ell \big\}=\mathbf{P}\big\{{x_{\ell }}={a_{\ell }}\mid {x_{{\mathcal{N}_{m}}(\ell )}}={a_{{\mathcal{N}_{m}}(\ell )}}\big\}.\]Definition 2.
For a fixed reference value $r\in \mathcal{A}$ and given m-order marginal conditional probabilities
the respective generalized logit ${\Lambda _{v}}(u)={\Lambda _{v|r}}(u)$ of a state $v\in \mathcal{A}$ versus r, given the neighbouring values $u\in {\mathcal{A}^{2m}}$, is defined as
where we set $\log (0/0)=0$ and $\log (p/0)=\infty $ for $p>0$.
Suppose that values of m-MRF x are fixed on the boundary $\partial \mathcal{N}$: ${x_{\partial \mathcal{N}}}=b$ a.s. for some $b\in {\mathcal{A}^{2m}}$. Denote
\[ {\mathcal{X}_{b}}:=\big\{w\in {\mathcal{A}^{n}}:\hspace{2.5pt}{w_{\partial \mathcal{N}}}=b\big\}.\]
From Hammersley–Clifford theorem (Besag, 1974), we obtain the following statement.Proposition 1.
Suppose the distribution of m-MRF x is positive on ${\mathcal{X}_{b}}$, i.e. $\mathbf{P}\{x=w\}>0$ for all $w\in {\mathcal{X}_{b}}$. Then the distribution of x is uniquely determined by the family of generalized logits ${\Lambda _{v|r}}(u)$, $r,v\in \mathcal{A}$, $u\in {\mathcal{A}^{2m}}$, which for $w\in {\mathcal{A}^{2m+1}}$, take the following form
and in general depend on $M=(|\mathcal{A}|-1)|\mathcal{A}{|^{m}}$ free scalar parameters. Here ${\lambda _{m}}$: ${\mathcal{A}^{m+1}}\to \mathbf{R}$ is an arbitrary function and ${w^{(r)}}$ is obtained from w by substituting r for ${w_{m+1}}$.
(8)
\[ {\Lambda _{{w_{m+1}}|r}}({w_{{\mathcal{N}_{m}}(m+1)}})={\sum \limits_{j=1}^{m+1}}\big[{\lambda _{m}}({w_{[j,m+j]}})-{\lambda _{m}}\big({w_{[j,m+j]}^{(r)}}\big)\big],\]The statement is well-known, it is just rewritten in the notation introduced above.
2.3 Local Symmetry: Definition and Characterization
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\to {w_{\ast }},w\in {\mathcal{A}^{2m+1}},$ defined in Section 2.1.
Definition 3 (See Židanavičiūtė, 2010).
A random sequence ${x_{[n]}}$ is m-order locally symmetric ($m<n/2$) iff
for all $\ell \in {\mathcal{N}^{\circ }},\hspace{2.5pt}v\in \mathcal{A},u\in {\mathcal{A}^{2m}}$.
(9)
\[ \mathbf{P}\{{x_{\ell }}=v\mid {x_{{\mathcal{N}_{m}}(\ell )}}=u\}=\mathbf{P}\big\{{x_{\ell }}={v^{\ast }}\hspace{0.1667em}\big|\hspace{0.1667em}{x_{{\mathcal{N}_{m}}(\ell )}}={u^{\ast }}\big\}\]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 ${\Lambda _{v}}(u)$.
Table 1
Nucleotide recoding rule.
| Purine | (Bonds) | Pyrimidine | s | |
| Weak (2 bonds) | A | $(\hspace{2.5pt}=\hspace{2.5pt})$ | T | $s=-1$ |
| Strong (3 bonds) | G | $(\hspace{2.5pt}\equiv \hspace{2.5pt})$ | C | $s=+1$ |
| y | $y=-1$ | $y=+1$ |
For characterization of local symmetry in terms of the generalized logits, it is convenient to change the initial alphabet $\mathcal{A}=\{\texttt{A},\texttt{C},\texttt{G},\texttt{T}\}$ of nucleotides v to ${\mathcal{A}_{1}^{2}}={\mathcal{A}_{1}}\times {\mathcal{A}_{1}}$, ${\mathcal{A}_{1}}:=\{-1,+1\}$, via mapping $v\to z:=(s,y)$ by the rule indicated in Table 1. The components $s=s(v)\in {\mathcal{A}_{1}}$ and $y=y(v)\in {\mathcal{A}_{1}}$ of a nucleotide $v\in \mathcal{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=({x_{1}},\dots ,{x_{k}})\in {\mathcal{A}^{k}}$ be a nucleotide sequence in the leading strand of DNA and let ${x^{\ast }}$ be its complement read from the left to the right but taken in the common direction. Set
Then $z({x_{\ast }})={z_{\ast }}(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):
(10)
\[\begin{aligned}{}& z=z(x):=(\overrightarrow{s},\overrightarrow{y})\in {\mathcal{A}_{1}^{2k}},\hspace{1em}\hspace{2.5pt}{(\overrightarrow{s})_{i}}:=s({x_{i}}),\hspace{1em}\hspace{2.5pt}\hspace{2.5pt}{(\overrightarrow{y})_{i}}:=y({x_{i}}),\hspace{1em}i=1,\dots ,k,\end{aligned}\]
\[\begin{aligned}{}& x=({x_{1}},\dots ,{x_{n}}):\hspace{1em}\dots \texttt{CGGATTTAGCTA}\dots \hspace{0.1667em},\\ {} & s=({s_{1}},\dots ,{s_{n}}):\hspace{1em}\dots \texttt{+++-----++--}\dots \hspace{0.1667em},\\ {} & y=({y_{1}},\dots ,{y_{n}}):\hspace{1em}\dots \texttt{+---+++--++-}\dots \hspace{0.1667em},\\ {} & {x_{\ast }}=({x_{n}^{\ast }},\dots ,{x_{1}^{\ast }}):\hspace{1em}\dots \texttt{TAGCTAAATCCG}\dots \hspace{0.1667em},\\ {} & {\overrightarrow{s}_{\ast }}=({s_{n}},\dots ,{s_{1}}):\hspace{1em}\dots \texttt{--++-----+++}\dots \hspace{0.1667em},\\ {} & {\overrightarrow{y}_{\ast }}=-({y_{n}},\dots ,{y_{1}}):\hspace{1em}\dots \texttt{+--++---+++-}\dots \hspace{0.1667em}.\end{aligned}\]
In what follows we identify $p(v\mid u)$ with $p(z(v)\mid z(u))$ and ${\Lambda _{v|r}}(u),\hspace{0.1667em}r=\texttt{A},$ with ${\Lambda _{z(v)}}(z(u))$.Let us introduce functions that are symmetric (antisymmetric) with respect to the complementary transformation $z\to {z_{\ast }},\hspace{0.1667em}z\in {\mathcal{A}_{1}^{2k}},$ defined in (10)–(12).
Definition 4.
A function $\psi :{\mathcal{A}_{1}^{2k}}\to \mathbf{R}$ is called symmetric (antisymmetric) with respect to the complementary transformation $w\to {w_{\ast }}$ iff $\psi (w)=\psi ({w_{\ast }})$ (respectively, $\psi (w)=-\psi ({w_{\ast }})$) for all $w\in {\mathcal{A}_{1}^{2k}}$.
Proposition 2.
Let $p(\eta \mid w),\eta \in {\mathcal{A}_{1}^{2}}$, $w\in {\mathcal{A}_{1}^{4m}}$, denote the m-order conditional probabilities of a bivariate random sequence $z(x)$ obtained from the nucleotide sequence $x\in {\mathcal{A}^{2m+1}}$ via z-transform (10). The following statements are equivalent:
-
(a) the sequence x and the marginal conditional probabilities $p(\cdot \mid \cdot )$ are m-order locally symmetric;
-
(b) there exist a symmetric function $\psi :{\mathcal{A}_{1}^{4m}}\to \mathbf{R}$ and two antisymmetric functions ${\psi _{-}}:{\mathcal{A}_{1}^{4m}}\to \mathbf{R}$ and ${\psi _{+}}:{\mathcal{A}_{1}^{4m}}\to \mathbf{R}$ such thatAnother form of (13)–(15) expressed in terms of the generalized logits ${\Lambda _{s,y}}(w)$:
(13)
\[\begin{aligned}{}& \log \bigg(\frac{p(-,+\mid w)}{p(-,-\mid w)}\bigg)={\psi _{-}}(w),\end{aligned}\]
Proof.
From the definition of generalized logits (7) and the recoding rule defined in Table 1 and (10), (11), we obtain, for all $w\in {\mathcal{A}_{1}^{4m}}$,
Consequently, for all $w\in {\mathcal{A}_{1}^{4m}}$,
Thus, ${\psi _{-}}(u)$ is antisymmetric. Analogously, for all $w\in {\mathcal{A}_{1}^{4m}}$,
and The proof is completed. □
\[\begin{aligned}{}{\Lambda _{-,+}}(w)& =\log \bigg(\frac{p(-,+\mid w)}{p(-,-\mid w)}\bigg)={\psi _{-}}(w),\\ {} {\Lambda _{+,+}}(w)-{\Lambda _{+,-}}(w)& =\log \bigg(\frac{p(+,+\mid w)}{p(+,-\mid w)}\bigg)={\psi _{+}}(w),\\ {} {\Lambda _{+,+}}(w)+{\Lambda _{+,-}}(w)-{\Lambda _{-,+}}(w)& =\log \bigg(\frac{p(+,+\mid w)\cdot p(+,-\mid w)}{p(-,-\mid w)\cdot p(-,+\mid w)}\bigg)=\psi (w).\end{aligned}\]
Let us check that the functions ${\psi _{-}}(w),{\psi _{+}}(w)$ and $\psi (w)$ possess the respective properties. By the definition of the local symmetry
(19)
\[ p(s,y\mid w)=p\big(s,-y\hspace{0.1667em}\big|\hspace{0.1667em}{w^{\ast }}\big),\hspace{1em}\forall w\in {\mathcal{A}_{1}^{4m}}.\]When estimating the generalized logits ${\Lambda _{\tau }}(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),\hspace{0.1667em}s,y\in {\mathcal{A}_{1}^{2m}},$ and hence in the sequel we deal with functions $\psi (s,y)$, $\psi :{\mathcal{A}_{1}^{k}}\times {\mathcal{A}_{1}^{k}}\to \mathbf{R}$, $k:=2m$.
Let $J\subset K:=\{1,\dots ,k\}$. Define the conjugate set ${J_{\ast }}$ of the set J by
For a given sequence $s=({s_{1}},\dots ,{s_{k}})\in {\mathcal{A}_{1}^{k}}$, denote
Any function $\psi :{\mathcal{A}_{1}^{k}}\times {\mathcal{A}_{1}^{k}}\to \mathbf{R}$ has the unique representation
where summation is over all subsets of K (including the empty set ∅), ${a_{{J^{\prime }}J}}={a_{{J^{\prime }}J}}(\psi ),\hspace{0.2778em}{J^{\prime }},J\subset K,$ are free parameters determining the function ψ. In general, there are ${4^{k}}$ free parameters.
(28)
\[ \psi (s,y)=\sum \limits_{{J^{\prime }},J\subset K}{a_{{J^{\prime }}J}}\hspace{0.2778em}{s^{{J^{\prime }}}}{y^{J}},\hspace{1em}s,y\in {\mathcal{A}_{1}^{k}},\]For a symmetric (antisymmetric) function ψ, we have
Consequently, in the case of the symmetric ψ, for all $s,y\in {\mathcal{A}_{1}^{k}}$,
and hence
If ${J_{\ast }}=J$ and ${J^{\prime }_{\ast }}={J^{\prime }}$ (i.e. the both subsets ${J^{\prime }}$ 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 ${a_{{J^{\prime }}J}}$ values in this case. Let ${k_{\ast }}={k_{\ast }}(k)$ denote the total number of the self-conjugate subsets of K.
(29)
\[ \psi (\stackrel{\gets }{s},-\stackrel{\gets }{y})=\psi (s,y)\hspace{1em}\big(\mathrm{respectively},\psi (\stackrel{\gets }{s},-\stackrel{\gets }{y})=-\psi (s,y)\big).\](30)
\[ \sum \limits_{{J^{\prime }},J\subset K}{a_{{J^{\prime }}J}}\hspace{0.2778em}{s^{{J^{\prime }}}}{y^{J}}=\sum \limits_{{J^{\prime }},J\subset K}{(-1)^{|J|}}{a_{{J^{\prime }}J}}\hspace{0.2778em}{s^{{J^{\prime }_{\ast }}}}{y^{{J_{\ast }}}}=\sum \limits_{{J^{\prime }},J\subset K}{(-1)^{|J|}}{a_{{J^{\prime }_{\ast }}{J_{\ast }}}}\hspace{0.2778em}{s^{{J^{\prime }}}}{y^{J}},\](31)
\[ {a_{{J^{\prime }}J}}={(-1)^{|J|}}{a_{{J^{\prime }_{\ast }}{J_{\ast }}}},\hspace{1em}{J^{\prime }},J\subset K.\]Let τ be some total order (enumeration of elements) in the class of pairs $({J^{\prime }},J)$ of the set K. Equations (31) imply that, for not self-conjugate pairs $({J^{\prime }},J)$, $({J^{\prime }},J)\ne ({J^{\prime }_{\ast }},{J_{\ast }})$, values of the coefficients ${a_{{J^{\prime }},J}},\hspace{0.2778em}\tau ({J^{\prime }},J)<\tau ({J^{\prime }_{\ast }},{J_{\ast }}),$ uniquely determine values of the remaining coefficients ${a_{{J^{\prime }},J}},\hspace{0.2778em}\tau ({J^{\prime }},J)>\tau ({J^{\prime }_{\ast }},{J_{\ast }})$. Define
From (28), (31), (32) and (33) we derive a general parametric form of a symmetric function ψ:
It has
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^{\prime }},J)$, these equations hold if and only if ${a_{{J^{\prime }},J}}=0$. Thus, the first summand in (34) and in (35) disappears giving the function
with
free parameters.
For $m=1$, ${k_{\ast }}=2$, thus ${k_{A}}=({4^{2}}-{2^{2}})/2=6$ and ${k_{S}}={k_{\ast }^{2}}+{k_{A}}=10$. Then symmetric (34) and antisymmetric (36) functions in a general form are given by
\[\begin{aligned}{}{\psi _{S}}(s,y)& ={a_{\varnothing \varnothing }}+{a_{\varnothing \{12\}}}{y_{1}}{y_{2}}+{a_{\{12\}\varnothing }}{s_{1}}{s_{2}}+{a_{\{12\}\{12\}}}{s_{1}}{s_{2}}{y_{1}}{y_{2}}\\ {} & \hspace{1em}+{a_{\varnothing \{1\}}}({y_{1}}+{y_{2}})+{a_{\{1\}\varnothing }}({s_{1}}-{s_{2}})\\ {} & \hspace{1em}+{a_{\{1\}\{1\}}}({s_{1}}{y_{1}}-{s_{2}}{y_{2}})+{a_{\{1\}\{2\}}}({s_{1}}{y_{2}}-{s_{2}}{y_{1}})\\ {} & \hspace{1em}+{a_{\{1\}\{12\}}}({s_{1}}{y_{1}}{y_{2}}-{s_{2}}{y_{1}}{y_{2}})+{a_{\{12\}\{1\}}}({s_{1}}{s_{2}}{y_{1}}+{s_{1}}{s_{2}}{y_{2}}),\\ {} {\psi _{A}}(s,y)& ={a_{\varnothing \{1\}}}({y_{1}}-{y_{2}})+{a_{\{1\}\varnothing }}({s_{1}}+{s_{2}})\\ {} & \hspace{1em}+{a_{\{1\}\{1\}}}({s_{1}}{y_{1}}+{s_{2}}{y_{2}})+{a_{\{1\}\{2\}}}({s_{1}}{y_{2}}+{s_{2}}{y_{1}})\\ {} & \hspace{1em}+{a_{\{1\}\{12\}}}({s_{1}}{y_{1}}{y_{2}}+{s_{2}}{y_{1}}{y_{2}})+{a_{\{12\}\{1\}}}({s_{1}}{s_{2}}{y_{1}}-{s_{1}}{s_{2}}{y_{2}}),\end{aligned}\]
respectively.Remark 1.
An ordered sequence of symbols $x=\overrightarrow{x}$ is said to be palindromic iff $\overrightarrow{x}=\stackrel{\gets }{x}$. We refer to the mapping $\overrightarrow{x}\to \stackrel{\gets }{x}$ as palindromic transformation. In particular, for a DNA sequence x, the sequence $(x,{x_{\ast }^{\ast }})$ (here ${x_{\ast }^{\ast }}={({x_{\ast }})^{\ast }}={({x^{\ast }})_{\ast }}$) is palindromic and for a palindromic DNA sequence x, we have $s(x)=s({x_{\ast }})$ and $y(x)=-y({x_{\ast }})$. Note that the mapping ${\mathcal{A}_{1}}\to \{+1,-1\}$ is a palindromic transform of the binary alphabet ${\mathcal{A}_{1}}$. Thus, the transform $y(x)\to y({x_{\ast }})$ is a superposition of two palindromic transforms: the transform of ordering of the sequence $y(x)$ elements and the transform of their alphabet.
Palindromic distributions are defined as those invariant under some palindromic operation. For instance, palindromic Bernoulli distributions (Marchetti and Wermuth, 2016) and palindromic Ising models (Marchetti and Wermuth, 2017) are invariant with respect to palindromic transforms of the alphabet. Formulas (34) and (36) are analogues of the characterization of palindromic Bernoulli distribution in terms of log-linear parameters of multivariate Bernoulli distribution given in Marchetti and Wermuth (2016).
3 Statistical Analysis
In this section, the first-order local symmetry of the longest non-coding sequences of bacterial genomes is tested by making use of its characterization in terms of generalized logits. A special regression-type probabilistic structure is imposed on the data.
3.1 Regression-Type Probabilistic Structure of the Data
Let us introduce the following data structure of the observed sequence $x\in {\mathcal{A}^{n}}$ with $n=({n_{m}}+1)\cdot (m+1)-1$, the quantity ${n_{m}}$ being an integer:
where ${v_{\ell }}:={x_{(m+1)\ell }}$ is a response variable and ${z_{\ell }}={x_{{U_{m}}((m+1)\ell )}}\in {\mathcal{A}^{2m}}$ is a vector of explanatory variables, $\ell \in S$.
(38)
\[ \mathcal{D}:=\big\{({v_{\ell }},{z_{\ell }}),\ell \in S\big\},\hspace{1em}S={S_{n,m}}=\{1,2,\dots ,{n_{m}}\},\]Assumption (Am):
Assumption (Am) ensures that usual conditions of the generalized logit model with the response variable $v\in \mathcal{A}$ and the vector $z\in {\mathcal{A}^{2m}}$ of explanatory variables are satisfied, see Agresti (1990), Stokes et al. (2001).
Remark 2 (Compatible evolutionary models).
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={\mathcal{N}^{\circ }}$ (to be precise, for all shifts $((m+1)S+\ell )\cap {\mathcal{N}^{\circ }}$ of the set of central nucleotides $(m+1)S$ by ℓ, $\ell =1,\dots ,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.
3.2 Testing of Local Symmetry
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.
Fig. 1
The length distribution density of the longest non-coding sequences of bacteria genomes plotted in a logarithmic scale.
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\to \mathrm{sgn}(S){\log _{2}}(1+|S|)$) for testing the significance of the coefficients of response functions ${\psi _{-}},{\psi _{+}}$ 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.
Fig. 2
Distribution of the logarithmized Student statistic of the 1st response function coefficients: the first 6 coefficients represent the antisymmetric part, the last 10 – the symmetric part (expected to be null).
Fig. 3
Distribution of the logarithmized Student statistic of the 2nd response function coefficients: the first 6 coefficients represent antisymmetric part, the last 10 – the symmetric part (expected to be null).
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).
Fig. 4
Distribution of the logarithmized Student statistic of the 3rd response function coefficients: the first 6 coefficients represent antisymmetric part (expected to be null), the last 10 – the symmetric part.
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]$.
3.3 Concluding Remarks
Elements of DNA sequences x are treated as random variables taking values from the alphabet $\mathcal{A}:=\{\texttt{A},\texttt{C},\texttt{G},\texttt{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.
-
2. The deviations from the local symmetry of the non-coding sequences are pretty symmetric: the sample distributions of estimates of the model parameters that should vanish in case of the local symmetry are very close to symmetric one. Apparently this symmetry is related to intra-strand frequency distribution parity noticed in Powdel et al. (2009).
-
3. As a by-product of the statistical analysis of the local symmetry, we show that distributions of adjacent nucleotides are not independent even for the non-coding sequences of bacteria genoms. Hence independent evolution models (see, e.g. Faith and Pollock, 2003; Marin and Xia, 2008) are not consistent with the data of bacteria genomes.
Further work
A natural next step is to study higher order asymmetry patterns. Under assumptions (Am) with $m=2$, for the statistical analysis of the second order local asymmetries the saturated generalized logit model can be applied. Then the analysis is based on 5-dimensional contingency tables (1024 cells). Hence for the data of the longest non-coding bacterial sequences, the average cell frequency in the contingency tables is less than 3, thus indicating their sparsity. A straightforward solution of the sparsity problem by joining all non-coding sequences of each genome seems to be inappropriate because of heterogeneity of DNA sequences (see, e.g. Cristadoro et al., 2018). Special statistical methods are needed to deal with both the sparsity and heterogeneity.
