On avoided words, absent words, and their application to biological sequence analysis

Introduction

The one-to-one mapping of a DNA molecule to a sequence of letters suggests that DNA analysis can be modelled within the framework of formal language theory [1]. For example, a region within a DNA sequence can be considered as a “word” on a fixed-sized alphabet in which some of its natural aspects can be described by means of certain types of automata or grammars. However, a linguistic analysis of the DNA needs to take into account many distinctive physical and biological characteristics of such sequences: The genome consists of coding regions that encode for polypeptide chains associated with biological functions as well as a plethora of regulatory and potentially functional non-coding regions, identified through multiple alignment of genomes of several organisms, and termed conserved non-coding elements (CNEs). In addition, it contains large non-coding regions most of which are not linked to any particular function. All these genomic components appear to have many statistical features in common with natural languages [2].

A computational tool oriented towards the systematic search for avoided words is particularly useful for in silico genomic research analyses. The search for absent words is already undertaken in the recent past and several results exist on the application and computation of such words [3–6]. However, words which may be present in a genome or in genomic sequences of a specific role (e.g., protein coding segments, regulatory elements, conserved non-coding elements etc.) but they are strongly underrepresented—as we can estimate on the basis of the frequency of occurrence of their longest proper factors—may be of particular importance. They can be words of nucleotides which are hardly tolerated because they negatively influence the stability of the chromatin or, more generally, the functional genomic conformation; they can represent targets of restriction endonucleases which may be found in bacterial and viral genomes; or, more generally, they may be short genomic regions whose presence in wide parts of the genome are not tolerated for less known reasons. The understanding of such avoidances is becoming an interesting line of research (for recent studies, see [7, 8]).

On the other hand, short words of nucleotides may be systematically avoided in large genomic regions or whole genomes for entirely different reasons, i.e. just because they play important signaling roles which confine their appearance only in specific positions: consensus sequences for the initiation of gene transcription and of DNA replication are well-known such oligonucleotides. Other such cases may be insulators, sequences anchoring the chromatin on the nuclear envelope like lamina-associated domains, short sequences like dinucleotide repeat motifs with enhancer activity, and several other cases. Again, we cannot exclude that this area of research could lead to the identification of short sequences of regulatory activities still unknown.

Brendel et al. in [9] initiated research into the linguistics of nucleotide sequences that focuses on the concept of words in continuous languages—languages devoid of blanks—and introduced an operational definition of words. The authors suggested a method to measure, for each possible word w of length k, the deviation of its observed frequency from the expected frequency in a given sequence. The values of the deviation, denoted by \(\textit{dev}(w)\), were then used to identify words that are avoided among all possible words of length k. The typical length of avoided (or of overabundant) words of the nucleotide language was found to range from 3 to 5 (tri- to pentamers). The statistical significance of the avoided words was shown to reflect their biological importance. This work, however, was based on the very limited sequence data available at the time: only DNA sequences from two viral and one bacterial genomes were considered. Also note that k might change when considering eukaryotic genomes, the complex dynamics and function of which might impose a more demanding analysis. The authors in [10–12] have studied the concept of unusual words—based on different definitions than the ones Brendel et al. use for expectation and variance—focusing on the factors of a string, whereas based on Brendel et al. definitions, we consider here any word over the alphabet.

Our contributions

The computational problem can be described as follows. Given a sequence x of length n, an integer k, and a real number \(\rho < 0\), compute the set of \(\rho \)-avoided words of length k, i.e. all words w of length k for which \(\textit{dev}(w) \le \rho \). We call this set the \(\rho \)-avoided words of length k in x. Brendel et al. did not provide an efficient solution for this computation [9]. Notice that such a word may be completely absent from x. Hence the set of \(\rho \)-avoided words can be naïvely computed by considering all possible \(\sigma ^k\) words, where \(\sigma \) is the size of the alphabet.

