Pangenome comparison via ED strings

Abstract

Introduction: An elastic-degenerate (ED) string is a sequence of sets of strings. It can also be seen as a directed acyclic graph whose edges are labeled by strings. The notion of ED strings was introduced as a simple alternative to variation and sequence graphs for representing a pangenome, that is, a collection of genomic sequences to be analyzed jointly or to be used as a reference.

Methods: In this study, we define notions of matching statistics of two ED strings as similarity measures between pangenomes and, consequently infer a corresponding distance measure. We then show that both measures can be computed efficiently, in both theory and practice, by employing the intersection graph of two ED strings.

Results: We also implemented our methods as a software tool for pangenome comparison and evaluated their efficiency and effectiveness using both synthetic and real datasets.

Discussion: As for efficiency, we compare the runtime of the intersection graph method against the classic product automaton construction showing that the intersection graph is faster by up to one order of magnitude. For showing effectiveness, we used real SARS-CoV-2 datasets and our matching statistics similarity measure to reproduce a well-established clade classification of SARS-CoV-2, thus demonstrating that the classification obtained by our method is in accordance with the existing one.

Keywords: SARS-CoV-2; elastic-degenerate string; intersection graph; matching statistics; pangenome comparison.

FIGURE 1

An example of an MSA (top left) and its corresponding (non-unique) ED string $T$ of length $n=7$ , cardinality $m=11$ and size $N=20$ (top right), and edge-labeled DAG for $T$ . Note that $\epsilon$ denotes the empty string. The DAG can also be viewed as an NFA with extended (multi-letter) transitions.

FIGURE 2

The two DAGs $G_1$ and $G_2$ for ED strings $T_1$ and $T_2$ . The filled black nodes are explicit states, while the orange empty nodes are implicit states.

FIGURE 3

Intersection graph $G$ for $T_1$ and $T_2$ , where $G_1$ and $G_2$ are shown at the left and on the top, respectively, to simplify the understanding of $G$ . A node $(i,j)$ in the intersection is represented by a square if both $i$ and $j$ are explicit nodes, and by a circle if only one of them is. The dashed edges of the intersection graph $G$ correspond to $\epsilon$ -transitions (namely, transitions such that no letter is read when traversed), while the solid edges correspond to the other extended transitions. A string in $\mathcal{L}(T_1)\cap \mathcal{L}(T_2)$ corresponds to a path from the starting node of $G$ to the accepting node. Here the intersection is nonempty and contains a single string $\mathtt{A}\mathtt{C}$ , which can be read on the red path.

FIGURE 4

Matching Statistics of $T_1$ and $T_2$ of our running example on their intersection graph $G$ , where, again - to simplify the understanding - we also draw $G_1$ and $G_2$ at the left and on the top, respectively. Note that this time, the pairs of implicit nodes that are reachable in a single extended transition from a pair that was previously computed are added. In the figure, there is only one such extra node, which is represented by a green open circle at the right of the graph. Here we highlight the paths that are relevant for computing the Matching Statistics array ${\mathsf{\text{MS}}}_{T_1,T_2}$ . To compute ${\mathsf{\text{MS}}}_{T_1,T_2}[1]$ , we look at the paths starting at nodes $(i,j)$ where $i$ is the explicit state one in the path-automaton of $T_1$ , and return the length of the longest label of such a path. These are the paths starting in one of the blue nodes (these are the nodes that correspond to the uppermost explicit node of $G_1$ paired with any node of $G_2$ , that is, they correspond to the uppermost dotted copy of $G_2$ ). The longest one of such paths (also drawn in blue) corresponds to the string $\mathtt{T}\mathtt{G}\mathtt{C}$ having length 3; therefore, ${\mathsf{\text{MS}}}_{T_1,T_2}[1]=3$ . For ${\mathsf{\text{MS}}}_{T_1,T_2}[2]$ we do the same but using as starting nodes those in red that correspond to the internal explicit node of $G_1$ paired with any node of $G_2$ (i.e., the nodes of the middle dotted copy of $G_2$ ). Here the longest path is drawn in red and it spells the string $\mathtt{C}\mathtt{A}$ , and therefore we set ${\mathsf{\text{MS}}}_{T_1,T_2}[2]=2$ . Computing ${\mathsf{\text{MS}}}_{T_2,T_1}$ can be performed in a dual manner on the same graph, but using as starting nodes those of the leftmost dotted copy of $G_1$ for ${\mathsf{\text{MS}}}_{T_2,T_1}[1]$ , and those of the middle dotted copy of $G_1$ for ${\mathsf{\text{MS}}}_{T_2,T_1}[2]$ .

