BWT for String Collections

With the advent of high-throughput sequencing technologies, the amount of genomic data available in public databases has undergone an exponential increase [43]. Similar data explosion can be witnessed in other areas based on textual data, such as webpages, versioned texts (such as Wikipedia), versioned software repositories (such as Github), digital books, and many others (for more details, see [93]).

However, not only has the amount or size of textual data increased, but it has also changed character. Most importantly, nowadays most data of interest consists of string collections rather than of individual strings. For the purpose of this survey, a string collection is any multiset of strings (or sequences, two terms we use interchangeably). Early examples include the famous 1000 Genomes Project [116], in the course of which almost $2000$ human genomes were sequenced, followed by the 10,000 Genomes Project [56], the 100,000 Human Genome Project [120], the 1001 Arabidopsis Project [117], the 3,000 Rice Genomes Project (3K RGP) [115], and others. Projects of the last few years include the COVID-19 Genomics UK Consortium’s [38] and the H3Africa Consortium’s work [34].

The Burrows-Wheeler Transform (BWT) [23] is one of the most fundamental methods for string processing and compressed data storage [48, 83, 54]; the tools built upon these BWT-based compressed data structures are among the most used in bionformatics [71, 76, 69, 78]. It is therefore natural that researchers and practicioners have gone ahead and applied BWT-based methods to this new kind of data. The BWT, however, was originally defined for one string, and it is not immediately clear how to extend it to several strings.

As it turns out, several different methods have been used to generalize the BWT to string collections. Some of these methods were carefully designed and consciously chosen, while others appear to have been simply implementation choices. In fact, in many publications it is not explicitly specified how the BWT of a string collection was computed. Underlying this is the implicit assumption – incorrect, as we will see – that all methods are equivalent.

The classic method, paralleling generalized suffix trees and generalized suffix arrays (see e.g.,  [65] or [95]), is to concatenate the strings in the collection, using distinct end-of-string symbols to separate consecutive strings. In effect, this means turning the string collection into one string, at which point one can apply the classic BWT in the usual way. Following established terminology, we will refer to these end-of-string-symbols as ’dollars’; so we have a distinct dollar (${\$}_i$) for each input string. In actual implementations, of course, the same dollar-symbol is used, because otherwise the alphabet would be blown up enormously, not a viable choice. However, the dollars remain distinct conceptually; one common way to handle this is to use the position of the dollars for breaking ties. This method, concatenating the strings with distinct dollars, is what we will refer to in this paper as multidollarBWT.