Here we present an \(\mathcal {O}(n)\)-time and \(\mathcal {O}(n)\)-space algorithm for computing all \(\rho \)-avoided words of length k in a sequence of length n over a fixed-sized alphabet. For words over an integer alphabet of size \(\sigma \), the algorithm requires time \(\mathcal {O}(\sigma n)\), which is optimal for sufficiently large \(\sigma \). We also present a time-optimal \(\mathcal {O}(\sigma n)\)-time algorithm to compute all \(\rho \)-avoided words (of any length) in a sequence of length n over an integer alphabet of size \(\sigma \). We provide a tight asymptotic upper bound for the number of \(\rho \)-avoided words over an integer alphabet and the expected length of the longest one. We also prove that the same asymptotic upper bound is tight for the number of \(\rho \)-avoided words of fixed length when the alphabet is sufficiently large.

As shown subsequently, the set of absent \(\rho \)-avoided words is a subset of the set of minimal absent words of a word. Hence the tight asymptotic bounds for \(\rho \)-avoided words are based on the proof we provide for the tightness of the known asymptotic bound on minimal absent words and the tightness of this bound for minimal absent words of fixed length over sufficiently large alphabets.

We make available an open-source implementation of our algorithm. Experimental results, using both real and synthetic data, show its efficiency and applicability. Specifically, using our method we confirm that restriction endonucleases which target self-complementary sites are not found in eukaryotic sequences [8]. In addition, we apply our algorithm in the case of CNEs, which are classes of sequences whose functions in our genomes remain largely enigmatic [13, 14]. We observe interesting patterns of occurring avoided words within CNEs compared to CNE-like sequences (surrogates) that are in accordance with their distinct sequence characteristics which classify them from other non-functional sequences [15, 16].

A preliminary version of this article has appeared in [17].

Terminology and technical background

Definitions and notation

We begin with basic definitions and notation generally following [18]. Let \(x=x[0]x[1] \cdots x[n-1]\) be a word of length \(n=|x|\) over a finite ordered alphabet \(\Sigma \) of fixed size \(\sigma \), i.e. \(\sigma = |\Sigma |=\mathcal {O}(1)\). We also consider the case of an integer alphabet; in this case each letter is replaced by its rank such that the resulting string consists of integers in the range \(\{1,\ldots ,n\}\). For two positions i and j on x, we denote by \(x[i \ldots j]=x[i]\cdots x[j]\) the factor (sometimes called subword) of x that starts at position i and ends at position j (it is empty if \(j < i\)), and by \(\varepsilon \) the empty word, word of length 0. We recall that a prefix of x is a factor that starts at position 0 (\(x[0\ldots j]\)) and a suffix is a factor that ends at position \(n-1\) (\(x[i \ldots n-1]\)), and that a factor of x is a proper factor if it is not x itself. A factor of x that is neither a prefix nor a suffix of x is called an \(\textit{infix}\) of x. We say that x is a power of a word y if there exists a positive integer k, \(k>1\), such that x is expressed as k consecutive concatenations of y; we denote that by \(x=y^k\).

Let \(w=w[0]w[1] \cdots w[m-1]\) be a word, \(0<m\le n\). We say that there exists an occurrence of w in x, or, more simply, that w occurs in x, when w is a factor of x. Every occurrence of w can be characterised by a starting position in x. Thus we say that w occurs at the starting position i in x when \(w=x[i \ldots i + m - 1]\). Further let f(w) denote the observed frequency, that is, the number of occurrences of a non-empty word w in word x. Note that overlapping occurrences are considered as distinct ones; e.g. \(f(\texttt {TT})=2\) in \(\texttt {TTT}\). If \(f(w) = 0\) for some word w, then w is called absent, otherwise, w is called occurring.

By \(f(w_p)\), \(f(w_s)\), and \(f(w_i)\) we denote the observed frequency of the longest proper prefix \(w_p\), suffix \(w_s\), and infix \(w_i\) of w in x, respectively. We can now define the expected frequency of word w, \(|w|>2\), in x as in Brendel et al. [9]:

