On the Construction of Elastic Degenerate Strings

Computing an optimal partition of an interval $[1..c]$, that is, a segmentation, is a problem that arises naturally when dealing with the transformation of sequences that are aligned into $c$ positions, after trivial and contradicting constraints are imposed on the desired features of the solution. Figure 1 gives an example of such segmentation on a multiple sequence alignment (MSA) of $6$ sequences. In an MSA, the $r$ input sequences are made to be of equal length $c$ by adding gaps “$-$” forming a matrix of $r$ rows and $c$ columns which we will denote $\mathrm{𝖬𝖲𝖠}[1..r,1..c]$. The quality of such MSA can be measured in many ways, but already the problem of placing gaps such that the number of columns containing a single symbol is maximized equals the problem of finding the length of the longest common subsequence (LCS) of the input sequences, which is NP-hard [19]. Nevertheless, MSAs play a vital role in bioinformatics, and practical heuristics have been developed to obtain high quality MSAs [23]. In this paper, we assume an MSA to be given as input, and we look for rigorous ways to convert the MSA into a sequence of sets of strings, that is, into an elastic degenerate string (EDS) [18] (also known as elastic-degenerate text). In Figure 1, we can observe that a segmentation of the MSA directly induces an EDS by interpreting the strings of each segment to be part of a set.

The transformation from an MSA to an EDS through segmentation can be seen as lossy compression: we merge together equal strings inside each segment and we ignore the exact connections between segments. Such lossy compression can be a desired feature, for example, for hiding the sensitive full sequence input by publishing only the EDS. A more rudimentary reason for such conversion is that EDSes and related pangenome graphs do compress well and enable fast processing and analysis of variation in populations [9]. A notable feature of this scheme is that if a query string occurs in a row of the MSA, then it will also occur in the EDS, but not the other way around. Various ways to optimize the transformation to minimize false positives while not sacrificing privacy could be considered. Minimizing a specific parameter is useful when the downstream analyses of the resulting EDS employ algorithms whose complexity depends from that parameter.

This paper is initiating the study of EDS optimization by looking at some selected optimization criteria (Section 3) for which we can modify algorithms from related literature. We give linear-time algorithms (Section 4) to minimize the total cardinality of the EDS sets (cardinality) and total length of the strings in the EDS sets (size) under specific assumptions: the algorithms directly work for gapless MSAs, and they serve as heuristics when gaps are interpreted as part of the alphabet. In the latter case, the optimization criteria should be adjusted accordingly: see Section 5 for the multiple problems that we leave open for general MSAs. We start by reviewing the related work in Section 2.

Founder reconstruction.

Ukkonen [25] studied the problem of explaining a given set of $r$ aligned sequences representing $c$ haplotype sites (thus no insertions or deletions) with the recombination of a few founder sequences, under the simple crossover model of recombination (sequence switching at aligned positions); minimizing the founder set is trivial if the recombination can happen at any position, but if the recombinations are allowed only at selected positions, then the resulting segmentation problems – minimizing the maximum cardinality of a segment (the number of founders in [25]) given a lower bound on the segment length, or given a target average segment length size to surpass – can be solved in polynomial time.

Norri et al. [22] solved one of the problems posed by Ukkonen – the segmentation of aligned sequences minimizing the maximum height given a segment-length lower bound $L$ – in linear $O(rc)$ time and $O(r+L)$ space. The solution proceeds in a dynamic programming fashion, from left to right, and for each aligned position $y$ it uses the positional Burrows–Wheeler Transform [10], enhanced by a few other arrays, to efficiently compute all the possible heights of a segment ending at $y$ and their optimal recursive values.

Cazaux et al. [8] studied the related problem of maximizing the minimum segment length or the average segment length of any segmentation that has height of at most a given $H$, motivated by the minimization of the size of the resulting pangenomic index based on founder sequences. To do so, they recognize that the approach by Norri et al. [22] fits into a more general left-to-right column-stream model of computation, and use range maximum queries to efficiently compute the recursive values. For this last task, they modify existing techniques for the semi-dynamic setting where an array is indexed for amortized constant-time range queries in an online fashion.

Algorithms on EDSes.

There is a wide literature on pattern search and comparison on (elastic) degenerate strings [18, 17, 5, 6, 1, 21, 4, 2, 3] where algorithmic results exhibit complexity that depend on both the cardinality and the size of input EDSes. Moreover, the latest algorithms for EDS-based pangenome comparison (intersection, matching statistics, similarity and distance measures) also depend on the cardinality and size metrics [15, 14, 16]. These results give further motivation for our focus on optimizing these measures.

Indexable Founder Graphs.

Elastic Founder Graphs (EFGs), introduced by Equi et al. [13], take the simplified recombination model of founder sequences (recombination at selected positions only) to transform an MSA of $r$ sequences and $c$ aligned positions – with insertions and deletions and thus gaps – into an acyclic Elastic Block Graph: the segmentation results into consecutive blocks of nodes. In the general case, the resulting graph is still hard to index just like more general graphs [12, 11], so Equi et al. [13] prove that if each chosen segment spells strings that exclusively occur from the segments’ starting column, then this property is conserved in the resulting graph, which is easy to index. This results in a framework of linearithmic (i.e. $O(rc\log (rc))$) time constructions of indexable EFGs, which support fast pattern matching.

Rizzo et al. [24] improved some construction solutions to linear time (i.e. $O(rc)$) while studying the problem of minimizing the maximum segment length or maximum height of a segmentation resulting in an indexable EFG.

MSA covers.

