HyDA-Vista: towards optimal guided selection of k-mer size for sequence assembly

Abstract

Motivation

Intimately tied to assembly quality is the complexity of the de Bruijn graph built by the assembler. Thus, there have been many paradigms developed to decrease the complexity of the de Bruijn graph. One obvious combinatorial paradigm for this is to allow the value of k to vary; having a larger value of k where the graph is more complex and a smaller value of k where the graph would likely contain fewer spurious edges and vertices. One open problem that affects the practicality of this method is how to predict the value of k prior to building the de Bruijn graph. We show that optimal values of k can be predicted prior to assembly by using the information contained in a phylogenetically-close genome and therefore, help make the use of multiple values of k practical for genome assembly.

Results

We present HyDA-Vista, which is a genome assembler that uses homology information to choose a value of k for each read prior to the de Bruijn graph construction. The chosen k is optimal if there are no sequencing errors and the coverage is sufficient. Fundamental to our method is the construction of the maximal sequence landscape, which is a data structure that stores for each position in the input string, the largest repeated substring containing that position. In particular, we show the maximal sequence landscape can be constructed in O(n + n log n)-time and O(n)-space. HyDA-Vista first constructs the maximal sequence landscape for a homologous genome. The reads are then aligned to this reference genome, and values of k are assigned to each read using the maximal sequence landscape and the alignments. Eventually, all the reads are assembled by an iterative de Bruijn graph construction method. Our results and comparison to other assemblers demonstrate that HyDA-Vista achieves the best assembly of E. coli before repeat resolution or scaffolding.

Availability

HyDA-Vista is freely available [1]. The code for constructing the maximal sequence landscape and choosing the optimal value of k for each read is also separately available on the website and could be incorporated into any genome assembler.

Introduction

The ability to accurately assemble genomes is a fundamental problem in bioinformatics that is vital to the success of many scientific projects, including the 10,000 vertebrate genomes (Genome 10K) [2], Arabidopsis variations (1001 genomes) [3], human variations (1000 genomes) [4], and Human Microbiome Project [5]. The genome assembly process aims to build contiguous sequences, called contigs, predominantly from short read sequencing data. Other sources of information have also been used to boost the accuracy, including genetic linkage data [6], optical mapping data [7], and longer sequencing reads (e.g. PacBio data) [8]. A potential source of information that has not been fully explored is the information contained in phylogenetically-close genomes. The genome of an individual of the same species or that of a phylogenetically-close species can potentially be used as an extra source of information, and increase the assembly quality. We argue that genome assemblers can benefit from using a reference genome to help guide the assembly process, particularly in those regions of the genome that are pervaded by repetitive sequences.

In Eulerian sequence assembly [9, 10], a de Bruijn graph is constructed with a vertex v for every (k − 1)-mer present in an input set of reads, and an edge $\left(v-v^{\prime }\right)$ for every observed k-mer in the reads with (k − 1)-mer prefix v and (k − 1)-mer suffix $v^{\prime }$. A contig corresponds to a non-branching path through this graph. SPAdes [11], IDBA [12], Euler-SR [13], Velvet [14], SOAPdenovo [15], ABySS [16] and ALLPATHS [17] all use this paradigm for assembly. The majority of de Bruijn graph based assemblers follow the same general outline: break the (possibly error corrected) reads into k-mers, construct the de Bruijn graph on the set of resulting k-mers, simplify the de Bruijn graph, resolve the repeated regions by using mate-pair information, and construct the contigs (simple paths in the de Bruijn graph). Therefore, the majority of assemblers require or allow the value of k to be specified by the user.

The problem of determining an appropriate value of k for the de Bruijn graph construction is important since it has a direct impact on assembly quality; stated very briefly, when k is too small the resulting graph is complicated by spurious edges and vertices, and when k is too large the graph becomes too sparse and possibly disconnected. Repetitive regions are especially problematic for genome assembly since they inadvertently result in spurious edges and vertices in the de Bruijn graph [18] and are very sensitive to the choice of k. There has been a significant effort in developing algorithms that will choose an ideal value for k by preprocessing the sequence reads, and thus, reduce the complexity of the de Bruijn graph [11, 12, 19].