$$\begin{aligned} E(w) = \frac{f(w_p) \times f(w_s)}{f(w_i)},\quad \text { if~ } f(w_i) >0; \text {~else~} E(w) = 0. \end{aligned}$$

(1)

The above definition can be explained intuitively as follows. Suppose we are given \(f(w_p)\), \(f(w_s)\), and \(f(w_i)\). Given an occurrence of \(w_i\) in x, the probability of it being preceded by w[0] is \(\frac{f(w_p)}{f(w_i)}\) as w[0] precedes exactly \(f(w_p)\) of the \(f(w_i)\) occurrences of \(w_i\). Similarly, this occurrence of \(w_i\) is also an occurrence of \(w_s\) with probability \(\frac{f(w_s)}{f(w_i)}\). Although these two events are not always independent, the product \(\frac{f(w_p)}{f(w_i)} \times \frac{f(w_s)}{f(w_i)}\) gives a good approximation of the probability that an occurrence of \(w_i\) at position j implies an occurrence of w at position \(j-1\). It can be seen then that by multiplying this product by the number of occurrences of \(w_i\) we get the above formula for the expected frequency of w.

Moreover, to measure the deviation of the observed frequency of a word w from its expected frequency in x, we define the deviation (\(\chi ^2\) test) of w as:

$$\begin{aligned} \textit{dev}(w) = \frac{f(w)-E(w)}{\max \{ \sqrt{E(w)}, 1\}}. \end{aligned}$$

(2)

For more details on the biological justification of these definitions see  [9].

Using the above definitions and a given threshold, we are in a position to classify a word w as either avoided or common in x. In particular, for a given threshold \(\rho < 0\), a word w is called \(\rho \)-avoided if \(\textit{dev}(w) \le \rho \). In this article, we consider the following computational problems.

Suffix trees

In our algorithms, suffix trees are used extensively as computational tools. For a general introduction to suffix trees, see [18].

The suffix tree \(\mathcal {T}(x)\) of a non-empty word x of length n is a compact trie representing all suffixes of x. The nodes of the trie which become nodes of the suffix tree are called explicit nodes, while the other nodes are called implicit. Each edge of the suffix tree can be viewed as an upward maximal path of implicit nodes starting with an explicit node. Moreover, each node belongs to a unique path of that kind. Then, each node of the trie can be represented in the suffix tree by the edge it belongs to and an index within the corresponding path.

We use \(\mathcal {L}(v)\) to denote the path-label of a node v, i.e., the concatenation of the edge labels along the path from the root to v. We say that v is path-labelled \(\mathcal {L}(v)\). Additionally, \(\mathcal {D}(v)= |\mathcal {L}(v)|\) is used to denote the word-depth of node v. Node v is a terminal node, if and only if, \(\mathcal {L}(v) = x[i \ldots n-1]\), \(0 \le i < n\); here v is also labelled with index i. It should be clear that each occurring word w in x is uniquely represented by either an explicit or an implicit node of \(\mathcal {T}(x)\). The suffix-link of a node v with path-label \(\mathcal {L}(v)= \alpha y\) is a pointer to the node path-labelled y, where \(\alpha \in \Sigma \) is a single letter and y is a word. The suffix-link of v exists if v is a non-root internal node of \(\mathcal {T}(x)\). We denote by Child \((v,\alpha )\) the explicit node that is obtained from v by traversing the outgoing edge whose label starts with \(\alpha \in \Sigma \).

