Metrics and Harmonic Analysis on DNA

TL;DR

This paper introduces a mathematical algorithm based on the crystal basis model to analyze DNA and RNA sequences, enabling sequence comparison and Fourier analysis to explore genetic similarities and differences.

Contribution

It presents a novel method to describe genetic sequences as integer strings, allowing for sequence comparison and harmonic analysis within the crystal basis framework.

Findings

Defined a new sequence distance measure.

Applied Fourier transform to genetic sequences.

Compared sequences from different species.

Abstract

A mathematical algorithm to describe DNA or RNA sequences of $N$ nucleotides by a string of $2N$ integers numbers is presented in the framework of the so called crystal basis model of the genetic code. The description allows to define a not trivial distance between two sequences and to perform the Discrete Fourier Transform on a sequence. Using the definition of distance the sequence of $\beta$-globin for Homo, Gallus and Opossum are compared. The Discrete Fourier Transform for some simple examples is computed.

Figures

Tables

Table 1. Table 1: Examples of labeling of sequences with the same content, but different position, of the 5 nucleotides.

Sequence	$J_{3,H}$	$J_H^5$	$J_H^4$	$J_H^3$	$J_H^2$
UUCCC	1/2	5/2	2	3/2	1
CUUCC	1/2	3/2	1	1/2	0
CCUCU	1/2	1/2	1	1/2	1
CCCUU	1/2	1/2	1	3/2	1
CCUUC	1/2	1/2	0	1/2	1
CUCUC	1/2	1/2	0	1/2	0

Table 2. Table 2: Distance between the 92 nucleotides sequence of the β − g ​ l ​ o ​ b ​ i ​ n 𝛽 𝑔 𝑙 𝑜 𝑏 𝑖 𝑛 \beta-globin of Homo, Gallus and Opossum. In the first column the considered couple of species.

	$d_H$	$d_K$	$d_{CB}$
Homo/Gallus	24 a	$10\alpha +14\beta$	$105{\alpha }_H+129{\alpha }_V$
Homo/Opossum	24 a	$11\alpha +13\beta$	$214{\alpha }_H+78{\alpha }_V$
Gallus/Opossum	28 a	$6\alpha +22\beta$	$313{\alpha }_H+125{\alpha }_V$
	$d_H^{st}$	$d_K^{st}$	$d_{CB}^n$
Homo/Gallus	$\approx 0,26$	$\approx 0,32$	$\approx$ 20,25
Homo/Opossum	$\approx 0,26$	$\approx 0,26$	$\approx$ 33,43
Gallus/Opossum	$\approx 0,30$	$\approx 0,10$	$\approx$ 43,20

Equations

$$\text{mbox}C\equiv (+\frac{1}{2},+\frac{1}{2})\text{mbox}T/U\equiv (-\frac{1}{2},+\frac{1}{2})\text{mbox}G\equiv (+\frac{1}{2},-\frac{1}{2})\text{mbox}A\equiv (-\frac{1}{2},-\frac{1}{2})$$

$$|J_1-J_2|\le J\le J_1+J_2$$

$$J_a=J_{(a,Max)},J_{(a,Max)}-1,J_{(a,Max)}-2,\dots \ge 0J_{(a,Max)}=N/2$$

$$J_{3,a}=-J_a,-J_a+1,\dots ,J_a$$

$$\text{mbox}\{N-nucleotides\}⟹\mid {\mathbf{J}}_H{\mathbf{J}}_V⟩=\mid J_{3,H},J_{3,V};J_H^N,J_V^N;\dots ;J_H^2,J_V^2⟩$$

$$\mid CGA⟩=\mid J_{3,H},J_{3,V};J_H^3,J_V^3;J_H^2,J_V^2⟩=\mid \frac{1}{2},-\frac{1}{2};\frac{1}{2},\frac{1}{2};1,1⟩$$

$$\Bigl{(}\frac{1}{2},\frac{1}{2}\Bigl{)}\ \otimes\ \Bigl{(}\frac{1}{2},\frac{1}{2}\Bigl{)}\ \otimes\ \Bigl{(}\frac{1}{2},\frac{1}{2}\Bigl{)}=\Bigl{(}\frac{1}{2},\frac{1}{2}\Bigl{)}\otimes\Biggl{[}\Bigl{(}\frac{1}{2},\frac{1}{2}\Bigl{)}\ \otimes\ \Bigl{(}\frac{1}{2},\frac{1}{2}\Bigl{)}\Biggl{]}$$

$$2J_{3,H}=N_C+N_G-N_U-N_A\equiv N_{CG}-N_{UA}2J_{3,V}=N_C+N_U-N_G-N_A\equiv N_{CU}-N_{GA}$$

$$J_{3,H}=N_{GC}-\frac{N}{2}$$

$$J_{3,V}=N_{CU}-\frac{N}{2}$$

$$J_{3,H}\equiv N_{CG}-\frac{N}{2}=3-4=-1J_{3,V}\equiv N_{CT}-\frac{N}{2}=3-4=-1$$

$$AGGACTTA⟶WSSWSWWW⟶WS\cdot \cdot WW⟶W\cdot \cdot \cdot W⟹J_H^8=1$$

$$AGGACTT⟶WSSWSWW⟶WS\cdot \cdot W⟶W\cdot \cdot \cdot ⟹J_H^7=1/2$$

$$AGGACT⟶WSSWSW⟶WS\cdot \cdot ⟹J_H^6=1$$

$$AGGAC⟶WSSWS⟶WS\cdot S⟹J_H^5=3/2$$