A more obvious combinatorial approach for building a simplified de Bruijn graph would be to allow the value of k to vary; having a larger value of k where the graph is more complex and a smaller value of k where the graph would likely contain fewer spurious edges and vertices. A major difficulty in implementing this approach is determining a practical method that makes this idea feasible assembling large genomes. Peng et al. [12] and Bankevich et al. [11] both introduced assemblers that use various values of k. IDBA builds the de Bruijn graph in an iterative manner from k = k _min to k = k _max ; these values of k are predetermined and (by default) do not change for different datasets or genomes. At iteration i, the de Bruin graph for k _i is built from the current set of reads and the contigs for that graph are constructed, then all the reads that align to at least one of those contigs are removed from the current set of reads. In the next iteration the graph is built by converting every edge from the previous graph to a vertex while treating contigs as edges. SPAdes [11] uses a similar approach but uses all the reads at each iteration.

While this method has been shown to greatly improve assembly quality [11, 12], it is not efficient since all the reads are assembled at each iteration. Thus, one challenge that remains to be addressed is how to efficiently determine which values of k should be used for this iterative assembly process and how to assign a k-mer value for each read. If this could be accomplished prior to assembly of the de Bruijn graph(s) then these iterative assembly methods could be made more efficient without degrading the assemblies quality.

Our contribution. We introduce an efficient algorithm for determining an optimal value of k for each read prior to constructing the de Bruijn graph, and implement this method into a modified version of HyDA, a de novo assembler developed by Movahedi et al. [20]. This modified assembler, which we refer to as HyDA-Vista, takes as input a phylogenetically-close genome and a set of paired-end reads. Imperative to HyDA-Vista is the construction of the maximal sequence landscape, which is a data structure that stores for every position in the input string, the longest repeat containing it. Prior to de Bruijn graph construction, HyDA-Vista constructs the maximal sequence landscape for the phylogenetically-close genome, and aligns the reads to the reference genome. The alignment and landscape allows the optimal value of k for each read be determined in linear time in the length of the read, provided the read is longer than the longest repeat. These values of k are "optimal" in the sense that for unchanged parts of the genome, the de Bruijn graph will have no spurious edges or vertices if there are no sequencing errors, and the length of the repeat is smaller than the read length. Unaligned reads are assigned a default value of k. After the assignment of values of k to each read, HyDA-Vista constructs the de Bruijn graph in an iterative manner.

Our approach for choosing values of k for each read takes into consideration the repeat structure of the genome, which enables us to avoid overly-complex regions of the graph since the assignment of values of k to reads is done prior to the assembly rather than during the assembly. We compare HyDA-Vista versus IDBA [12], SPAdes [11], SOAPdenovo [15], ABySS [16] and HyDA [20]. Our results demonstrate that HyDA-Vista produces the best assembly of E. coli before repeat resolution or scaffolding. We aim to achieve the best assembly without repeat resolution and scaffolding and note that such methods could be applied to all these initial assemblies. Lastly, we demonstrate that this method improves the efficiency of iterative assembly.

Roadmap. We review related tools for the problem in the remainder of this section. Section then sets notation and formally lays down the problem and the data structures that will be used for the construction of the maximal sequence. We formally define the maximal sequence landscape in Section. Section gives details of HyDA-Vista. In Section we give results that demonstrate how HyDA-Vista compares against competing assemblers. Finally, Section offers reflections and some future areas of research that warrant investigation.

Related work. The re-sequencing methods include LOCAS and SUPERLOCAS [21, 22], e-RGA [23], Colombus module of Velvet [24] and IDBA-Hybrid [25]. Since these methods aim to determine structural variations between species, and require extremely high sequence similarity to produce reasonable results, they have only been applied to individuals of the same species. Our focus is to produce high quality de novo assemblies using homology information contained in the reference genome of the same species or phylogenetically-close species. Gnerre et al. [26] also consider how to improve assembly quality by using the alignment of reads to a reference genome. Their method simultaneously builds a de novo assembly from the reads and aligns these same reads to one or more related genomes. The alignment is then used to improve the assembly quality, e.g., reads that were not used in the assembly are incorporated into the assembly using the alignment.