Another method that is commonly used is to concatenate the input strings but drop the distinction between the dollars; this is a much simpler choice on an implementation level because one does not have to take care of distinguishing between the different dollars. In order to ensure correctness (of BWT construction, LF-mapping, backward search, etc.), it is necessary to append one more symbol at the very end of the concatenated string (often termed ’hash’ and denoted $\mathrm{\#}$), which is smaller than all other characters, including the dollar. We will call this method concatBWT.

In 2007, a mathematically clean extension of the BWT to string collections was proposed by Mantaci et al. [84], which the authors called extended BWT ($\mathrm{EBWT}$). The $\mathrm{EBWT}$ is essentially the inverse of the Gessel-Reutenauer bijection [57], known in combinatorics on words, but hitherto unknown in the research community working on textual data structures and string processing. The $\mathrm{EBWT}$ is a generalization of the BWT from one string to a multiset of primitive strings. Similarly to the BWT, it is a permutation of the characters of the strings in input and therefore has length equal to the total length of the multiset in input. It also allows backward search and has the same ’clustering effect’ as the BWT, i.e., it tends to produce long same-character runs on repetitive inputs. We refer to this method simply as $\mathrm{EBWT}$.

Applying the classic BWT to a string which terminates with a dollar is different in many ways from applying it to one without. (For example, $\mathrm{BWT}(\mathrm{𝚋𝚊𝚗𝚊𝚗𝚊}\$)=\mathrm{𝚊𝚗𝚗𝚋}\$\mathrm{𝚊𝚊}$, while $\mathrm{BWT}(\mathrm{𝚋𝚊𝚗𝚊𝚗𝚊})=\mathrm{𝚗𝚗𝚋𝚊𝚊𝚊}$.) Similarly, appending a dollar to the ends of the input strings changes the properties of the $\mathrm{EBWT}$. Again, one possibility is to append the same dollar to every input string; this is what happens, for example, reading in strings line-by-line, since each line is terminated by a newline symbol. We will call this method dollarEBWT.

Finally, if one appends different dollars to the input strings (or the same dollar but considering them as different), then we get what we call multidollarEBWT. Regarding the resulting transform, this method is similar to the multidollarBWT method (concatenate the strings with distinct dollars separating them), except that the indices of the dollars in the final transform are shifted by one.

In Table 1, we give a classification of existing BWT construction tools, geared specifically to string collections, according to the BWT variant they produce. As it is not always possible to distinguish the underlying method, our categorization is based on a combination of the information given in the accompanying papers and on reverse engineering the final output of the tools.

In this paper, we survey the five different methods sketched above, giving their exact definitions and combinatorial properties. As we will see, which method is chosen has important consequences for the resulting transform. This includes issues such as input order dependence (giving the string collection in a different order may or may not result in different outputs), dynamicity (how does adding or removing a string from the collection impact on the transform), and number of runs of the BWT. This last parameter, often denoted $r$, is central in BWT-based data structure design, since modern compressed data structures such as the $r$-index [54] require storage space linear in $r$.

We will in fact show that if the goal is to reduce the number of runs, then either the multidollarEBWT or the multidollarBWT should be chosen, because methods leading to few runs – the exact method, called optimalBWT, or heuristics like colexBWT or plusBWT – only work with these variants. Moreover, both allow adding or removing strings from the collection without having to recompute the BWT. On the other hand, if input order independence is desired, then the $\mathrm{EBWT}$ or the dollarEBWT are the right choice.

We will also see that the concatBWT, the variant where the strings are concatenated using the same dollar, has several undesirable properties.

The number $r$ of runs of the BWT appears increasingly as a parameter characterizing the input data: $r$, or equivalently, the average runlength $n/r$ of the BWT, is often used as a measure of how repetitive the data is. This is well justified in the case of one string, as many recent publications on the theory of $r$ as a repetitiveness measure attest (see [93] for a recent survey). However, for string collections, the parameter $r$ is not well defined, given that different methods on the same input collection can give different $r$’s, and even the same method on the same collection, if presented in a different order, can give very different $r$’s (up to a multiplicative factor of over $31$ on real data [30]). We argue that this parameter should be standardized, to $r_{\mathrm{opt}}$, the minimum number of runs among all possible transforms on a given string collection.

We round off the paper with a section about applications in bionformatics that make use of string collections. These are highly repetitive data sets that are the ideal object for data structures based on the BWT. We give an overview of the different areas and problems using highly repetitive sequence collections, which are at the center of today’s bioinformatics research.

Let $\mathrm{\Sigma }=\{a_1,a_2,\mathrm{\dots },a_{\sigma }\}$ be a finite ordered alphabet $\mathrm{\Sigma }$ with . A string (or text) $T=T[1..n]$ is a sequence of symbols $T[1]\mathrm{\cdots }T[n]$ over the alphabet $\mathrm{\Sigma }$. We denote by $ϵ$ the empty string and by $|T|$ the length of string $T$. For $1\le i\le n$, the substring $T[1..i]$ is called the $i$th prefix of $T$ and $T[i..n]$ the $i$th suffix of $T$. The string $T^{\mathrm{rev}}=T[n]\mathrm{\cdots }T[1]$ is the reverse string of $T$. Let denote the standard lexicographic order. The colexicographic order (sometimes also referred to as reverse lexicographic order) is defined as if and only if .

Given a string $T$, the string $T^k=TT\mathrm{\cdots }T$ is the string obtained by concatenating $k$ copies of $T$. A string $T$ is called primitive if $T=U^k$ implies $T=U$ and $k=1$. For every string $T$, there exists a unique primitive string $U$ and a unique integer $k$ such that $T=U^k$. The string $U$ is called $\mathrm{root}(T)$ and $k$ is called $\exp (T)$; thus, $T=\mathrm{root}{(T)}^{\exp (T)}$.

For a string $T$, we denote by $T^{\omega }=TT\mathrm{\cdots }$ the infinite string obtained by concatenating an infinite number of copies of $T$. The omega-order on ${\mathrm{\Sigma }}^{\ast }$ is defined as follows: $S{\prec }_{\omega }T$ if (a) , or (b) $S^{\omega }=T^{\omega }$ and (note that $S^{\omega }=T^{\omega }$ implies $\mathrm{root}(S)=\mathrm{root}(T)$). It can be verified that the $\omega$-order relation coincides with the lexicographic order if neither of the two strings is a proper prefix of the other. Otherwise, the two orders can be different. For instance, but $\mathrm{𝙲𝙶𝙰𝙲𝙲}{\prec }_{\omega }\mathrm{𝙲𝙶𝙰}$.

Given string $T$, the alphabet of $T$ is defined as $\mathrm{alph}(T)=\{T_i\mid 1\le i\le |T|\}$. A run in a string $T$ is a maximal substring $U$ such that $\mathrm{alph}(U)$ is a singleton. Given a string $T$, we denote by $\rho (T)$ the number of runs of $T$, e.g., $\rho (\mathrm{𝚋𝚊𝚊𝚊𝚋𝚋𝚌})=4$. Moreover, we denote by $r(T)$ the number of runs of the BWT of $T$ (defined below).

The string $S$ is a conjugate (or cyclic rotation, or simply rotation) of the string $T$ if $S=T[i..n]T[1..i-1]$, for some $i\in \{1,\mathrm{\dots },n\}$. In this case, $S$ is also called the $i$th conjugate of $T$, denoted ${\mathrm{conj}}_i(T)$. The conjugate array $\mathrm{CA}$ of a string $T$ of length $n$ is the permutation of $\{1,\mathrm{\dots },n\}$ such that $\mathrm{CA}[j]=i$ if ${\mathrm{conj}}_i(T)$ is the $j$th conjugate of $T$ with respect to the lexicographic order, with ties broken according to string order, i.e. if $\mathrm{CA}[j]=i$ and $\mathrm{CA}[j^{\prime }]=i^{\prime }$ for some , then either , or ${\mathrm{conj}}_i(T)={\mathrm{conj}}_{i^{\prime }}(T)$ and . If $T$ is a primitive string, a conjugate $S$ of $T$ is called Lyndon rotation of $T$, denoted as $Lynd(T)$, if it is lexicographically the smallest among all conjugates of $T$.

To define the suffix array $\mathrm{SA}$, it is common to assume that the string has an end-of-string symbol $\$$ appended at the end, such that $\$$ is smaller than all other characters and does not occur elsewhere in the string. Then, the $\mathrm{SA}$ is a permutation of the indices $\{1,2,\mathrm{\dots },n+1\}$ such that $\mathrm{SA}[j]=i$ if suffix $T[i..n+1]$ is the $j$th suffix of all suffixes in lexicographic order. Note that in this case, there can be no ties, since two distinct suffixes necessarily have different lengths. Since the string $T\$$ has exactly one occurrence of $\$$, the conjugate array and the suffix array of $T\$$ coincide.

The Burrows-Wheeler Transform ($\mathrm{BWT}$) [23] is a reversible transformation that permutes the characters of the input string. Given a string $T$, the BW-matrix of $T$ is defined as the matrix whose rows are the conjugates of $T$, given in lexicographic order. We denote by $\mathrm{BWT}(T)$ the last column of the matrix. The transformation also outputs an integer $i$, which is the position of $T$ within the matrix. It follows from the definition that if a string $S$ is a conjugate of $T$, then $\mathrm{BWT}(S)=\mathrm{BWT}(T)$. The index $i$ is fundamental for recovering $T$ from $\mathrm{BWT}(T)$, otherwise the original string can be recovered only up to conjugates. When the reconstruction of the original string is not relevant, the index is omitted.

The standard permutation of a word $S=S[1..n]$ over $\mathrm{\Sigma }$ is the permutation ${\pi }_S$ of $\{1,2,\mathrm{\dots },n\}$ such that if and only if or $S[i]=S[j]$ and . The standard permutation ${\pi }_S$ of a string $S$ is also called LF-mapping if $S=\mathrm{BWT}(T)$, for some string $T$. It is known that $T$ is a primitive word if and only if ${\pi }_S$ is a cycle of length $n$. In this case, given $S$, ${\pi }_S$ and an index $i$, the $i$th conjugate $U$ of $T$ in lexicographic order can be recovered as follows: $U[n]=S[i]$ and $U[n-j]=S[{\pi }_S^j(i)]$, for $j=1,\mathrm{\dots },n-1$. The conjugate $Lynd(T)$ can be recovered with $i=1$.

The permutation $\mathrm{BWT}(T)$ can be computed from the conjugate array $\mathrm{CA}$, since for all $j=1,\mathrm{\dots },n$, $\mathrm{BWT}(T)[j]=T[\mathrm{CA}[j]-1]$, assuming that $T[0]=T[n]$. Alternatively, if the input string $T$ terminates with an end-of-string symbol $\$$, $\mathrm{BWT}(T\$)$ can be computed from the suffix array $\mathrm{SA}$. Moreover, recovering the first row of the BW-matrix will yield $\$T$, because $\$$ is smaller than all other symbols, and since $\$$ occurs exactly once, unique recovery is ensured also in this case.

Appending an end-of-string symbol to the string $T$ has traditionally allowed $\mathrm{BWT}(T\$)$ to be computed using linear-time algorithms for suffix array construction. However, it is possible to compute $\mathrm{BWT}(T)$ in linear time without appending any end-of-string symbol $\$$ to $T$. This can be done by computing the suffix array of $U=Lynd(\mathrm{root}(T))$ [58] and then obtaining $\mathrm{BWT}(T)=\mathrm{BWT}(U^k)$, where $k=\exp (T)$, from $\mathrm{BWT}(U)$, applying a simple relationship between the two transforms: $\mathrm{BWT}(U^k)$ can be obtained by concatenating $k$ copies of each character of $\mathrm{BWT}(U)$ [85]. Alternatively, it is also possible to compute $\mathrm{BWT}(T)$ directly in linear time with the $\mathrm{𝚌𝚊𝚒𝚜}$ algorithm, which is based on induced sorting [18].

Now let $ℳ=\{T_1,T_2,\mathrm{\dots },T_m\}$ be a multiset of strings (often referred to as a string collection), with $|T_d|=n_d$ for $1\le d\le m$. We denote by $|ℳ|=m$ the number of strings in $ℳ$ and by $\Vert ℳ\Vert ={\sum }_{d=1}^m{n}_d$ their total length. Further, we will denote by ${ℳ}_{\$}=\{T_1{\$}_1,T_2{\$}_2,\mathrm{\dots },T_m{\$}_m\}$ the set of strings with end-of-string markers appended to all strings in $ℳ$, of total length $m+{\sum }_{d=1}^m{n}_d$.

The generalized conjugate array $\mathrm{GCA}$ of the collection $ℳ=\{T_1,T_2,\mathrm{\dots },T_m\}$ is an array of total length $\Vert ℳ\Vert$, where $GCA[j]=(d,i)$ if $con{j}_i(T_d)$ is the $j$th string in the $\omega$-sorted list of the conjugates of the strings in $ℳ$, with ties broken first w.r.t. the index of the string (in case of identical strings), and then w.r.t. the position in the string (in case of identical conjugates of the same string).

The generalized suffix array $\mathrm{GSA}$ of the collection ${ℳ}_{\$}$ is an array of length $\Vert {ℳ}_{\$}\Vert =m+{\sum }_{d=1}^m{n}_d$, containing pairs of integers $(d,i)$, corresponding to the lexicographically sorted suffixes of the strings in ${ℳ}_{\$}$. In particular, $\mathrm{GSA}[j]=(d,i)$ is the pair corresponding to the $j$th smallest suffix of the strings in ${ℳ}_{\$}$, with ties broken according to .

For string collections, several variants of the BWT have been proposed. Each of these variants aims to extend the features of the BWT for single strings to collections of strings. In this section, we review the most significant variants, which form the foundation of the most widely used tools for indexing and compressing string collections based on the BWT. We use the classification introduced in [32, 33], and in addition group these variants into two families, depending on the strategy used to process the strings of the collection.²²2In difference to [32, 33], here we include a new separator-based BWT variant, the multidollarEBWT, which was previously subsumed under multidollarBWT. The first family includes transformations that explicitly or implicitly apply the Extended BWT ($\mathrm{EBWT}$) introduced in [84] to the multiset of strings, to which end-of-string symbols may or may not have been appended. The second family, on the other hand, consists of those transformations that involve concatenating the strings of the collection, and then applying the traditional $\mathrm{BWT}$ to the concatenated string. We will see that these BWT variants, while preserving the general functionality of the $\mathrm{BWT}$, can produce significantly different outputs.

Let us denote by $ℳ=\{T_1,T_2,\mathrm{\dots },T_m\}$ the input string collection, i.e., a multiset of $m$ strings, with $|T_d|=n_d$ for $1\le d\le m$. Even though the input is a multiset of strings, in real applications, this input is always given in some order. As we will see, the input order plays an important role in the resulting transform for several of the BWT variants we will review. For this reason, in slight abuse of notation, whenever the definition depends on the input order, we will treat $ℳ$ in the following as an ordered collection and will write accordingly $ℳ=(T_1,T_2,\mathrm{\dots },T_m)$. When it does not play a role, in the sense that the output is always the same whatever the input order, in those cases we will stick to the multiset notation.

In Table 1 we list, for each BWT variant, the tools that compute it (up to renaming of dollars).

In this section, we study combinatorial properties of the five BWT variants introduced in Section 3. We give an overview in Table 4. Full proofs for Sections 4.1 and 4.2 can be found in the Supplementary Material of [33]. Section 4.3 is new, and Section 4.4 is a generalization of the respective results in [32].

First notice that all BWT variants except the $\mathrm{EBWT}$ append an end-of-string symbol (resp. separator-symbol) at the end of each input string. For this reason, we will in the following refer to these BWT variants as separator-based BWT variants.

Recall that $L\widehat{=}{L}^{\prime }$ if the $L$ and $L^{\prime }$ are equal up to renaming the dollar symbols. In Section 4, we showed that all separator-based BWT variants can be simulated by mdolEBWT or by mdolBWT (Thm. 17): if $L$ is the output of one of these variants, then the input collection can be permuted in such a way that mdolEBWT and mdolBWT will produce a transform which differs from $L$ only in the denomination of the dollar signs. This justifies the following definition.

Given a string collection $ℳ$, let ${𝔖}_{ℳ}$ denote the string family obtained by applying mdolEBWT to all possible orders of the strings in $ℳ$. Formally, for $ℳ=(T_1,\mathrm{\dots },T_m)$,

This string family coincides with the one introduced in [14] (albeit there given with a different definition).

In this section, we focus on reducing the number of runs of the BWT, commonly referred to as $r$. This is a fundamental parameter, since the BWT is primarily used for compressed text indexes, such as the FM-index [48], the $r$-index [54], or the extended $r$-index [19, 20], whose storage space depends on $r$. Moreover, the value $n/r$, the average runlength of the BWT, is increasingly being used as a parameter describing the input data. This is because it is well known that the more repetitive the input data, the fewer runs the BWT tends to have.

Therefore, we are interested in the following:

As we saw in Section 4.2, any differences between two transforms $L$ and $L^{\prime }$ on the same collection $ℳ$, in whatever order, necessarily lie within interesting intervals. This ignores the naming of the separators, which is justified, given that no tool outputs the transform with different dollar-symbols, even if it computes the mdolBWT or mdolEBWT, which means that internally it handles the dollar-symbols as different symbols.⁵⁵5Some tools are able to output the indices corresponding to the dollars, such as BCR and BEETL [6, 7, 8], but this is done in a separate file, or the indices can be extracted from the DA or the GSA, such as eGAP [46, 45]. Thus, for the sake of compression, it is appropriate to treat the output independently of the names of the dollar-symbols.

Bentley et al. [9] gave a linear-time algorithm for minimizing the number of runs but they gave no implementation. An implementation called optimalBWT [99] was presented by Cenzato et al. [30] which computes an optimal transform $L$. To understand how the algorithm works, let us consider an interesting interval $[b,e]$. It is clear that within this interval, the number of runs can be minimized if all occurrences of same characters are grouped together. Therefore, if there are $k$ distinct characters in $L[b..e]$, then $k$ is the minimum number of runs achievable in this interval. Now let’s say we have grouped all characters into $k$ runs and $L[b..e]=a_1^{r_1}{a}_2^{r_2}\mathrm{\cdots }{a}_k^{r_k}$. If the previous character $L[b-1]\ne a_1$ then it may be possible to reduce the number of runs by swapping some other run, say $a_i^{r_i}$ with the first run $a_1^{r_1}$, namely if $L[b-1]=a_i$. Similarly, one may be able to reduce the number of runs by swapping some other run with the last one $a_k^{r_k}$. The problem is that this choice may not be unique, and when there are consecutive interesting intervals then it may not be possible to decide until these interesting intervals end. Bentley et al. give an algorithm for making these choices, modeling the problem as what they refer to as a tuple ordering problem, which can be turned into a shortest path problem on a DAG and solved optimally in linear time.

Now consider our toy example $ℳ=(\text{CTGA},\text{TG},\text{GTCC},\text{TCA},\text{CGACC},\text{CGA})$. In Table 5, we give different transforms from ${𝔖}_{ℳ}$, with the number of runs between $19$ (top row, input order as given above), and $14$ (bottom row), which is the minimum $r_{\mathrm{opt}}$. In the middle rows we give two heuristics, which we explain next.

Recall that the colex order (also referred to as reverse lexicographic order (rlo)) is defined as if . It is easy to see that if the input strings are in colex order, then every interesting interval will have the minimum number of runs. On our example, the colex results in $18$ runs, one less than the input order. In particular, several existing tools (BCR [6, 7], ropeBWT2 [107, 74], ropeBWT3 [108, 75]) have the option to output the BWT with respect to the colex order. This was already shown in [39] to significantly reduce the number of runs.

Several other heuristics were introduced in [14], all based on the idea of grouping all occurrences in interesting intervals into one single run and possibly swapping the runs within the interesting interval. The one that performed best was plusBWT, which during the construction greedily puts as first run the character immediately before the interesting interval and as last run the character immediately after, if these characters are present. In Table 7, we report experimental results, over benchmarks that are described in Table 6, that show that both of these heuristics get very close to the optBWT. Interestingly, even randomly permuting the runs within interesting intervals results in a slight improvement over the colex order. For details, see [14, 6].

What is also remarkable is the improvement that can be achieved by the optimalBWT as compared to the input order, namely up to a multiplicative factor of over $31$, on real biological data [30, 14], see Table 7. In other words, a resulting data structure built on top of the BWT would be far smaller if the optBWT was used, rather than just any BWT transform. Indeed, in [31], first experimental results on the size of the resulting $r$-index on real datasets are presented.

Another interesting aspect is the following: in order to reduce the number of runs, a necessary step is grouping all equal characters together within the interesting intervals. The larger the interesting intervals, the fewer the fraction of permutations of the string collection that do this. This implies that the probability that some arbitrary input order – as given by the order in which the collection comes, e.g. as downloaded from a database – will group these characters together is very small. This explains why the potential gain from applying the optBWT (or, indeed, one of the heuristics given above) is highest on very large sets of very similar short strings: many suffixes will be shared by many input strings, which produces long interesting intervals. (See [32, 33] for a more detailed analysis.)

Now consider the following example, which shows that with the concatBWT, it can be impossible to minimize the number of runs of the BWT. In other words, on some data sets, whatever the input order, $\text{concatBWT}(ℳ)$ will always have more runs than necessary.

Recent advances in DNA and RNA sequencing technologies led to the explosion of the amount of in silico genomic data to be assembled, investigated, stored and compared. In this final section we will overview the applications in bioinformatics that benefit from the efficiency of BWT-based indexing of string collections.

Both DNA and RNA sequences can be seen as strings over an alphabet of size four: the set of nucleotides that compose them. Furthermore, as we will see below, in all the applications we will mention, sequences are highly repetitive and the different strings of the collection encode redundant data: an ideal target for BWT based tools that are specifically designed for string collections.

When the string collection consists of raw sequencing data, that is, a large amount of fragments that directly result from the sequencing process, the most requested tasks can be:

Assembling.

Size and accuracy details have changed across the time and also change from a technology to another, but all sequencing machines are not able to process molecules of arbitrary size. As a consequence, the typical procedure is to: first, create a redundant number of copies of the sequence while still in vitro; second, cut (this can be done using specific enzymes) the sequences into small enough fragments randomly enough to ensure that distinct copies are cut at different sites; third, perform the actual sequencing that turn the sequence into in silico data; fourth, exploiting the redundancy of the original copies and hence fragments overlap, build a layout for the fragments that enables to reconstruct the original sequence. This last task is named Fragment Assembly, and its efficient solution played a crucial role in the historical human genome project. For decades the problem of fragment assembly has been the Holy Grail of algorithmic tasks for bioinformatics, having to deal with millions of DNA fragments, each of a few hundreds of length, that had to be jointly analysed, looking for local similarities among them [111].

Mapping.

When the newly sequenced genome belongs to species for which a fully assembled reference genome already exists, then one does not need to perform a de novo assembly out of the fragments, as the reference can be exploited. In this case, the strings of the collection are, in general, not as many as in the case of the assembly task, because redundancy is not required as much as there. However, a huge amount of strings must be mapped onto the reference genome, dealing with sequencing errors, repetitions inside the genome that make the mapping loci ambiguous, as well as the physiological differences among the query genome and the reference [71, 76].

Haplotype Matching.

A Haplotype is a set of genetic markers or DNA variations that are inherited together from a single parent. Haplotype Matching asks for finding long matches between sequences within a given large collection of aligned genetic sequences (that are, indeed, the haplotypes). The task is relevant because these long matches are candidates to be regions that are identical by descent (IBD) from a common ancestor. Even for this task BWT based methods have been designed: the positional BWT [44, 113, 114, 110] for linear strings, and more recently the GBWT for Pangenome Graphs [109].

Assembly free (and reference free) variant calling.

Variant calling is the fundamental task that has to be performed with a genome, and it consists of detecting individual traits and genetic diversity of the query genome to discover genetic based diseases, perform personalised medicine, as well as to understand evolution (when variants are annotated at population scale). When a reference genome is not available, this task is performed directly on raw data without a preliminary assembling phase, and hence this requires the investigation and comparison of (pairs of) large collections of strings [104, 102, 103, 70].

When, instead, the string collections are multiple genomes, or several RNA isoforms, or belong to possibly many different genomes, we have a (possibly large) set of strings that are as long as billions of letters. In these cases, the downstream analyses tasks are different, and they include:

Comparative genomics and variant calling.

Here two or more (whole) genomes must be compared, one of them possibly being the reference. Like for the above mentioned tasks, criticalities arise because of the highly repetitive nature of genomes, and by the different nature of variants that one needs to detect: the so called point mutations or the SNPs that involve single letters, or the short INDELs that are the insertion or the deletion of a short fragment, up to the fragment level mutations: the deletions or insertion or duplication of a possibly large DNA fragment, or the translocation (that is, its position in the genome changes), or the copy number variations (fragments that are commonly present in multiple copies show a different number of them in different individuals). Moreover, when a DNA mutation event involves a whole fragment, there is in general a $50\%$ probability that a copy is inverted, that is it occurs in the reversed direction. The comparison tasks are very complex and cannot be made by means of global alignments due to the fragment level mutations. Hence efficient indices of the strings are a basic toolkit for the various possible algorithmic strategies [20, 44].

Transcriptomics.

With the so called new generation sequencing technologies, also RNA-Seq became a common task, that is to move to in silico also RNA the sequences. In particular, mRNA (messenger RNA) isoforms are very interesting to investigate as they are RNA fragments that encode the information that leads a coding part of the DNA (that is, a gene) to the synthesis of a protein. The transcription of this coding part of DNA into RNA, called gene expression, is a very interesting topic for molecular biologists and genetists because the same gene can express in different ways according to the cell it is (say, in the liver rather than in the skin), or according to external conditions (say, in a skin cell whether it is cold or it is being burned). The set of RNA isoforms that are in a certain cell under specific conditions is called the transcriptome, and this is a sort of a picture of gene expression. Investigating or comparing transcriptome(s) requires counting the abundance of each kind, being each of them a different concatenation of a subset of the fragments of the gene: again several strings, and again highly redundant, and hence a perfect target for a string collection index based on the BWT [119, 18, 20].

Metagenomics.

Metagenomics is the (comparative) study of genetic material found in environmental samples (say, into the soil, or in sea water, or in an animal’s gut), with the purpose of detecting and investigating (the distribution of) genomes of different species of microorganisms such as bacteria, viruses, and fungi. The ultimate goals of metagenomic is the study of microbial diversity, and the interactions: both mutual and with the ecosystem it belongs to. This requires various tasks that include profiling, assembling, comparing and classifying data consisting of a huge amount of fragments belonging to different individuals and different species. As such, the redundancy will be less than in the other applications; however, indexing these string collections with BWT is still fundamental to make efficient compression and analyses [67, 64, 14, 62, 63, 40].

Pangenomics.

Using a single genome as the reference clearly introduces a bias for the discovery of genomic variations. Therefore, in the literature there came a new challenge that simultaneously deals with, and exploits, the widespread availability of sequencing data: using a pangenome rather than a genome. In [118], a pangenome was defined as “any collection of genomic sequences to be analyzed jointly or to be used as a reference”. In contrast to a linear reference, a pangenome reference aims to compactly represent the variation within a population by encoding the commonalities and differences among the underlying sequences. In the literature there are many suggestions for pangenome representations [77, 4, 28, 106, 55, 100, 79, 13, 12, 92, 26, 72] whose construction may benefit from efficient indexing of data as a starting step [81, 46, 53, 89, 98, 110, 18, 109, 35]. Some of these representations support efficient and accurate analyses tools such as comparison [1, 2, 51, 52, 50] or alignments [90, 11, 91, 10, 3, 36], that could be sped up with string collections indices. However, by now none of these graph structures has been acknowledged as the standard yet. The computational challenges posed by pangenomes often result in a trade-off between the efficiency and accuracy of the methods and the information content of the chosen representation. In this view, a valid option for certain tasks turns out to be to just index and investigate the collection of genomes, disregarding loci and alignment information and rather choosing speed and low memory usage [21, 68, 22, 44].

In this paper, we studied different methods for constructing the BWT of string collections (as opposed to individual strings). We showed that the methods in use are non-equivalent, in the sense that they produce different transforms, which differ not only in the denomination of separator-symbols (“dollars”). We analyzed precisely where and how these differences occur (“interesting intervals”). The differences between the transforms extend to the number of runs $r$ of the BWT, a central parameter in compressed text indexing.

We showed that several of these methods are input-order dependent, and that this can be exploited to significantly reduce the number of runs: by up to a factor of over $31$, on real biological data. This is the multiplicative gain of the optBWT as compared to computing the transform on the collection in the order it comes in. We also saw that several simple heuristics come close to this optimum.

On the downside, we saw that one of the most common methods, which we refer to as concatBWT (concatenating the input strings and separating them with the same dollar), has several undesirable properties. First, it cannot produce all possible transforms, and can therefore in general not be used as a method for run minimization. Second, it is not robust with respect to addition or removal of strings from the collection (dynamicity). Finally, it is input-order dependent but in a fairly opaque manner (its order being in essence the BWT of the input order), which makes it difficult to control the output in any way.

Based on our findings, we make the following suggestions:

1. 1.

   We advise to exercise care when using the method we term concatBWT, due to the issues listed above. It is a good method because it is simple and efficient, but the user (or programmer) should be aware of its shortcomings.

2. 2.

   We believe that the definition of the number of runs of a string collection should be standardized to $r_{\mathrm{opt}}$, the minimum number of runs among all possible separator-based transforms. This number is also fairly easy to approximate with heuristics, the simplest of which is colexBWT.