While segmentation induces a limited class of acyclic pangenome graph representations, there is an analogous approach by Cartes et al. [7] to cover MSAs and generate arbitrary acyclic graphs; the related optimization problems are much harder, although small instances can be solved using integer linear programming.

We denote integer interval $\{1,\mathrm{\dots },n\}$ as $[1..n]$. Let $\mathrm{𝖬𝖲𝖠}[1..r,1..c]\in {(\mathrm{\Sigma }\cup \{-\})}^{r\times c}$ be a multiple sequence alignment of $r$ rows and $c$ columns representing the input sequences (or strings) and the aligned positions, respectively. It is built from a finite integer alphabet $\mathrm{\Sigma }=[1..\sigma ]$ ($\mathrm{\Sigma }=\{𝙰,𝙲,𝙶,𝚃\}$ in our examples) augmented with the gap symbol $-\notin \mathrm{\Sigma }$ to represent insertions and deletions. It is immediate to see that the alignment takes at most $O(rc)$ words of space, or $O(rc\log |\mathrm{\Sigma }|)$ bits of space. We denote the $i$-th row of the MSA as $\mathrm{𝖬𝖲𝖠}[i,1..c]$, and the concatenation of its characters from position $x$ to $y$, with $x\le y$, as $\mathrm{𝖬𝖲𝖠}[i,x..y]$. We say that $\mathrm{𝖬𝖲𝖠}[1..r,1..c]$ is gapless if no gap symbol is used. The operator $\mathrm{spell}(T)$ removes all the gap symbols from a given string $T\in {(\mathrm{\Sigma }\cup \{-\})}^c$: for example, $\mathrm{spell}(\mathrm{𝖬𝖲𝖠}[i,1..c])$ is the $i$-th (unaligned) sequence. If $T$ contains only gaps then $\mathrm{spell}(T)=\epsilon$, the empty string. Consider the following definition alongside Figure 1.

Consider the problem of finding a segmentation $S$ of minimum cardinality $\mathrm{m}(S)$. If the MSA contains few but very long sequences, that is, $r\ll c$, the trivial segmentation $S^{\equiv }=[1..c]$ might be optimal. We can encourage a certain level of recombination in the corresponding EDS by setting an upper bound $U$ on the length of the accepted segments.

On the other hand, consider finding a segmentation $S$ of minimum gap-aware size ${\mathrm{N}}_{\epsilon }(S)$. In Section 4.1 we show that the trivial segmentation $S^{\equiv }=[1],\mathrm{\dots },[c]$ has always minimum size in the gapless setting, resulting in a highly recombinant EDS. We can symmetrically discourage recombination in regions of low sequence similarity by fixing a lower bound $L$ on the length of the accepted segments.

One crucial data structure that we use in Section 4 is based on sorting iteratively the MSA rows while reading the alignment from left to right, column by column. Let $\mathrm{\Sigma }$ have an implicit total order $\le$ on its characters. Given strings $S,T\in {\mathrm{\Sigma }}^c$, we say that $S$ is lexicographically smaller than $T$, in symbols $S{\prec }_{\mathrm{lex}}T$, if $S[1..x]=T[1..x]$ and for some $x\in [0..c-1]$ (where $S[\mathrm{1..0}]$ is equal to the empty string $\epsilon$). Symmetrically, $S$ is co-lexicographically smaller than $T$, in symbols $S{\prec }_{\mathrm{colex}}T$, if $S[x..c]=T[x..c]$ and for some $x\in [1..c+1]$ (where $S[c+1..c]=\epsilon$). Moreover, we call $S[1..x]$ a prefix of $S$ and $S[1..x]$ a suffix of $S$. They are non-empty if $x\ge 1$ and $x\le c$, respectively.

Finally, one last data structure that we use in this paper indexes an integer array to later find out the minimum value inside a given range of array positions.

In Section 4.1 we remark useful properties of gapless MSAs. Then, in Sections 4.2 and 4.3 we provide linear-time solutions to min-$U$-cardinality and min-$L$-size for MSAs under this setting.

This initial study on EDS construction opens many interesting directions for further study:

Gaps as symbols yield upper bounds

Even though Equations 2 and 3 (but not Equation 4) hold for MSAs with gaps, Theorems 13 and 14 work only on MSAs without gaps. A straightforward way to extend them is to treat them as normal symbols in the alphabet. After segmentation, the gap symbols can be removed to obtain the final EDS. With all the quality measures we proposed, this works as an upper bound: for example, consider the height of a segment that contains A- and -A; height is reduced by one after gaps are removed. Thus, the removal of gaps cannot increase the height, and it is of interest to study experimentally how much the measures change after gap removal.

Gaps as symbols yield exact solutions?

For height optimization, it seems possible to characterize the class of MSAs with gaps where the approach above yields an optimal solution, and not just an upper bound. Such characterization should be possible through forbidding local sub-optimal alignments like A- and -A. This raises several natural questions: What is the exact form of the characterization? How to efficiently detect if a given MSA is part of this class? Do all optimal alignments under some column-wise scoring function belong to this class? Are there efficient algorithms to convert an MSA into one in this class without sacrificing the alignment score?

Tailored algorithms with gaps

The presented Algorithms 1 and 2 do not find the optimal segmentation of arbitrary MSAs, so it is natural to ask if there are efficient algorithms to directly optimize the studied measures on general MSAs. Alternatively, there may be slight adjustments to the measures that are amenable to efficient algorithms: in a different context, Rizzo et al. [24] introduced the metric of prefix-aware height, which was optimized in linear time, and also optimized in linear time the metric of maximum segment length.