Complementary to the work of Chikhi and Medvedev [19], Peng et al. [12], and Bankevich et al. [11], there has been an effort in developing methods that use paired-end data to constrain the construction of the de Bruijn graph [21, 27–30]. Medvedev et al. [29] introduced the concept of a paired de Brujin graph. Since the insert size is variable among mate pairs, this method requires that all the paths within some threshold be considered in order to ensure an edge is not missed. Thus, Bankevich et al. [11] improve upon this idea by developing the rectangle graph, which eliminates the need to consider all paths. Vyahhi et al. [30] furthered this study of rectangle graphs for genome assembly. These methods merit mentioning the goal of these methods is the same as the goal of increasing the value of k in certain regions; both aim to minimize spurious edges and branching in the graph but in a different manner.

Determining all maximal exact repeats in a string has been previously studied [31–33]. It has been shown that all maximal repeats of a string can be found and stored in O(n)-time and O(n)-space using a suffix tree (although the output maybe of size Θ(n²)) or directed acrylic graph [32]. Therefore, the maximal sequence landscape, which we define in this paper, can be constructed from either a suffix tree or a directed acrylic graph in O(n)-time and O(n)-space using these algorithms directly or adapting them. However, the constant in the order notation of the space complexity of these constructions is relatively large. The algorithm we present uses a suffix array and thus, requires linear space with a smaller constant and O(n log n)-time. Hence, we pay a log n cost in time to remove the large constant from the linear space time. We also note that the related problem of finding inexact maximal repeats also has been previously studied [34–37].

Background

Strings. Consider a fixed alphabet Σ = {σ₁,..., σ _m } and a total order ≤ _L defined over Σ = Σ ∪ {$} where $ ∉ Σ, and for all σ ∈ Σ we have $ ≤ _L σ. We denote a finite string s as s₁s₂ ... s _n , where s _i ∈ Σ. We use s _ij , where 1 ≤ i ≤ j ≤ n, to indicate substring s _i s _i+1 
 s _j of string s. We call substrings $s_i^{\mathsf{\text{Pre}}}=s_{1i}$ and $s_i^{\text{Suf}}=s_{in}$ with i ∈ {1, . . . , n} the i^th prefix and i^th suffix of s respectively. Based on the total order ≤ _L , we define a lexicographical total order on the strings in Σ^∗.

Arrays. We denote arrays of integers by all capital letter strings like A, SP, LCP, etc. A[i], with 1 ≤ i ≤ |A|, stands for the integer in the i^th cell of array A. Also, A[i, j] indicates the projection of A onto indices i to j, inclusive of both ends. For an array A, with |A| = n, that holds a permutation of integers {1, ..., n}, index array of A is another array I(A) with |I(A)| = n such that I(A)[i] = j if and only if A[j] = i.

Suffix and longest common prefix arrays. SA _s , for some string s, denotes the suffix array associated with s [38]. SA _s [i] = l for i ∈ {1, ... , n} if and only if $s_l^{\text{Suf}}$ is the i^th string in the lexicographically sorted list of all suffixes of s. We also indicate the longest common prefix array of some string s with LCP _s , and LCP _s [i] = l for i ∈ {1, . . . , n − 1} if and only if the length of the longest common prefix between $s_{\mathsf{\text{S}}{\mathsf{\text{A}}}_{s\left[i\right]}}^{\mathsf{\text{Suf}}}$ and $s_{\mathsf{\text{S}}{\mathsf{\text{A}}}_{s\left[i\right]}+1}^{\mathsf{\text{Suf}}}$ equals to l.

Approach

The maximal sequence landscape

We formally define the maximal sequence landscape in this section. Clift et al. [32] introduced the concept of a sequence landscape, which is a data structure that stores the occurrences of any substring from a source string s in a target string t. In set representation, the sequence landscape L _t|s of a target string t with respect to a source string s is defined as a set of tuples {m₁, m₂, ..., m _l }, where m _i = (b _i , e _i , f _i ) corresponds to the occurrence of substring sb _i e _i = sb _i sb _{i +1} 
 se _i from s in t with frequency f _i . If s and t are equal then the sequence landscape categorizes all repeated substrings in the source string s. We define to this special case where s = t as the self sequence landscape. Figure 1 illustrates an example of a self sequence landscape and a sequence landscape. Given a position i of the input string s, all the repeated substrings containing s _i can be recovered from the self sequence landscape in linear-time in the number of different repetitions.

