Specific Patterns Against Reference Sequences

The goal of this article is to design an alignment-free technique for comparing a set of sequences or words, called a target, against a set of words, called a reference.

The motivation comes from the analysis of genomic sequences as done, for example, by Khorsand et al. in [23] in which the authors introduce the notion of sample-specific strings. To avoid alignments but to extract interesting elements that differentiate the target from the reference, or in general two words, the chosen specific fragments are minimal absent words, also called minimal forbidden words. Target-specific words are factors of the target that are minimal absent words of the reference.

These types of factors associated with absent patterns are rather commonly used to efficiently compare sequences by avoiding complete alignments of full sequences, see, for example, [11] and references therein. In bioinformatics, target-specific words serve as signatures for newly sequenced biological molecules, helping to identify their characteristics. In the domain of molecular biology they allow the discovery of remarkable patterns in some genomic sequences, such as persitent patterns in the analysis of SARS-CoV-2 genomes that are absent in the human genome [31] and minimal sequences in Ebola virus also absent in the human genome [33]), or to build phylogenies of biological molecular sequences using a distance based on absent words (see [10, 32]). These patterns are also helpful to improve both pattern matching methods (see [15]) and text compression with the concept of antidictionaries (see for example [17, 1]). They are also central in some approaches related to the sanitisation of data, that is, the process of hiding confidential information [5, 6], to quote only a few applications.

The notion of a minimal absent word was introduced by Mignosi et al. [27] (see also [4]) in relation to combinatorial aspects of certain sequences. The first linear-time computation of the minimal absent words of the factors of a single sequence is described in [16] (see also [14]). Its time complexity is $O(n)$ on a fixed-size alphabet for a sequence of size $n$. The algorithm uses the computation of the directed acyclic word graph (DAWG), also called suffix automaton (see [8, 14]), of a single sequence. The same time complexity holds for an integer alphabet of polynomial size (see [19], [20]). This result is obtained using the $O(n)$ time complexity for computing the DAWG of a single sequence, in the case of an integer alphabet of polynomial size. For a general ordered alphabet $A$, the running time becomes $O(n\log |A|)$.

These algorithms extend to the computation of the minimal absent words of the set of factors of a finite set of sequences of total size $n$. This is done in [3] in $O(n)$ for a fixed-size alphabet. It becomes $O(n\log |A|)$ for a general ordered alphabet. It is done in $O(n)$ in [29] for an integer alphabet of polynomial size, extending the linear time computation of the DAWG in [19] and [20] to a finite set of sequences. It is also mentioned in [2] that the DAWG of a finite set of sequences can be computed in $O(n)$ time using sparse matrices for an integer alphabet of polynomial size.

Due to the significant role of the notion, the efficient computation of minimal absent words has attracted considerable attention (see, for example, [30] and references therein).

In this article, we continue exploring the approach of target-specific words as in [23] by introducing new algorithmic techniques to detect them. A more general view of the usefulness of formal languages in analyzing a series of genomes using pangenomic graphs is described by Bonizzoni et al. in [9].

A preliminary version of this paper was presented at the DLT 2023 conference [2]. In the present article, the motivation is strengthened due to additional references, the presentation of algorithms has been improved and complete proofs have been added. In addition, Section 5 reveals a surprising phenomenon on the tight link between DAWG matching and the search for minimal absent words, a technique at the core of our solution in Section 4.

In this section, we recall the definition and sketch the construction of the directed acyclic word graph of a finite set of words. This description already appears in [3].

Let $P=\{x_1,x_2,\mathrm{\dots },x_r\}$ be a finite set of $r$ finite words. A linear-time construction of a deterministic finite state automaton recognizing $\mathrm{Fact}(P)$ has been obtained by Blumer et al. in [8, 7] (see also [28]). Their construction is an extension of the well-known incremental construction of the suffix automaton of a single word (see for instance [12, 14]). The words are added one by one to the automaton. In the sequel, we call this algorithm the Dawg algorithm since it outputs a deterministic automaton called a directed acyclic word graph. Let us denote it by $\text{DAWG}(P)=(Q,A,i,Q,\delta )$. Let $\text{Suff}(v)$ denote the set of suffixes of a word $v$ and $\text{Suff}(P)$ the union of all $\text{Suff}(v)$ for $v\in P$.

The states of $\text{DAWG}(P)$ are the equivalence classes of the right invariant equivalence ${\equiv }_{\text{Suff}(P)}$ defined as follows. For $u,v\in \mathrm{Fact}(P)$,

There is a transition labelled by $a$ from the class of a word $u$ to the class of $ua$. The automaton $\text{DAWG}(P)$ has a unique initial state, which is the class of the empty word, and all its states are final. Note that the (syntactic) congruence $\sim$ defining the minimal automaton of the language is

and is not the same as the below equivalence. In other words, $\text{DAWG}(P)$ is not always a minimal automaton.

The construction of $\text{DAWG}(P)$ is performed in time $O(\mathrm{size}(P)\times \log |A|)$ in [8]. A time complexity of $O(\mathrm{size}(P))$ can be obtained for an integer alphabet of polynomial size in [19] and [20]. It can also be obtained with an implementation of automata with sparse matrices (see [14, Exercise 1.15]).

An essential element of the efficient construction is the notion of suffix links between states, denoted by $s$. We first define the function $s^{\prime }$ from $\mathrm{Fact}(P)\setminus \{\epsilon \}$ to $\mathrm{Fact}(P)$ as follows: for any $v\in \mathrm{Fact}(P)\setminus \{\epsilon \}$,

