{'fraction_non_alphanumeric': 0.04425470761037259, 'fraction_numerical': 0.011222182302343767, 'mean_word_length': 4.1259885169537425, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FSG: Fast String Graph Construction for De Novo Assembly of Reads Data

The string graph for a collection of next-generation reads is a lossless data representation that is fundamental for de novo assemblers based on the overlap-layout-consensus paradigm. In this paper, we explore a novel approach to compute the string graph, based on the FM-index and Burrows-Wheeler Transform. We describe a simple algorithm that uses only the FM-index representation of the collection of reads to construct the string graph, without accessing the input reads. Our algorithm has been integrated into the SGA assembler as a standalone module to construct the string graph. The new integrated assembler has been assessed on a standard benchmark, showing that FSG is significantly faster than SGA while maintaining a moderate use of main memory, and showing practical advantages in running FSG on multiple threads.


Introduction
De novo sequence assembly continues to be one of the most fundamental problems in Bioinformatics. Most of the available assemblers (Simpson et al., 2009;Peng et al., 2012;Bankevich et al., 2012;Chikhi and Rizk, 2013;Salikhov et al., 2014;Chikhi et al., 2015) are based on the notions of de Bruijn graphs and k-mers (short k-long substrings of input data). Currently, biological data are produced by different Next-Generation Sequencing (NGS) technologies which routinely and cheaply produce a large number of reads whose length varies according to the specific technology. For example, reads obtained by Illumina technology (which is the most used) usually have length between 50 and 150 bases (Salzberg et al., 2012).
To analyze datasets coming from different technologies, hence with a large variation of read lengths, an approach based on same-length strings is likely to be limiting, as witnessed by the recent introduction of variable order de Bruijn graphs (Boucher et al., 2015). The string graph (Myers, 2005) is an alternative approach that does not need to break the reads into k-mers (as in the de Bruijn graphs), with the advantage of immediately distinguishing repeats longer than k but contained in a read-when using a de Bruijn graph, those repeats are resolved only at later stages. The string graph is the main data representation used by assemblers based on the overlap-layout-consensus paradigm. Indeed, in a string graph the vertices are the input reads and the arcs corresponds to overlapping reads, with the property that contigs are paths of the string graph.
Even without repetitions, analyzing only k-mers instead of the full-length reads can result in some information loss, since two bases that are k + 1 positions apart can belong to the same read, but are certainly not part of the same k-mer. Indeed, differently from de Brujin graphs, any path of a string graph is a valid assembly of reads. String graphs are more computationally intensive to compute (Simpson and Durbin, 2012), justifying our search for faster algorithms. The most widely used string graph assembler is SGA (Simpson and Durbin, 2010), which first constructs the BWT (Burrows and Wheeler, 1994) and the FM-index of a set of reads, and then uses those data structures to efficiently compute the arcs of the string graph (connecting overlapping reads). Another string graph assembler is Fermi (Li, 2012) which implements a variant of the original SGA algorithm (Simpson and Durbin, 2010) that is tailored for SNPs and variant calling.
Several recent papers face the problem of designing efficient algorithmic strategies or data structures for building string graphs. Among those works, we can find a string graph assembler (Ben-Bassat and Chor, 2014), based on a careful use of hashing and Bloom filters, with performance comparable with the first SGA implementation (Simpson and Durbin, 2010). Another important alternative approach to SGA is Readjoiner (Gonnella and Kurtz, 2012) which is based on an efficient computation of a subset of exact suffix-prefix matches, and by subsequent rounds of suffix sorting, scanning, and filtering, obtains the non-redundant arcs of the graph.
All currently available assemblers based on string graphs (such as SGA) need to both (1) query an indexing data structures (such as an FM-index), and (2) access the original reads to detect prefix-suffix overlaps between the elements. Since the self-indexing data structures, such as FM-index, represent the whole information of the original dataset, an interesting problem is to design efficient algorithms for the construction of string graphs that only require to keep the index, while discarding the original reads.
Improvements in this direction have both theoretical and practical motivations. Indeed, detecting prefix-suffix overlaps only by analyzing the (compressed) index is an almost unexplored problem.
The information contained in the indexing data structure can be analyzed with different and almost orthogonal approaches. A natural and straightforward goal that we have explored previously is to minimize the amount of data maintained in RAM (Bonizzoni et al., 2016a). In this paper we will focus instead on reducing the running time, by introducing a method that is able to build the whole string graph, via a limited number of sequential scans of the index. This property leads to the design of an algorithm that can exploit some features of modern processors. Moreover, since our new algorithm computes the string graph, we have a memory conscious and time efficient tool that may be directly integrated in a pipeline for assembling DNA reads.
We propose a new algorithm, called FSG, to compute the string graph of a set R of reads with O(nm m ) worst-case time complexityn is the number of reads in R and m is the maximum read length. To the best of our knowledge, it is the first algorithm that computes a string graph using only the FM-index of the input reads. The vast literature on BWT and FM-index hints that this approach is amenable to further research. Our algorithm is based on a characterization of the string graph given in (Bonizzoni et al., 2016a), but we follow a completely different approach.
An important observation is that, to compute the arcs outgoing from each read r, SGA queries the FM-index for each character of r. While this approach works in linear, i.e., O(nm), time, it can perform several redundant queries, most notably when the reads share common suffixes (a very common case). Our algorithm queries the FM-index in a specific order, so that each distinct string is processed only once, while SGA might process more than once each repeated string.
It is important to notice that our novel algorithm uses a characterization of a string graph  that is different, but equivalent, to the one in (Myers, 2005). We have implemented FSG and integrated it with the SGA assembler, by replacing the procedure to construct the string graph. Our implementation follows the SGA guidelines, that is we use SGA's read correction step before computing the overlaps without allowing mismatches (which is also SGA's default choice). Indeed, the guidelines to reproduce the assembly of the dataset NA12878 included in the SGA software package set the error-rate parameter to 0, the default value. Therefore, it is fair to compare the performances of the two tools. Also the assembly phases of SGA can be applied without any modification. These facts guarantees that the assemblies produced by our approach and SGA are the same, except for the unusual case when two reads have two different overlaps. In that case, SGA considers only the longer overlap, while we retain all overlaps. While it is trivial to modify our approach to guarantee the the assembly is the same, we have decided that considering all overlaps is more informative. We want to point out that the FSG algorithm is relatively simple and could be useful also for different assembly strategies.
We have compared FSG with SGA -a finely tuned implementation that has performed very nicely in the latest Assemblathon competition (Bradnam et al., 2013) -where we have used the latter's default parameter (that is, we compute overlaps without errors). Our experimental evaluation on a standard benchmark dataset shows that our approach is 2.3-4.8 times faster than SGA in terms of wall clock time (1.9-3 times in terms of user time). requiring only 2.2 times more memory than SGA.