Figure 1

(left)The self sequence landscape for CATCATTTG, and (right) the sequence landscape of another string GGCATCATTGGGTATAACC with respect to CATCATTTG. The mountains (light grey or black) demonstrate occurrences of substrings of source string in the target string. Numbers at the peak of the mountains denote the frequency of occurrence. The maximal sequence landscape is highlighted in light grey, and the arrows demonstrate the ascent and descent of the landscapes.

The occurrences of the substrings in the source string are defined as mountains. This terminology reflects the visual representation that was first introduced by Clift et al. [32] that illustrates each occurrence as a mountain having height equal to the length of the substring, i.e. the height of mountain mi of L _t|s is denoted as h(m _i ) and equal to e _i − b _i + 1. The peak of each mountain is labelled with the frequency of the substring corresponding to it. In Figure 1 (left), the substring CAT is represented as two mountains each of which has a height equal to three and frequency equal to two.

We say that a mountain m _j = (b _j , e _j , f _j ) in a landscape L _t|s covers index i and denote it by m _j Δ i if and only if i ∈ {b _j ..., e _j }. Hence, the cover set of a specific index i of the sequence landscape L _t|s is the set of all the mountains that covers i. We denote the cover set as $C_{L_{t|s}}\left(i\right)$ and define it as follows:

$$C_{L_{t|s}}\left(i\right)=\left\{m_j|m_j\Delta i,h\left(m_j\right)>1.f_j>1\right\}.$$

(1)

Lastly, we define the summit of index i as the highest mountains in its cover set. We denote the summit of i by SLt s (i) and define it as follows:

$$S_{L_{t|s}}\left(i\right)=\left\{m_j|h\left(m_j\right)\ge h\left(m_k\right)\forall m_k\in C_{L_{t|s}}\left(i\right)\right\}.$$

(2)

Please note that the summit of index i can be empty or non-unique so the height of summit of index i is defined to be zero for empty set.

Definition 1 The maximal sequence landscape, which we denote as $L_{t|s}^{*}$, is the set of the summits of all positions in s that have frequency greater than one. $L_{t|s}^{*}$can be formally defined as follows: $L_{t|s}^{*}=\left\{S_{L_{t|s}}\left(i\right)|i=1,\dots ,n\right\}$The maximal sequence landscape is highlighted in light grey in Figure 1.

The maximal sequence landscape is obtained from the sequence landscape by removing all mountains except those that are highest and have frequency greater than one at each position. In the case of the maximal sequence landscape constructed from a self sequence landscape, this results in a data structure containing the longest repeat at each position of the input string. In Subsection we give an algorithm that builds the maximal sequence landscape and returns an array containing the length of the longest repeat at each position of the input string. Therefore, given a position i in s, we can determine the length of the longest repeat in s containing that position in constant time by simply indexing the maximal sequence landscape at position i. By choosing a value for k that is larger than the length of this repeat it can guaranteed that there will be no branching in the corresponding vertices of the de Bruijn graph, if the same substring is not repeated in changed parts of the genome that is being assembled. This is our idea based on which we determine the optimal value of k for each read. We consider the maximal sequence landscape constructed from the self sequence landscape for the remainder of this paper since it is what is used by HyDA-Vista.

Methods

Algorithm 1 gives an overview of HyDA-Vista algorithm. We explain each of these steps in detail in the subsequent subsections.

Algorithm 1 An overview of HyDA-Vista

1: Build the maximal sequence landscape for the reference genome.

2: Align all reads to the reference using BWA.

3: For each aligned read: assign a value of k using the maximal sequence landscape.

4: Unaligned reads are assigned a value of k using a heuristic.

5: The de Bruijn graph is constructed in an iterative manner, as shown in Algorithm 2.

Construction of the maximal sequence landscape

