Subsequence-Based Indices for Genome Sequence Analysis

Sequence analysis is one of the core branches of bioinformatics, and it is arguably one of the most fundamental tasks due to the abundance of genomic sequences enabled by advancements in sequencing technologies. Sequence analysis methods have a huge advantage compared to “in vitro” methods: once a dataset is available, it can be instantly and easily shared with anyone, it does not deteriorate or deplete, and can be analysed repeatedly with just regular computers, without the need of expensive ad-hoc machines. This is one of the reasons why much effort is dedicated to quickly produce new and more powerful sequencers, resulting in larger datasets available to the scientific community, ranging from bacteria and viruses to humans.

The size of genomic data can be daunting, as the length of genomic sequences ranges from thousands (for some viruses, or proteins) to millions (e.g., E. Coli genome) or billions (e.g., human genome). An important trade-off is immediately evident: more complex and refined approaches can extract better-quality information from the data, but require more computational resources to be executed, and may be not be applicable to complex organisms. The research community has been advancing in two main directions: developing more sophisticated algorithms, and extending their applicability to increasingly complex data. A key task in the latter direction is optimizing the way data is represented and handled, since just storing sequences in an uncompressed format may already require tens of GB of space. In this scenario, compact indices are extremely valuable tools, that typically employ algorithmic tricks to provide a good trade-off between the size of the index, and ease of access to the data (as well as support for specific queries).

In this paper, our aim is to give an overview of two recent research results concerning subsequence-based indices in bioinformatics, which involve Roberto Grossi both in their past development and in their current investigation of future directions. We experimentally highlight how these indices do more than just represent genomic data: their clever processing of the input enables the extraction of salient information that can be used for sequence analysis application tasks.

Specifically, the first research result we review is a method, called phyBWT [25, 26], that addresses the problem of phylogenetic inference employing the Burrows-Wheeler Transform (BWT) [11]. The BWT is a text transformation with the remarkable feature of clustering together repeated sequences, a property originally intended to enhance compression, but now widely used in genome indexing algorithms [32]. Phylogenetic inference refers to the process of reconstructing the evolutionary relationships among species, or more generally, among taxa. The phyBWT methodology for phylogeny reconstruction uses the BWT of a string collection [7, 35], precisely of a group of sequences representing different taxa, to group related taxa together and to suggest evolutionary relationships. The main features of the phyBWT tool are its ability to work directly on raw sequencing data in an assembly-free scenario, and the fact that it does not rely on pairwise sequence comparisons, and thus on a distance matrix, but rather compares all the sequences simultaneously and efficiently.

The subsequent research result reviewed in this paper concerns the construction of a deterministic finite automaton (DFA) to efficiently index all maximal common subsequences (MCS) of two (or more) input strings. A common subsequence is a sequence of characters that occurs in the same order in all input strings, albeit not necessarily consecutively. An MCS is a common subsequence that is not a subsequence of any other common subsequence. The DFA, implemented as a labelled directed acyclic graph (DAG), is called McDag, and was presented in [12, 13]. Even if some previous works provide indices with better worst-case bounds [18, 28], the McDag index has been experimentally shown to be more efficient in practice for real-world genomic data, as it is significantly faster to build and typically smaller in size. Preliminary experiments (see Figure 3 from [12]) showed that the distribution of MCS lengths appears to behave differently when comparing very similar or dissimilar genomes.

Aiming to illustrate the potential of the two methods for genomic data analysis, we experimentally observe the MCS lengths distributions on genomic sequences within the same taxonomic group, to get a picture of how these distributions behave on closely related taxa versus more distant ones, using the phylogeny produced by phyBWT as a guide.

