DNA Sequences as Numerical Sequences

Formally, a DNA sequence ${\displaystyle s[i]}$ is an ordered list of symbols from a dictionary of nucleotides ${\displaystyle (A,C,T,G)}$. There are multiple ways to view this as a numerical sequence, and they will all be applied to the same example:

${\displaystyle s[i]=\{A,T,G,C,A,G,C\}}$

• Binary Projection

Define a binary signal for each nucleotide, ${\displaystyle U_{\alpha }[i]}$, which is 1 when the i-th position in the sequence is the nucleotide ${\displaystyle \alpha }$ and 0 otherwise. Formally,

${\displaystyle u_{\alpha }[i]={\begin{array}{cc}\{&{\begin{array}{cc}1&s[i]\,{\text{is}}\,\alpha \\0&{\text{o/w}}\end{array}}\end{array}}}$

This creates four signals which encode the position of the four nucleotides in the sequence. For the above example, the projected signals would be

ATGCAGC, ${\displaystyle u_{a}}$ = 1000100, ${\displaystyle u_{c}}$ = 0001001, ${\displaystyle u_{t}}$ = 0100000, ${\displaystyle u_{g}}$ = 0010010

• Correlation Function

Let ${\displaystyle E(i,j)}$ be a function that is 1 if ${\displaystyle s[i]=s[j]}$ and ${\displaystyle 0}$ otherwise. Let ${\displaystyle c[t]=\sum _{i=t}^{N}E(i,i-t)}$. In other words, it counts the number of times that the nucleotide at a position is equal to the nucleotide t-away on the sequence. It is invariant to changing the labeling of the bases. ${\displaystyle c[0]}$ will always be equal to the length of the sequence. For the above example, the projected signals would be

${\displaystyle c(0)=7}$, ${\displaystyle c(1)=0}$, ${\displaystyle c(2)=0}$, ${\displaystyle c(3)=2}$, ${\displaystyle c(4)=1}$, ${\displaystyle c(5)=0}$, ${\displaystyle c(6)=0}$

• Linear Map

For mathematical purposes, a gene sequence could be viewed as a signal by mapping each nucleotide into the range [1,4]. For the example above and map ${\displaystyle \{A->1,T->2,C->3,G->4\}}$, the sequence would manifest itself as ${\displaystyle 1243143}$.

This process generates a single signal for the sequence, but raises questions such as which out of 24 (4!) maps should be used and what effect does this map have on our analysis. Unlike the other two mapping methods, this one is not invariant to changing the labeling of the bases but having the same structure (i.e. AACA -> TTGT), which could be detrimental for some applications.

Spectral Analysis of DNA Sequences

Once the DNA sequence has been converted into a numerical sequence, spectral analysis can be performed on that sequence. Recall the DFT is defined by the analysis equation

${\displaystyle X[k]=\sum _{n=0}^{N-1}x[n]\cdot e^{-2\pi ikn/N},}$

and produces a N-long complex signal ${\displaystyle X[k]}$ that represents the spectrum of the signal. The index ${\displaystyle k}$ of the frequency response corresponds to an angular frequency ${\displaystyle \omega }$ of ${\displaystyle 2\pi k/N}$. For the 3-base periodicity property, the frequency content at ${\displaystyle k=N/3}$ should be analyzed as it corresponds to a period of 3.

Recall that the power spectra of a sequence is equal to the magnitude of the frequency vector squared,

${\displaystyle P[k]=|X[k]|^{2}}$

Tiwari applied the DFT to analysis of DNA sequences using the binary projection operator.^[6] They calculated the spectra of the four projected nucleotide sequences using the DFT. Call them ${\displaystyle C_{A}[k],C_{C}[k],C_{T}[k],C_{G}[k]}$. Then, the power spectra of all the signals were combined,

${\displaystyle C[k]=|C_{A}[k]|^{2}+|C_{C}[k]|^{2}+|C_{T}[k]|^{2}+|C_{G}[k]|^{2}.}$

Then, they calculated the Signal-to-Noise ratio of this signal and performed a threshold test on that value to determine whether or not that stretch of DNA is coding.

Anti-Notch Filter

Instead of computing a Discrete Fourier Transform on different segments of the signal, the analysis can be performed in the time-domain through the use of an anti-notch filter at frequency ${\displaystyle N/3}$, which is much faster and performs a similar operation. An anti-notch filter applied to the signals essentially approximates the magnitude of the DFT at ${\displaystyle N/3}$, and a similar method is used to threshold this value to determine whether or not the window is coding or non-coding. Vaidyanathan and Yoon implemented and evaluated this method.^[8] An ideal anti-notch filter at frequency ${\displaystyle N/3}$ has an impulse response of

${\displaystyle h[n]=e^{2\pi nj/3}}$

and thus the frequency content of the output of an anti-notch filtered signal is ${\displaystyle X[k]=A[N/3]}$ if ${\displaystyle k=N/3}$ and 0 otherwise. By taking the magnitude of the time-domain signal, and invoking Parseval's Theorem, we get the magnitude of the frequency response. By the above logic, this is just ${\displaystyle |A[N/3]|^{2}}$ which is the quantity of interest.

Spectrograms

Spectrograms are a good way to view how the frequency content of a signal changes over time. The most common way to compute a spectrogram is to compute a Fourier transform over different segments of the signals, convert the frequency magnitude plot into an image, and concatenate those images. This is a useful way to visually identify coding and non-coding regions of DNA and to inspect other patterns that might exist.^[9]

As a Feature in Gene Identification

Identifying genes in a DNA sequence is harder than just finding what segments are coding and near impossible to identify by visually inspecting spectrograms. Genes are made up of both coding and non-coding regions, called introns and exons. Thus, the transition between coding and non-coding regions must be examined and analyzed properly to identify genes. Computing the "level" of 3-periodicity over different (possibly overlapping) windows of the sequence generates a plot of 3-periodicity over time.^[10]

These long stretches of coding vs. non-coding can then be classified as introns or exons and the entire segment heuristically labeled as gene or non-gene.

Why Can it Discriminate Coding from Non-Coding?

Binary signals can be parsed into something called a position count function (PCF), which counts the number of one's at phase ${\displaystyle s}$ at positions that are multiples of ${\displaystyle \omega }$. This is expressed mathematically for a binary signal ${\displaystyle A}$ as

${\displaystyle C_{w}^{A}(s)=\sum _{i=0}^{\frac {N-1}{\omega }}A[\omega i+s].}$

For ${\displaystyle \omega =3}$, there are three PCFs ${\displaystyle (C_{3}^{A}(0)}$, ${\displaystyle C_{3}^{A}(1)}$ and ${\displaystyle C_{3}^{A}(2))}$ which each count the number of one's at positions ${\displaystyle \{0,3,6,...\}}$, ${\displaystyle \{1,4,7,...\}}$, ${\displaystyle \{2,5,8,...\}}$ respectively. It can be shown that the magnitude of the DFT at ${\displaystyle N/3}$ is in fact directly related to the PCF and is equal to

${\displaystyle |{\tilde {A}}[N/3]|^{2}={\frac {1}{2}}[(C_{3}^{A}(0)-C_{3}^{A}(1))^{2}+(C_{3}^{A}(1)-C_{3}^{A}(2))^{2}+(C_{3}^{A}(2)-C_{3}^{A}(0))^{2}].}$ ^[11]

In other words, the spectral power at ${\displaystyle N/3}$ in the DFT depends only on the difference between the PCFs at different phases. This relation here is the fundamental to the explanation for why this method can identify coding regions. This is because the spectral power at ${\displaystyle N/3}$ will be determined by how the codons that make up the sequence are sampled.

If the codons are sampled uniformly at random, as they would be in a noncoding segment of DNA, there is a high chance that the PCFs would not differ by a significant amount and the power at that frequency will be low.

However, in a protein-coding sequence, the DNA sequence is made up of a string of codons which correspond to amino acids. Because the genetic code is degenerate (more than one codon map to a single amino acid) and samples from the amino acids rather than the codons, the codons are not sampled uniformly thus leading to differences in the PCFs.

There is also empirical evidence for why this method works. In other words, over multiple studies this method has been able to discriminate coding vs. non-coding DNA segments. These are discussed in the next section.

Pros

• Requires no training data

The method can be run on any DNA sequence, where as other methods such as BLAST, FASTA and Smith-Waterman require empirical data.

• Independent of variations in base compositions

This is because the total spectral power is

${\displaystyle |S[k]|^{2}=|A[k]|^{2}+|T[k]^{2}+|C[k]|^{2}+|G[k]|^{2},}$

and no one base contributes more than another.

• Low false-positive rate
• Robust to sequencing errors resulting in frameshifts

This is due to the property that shifting a sequence does not change the magnitude of its Discrete Fourier Transform.

${\displaystyle x[n-n_{0}]\leftrightarrow X[k]e^{-j2\pi kn_{0}/N}}$

${\displaystyle |X[k]e^{-j2\pi kn_{0}/N}|=|X[k]|}$

Cons

• Three-base periodicity found to be lacking in some genes
• Outperformed by modern empirical and specialized gene prediction methods