In this section we demonstrate that the maximal sequence landscape for an input string s can be built in O(n + n log n)-time and O(n)-space using a simple algorithm, where n is the length of s. Our method relies on the use of suffix array and longest common prefix array and thus, begins by building the suffix array (SA _s ) and the longest common prefix array (LCP _s ). This construction can be done in O(n)-space and O(n)time [39]. Two other auxiliary data structures are constructed at the beginning of the algorithm. However, we delay the definition of these to later in this section. The algorithm then iterates through each position of s and finds the longest repeated substring in s that contains it using SA _s , LCP _s , and the auxiliary data structures. An important aspect of our algorithm that allows us to achieve O(n + n log n)-time is that we only search the interval of SA _s that is between the indices SA_min and SA_max at each iteration; not the entire array. In other words this invariant holds at each iteration of our algorithm: [SA_min, SA_max] holds the the interval in the suffix array that corresponding suffixes share the same prefix. This prefix is the longest repeat that has been seen so far and covers that position. Thus, each time the largest repeated substring is found for a particular position, the maximal sequence landscape, SA_min, and SA_max are updated for search at the next iteration. Two possibilities exist at each iteration i of the algorithm when we are processing si;

• (a) The longest repeated substring at position i - 1 can be extended by appending s _i . The maximal sequence landscape, SA_min, and SA_max are updated.

• (b) The longest repeated substring at position i - 1 cannot be extended by appending s _i (either the extended string does not occur or it does occur but its frequency is one). Let p = s _j ⋯ s _i−1 be the longest repeated substring yet found that contains s _i−1 , and $p^{\prime }=s_{j+1}\cdots s_i$ be the string obtained by removing the first letter of p and appending s _i . If the frequency of occurrence of the $p^{\prime }$ is greater than one, then the maximal sequence landscape, and the search interval is updated as in (a). Otherwise, the search for the longest repeated substring continues by eliminating the first character of $p^{\prime }$ each time until a repeating match is found or the null string is reached. If the null string is reached then the maximal sequence landscape is empty at that position and the search interval is updated to [1, n].

The search interval contains all indices in SAs for which the corresponding suffixes have the current longest repeated substring as a prefix. In (a), the interval is updated by performing binary search. In (b) the search interval is no longer valid since we removed a letter from the beginning of the current longest repeated substring and we need (a more complicated) scheme to efficiently find the correct search interval. To accomplish this we need two auxiliary data structures that are constructed at the beginning of the algorithm: the SP _s array, and an ordered binary search tree containing all consecutive intervals of LCP _s . SP _s [j] holds the index in SAs that is obtained by removing the first letter from the beginning of $s_j{s}_{j+1}\cdots s_n$ in order to obtain $s_{j+1}\cdots s_n$. This array can be built in linear time by scanning the index array of SA _s . Thus, to find the correct interval in (b), we locate an index of SA _s (denoted as sp) where the corresponding suffix contains $s_{j+1}\dots s_{i-1}$ as a prefix, and find the largest interval around sp where all the suffixes in the interval have $s_{j+1}\dots s_{i-1}$ as prefix. This is the new search interval. The sp index can be found in constant time by correctly indexing SP _s . The second step is equivalent to finding the largest interval [d, u] around sp that for all j ∈ [d, u] we have $\mathsf{\text{LC}}{\mathsf{\text{P}}}_s\left[j\right]\ge \left|p^{\prime }\right|-1$. This can be done in O(log n)-time and O(n)-space using an ordered binary search tree.

The time complexity can be seen by first noting that at each iteration of the algorithm the maximal sequence landscape either ascends by one after a number of descents (possible zero) or it is undefined after a number of nonzero descents, and each of these ascents or descents require O(log n)-time. Note that in (a) the maximal sequence landscape is ascending, and in (b) the maximal sequence landscape is descending, and the frequency of a substring in s can be determined in constant time using LCP _s . Second, since each time it ascends one character from s is processed and the number of ascents equals the number of descents, the total number of ascents and descents is 2n. Therefore, since the data structures are constructed in O(n)-time, and since there are at most 2n ascents or descents which take O(log n)-time, the running time of the algorithm is O(n + n log n).

Assignment of k-mer sizes to reads

