Received 23 June 2025, accepted 12 July 2025, date of publication 11 August 2025, date of current version 8 September 2025.

Digital Object Identifier 10.1109/ACCESS.2025.3597547

RESEARCH ARTICLE

# Dynamic Programming Alignments With Skips

NJAGI MOSES MWANIKI ID AND NADIA PISANTI ID

Department of Computer Science, University of Pisa, 56127 Pisa, Italy

Corresponding author: Nadia Pisanti (nadia.pisanti@unipi.it)

This work was supported in part by European Union's Horizon 2020 Research and Innovation Program under the Marie Skłodowska-Curie Grant through the Algorithms for Pangenome Computational Analysis (ALPACA) Project under Grant 956229 and through the Pangenome Graph Algorithms and data Integration (PANGAIA) Project under Grant 872539, in part by Miistero Universita' e Ricerca (MUR) Progetto di Ricerca di Interesse Nazionale (PRIN) through the Project Pangenome Informatics: from Theory to Applications (PINC) under Grant 2022YRB97K, and in part by the NextGeneration European Union (EU) Program Tuscany Health Ecosystem under Grant PNRR ECS00000017.

**ABSTRACT** The outcome of a Multiple Sequence Alignment (MSA) can be compactly represented by means of Elastic-Degenerate Strings (ED-strings) by collapsing conserved fragments into standard linear strings, while representing gaps and variants as sets of alternative strings. These alternative *variants* can differ in size and can possibly include the empty string. In 2022, Lee et al. introduced Partial Order Alignment (POA) to enable the alignment of a string against a graph-like structure derived from an MSA. However, the POA edit transcript (the sequence of edit operations that describe the alignment) does not reflect the possible elasticity of the MSA (such as different gaps sizes in the aligned string), leaving room for possible misalignment and its propagation in progressive MSA strategies. In this paper, we propose a dynamic programming based method that optimally aligns a string to an ED-string, the latter compactly representing an MSA, overcoming the ambiguity in the POA edit transcript while maintaining its time and space complexity. Moreover, since pangenomes can also be represented using ED-strings, our algorithm paves the way to a new class of sequence to graph alignment methods capable of taking into account possible gaps in the pangenome representation, thus offering a richer and more flexible model for pangenomic analysis.

**INDEX TERMS** Computational pangenomic, dynamic programming, ED string, partial order alignment, sequence alignment.

## I. INTRODUCTION

The analysis of genomic data is currently undergoing a paradigm shift: from relying on a single genome sequence as a reference - the reference genome - to using a set of aligned related sequences jointly represented in a compact structure: the reference pangenome.<fnref n="1" /> As opposed to a linear reference genome, a pangenome reference allows a simultaneous representation of variations and commonalities among the underlying sequences, to be used as a baseline for variants analyses of individuals as well as populations, without incurring in a single reference genome bias related to the choice of the specific individual genome.

In particular, for species like humans that enjoy a widespread availability of sequencing data as well as a

The associate editor coordinating the review of this manuscript and approving it for publication was Yilun Shang.

growing awareness of genomic variants, many in silico biomedical analysis tasks using pangenomes can heavily improve the accuracy of significant variants detection [2].

This promising opportunity also introduces significant computational challenges: choosing a representation that effectively balances expressive power and computational efficiency for structures that have to be stored, searched, updated, and investigated [3]. This is one of the key reasons why the aforementioned paradigm shift has led to the emergence of the new *-omic* discipline *pangenomics*. More specifically, *computational pangenomics* addresses new algorithmic tasks where the traditional linear genome sequence model is replaced by more complex graph-like structures [2], [4]. As a result, recent years have seen a growing attention of the scientific bioinformatics community towards graph-based representations that can model a pangenome. The most general of these representations are edge- and/or node-labeled directed graphs [3], such as the *variation graphs* [5],

<footnote n="1">¹A *pangenome* is “any collection of genomic sequences to be analyzed jointly or to be used as a reference” [1].</footnote>
[6] (with their haplotype aware version of [7]) and sequence graphs [8]. Although all such representation are conceptually related, sharing the feature of modeling variants as alternative paths in the graph, and common/conserved regions as nodes or edges that can be reached by distinct paths, there are significant differences in the constraints of the graph topology as well as in what kind of sequence a single node or edge can represent. As a consequence, there is currently no universally accepted standard for a “(pan)genome graph” representation. Some methods are designed to operate across several graph types, while others are tailored specifically, or are particularly efficient, for a subset of them.

Simpler acyclic alternatives to the graph-based representations discussed above are *Elastic Founder graphs* [9], [10] (EF graphs), and *Elastic-Degenerate strings* [11], [12], [13] (ED strings). These build upon earlier non-elastic versions *Founder graphs* [14] (F graphs) and *Degenerate strings* [15], [16], [17] (D strings), respectively. Both ED strings and EF graph are sort of *linearized* versions of pangenome graphs in that they are both directed and acyclic (they do not allow back edges), and they both explicitly list local variants. As opposed to ED strings, EF graphs are slightly more haplotype-aware (that is, they partially track information of which of the original sequences the variant is from), at the cost of storing - for each pair of adjacent variants - explicit edges from one variant to the next, and also at the cost of checking such edges when traversing the graph. On the other hand, ED strings include the empty string among the possible variants, whereas EF graphs do not. As a consequence, with ED strings one can explicitly represent insertions and deletions of genomic fragments (a.k.a. *INDELs*), that are one of the most common type of variants.

Degenerate Strings were first introduced to deal with a common issue in experimental data [18] or, as it happens with genomic data coming from sequencing, when one needs to represent symbols whose identity is *unknown* or *uncertain*. In computational biology, for example, the IUPAC notation [19] extends the standard {A, C, G, T} genome alphabet to represent ambiguous positions in a DNA sequence where multiple alternative bases or fragments can be observed. Similarly to how pangenomes represent variations, uncertain sequences can be represented by means of degenerate string (or *indeterminate string* as they had first been named when variants where bound to be made of single letters [20]). Unlike a standard string, which is a sequence of single letters, a degenerate string is a sequence of finite *sets of letters*. When every set contains exactly one letter, then we are in the special case of a *standard string* (the traditional linear deterministic string).

Due to the challenges of pangenomic data analysis, in recent years there has been a growing interest in algorithmic tasks that involve the various structures, suggesting new tools that have good time and accuracy performance for the different graphs types, even if they do not always guarantee exact solutions (see [6], [21], [22], [23] just to

mention a few; for a comprehensive survey one can read [24] or [3]). In particular, many algorithms and tools [12], [17], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], as well as lower bounds [35], [36], [37], [38] have been designed in the last few years for the problem of aligning a sequence to a graph. Efforts on comparing pangenomes (that is, graph-to-graph alignment or comparison) are also worth being mentioned, even if this is more sophisticated and less developed task [15], [16], [39], [40], [41]. The applications of aligning these graphs to sequences (or among themselves) are many and include variant calling, progressive multiple sequence alignments, pangenome construction, and more in general pangenomic data analysis.

Compared to more expressive representations among those mentioned above, ED strings have a few algorithmic advantages that are basically the consequence of not allowing back edges and loops on their graph structure. First of all, they support fast and simple on-line string matching [12], [42]. Moreover, they allow deterministic subquadratic string matching [37], [43]. Also, and very importantly for the purposes of this paper, while the above mentioned more general pangenomes representations are more suitable than ED-strings to represent structural variants and large repeats, they all lack an explicit representation of skips.

The notion of ED-strings and, over them, the exact (that is, no edit operations other than exact match are allowed) matching problem *Elastic-Degenerate String Matching* (EDSM) problem has attracted some attention in the combinatorial pattern matching community in recent years. Since its introduction in [11], a series of results have been published [12], [36], [37], [43], and a few tools produced [42], [44]. In [45], [46], [47], [48], and [49], the approximate version of the problem has also been addressed allowing non-zero Hamming or edit distance between the pattern and its occurrence inside the ED-string. In particular, the current state of the art for computing the edit distance between a pattern *P* of length *m* and an ED-string $\tilde{T}$ of size<fnref n="2" /> *N* takes time in $O(k \cdot m \cdot G + k \cdot N)$, where $G \in O(N)$ is the total number of strings in $\tilde{T}$ and *k* is the computed edit distance [47]. Moreover, it can be efficiently decided whether two ED strings share a string [39] (that is, a boolean answer to an *intersection query*). Finally, D strings support fast dynamic-programming-based alignment [17] for approximate matching with a solid string, and linear-time answer to the intersection query [15], [16].

Sequence comparison is a well studied problem in computer science, with a central role as a fundamental task for many applications in computational biology [50], and therefore in computational pangenomics. Given two or more sequences and a distance or similarity function, the

<footnote n="2">²In [45] and [46] the size of the ED-string $\tilde{T}$ is defined as the actual number of letters or empty strings that appear in $\tilde{T}$, whereas here (see Section II) sizes of degenerate letters are the number of variants times the length of the longest among them.</footnote>
task is to compare the sequences in order to infer their (dis)similarities. For pairwise comparison, a typical such function in bioinformatics enumerates and weights possible edit events that can turn a sequence into another, and hence seeks for the transcript of such events (the *edit transcript*) whose total weight minimizes (resp. maximizes) the distance (resp. similarity). These edit events are typically insertions, deletions, matches, and mismatches. An edit transcript is thus traditionally a sequence in the letters M, X, I, and D, where M denotes a match, X is a mismatch, I is an insertion, and D denotes a deletion<fnref n="3" />

$$T_1 \rightarrow \text{GAATTATAC}$$

$$T_2 \rightarrow \text{G-G-TA--C}$$

Insertions and deletions in sequence alignment are also generalized as *gaps* (when referred to the single letter insertion or deletion), or *INDELs* (when referred to the biological mutation event) that can involve whole fragments as well as single letters, raising issues on how, and how much, to account for their costs. One sort of obligated choice is to weight insertions as much as deletions in order to guarantee a symmetry in the distance or similarity measure between the two sequences. Besides this, for example, a gap cost purely proportional to its length can erroneously discourage long runs of gaps, and a low cost for gaps can fragment the alignment into too many short gaps that hardly correspond to as many mutation events that might have taken place.

For example, the two sequences $T_1$ and $T_2$ above could as well be aligned as follows

$$T_1 \rightarrow \text{GAATTATAC}$$

$$T_2 \rightarrow \text{GG---TAC}$$

leading to the transcript is MXDDDDMMM, which is just a permutation of the other one. If a gap cost is purely proportional to its length, then a run of 4 consecutive deletions costs the same as 4 scattered deletions, and hence the two alignments would be equally optimal. If, however, the alignment is supposed to highlight actual mutation events that have turned one sequence into the other, then the run of consecutive gaps better fulfills the task because in real DNA a single deletion of a short fragment from $T_1$ to $T_2$ (or an insertion from $T_1$ to $T_2$) is more likely than several single deletions so close to each other. This has motivated the definitions of gap penalty functions such as the *affine gap penalty* score, specifically introduced for genomic sequences [51], that we will use in this paper. Given the two parameters o and e, the affine gap penalty function assigns to g consecutive gaps a score equal to o + g * e: a fixed opening cost o plus as many as g times the extension cost

e. With this scoring scheme, the first alignment above would require 4 times the opening cost o, resulting in a worst scoring alignment with respect to the second one.

When one has to compare more than just two sequences, the alignment task turns into *Multiple Sequence Alignment* (MSA). In genomic and pangenomic data analyses, as well as in other applications, MSAs are typically performed in a progressive manner by iteratively adding a new linear sequence to a “preliminary” MSA. To this purpose, the MSA can be sort of linearized into a lossy representation of the actual similarities and differences of the aligned strings, and then the iterative - crucial - step of progressive MSA boils down to a pairwise alignment of the new string against such MSA. This alignment involves two strings of different type, one linear and a more complex one, the latter representing the MSA, and therefore possibly carrying a whole set of variants, or a gap, in a single location, rather than a single value. This breaks up the traditional above mentioned symmetry of pairwise alignments, therefore suggesting new challenges that include re-opening the issue of how to account for gaps. Indeed, whereas in traditional pairwise sequence alignments an insertion is just the dual of a deletion and their cost is the same, here the insertion can mean two events that are substantially different: wherever the MSA has already located, and already accounted, a gap, then this insertion cost should not be accounted again and again to each new string iteratively added in a progressive alignment; on the other hand, when the gap is an insertion that had not occurred yet in the MSA, or a deletion in the new string, then this is indeed a new event, and it has to be penalized as such [52], [53], and [54]. In this paper we will formalize the distinction of these two events (the already accounted gap, and the new one) the former being denoted as *skip*. To the best of our knowledge, ours is the first attempt to manage skip events in a sequence alignment algorithm.