$$AGGA⟶WSSW⟶WS\cdot ⟹J_H^4=1$$

$$AGG⟶WSS⟹J_H^3=3/2AG⟶WS⟹J_H^2=1$$

$$AGGACTTA⟶RRRRYYYR⟶RRRRYY\cdot ⟹J_V^8=3$$

$$AGGACTT⟶RRRRYYY⟶RRRRYYY⟹J_V^7=7/2$$

$$⟹J_V^6=3⟹J_V^5=5/2⟹J_V^4=2⟹J_V^3=3/2⟹J_V^2=1$$

$$\mid AGGACTTA⟩=\mid -1,-1;1,3;\frac{1}{2},\frac{7}{2};1,3;\frac{3}{2},\frac{5}{2};1,2;\frac{3}{2},\frac{3}{2};1,1⟩$$

$$\mid ATGGTGCACC⟩=\mid 1,0;3,2;\frac{5}{2},\frac{3}{2};2,1;\frac{5}{2},\frac{3}{2};2,1;\frac{3}{2},\frac{3}{2};2,1;\frac{3}{2},\frac{1}{2};1,1⟩$$

$${\lambda }_{0,a}\equiv J_{3,a}{\lambda }_{k,a}\equiv J_{\alpha }^{N-k+1}k=1,2,\dots ,N-1a=H,V$$

$$ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGT$$

$$GGGGCAAGGTGAACGTGGATGAAGTTGGTGGTGAGGCCCTGGGCAG$$

$$d^1(1,2)=\sum _{i=H,V}{\alpha }_i\left(\sum _{k=0}^{N-1}|{\lambda }_{k,i}^1-{\lambda }_{k,i}^2|\right)$$

$${\Psi }_N\equiv \sum _{i=H,V}{\alpha }_i\left(\sum _{k=0}^{N-1}{\lambda }_{k,i}\right)$$

$$d^E(1,2)=\sum _{i=H,V}{\alpha }_i\sqrt{\sum _{k=0}^{N-1}\left[{\left({\lambda }_{k,i}^1-{\lambda }_{k,i}^2\right)}^2\right]}$$

$${\lambda }_k^c={\lambda }_{k,H}+i{\lambda }_{k,V}$$

$$d_{stat}(1,2)=\sum _{p_1}^{p_2}\frac{p_x}{\sqrt{X}}$$

$$d_H(1,2)=3a{d}_H(3,4)=2a{d}_H(5,6)=a$$

$$d_K(1,2)=\alpha +2\beta d_K(3,4)=2\beta d_K(5,6)=\beta$$

$$d_{CB}(1,2)={\alpha }_H+2{\alpha }_V{d}_{CB}(3,4)=4{\alpha }_H+{\alpha }_V{d}_{CB}(5,6)={\alpha }_V$$

$$CAUC=\mid 0,1;1,1;\frac{1}{2},\frac{1}{2};0,0⟩$$

$$AUGU=\mid -1,0;1,1;\frac{3}{2},\frac{1}{2};1,1⟩$$

$$CUGU=\mid 0,1;0,1;\frac{1}{2},\frac{1}{2};0,1⟩$$

$$d_H^{st}\equiv \frac{d_H}{aN}$$

Full text

Keywords: DNA, RNA, nucleotide sequences, crystal basis model, genetic code, distance, discrete Fourier transform

1 Introduction

In the last years, complete genomes have been identified and the expressions of tens of thousands of genes for many organisms have been identified.

As more genomes become sequenced and more genes are located, the identification and characterization of gene functions are the great new challenges of the functional genomics.

The advent of large-scale “transcriptomic” and “proteomic” technologies is providing each day huge amounts of data, pushing to search new approaches to discover and characterize gene expression patterns and, more generally, nucleotide sequences in DNA as well as in RNA. Moreover, with the advent of computers, new scientific branches have arisen, the bioinformatics and next the bioengineering. As emblematic example, let us quote the MIT programming language to design DNA circuit Nielsen and al. [2016], which very recently has gathered a very wide appeal in the world media for its implication on the genetic manipulations. So it appears evident the need to translate any DNA or RNA sequence in a mathematical language, in order to allow its study by means of the tools of the mathematics and/or of the artificial intelligence.

The standard representation of nucleotide sequences is by letter-series using the four letter “alphabet” {C,T/U,G,A} of the DNA/RNA. However, the availability of DNA and RNA data, each day extremely fast growing, has prompted in the last decades, a research toward mathematical and graphical representation of nucleotide series in order to quickly extract relevant information. At our knowledge one of the first alternative representation has been proposed in Hamori and Ruskin [1983], where a method has been described through which any sequence is represented by a three dimensional curve. In Nandy et al. [2006] several graphical and numerical methods are reviewed and compared.

Another very relevant question arising when nucleotide sequences are compared is the characterization of their differences, in other words their “distance”.

The first and commonly used measure between two sequences is the Hamming distance Hamming [1950]. The Hamming distance between two sequences of equal length is the number of positions at which the corresponding symbols are different. In other words it measures the minimum number of substitutions required to change one sequence into the other. The most simple generalization of the Hamming distance is the Kimura 2-parameter distance Kimura [1980] assigning different parameter if a pyrimidine Y (Y=C,T/U) (a purine R (R=A,G)) nucleotide is replaced by the other one, transition, or by a purine R, transversion, (resp. pyrimidine Y).