Then, we define the partial function $s$ from $Q$ to $Q$ as follows. When $p=\delta (i,v)$ for $v\ne \epsilon$, $s(p)$ is the state $\delta (i,s^{\prime }(v))$. The function $s$ is not defined on the initial state $i$.

In this section, we assume that the reference $R$ and the target $T$ are two finite sets of words and our goal is to compute the set of T-specific factors of $T$ against $R$. To do so, we first compute the directed acyclic word graph $\text{DAWG}(R\cup T)=(Q,A,i,Q,\delta )$ of $R\cup T$. Further, we compute a table, called $\text{mark}$, indexed by the set of states $Q$ and that satisfies: for each state $p$ in $Q$, $\text{mark}[p]$ is one of the three values $r$, $t$ or both $r,t$ according to the fact that each word labeling a path from $i$ to $p$ is a factor of some word in $R$, in $T$ or in both. This information can be obtained during the construction of the directed acyclic word graph without increasing the running time.

The following algorithm outputs the set of $T$-specific words occurring in $T$ with respect to $R$ in the form of a trie (digital tree, see [14]).

In this section, we consider that $T$ is a single word and $R$ is a finite set of words as before. The goal of the section is to design an algorithm that computes all the occurrences in $T$ of words that are $T$-specific with respect to $R$. To do so, we define the $T$-specific table associated with the pair $R,T$ of words of the problem.

A letter of $T$ at position $k$ is denoted by $T[k]$ and $T[i..j]$ denotes the factor $T[i]T[i+1]\mathrm{\cdots }T[j]$ of $T$. Then, the T-specific table $\mathrm{𝑇𝑠}$ is defined, for $i=0,\mathrm{\dots },|T|-1$, by

Note that the set of $T$-specific factors is prefix-free, that is, no proper prefix of an element of the set is also in the set. (The set of $T$-specific factors is also suffix-free.) Therefore, for each position $k$ on $T$, there is at most one $T$-specific factor of $T$ starting at $k$ (and for each position $j$ on $T$ there is at most one $T$-specific factor of $T$ ending at $j$).

Instead of computing the T-specific table $\mathrm{𝑇𝑠}$, in a straightforward way, the algorithm below can be transformed to compute the list of pairs $(i,j)$ of positions on $T$ for which $\mathrm{𝑇𝑠}[i]=j$ and $j\ne -1$.

To compute the table we use $ℛ$, the Suffix automaton or rather the DAWG of $R$. The former is the minimal automaton accepting the suffixes of $R$ (see [14, Section[5.4]) and the latter has the same structure but with all states as terminal states instead (see [8]). As such, it accepts $\mathrm{Fact}(R)$ the set of factors of $R$.

The automaton is given with its transition function $\delta$, its initial state $i$ and and is equipped with both the suffix link $s$ (used here as a failure link) and the length function $\mathrm{\ell }$ defined on states. The function $\mathrm{\ell }$ is defined by: $\mathrm{\ell }[p]=\max \{|z|\in A^{\ast }\mid \delta (i,z)=p\}$. The functions $s$ and $\mathrm{\ell }$ transform the automaton into a search machine (see [14, Section 6.6]).

Figure 3 illustrates the principle of Algorithm TsTable. Let us assume that the factor $u=T[k..j-1]$ is a factor of $R$ but $ub$ is not for some letter $b$. Then, let $v$ be the longest suffix of $u$ for which $vb$ is a factor of $R$. If it exists, then clearly $avb$, with the letter $a$ preceding $v$, is $T$-specific. Indeed, $av,vb\in \mathrm{Fact}(R)$ and $avb\notin \mathrm{Fact}(R)$, which means that $avb$ is a minimal absent word of $R$ while occurring in $T$. Therefore, setting $q^{\prime }=\delta (i,v)$ and since $|v|=\mathrm{\ell }[q^{\prime }]$, the minimal absent word $avb$ is identified by setting $\mathrm{𝑇𝑠}[j-\mathrm{\ell }[q^{\prime }]-1]=j$. If there is no suffix of $u$ satisfying the condition, the letter $b$ alone is $T$-specific and is identified by setting $\mathrm{𝑇𝑠}[j]=j$.

Searching for absent factor occurrences as done in the previous section is based on the DAWG matching technique [13] (see also [14, Section 6.6]) in the domain of string matching algorithms. As such, Algorithm TsTable looks like a side product of string matching.

The goal of string matching algorithms is to locate a given pattern $p$ in a longer text $T$. Figure 3 displays the generic situation when a part $u$ of the pattern has been found in $T$ and is about to be appended with the letter $b$. The extension of $u$ is successful if $b$ matches its aligned letter of the pattern, which eventually can lead to the detection of an occurrence of the whole pattern.

The situation in Figure 3 is also the generic situation when searching for minimal absent words. However, in contrast, this is an unsuccessful match of the letter $b$ that immediately yields the detection in $T$ of a minimal absent word of $p$.

The role of the DAWG matching technique is essential here to detect minimal absent words of $p$ occurring in $T$ because the DAWG of $p$ stores and provides direct access to all factors of $p$. Instead, if an online string matching algorithm is used, like KMP [24] or Simon-Hancart [34, 22] algorithms (see also [14, Chapter 2]), the algorithm will detect only some absent words. That is, those of the form $aub$ in which $au$ is a factor of $p$ but $ub$ is only a prefix of it. This does not produce all the minimal absent words. However, other types of string matching based on text indexes can certainly be used for the same purpose.

As a conclusion, the algorithm designed in the previous section reveals a surprising phenomenon: the very tight link between the existence of minimal absent words that pop up naturally in the analysis of pattern searching based on a Suffix automaton or a DAWG.