Degenerate strings can also be used to encode the consensus of a population of DNA sequences in a gapless multiple sequence alignment. Nevertheless, since real multiple alignments do contain gaps, in order to also encode insertions and deletions that occur in MSAs, the notion of Elastic Degenerate String is again better suited: a sequence of finite sets of *strings* that may have different size and that may possibly be the empty string: a gap.

As a result, even for the purpose of compactly represent the outcome of a MSA by collapsing conserved fragments into linear strings, and that of highlighting gaps and variants as list of strings that can have different size and that can include the empty string, we suggest to use the notion of ED-String.

Consider, for example, for the following MSA of three closely-related sequences $T_1$, $T_2$, and $T_3$:

$$T_1 \rightarrow \text{TAGTCTGTCCCCAT}$$

$$T_2 \rightarrow \text{TAAT-TTTCCCCAT}$$

$$T_3 \rightarrow \text{TA---TTGCCCCAT}$$

<footnote n="3">³We here anticipate that this paper will introduce a new edit transcript event S representing a skip. For example, the *edit transcript of an optimal alignment of the two sequences* $T_1 = \text{GAATTATAC}$ and $T_2 = \text{GGTAC}$ is *MXDMMDDM* corresponding to the alignment.</footnote>
whose ED-string representation $\tilde{T}$ is as follows

$$ \tilde{T} = \underline{TA} \left\{ \begin{array}{c} GTC \\ AT \\ \epsilon \end{array} \right\} \underline{T} \left\{ \begin{array}{c} \underline{GT} \\ \underline{TT} \\ \underline{TG} \end{array} \right\} \underline{CCCCA}T $$

where $\epsilon$ denotes the empty string, and where those in curly brackets are called degenerate positions and the strings listed therein are the variants.

It is important to remark that the variants represented by an ED-String introduce a partial order in its letters (letters are ordered from left to right as usual, but as opposed to the total order of letters in a linear string, here letters that belong to distinct variants of the same degenerate positions are not comparable). In this case, its edit distance can be computed using Partial Order Alignment (POA) alignments [55], as shown in [17] for the simpler case of D-strings. However, POA alone does not yield the most informative edit transcript as the empty string is not properly managed. Consider for example the simple ED-string $\tilde{T} = A\{C\}T$ and the string $P = AT$: we want the edit transcript of the optimal alignment to exhibit the skip operation, and thus to be MSM and not just MM. In this way, should the MSA be used as an intermediate step in an iterative progressive alignment, when adding a possible new string there would be a record of an existing skip at that location. In this paper we solve this problem by adding a new edit event in the edit transcript: the skip operation denoted by $S$, and we suggest a dynamic programming algorithm to optimally align a linear string $P$ and an ED-String $\tilde{T}$ with (free) skips operations.

Our aim is to compute the correct edit distance with skips (the lowest number of operations needed to convert $P$ into one of the strings represented by $\tilde{T}$ with free cost of skips) and its optimal edit transcript that includes skips events.

For example, the optimal alignment of the pattern $P = TATGCCCCCTTT$ to the ED-string $\tilde{T}$ above involves the letters of $\tilde{T}$ shown with: (i) underlined letters being matches, (ii) a skip operation in that the $\epsilon$ represent an existing gap locus in the MSA, and hence a free insertion, also underlined (iii) the overlined $T$ is instead a new insertion in $P$ and therefore it has a cost, (iv) the boldface letter $A$ is a mismatch (at its cost), and (v) the last letter of $P$ has been deleted (at a cost). Assuming unitary cost for all edit operations and no cost for matches and skips, then the distance between $P$ and $\tilde{T}$ above is 3 (1 insertion, 1 mismatch, 1 deletion), and the edit transcript representing an optimal alignment is MMSMMIMMMMXMD.

Furthermore, we also extend the traditional dynamic programming alignment algorithm with affine gap penalty function to work with ED-strings and skips. Our new edit operation skip, together with (i) a suitably designed rewriting $\tilde{T}'$ of the ED-string $\tilde{T}$ that conserves its size and its edit distance with any linear string, and (ii) the recursive formula that we consequently introduce for the dynamic programming computation of the edit distance between $\tilde{T}'$ and a string $P$, lead to our main result: a dynamic programming algorithm that computes an optimal alignment of $\tilde{T}$ and $P$ in time and

space complexity in $O(m \cdot N)$ (where $m = |P|$ and $N$ is the size of $\tilde{T}$).

The algorithm we describe here is for a global alignment (the natural choice for the progressive MSA application mentioned above), but its extension to local and semiglobal alignment is straightforward. A semiglobal alignment could be used, for example, in computational pan-genomics, with the ED-string representing a pan-genome [1], [15], [16], [24] and the pattern being a read to be mapped therein: this is currently one of the most relevant problems in bioinformatics [6], [29], [30], [42], [44], [56].

In Section III we explain how we propose to manage skips using as an example the basic case of an elastic but not degenerate strings (which we define in Section II with other preliminary definitions), and in Section IV we show how to extend the dynamic programming framework to ED-strings. This work is an extended version of part of [31].

## II. PRELIMINARY DEFINITIONS

In this section we formally define some of the notions that have been mentioned in Section I.

A string $X$ is a sequence of elements on an alphabet $\Sigma$, where the alphabet $\Sigma$ is a non-empty finite set of letters of size $|\Sigma|$: we say that $X \in \Sigma^*$, and specifically that $X \in \Sigma^n$ if $n$ is the length of $X$ ($n = |X|$). We denote with $\epsilon$ the empty character (a character that does not change the string it is in) which we will name skip, and with $\epsilon$ the empty string: a string made up of empty characters only. A string $\Upsilon$ containing one or more empty characters is an elastic string (E-string) ($\Upsilon \in (\Sigma \cup \{\epsilon\})^*$) where $\epsilon$ is the skip. The length $|\Upsilon|$ of $\Upsilon$ is computed including its skips as well.

For any string $X$, we denote by $X[i, j]$ the substring of $X$ that starts at position $i$ and ends at position $j$. In particular, $X[0, j]$ is the prefix of $X$ that ends at position $j$ (it is the empty prefix if $j = 0$), and $X[i, |X|]$ is the suffix of $X$ that starts at position $i$ (being the empty suffix if $i = 0$).

**Definition 1:** An elastic degenerate string (ED-string) $\tilde{S} = \tilde{S}\{1\}\tilde{S}\{2\}...\tilde{S}\{n\}$ of length $n$ over an alphabet $\Sigma$ is a finite sequence of $n$ positions $\tilde{S}\{i\}$'s. Each position $\tilde{S}\{i\}$ is a finite non-empty set of $s_i$ strings (also called variants) that can possibly properly include the empty string. We denote with $\ell_i$ the length of the longest string in $\tilde{S}\{i\}$.

When $s_i = 1$, then $\ell_i = 1$ as well, and in that case $\tilde{S}\{i\}$ is just a simple letter of $\Sigma$ and in this case we say that $\tilde{S}\{i\}$ is solid, while if $s_i > 1$ then we say that it is a degenerate position. The following parameters measure the degree of degeneracy of an ED-string: the (total) size $N$ and (total) width $W$ of an ED-string $\tilde{S}$ are respectively defined as $N = \sum_{i=1}^n s_i \cdot \ell_i$ and $W = \sum_{i=1}^n \ell_i$. The “total” in the definition of $W$ and $N$ above, which we will often omit, refers to the fact that these are computed along the whole length of $\tilde{S}$.

For example, $\tilde{T}$ of Section I has length $n = 11$, width $W = 14$, and size $N = 24$. Moreover, there we have that all positions are solid except for 3 and 5, and hence all $\ell_i$'s
are equal to 1 except $\ell_3 = 3$ and $\ell_5 = 2$, and therein we have $s_3 = s_5 = 3$.

### III. EDIT DISTANCE WITH SKIPS

In this section we describe a dynamic programming algorithm to optimally compute the edit distance with affine gap penalty cost with skips between a string $P \in \Sigma^m$ and an E-string $\Upsilon \in (\Sigma \cup \{\epsilon\})^n$. A skip is an insertion of an empty letter $\epsilon$ of $\Upsilon$ whose cost $s$ will be 0 in the rest of this paper.<fnref n="4" />

The sequence of edit events that optimally turn a string into another, and that determine the alignment between the two strings, is described by an *edit transcript*: a sequence of events among $M, X, I, D, S$ (where $M$ denotes a match, $X$ is a mismatch, $I$ is an insertion, $D$ is a deletion, and $S$ denotes a skip), each one having its own cost, that transform a string into the other. The alignment is *optimal* when such cost is minimized. The affine gap penalty functions assigns a unique cost $o + k \cdot e$ to a gap of length $k$, that is a maximal run of $k$ consecutive insertions $I$ (resp. deletions $D$), with $o$ being for gap opening and $e$ for gap extension.

In an edit distance computation framework, the match has typically null cost in order to fulfill the mathematical metric requirement that a string has zero-distance to itself; for this reason in our examples we assume $a = 0$. However, since the dynamic programming method we design also works when one wants to compute a similarity score rather than a distance (it suffices to adapt the penalty scores and seek the maximum instead of the minimum), then in our problem statement as well as in the recurrence formula that describe our algorithm, we parametrize the score of a match with $a$.

**Example 1:** Let $\Upsilon = \text{GTCAAATGC}$, and let $a = 0$ be the cost of a match, $x = 1$ for a mismatch, and $o = 2$ (resp. $e = 1$) for gap opening (resp. extension). The edit transcript *MSMIIMMDDMSMX* is an alignment of $\Upsilon$ and $P$ with cost 9 (twice $o + 2e = 4$ for the two gaps of size two, and 1 for the final mismatch), while *MSMMXXMMSMX* is an optimal alignment of cost 3.

The solution of this problem, which we will name STOES, has the purpose of showing how to manage skips with affine gap penalty costs without adding the degenerate strings framework yet. In Section IV we will complete the picture showing how we will handle degenerate position in the dynamic programming. The formal statement of the problem STOES we solve here is the following:

**String Alignment to E-string (STOES)**

**Input:** An E-string $\Upsilon$ of length $n$, a pattern $P$ of length $m$, and penalty scores $a, x, o, e$.

**Output:** An optimal alignment between $P$ and $\Upsilon$ using scores: $a$ for match, $x$ for mismatch, $o$ for gap opening, and $e$ for gap extension.

Like in the dynamic programming (DP) framework for traditional strings, we incrementally compute the optimal alignments of prefix $P[0, j]$ of $P$ against prefix $\Upsilon[0, i]$ of $\Upsilon$, for growing values of $i$ and $j$ starting with $i$ or $j$ equal to 0,

up to the final lengths $n$ and $m$. As it is well known [50], the problem of computing the optimal alignment with such gap penalty function does not directly fulfill the requirement of having optimal substructure (which is needed to obtain an optimal solution with dynamic programming), and to overcome this limitation, the alignment is conceptually split in blocks grouping runs of edit events. To this purpose, the traditional computation of an optimal alignment with affine gap penalty function requires to use three matrices: $I$ for insertions, $D$ for deletions, $M$ for (mis)matches [57]. Since the skip operation we want to implement raises a similar issue, we here add an additional matrix $S$ for skips. As a consequence, our solution of STOES requires to compute the DP matrices $M, S, I, D$ of size $(n+1) \times (m+1)$ each, and whose entries for $0 \le i \le n$ and $0 \le j \le m$ have the following meaning:

$M[i, j]$ =score of best alignment of $P[0, j]$ and $\Upsilon[0, i]$, ending with a match or a mismatch of the letters $P[j]$ and $\Upsilon[i]$,

$S[i, j]$ = score of best alignment of $P[0, j]$ and $\Upsilon[0, i]$, ending with a skip at $\Upsilon[i]$,

$I[i, j]$ = score of best alignment of $P[0, j]$ and $\Upsilon[0, i]$, ending with an insertion of $\Upsilon[i]$ in $P$,