Preliminaries
We briefly recall some standard definitions that will be used in the following. Let Σ be a constant-sized alphabet and let S be a string over Σ. We denote by S [i] the i-th symbol of S , by ł = |S | the length of S , and by S [i : j] the substring respectively. Given two strings (S i , S j ), we say that S i overlaps S j iff a nonempty suffix β of S i is also a prefix of S j , that is S i = αβ and S j = βγ. In that case we say that that β is the overlap of S i and S j , denoted as ov i, j , that γ is the right extension of S i with S j , denoted as rx i, j , and α is the left extension of S j with S i , denoted as lx i, j .
In this paper we consider a set R of n strings over Σ that are terminated by the sentinel $, which is the smallest character. To simplify the exposition, we will assume that all input strings have exactly m characters, excluding the $. The overlap graph of a set R of strings is the directed graph G O = (R, A) whose vertices are the strings in R. For each three strings α, β, and γ such that r i = αβ and r j = βγ are two strings, there is the arc (r i , r j ) ∈ A. In this case β is called the overlap of the arc.
Observe that the notion of overlap graph originally given by (Myers, 2005) is defined by labeling with the right extension rx i, j = γ the arc (r i , r j ) ∈ A. The assembly string related to (r i , r j ) is given by r i γ. More in general, given a path π = r 1 , r 2 , · · · , r n in the overlap graph, the assembly string is r 1 rx 1,2 · · · , rx n−1,n .
The notion of a string graph derives from the observation that in a overlap graph the label of an arc (r, s) may be equal to the assembly string of a path r, . . . , s : in this case the arc (r, s) is called called redundant and it can be removed from the overlap graph without loss of information, since all paths resulting in a valid assembly are still in the graph, even after the removal of such redundant arcs (r, s). In (Myers, 2005) redundant arcs are those arcs (r, s) labeled by γ, for γ containing as prefix the label of an arc (r, t). In (Bonizzoni et al., 2016b) we state an equivalent characterization of string graphs (given below) which is a direct consequence of the fact that an arc (r i , r j ) is labeled by the left extension α and its assembly is αr j . An arc e 1 = (r i , r j ) of G O labeled by α is transitive (or reducible) if there exists another arc e 2 = (r k , r j ) labeled by δ where δ is a suffix of α. Therefore, we say that e 1 is non-transitive (or irreducible) if no such arc e 2 exists. The string graph of R is obtained from G O by removing all reducible arcs. This definition allows to use directly the FM-index to compute the labels of the overlap graph since the labels are obtained by backward extensions on the index. The Generalized Suffix Array (GSA) (Shi, 1996) of R is the n-long array SA where each element SA[i] is equal to (k, j) if and only if the k-long suffix r j [|r j | − k + 1 :] of the string r j is the i-th smallest element in the lexicographic ordered set of all suffixes of the strings in R. The Longest Common Prefix (LCP) array of R, is the n-long array L such that L[i] is equal to the length of the longest prefix shared by the the k i -suffix of r j i and the k i−1 -suffix of is the symbol that precedes the k-long suffix of a string r j where such suffix is the i-th smallest suffix in the ordering given by SA.
The i-th smallest (in lexicographic order) suffix is denoted by . Given a string ω, all suffixes of R whose prefix is ω appear consecutively in LS . We call ω-interval (Bauer et al., 2013) the maximal interval [b, e] such that ω is a prefix of LS [i] for each i, b ≤ i ≤ e. We denote the ω-interval by q(ω). The width e − b + 1 of the ωinterval is equal to the number of occurrences of ω in some read of R. Since LS , the BWT B and SA are closely related, we also say that [b, e] is a ω-interval on all those arrays. Given a ω-interval and a character c, the backward c-extension of the ω-interval is the (possibly empty) cω-interval. We recall that the FMindex (Ferragina and Manzini, 2005) is essentially made of the two functions C and Occ, where C(c), with c a character, is the number of occurrences in B of characters that are alphabetically smaller than c, while Occ(c, i) is the number of occurrences of c in the prefix B[: i − 1]. Given a string α and a character c, the backward c-extension of q(ω) = [b, e] is q(cω) = [C(c) + Occ(c, b) + 1, C(c) + Occ(c, e + 1)] (Ferragina and Manzini, 2005).