Given two symbolic sequences of different length, their relative distance, that is a measure of the amount of difference between them, is commonly computed by the Levenshtein distance Levenshtein [1966], which is defined as the minimum number of “edit operation” needed to transform one sequence into the other, the “edit operation” being insertion, deletion, or substitution of a single symbol.111The Levenshtein distance between two strings of the same length is strictly less than the Hamming distance.

As the measure of the distance between two sequences is extremely relevant in the theory of the evolution, many and more sophisticated measures have been introduced, whose discussion is beyond the aim of the present paper. However these measures, some of which are not strictly measures, depend only on the positions or, at most, by some global property, e.g. the GC content.222In Hosangadi [2012] several distance measures for sequences have been examined and compared.

In the present work we introduce a parametrization to describe a DNA or RNA sequence, which provides information on the distribution of the nucleotides in the sequence. The mathematical scheme establishes a correspondence between any sequence of N nucleotides and a set of 2N half-integer and integer numbers333Clearly one can deal only with integer numbers by multiplying our parameters by 2., so DNA/RNA can be considered as a discrete function with codomain in the set of integer number.

Let us point out some features of our parametrization.

it allows to introduce a well defined distance between any two sequences; 2. 2.

it allows to perform a Discrete Fourier Transform (DFT) and, consequently, to exploit all the power of the harmonic analysis; 3. 3.

a transition, a transversions and, more generally, any mutation can be defined as tensor operator, so allowing to have a mathematical description where the mutation of a nucleotide is position dependent.

The aim of the paper is not to discuss specific problems, but to list, surely in a not exhaustive way, the possibilities that our mathematical description provides. The paper is organized as follows: in Sect. 2 we concisely recall the mathematical scheme which is at the base of the parametrization; in Sect. 3 we describe the algorithm which allows to represent the nucleotide sequence in terms of a discrete function taking value in the set of half-integer and integer numbers, discussing the main features of the parametrization, presenting also a few simple examples, and in Sect. 4 we give a short summary and some highlights on future developments, in particular on the above quoted 3-rd point.

2 Mathematical parametrization of DNA/RNA sequences

2.1 Basic mathematical structure of the crystal basis model