$D[i, j]$ = score of best alignment of $P[0, j]$ and $\Upsilon[0, i]$, ending with a deletion of $P[j]$ after $\Upsilon[i]$.

The initialization of the matrices is as follows:

**Matrix M:** Like in traditional dynamic programming with affine gap penalty, we set $M[0, 0] = 0$ and $M[i, 0] = M[0, j] = \infty$ for $0 < i \le n$ and $0 < j \le m$ because no alignment with an empty prefix can end with a match or with a mismatch.

**Matrix D:** Again like in traditional dynamic programming with affine gap penalty, we set $D[0, 0] = 0$, while for $0 < i \le n$ and $0 < j \le m$ we have $D[0, j] = o + e \cdot j$ because we align an empty prefix of $\Upsilon$ with $P[0, j]$, while $D[i, 0] = \infty$ because no alignment of $\Upsilon[0, i]$ and the empty prefix of $P$ ends with a deletion in $P$.

**Matrix I:** Dually, $I[0, 0] = 0$ and $I[0, j] = \infty$ because no alignment of $P[0, j]$ and an empty prefix of $\Upsilon$ can end with an insertion of a letter of $\Upsilon$. In order to set the values $I[i, 0]$ of the first column, instead, we need to take into account whether $\Upsilon[i] = \epsilon$ and, if not, then we need to pay the gap opening only the first time this is the case, and pay the gap extension only for successive insertions: let $i'$ be the smallest $i$ from 1 to $n$ such that $\Upsilon[i] \ne \epsilon$. Therefore we set:

$$ I[i, 0] = \min \begin{cases} \infty & \text{if } \Upsilon[i] = \epsilon \\ I[i - 1, 0] + o + e & \text{if } \Upsilon[i] \ne \epsilon \text{ & } i = i' \\ I[i - 1, 0] + e & \text{if } \Upsilon[i] \ne \epsilon \text{ & } i \ne i' \\ S[i - 1, 0] + e & \text{if } \Upsilon[i] \ne \epsilon \text{ & } i \ne i' \end{cases} \quad (1) $$

because we need to align $\Upsilon[0, i]$ and the empty prefix of $P$ in such a way to end with an insertion of $\Upsilon[i]$, and this is only possible if $\Upsilon[i] \ne \epsilon$, and else this either opens a new gap (when $i = i'$), or it extends an existing gap that was already

<footnote n="4">⁴The method can be generalized to the case of cost $s > 0$ for skips.</footnote>
open at row $i-1$ whose cost is either in the matrix $I$ itself, or in $S$ if $\Upsilon[i-1] = \varepsilon$, in which case the gap has included a skip at no cost.

**Matrix S:** We have that $S[0, j] = \infty$ because no alignment of an empty prefix of $\Upsilon$ and $P[0, j]$ can end with a skip of $\Upsilon$, while:

$$S[i, 0] = \begin{cases} \min\{S[i-1, 0], I[i-1, 0]\} & \text{if } \Upsilon[i] = \varepsilon \\ \infty & \text{if } \Upsilon[i] \neq \varepsilon \end{cases} \quad (2)$$

because we need to align $\Upsilon[0, i]$ and the empty prefix of $P$ ending with a skip. Being this not possible when $\Upsilon[i] \neq \varepsilon$, we set the value as $\infty$ when this is the case. When, instead, we have $\Upsilon[i] = \varepsilon$, then either a skip had occurred already at the previous row (in which case we can inherit its score), or we have to inherit the score of an insertion gap.

We remark that, besides the apparent circular definition of the first columns of $I$ and $S$, they can just both filled in a non ambiguous way once the index $i'$ has been located beforehand (which can be done in linear time).

We can then progressively fill in the matrices for $i, j > 0$ as follows:

$$M[i, j] = \min \begin{cases} \infty & \text{if } \Upsilon[i] = \varepsilon \\ f(i, j) + M[i-1, j-1] & \text{if } \Upsilon[i] \neq \varepsilon \\ f(i, j) + S[i-1, j-1] & \text{if } \Upsilon[i] \neq \varepsilon \\ f(i, j) + D[i-1, j-1] & \text{if } \Upsilon[i] \neq \varepsilon \\ f(i, j) + I[i-1, j-1] & \text{if } \Upsilon[i] \neq \varepsilon \end{cases} \quad (3)$$

where $f(i, j) = 0$ is only defined when $\Upsilon[i], P[j] \in \Sigma$, as follows:

$$f(i, j) = \begin{cases} a & \text{if } \Upsilon[i] = P[j] \\ x & \text{otherwise} \end{cases} \quad (4)$$

as in traditional dynamic programming with affine gap penalty, except that no (mis)match is defined in rows for which $\Upsilon[i] = \varepsilon$ because in that case no alignment between $P[0, j]$ and $\Upsilon[0, i]$ can end with a match or a mismatch as $\varepsilon$ is not a letter of the alphabet.

The deletion matrix $D$ is the very same as in the traditional alignment because deletions remove letters of $P$ which does not contain skips:

$$D[i, j] = \min \begin{cases} D[i, j-1] + e \\ M[i, j-1] + o + e \\ S[i, j-1] + o + e \\ I[i, j-1] + o + e \end{cases} \quad (5)$$

The DP recurrence relations to fill in the insertion matrix $I$ is as follows:

$$I[i, j] = \min \begin{cases} \infty & \text{if } \Upsilon[i] = \varepsilon \\ D[i-1, j] + o + e & \text{if } \Upsilon[i] \neq \varepsilon \\ M[i-1, j] + o + e & \text{if } \Upsilon[i] \neq \varepsilon \\ S[i-1, j] + e & \text{if } \Upsilon[i] \neq \varepsilon \\ I[i-1, j] + e & \text{if } \Upsilon[i] \neq \varepsilon \end{cases} \quad (6)$$

because when a skip is encountered at row $i$, then no costly insertion can take place and hence no alignment of $P[0, j]$ and $\Upsilon[0, i]$ can end with an insertion when $\Upsilon[i] = \varepsilon$. On the other hand, when $\Upsilon[i] \neq \varepsilon$, then we are basically in the same (actually dual) situation as with the filling of the deletion matrix $D$. Finally:

$$S[i, j] = \min \begin{cases} \infty & \text{if } \Upsilon[i] \neq \varepsilon \\ S[i-1, j] & \text{if } \Upsilon[i] = \varepsilon \\ M[i-1, j] & \text{if } \Upsilon[i] = \varepsilon \\ I[i-1, j] & \text{if } \Upsilon[i] = \varepsilon \\ D[i-1, j] & \text{if } \Upsilon[i] = \varepsilon \end{cases} \quad (7)$$

because if $\Upsilon[i] \neq \varepsilon$, then no skip operation can be performed. On the other hand, when $\Upsilon[i] = \varepsilon$, we can copy in $S[i, j]$ the best score in the cell above from any of the four matrices.

The computation can end when all matrices have been filled in up to the last row for $i = n$ and the last column for $j = m$, and then we have that the optimal score of the solution for STOES is $min(M[n, m], I[n, m], D[n, m], S[n, m])$.

An optimal edit transcript can be recovered by tracing back the computation of the optimal cost suitably jumping among the four matrices like in the traditional dynamic programming algorithm with affine gap score [50], [51], [58], with the only difference that here we have matrix $S$ as well. We remark that matrices $I$ and $S$ could possibly be merged into one to save space, but we maintain the distinction for the sake of clarity of their distinct meaning.

The correctness of the algorithm derives from the observation that any alignment either ends with a match, or a mismatch, or an insertion, or a deletions, or a skip, and that for each entry $(i, j)$ all these options have been correctly taken into account in the matrices as we defined. Also, with the incremental computation of optimal alignments of prefixes of $\Upsilon$ and $P$ using the recurrence relations that correctly refer to the optimal solution of all sub-problems of the possible sub-alignments that can precede them and according to the cost definition. In particular, our method has shown how to realize a null cost for skips operations. The algorithm we just described optimally solves the problem STOES in time and space $O(|\Upsilon| \cdot |P|)$ as this is the size of the matrices and each entry can be computed in constant time.

## IV. EDIT DISTANCE WITH SKIPS ON ED-STRINGS

In this section we show how to compute an optimal alignment with skips for a string $P$ and an ED-string $\tilde{T}$ with affine gap penalty score using the free skips rules we introduced in Section III. We can view an ED-string as the complete set of the linear strings with skips it represents: let us denote such set with $S(\tilde{T})$. Seeking an optimal alignment between a linear string $P$ and $\tilde{T}$ means to find the E-string $P' \in S(\tilde{T})$ that minimizes the edit distance with skips between $P$ and $P'$. For example, for the ED-string $\tilde{T}$ shown in Section I and $P = TATGCCCCCTTT$, we have that $P' = TA\epsilon TGTCCCCAT \in S(\tilde{T})$.
The problem we want to solve is defined as follows:

String Alignment to ED-string (STOED)

**Input:** An ED-string $\tilde{T}$ of length $n$, size $N$, and width $W$, a pattern $P$ of length $m$, and penalty scores $a, x, o, e$.

**Output:** An optimal alignment between $P$ and $\tilde{T}$ using scores: $a$ for match, $x$ for mismatch, $o$ for gap opening, and $e$ for gap extension.

The first step of our algorithm is a rewriting of the input ED-string $\tilde{T}$ that modifies degenerate positions $\tilde{T}\{i\}$'s in order to let all strings in $\tilde{T}\{i\}$ reach length $\ell_i$ by padding shorter strings with $\varepsilon$'s (we recall that $\ell_i$ is the length of the longest string in $\tilde{T}\{i\}$). For example, the ED-string $\tilde{T}$ of Section I becomes:

$$ \tilde{T}' = TA \begin{Bmatrix} G & T & C \\ A & T & \varepsilon \\ \varepsilon & \varepsilon & \varepsilon \end{Bmatrix} A \begin{Bmatrix} GT \\ TT \\ TG \end{Bmatrix} CCCCAT $$

We will refer to $\tilde{T}'$ as the padded version of $\tilde{T}$. Observe that the width $W$, the length $n$, and the size $N$ of $\tilde{T}'$ are the same as that of $\tilde{T}$. We can then build dynamic programming matrices with $W+1$ rows and $m+1$ columns, associating matrices rows $i'$ with the columns of $\tilde{T}'$: each one is one base wide, but some have several variants, and therefore refer to tuples of letters. In $\tilde{T}'$ above we have $W=14$ and the numbers below $\tilde{T}'$ show the (set of) letters that will be associated with rows of the dynamic programming; for example, $\tilde{T}[1] = A$, $\tilde{T}[2] = C$, $\tilde{T}[3]$ is the tuple $\langle \varepsilon, C, T \rangle$, $\tilde{T}[4] \langle \varepsilon, A, A \rangle$, etc. We denote with $\tilde{T}[i][h]$ the $h^{th}$ variant of the tuple $\tilde{T}[i]$.

The cells of the DP matrices within degenerate positions will thus contain tuples of values: one possible different edit distance per variant. Namely, we have now matrices $M'$, $S'$, $D'$, $I'$ whose entries are in general tuples. We denote with $M'[i,j][h]$ (with $1 \le h \le s$) the $h^{th}$ element of a tuple of size $s$ in $M'[i,j]$ (resp. in $S'$, $I'$, $D'$). Then the entries of the matrices are defined as follows:

$M'[i, j][h] =$ score of best alignment of $P[0, j]$ and $\tilde{T}[0, i][h]$ ending with a match or a mismatch of $P[j]$ and $\tilde{T}[i][h]$,

$I'[i, j][h] =$ score of best alignment of $P[0, j]$ and $\tilde{T}[0, i][h]$ ending with an insertion of $\tilde{T}[i][h]$ in $P$,

$D'[i, j][h] =$ score of best alignment of $P[0, j]$ and $\tilde{T}[0, i][h]$ ending with a deletion of $P[j]$ after $\tilde{T}[i][h]$,

$S'[i, j][h] =$ score of best alignment of $P[0, j]$ and $\tilde{T}[0, i][h]$ ending with a skip at $\tilde{T}[i][h]$.

We remark that we must ensure that, when an alignment of prefixes of $p$ and $\tilde{T}$ has chosen a variant $h$ in a degenerate position of $\tilde{T}$, then as long as rows involve variants of the same degenerate position, it must hold that the $h^{th}$ entry of all tuples always refers to the corresponding $h^{th}$ entry at previous rows and/or columns of the matrices.

A. THE MATCH MATRIX

We now show how the matrix $M'$ is filled in. The novelty with respect to the alignment to E-strings is the presence of

tuples, and hence the notion of match or mismatch is modified accordingly, as it may now need to refer to a specific item $\tilde{T}'[i][h]$ of the tuple $\tilde{T}'[i]$.

The function $f'(i, j)$ that checks whether we have a match or a mismatch when both the pattern $P$ and the text $T$ have letters from $\Sigma$ must then be modified for ED-strings from equation (4) as follows:

$$ f'(i, j)[h] = \begin{cases} a & \text{if } \tilde{T}'[i][h] = P[j] \\ x & \text{otherwise} \end{cases} \qquad (8) $$

to be used to compute the entries of $M'$ as follows when both $\tilde{T}'[i]$ and $\tilde{T}'[i-1]$ belong to the same degenerate position (and hence their tuples have the same size and letters of $\tilde{T}'[i][h]$ and $\tilde{T}'[i-1][h]$ belong to the same variant).

$$ M'[i, j][h] = \min \begin{cases} \infty & \text{if } \tilde{T}'[i][h] = \varepsilon \\ f'(i, j)[h] + M'[i-1, j-1][h] & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ f'(i, j)[h] + S'[i-1, j-1][h] & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ f'(i, j)[h] + D'[i-1, j-1][h] & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ f'(i, j)[h] + I'[i-1, j-1][h] & \text{if } \tilde{T}'[i][h] \neq \varepsilon \end{cases} \qquad (9) $$

Similarly, if $\tilde{T}'[i]$ and $\tilde{T}'[i-1]$ are both solid positions, then the formula (9) applies just ignoring $h$ because none of the involved entries is a tuple (basically formula (3) of the previous section can be used). More in general, we must take into account also the cases in which $\tilde{T}'[i]$ and $\tilde{T}'[i-1]$ are such that:

1) $\tilde{T}'[i-1]$ is solid and $\tilde{T}'[i]$ is degenerate

2) $\tilde{T}'[i-1]$ is degenerate and $\tilde{T}'[i]$ is solid