We consider a string $X=X[1]\mathrm{\dots }X[|X|]$ as a sequence of characters from a finite and ordered alphabet $\mathrm{\Sigma }$, where $X[i]\in \mathrm{\Sigma }$ denotes the character at position $i$ in $X$ and $|X|$ denotes the total number of characters in $X$. Let $X[i,j]$ denote the substring $X[i]\mathrm{\dots }X[j]$, for $1\le i\le j\le |X|$: a substring $X[i,|X|]$ (resp. $X[1,i]$) is called suffix (resp. prefix) of $X$, for any $1\le i\le |X|$. We use special characters $\{\mathrm{\#},\$\}$ as markers delimiting input strings.

We say that string $Z$ is a subsequence of $X$ if there exist indices such that $X[i_k]=Z[k]$ for . If a subsequence can be mapped on $X$ contiguously, i.e., for all , $i_k=i_1+k-1$, then $Z$ is a substring of $X$. Moreover, $Z$ is a common subsequence of strings $X$ and $Y$ if $Z$ is a subsequence of both $X$ and $Y$(see Figure 1). More specifically, let us call a pair $(i,j)$ a match when $X[i]=Y[j]$: letting be the indices such that $Y[j_k]=Z[k]$ for , then each pair $(i_k,j_k)$ is a match (as $X[i_k]=Y[j_k]$) and we say that the pairs $(i_1,j_1),\mathrm{\dots },(i_{|Z|},j_{|Z|})$ form a matching in $X$ and $Y$, whose corresponding string is $Z$. We observe that matches induce a partial order, defined as iff and , which is total if the pairs belong to the same matching; $(i,j)\le (i^{\prime },j^{\prime })$ is analogously defined. A string $Z$ is a maximal common subsequence (MCS) of $X$ and $Y$ if there is no string $W\ne Z$ that satisfies both conditions: (i) $W$ is a common subsequence of $X$ and $Y$, and (ii) $Z$ is a subsequence of $W$. The set of all strings that are maximal common subsequences is denoted by $MCS(X,Y)$.

We next introduce some graph notions. A directed graph $G=(V,E)$ consists of a set of nodes $V$ and a set of edges $E\subseteq V\times V$, where each edge $(u,v)$ is an ordered pair of nodes that specifies a direction from $u$ to $v$. Two edges $(u,v)$ and $(w,z)$ are said to be adjacent if $v=w$. A path in $G$ is a sequence of distinct edges, each adjacent to the next. If the path starts at node $s$ and ends at node $t$, it is called an $st$-path; it is a cycle when $s=t$. A DAG is a directed acyclic graph. Given a node $u$, the set $N^{+}(u)$ indicates the out-neighbor nodes $v$ such that $(u,v)\in E$, and the set $N^{-}(u)$ indicates the in-neighbor nodes $v$ such that $(v,u)\in E$. The out-degree of $u$ is $d^{+}(u)=|N^{+}(u)|$, and its in-degree is $d^{-}(u)=|N^{-}(u)|$; $u$ is a source if $d^{-}(u)=0$, and a sink if $d^{+}(u)=0$. In a labeled DAG $G=(V,E,l)$ each node $u$ is associated with a character $l(u)\in \mathrm{\Sigma }\cup \{\mathrm{\#},\$\}$. In Section 4 we also consider labeled DAGs in which each node $u$ is associated with a match $m(u)=(i_u,j_u)$.

In this section we review phyBWT, a methodology first introduced in [25] and subsequently refined in [26], for reconstructing a phylogenetic tree bypassing the standard computationally expensive steps of sequence alignment and de novo assembly. It can take as input any type of sequence data representing taxa, such as whole-genome sequences and raw sequencing reads.

The approach exploits the inherent combinatorial properties of the extended Burrows-Wheeler Transform (eBWT) [7, 35] to index and detect relevant common substrings of varying length. The common substrings then play a crucial role in building partition trees without performing pairwise comparisons between sequences. This is the primary feature of phyBWT: the tree structure is inferred by comparing all the sequences simultaneously and efficiently, without resorting to a distance matrix.

The second remarkable feature of phyBWT is that, to the best of our knowledge, it is the first approach to apply the properties of the eBWT to the idea of decomposition for phylogenetic inference. By indexing the sequences in the eBWT, we can identify maximal common substrings of varying lengths that are used to group sequences together and to partition groups of taxa based on their shared substrings.

Finally, the worst-case running time of phyBWT is $O(N\mathrm{\ell })$, where $\mathrm{\ell }$ is the number of taxa and $N$ is the total length of all the taxa sequences, using $O(N+{\mathrm{\ell }}^2)$ space.

In the following, we first provide an overview of the preliminary notions. These include the extended Burrows-Wheeler Transform employed for detecting common substrings and the positional clustering framework, which overcomes the limitation of a priori fixing the length of the common substrings. Subsequently, we delineate the tree reconstruction methodology of phyBWT according to [26].

In this section we describe the main ideas of the practically efficient compact index McDag introduced in [12]. In the present paper, for the sake of simplicity, we describe the results for two input strings, even if they have also been extended to handle an arbitrary number of strings in [13].

Let us first formalize the definition of an MCS index as follows:

Let us note that a naive construction of such an MCS index (e.g., through a trie) could potentially require exponential time and space, as the number of nodes may be proportional to the number of MCS, which can in turn be exponential in the input size. Constructing such a compact MCS index in an efficient way is therefore not trivial.

We start by giving a high-level idea of how this problem is solved by McDag in Section 4.1, and then, in Section 4.2, we focus on explaining how to efficiently compute the frequency distribution of MCS lengths from the McDag index.

In this section, we experimentally review that the BWT-based tool phyBWT can achieve benchmark-level accuracy in phylogenetic reconstruction by exploiting the common substrings among taxa. Furthermore, we provide experimental evidence suggesting that the tool for MCSs indexing could offer valuable insights for inferring evolutionary relationships.

Specifically, the plots showing the distribution of the MCS lengths reveal a notable correlation between sequences associated with taxa that are close in the phylogenetic tree.

We have reviewed two recent results that use compact string indices to naturally highlight relevant information in a genomic context. The BWT-based approach phyBWT infers evolutionary links by clustering similar substrings, and the DAG-based index McDag can be used to show the distribution of Maximal Common Subsequences, which exposes similarities among strings. We have also shown experimentally that MCS length distributions vary among closely related and more distantly related taxa, using the phylogeny generated by phyBWT as a reference. Further work in visualizing and analysing the information emerging from these indices, as well as extending the analysis to new indices, is an interesting direction to explore and may yield positive results in phylogeny and, more generally, in the analysis of genomic sequences.