In any standard implementation of the suffix tree, we assume that each node of the suffix tree is able to access its parent. Note that once \(\mathcal {T}(x)\) is constructed, it can be traversed in a depth-first manner to compute the word-depth \(\mathcal {D}(v)\) for each node v. Let u be the parent of v. Then the word-depth \(\mathcal {D}(v)\) is computed by adding \(\mathcal {D}(u)\) to the length of the label of edge (u, v). If v is the root then \(\mathcal {D}(v) = 0\). Additionally, a depth-first traversal of \(\mathcal {T}(x)\) allows us to count, for each node v, the number of terminal nodes in the subtree rooted at v, denoted by \(\mathcal {C}(v)\), as follows. When internal node v is visited, \(\mathcal {C}(v)\) is computed by adding up \(\mathcal {C}(u)\) of all the nodes u, such that u is a child of v, and then \(\mathcal {C}(v)\) is incremented by 1 if v itself is a terminal node. If a node v is a leaf then \(\mathcal {C}(v) = 1\).

Example 1

Consider the word \(x=\texttt {AGCGCGACGTCTGTGT}\). Fig. 1 represents the suffix tree \(\mathcal {T}(x)\). Note that word \(\texttt {GCG}\) is represented by the explicit internal node v; whereas word \(\texttt {TCT}\) is represented by the implicit node along the edge connecting the node labelled 15 and the node labelled 9. Consider node v in \(\mathcal {T}(x)\); we have that \(\mathcal {L}(v) = \texttt {GCG}\), \(\mathcal {D}(v) = 3\), and \(\mathcal {C}(v)=2\).

Tight bounds on minimal absent words

We first show that the known asymptotic upper bound on the number of minimal absent words of a word is tight.

In the following lemma we show that, for sufficiently large alphabets, \(\mathcal {O}(\sigma n)\) is a tight asymptotic bound for the number of minimal absent words of fixed length.

Useful properties of avoided words

In this section, we provide some useful insights of combinatorial nature which were not considered by Brendel et al. [9]. By the definition of \(\rho \)-avoided words it follows that a word w may be \(\rho \)-avoided even if it is absent from x. In other words, \(\textit{dev}(w) \le \rho \) may hold for either \(f(w) > 0\) (occurring) or \(f(w) = 0\) (absent).

This means that a naïve computation should consider all possible \(\sigma ^k\) words. Then for each possible word w, the value of \(\textit{dev}(w)\) can be computed via pattern matching on the suffix tree of x. In particular, we can search for the occurrences of w, \(w_p\), \(w_s\), and \(w_i\) in x in time \(\mathcal {O}(k)\) [18]. In order to avoid this inefficient computation, we exploit the following crucial lemmas.

Avoided words algorithm

In this section, we present Algorithm AvoidedWords for computing all \(\rho \)-avoided words of length k in a given word x. The algorithm builds the suffix tree \(\mathcal {T}(x)\) for word x, and then prepares \(\mathcal {T}(x)\) to allow constant-time observed frequency queries. This is mainly achieved by counting the terminal nodes in the subtree rooted at node v for every node v of \(\mathcal {T}(x)\). Additionally during this pre-processing, the algorithm computes the word-depth of v for every node v of \(\mathcal {T}(x)\). By Lemma 3, \(\rho \)-avoided words are classified as either occurring or (minimal) absent, therefore Algorithm AvoidedWords calls Routines AbsentAvoidedWords and OccurringAvoidedWords to compute both classes of \(\rho \)-avoided words in x. The outline of Algorithm AvoidedWords is as follows.

Computing absent avoided words

In Lemma 3, we showed that each absent \(\rho \)-avoided word is a minimal absent word. Thus, Routine AbsentAvoidedWords starts by computing all minimal absent words in x; this can be done in time and space \(\mathcal {O}(n)\) for a fixed-sized alphabet or in time \(\mathcal {O}(\sigma n)\) for integer alphabets [4, 5]. Let \(< (i,j), \alpha>\) be a tuple representing a minimal absent word in x, where for some minimal absent word w of length \(|w| > 2\), \(w = x[i \ldots j]\alpha \), \(\alpha \in \Sigma \); this representation is clearly unique.