3) $\tilde{T}'[i-1]$ and $\tilde{T}'[i]$ are both degenerate but for two distinct consecutive degenerate positions.

In case 1. above, we have that $\tilde{T}'[i-1]$ is a letter of $\Sigma$ and $M'[i-1, j]$ is a single entry, while $\tilde{T}'[i]$ is a tuple of $\ell$ letters and $M'[i, j]$ is also a tuple of size $\ell$. In this case, whenever the computation of $M'[i, j]$ needs to refer to the previous row, then all variants $h$ of the tuple have to consider the same entry at row $i-1$, that is:

$$ M'[i, j][h] = \min \begin{cases} \infty & \text{if } \tilde{T}'[i][h] = \varepsilon \\ f'(i, j)[h] + M'[i-1, j-1] & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ f'(i, j)[h] + S'[i-1, j-1] & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ f'(i, j)[h] + D'[i-1, j-1] & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ f'(i, j)[h] + I'[i-1, j-1] & \text{if } \tilde{T}'[i][h] \neq \varepsilon \end{cases} \qquad (10) $$

In case 2. we have the inverse situation: $\tilde{T}'[i-1]$ is a tuple of letters and $M'[i-1, j]$ is tuple of values, while $\tilde{T}'[i]$ is a letter of $\Sigma$ and $M'[i, j]$ is a single entry. In this case, whenever the computation of $M'[i, j]$ needs to refer to the previous row, then all variants $h$ of the tuple at row $i-1$ have to be taken into account, and among them the dynamic programming
framework has to pick the value that minimizes the distance. That is:

$$ M'[i, j] = \min \left\{ \begin{array}{ll} \infty & \text{if } \tilde{T}'[i] = \varepsilon \\ f'(i, j) + \min_h \{M'[i-1, j-1][h]\} & \text{if } \tilde{T}'[i] \neq \varepsilon \\ f'(i, j) + \min_h \{S'[i-1, j-1][h]\} & \text{if } \tilde{T}'[i] \neq \varepsilon \\ f'(i, j) + \min_h \{D'[i-1, j-1][h]\} & \text{if } \tilde{T}'[i] \neq \varepsilon \\ f'(i, j) + \min_h \{I'[i-1, j-1][h]\} & \text{if } \tilde{T}'[i] \neq \varepsilon \end{array} \right. \quad (11) $$

Finally, In case 3. both $\tilde{T}'[i-1]$ and $\tilde{T}'[i]$ are tuples of letters (say of size $\ell_1$ and $\ell_2$, respectively), and consequently $M'[i-1, j]$ is a tuple of $\ell_1$ values, while $M'[i, j]$ is a tuple of $\ell_2$ values. Since they belong to distinct degenerate positions, whenever the computation of $M'[i, j]$ needs to refer to the previous row, there are no restrictions on which preceding alignment must be picked, and therefore it is again the minimum that must drive the choice. The recurrence formula is then as follows, where $h_2$ is the index on the tuple of $M'[i, j]$, and $h_1$ is that of $M'[i-1, j]$:

$$ M'[i, j][h_2] = \min \left\{ \begin{array}{ll} \infty & \text{if } \tilde{T}'[i][h_2] = \varepsilon \\ f'(i, j)[h_2] + \min_h \{M'[i-1, j-1][h_1]\} & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \\ f'(i, j)[h_2] + \min_h \{S'[i-1, j-1][h_1]\} & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \\ f'(i, j)[h_2] + \min_h \{D'[i-1, j-1][h_1]\} & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \\ f'(i, j)[h_2] + \min_h \{I'[i-1, j-1][h_1]\} & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \\ f'(i, j)[h_2] + \min_h \{I'[i-1, j-1][h_1]\} & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \end{array} \right. \quad (12) $$

**Example 2:** Figure 1 shows a sketch of matrix $M'$ with a zoom on the entries involved in a case in which a tuple is preceded and followed by a letter in $\tilde{T}$. For example, the computation of row $i$ of Figure 1 requires using equation (10) with $h=1$ and with $h=2$, meaning that two values (because there are two variants) have to be computed in all entries of row $i$ of $M'$ and that, in particular, entry $M'[i, j]$ is also a tuple of size 2. If $\tilde{T}[i][1]$ (resp. $\tilde{T}[i][2]$) is $\varepsilon$, then the corresponding entry $M'[i, j][1]$ (resp. $M'[i, j][2]$) is $\infty$. Else, letters $P[j]$ and $\tilde{T}[i][1]$ (resp. $\tilde{T}[i][2]$) are compared in order to compute $f'(i, j)[1]$ (resp. $f'(i, j)[2]$); the resulting value has to be added to the best one among the four entries $[i-1, j-1]$ of the four matrices $M'$, $S'$, $D'$, $I'$. Furthermore, since $\tilde{T}[i+1]$ represents the continuation of the same degenerate position that started in $\tilde{T}[i]$, then the computation of row $i+1$ of Figure 1 requires using equation (9) again with $h=1$ and with $h=2$, and checking entries $[(i+1)-1, j-1][h]$ of all matrices; the choice of $h$ for $\tilde{T}[i+1]$ being forced to be the same as that in the entries that have been computed for $\tilde{T}[i]$ ensures the

alignment staying on the same variant within the degenerate position. Finally, since $\tilde{T}[i+2]$ is again a solid letter after the degenerate position of $\tilde{T}[i, i+1]$ has ended, then the computation of row $i+2$ of $M'$ uses equation (11): entries of $M'$ at row $i+2$ are no longer tuples but rather single numbers making their best choice among entries $[(i+2)-1, j-1][h]$ for both $h=1$ and $h=2$.

## B. THE INDEL MATRICES

Also the Deletions and Insertions matrices $D'$ and $I'$ must be filled in using recurrences formula obtained modifying formulae (5) and (6) to take into account the case in which $\tilde{T}[i]$ is a tuple. Let us begin with $D'$, for which we do not even have to check whether $\tilde{T}[i, j]$ or $\tilde{T}[i, j][h]$ is equal to $\varepsilon$ because the edit operation there deletes a letter of $P$ whatever the letter of $\tilde{T}$ is. Hence, we have that

$$ D'[i, j][h] = \min \left\{ \begin{array}{l} D'[i, j-1][h] + e \\ M'[i, j-1][h] + o + e \\ S'[i, j-1][h] + o + e \\ I'[i, j-1][h] + o + e \end{array} \right. \quad (13) $$

because we just need to ensure that the same variant $h$ is read from the cell on the left when $\tilde{T}'[i]$ is a degenerate position. Else, if $\tilde{T}'[i]$ is a letter, we just apply the very same formula (13) disregarding any $"[h]"$, which is basically the same as applying equation (5).

In order to fill matrix $I'$ in, instead, we do need to check whether $\tilde{T}[i, j]$ or $\tilde{T}[i, j][h]$ is a letter or is equal to $\varepsilon$ because in it the content of $\tilde{T}[i]$ that is actually involved in the insertion. Moreover, we also need to check whether or not $\tilde{T}[i-1, j]$ is in the same degenerate position as $\tilde{T}[i, j]$ because the reference to the previous row of matrix $I'$ will have to change accordingly, as it was the case for $M'$. Therefore, we again need to distinguish all the cases that concern the relation between $\tilde{T}'[i][h]$ and $\tilde{T}'[i-1][h]$.

If both $\tilde{T}'[i]$ and $\tilde{T}'[i-1]$ belong to the same degenerate position, then row $i$ of $I'$ is a tuple whose values can be filled in as follows, for all $h$'s:

$$ I'[i, j][h] = \min \left\{ \begin{array}{ll} \infty & \text{if } \tilde{T}'[i][h] = \varepsilon \\ I'[i-1, j][h] + e & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ M'[i-1, j][h] + o + e & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ D'[i-1, j][h] + o + e & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ S'[i-1, j][h] + o + e & \text{if } \tilde{T}'[i][h] \neq \varepsilon \end{array} \right. \quad (14) $$

because, when $\tilde{T}'[i][h] = \varepsilon$ we have no letter of $\tilde{T}'[i]$ to insert. Otherwise, we insert the letter $\tilde{T}'[i, j][h]$ with the cost it takes, possibly extending an existing gap that was already open at $I'[i-1, j][h]$ with the very same variant $h$, or opening a new one inheriting the score from the entry $[i-1, j]$ of either of the other matrices, still provided this refers to the same variant $h$. Once again, the case in which both $\tilde{T}'[i]$ and $\tilde{T}'[i-1]$ are solid letters and not tuples, all we need to do
<figure><img src="image_1.png" /><figcaption>FIGURE 1. A sketch of matrix M' showing a case in which a tuple of size 2 is preceded and followed by a letter in T showing the consequent tuples in M' and the dependencies among its entries.</figcaption></figure>

is to apply equation (14) disregarding “[h]”, which actually means applying formula (6).

If we are in the case in which $\tilde{T}'[i-1]$ is solid whereas $\tilde{T}'[i]$ is degenerate, then we have to compute the tuple for entry $I'[i, j]$ for as many $h$'s as the number of variants at $\tilde{T}'[i]$, as follows:

$$I'[i, j][h] = \min \left\{ \begin{array}{ll} \infty & \text{if } \tilde{T}'[i][h] = \varepsilon \\ I'[i-1, j] + e & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ M'[i-1, j] + o + e & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ S'[i-1, j] + o + e & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ D'[i-1, j] + o + e & \text{if } \tilde{T}'[i][h] \neq \varepsilon \end{array} \right. \quad (15)$$

because there is no variant in $\tilde{T}'[i-1]$ to be coherent with, as $\tilde{T}'[i-1]$ is a single letter and hence entries $[i-1, j]$ contain a single value in all matrices. Therefore, which matrix the $I'[i-1, j]$ entry is inherited from when $\tilde{T}'[i][h] \neq \varepsilon$, will depend from what is their best score at entry $[i-1, j]$, or from whether or not an insertion gap had already been opened at $I'[i-1, j]$.

Conversely, if $\tilde{T}'[i-1]$ is degenerate and $\tilde{T}'[i]$ is solid, then the entry $I'[i-1, j]$ is a tuple of values (and so are $M'[i-1, j]$, $D'[i-1, j]$, and $S'[i-1, j]$), while $I'[i, j]$ will be a single value (like $M'[i, j]$, $D'[i, j]$, and $S'[i, j]$). Since the computation of $I'[i, j]$ needs to refer to the previous row, then all variants $h$ of the tuple at row $i-1$ have to be considered and, among them, we will choose the one that minimizes the distance. Again, this will also depend from whether or not an insertion gap had already been opened at $I'[i-1, j]$. Formally:

$$I'[i, j] = \min \left\{ \begin{array}{ll} \infty & \text{if } \tilde{T}'[i] = \varepsilon \\ \min_{h} \{I'[i-1, j][h]\} + e & \text{if } \tilde{T}'[i] \neq \varepsilon \\ \min_{h} \{M'[i-1, j][h]\} + o + e & \text{if } \tilde{T}'[i] \neq \varepsilon \\ \min_{h} \{S'[i-1, j][h]\} + o + e & \text{if } \tilde{T}'[i] \neq \varepsilon \\ \min_{h} \{D'[i-1, j][h]\} + o + e & \text{if } \tilde{T}'[i] \neq \varepsilon \end{array} \right. \quad (16)$$

The last case to be described is when both $\tilde{T}'[i-1]$ and $\tilde{T}'[i]$ are degenerate, but they are so for two distinct degenerate positions that happen to be consecutive (one is ending at $\tilde{T}'[i-1]$ and has size say $-\ell_1$, and another is beginning at $\tilde{T}'[i]$ and has size say $-\ell_2$). When this is the case, then $I'[i-1, j]$ is a tuple of $\ell_1$ values (and so are $M'[i-1, j]$, $D'[i-1, j]$, and $S'[i-1, j]$ from which we will have to read), while $I'[i, j]$ is a tuple of $\ell_2$ values, and we have to set them all. Let $1 \le h_1 \le \ell_1$ (resp. $1 \le h_2 \le \ell_2$). When $\tilde{T}'[i][h_2]$ is a letter, then the computation of $I'[i, j][h_2]$ will refer to the previous row of possibly any of the other matrices, selecting the best choice in terms of the minimum score. We remark that, even if both $\tilde{T}'[i-1]$ and $\tilde{T}'[i]$ are degenerate, and hence a variant ended at row $i-1$ of $I'$ and a new one begins at row $i$ of matrix $I'$, nothing forbids a gap to span over both variants, and therefore we also need to consider the case in which a gap is just going to be extended. The following recurrence formula computes $I'[i, j][h_2]$ for all $1 \le h_2 \le \ell_2$, where $h_1$ spans over the variants of $I'[i-1, j]$ (and hence $1 \le h_1 \le \ell_1$):

$$I'[i, j][h_2] = \min \left\{ \begin{array}{ll} \infty & \text{if } \tilde{T}'[i][h_2] = \varepsilon \\ \min_{h_1} \{I'[i-1, j][h_1]\} + e & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \\ \min_{h_1} \{M'[i-1, j][h_1]\} + o + e & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \\ \min_{h_1} \{S'[i-1, j][h_1]\} + o + e & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \\ \min_{h_1} \{D'[i-1, j][h_1]\} + o + e & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \end{array} \right. \quad (17)$$

### C. THE SKIP MATRIX

We finally show how the skip matrix $S'$ is filled in: this is realized by replacing equation (7) with equations (18),(19),(20) as shown below.

When both $\tilde{T}'[i]$ and $\tilde{T}'[i-1]$ belong to the same degenerate position (and hence their tuples have the same size and letters of $\tilde{T}'[i][h]$ and $\tilde{T}'[i-1][h]$ belong to the same variant), row $i$ of $S'$ is a tuple whose values can be filled in as follows, for all $h$'s:

$$S'[i, j][h] = \min \left\{ \begin{array}{ll} \infty & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ M'[i-1, j][h] & \text{if } \tilde{T}'[i][h] = \varepsilon \\ S'[i-1, j][h] & \text{if } \tilde{T}'[i][h] = \varepsilon \\ D'[i-1, j][h] & \text{if } \tilde{T}'[i][h] = \varepsilon \\ I'[i-1, j][h] & \text{if } \tilde{T}'[i][h] = \varepsilon \end{array} \right. \quad (18)$$

because, when $\tilde{T}'[i][h] = \varepsilon$ we can skip from $\tilde{T}'$ without incurring in any additional cost, and hence inherit the score from the entry $[i-1, j]$ above just copying its value, provided it refers to the same variant $h$; this holds true whatever matrix it comes from, and we just pick the minimum among all other matrices. If, instead, $\tilde{T}'[i][h] \neq \varepsilon$, then no skip can take place at this row, and therefore we set $S'[i, j][h] = \infty$.

When, instead, $\tilde{T}'[i-1]$ is solid and $\tilde{T}'[i]$ is degenerate, then we have to compute the tuple for entry $S'[i, j]$ as follows
for as many h's as the number of variants at $\tilde{T}'[i]$:

$$ S'[i, j][h] = \min \begin{cases} \infty & \text{if } \tilde{T}'[i][h] \neq \varepsilon \\ M'[i-1, j] & \text{if } \tilde{T}'[i][h] = \varepsilon \\ S'[i-1, j] & \text{if } \tilde{T}'[i][h] = \varepsilon \\ D'[i-1, j] & \text{if } \tilde{T}'[i][h] = \varepsilon \\ I'[i-1, j] & \text{if } \tilde{T}'[i][h] = \varepsilon \end{cases} \quad (19) $$

because the inherited entry can still come from any of the four matrices, but now there is no variant to be coherent with as $\tilde{T}'[i-1]$ was a single letter and hence entries $[i-1, j]$ contain a single value in all matrices.

When $\tilde{T}'[i-1]$ is degenerate and $\tilde{T}'[i]$ is solid, then we just set $S'[i, j] = \infty$ because it must be $\tilde{T}'[i] \neq \varepsilon$ and hence no alignment can end at $\tilde{T}'[i]$ with a skip event.

Finally, when both $\tilde{T}'[i-1]$ and $\tilde{T}'[i]$ are degenerate, but they are so for two distinct degenerate positions that are consecutive (say they are of size $\ell_1$ and $\ell_2$, respectively), then $S'[i-1, j]$ is a tuple of $\ell_1$ values, while $S'[i, j]$ is a tuple of $\ell_2$ values. Since they belong to distinct degenerate positions, whenever the computation of $S'[i, j]$ needs to refer to the previous row, there are no restrictions on which preceding alignment should be picked to add up its cost, and therefore it is again the minimum score that has to drive the choice. The following recurrence formula computes the entry $S'[i, j][h_2]$ for all $1 \le h_2 \le \ell_2$, and refers to $h_1$ that spans over all the variants of any of the matrices at row $i-1$ (and hence $1 \le h_1 \le \ell_1$):

$$ S'[i, j][h_2] = \min \begin{cases} \infty & \text{if } \tilde{T}'[i][h_2] \neq \varepsilon \\ \min_{h_1}\{M'[i-1, j][h_1]\} & \text{if } \tilde{T}'[i][h_2] = \varepsilon \\ \min_{h_1}\{S'[i-1, j][h_1]\} & \text{if } \tilde{T}'[i][h_2] = \varepsilon \\ \min_{h_1}\{D'[i-1, j][h_1]\} & \text{if } \tilde{T}'[i][h_2] = \varepsilon \\ \min_{h_1}\{I'[i-1, j][h_1]\} & \text{if } \tilde{T}'[i][h_2] = \varepsilon \end{cases} \quad (20) $$

As an example, Figure 2 shows a sketch of matrix $S'$ showing, for the same case as that of Figure 1, a zoom over two entries that belong to a degenerate position, as well as the rows in $S'$ that precedes and follows them, with the dotted line highlighting dependencies. In particular, we can see that also in $S'$ when $\tilde{T}'[i]$ we have that row $i$ contains tuples. Then, According to equations (18),(19), and (20), in the matrix of Figure 2 we will have $S'[i, j][h] = \infty$ (resp. $S'[i+1, j][h] = \infty$) when $\tilde{T}'[i][h] \neq \varepsilon$ (resp. $\tilde{T}'[i+1][h] \neq \varepsilon$). On the other hand, if $\tilde{T}'[i][h] = \varepsilon$, then $S'[i, j][h]$ is computed for $h=1$ and $h=2$ according to equation (19) choosing to inherit from the best among the entries $[i-1, j]$ of the four matrices. As for row $i+1$ of $S'$, since $\tilde{T}[i+1]$ represents the continuation of the same degenerate position that started in $\tilde{T}[i]$, then the computation of $S'[i+1, j][h]$ requires using equation (18) again with $h=1$ and with $h=2$, and checking entries $[(i+1)-1, j-1][h]$ of all matrices; the choice of $h$ for $\tilde{T}[i+1]$ is thus forced to be the same as that in the

<figure><img src="image_2.png" /><figcaption>FIGURE 2. Sketch of matrix $S'$ showing, for the same case as Figure 1, the consequent tuples in $S'$ and the dependencies among its entries.</figcaption></figure>

entries that have been computed for $\tilde{T}[i]$, making sure the alignment stays on the same variant within the degenerate position. Finally, since $\tilde{T}[i+2]$ is again a solid letter after the degenerate position of $\tilde{T}[i, i+1]$ has ended, then the computation of row $i+2$ of $S'$ uses equation (20): entries of $S'$ at row $i+2$ are no longer tuples but rather single numbers making their best choice among entries $[(i+2)-1, j-1][h]$ for both $h=1$ and $h=2$.

The computation ends when all matrices have been filled in up to the last row for $i=W$ and the last column with $j=m$, and then we have that the final result is $\min(M'[W, m], I'[W, m], D'[W, m], S'[W, m])$, and an optimal edit transcript can be recovered by tracing back the computation of the optimal cost suitably jumping among the four matrices.

The algorithm above correctly solves problem STOEDS in time and space in $O(m \cdot N)$ where $m = |P|$ and $N$ is the size of $\tilde{T}$, and its trace-back reconstruct an explicit edit transcript with skips.

We conclude this section showing, with a very simple example of the difference between an edit transcript with, or without, the skip operation made explicit.

**Example 3:** Let us consider the simple ED string $\tilde{T} = A\{C\}T$ and the string $P = AT$. For the sake of simplicity, here we ignore the requirement of affine gap penalty function and hence that of computing four matrices (the optimal alignment contains no gaps anyhow). In Fig. 3, on the left we show the alignment obtained with just POA (the partial order alignment is represented by the arrows on the left) producing the edit transcript MM that correctly records the match between the pattern $P$ and its occurrence $P' = AT$ in $\tilde{T}$ disregarding the reason why the match score has been propagated from an entry in earlier rows of the matrix. In Fig. 3, on the right
<figure><img src="image_3.png" /><figcaption>**FIGURE 3.** Alignments of ED string $\tilde{T} = A\{\epsilon, C\}T$ and string $P = AT$ with scoring scheme 0 for match and skips, and 1 for mismatch and gap. The alignment on the left does not make use of the skip operation, and uses only partial order alignment (and hence the last row has the two rows above as possible preceding line). The alignment on the right makes use of the skip operation.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**FIGURE 4.** An example of Wave Front Alignment of strings ACGTT and ACCTA where any edit operation (mismatch, insertion, or deletion) scores 1. In this case no INDELs is required and the final cell is reached by WF₂ (and hence the value 2 is reported). Light grey squares correspond to entries of the matrix whose values do not need to be computed.</figcaption></figure>

we show a single matrix<fnref n="5" /> for the alignment we compute with skips. Since the optimal alignment makes use of a skip operation (as the vertical dotted line shows the propagation of score 0 because of $\epsilon$ at that row), the resulting edit transcript is MSM.

## V. WAVE FRONT ALIGNMENT WITH AFFINE GAP PENALTY AND SKIPS

When two strings are similar, their alignment can be performed in an efficient way by means of the so-called *Wave Front Alignment* (WFA). This powerful technique was introduced in [59] as a linear time and space dynamic programming method to compute the edit distance of similar strings. In a few words, wave front alignment works as follows: the input are two strings of size *N* whose distance is known to be limited. Assuming for simplicity that each edit operation (among insertion, deletion, and mismatch) costs 1 (that is, no affine gap penalty is used for now), and assuming that their distance is upper-bounded by *D*, this information is exploited by avoiding to compute the whole dynamic programming table of size *N*², and actually computing only a stripe of *O*(*D*) diagonals around the main one. Indeed, the similarity guarantees that the optimal alignment lays therein: if at most *D* edit operation are required to turn a string into the other, then the optimal alignment will move at most *D* diagonals away from the main one. One begins with trying to advance along the main diagonal pushing a null distance as far as possible as long as there are only matching characters: this is WF₀, the wave front for distance 0. When this is no longer possible, WF₁ proceeds in the very same way along the main diagonal after having payed the cost of a mismatch, as well as along its adjacent diagonals after having payed the cost of an INDEL. This procedure is incrementally iterated from WFᵢ to WFᵢ₊₁ for i ≥ 1 until the bottom-right cell of the matrix is reached. When this happens for, say, WFₐ the computation ends and the value of *d* is returned. The result is a *O*(*DN*) time and space algorithm replacing the quadratic one, and an algorithm that, rather than filling in a matrix row by row, proceeds with wave fronts along diagonals using a

coordinate system of pairs < diagonal, offset > instead of the traditional < row, column > to navigate the matrices.

A simple example for linear strings and disregarding the affine gap penalty score, is shown in Figure 4, where WF₀ reaches offset 2 opening three possible paths for WF₁ along as many diagonals. The furthest reaching offset of WF₁ (along diagonal 0) is thus extended by WF₂ up to the final cell. Once this is reached, no further wave fronts are explored and the computation ends reporting the distance of the latest wave front.

In [60], the linear time and space optimization of Myers was extended to the exact computation of an optimal pairwise alignment of strings using the affine gap penalty score function. Hence, the wave front proceeds iterating over a similarity score rather than over a distance and, moreover, both opening and extension scores are managed for the gaps. Roughly, for the alignment of two strings of size $N_1$ and $N_2$, the matrices *I*, *D*, *M* of size $N_1 \times N_2$ are replaced by *wave front* records $I_{d,k}$, $D_{d,k}$, $M_{d,k}$ that store, for each score/distance *d* and diagonal *k*, the furthest-reaching offset in *k* that scores *d*. The dynamic programming recurrence is then operated on increasing values of *d* and adjacent diagonals *k*, up to the final distance: the one whose wave front reaches the end of both input strings in any of the records *I*, *D*, or *M*.

In [17], the Wave Front Alignment with Affine Gap Penalty scores was further extended to work with Degenerate Strings (D strings) that can be viewed as ED strings without the empty string and with the constraint of having all variants of the same length $l_i$ within $\tilde{T}\{i\}$. As both theoretical and experimental analyses reported in [17] demonstrate, this results in a very fast tool that finds the optimal alignment of a linear string and a degenerate string in a time which is linear with the width *W* of the degenerate string. The method, besides being fast, enjoys the property of not being an heuristic, thus giving a guarantee of the optimal solution and its corresponding alignment and transcript.

As we have shown in Section IV, the first step of our algorithm that solves STOEDS is a padding that rewrites the ED-string $\tilde{T}$ by modifying its degenerate positions $\tilde{T}\{i\}$'s in

<footnote n="5">⁵It is a matrix unifying matrices $M'$ and $S'$ showing in each entry the value of the best score among them.</footnote>
order to let all strings in $\tilde{T}\{i\}$ reach length $\ell_i$ by appending $\epsilon$'s on the right of possible shorter strings. This basically results into a special kind of Degenerate String. The strings we turn $\tilde{T}$ into after the padding step is basically a D string where the symbol $\epsilon$ might also occur in one or more variants after a symbol in $\Sigma$. The efficient Wave Font Algorithm introduced in [17] for D strings can be extended in a straightforward way to the padded ED strings for which we described here the alignment algorithms with skips. As a result, the technique suggested here can be combined with that of [17] to obtain a fast dynamic programming algorithm that uses wave front alignment as well as partial order to compute the optimal alignment of an ED string and a linear string, with affine gap penalty.

In particular, it is enough to keep four wave front records $I_{d,k}$, $D_{d,k}$, $M_{d,k}$, and $S_{d,k}$ to represent the current furthest reaching wave fronts on all matrices that here include the skip matrix too. These records can be updated using the algorithm we have shown in Section IV for as much as skips concerns, and maintaining the combined use of affine gap penalty score and wave fronts shown in [17] for degenerate strings. The only modification that should be done to the formula od Section IV are those that realize the change of coordinate system using pairs < *diagonal*, *offset* > that represent furthest reaching points of waves along diagonals, rather than pairs < *row*, *column* > to navigate the matrices. This can again be done in the very same way as in [17].

Observing that wave fronts proceed in general along the records $I_{d,k}$, $D_{d,k}$, $M_{d,k}$, and $S_{d,k}$, and that an entry $\tilde{T}\{i\}[h]$ for each $1 \le h \le \ell_i$ (resp. $\tilde{T}\{i\}$ if this is a solid letter) has a unique wave front that reaches it before any other, we have that the problem STOEDS can be solved in time $O(Nd)$ where $N$ is the size of $\tilde{T}$ and $d$ is the returned distance.

As for space requirement, the complexity coincides with the largest amount of furthest reaching wave front records that are simultaneously active for any current distance $d$, which is again upper bounded by $O(Nd)$. Concerning space complexity, we additionally, observe that matrix $M$ and matrix $S$ can actually be merged into a single matrix because each entry $M[i,j]$ (resp. $S[i,j]$) therein, as well as each “sub-entry” $M[i,j][h]$ (resp. $S[i,j][h]$) for each $1 \le h \le \ell_i$ when $\tilde{T}\{i\}$ is not a solid letter, is either a letter of $\Sigma$ or it is an $\epsilon$. In the former case the (sub)entry would not contain a finite number in $S$, and in the latter it would not contain a finite number in $M$, and therefore the useful entries of the two matrices are somewhat complementary, and their contents can be merged into a single matrix or wave front record.

As an example, Figure 5 shows the wave fronts along a matrix that simultaneously contains entries of $M$ and of $S$: the dark blue dot at the left on the tuple in the fifth row and fourth column is the only entry that refers to a $\epsilon$, and hence would be recorded only in the wave front record of $S$ for $d = 2$. All other colored dots correspond to entries visited by wave fronts that refer to solid letters in $\tilde{T}$, and hence would instead be all recorded in the wave front record of $M$ (for $d = 0, 1, 2$)

<figure><img src="image_5.png" /><figcaption>**FIGURE 5.** An example of Wave Front Alignment of $\tilde{T} = AS\{T, GA\}A$ and $P = AAGA$ where we also show wave fronts for $d = 0, 1, 2$ with different colors, as well as the final optimal alignment and edit transcript. Here mismatches and gaps extension cost 1, a gap opening costs 2, while matches and skips are free. The wave fronts of matrix $M$ and of matrix $S$ are shown here all in a single table.</figcaption></figure>

because they are all reached with a match or with a mismatch of that letter.

As illustrated in Fig. 5, the (pink) wave front for distance $d = 0$ halts after matching the two A's, getting to a single furthest reaching point in the main diagonal. When $WF_0$ stops due to a mismatch, $WF_1$ is activated for distance $d = 1$ (the cost of the mismatch) along the same diagonal. The (light blue) wave front $WF_1$ advances along the same diagonal as the right variant of the degenerate position matches the remaining letters GA of $P$. However, $WF_1$ will stop at the fifth row and fifth column, because the only way to complete the alignment would be to insert the last letter of $\tilde{T}$, but opening the gap would cost 3 (a gap opening cost 2 plus a gap extension that costs 1), and hence this would be deferred to $WF_4$ that would activate the wave front record in the insertion matrix. Additionally, the wave front $WF_1$ has another furthest reaching point left after the mismatch which opens now $WF_2$: incrementing the distance to 2 with a mismatch choosing the other variant, and then proceeding down with the skip, and closing with a final match. At this point $WF_2$ has reached the bottom right cell of the matrix, and hence the computation ends returning 2. The optimal alignment is highlighted in yellow on the right with its edit transcript. Notably, this alignment involves no deletions nor insertions, and hence these matrices/records will remain untouched. This is why Fig. 5 displays only the relevant entries in the (Mis)Match and Skip matrices. Finally, all empty cells in the matrix of Fig. 5 represent entries that are never explored because they would result in scores exceeding that of the optimal alignment.

Algorithm 1 shows the pseudocode of WFSALIGN (where WFS stands for Wave Front with Skips) that takes as input an ED string $\tilde{T}$ and a string $P$, and computes their optimal alignment using a wave front algorithm within the affine gap penalty function framework, following the structure of [60] using the wave front coordinates system of diagonals and
offsets rather than rows and columns indices. However, our approach works with ED strings (and hence in particular with the empty string and skips) and their tuples-based partial order alignment, as shown in Section IV. Since the padding does not affect the width nor the size of $\tilde{T}$, we assume without loss of generality that $\tilde{T}$ is already padded.

We denote with $\text{offset}_{d,k}^M$ (resp. $\text{offset}_{d,k}^S$, $\text{offset}_{d,k}^I$, $\text{offset}_{d,k}^S$) the offsets of the wave front of matrix $M$ (resp. matrix $D, I, S$) along diagonal $k$ for distance $d$. Observe that, for any value of $d$, in general the wave front $WF_d^M$ of $M$ (resp. $WF_d^S, WF_d^I, WF_d^S$ of $D, I, S$) may span over multiple diagonals.

Note that only a match or a skip can “consume” letters at zero cost. As a result, the following holds:

* the only wave front offset that can be extended within the same distance $d$ and diagonal $k$ is $\text{offset}_{d,k}^M$ because, in case of matches along $k$, the wave front of $M$ can increase the offset maintaining distance $d$;

* still within the same distance $d$, whenever an offset encounters a skip along $\tilde{T}$, then a new wave front $\text{offset}_{d,k-1}^S$ in the diagonal $k-1$ below (assuming $\tilde{T}$ is along rows and $P$ along columns) can be opened up.

Both $WF_d^M$ and $WF_d^S$ should be extended in this way whenever, and as long as, any of the above conditions holds. The function WFSROLL realizes these free extensions up to the furthest reaching offsets of $WF_d^M$ and $WF_d^S$, starting from all opened wave fronts for $d$. Observe that a letter of $\tilde{T}$ is either $\epsilon$ or a character in $\Sigma$, and hence $WF_d^M$ and $WF_d^S$ may alternate each other to extend $WF_d$, but they will never cross nor stop at the same spot. For example, a run of matches can be interrupted at a certain $\text{offset}_{d,k}^M$ because it encounters a skip instead of a letter: when that is the case $\text{offset}_{d,k}^M$ will remain set there for the records (for possible wave fronts with distances larger than $d$), while $WF_d$ will get extended downward with $\text{offset}_{d,k-1}^S$.

Notice that, to avoid clatter, in the explanation above for WFSROLL as well as below for WFSNEXT, we disregard the fact that $\text{offset}_{d,k}^*$ can actually be tuples and that in that case the values will have to be computed accordingly as shown in Section IV.

When WFSROLL can no longer extend $WF_d^M$ nor $WF_d^S$, the starting points of $WF_{d+1}$ must be computed. This is done by first incrementing $d$, and then calling the function WFSNEXTM (resp. WFSNEXTD, WFSNEXTI) that computes $WF_d^M$, (resp. $WF_d^D, WF_d^I$), setting all possible $\text{offset}_{d,k}^M$ (resp. $\text{offset}_{d,k}^S, \text{offset}_{d,k}^I$) in all diagonals $k$ where this is possible.

Let us consider for example WFSNEXTM: scores referring to alignment ending with matches have been computed by WFSROLL before $d$ was incremented, and hence WFSNEXTM sets the baseline for $WF_d^M$ in the hypothesis that the latest event was a mismatch. An alignment scoring $d$ that ends with a mismatch can be obtained whenever in $\tilde{T}$ there is a letter mismatching that of $P$ and, along the same diagonal, the furthest reaching offset of any $WF_{d-x}$ had been reached, where by “any” we mean in any wave front among $WF_{d-x}^M$,

$WF_{d-x}^D$, $WF_{d-x}^I$, or $WF_{d-x}^S$. Specifically, for each possible $\text{offset}_{d-x,k}^* = f$ in $WF_{d-x}^*$ (for $* = M, D, I, S$), the function WFSNEXTM will add $\text{offset}_{d,k}^M = f+1$ to $WF_d^M$.

Summing up, WFSALIGN shown in Algorithm 1 works as follows.

An initialization phase sets:

* [line 1] the width $W$ of $\tilde{T}$ that would correspond to the number of rows of the matrices,

* [line 2] the central diagonal $A_k$ along which the alignment must begin and end,

* [line 3] the offset end that, if along $A_k$, corresponds to the bottom right entry of the matrices,

* [line 4] the initial null distance,

* [line 5] the initial tuple of size $\ell_1 = |\tilde{T}\{1\}|$ for $\text{offset}_{0,A_k}^M$. The offsets of the other wave fronts $WF^S, WF^D$, and $WF^I$ will be generated and initialized by WFSROLL (for $WF^S$), WFSNEXTD (for $WF^D$), and WFSNEXTDI (for $WF^I$).

An iterative phase that, for increasing values of $d$ until the alignment is complete, computes $WF_d$ as follows:

* [line 7] the function WFSROLL extends $WF_d$ as much as possible by possibly “rolling” $WF_{d,k}^M$ along $k$ and $WF_{d,k}^S$ down to $k' < k$ as explained above,

* [line 9] since $WF_d$ has reached its furthermost boundaries, $d$ is incremented,

* [line 10] the function WFSNEXTM computes the new wave front $WF_d^M$ from $WF_{d-x}$ by paying the cost of a mismatch: every $\text{offset}_{d-x,k}^*$ (for $* = M, D, I, S$) generates a new $\text{offset}_{d,k}^M$,

* [line 12] the function WFSNEXTD computes the new wave front $WF_d^D$, which can be either

- from $WF_{d-e}^D$ by paying the cost $e$ of extending an existing deletion gap: every $\text{offset}_{d-e,k}^D$ can generate a new $\text{offset}_{d,k+1}^D$ on the diagonal above, or

- from $WF_{d-o-e}^*$ (for $* = M, I, S$) by paying the cost $o + e$ of opening a new deletion gap: every $\text{offset}_{d-o-e,k}^*$ can generate a new $\text{offset}_{d,k+1}^D$,

* [line 14] dually, the function WFSNEXTI computes the new wave front $WF_d^I$ either from $WF_{d-e}^I$ extending an existing insertion gap, or from any $WF_{d-o-e}^*$ (for $* = M, D, S$) opening a new one.

Lines 8, 11, 13, 15 and 16 manage the conditions to terminate WFSALIGN. Observe that the alignment is complete whenever a wave front offset reaches the value end along diagonal $A_k$. At lines 8, 11, 13, and 15 we set the furthest reaching value along the main diagonal $A_k$ for $WF_d^S, WF_d^M, WF_d^D$, and $WF_d^I$, respectively. This is done finding the maximum value in a possible tuple, and right after the latest updated done to the corresponding wave front. Finally, line 16 checks whether any of them has reached the ending value.

## VI. PROGRESSIVE MSA WITH SKIPS

As we mentioned, a progressive multiple sequence alignment reduces the MSA to a linear profile and then iterates a pairwise alignment task between the linearized profile and
**Algorithm 1** The Wave Front Alignment Algorithm With Skips

0: **function** WFSALIGN($\tilde{T}$, P)

1: W ← width of $\tilde{T}$

2: $A_k ← W - |P|$ {central diagonal}

3: end ← max{$W, |P|$} {end of alignment index}

4: d ← 0 {initial distance}

5: offset$^{M}_{0,A_k}$ ←<0, ..., 0> {initialize tuple of length $\ell_1$}

6: repeat

7: WFSROLL($WF_d$) {extends the wave front within $d$}

8: $f_S ← \max \text{offsets}^S_{d,A_k}$ {furthest skip on $A_k$}

9: $d ← d + 1$ {increases current distance} {baselines of the next wavefront computed below:}

10: WFSNEXTM($WF^M_{d-x}, WF^D_{d-x}, WF^I_{d-x}, WF^S_{d-x}, d$)

11: $f_M ← \max \text{offsets}^M_{d,A_k}$ {furthest (mis)match on $A_k$}

12: WFSNEXTD($WF^D_{d-e}, WF^M_{d-o-e}, WF^I_{d-o-e}, WF^S_{d-o-e}, d$)

13: $f_D ← \max \text{offsets}^D_{d,A_k}$ {furthest deletion on $A_k$}

14: WFSNEXTIWF^I_{d-e}, WF^M_{d-o-e}, WF^D_{d-o-e}, WF^S_{d-o-e}, d

15: $f_I ← \max \text{offsets}^I_{d,A_k}$ {furthest insertion on $A_k$}

16: **until** ($f_M \ge \text{end}$) OR ($f_S \ge \text{end}$) OR ($f_D \ge \text{end}$) OR ($f_I \ge \text{end}$)

the new sequence. Partial Order Alignment (POA) [55] overcomes problems from earlier progressive methods by avoiding to reduce partial MSA to a linear profile, but rather aligning against a graph-like text and then more correctly computing the edit distance. However, partial order alignment does not yield the optimal edit transcript because the traceback does not reflect the distance introduced by variants with gaps of different lengths. This problem is solved here with the skip operation.

In this section we show how progressive multiple alignment can be improved using our optimal solution for STOEDS by performing the following steps, assuming that the order in which the sequences are aligned is known, e.g. guided by a tree that is given in input or that can be computed as a preprocessing step [53], [54], [61]:

1) Performing a first pairwise global alignment that produces a partial MSA that is turned into an ED-string.

2) Progressively adding a new string to the MSA.

3) Producing a new richer MSA as a result of Step 2. and hence a new ED-string.

4) Iterate steps 2. and 3. until all input sequences have been included.

Step 1 can be performed turning an edit transcript of the pairwise alignment into an ED-string as follows: runs of matches will correspond to solid positions; a mismatch can be expressed with a degenerate position with two variants and unitary width; gaps can be represented by degenerate positions that include an empty string among its options. Step 2. can be performed with our optimal solution for STOEDS. Step 3. can be done similarly to Step 1, with new degenerate positions that include the empty string being added when insertions and deletions appear in the edit transcript, while skips in the ED-string will just confirm existing ones in the multiple alignment.

## VII. CONCLUSION AND FURTHER WORK

We remark that with our dynamic programming method, it is very natural to add weights to letters and to use a weighted modification of the actual score as objective function. This could be exploited by using two kinds of metadata that are typically available in genomic sequences: (i) confidence score, and (ii) variants abundance.

**Confidence score:** Given the low but unavoidable rate of sequencing errors, all DNA sequencing machines actually output metadata associating a confidence score to each nucleic base (that is, the letter in the DNA sequence), indicating how much one should trust that specific base call. Our algorithm could associate this possibly useful metadata to the query and/or the ED string.

**Variants abundance:** In ED strings, like in pangenome graph in general, one may add abundance information to the strings that represent the variants (that is, how many of the input linear alleles exhibit a specific variant). Tracking this information in alignments could give a measure of how frequent is the aligned variants of the pangenome. Also, abundance information in a MSA could also be used to account the frequency of the variants in the consensus ED string, and our DP approach could naturally be updating this information for the next iteration.

In both cases managing the weights does not affect the theoretical complexity, but their treatment would slow down the running time by a constant factor, and therefore the use of this feature should be optional and its use evaluated in specific applications trading off accuracy and speed.

## REFERENCES

[1] T. Marschall, "Computational pan-genomics: Status, promises and challenges," *Briefings Bioinf.*, vol. 19, no. 1, pp. 118-135, 2018.

[2] W. W. Liao et al., "A draft human pangenome reference," *Nature*, vol. 617, no. 7960, pp. 312-324, 2023. [Online]. Available: https://www.nature.com/articles/s41586-023-05896-x

[3] J. A. Baaijens, P. Bonizzoni, C. Boucher, G. D. Vedova, Y. Pirola, R. Rizzi, and J. Sirén, "Computational graph pangenomics: A tutorial on data structures and their applications," *Natural Comput.*, vol. 21, no. 1, pp. 81-108, Mar. 2022.

[4] B. Paten, A. M. Novak, J. M. Eizenga, and E. Garrison, "Genome graphs and the evolution of genome inference," *Genome Res.*, vol. 27, no. 5, pp. 665-676, May 2017.

[5] E. Garrison, J. Sirén, A. M. Novak, G. Hickey, J. M. Eizenga, E. T. Dawson, W. Jones, S. Garg, C. Markello, M. F. Lin, B. Paten, and R. Durbin, "Variation graph toolkit improves read mapping by representing genetic variation in the reference," *Nature Biotechnol.*, vol. 36, no. 9, pp. 875-879, Sep. 2018.

[6] J. M. Eizenga, A. M. Novak, E. Kobayashi, F. Villani, C. Cisar, S. Heumos, G. Hickey, V. Colonna, B. Paten, and E. Garrison, "Efficient dynamic variation graphs," *Bioinformatics*, vol. 36, no. 21, pp. 5139-5144, Jan. 2021.

[7] J. Sirén, E. Garrison, A. M. Novak, B. Paten, and R. Durbin, "Haplotype-aware graph indexes," in *Proc. 18th Int. Conf. Algorithms Bioinf. (WABI)*, vol. 36, 2019, pp. 400-407.

[8] G. Rakocevic et al., "Fast and accurate genomic analyses using genome graphs," *Nature Genet.*, vol. 51, no. 2, pp. 354-362, Feb. 2019.

[9] N. Rizzo and V. Mäkinen, "Linear time construction of indexable elastic founder graphs," in *Proc. 33rd Int. Workshop Combinat. Algorithms (IWOCA)*, 2022, pp. 480-493.
[10] N. Rizzo, M. Equi, T. Norri, and V. Mäkinen, "Elastic founder graphs improved and enhanced," *Theor. Comput. Sci.*, vol. 982, Jan. 2024, Art. no. 114269.

[11] C. S. Iliopoulos, R. Kundu, and S. P. Pissis, "Efficient pattern matching in elastic-degenerate texts," in *Proc. 11th Int. Conf. Language Automata Theory Appl. (LATA)*, vol. 10168, 2017, pp. 131–142.

[12] R. Grossi, C. S. Iliopoulos, C. Liu, N. Pisanti, S. P. Pissis, A. Retha, G. Rosone, F. Vayani, and L. Versari, "On-line pattern matching on similar texts," in *Proc. 28th Annu. Symp. Combinat. Pattern Matching (CPM)*, vol. 78, 2017, p. 14.

[13] C. S. Iliopoulos, R. Kundu, and S. P. Pissis, "Efficient pattern matching in elastic-degenerate strings," *Inf. Comput.*, vol. 279, Aug. 2021, Art. no. 104616.

[14] V. Mäkinen, B. Cazaux, M. Equi, T. Norri, and A. I. Tomescu, "Linear time construction of indexable founder block graphs," in *Proc. 20th Int. Conf. Algorithms Bioinf. (WABI)*, vol. 172, 2020, p. 18.

[15] M. Alzamel, L. A. K. Ayad, G. Bernardini, R. Grossi, C. S. Iliopoulos, N. Pisanti, S. P. Pissis, and G. Rosone, "Degenerate string comparison and applications," in *Proc. 18th Int. Workshop Algorithms Bioinf. (WABI)*, 2018, pp. 21:1–21:14.

[16] M. Alzamel, L. A. K. Ayad, G. Bernardini, R. Grossi, C. S. Iliopoulos, N. Pisanti, S. P. Pissis, and G. Rosone, "Comparing degenerate strings," *Fundamenta Informaticae*, vol. 175, nos. 1–4, pp. 41–58, Sep. 2020.

[17] N. Mwaniki, E. Garrison, and N. Pisanti, "Fast exact string to D-texts alignments," in *Proc. 16th Int. Joint Conf. Biomed. Eng. Syst. Technol.*, 2023, pp. 70–79.

[18] G. Bernardini, A. Conte, G. Gourdel, R. Grossi, G. Loukides, N. Pisanti, S. P. Pissis, G. Punzi, L. Stougie, and M. Sweering, "Hide and mine in strings: Hardness and algorithms," in *Proc. IEEE Int. Conf. Data Mining (ICDM)*, Nov. 2020, pp. 924–929.

[19] IUPAC-IUB Commission on Biochemical Nomenclature, "Abbreviations and symbols for nucleic acids, polynucleotides, and their constituents," *Biochemistry*, vol. 9, no. 20, pp. 4022–4027, Sep. 1970.

[20] N. Pisanti, H. Soldano, M. Carpentier, and J. Pothier, "A relational extension of the notion of motifs: Application to the common 3D protein substructures searching problem," *J. Comput. Biol.*, vol. 16, no. 12, pp. 1635–1660, Dec. 2009.

[21] É. Birmelé, P. Crescenzi, R. Ferreira, R. Grossi, V. Lacroix, A. Marino, N. Pisanti, G. Sacomoto, and M. Sagot, "Efficient bubble enumeration in directed graphs," in *Proc. 19th Int. Symp. String Process. Inf. Retr. (SPIRE)*, 2012, pp. 118–129.

[22] J. M. Eizenga, A. M. Novak, J. A. Sibbesen, S. Heumos, A. Ghaffaari, G. Hickey, X. Chang, J. Seaman, R. Rounthwaite, J. Ebler, M. Rautiainen, S. Garg, B. Paten, T. Marschall, J. Sirén, and E. Garrison, "Pangenome graphs," *Annu. Rev. Genomics Hum. Genet.*, vol. 21, no. 1, pp. 139–162, Aug. 2020.

[23] B. Letcher, M. Hunt, and Z. Iqbal, "Gramtools enables multiscale variation analysis with genome graphs," *Genome Biol.*, vol. 22, no. 1, pp. 1–27, Dec. 2021.

[24] V. Carletti, P. Foggia, E. Garrison, L. Greco, P. Ritrovato, and M. Vento, "Graph-based representations for supporting genome data analysis and visualization: Opportunities and challenges," in *Proc. 12th IAPR-TC-15 Int. Workshop Graph-Based Represent. Pattern Recognit. (GbRPR)*, Jun. 2019, pp. 237–246.

[25] L. Denti, R. Rizzi, S. Beretta, G. D. Vedova, M. Previtali, and P. Bonizzoni, "ASGAL: Aligning RNA-seq data to a splicing graph to detect novel alternative splicing events," *BMC Bioinf.*, vol. 19, no. 1, pp. 444:1–444:21, Dec. 2018.

[26] J. Avila, P. Bonizzoni, S. Ciccolella, G. D. Vedova, L. Denti, D. Monti, Y. Pirola, and F. Porto, "Recgraph: Adding recombinations to sequence-to-graph alignments," *Biorxiv*, Oct. 2022. [Online]. Available: https://www.biorxiv.org/content/10.1101/2022.10.27.513962v1

[27] C. Jain, S. Misra, H. Zhang, A. Dilthey, and S. Aluru, "Accelerating sequence alignment to graphs," in *Proc. IEEE Int. Parallel Distrib. Process. Symp. (IPDPS)*, May 2019, pp. 451–461.

[28] G. Chandra and C. Jain, "Sequence to graph alignment using gap-sensitive co-linear chaining," in *Proc. 27th Annu. Int. Conf. Res. Comput. Mol. Biol. (RECOMB)*, Cham, Switzerland: Springer, 2022, pp. 58–73.

[29] M. Rautiainen and T. Marschall, "GraphAligner: Rapid and versatile sequence-to-graph alignment," *Genome Biol.*, vol. 21, no. 1, p. 253, Dec. 2020.

[30] H. Li, X. Feng, and C. Chu, "The design and construction of reference pangenome graphs with minigraph," *Genome Biol.*, vol. 21, no. 1, p. 265, Dec. 2020.

[31] N. M. Mwaniki and N. Pisanti, "Optimal sequence alignment to ED-strings," in *Proc. 18th Int. Symp. Bioinf. Res. Appl. (ISBRA)*. Cham, Switzerland: Springer, 2022, pp. 204–216.

[32] P. Ivanov, B. Bichsel, H. Mustafa, A. Kahles, G. Rätsch, and M. Vechev, "AStarix: Fast and optimal sequence-to-graph alignment," in *Proc. 24th Annu. Int. Conf. Res. Comput. Mol. Biol. (RECOMB)*. Cham, Switzerland: Springer, 2020, pp. 104–119.

[33] P. Ivanov, B. Bichsel, and M. Vechev, "Fast and optimal sequence-to-graph alignment guided by seeds," in *Proc. 26th Annu. Int. Conf. Res. Comput. Mol. Biol. (RECOMB)*. Cham, Switzerland: Springer, 2021, pp. 306–325.

[34] R. G. Koerkamp, "A*PA2: Up to 20 times faster exact global alignment," *Biorxiv*, Mar. 2024. [Online]. Available: https://www.biorxiv.org/content/10.1101/2024.03.24.586481v1

[35] M. Equi, R. Grossi, V. Mäkinen, and A. I. Tomescu, "On the complexity of string matching for graphs," in *Proc. 46th Int. Colloq. Automata*, 2019, p. 15.

[36] G. Bernardini, P. Gawrychowski, N. Pisanti, S. P. Pissis, and G. Rosone, "Even faster elastic-degenerate string matching via fast matrix multiplication," in *Proc. 46th Int. Colloq. Automata*, 2019, pp. 21:1–21:15.

[37] G. Bernardini, P. Gawrychowski, N. Pisanti, S. P. Pissis, and G. Rosone, "Elastic-degenerate string matching via fast matrix multiplication," *SIAM J. Comput.*, vol. 51, no. 3, pp. 549–576, Jun. 2022.

[38] R. Ascone, G. Bernardini, A. Conte, M. Equi, E. Gabory, R. Grossi, and N. Pisanti, "A unifying taxonomy of pattern matching in degenerate strings and founder graphs," in *Proc. 24th Int. Conf. Algorithms Bioinf. (WABI)*, 2024, pp. 14:1–14:21.

[39] E. Gabory, M. N. Mwaniki, N. Pisanti, S. P. Pissis, J. Radoszewski, M. Sweering, and W. Zuba, "Comparing elastic-degenerate strings: Algorithms, lower bounds, and applications," in *Proc. 34th Annu. Symp. Combinat. Pattern Matching (CPM)*, 2023, pp. 11:1–11:20.

[40] E. Gabory, M. N. Mwaniki, N. Pisanti, S. P. Pissis, J. Radoszewski, M. Sweering, and W. Zuba, "Pangenome comparison via ED strings," *Frontiers Bioinf.*, vol. 4, Sep. 2024, Art. no. 1397036.

[41] E. Gabory, M. N. Mwaniki, N. Pisanti, S. P. Pissis, J. Radoszewski, M. Sweering, and W. Zuba, "Elastic-degenerate string comparison," *Inf. Comput.*, vol. 304, May 2025, Art. no. 105296.

[42] A. Cislak, S. Grabowski, and J. Holub, "SOPanG: Online text searching over a pan-genome," *Bioinformatics*, vol. 34, no. 24, pp. 4290–4292, Dec. 2018.

[43] K. Aoyama, Y. Nakashima, I. Tomohiro, S. Inenaga, H. Bannai, and M. Takeda, "Faster online elastic degenerate string matching," in *Proc. 29th Annu. Symp. Combinat. Pattern Matching (CPM)*, vol. 105, Nov. 2018, p. 10.

[44] A. Cislak and S. Grabowski, "SOPanG 2: Online searching over a pan-genome without false positives," 2020, arXiv:2004.03033.

[45] G. Bernardini, N. Pisanti, S. P. Pissis, and G. Rosone, "Pattern matching on elastic-degenerate text with errors," in *Proc. 24th Int. Symp. String Process. Inf. Retr. (SPIRE)*, 2017, pp. 74–90.

[46] G. Bernardini, N. Pisanti, S. P. Pissis, and G. Rosone, "Approximate pattern matching on elastic-degenerate text," *Theor. Comput. Sci.*, vol. 812, pp. 109–122, Apr. 2020.

[47] S. P. Pissis, J. Radoszewski, and W. Zuba, "Faster approximate elastic-degenerate string matching—Part A," in *Proc. 36th Annu. Symp. Combinat. Pattern Matching (CPM)*, vol. 331, 2025, pp. 28:1–28:19.

[48] P. Gawrychowski, A. Górkiewicz, P. Marciniak, S. P. Pissis, and K. Pokorski, "Faster approximate elastic-degenerate string matching—Part B," in *Proc. 36th Annu. Symp. Combinat. Pattern Matching (CPM)*, vol. 331, Aug. 2025, pp. 29:1–29:21.

[49] G. Bernardini, E. Gabory, S. P. Pissis, L. Stougie, M. Sweering, and W. Zuba, "Elastic-degenerate string matching with 1 error or mismatch," *Theory Comput. Syst.*, vol. 68, no. 5, pp. 1442–1467, Oct. 2024.

[50] D. Gusfield, *Algorithms on Strings, Trees, and Sequences: Computer Science and Computational Biology*. Cambridge, U.K.: Cambridge Univ. Press, 1997.

[51] T. F. Smith and M. S. Waterman, "Identification of common molecular subsequences," *J. Mol. Biol.*, vol. 147, no. 1, pp. 195–197, Mar. 1981.

[52] A. Löytynoja and N. Goldman, "An algorithm for progressive multiple alignment of sequences with insertions," *Proc. Nat. Acad. Sci. USA*, vol. 102, no. 30, pp. 10557–10562, Jul. 2005.
[53] D.-F. Feng and R. F. Doolittle, "Progressive sequence alignment as a prerequisitetto correct phylogenetic trees," *J. Mol. Evol.*, vol. 25, no. 4, pp. 351-360, Aug. 1987.

[54] D. G. Higgins and P. M. Sharp, "CLUSTAL: A package for performing multiple sequence alignment on a microcomputer," *Gene*, vol. 73, no. 1, pp. 237-244, Dec. 1988.

[55] C. Lee, C. Grasso, and M. F. Sharlow, "Multiple sequence alignment using partial order graphs," *Bioinformatics*, vol. 18, no. 3, pp. 452-464, Mar. 2002.

[56] C. A. Darby, R. Gaddipati, M. C. Schatz, and B. Langmead, "Vargas: Heuristic-free alignment for assessing linear and graph read aligners," *Bioinformatics*, vol. 36, no. 12, pp. 3712-3718, Jun. 2020.

[57] J. Setubal and J. Meidanis, "*Introduction to Computational Molecular Biology*. Boston, MA, USA: PWS Publishing Company, 1997.

[58] O. Gotoh, "An improved algorithm for matching biological sequences," *J. Mol. Biol.*, vol. 162, no. 3, pp. 705-708, Dec. 1982.

[59] E. W. Myers, "AnO(ND) difference algorithm and its variations," *Algorithmica*, vol. 1, nos. 1-4, pp. 251-266, Nov. 1986.

[60] S. Marco-Sola, J. C. Moure, M. Moretó, and A. Espinosa, "Fast gap-affine pairwise alignment using the wavefront algorithm," *Bioinf.*, vol. 37, no. 4, pp. 456-463, Sep. 2020.

[61] C. Grasso and C. Lee, "Combining partial order alignment and progressive multiple sequence alignment increases alignment speed and scalability to very large alignment problems," *Bioinformatics*, vol. 20, no. 10, pp. 1546-1556, Jul. 2004.

<figure><img src="image_6.png" /></figure>

**NJAGI MOSES MWANIKI** received the B.S. degree in information science from Moi University, Kenya, in 2018, the M.S. degree in bioinformatics from Pwani University, Kenya, in 2021, and the Ph.D. degree in informatics from the University of Pisa, Italy, in 2025, funded by the EU ITN Project ALPACA. His research interests include developing methods and tools for computational pangenomics.

<figure><img src="image_7.png" /></figure>

**NADIA PISANTI** received the M.S. and Ph.D. degrees in computer science from the University of Pisa, Italy, in 1996 and 2002, respectively. From 2001 to 2005, she was a Postdoctoral Researcher at INRIA, Pisa, France; and CNRS, France. Since 2006, she has been an Assistant Professor with the Computer Science Department, University of Pisa, where she has been an Associate Professor, since 2016. She is the author of more than 100 scientific publications. She has

been in the technical or scientific PC of more than 100 conferences and seven times as the chair. She is in the steering committee of three international conferences and on the editorial board of two international journals. Her research interests include algorithmic methods and tools design for the analysis of genomic data.

...