FIGURE 5

Breakpoint Matching Statistics computation in the intersection graph $G$ of $T_1$ and $T_2$ . To compute ${\mathsf{\text{BMS}}}_{T_1,T_2}[1]$ , the candidate starting nodes of the match in $G$ are those in blue: nodes $(i,j)$ where $i$ is an explicit state of $T_1$ in the uppermost dotted copy of $G_2$ , and $j$ is either an explicit state of $T_2$ (squared blue nodes) or an implicit one (circled blue nodes). Note that $\mathtt{T}\mathtt{G}\mathtt{C}$ is the longest match that starts at the first set of $T_1$ but it does not fulfill the conditions for the Breakpoint Matching Statistics because it does not end at any breakpoint; for the same reason, $\mathtt{T}\mathtt{G}$ is also not a good candidate match. The occurrence of $\mathtt{A}\mathtt{C}$ corresponding to the blue edge starts at a blue square node; hence it is reachable from the node itself that corresponds to a pair of explicit states, and it ends at a node that is again a pair of explicit states, and hence a breakpoint for both $T_1$ and $T_2$ . There is no longer match satisfying these conditions; therefore we set ${\mathsf{\text{BMS}}}_{T_1,T_2}[1]=2$ . For ${\mathsf{\text{BMS}}}_{T_1,T_2}[2]$ we do the same but use as starting nodes those in red that correspond to the internal explicit node of $G_1$ paired with any node of $G_2$ (i.e., the nodes of the middle dotted copy of $G_2$ ). The red path spelling $\mathtt{C}$ : (i) is a prefix in $T_1[2]$ starting at an explicit node of $T_1$ ; (ii) is reachable from a square node in $G$ by spelling $\mathtt{A}$ in both strings (curved brown red edge labeled with $\mathtt{A}$ ); and (iii) ends where $T_2[2]$ does, that is, at a breakpoint. Since this is the longest such path in $G$ , we set ${\mathsf{\text{BMS}}}_{T_1,T_2}[2]=1$ . Note, for example, that the match $\mathtt{C}\mathtt{A}$ that occurs in $T_1[2]$ and inside $T_2[2]$ cannot be used for ${\mathsf{\text{BMS}}}_{T_1,T_2}[2]$ because it starts at a node that is not reachable from a pair of explicit nodes, meaning that it is not upperbounded by a breakpoint in $T_2$ . Computing ${\mathsf{\text{BMS}}}_{T_2,T_1}$ , which is of size $n_2=2$ , can be done in a dual manner on the very same graph, using as starting nodes those of the leftmost dotted copy of $G_1$ for ${\mathsf{\text{BMS}}}_{T_2,T_1}[1]=2$ (obtained by traversing an $\epsilon$ -transition and then $\mathtt{A}\mathtt{C}$ ), and those of the middle dotted copy of $G_1$ for ${\mathsf{\text{BMS}}}_{T_2,T_1}[2]=2$ ( $\mathtt{A}\mathtt{C}$ again).

FIGURE 6

SARS-CoV-2 clades pairwise similarity graph generated according to average Breakpoint Matching Statistics. The annotation (all non grey nor black graphics and text) highlights similarities with Figure 7.

FIGURE 7

Phylogeny of 3357 SARS-CoV-2 genomes samples. The figure is generated and downloaded from Nextstrain https://nextstrain.org/ncov/open/global/all-time (2024), and some annotation is added here to highlight similarities with the graph of Figure 6.