We assign values of k to the reads using the maximal sequence landscape constructed for the reference genome by first aligning the reads to the reference genome using BWA (version 0.7.4) [40] in pairedend mode. We consider all forward and reverse alignments of every read. Let p be the position in the reference genome where a read of length l aligns, and let k^∗ be the maximum of {MSL[p] + 1, MSL[p + 1] + 1, ... , MSL[p + l] + 1}, where MSL is the maximal sequence landscape array that contains the height of the maximal sequence landscape at each position. We compute k^∗ for each forward alignment and let K^∗ be the set of all these values. The optimal value of k for the (forward) read is equal to the maximum value in K^∗. We follow the same procedure for the reverse alignments with the exception that we compute the reverse complement of the read. Thus, the optimal values of k can be computed in linear time in the length of the read.

If the computed k-mer size (maximum of all maximal sequence landscape heights of all aligned nucleotides) is larger than the read length, then a default value (k = 77 is the default) is used instead. Unaligned reads are also assigned a default k-mer value (k = 55 is the default).

The de Bruijn graphs

Let R = {r₁,..., r _N } denote the set of reads. We also denote the k-mer size assigned to r _i in the previous section by K(r _i ). In the first step of constructing the assembly de Bruijn graphs, we partition R into R _k := {r ∈ R | K(r) = k}, in which k ranges from k_min = min _r∈R K(r) to k_max = max _r∈R K(r). The HyDA-Vista assembly procedure, shown in Algorithm 2, iteratively builds de Bruijn graphs $G_{k_{\mathsf{\text{min}}}},\cdots ,G_{k_{\mathsf{\text{max}}}}$ with k = k_min,..., k_max respectively and obtains $A_{k_{\mathsf{\text{min}}}},\dots ,A_{k_{\mathsf{\text{max}}}}$ assembly contig sequences after iterative graph condensation and error removal. Each G _k is constructed from the reads whose assigned k-mer size is not more than k and the contigs resulting from G _k−1 constructed in the previous iteration,

$$\underset{j=k_{\mathsf{\text{min}}}}{\overset{k}{\cup }}{R}_j\cup A_{k-1}.$$

(3)

The rationale behind this idea is that those reads that have an assigned k-mer size not more than k should ideally not create any repeats when they are assembled with the k-mer size k. The iterative inclusion of contigs from previous rounds, first introduced in IDBA [12] and later adopted by others [11], is an idea that has already shown merit in improving assembly quality. In Algorithm 2, HyDA is a function that accepts a set of input sequences and an integer k, and returns a set of contigs which are obtained from assembling the input sequences with a k de Bruijn graph.

Algorithm 2 Construction of the de Bruijn graphs

1: function HyDA-Vista(R, K)

2:     $k_{\mathsf{\text{min}}}\leftarrow \underset{r\in R}{\mathsf{\text{min}}}K\left(r\right),k_{\mathsf{\text{max}}}\leftarrow \underset{r\in R}{\mathsf{\text{max}}}K\left(r\right)$

3:     for all k_min ≤ k ≤ k_max do

4:         $R_k\leftarrow \mathrm{0̸}$

5:     end for

6:     for all r ∈ R do

7:         k ← K(r)

8:         R _k ← R _k ∪ {r}

9:     end for

10:     $R^{\prime }\leftarrow \mathrm{0̸}$

11:     $A_{k_{\mathsf{\text{min}}}-1}\leftarrow \mathrm{0̸}$ ▷ assembly contigs

12:     for k := k_min to k_max do

13:         $R^{\prime }\leftarrow R^{\prime }\cup R_k$

14:         $A_k\leftarrow \mathsf{\text{HyDA}}\left(R^{\prime }\cup A_{k-1},k\right)$ ▷ contigs resulting from assembly with HyDA

15:     end for

16:     return $A_{k_{\mathsf{\text{max}}}}$

17: end function

Results

Improved efficiency due to maximal sequence landscape