The Algorithm
Our algorithm is based on two steps: the first is to compute the overlap graph, the second is to remove all transitive arcs. Given a string ω and R a set of strings (reads), let R S (ω) and R P (ω) be respectively the subset of R with suffix (resp. prefix) ω. As usual in string graph construction algorithms, we will assume that the set R is substring free, that is no string is a substring of another one. A fundamental observation is that the list of all nonempty overlaps β is a compact representation of the overlap graph, since all pairs in R S (β) × R P (β) are arcs of the overlap graph. Moreover, each arc (r i = αβ, r j = βγ) of the overlap graph can be represented by the triple (α, β, γ).
Our approach to compute all overlaps between pairs of strings is based on the notion of potential overlap, which is a nonempty string β * ∈ Σ + that is a proper suffix of an input string r i = αβ * (α ǫ) and such that there exists an input string r j = γβ * δ, with δ ǫ, containing β * as a substring, but not a suffix.
A simple relation between overlaps and potential overlaps is given in Proposition 1, which is a direct consequence of the definition of potential overlap.
Proposition 1. Let β be an overlap. Then all suffixes of β are potential overlaps.
We can now briefly sketch our algorithm which consists of two main parts. The first part computes all potential overlaps, starting from those of length 1 and extending the potential overlaps by adding a new leading character. Each potential overlap is also checked to determine whether it is also an actual overlap, to compute the set of all overlaps.
The second part of our algorithm, that is to detect all transitive arcs, starts from the sets ARC(α = ǫ, ǫβ, X = R p (β)) (a set for each one of the overlaps β) that can be immediately obtained from the overlaps β computed in the first step, where in general a set ARC(α, αβ, X) consists of arcs with overlap β, a left extension that has α as a suffix, and are incoming into a read in X. Observe that ARC(α = ǫ, αβ, X = R p (β)) is the set of the arcs of the overlap graph having overlap β. During the second part of the algorithm, the transitive arcs of the overlap graph are removed by iteratively computing arc-sets of increasing extension length ℓ (starting from α = ǫ) by adding a leading character to the extensions α of length ℓ − 1 (of the sets computed at same previous iteration) and deleting reads from X. All the computed sets ARC(α, αβ, X), where αβ is a read (that is, α is the complete left extension of the arcs), contains only irreducible arcs. The following definition formalizes the previously discussed notions.
Definition 1. Let α be a string, let β be a nonempty string, and let X be a subset of R P (β). The arc-set ARC(α, αβ, X) is the set {(r 1 , r 2 ) : αβ is a suffix of r 1 , β is a prefix of r 2 , and r 1 ∈ R, r 2 ∈ X}. The strings α and β are called the extension and the overlap of the arc-set. The set X is called the destination set of the arc-set.
An arc-set ARC(α, αβ, X) is terminal if there exists r ∈ R s.t. r = αβ, while an arc-set is basic if α = ǫ (that is the empty string). Since the arc-set ARC(α, αβ, X) is uniquely determined by strings α, αβ, and X, the triple (α, αβ, X) wil be used in our algorithm to encode the arc-set ARC(α, αβ, X).
Observe that the string α in Definition 1 is a suffix of the left extension (label) of the arcs in ARC(α, αβ, X). When the arc-set is terminal then the extension α of the arc-set is also the label of its arcs.
Moreover, the arc-set ARC(α, αβ, X) is correct if X ⊇ {r 2 ∈ R P (β) : r 1 ∈ R S (αβ) and (r 1 , r 2 ) is irreducible}, that is all irreducible arcs whose overlap is β and whose (left) extension has a suffix α are incoming in a read of X.
Observe that our algorithm outputs only correct arc-sets (hence all irreducible arcs are preserved). Moreover terminal arc-sets, computed by our algorithm, only contain irreducible arcs (see Lemma 5), hence all transitive arcs are removed.
Our algorithm is based on an extension of the definition of transitive arc to arc-sets. We present Algorithm 1 which computes the overlap graph of a set of input strings (the overlap graph is represented by the set of all overlaps), and Algorithm 2 which receives the overlaps computed by Algorithm 1 and outputs the string graph. In our description we assume that, given a string ω, we can compute in constant time (1) the number suff(ω) of input strings whose suffix is ω, (2) the number pref(ω) of input strings whose prefix is ω, (3) the number substr(ω) of occurrences of ω in the input strings. Moreover, we assume to be able to list the set listpref(ω) of input strings with prefix ω in O(|listpref(ω)|) time.
In the next section we will describe how to compute such a data structure.
We recall that Algorithm 1 exploits Proposition 1 to compute all overlaps. More precisely, given a k-long potential overlap β * , the (k + 1)-long string cβ * , for c ∈ Σ, is a potential overlap if and only if suff(cβ * ) > 0 and substr(cβ * ) > suff(cβ * ). We construct incrementally all potential overlaps, first by determining if each string β * consisting of a single character is a potential overlap. Then, starting from the potential overlaps of length 1, we iteratively compute the potential overlaps of increasing length by prepending each character c ∈ Σ to each k-long potential overlap β * (stored in the list Last), and we determine if cβ * is an (k + 1)long potential overlap: in this case we store the potential overlap in the list New.
The lists Last and New store the potential overlaps computed at the previous and at the current iteration respectively. By our previous observation, this procedure computes all potential overlaps. Observe that a potential overlap β * is an overlap iff pref(β * ) > 0. Since each potential overlap is a suffix of some input string, there are at most nm distinct suffixes, where m and n are the length and the number of input strings, respectively. Each query suff(·), pref(·), substr(·) requires O(1) time, thus the time complexity of all of such queries is O(nm). Given two distinct strings β 1 and β 2 , when |β 1 | = |β 2 | then no input string can be in both listpref(β 1 ) and listpref(β 2 ). Since each overlap is at most m long, the overall time spent in the listpref(·) queries is O(nm). The first phase produces (line 7) the set of disjoint basic arc-sets ARC(ǫ, ǫβ, R p (β)) for each overlap β, whose union is exactly the set of arcs of the overlap graph. Recall that listpref(β) gives the set of reads with prefix β, which has been denoted by R p (β).
The following definitions and lemma are fundamental in the design of Algorithm 2, since they allow to restrict the search that determines whether a certain arc is transitive. We recall that each arc-set ARC(α, αβ, X) is actually encoded as the triple (α, αβ, X).
Definition 3. Let A = ARC(αγ, αγβ, X 1 ), B = ARC(γ, γβδ, X 2 ) be two arc-sets, such that B is terminal and X 2 ⊆ X 1 . Then B reduces A and the tuple (αγ, αγβ, Based on the previous definition, we say that the arcs of (αγ, αγβ, X 1 ∩ X 2 ) are removed by B, or B removes those arcs.
In Figure 1 and Table 1 an example of arc-sets is presented for a set of five reads. Lemma 2. Let (r 1 , r 2 ) be an arc with overlap β. Then (r 1 , r 2 ) is transitive if and only if (i) there exist α, γ, δ, η ∈ Σ + such that r 1 = γαβ, r 2 = βδη, and (ii) there exists an input read r 3 = αβδ such that (r 3 , r 2 ) is an irreducible arc of a nonempty correct arc-set ARC(α, αβδ, X) for all X such that r 2 ∈ X.
Assume now that such an input string r 3 exists, we will prove that the arc (r 3 , r 2 ) reduces (r 1 , r 2 ). First we prove that (r 3 , r 2 ) is irreducible. Assume to the Table 1: Example of arc-sets corresponding to the overlap graph of Figure 1. Among the five reported arc-sets, B, C and E are terminal, and only D is not correct (and non-terminal) since it represents reducible arcs of the overlap graph.
Observe that E reduces D (by Definition 3) and the residual arc-set D\E represents an empty set of arcs (all the four reducible arcs are removed).
A direct consequence of Lemma 2 is that a nonempty correct terminal arcset ARC(α, αβδ, X) implies that all arcs of the form (γαβ, βδη), with γ, η ǫ are transitive.
Algorithm 2 classifies the arcs of the overlap graph into reducible or irreducible by iteratively computing arc-sets of increasing extension length ℓ, starting from the basic arc-sets ARC(ǫ, ǫβ, R p (β)) of extension length ℓ = 0 obtained in the previous phase. By Lemma 2, we compute all correct arc-sets of a given extension length and we discard (according to Definition 3) all arcs of non-terminal arc-sets that are removed by terminal arc-sets. The set D is used to store the reads of the destination sets X of the computed terminal arc-sets. Notice that if ARC(α, αβ, X) is terminal, then all of its arcs have the same origin r = αβ, that is The processed arc-sets (that have the same extension length) are partioned into clusters C(·). We denote with C(α) the set of the arc-sets ARC(α, ·, ·) with a given extension α which are contained in the stack Clusters at a certain point of the execution of Algorithm 2. Since the arc-sets pushed to Clusters (lines 16,19) have an extension cα, then for each α there can be at most a cluster C(α) during the entire execution of Algorithm 2. Observe that at the first iteration (when the processed arc-sets are basic and have the same extension ǫ) there is only one cluster C(ǫ) that is the output of Algorithm 1.
Each cluster C(α) is processed independently from the other ones, and the set D is used to store the reads of the destination sets X of its terminal arc-sets ARC(α, αβ, X). A consequence of Lemma 2 is that, for each one of the computed terminal arc-sets ARC(α, αβ, X), all arcs in C(α) with a destination in X ∈ D and with an origin different from r = αβ are transitive and can be removed simply by removing X from all destination sets in the non-terminal arc-sets of C(α). Another application of Lemma 2 is that, when we find a terminal arc-set, then all of its arcs are irreducible, that is it is also correct. In fact by Lemma 2, an arc (α * αβ * , βδ), where β * (overlap of the arc) is a prefix of β, is classified as transitive in relation to the existence of a read r = αβ that is the origin of an arc (αβ, βδ) with (left) extension α. Since the algorithm considers arc-sets by increasing extension length, all arcs that have extensions shorter than |α| have been reduced in a previous step of the algorithm and thus terminal arc-sets computed by previous iterations contain only irreducible arcs. More precisely, the test at line 8 is true iff the current arc-set is terminal. In that case, at line 10 all arcs of the arc-set are output as arcs of the string graph, and at line 11 the reads in the destination set X is added to the set D that contains the destinations of C(α) that must be removed from the destination sets of non-terminal arc-sets.
For each cluster C(α), we read twice all arc-sets that are included in C(α). The first time to determine which arc-sets are terminal and, in that case, to determine the reads (see the set D) that must be removed from all destinations of the non-terminal arc-sets included in C(α). The second time to compute, from the non-terminal arc-sets ARC(α, β, X), the clusters C(cα), for c ∈ Σ, that will contain the nonempty arc-sets ARC(cα, β, X \ D) with extension cα consisting of the arcs that we still have to check if they are transitive or not. Notice that, in Algorithm 2, the cluster C(α) that is currently analyzed is stored in Current-Cluster, that is a list of the arc-sets included in the cluster. Terminal arc-sets are removed from CurrentCluster before computing the extended clusters C(cα) (see line 12). Moreover, the clusters that still have to be analyzed are stored in the stack Clusters. We use a stack to guarantee that the clusters are analyzed in the correct order, that is the cluster C(α) is analyzed after all the clusters C(α[i :]) where α[i :] is a generic suffix of α -Lemma 3 will show that a generic irreducible arc (r 1 , r 2 ) with extension α and overlap β belongs exactly to the clusters C(ǫ), . . . , C(α[3 :]), C(α[2 :]), C(α). Moreover, r 2 does not belong to the set D when considering C(ǫ), . . . , C(α[3 :]), C(α[2 :]), hence the arc (r 1 , r 2 ) is correctly output by the algorithm when considering the cluster C(α).
Each cluster C(α) is analyzed separately, and each arc-set in any given cluster is tested to determine if the arc-set is terminal. In that case (see the test condition at line 8), the arcs having origin in r = αβ and destination in a read of X are produced in output (see line 10). All such arcs have label α. Moreover the reads in the destination set X are added to the set D, initially empty for each cluster, that will contain the destination sets of all the terminal arc-sets computed at the current iteration. We will use this information in the next step in order to remove all the arcs that have α as suffix of the label.
After having analyzed all arc-sets in C(α), those arc-sets are scanned again. During this second scan, for each arc-set with overlap β and destination set X, we test if there is at least a read in X \ D, that is at least a read of the destination set has not been reduced. In that case, we split C(α) into the non-empty cα-cluster (recall that such a cluster is empty iff suff(cαβ) = 0).
Corollary 4. The set of arcs computed by the algorithm is a superset of the irreducible arcs of the string graph.
Lemma 5. Let ARC(α, αβ, X) be an arc-set inserted into a cluster by Algorithm 2. Then such an arc-set is correct.
Proof. Let e 1 be an irreducible arc (r 1 , r 2 ) of ARC(α, αβ, X), and let α 1 be respectively the extension and the overlap β of e 1 . Since e 1 ∈ ARC(α, αβ, X), then α is a suffix of α 1 , therefore we can apply Lemma 3 which implies that e 1 ∈ C(α). Since the only arc-set contained in C(α) to which e 1 can belong is ARC(α, αβ, X), then r 2 ∈ X which completes the proof.
We can now prove that no transitive arc is ever output.
Lemma 6. Let e 1 be a transitive arc (r 1 , r 2 ) with overlap β. Then the algorithm does not output e 1 .
Proof. Since e 1 is transitive, by Lemma 2 the two reads r 1 , r 2 are r 1 = γαβ, r 2 = βδη, and there exists an input string r 3 = αβδ such that the arc e 2 = (r 3 , r 2 ) with overlap βδ is irreducible. Moreover all correct arc-sets of the form ARC(α, αβδ, X) with r 2 ∈ X are nonempty and terminal.
Assume to the contrary that e 1 is output by Algorithm 2, and notice that such arc can be output only when the current cluster is C(α) and the current arc-set is ARC(γα, γαβ, X) with r 2 ∈ X.
By the construction of our algorithm, since the cluster C(γα) is nonempty, also C(α) is nonempty: let us consider the iteration when the current cluster is C(α). By Lemma 5 the arc-set ARC(α, αβδ, X 1 ) is correct, hence it contains the arc e 2 . But such arc-set satisfies the condition at line 8, hence r 2 ∈ D at that iteration. Consequently, C(α) cannot contain an arc-set with destination set with r 2 .
Theorem 7 is a direct consequence of Corollary 4 and Lemma 6.
Theorem 7. Given as input a set of strings R, Algorithm 2 computes exactly the arcs of the string graph. A similar argument shows that the number of arc-sets managed by Algorithm 2 is at most O(nm 2 ), since each suffix αβ can be considered for different extensions α. Moreover, differently from Algorithm 1, the time spent by Algorithm 2 for each string αβ is not constant. More precisely, for each cluster, besides computing substr(·), pref(·), suff(·), Algorithm 2 computes a union D (line 11) and the difference X \ D (line 16) -a direct inspection of the pseudocode shows that all other operations require constant time for each string αβ.
The union at line 11 and the difference X \ D at line 16 are computed for each string suffix αβ where β is a potential overlap. Let d(n) be the time complexity of those two operations on n-element sets (the actual time complexity depends on the data structure used). Therefore, the time complexity of the entire algorithm is O(nmd(n)). We point out that a representation based on an n-long bitvector (like the one we will discuss in the following) implies an O(n) time complexity for each union and difference.
As a consequence, the time complexity of our algorithm is O(nm(n + m)). On the other hand, we conjecture that this time complexity is highly pessimistic, since usually the number of potential overlaps is smaller than the worst case.
One of the main features of our approach is that we operates only on the (po-Algorithm 2: Compute the string graph Input : The set Basics of the basic arc-sets ARC(ǫ, ǫβ, R p (β)) Output: The arcs of the string graph of R 1 Clusters ← empty stack; 2 Push Basics to Clusters; 3 while Clusters is not empty do Push E xtendedClusters[c] to Clusters; tentially compressed) FM-index of the collection of input reads. To achieve that goal we cannot use the naïve representation of a string ω, but we must employ a BWT-based representation. More precisely, we represent a string ω with the ω-interval q(ω) = [b ω , e ω ] and the fact that ω is suffix in some read of R with the ω$-interval q(ω$) = [b ω$ , e ω$ ] on the BWT, hence using four integers for each string.
We need to show how to compute efficiently this representation of ω as well as pref(ω), suff(ω), and substr(ω). The following proposition is instrumental and can be verified by a direct inspection of Algorithms 1 and 2: both algorithms compute strings by prepending characters, that is they need to obtain the (representation of the) string cω from the (representation of the) string ω that has been processed previously.
Proposition 8. Let ω be a string processed by Algorithms 1 or 2, and such that |ω| > 1. Then the string ω[2 :] is processed by the same algorithm before ω.
Since each string ω considered by the algorithm is a substring of some input read, we can represent ω in constant space by the boundaries (i.e., the first and the last index) of q(ω) and q(ω$), instead of using the naïve representation with O(|ω|) space. Consequently, the algorithm operates only on the (potentially compressed) FM-index of the collection of input reads.
This representation allows to answer each query pref(ω), suff(ω) and substr(ω) in O(1) time. In fact, given q(ω) = [b ω , e ω ], then substr(ω) = e ω − b ω + 1 and Moreover, listpref(ω) corresponds to the set of reads appearing in the interval q($ω) of the GSA. Notice that no read can appear twice in such interval, hence a linear scan of the interval q($ω) suffices.
A further optimization of the representation is possible. Recall that q(ω) = [b ω , e ω ] and q(ω$) = [b ω$ , e ω$ ]. Notice that b ω = b ω$ , since the sentinel $ is lexicographically smaller than all characters in Σ. Hence the two intervals q(ω) and q(ω$) can be represented by the three integers b ω , e ω , e ω$ .
Instead of storing directly the three integers b ω , e ω , e ω$ , we use two n(m + 1)long bitvectors, requiring 2n(m+1)+o(nm) bits and allowing to answer in constant time to rank and select queries (Clark and Munro, 1996;Jacobson, 1989).
Algorithm 1 mainly has to represent the set of potential overlaps (i.e., the lists Last and New). At each iteration, the potential overlaps in Last (in New, resp.) have the same length, hence their corresponding intervals on the BWT are disjoint. For each potential overlap β ∈ Last (in New, resp.) represented by the triple (b β , e β , e β$ ), the first bitvector has 1 in position b β and the second bitvector has 1 in positions e β$ and e β -since each q(β) is disjoint from all other intervals, the i-th interval is represented by the i-th 1 in the first bitvector, and the 2i-th and (2i + 1)-th 1s of the second bitvector.
Algorithm 2 mainly has to represent clusters and, for each cluster, the set D containing the reads that must be deleted from all destination set of the nonterminal arc-sets of the cluster. A cluster groups together arc-sets whose overlaps are either pairwise different or one is the prefix of the other. Thus, the corresponding intervals on the BWT are either disjoint or one contained in the other (i.e., partial overlap of the intervals cannot happen). Moreover, also the destination set of the basic arc-sets can be represented by a set of pairwise disjoint or contained intervals on the BWT (since listpref(β) of line 7 corresponds to the reads of the interval q($β) on the GSA). The following proposition, which can be proved by a similar argument as Lemma 3, describes the relation between destination sets.
Given a cluster C(α), let γ be the longest common prefix of all strings αβ such that ARC(α, αβ, X) is an arc-set of the cluster. Then the numbers that represent the strings αβ are all contained in the γ-interval: therefore it suffices to consider numbers in the interval [b γ , e γ ] instead of the interval [1, n(m + 1)]. A more compact representation of the data that we have to manage consists of two (e γ − b γ + 1)long vectors V b , V e of integers, and a bitvector B x of length e $γ − b $γ + 1 -where [b $γ , e $γ ] is the $γ-interval. Each entry V b [i] (V e [i], resp.) is the number of arcsets ARC(α, αβ, ·) in C(α) with initial (final, resp.) boundary b γ + i. Moreover, the destination sets that are considered in the cluster are contained in the set of reads in q($γ), hence the bitvector B x is to encode the set D: in fact B x [i] is 1 iff the (b $γ + i)-th read (in lexicographic order), belongs to D.
Notice that the arc-sets ARC(α, αβ, X) in each cluster C(α) are sorted according to the lexicographic order of β: this fact allows to coordinate the sequential scan of V b and V e that corresponds to the for loop at lines 13-16 with a scan of the list of destination sets (one for each arc-set) that must be updated at line 16. Also notice that V b , V e , and B x are constructed at lines 11-16 and exploited at lines 13-16, hence the data structures necessary to answer rank and select queries in constant time must be built once for each cluster just before line 13.

Experimental Analysis
A C++ implementation of our approach, called FSG (short for Fast String Graph), has been integrated in the SGA suite and is available at http://fsg.algolab.eu under the GPLv3 license. Our implementation uses the Intel ® Threading Building Blocks library in order to manage the parallelism. The software is conceptually divided in the two phases illustrated in the previous section, and each phase has been implemented as a computational pipeline using the pipeline construct made available by the library.
We have evaluated the performance of FSG with three experiments: the first experiment compares FSG and SGA on a standard benchmark of 875 million 101bp-long reads sequenced from the NA12878 individual of the International HapMap and 1000 genomes project (extracted from ftp://ftp-trace.ncbi.nih.gov/1000genomes/ The second experiment, on an Escherichia coli dataset, aims at investigating the cause of the speedup obtained by FSG with respect to SGA. Finally, the third experiment (which is on a synthetic dataset obtained from the Human chromosome 1), studies the effect of coverage on the performance of FSG and SGA. All experiments have been performed on an Ubuntu 14.04 server with four 8-core Intel ® Xeon E5-4610v2 2.30GHz CPUs (hyperthreading was enabled for a total of 16 threads per processor). The server has a NUMA architecture with 64GiB of RAM for each node (256GiB in total). To minimize the effects of the architecture on the executions, we used numactl to preferably allocate memory on the first node where also the threads have been executed (with 32 threads also the second node was used).
We have run SGA with its default parameters, that is SGA has computed exact overlaps after having corrected the input reads. We could not compare FSG with Fermi, since Fermi does not split its steps in a way that allows to isolate the running time of the string graph construction-most notably, it includes reads correction and scaffolding. Since the string graphs computed by FSG and SGA are essentially the same, we have not focused our analysis on the quality of the resulting assemblies, but we give a brief analysis hinting that it is not possible to determine which assemblies are better (see Table 2).
For genome assembly purposes, only overlaps whose length is larger than a user-defined threshold are considered. The value of the minimum overlap length threshold that empirically showed the best results in terms of genome assembly quality is around the 75% of the read length (Simpson and Durbin, 2012). In order to assess how graph size affects performance, different values of minimum overlap length (called τ) between reads have been used (clearly, the lower this value, the larger the graph). The minimum overlap lengths used in this experimental assessment are 55, 65, 75, and 85, hence the chosen values test the approaches also on larger-than-normal (τ = 55) and smaller-than-normal (τ = 85) string graphs.
Another aspect that we wanted to measure is the scalability of FSG for a differ- In all cases, we have measured the elapsed (wall-clock) time and the total CPU time (the time a CPU has been working).
In terms of memory, SGA does not maintain the computed string graph in memory, hence its peak memory usage is only dependent on the input size and in these experiments was always about 63GiB. Also the peak memory usage of our approach was approximately equal to 138GiB for all the configurations. As a consequence, the memory usage of our approach is practically only dependent on the input size since it compactly stores the arc-sets of the first phase and the stack maintaining the clusters to be processed in the second phase does not grow to have more than |Σ| · (m − τ) elements. Table 3 summarizes the running times of both approaches on the different configurations of the parameters. Notice that FSG approach is from 2.3 to 4.8 times faster than SGA in terms of wall-clock time and from 1.9 to 3 times in terms of CPU time. On the other hand, FSG uses approximately 2.2 times the memory used by SGA -on the executions with at most 8 threads.
While FSG is noticeably faster than SGA on all instances, there are some other interesting observations. The combined analysis of the CPU time and the wallclock time on at most 8 threads (which is the number of physical cores of each CPU on our server) suggests that FSG is more CPU efficient than SGA and is able to better distribute the workload across the threads. The latter value of 8 threads seems to be a sweet spot for the parallel version of FSG.
On a larger number of threads, and in particular the fact that the elapsed time of FSG on 32 threads is larger than that on 16 threads suggests that, in its current form, FSG might not be suitable for a large number of threads. However, since the current implementation of FSG is almost a proof of concept, future improvements to its codebase and a better analysis of the race conditions of our tool will likely lead to better performances with a large number of threads. Furthermore, notice that also the SGA algorithm, which is (almost) embarrassingly parallel and has a stable implementation, does not achieve a speed-up better than 6.4 with 32 threads. As such, a factor that likely contributes to a poor scaling behaviour of both FSG and SGA could be also the NUMA architecture of the server used for the experimental analysis, which makes different-unit memory accesses more expensive (in our case, the processors in each unit can manage at most 16 logical threads, and only 8 on physical cores).
Notice that, FSG uses more memory than SGA. The reason is that genome assemblers have to correctly manage reads extracted from both strands of the genome. In our case, this fact has been addressed by adding each reverse-andcomplement read to the set of strings on which the FM-index has been built, hence immediately doubling the size of the FM-index. Moreover, FSG needs some additional data structures to correctly maintain potential overlaps and arc-sets, as described in Section 4. The main goal of FSG is to improve the running time, and not necessarily to decrease memory usage.
We also wanted to estimate the effect of the optimizations discussed in Section 4, by measuring the number of backward extensions performed. We have instrumented FSG and SGA to count the number of accesses to the FM-index (which can happen only when computing a backward extension) and we have run both program on the NA12878 dataset with τ = 85: FSG has made 877 · 10 9 accesses, while SGA has made 947 · 10 9 accesses, i.e., SGA has made 8% more backward extensions than FSG. This result confirms that clustering together arcsets and avoiding unnecessary backward extensions actually improves on SGA's strategy. On the other hand, such difference on the number of backward extensions is unable to fully explain the different running times of FSG and SGA - Table 3 shows that, on this specific instance, SGA needs 2.5x the time used by SGA. Therefore we have designed a second experiment to investigate the main causes of the different performances of FSG and SGA.
In our second experiment, we have run both FSG and SGA on the Escherichia coli dataset downloaded from http://www.ebi.ac.uk/ena/data/view/CP009789 under valgrind (Nethercote and Seward, 2007) to measure the instruction and memory access patterns. Since running a program under valgrind increases the running time by two orders of magnitude, it was not feasible to use the same dataset as the first experiment, but we had to use a much smaller one. The results of this experiment are summarized in Table 4.
Our initial conjecture was that the better efficiency achieved by FSG originated from operating only on the FM-index of the input reads and by the order on which extension operations (i.e., considering a new string cα after α has been processed) are performed. These two characteristics of our algorithm allow to eliminate the redundant queries to the index which, instead, are performed by SGA. In fact, FSG considers each string that is longer than the threshold at most once, while SGA potentially reconsiders the same string once for each read in which the string occurs. A consequence of our conjecture should have been fewer memory accesses and cache misses. Unfortunately the results we have obtained on memory accesses are inconclusive: the number of memory accesses made by SGA are almost twice as many as those made by FSG. On the other hand, the number of cache misses that result in RAM accesses hints that SGA is much more efficient in that regard. We have estimated the total time spent in memory accesses using the values suggested by Valgrind: 1 cycle for each cache access, 10 cycles for each level 1 cache miss, and 100 cycles for each cache miss that results in a RAM access. Overall, FSG is more efficient than SGA, but definitely less than 10% more efficient: hence this reason alone cannot justify the difference in running times of the two programs.
The fact that FSG consists of several linear scans suggests that it should be able to better exploit the superscalar features of modern CPUs. In fact, our analysis of branch mispredictions confirms this fact: the ratio of the total number of branch predictions made by FSG is about 75% of those made by SGA, and the total number of branch mispredictions made by FSG is less than 50% of those made by SGA.  We have compared the string graphs produced by FSG and SGA, since their respective notions of string graphs are slightly different. More precisely, FSG keeps multiple arcs with distinct labels between the same pairs of vertices, while SGA retains only the arc having the shortest label. When the minimum overlap length is 65bp, the string graph computed by FSG had ∼ 3.5% more arcs than the one computed by SGA, but the impact on the actual assembly is not relevant (see Table 2).  A third experiment has been performed with the goal of studying the scalability of FSG on larger coverage values. We have extracted random reads from the human chromosome 1, with different average coverage values and different read lengths. Since the goal of this experiment is analyzing the running time and memory usage, while we are not interested into the accuracy of the predictions, the reads do not contain errors (i.e., they are substrings extracted from the reference genome). The results of these experiments are reported in Tables 5 and 6. The comparison between FSG and SGA confirms the results of the first experiment, and shows that for large coverage values, FSG becomes even faster than SGA but uses an increasing amount of memory.

Conclusions and future work
We present FSG: a tool implementing a new algorithm for constructing a string graph that works directly querying a FM-index representing a collection of reads, instead of processing the input reads. Our main goal is to provide a simpler and fast algorithm to construct string graphs, so that its implementation can be easily integrated into an assembly pipeline that analyzes the paths of the string graph to produce the final assembly. Indeed, FSG could be used for related purposes, such as transcriptome assembly (Lacroix et al., 2008;Beretta et al., 2014), and haplotype assembly (Bonizzoni et al., 2003), and variant detection via aligning paths of the string graph against a reference genome. These topics are some of the research directions that we plan to investigate. More precisely, our algorithm uses string queries that are efficiently implemented using the information provided by the index and takes advantage of a lex-icographic based-ordering of string queries that allows to reduce the total number of such queries to build the string graph. Since FSG reduces the total number of queries and does not process the input to compute the transitive reduction as done by SGA, the current state-of-art tool for computing the string graph, we are able to show that FSG is significantly faster than SGA over genomic data. It would be interesting to test our implementation to compute the string graph for a large collection of strings with different characteristics than genomic reads, such as for example when the alphabet is of larger size or the input data consists of strings of variable length.