Intuitively, the idea is to check the length of every minimal absent word. If a tuple \(< (i,j), \alpha>\) represents a minimal absent word w of length \(k = j-i+2\), then the value of \(\textit{dev}(w)\) is computed to determine whether w is an absent \(\rho \)-avoided word. Note that, if \(w = x[i \ldots j]\alpha \) is a minimal absent word, then \(w_p= x[i \ldots j]\), \(w_i= x[i+1 \ldots j]\), and \(w_s = x[i+1 \ldots j]\alpha \) occur in x by Definition 1. Thus, there are three (implicit or explicit) nodes in \(\mathcal {T}(x)\) path-labelled \(w_p\), \(w_i\), and \(w_s\), respectively.

The observed frequencies of \(w_p\), \(w_i\), and \(w_s\) are already computed during the pre-processing of \(\mathcal {T}(x)\). For an explicit node v of \(\mathcal {T}(x)\), path-labelled \(w'= x[i' \ldots j']\), the value \(\mathcal {C}(v)\), which is the number of terminal nodes in the subtree rooted at v, is equal to the number of occurrences (observed frequency) of \(w'\) in x. For an implicit node along the edge (u, v) path-labelled \(w''\), the number of occurrences of \(w''\) is equal to \(\mathcal {C}(v)\) (and not \(\mathcal {C}(u)\)). The implementation of this procedure is given in Routine AbsentAvoidedWords.

Computing occurring avoided words

Lemma 4 suggests that for each occurring \(\rho \)-avoided word w, \(w_p\) is a path-label of an explicit node v of \(\mathcal {T}(x)\). Thus, for each internal node v such that \(\mathcal {D}(v)= k-1\) and \(\mathcal {L}(v)= w_p\), Routine OccurringAvoidedWords computes \(\textit{dev}(w)\), where \(w =w_p \alpha \), \(\alpha \in \Sigma \), is a path-label of a child (explicit or implicit) node of v. Note that if \(w_p\) is a path-label of an explicit node v then \(w_i\) is a path-label of an explicit node u of \(\mathcal {T}(x)\); node u is well-defined and it is the node pointed at by the suffix-link of v. The implementation of this procedure is given in Routine OccurringAvoidedWords.

Optimal computation of all ρ-avoided words

Although the biological motivation is yet to be shown for this, we present here how we can modify Algorithm AvoidedWords so that it computes all \(\rho \)-avoided words (of all lengths) in a given word x of length n over an integer alphabet of size \(\sigma \) in time \(\mathcal {O}(\sigma n)\). We further show that this algorithm (AllAvoidedWords) is in fact time-optimal.

Based on Lemma 1 and similarly to the proof of Lemma 5 we obtain the following result.

It is clear that if we just remove the condition on the length of each minimal absent word in Line 2 of AbsentAvoidedWords we then compute all absent \(\rho \)-avoided words in time \(\mathcal {O}(\sigma n)\). In order to compute all occurring \(\rho \)-avoided words in x it suffices by Lemma 4 to investigate the children of explicit nodes. We can thus traverse the suffix tree \(\mathcal {T}(x)\) and for each explicit internal node, check for all of its children (explicit or implicit) whether their path-label is a \(\rho \)-avoided word. We can do this in \(\mathcal {O}(1)\) time as described. The total number of these children is at most \(2n-1\), as this is the bound on the number of edges of \(\mathcal {T}(x)\) [18]. This modified algorithm is clearly time-optimal for fixed-sized alphabets as it then runs in time \(\mathcal {O}(n)\). The time optimality for integer alphabets follows directly from Lemma 8. Hence we obtain the second result of this article.

Experiment I