HyDA-Vista uses the maximal landscape to break the reads into groups by assigning each a value of k. It then uses these groups to build the graph iteratively. This is in contrast to other methods that also iteratively build of the graph; SPAdes [11] uses all the reads at each iteration, and IDBA [12] uses a more complicated approach to remove some subset of reads at each iteration. Thus, one of the main advantages of using the maximal sequence landscape is that it increases the efficiency of building the assembly graph iteratively without degrading assembly quality (see the next subsection for a comparison of the different assemblers). To demonstrate this efficiency experimentally we ran HyDA-Vista with and without the maximal sequence landscape on multicell E. coli (substr. K-12) Illumina data and the E. coli (substr. K-12) reference genome. See Subsection for a description of this dataset. Without the maximal sequence landscape the assembly took 1,414 minutes, and with the maximal sequence landscape the assembly took 822 minutes with 42 number of minutes for building the maximal sequence landscape and assigning the values of k to the reads.

Comparison between competing assemblers and HyDA-Vista on E.Coli

The first data set consists of approximately 27 million paired-end 100 bp reads from multicell E. coli (substr. K-12), generated by the Illumina Genome Analayzer (GA) IIx platform. It was obtained from the NCBI Short Read Archive (accession ERA000206, EMBL-EBI Sequence Read Archive). To assess assembly quality, we aligned the reads to the E. coli reference genome (substr. K-12) using BWA (version 0.5.9) [40] with default parameters. We call a read mapped if BWA outputs an alignment for it and unmapped otherwise. Analysis of the alignments revealed that 98% of the reads mapped to the reference genome, representing an average depth of approximately 600×; An analysis using BLAST against known contaminants revealed that the unmapped reads are attributed to minor contamination of the sample [41]. All reads were error corrected using BayesHammer [42] with default parameters.

KmerGenie [19] predicted 41 to be the optimal k-mer value for this dataset. Therefore, for the assemblers that require a single value of k to be specified (SOAP-denovo, ABySS, HyDA) we used k = 41. HyDA and HyDA-Vista were ran with a cut off of five. All other parameters of SOAPdenovo and ABySS were kept at their default. SPAdes and IDBA were run with their default parameters in single-end mode. Since the input reads were corrected prior to assembly, the reported data for SPAdes is from the "only assembly" stage. IDBA was run with and without error correction, yielding the same statistics as expected.

Table 1 gives the standard assembly statistics of all the assemblies. All statistics in Table 1 were computed by QUAST in default mode [43]. The results demonstrate that HyDA-Vista achieves the best assembly prior to repeat resolution or scaffolding. Note that upon determination of k-mer sizes, all assemblers were run in single-end mode, i.e., ignoring the pairing information, to study only the effect of our k assignment on contigs. HyDA provides a skeletal de Bruijn graph implementation on which other technologies can be developed. HyDA alone does not compete with some of the state of the art assemblers such as SPAdes and IDBA that employ multiple sophisticated technologies. However, empowered by only the maximum landscape information without any other sophisticated technology particularly pairing information, HyDA-Vista increases the N50 and NG50 more than twice (in comparison to HyDA) and outperforms the competing assemblers in all measures (the NG50 of 36 kbps obtained by HyDA increases to 82 kbps by HyDA-Vista).

Table 1 The performance comparison between major assembly tools and HyDA-Vista on the error-corrected standard multicell E.coli dataset (6.2 Gbps, 28 million reads, 100 bp, treated as single-end) using QUAST in default mode [43].

Conclusion

We demonstrated that HyDA-Vista achieves superior performance with respect to standard assembly statistics for ecoli genome before repeat resolution and scaffolding. A crucial aspect of our method is the construction of the maximal sequence landscape for the phylogenetically-close genome, which allows for the optimal k-mer value to be computed for each read. The maximal sequence landscape requires that a pair of substrings be an exact match in order for them to be considered to be the same repetition. An area that warrants future investigation is determining if there is an efficient algorithm for computing the maximal sequence landscape for inexact matches, i.e., the maximal sequence landscape with the condition that two substrings x and y are considered to represent the same repetition (mountain) if the Hamming distance or edit distance between the two is at most d, where d is a parameter in the problem. Another open problem is determining whether the maximal sequence landscape could be constructed with only a suffix array. Since all the data structures (including the suffix array) are constructed in O(n)-time and O(n)-space the order notation of the running time of such a maximal sequence landscape construction algorithm would likely not improve the running time of the existing algorithm by more than a small constant factor. However, the removal of the auxiliary data structures may simplify the algorithm and would be of theoretical interest.