In order to label a sequence of N nucleotides with 2N half-integer and integer numbers we have to briefly recall the model of genetic code introduced in Frappat et al. [1998] and called “crystal basis model”. In the model the 4 nucleotides, denoted by the four letter set {C,T/U,G,A} in DNA/RNA, are assigned to the 4-dim fundamental irreducible representation (irrep.) of $\mathcal{U}_{q\to 0}(sl(2)\oplus sl(2))$, the two $sl(2)$ being denoted as $sl_{H}(2)$ and $sl_{V}(2)$.444The special basis of $\mathcal{U}_{q}(sl(n)$ in the limit $q\to 0$ is called a crystal basis, see Kashiwara [1990]. This explains the name of the model. The fundamental irrep., of dimension 2, is labelled by $(J_{H}=1/2,J_{V}=1/2)$ and the states are identified by the following assignment for the value of the third component of the operator $J_{H}$ ($J_{V}$) denoted by $J_{3,H}=\pm{\textstyle{\frac{1}{2}}}$ (resp. by $J_{3,V}=\pm{\textstyle{\frac{1}{2}}}$):

[TABLE]

In Kashiwara [1990] it has been shown that the tensor product of two crystal bases labelled by $J_{1}$ and $J_{2}$ can be decomposed into a direct sum of crystal bases labelled, as in the case of the tensor product of two standard $sl(2)$ irreducible representations, by a set of integer ($J_{1}+J_{2}$ integer) or half-integer ($J_{1}+J_{2}$ half-integer) numbers $J$ such that

[TABLE]

The new peculiar and crucial feature is that, in the limit $q\to 0$, the basis vectors of the $J$-space are pure states, that is they are the product of a state belonging to the $J_{1}$-space and of a state belonging to the $J_{2}$-space, while in the case of $sl(2)$ or of ${\cal U}_{q}(sl(2))$ they are linear combinations of states belonging to the $J_{1}$-space and to the $J_{2}$-space with coefficients called respectively Clebsch-Gordan coefficients and $q$-Clebsch-Gordan coefficients. It should be noted that in the case of $\mathcal{U}_{q}(sl(2))$, then in particular for $\mathcal{U}_{q\to 0}(sl(2))$, the tensor product is not symmetric and one has to specify which is the first irrep. in the tensor product, i.e. an order has to be fixed. In the framework of DNA/RNA the order appears in a very natural way as the reading frame moves toward the $5\prime\rightarrow 3\prime$ direction.

It follows that a state of an irrep. contained in the tensor product of $N$ fundamental irreps. is an N-ordered sequence of states belonging to the fundamental irreps., i.e. an ordered sequence of N letters belonging to the four letter “alphabet” {C,T/U,G,A}, and is labelled by $J_{a}$ and $J_{3,a}$, $a=H,V$:

[TABLE]

As for $N>2$ an irrep. labelled by $J_{a}\neq J_{(a,Max)}$ is degenerate, i.e. there are several irreps. labelled by the same values of $J_{a}$, we need to identify an irrep. by other $2(N-2)$ labels, which we choose to be the labels of the tensor product of $k$ ($k=N-1,N-2,\ldots,2$) fundamental irreps., obtained by removing $(N-k)$ irreps. in the tensor product, starting from the right. In the following we denote the J labels by un upper index $k$.

Let us point out that the above given prescription to remove $(N-k)$ irreps. starting from the right is just a convention. One may equally well choose to remove $(N-k)$ irreps. starting from the left. However the results of the two prescriptions are not, in general, equal and the choice of the prescription modifies, e.g., the distance between sequences, see next Section. Our chosen prescription is suggested by the direction of the transcription process from DNA to RNA.

Summarizing, making use of the properties of crystal basis, any string of $N$ binary label (spin) $\;x\in\{\pm\equiv R,Y\;\mbox{or}\;\equiv S,W\;\mbox{or}\;\equiv C,U\;\mbox{or}\;\equiv G,A\;\mbox{or}\;\equiv C,G\;\mbox{or}\;\equiv U,A\}\;$ can be seen as a state of an irrep. contained in $N$-fold tensor product of the the 2-dim fundamental irrep. (labelled by $J=1/2$) of $\;sl_{q\to 0}(2)\;$ whose state are labelled by $\;J_{3}=\pm 1/2\equiv\pm\,$. In the case of the nucleotides, due to the presence of two $sl(2)$, eq.(1), the number of needed labels has to be doubled. Therefore, in the most general case, a sequence of N nucleotides can be identified by the following $2N$ labels.

[TABLE]

where $J^{N}_{a}$ ($a=H,V$) labels the irrep. of of $\mathcal{U}_{q\to 0}(sl_{a}(2))$ which the state belongs to, $J_{3,a}$ is the value of the 3rd component of the operator $J_{a}^{N}$. Note that $J^{k}$ and $J_{3}$ are integer for even $k$ and half-integer for odd $k$. Obviously a H spin flip, i.e. the change by $\pm 1$ of $J_{H}$, (V spin flip) corresponds, respectively, to a biological transition (transversion).

For example, the codon CGA is represented as

[TABLE]

The labels have been computed considering the codon as a composite state resulting from the 3-fold tensor product of the fundamental irrep., tensor product computed following the below indicated order

[TABLE]

2.2 Labelling of an ordered sequence of N nucleotides

In the following we give the algorithm for the computation of the labels $J_{a}^{k}$ appearing in eq.(5) for a sequence of $N$ nucleotides.

Denoting by $N_{X}$ the number of nucleotides of type $X\in\{C,U,G,A\}$555Clearly for DNA U should be replaced by T., we have:

[TABLE]

The GC content, i.e. the number $N_{GC}=N_{C}+N_{G}$ of G and C nucleotides is clearly related to the value of $J_{3,H}$ by

[TABLE]

while the asymmetry, i.e. the difference of the number of the complementary bases or the excess of purines, is related to the value of $J_{3,V}$

[TABLE]

Let us recall that we have $|J_{3,H}|\leq J_{H}^{N}\leq N/2$ and $|J_{3,V}|\leq J_{V}^{N}\leq N/2$, where the equality for H (V) labels holds for a sequence of only C and G or U and A (resp. C and U or G and A) nucleotides.

• •

The algorithm - From eqs.(9)-(10) we can compute the values of $J_{3,H}$ and $J_{3,V}$ corresponding to a fixed sequence of N nucleotides. In the following we describe, in a very didactic way, the algorithm to compute the values of $J_{H}^{N},J_{V}^{N};\ldots;J_{H}^{2},J_{V}^{2}$ corresponding to the sequence of N nucleotides. The value of $J^{N}_{H}$ is computed replacing, in the sequence, the nucleotides C, G by the symbol $S$ and A, U/T by the symbol W (the standard symbols denoting the two couples of nucleotides with respectively Strong (Weak) hydrogen bond), then taking away the $W$ elements, which are at right of the $S$ elements, contracting each of them with a $S$ element on the left, and then summing the numbers of left $S$ and $W$, which of course are, respectively, at the right of $W$ and at the left of $S$. In other words this value is computed deleting all the couples $SW$ of first neighbours which are present in the sequence and iterating the procedure, on the sequence obtained by deleting the first couples $SW$, as long as no $W$ appears on the right of any $S$. We refer to the elements which are deleted in the procedure as contracted. The value $J^{N}_{V}$ is computed replacing in the sequence the nucleotide C, U/T by $Y$ and A,G by $R$ (the standard symbols for pyrimidines and purines) and then proceeding in the same way as for the computation of $J^{N}_{H}$, but now deleting all the couples $YR$, that is $R$ ($Y$) takes the role of $W$ (resp. $S$). Then we proceed to compute the value of $J_{H}^{k}$ and of $J_{V}^{k}$ ($k=N-1,\ldots,2$) deleting in the original sequence $N-k$ nucleotides from the right and applying to the obtained $k$-nucleotide sequence the above given prescriptions. The algorithm, nevertheless its apparent complexity, can be easily implemented in computers and is in detail illustrated in the example below.

Example: Let us consider the sequence AGGACTTA

Calculation of $J_{3,H}$ and $J_{3,V}$. From eqs.(9)-(10), we find

[TABLE] 2. 2.

Calculation of $J_{H}^{k}$, ($k=8,7,\ldots,2$). We replace a contracted couple $SW$ by $\cdot$.

[TABLE] 3. 3.

Calculation of $J_{V}^{k}$, ($k=8,7,\ldots,2$). We replace a contracted couple $YR$ by $\cdot$.

[TABLE]

Summarizing, the sequence is labelled, the order of the labels being specified in eq.(5), by:

[TABLE]

Remarks:

• •

$J_{a}^{n}$ differs by $J_{a}^{n\pm 1}$ by $\pm{\textstyle{\frac{1}{2}}}$, ($n=2,\ldots,N-1$. $a=H,V$);

• •

if $J_{a}^{n}=J_{a}^{n+2k}$ ($n\geq 2,k=1,2,\ldots$) the numbers of S and W nucleotides, between the n-th and (n+2k)th-position, are equal and the W nucleotides are arranged on the left of the S nucleotides (a = H) or the numbers of Y and R nucleotides are equal and the R nucleotides are arranged on the left of the Y nucleotides (a = V);

• •

the values of the labels take somehow into account the arrangement of the nucleotides. To further illustrate this point we report in Table 1 the value computed with a string of 5 pyrimidines formed by 2 U and 3 C, in different position. For all the strings the values of the V-labels are: $J_{3,V}=5/2,J_{V}^{5}=5/2,J_{V}^{4}=2,J_{V}^{3}=3/2,J_{V}^{2}=1$.

Let us apply the above described algorithm to two more simple, but may not irrelevant, cases.

let us consider the $N=10$ sequence $ATGGTGCACC$ which corresponds to the first 10 nucleotides of the exon 1 in Human $\beta-globin$ gene, see below.666We use this sequence as it has been used in many examples in Nandy et al. [2006] in order to allow an immediate comparison between our and other descriptions existing in the literature. The sequence is labelled, the order of the labels being specified in eq.(5), by:

[TABLE]

We report in Fig. 2 the diagram of the labels of the above string. In the following, to simplify the notation, we denote

[TABLE] 2. 2.

let us consider the $N=92$ coding sequence of the human $\beta-globin$ gene, see Table 2 of Qi and al [2012] (http://www.hindawi.com/journals/tswj/2012/104269/tab2/)

[TABLE]

We report in Fig. 4 the labels of the above sequence. In Figg. 6-6 the labels for the $N=92$ sequence of the $\beta-globin$ for Gallus and Opossum are reported, the nucleotide sequences being omitted to save space, see Table 2 of Qi and al [2012].

Once a discrete function is associated to a DNA/RNA sequence, one may study this function from many points of view. Let us cite just a few: one can

• •

perform discrete circular convolution looking for symmetry, e.g. for function even or odd under $z$ ($z\in{\bf Z}$) circular symmetry,

• •

look for other kinds of symmetry, in the whole sequence or in some subsequence, for example:

In Fig. 2, neglecting the point n = 0, and subtracting 2 to the values of the H labels, one remarks that the function shows an inverse specular symmetry, i.e. is odd, around the point n = 5, that is $\lambda_{mp-n,H}=-\lambda_{mp+n,H}$, where $\lambda_{mp,H}$ is the value of the H label in the middle position ($mp=5$). The meaning of this symmetry is that the number of not contracted S and W nucleotides between the positions $N+1-mp+n=6+n$ and $N+1-mp-n=6-n$ is given by $\lambda_{mp-n,H}-\lambda_{mp+n,H}$. Of course this information can be obtained out from an analysis of the labels $\lambda_{k,H}$, but the graphical representation may help to more quickly extract the information; 2. 2.

In Fig. 4, (human $\beta-globin$) one can also notice that: in the case of V labels, considering approximately the range of labels from 7 to 70 and subtracting 5 to the values of $\lambda_{k,V}$ and in the case of H labels, considering approximately the range of labels from 3 to 50 and subtracting 7 to the values of $\lambda_{k,H}$, an approximate specular symmetry does show around the middle point of each series of labels; 3. 3.

In Fig. 6, ($\beta-globin$ of Gallus) the 9 labels on the right of the label $\lambda_{28,V}$ have the same value of the 9 labels on the left. The meaning of this symmetry can be understood on the light of the above remark; 4. 4.

In Fig. 6, ($\beta-globin$ of Opossum) the 8 labels on the right of the label $\lambda_{73,H}$ have the opposite value of the 8 labels on the left subtracting to all labels the value 2,5; 5. 5.

In Fig. 6, the first 45 positive H labels fit a simple linear regression with a factor $R^{\approx}0,9$. This means that first 45 nucleotides are almost or purine or yrimidine. Indeed a linear segment in a graphic of the H (V) labels means that the nucleotides are all S or W nucleotides (resp. Y or R).

• •

study the behavior under shifts of a sequences of length $M<N$;

• •

analyze the “time reversal” of the function, that is the function obtained by reading the DNA frame in the inverse direction $3\prime\rightarrow 5\prime$ direction;

• •

sample the sequence with windows of variable length, suggested by preliminary analysis of the discrete function or by the interest to focus on some particular sub-sequences.

Let us remark that from Fig. 4 a some kind of correlation does show between the H and V labels, the correlation factor between the two sets of labels (not including $J_{3,H}$ and $J_{3,V}$) being $r\approx 0,73$. However, at the present stage of our analysis, we cannot state the this weak correlation is significant. Indeed the correlation between the $J_{H}$ and $J_{V}$ labels for the $\beta-globin$ sequence for Gallus is $r\approx 0,81$ and for Opossum $r\approx 0,53$.

Let us emphasize once more that the knowledge of all the labels $\;J_{3,H},J_{3,V};J_{H}^{N},J_{V}^{N};\ldots;J_{H}^{2},J_{V}^{2}\;$ allows to reconstruct the sequence in terms of the four letter nucleotide alphabet and conversely.

3 Metrics and harmonic analysis

Once established a correspondence between any DNA/RNA sequence and the value of a discrete function in ${\bf Z}$, according to the procedure illustrated in the previous Section, it is possible to go further, in particular to define a distance between two sequences and to perform the Discrete Fourier Transform (DFT) of a sequence.

3.1 Distance between two sequences of equal length

The most simple definition of distance, based on two parameters $\alpha_{i}$ is the following one (1-norm): given two $N$-elements strings the distance is defined by

[TABLE]

where 1,2 denotes the two sequences and $\alpha_{i}$ is a real parameter denoting the “scale” of the distance in the $i-th$-direction. Rigorously speaking eq.(18) is the 1-norm metrics on the function

[TABLE]

computed for the two sequences. Eq.(19) is the natural distance in the “Taxicab geometry” and is also called the “Manhattan” or “City block distance” , up tp the parameters $\alpha_{i}$.

Let us remark that this distance can be seen as the minimum number of different states between the states corresponding to the first and second sequence, assuming the change of the parameters $\lambda_{k,i}$ by step of one H or V spin flip.

Of course one may also introduce an Euclidean distance

[TABLE]

which can be read as the 2-norm metrics on the function eq.(19). However

One may consider $\lambda_{H}$ and $\lambda_{V}$ as respectively the real an imaginary part of a complex number so establishing a correspondence between a nucleotide sequence and a set of N complex numbers

[TABLE]

With this generalization it appears more natural to define a distance as the modulus of the complex number, that is the Euclidean distance eq.(20). However one should give some arguments to motivate the definition eq.(21), for example an interpretation of the complex conjugate $\lambda_{k}^{c*}$.

As we are speaking of distance between states, one may wonder if there is a relation with the distance between Hilbert states in quantum mechanics. With the interpretation above sketched, the measure eq.(19) is, in some sense, a statistical distance between the sequences, close to the distance between quantum states defined in Wootters [1981],

[TABLE]

where $p_{z}$ is the probability for a state $x$ between the states 1 and 2 and $X$ is the total number of different intermediate states. In our case, $X$ should be read as the number of different sequences in the path from the first to the second sequence, going step by step, each step corresponding to a transition or to a transversion. In our case $p_{x}$ is equal $\forall\;x$.

Let us discuss, with a few simple examples, the difference between the Hamming ($d_{H}$) and Kimura ($d_{K}$) measure, respectively with parameter $a$, parameter $\alpha$ (transition) and $\beta$ (transversion), and the above introduced one eq.(19), which we call crystal basis measure $d_{CB}$

This is a very simple didactic case. Let us compute the distance between the following 4-strings $1\equiv CAUC,2\equiv CUCG,3\equiv AUGU,4\equiv CUGA,5\equiv CUGU,6\equiv CAGU$ we have

[TABLE]

where the sequences are labelled as:

[TABLE]

The above defined distances are very simple, but they show the drawbacks to be explicitly dependent from somehow arbitrary parameters and implicitly from the length of the sequence. It is convenient do define normalized distance in the following way:

1. (a)

the standard Hamming distance, $d^{st}_{H}$, the percentage of different symbols, that is the number of different symbols divided by the sequence length

[TABLE] 2. (b)

the standard 2-parameters Kimutra distance $d^{st}_{K}$ Kimura [1980]

[TABLE]

where $p$ ($q$) is the percentage of positions which are transitionally (resp. transversionally) different, that is ($n_{tns}$ and $n_{tnv}$ denoting resp. the number of transitions and transversion)

[TABLE] 3. (c)

we can define the normalized crystal basis distance as it follows

[TABLE]

Eq.(28) is known in the literature as the “Canberra distance”. Note the eq.(28) is not rigorously a distance as it does not satisfy the triangle assiom. In order to get a normalized rigorous distance one should use the normalized form of the Euclidean distance eq.(20).

Eq.(23), using the above defined measures, becomes

[TABLE] 2. 2.

Let us consider the 6 sequences of Table 1 as obtained from the ancestral sequence $\,CCCCC$, by two point mutations $C\rightarrow U$. The Hamming and Kimura distance from the ancestral sequence is, for all the 6 sequences,

[TABLE]

The crystal basis measure between the ancestral sequence $d_{CB,i}$ (labelled by $(k=0,\ldots,4)\,\lambda_{k,H}\equiv\mid{\textstyle{\frac{5}{2}}};{\textstyle{\frac{5}{2}}};2;{\textstyle{\frac{3}{2}}};1\rangle$) and the i-th sequence (the sequences being ordered as in Table 1) is almost different for any sequence and can be computed as

[TABLE]

From eq.(31) we read that the most distant sequence from the ancestral one is the one where the mutations $C\rightarrow U$ have occurred in the 2nd and 4th position while the most close is the one where the nucleotides have mutated in the first two positions. This example shows clearly that the distance $d_{CB}$ is not only depending from the number of the substitutions in the nucleotides in the sequences, but from the positions where the substitutions have occurred, a peculiar feature which may be of biological interest. 3. 3.

Let us compute the distance between the $N=92$ coding sequence of $\beta-globin$ gene for Homo sapiens, Didelphis virginiana (Opossum) and Gallus gallus (see web site

http://www.hindawi.com/journals/tswj/2012/104269/tab2/ for the nucleotide sequences and Figg. 6-6 for the values of the $J_{H}$ and $J_{V}$).

From Table 2 we remark that the most distant (the most close) $\beta-globin$ sequences, according to our parametrization, are those between Gallus and Opossum (resp. between Homo and Gallus) in agreement with estimates of Yu [2011] and, partially, with those of Todeschini and al. [2008]. Differently the values of the Hamming and Kimura distance are not at all reliable.

Let us point out once more that the value of $d_{H}(1,2)$ ($d_{K}(1,2)$) depends only on the number of different nucleotides (resp. on the number of transitions and transversions) between the sequence 1 and 2, not from their position, while $d_{CB}(H,G)$ depends from the number of different nucleotides, from the nature of the nucleotides as well as from their position; $n_{1}$ transitions and $n_{2}$ transversions in different positions will in general lead to a different value of $d_{CB}$. Indeed a mutation of nucleotide in the m-th position changes the values of $J^{m}_{i}$ and of $J_{3,i}$ and may change the value of $J^{k}_{i}$,  $k>m$.

Remarks:

• •

The $d_{CB}$ distance is naturally defined for sequences of the same length. It may be extended to a distance between two sequences of different length by some suitable scaling procedure, but care must be paid for preserving the nature of integer or half-integer numbers of the labels $\lambda_{i}$.

• •

As $\lambda_{0,i}=J_{3,i}$ labels different states in the same irrep., one may introduce a different parameter for this label, that is to replace eq.(19) by

[TABLE]

However this choice is irrelevant if one makes use of the normalized distance.

3.2 Discrete Fourier Transform of a sequence

The DFT of a numerical sequence is extremely useful because it allows to reveal periodicities in the input data as well as the relative strength of the periodic components and is defined by777See Wikipedia, https://en.wikipedia.org/wiki/Discrete_Fourier_transform.

[TABLE]

The inverse DFT is obtained as

[TABLE]

For real $\lambda_{n}$, which is our case, the above equation takes the simpler form

[TABLE]

where

[TABLE]

Eq.(35) can be rewritten, N even, as

[TABLE]

and, N odd, as

[TABLE]

For real $\lambda_{n}$ it can be useful to rewrite eq.(33) in the following form, which helps to understand the features of the Figg. 8-10

[TABLE]

The first term of eq.(37)-(38) ($Re\,\widehat{\lambda_{0}}=\widehat{\lambda_{0}}=\sum_{k=0}^{k=N-1}\lambda_{k}$, $Im\,\widehat{\lambda_{0}}=0$) can be read as the average value of a constant real component, the terms in the sum are sine waves of frequency k/N, amplitude $2|\widehat{\lambda_{k}}|/N$ and phase shift $\Phi_{k}$ and the last term is the (real) coefficient of the highest frequency (= 1/2), not present for N odd.

The challenge to define Fourier transforms for DNA sequences, more generally for strings of symbolic data, has already been faced in the literature, one of the main aim being the detection of periodicities in DNA or protein sequences. The simplest solution is to map each symbol to a number, but, without an underlying algebraic framework, in general this choice is rather arbitrary and the result depends from the numerical labelling of the symbols, the nucleotides in the DNA/RNA case. Therefore several distinct and unrelated approaches exist. At our knowledge the first paper on this subject was published in 1986 Silverman and Linsker [1986] assigning to each letter of the DNA alphabet a vertex of a regular tetrahedron. Since then, many different DFT methods have been proposed and extensively used to study periodicities and repetitive elements in DNA sequences, genomes and protein structures. Indeed a few years later in S.Tavare and Giddings [1989] a numerical sequence, called the indicator sequence, of numerical values [math] and $1$ has been associated to a symbolic sequence and in Coward [1997] the two methods have been compared. Other approaches have been proposed and in Afreixo et al. [2004] the authors discuss the equivalence and connections between several different and distinct approaches. In general the DFT are different for different approaches, even if sometimes they share common features. Moreover in Yin et al. [2014] an euclidean distance between the power spectra of the DFT of two sequences, even of different length, has been proposed and applied in particular for hierarchical clustering. However all these approaches are based on the assignment of a numerical indicator, different for any nucleotide, to the nucleotide in the k-th position.

In the present paper, making explicit use of the underlying algebraic structure definite in the crystal basis model, the DNA symbolic sequence is transformed in a numerical sequence, build up from the whole sequence composition, on which the DFT can be defined in the standard way. Once introduced the DFT one may apply all the tools of this transform. The DNA/RNA sequence is seen as a signal, which can be decomposed, through the DFT, into different “frequency” components. Various filtering processing can be carried out as needed for the specific application, making more easily the manipulation of the coefficients for different frequency components. When we try to determine the “frequency content”, we are trying to decompose a complicated signal into simpler parts. Many physical signals are best described as a sum of many individual frequencies components instead of the time domain description. The job of the DFT is exactly to determine which frequencies a complicated signal is composed of. This translates for DNA/RNA signal to a description in terms of the “frequency content” rather than a description in terms of the nucleotides content. This type of analysis may reveal a powerful tool, for example, in the analysis of genes, in the comparison of coding and not coding sequences and for detecting periodicities.

Of course the DFT can be directly performed on the the sequence of N nucleotides using the labels of eq.(1), but in this case the set of half-integer number contains a limited amount of information, for example the Parseval theorem ($\lambda_{k}\in\{\pm 1/2\}$)

[TABLE]

for any sequence of fixed length gives the same value, irrespective of the arrangement of the nucleotides in the sequence.

We report in Fig.2 and Fig.4 the real, respectively the imaginary part, of the DFT of the labels of sequence eq.(15). Let us note that the symmetry, which shows in the two figures around the point 6, neglecting the point 0, is just a consequence of definition eq.(33), that is for $\lambda_{k}$ real the DFT shows the circular hermitian symmetry, i.e. $\widehat{\lambda_{n,a}}=\widehat{\lambda_{N-n,a}}^{*}$, where ∗ denotes the complex conjugate.

In Fig.8 and Fig.8 (Fig.10 and Fig.10) we draw the real and the imaginary part of the DFT of the H labels (resp. of V labels) of human $\beta-globin$ gene. We do not draw $\widehat{\lambda_{0,H}}=512,5$ and $\widehat{\lambda_{0,V}}=384,5$. In the figures we remark that, due to the circular hermitian symmetry, the real part (the imaginary part) of DFT coefficient is symmetric (resp. anti-symmetric) with respect to the label 46. This property reduces the analysis of the DFT coefficients of the labels of very long sequences of a one-half factor. Moreover we remark the presence of a peak for $k=2$ ($|\widehat{\lambda_{2,H}}|\approx 23$, $|\widehat{\lambda_{2,V}}|\approx 89$), corresponding to a frequency $k/N=1/46$, and a secondary peak for the case of V labels for $k=1$ ($|\widehat{\lambda_{1,V}}|\approx 57$), corresponding to a frequency $k/N=1/92$. Due to the circular symmetry these peak appear also for $k\prime=N-k$ with shift of equal sign for H labels (of opposite sign for V labels).

The edge effects which appear in Figg. 8-10 are due to the fact that the $\lambda_{k}$ in eq.(33) are real, with a specific pattern.

The real part of the DFT of a motif of the same nucleotide ($NNN..N$,$N\in\{C,G,U/T,A\}$), e.g. an isochore motif, is an horizontal segment, while the imaginary part is an increasing straight segment.

4 Conclusions

We propose a mathematical modeling in which a DNA/RNA sequence is represented by a numerical sequence, which somehow reflects the distribution of nucleotides in the original sequence.

It is commonly believed that all the information for protein folding is contained in the coding RNA sequences. May the crystal basis parametrization of a sequence be useful for the investigation of this point ?

From the proposed parametrization, it follows that a metrics can be defined between DNA/RNA sequences which may be useful to compute, e.g., which mutations push away more distantly a mutated sequence from the ancestral one.

The description allows, e.g., to look for symmetry of subsequences from a completely new point of view and in Subsection 2.2 some remarks and comments on this point are given.

Moreover, as a sequence takes the form of a signal, it can be decomposed into its elementary components performing the DFT of the signal. The “frequency content” of a sequence should be better interpreted to look for a possible deeper biological meaning than its mathematical definition. It seems that there is a dominant frequency $\nu=2/N$, depending on the length of the sequence.

Let us emphasize once more that the aim of this paper is not to make specific applications of this mathematical description for DNA or RNA sequences, but to try to illustrate the method and to enlighten, without any claim of completeness, its potentiality and versatility. Dealing with specific problems may possibly validate some of the possibilities here outlined.

The algebraic structure, underlying the model, illustrated in Subsection 2.1, which provides the framework for the mathematical parametrization, has been exploited, in the present paper, only to represent a nucleotide sequence as a state in the N-fold tensor product of fundamental irreps.. The wealth of the algebraic structure can be further used to model a mutation or deletion or insertion of a nucleotide in a sequence as the action of suitable tensor operator. Let us recall that tensor operators $\tau_{m}^{j}$ of $\mathcal{U}_{q\to 0}(sl(2))$, with the consequent extension of the Wigner-Eckart theorem, have been defined in Marotta and Sciarrino [2000]. An explanation of the organization of the genetic code in this framework has been given in Sciarrino [2003] and a preliminary application to deletion from 4-nucleotide (3-nucleotide) to 3-nucleotide (resp. 2-nucleotide) has been discussed Frappat et al. [2001].

The ($q\to 0$)-Wigner-Eckart theorem has the peculiar feature that the selection rules, that is the transitions between states which are permitted or forbidden, do not depend only on the rank of the tensor operator and on the initial state, but in a crucial way from the specific component of the tensor in consideration.

For example, the modeling of a pure negative trannsversion ($\Delta J_{H}=0,\Delta J_{V}=-1$, i.e. $C\rightarrow G$ or $T\rightarrow A$) as the tensor product of the identity operator on the H spin and a vector operator on the V spin (${\bf 1}_{H}\otimes\tau_{-1}^{1_{V}}$) acts on the state of human $\beta-globin$ as it follows (we indicate only the values of the labels $J_{3,i},J^{N}_{i},i=H,V$)

[TABLE]

and, in the considered case, does not give any restriction on the mutation. On the contrary the modeling of a pure positive trannsversion ($\Delta J_{H}=0,\Delta J_{V}=1$, i.e. $G\rightarrow C$ or $A\rightarrow T$) as the tensor product of the identity operator on the H spin and a vector operator on the V spin (${\bf 1}_{H}\otimes\tau_{1}^{1_{V}}$) acts on the state of human $\beta-globin$ as it follows

[TABLE]

and it implies that the mutation can only occur in a nucleotide in the $n-th$ position, $n\leq 81$. This is a relatively weak bound, further constraints arise by requiring the tensor operator to behave in a fixed way for the irreps. labelled by $J^{K}_{i}$, $K<N$, that is considering a kind of “russian doll tensor operator”.

We are hopeful that the presently proposed analysis may prove to be an useful tool for elucidating the topology of genetic regulatory networks, for putting into evidence possible existence of symmetry, for pointing out possible DNA “fundamental frequencies”, for analyzing exon and intron regions, cis-regulatory regions and so on.

Let us repeat that the approach here presented has to be further developed and, above all, its usefulness has to be validated by application to specific problem.

Acknowledgments It is a pleasure to thank Paul Sorba for useful discussions and suggestions.