To evaluate the time performance of our implementation, synthetic DNA (\(\sigma =4\)) and protein (\(\sigma =20\)) data were used. The input sequences were generated using a randomised script. In the first experiment, our task was to establish that the performance of the program does not essentially depend on k and \(\rho \); i.e., the elapsed time of the program remains unchanged up to some constant with increasing values of k and decreasing values of \(\rho \). As input datasets, for this experiment, we used a DNA and a protein sequence both of length 1M (1 Million letters). For each sequence we used different values of k and \(\rho \). The results, for elapsed time are plotted in Fig. 2. It becomes evident from the results that the time performance of the program remains unchanged up to some constant. The longer time required for the protein sequences for some value of k is explained by the increased number of branching nodes in this depth in the corresponding suffix tree due to the size of the alphabet (\(\sigma =20\)). To confirm this we counted the number of nodes considered by the algorithm to compute the \(\rho \)-avoided words for \(k=4\) and \(\rho =-10\) for both sequences. The number of considered nodes for the DNA sequence was 260 whereas for the protein sequence it was 1,585,510.

Experiment II

In the second experiment, our task was to establish the fact that the elapsed time and memory usage of the program grow linearly with n, the length of the input sequence. As input datasets, for this experiment, we used synthetic DNA and proteins sequences ranging from 1 to 128 M. For each sequence we used constant values for k and \(\rho \): \(k=8\) and \(\rho =-10\). The results, for elapsed time and peak memory usage, are plotted in Fig. 3. It becomes evident from the results that the elapsed time and memory usage of the program grow linearly with n. The longer time required for the protein sequences compared to the DNA sequences for increasing n is explained by the increased number of branching nodes in this depth (\(k=8\)) in the corresponding suffix tree due to the size of the alphabet (\(\sigma =20\)). To confirm this we counted the number of nodes considered by the algorithm to compute the \(\rho \)-avoided words for \(n=64\)M for both the DNA and the protein sequence. The number of nodes for the DNA sequence was 69,392 whereas for the protein sequence it was 43,423,082.

Experiment III

In the next experiment, our task was to evaluate the time and memory performance of our implementation with real data. As input datasets, for this experiment, we used all chromosomes of the human genome. Their lengths range from around 46M (chromosome 21) to around 249M (chromosome 1). For each sequence we used \(k=8\) and \(\rho =-10\). The results, for elapsed time and peak memory usage, are plotted in Fig. 4. The results with real data confirm that the elapsed time and memory usage of the program grow linearly with n.

Experiment IV

Experiment V

Then, we proceeded to the examination of several collections of CNEs obtained through multiple sequence alignment between the human and other genomes. The detailed description of how those CNEs were identified could be found in [15]. For each CNE of these datasets, a sequence stretch (surrogate sequence) of non-coding DNA of equal length and equal GC content was taken at random from the repeat-masked human genome. The CNEs of each collection were concatenated into a single long sequence and the same procedure was followed for the corresponding surrogates. Seven CNEs concatenates and the corresponding surrogate datasets have been formed and used in this experiment. We have determined through the proposed algorithm the avoided words for \(k=10\) (decamers) and \(\rho =-2\) for these fourteen datasets and the results are presented in Table 1. In Table 2, we show likewise for \(k>2\) (all avoided words) and \(\rho =-2\).

The first five CNEs collections have been composed through multiple sequence alignment of the same set of genomes and they differ only in the thresholds of sequence similarity applied between the considered genomes: from 75 to 80 (the least conserved CNEs, which thus are expected to serve less demanding functional roles) to 95–100 which represent the extremely conserved non-coding elements (UCNEs or CNEs 95–100) [15]. The remaining two collections have been composed under different constraints and have been derived after alignment of genomes belonging to the Mammalian and Amniotic groups. In Tables 1 and 2, the last line shows the ratios formed by the numbers of avoided words of each concatenate of surrogates divided by the numbers of avoided words of the corresponding CNE dataset.

Our second observation, i.e. the positive correlation between (1) the paucity of avoided decamers in CNEs collections and (2) the similarity thresholds used for their identification comes in accordance with the above argument. CNEs extracted under a stricter requirement of sequence similarity between evolutionary distant species are CNEs whose functionality is less tolerant to alterations due to random mutations in general. Hence, they also tolerate less the propagation within their sequence of parasite polypyrimidine/polypurine tracts too.