The use of short reads from High Throughput Sequencing (HTS) techniques is now commonplace in de novo assembly. Yet, obtaining contiguous assemblies from short reads is challenging, thus making scaffolding an important step in the assembly pipeline. Different algorithms have been proposed but many of them use the number of read pairs supporting a linking of two contigs as an indicator of reliability. This reasoning is intuitive, but fails to account for variation in link count due to contig features.
We have also noted that published scaffolders are only evaluated on small datasets using output from only one assembler. Two issues arise from this. Firstly, some of the available tools are not well suited for complex genomes. Secondly, these evaluations provide little support for inferring a software’s general performance.
Results
We propose a new algorithm, implemented in a tool called BESST, which can scaffold genomes of all sizes and complexities and was used to scaffold the genome of P. abies (20 Gbp). We performed a comprehensive comparison of BESST against the most popular standalone scaffolders on a large variety of datasets. Our results confirm that some of the popular scaffolders are not practical to run on complex datasets. Furthermore, no single standalone scaffolder outperforms the others on all datasets. However, BESST fares favorably to the other tested scaffolders on GAGE datasets and, moreover, outperforms the other methods when library insert size distribution is wide.
Conclusion
We conclude from our results that information sources other than the quantity of links, as is commonly used, can provide useful information about genome structure when scaffolding.
Recent highthroughput sequencing (HTS) technologies are attractive for de novo assembly projects since they produce millions of short DNAsequences (referred to as reads) at low cost. However, these reads are only a couple of hundred base pairs long making it difficult for an assembler (e.g., [1, 2]) to reconstruct the genome. As a result, the output of an assembly often consists of contigs, i.e., subsets of reads assembled into longer fragments of genomic sequence.
However, HTStechnologies provide protocols for creating read pairs that can be used to increase the contiguity of an assembly. We define a read pair as two reads that are sequenced at a known distance and orientation where the distance between the reads, is referred to as insert size. If the two reads within a read pair belong to different contigs c_{
a
} and c_{
b
}, a link is created between c_{
a
} and c_{
b
}, see Figure 1a. From this link, we can infer a relative order, orientation and distance between c_{
a
} and c_{
b
}.
The process of linking and ordering contigs is called scaffolding. In addition to paired reads, information such as reference sequences of related organisms [3], restriction maps [4] and RNAseq data [5], can be used for contig linking. However, reference based assembly is not applicable to most de novo sequencing projects, restriction maps are often not available, and RNAseq data only have coverage over genes and contains no information about distance between reads which makes contig placement ambiguous. This makes read pair information the most commonly used (and often also the only applicable) source of information for scaffolding.
Unfortunately, scaffolding with read pairs poses challenges: reads may create spurious links because of read errors, heterozygosity and the repeated nature of the genome, and these spurious links make ordering and orientations among the contigs ambiguous. Hence, the scaffolding problem can be summarized as detecting and utilizing the correct links in order to find a consistent ordering and orientation of the contigs. The existing formalizations of scaffolding have been proven to be NPcomplete, but it is still unclear if these formulations, even when finding the optimal solution with respect to the objective, solves the real (i.e. biological) problem. These approaches have focused on structural properties of the graph induced by contig links, with little emphasis on assessing correctness of individual links. Our approach focuses on removing incorrect links and employing sophisticated statistics to evaluate whether linking reads come from the underlying library distribution, or from misalignments. Only in a second step are structural properties used.
The following section discusses the formalization of scaffolding and related work, as well as gives an outline and motivation for our work. Our algorithm, realized in an implementation called BESST (Bias Estimating Stepwise Scaffolding Tool), is presented in detail in the Methods section. The algorithm scales well and is practical on very large and complex genomes, as proved by its use in the Picea abies genome project (20 Gbp) [6]. Furthermore, it excels at scaffolding with wider insert size distributions.
We present an evaluation of BESST against other popular standalone scaffolders on a large variety of datasets from GAGE [7]. Compared to previous assessments of novel scaffolding methodologies, the results obtained from our evaluation allows us to draw conclusions about the general performance of standalone scaffolders to a much higher extent. Another recent extensive evaluation of scaffolding tools is given in [8]. In our study we primarily compare standalone scaffolders because they have access to the same amount of information and are applicable in the same contexts (e.g. scaffolding with mate pair libraries that was not use in the original assembly). Nonetheless, we also include GAGE results on integrated scaffolders.
Our results indicate that no single scaffolder outperforms the others on all datasets although in total, BESST shows the most favorable results among standalone scaffolders. Furthermore, our algorithm outperforms other standalone scaffolders when the library insertsize distribution has a high standard deviation. Although there is wide performance variation around integrated scaffolders, overall, GAGE results demonstrate that AllpathsLG’s assemblies scaffolded with its integrated scaffolder have the highest quality.
The problem
Formalizing Scaffolding
As input for scaffolding we assume a set of contigs produced by an assembler and a number of read pairs from a read pair library that have been aligned to the contigs. These read pairs have an insert size distribution with mean μ and standard deviation σ. By aligning all reads in on we can define the graph as follows: Each contig gives rise to precisely two vertices c_{
i,L
} and c_{
i,R
} in where c_{
i,L
} denotes it’s 5’ end and c_{
i,R
} denotes it’s 3’ end (see Figure 1b). In a read pair, if aligns to precisely one contig c_{
k
} and aligns to precisely one contig c_{
m
}, with k ≠ m, this read pair induces a relative orientation and an approximate distance between c_{
k
} and c_{
m
}. This relationship is represented as an edge e, see Figure 1b. We let V and E denote the set of vertices and edges respectively in . Given , several formulations and methods have been proposed for scaffolding. We will discuss some of them below.
Problem formulations in related works
The scaffolding problem (SP) defined by Huson et al.[9] is a formulation that is commonly referred to. Using their notation, let be defined as above and let n links between two contigs induce a weight n on the edge e between these two contigs. Furthermore, let Φ : V → N be an ordering, orientation and distance map of , that is, an assignment of non negative integer coordinates to the vertices V in that preserves the contig lengths. Given such a mapping instance ϕ, [9] states that an edge e between c_{
i
} and c_{
j
} is consistent if c_{
i
} and c_{
j
} have the correct relative orientation (induced by aligned read pairs), and the distance between c_{
i
} and c_{
j
} is approximately correct. Here, approximately correct means that e suggests a distance between c_{
i
} and c_{
j
} that is less than μ + 3 σ, a heuristically chosen bound. If an edge does not satisfy these constraints, it is called inconsistent. Huson et al. [9] define SP to be the problem of finding a maximum weight consistent edge subset. SP has been used as foundation for a number of other works discussing scaffolding and proposed heuristics for solving it can generally be categorized as either “greedy” or “graphstructure” optimization algorithms.
Greedy algorithms proposed to solve SP include SSPACE and Bambus [10, 11]. SSPACE extends scaffolds in a greedy fashion applying a heuristic stopping criterion. Bambus builds scaffolds greedily with heuristics to remove inconsistent link constraints.
Graphstructure optimization algorithms that have been proposed to solve the SP are for instance: SOPRA [12] formulates a global optimization problem for solving relative contig orientation (exact for simple regions while a simulated annealing approach is employed in more complex regions of the graph). In a second step, readpair distribution is used to determine the relative positions of contigs within a scaffold. If an inconsistency is found in the positioning step, the link causing the inconsistency is removed and the algorithm restarts at the orientation step. OPERA [13] builds scaffolds using the number of inconsistent edges p in a subgraph as a design criterion (the subgraph represents a potential scaffold). By treating p as fixed, they can obtain a polynomial time algorithm to find an optimal (with respect to a given p) solution to their slightly modified version of SP. The algorithm then tries all p starting from p = 0 and stops when a scaffold can be constructed. SLIQ [14] formulates a set of linear inequalities together with majority voting to predict placements of contigs. MIP Scaffolder [15] and GRASS [16] formulate SP as a mixed integer programming problem, but uses different techniques to find a solution. MIP Scaffolder resolves conflicting regions in the obtained MIP solution using heuristics such as removing edges that are stretched or contracted more than a given threshold. GRASS uses an ExpectationMaximization algorithm. The maximization step obtains degrees of penalties on contig links given fixed contig orientations. The penalties are set according to what magnitude the constraints for a link is violated. If a penalty is higher than a given threshold, the penalty of the link is “deactivated”, i.e., its constraints are not considered. The expectation step is used to obtain the expected contig orientation of links given (the “activated”) penalties set in the maximization step. Links that are activated in the final solution are used for scaffolding.
There are advantages and disadvantages with these two classes of methods. Algorithms that are solving a local problem using a greedy approach often have better runtime and scale well on larger genomes but use oversimplified methods to find a solution which may only work for some genomes. Graphstructure optimization methods are instead hindered by their time complexity for finding a solution. The runtime scales poorly and it is difficult to predict if such an algorithm will ever finish on a larger dataset (see section Results).
Additionally, current methods that use insert sizes of paired reads for contig placement are built on false assumptions as we have previously shown [17]. This can complicate scaffolding when libraries with large insertsize variation are used.
Link inconsistency detection
The methods previously described define SP similar to [9] with modifications on how to define a consistent edge. Different heuristics are used between the methods to obtain a solution to SP. Yet, a common denominator is that the number of links supporting an edge is used as an indicator of reliability; edges with many links are preferred and those with few links are avoided. This reasoning is intuitive, but fails to account for variation in link count due to contig features. Firstly, the number of links between two contigs depends on the real (i.e., biologically) distance between the two contigs and on their size [17]. Secondly, in SP we face structural features such as repeated regions, heterozygosity, and chimeric contigs. These features create clusters with reads being misaligned which cannot be seen as individual random events. It is our assumption that the number of random, nonstructural, misalignments caused by, e.g., sequencing errors are almost negligible compared to the structural misalignments. Link count is therefore a poor indicator of edge reliability.
We take a different approach to SP and, instead of link count, evaluate edges based on link statistics. When read pairs are mapped to contigs, are they placed on and connecting contigs in a reasonable way? In other words, we want to answer the question: given an edge e, is the cluster of readpairs forming e coming from the readpair library, or are they a consequence of a structural feature? If these reads together show similar properties as the read pair library we are scaffolding with (e.g., mean, standard deviation), the edge is more likely to be correct.
We propose an algorithm, BESST (BiasEstimating Stepwise Scaffolding Tool), that puts focus on analyzing the scaffold graph in local regions using statistics to filter out spurious edges created by structural errors. BESST starts scaffolding with contigs that meet a length criterion for the library (definition given in section Methods). It then continues with smaller contigs in an optional step. If several different pairedread libraries are used, BESST scaffolds with one library at a time in an increasing order of insert size of the library. Separating contigs with respect to size is mainly due to two reasons: (i) Links between large contigs make gap size estimation more stable (see [17]) giving a more robust statistical analysis. (ii) The gain in speed is significant since correct regions are simple path components in which are found by visiting each edge once, thus, the time complexity is O (E).
Results and discussion
De novo assembly validation is a task as difficult as de novo assembly itself. Recent evaluation efforts like GAGE [7] and Assemblathon [18] encountered several problems in identifying the best assembler. GAGE clearly demonstrated how the same assembler can have completely different performances (e.g., quality) even on similar datasets (e.g., bacterial genomes). This predicament was also supported in recent evaluation efforts [19, 20]. Despite this, as noted by [21], all new published assemblers and scaffolders have been compared to the thenexisting tools highlighting better performances on a specific dataset using some specific metrics. We argue that evaluation of tools should be performed on multiple datasets and/or scenarios to avoid overgeneralization and confirmation bias. For standalone scaffolders without stated dependencies, it is advisable to test on output from several assemblers to investigate overall performance.
We have tried to address the above issues in our evaluation of BESST, using a wide range of different datasets and assemblers. BESST has been compared with three other stateoftheart scaffolders: OPERA, SOPRA, and SSPACE.
Datasets
Evaluation has been performed using the three GAGE datasets [7] which gave us the possibility to evaluate scaffolders on three highly different genomes: Staphylococcus aureus, Rhodobacter spaeroides, and Human chromosome 14 (hereafter referred to as Hs14). All three datasets have been sequenced with high coverage Illumina pairedend (i.e., PEreads) and matepairs (i.e., MPreads) libraries. Moreover each organism has been assembled with up to 8 different assemblers.
GAGE provides high quality MPlibraries with narrow insert size distributions with standard deviation lower than 10% of the mean. However, narrow insert size libraries cannot be obtained in assembly projects where only small amounts of DNA are available. The MP libraries obtained in these cases are wide and the standard deviation can be up to 50% of the mean. BESST uses a technique that works well for larger uncertainties in insert size as this was one of the design assumptions. Therefore we have included the MP library provided in [22] which is characterised by a large variation in insert size. We used picard [23] to estimate the mean and standard deviation of insert size to 2600 and 1250 base pairs respectively. This library will from now on be referred to as the “wide MP” library. An insert size histogram of this distribution is available in Additional file 1: Figure S2.
Evaluation
We scaffolded all 23 available (contig level) GAGEassemblies with BESST v1.0.4.2, and the standalone scaffolders OPERA v1.2, SOPRA v1.4.6, and SSPACEbasic v2.0 using both PE and MP libraries provided by GAGE. Results for assemblerintegrated scaffolders, as computed by GAGE, are also presented, but we primarily compare with the standalone scaffolders because they have access to the same amount of information as BESST and are applicable in similar situations. Note that in GAGE evaluation, Bambus2 was used both for contig and scaffold assembly (with unitigs provided by Celera Assembler).
All scaffolders were run with default parameters (see Additional file 1 for details) on a 1 TB RAM machine equipped with 24 CPUs. Read pairs were mapped to contigs using BWA v0.6.1 for BESST, OPERA, and SOPRA. SSPACEbasic is distributed with Bowtie, thus we used the included version of Bowtie (v0.12.5) for alignments with SSPACE. SSPACE also have a commercial version that supports alignments with BWA. The difference in read alignment method may have an impact on the scaffolding result but we did not investigate this. Out of the 124 scaffolding experiments, 117 successfully terminated within our runtime limit of 48 hours (OPERA and SOPRA were not able to scaffold the Hs14 dataset within this time limit in 3 and 4 cases respectively). Moreover, for the Rhodobacter genome, we also scaffolded the 8 available contiglevel assemblies employing the wide MP library. To summarize, a total of 156 scaffolding experiments have been run, and of these, 149 terminated within the runtime limit and were evaluated.
Each of the 149 results have been evaluated with GAGE validation scripts http://gage.cbcb.umd.edu/results/gagevalidation.tar.gzfor scaffolds, using the available reference sequence. For each assembly, we used GAGE evaluation scripts to compute:
Scaffold errors: number of indels, inversions, relocations, and translocations (as defined by [7]).
Scaffold NG50: size of the longest scaffold such that the sum of the lengths of all scaffolds longer than it is at least half of the (known) reference genome size.
Scaffold Esize: The expected scaffold size at a randomly chosen position on the genome (introduced and defined by [7]). The Esize is calculated as where L_{
c
} is the length of scaffold c and G is the genome length estimated by the sum of all scaffold lengths. Esize is computed similarly for contigs.
Scaffold corrNG50: NG50 after scaffolds have been broken at every position a scaffold error is found.
Scaffold corrEsize: Esize after scaffolds have been broken at every position a scaffold error is found.
Moreover, for each entry, we also compute:
Number of initial contigs and number of produced scaffolds.
Time used by the scaffolder (without considering time required to align reads).
NG50 is a common metric to evaluate an assembly, often offering a good indication of the connectivity as it gives the weighted median contig length. However, the size of one scaffold can be misleading as a measure of the general connectivity of an assembly (as discussed in [7]) Consider, for example, a simple case of two error free assemblies a and b of a 1000 bp genome. If assembly a has one contig of 499 bp and 5 contigs of 100 bp while assembly b has 10 contigs of 100 bp, both will have an NG50 of 100 bp. The measure will therefore not expose the difference in quality between a and b. However, the respective Esizes for assembly a and b are 299 and 100, and thus better capturing the average assembly fragmentation.
Results
Tables 1, 2 and 3 presents the scaffolding performances for high quality libraries provided by GAGE. With the evaluation metrics provided here, no standalone scaffolder is a clear winner (as expected [7, 20]). In general, BESST produces favorable results on all of the organisms. Contrary to the results in [8], SOPRA does not perform well on the metrics provided by GAGE. The results for assemblerintegrated scaffolders, as computed by GAGE, are presented alongside the standalone scaffolders results. There is a large variation in performance of integrated scaffolders but in general, BESST fares well also here. We note that only AllpathsLG has better scaffolded assembly on all three GAGE datasets. Scaffolds from Bambus2 on S. aureus and SGA on Hs14 are two other instances where the integrated scaffolder outperforms the standalone ones.
Table 1
Staphylococcus aureusGAGE data
BESST
OPERA
SOPRA
SSPACE
Integrated scaffolder
Unscaffolded
CorrEsize
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
ABySS
263,4
1
316,7
12
103,4
2
126,3
5
35,3
1
31,4
AllpathsLG
436,4
0
607,4
12
295,5
0
1030,0
1
1136,2
0
90,0
Bambus2
827,3
1
560,0
4
125,2
2
665,7
2
1119,5
0
19,6
MSRCA
744,7
3
302,4
11
117,4
0
781,6
2
999,9
3
50,3
SGA
75,1
0
920,1
3
239,9
6
32,6
2
162,9
1
4,7
SOAPdenovo
346,9
0
333,1
7
227,2
0
286,7
5
229,3
0
68,0
Velvet
204,2
4
236,8
5
154,4
1
162,2
12
194,6
17
48,5
SUM
2898,1
9
3276,6
54
1263,0
11
3085,1
29
The numbers in bold face style indicate the best corrected Esize and number of errors among the standalone scaffolders for each assembly.
Table 2
Rhodobacter sphaeroides, GAGE data
BESST
OPERA
SOPRA
SSPACE
Integrated scaffolder
Unscaffolded
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
ABySS
70,2
13
65,8
20
44,9
17
34,7
4
73,4
3
6,9
AllpathsLG
2005,7
0
852,1
4
425,4
2
1271,9
1
2401,7
0
35,9
Bambus2
1426,0
4
1446,0
8
1469,0
3
789,9
1
1348,4
2
16,2
CABOG
474,0
2
362,6
7
293,4
2
419,1
4
211,3
5
21,5
MSRCA
1757,5
3
573,5
8
138,2
1
1579,8
2
2001,1
5
21,6
SGA
100,5
6
148,3
5
105,7
41
44,9
9
48,0
1
3,2
SOAPdenovo
1551,2
0
841,5
7
1477,1
3
1500,6
3
687,6
0
18,6
Velvet
332,9
2
336,1
10
175,6
11
329,6
6
348,1
19
16,7
SUM
7718,1
30
4626,0
69
4129,2
80
5970,5
30
The numbers in bold face style indicate the best corrected Esize and number of errors among the standalone scaffolders for each assembly.
Table 3
Hs14, GAGE data
BESST
OPERA
SOPRA
SSPACE
Integrated scaffolder
Unscaffolded
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
ABySS
21,6
13
15,8
200


15,3
47
2,8
9
3,1
AllpathsLG
513,6
32
311,0
104
194,9
17
559,0
22
4652,3
45
27,1
Bambus2
88,2
75
61,7
331


99,0
109
157,6
143
6,3
CABOG
421,9
31
349,1
77
234,0
19
411,0
23
347,7
597
30,7
MSRCA
51,3
95




51,9
146
111,9
1068
5,9
SGA
57,2
58
3,5
39
22,2
2253
24,8
42
89,9
19
3,7
SOAPdenovo
94,1
211




75,3
205
99,2
268
9,8
Velvet
35,7
52


75,4
734
22,6
140
26,6
9156
3,0
SUM
1283,6
567
741,0
751
526,5
3023
1259,0
734
The numbers in bold face style indicate the best corrected Esize and number of errors among the standalone scaffolders for each assembly.
High quality datasets where insert size variation does not deviate much from the mean are not always available (see Additional file 1 for a discussion of this). The wide MP library has higher variation, and thus, increases the difficulty in scaffolding by introducing more uncertainty for contig placement. Table 4 shows scaffolding results for the wide MP library. Considering this scenario, BESST is outperforming the other standalone scaffolders having the highest total connectivity whilst giving the fewest errors. OPERA shows a slightly higher connectivity in some cases yet produces 17 times more errors in total. Withstanding the SGA assembly, SOPRA shows few errors in all cases. Yet in most cases, SOPRA also shows extremely low connectivity close to the original contig assembly. Similarly, SSPACE is shown to produce few errors but also struggle with connectivity. As mentioned, larger variation in library insertsize introduces more uncertainty of distance estimates and placement of contigs. Thus, scaffolding becomes more error prone. We believe that our performance here is a consequence of BESST’s ability to infer correct linkstatistics despite wide library distribution.
Table 4
Rhodobacter sphaeroideson GAGE contig assemblies using the wide MP library
BESST
OPERA
SOPRA
SSPACE
Unscaffolded
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
err
CorrEsize (kbp)
ABySS
17,6
5
19,7
51
6,7
4
9,9
1
6,9
AllpathsLG
318,1
0
314,5
1
44,4
0
70,3
2
35,9
Bambus2
460,6
0
267,1
0
93,9
0
137,9
0
16,2
CABOG
199,8
3
97,3
4
22,3
0
33,5
0
21,5
MSRCA
192,1
1
203,7
6
24,0
0
38,1
1
21,6
SGA
4,3
0
13,3
111
3,3
12
6,5
3
3,2
SOAPdenovo
756,7
0
720,5
10
156,2
0
206,4
2
18,6
Velvet
350,2
2
62,3
4
17,8
3
30,9
3
16,7
SUM
2299,3
11
1698,4
187
368,7
19
533,6
12
The numbers in bold face style indicate the best corrected Esize and number of errors among the standalone scaffolders for each assembly. Note that the corrected Esize for SOPRA is slightly less than the corrected contig size in the ABySS assembly. This can happen for low contiguity scaffolded assemblies that contains more bases than the contig assemblies (5,0Mbp and 4,5Mbp respectively on this instance). The difference in number of bases is due to the facts that GAGE evaluation script only compute statistics on contigs and scaffolds that are longer than 200bp. GAGE evaluation script returned an error when computing statistics for seven, two and one scaffolds on SSPACE results of ABySS, CABOG and MSRCA assemblies respectively. On SOPRA and OPERA results, 1 respectively 3 scaffolds of the SGA assembly returned this error. We removed these scaffolds from the evaluation in order to compute the results. In all cases, the scaffolds removed summed up to a total length of less than 110 kbp. Thus, this has a negligible (either positive or negative) effect on Esize computation and an eventual positive effect on the number of errors. Results for BESST contained no scaffolds giving this error.
Tables 5, 6, 7 and 8 illustrates the types of errors that the standalone scaffolders make on the different data sets and Tables 9, 10, 11 and 12 shows runtime of the scaffolders excluding alignment time. An upper bound on runtime was set to 48 hours. BESST and SSPACE present good time scalability in contrast to SOPRA and OPERA which did not finish after 48 hours on 3 and 4 Hs14 instances respectively.
Table 5
Types of errors onS. aureussummed over all assemblies
Assembly
BESST
OPERA
SOPRA
SSPACE
Indels
2
17
2
6
Inversions
6
30
6
16
Translocations
0
0
0
1
Relocations
1
7
3
6
Table 6
Types of errors onRhodobacter sphaeroidessummed over all assemblies
Assembly
BESST
OPERA
SOPRA
SSPACE
Indels
16
28
13
7
Inversions
4
9
7
3
Translocations
2
26
21
14
Relocations
8
6
39
6
Table 7
Types of errors on Hs14 summed over all assemblies
Assembly
BESST
OPERA
SOPRA
SSPACE
Indels
398
149
1062
442
Inversions
163
383
154
289
Translocations
0
0
0
0
Relocations
6
219
1807
3
Table 8
Types of errors onRhodobacter sphaeroideswith wide MP library summed over all assemblies
Assembly
BESST
OPERA
SOPRA
SSPACE
Indels
5
87
4
3
Inversions
3
5
1
0
Translocations
1
9
1
7
Relocations
2
86
13
2
Table 9
Runtime for scaffolders onStaphylococcus aureus
Runtime (hh:mm:ss)
Assembly
BESST
OPERA
SOPRA
SSPACE
ABySS
00:00:40
00:28:47
01:18:24
00:00:26
Allpaths
00:00:25
00:00:47
00:11:56
00:00:20
Bambus2
00:00:26
00:00:49
00:22:11
00:00:21
MSRCA
00:00:26
00:01:05
00:19:21
00:00:21
SGA
00:01:03
00:05:58
04:30:08
00:00:51
SOAPdenovo
00:00:25
00:00:50
00:26:51
00:00:19
Velvet
00:00:27
00:00:53
00:39:25
00:00:21
Table 10
Runtime for scaffolders onRhodobacter sphaeroideswith GAGE data
Runtime (hh:mm:ss)
Assembly
BESST
OPERA
SOPRA
SSPACE
ABySS
00:01:22
00:12:38
01:17:49
00:00:40
AllpathsLG
00:00:32
00:01:13
00:10:35
00:00:27
Bambus2
00:00:33
00:01:38
00:10:42
00:00:25
CABOG
00:00:35
00:00:59
00:11:13
00:00:27
MSRCA
00:00:38
00:01:12
00:18:10
00:00:29
SGA
00:01:45
00:01:35
01:30:44
00:00:45
SOAPdenovo
00:00:33
00:03:48
00:18:08
00:00:27
Velvet
00:00:49
00:01:16
00:36:54
00:00:27
Table 11
Runtime for scaffolders on Hs14 with GAGE data (upper bound time requirement was set to 48 hours)
Runtime (hh:mm:ss)
Assembly
BESST
OPERA
SOPRA
SSPACE
ABySS
00:19:37
00:58:22

00:32:55
AllpathsLG
00:05:06
00:53:24
22:02:09
00:12:55
Bambus2
00:07:43
01:18:06

00:14:26
CABOG
00:04:16
00:16:16
11:50:23
00:08:33
MSRCA
00:11:22


00:15:38
SGA
00:53:42
00:23:18
16:15:22
00:38:42
SOAPdenovo
00:07:50


00:10:46
Velvet
00:10:07

01:27:16
00:15:35
Table 12
Runtime for scaffolders onRhodobacter sphaeroideswith wide mate pair library
Runtime (hh:mm:ss)
Assembly
BESST
OPERA
SOPRA
SSPACE
ABySS
00:00:52
00:04:33
00:14:49
00:00:50
AllpathsLG
00:00:31
00:04:03
00:08:48
00:00:36
Bambus2
00:00:25
00:03:50
00:09:08
00:00:33
CABOG
00:00:31
00:03:10
00:07:17
00:00:37
MSRCA
00:00:30
00:03:46
00:08:29
00:00:39
SGA
00:00:36
00:04:06
00:16:29
00:00:49
SOAPdenovo
00:00:27
00:04:35
00:09:54
00:00:35
Velvet
00:00:33
00:03:58
00:09:38
00:00:40
Conclusion
We proposed a new algorithm, BESST, for the scaffolding problem. This algorithm works well on both small and large datasets. Moreover, we performed a large evaluation of our software against other popular standalone scaffolders. BESST places favorably compared to the other scaffolders on GAGE datasets and outperforms the other methods on libraries with a wide insertsize distribution.
Methods
Scaffolding of larger contigs
BESST works on a graph structure (as defined under Formalizing scaffolding in the Background section). We apply statistics to assess similarity of observed link distribution to the expected link distribution between contigs larger than μ+4σ: a value chosen so that it is very unlikely that a properly mapped read pair will span over a contig of such size. This means that a correctly assembled contig that is not a perfect repeat will have at most one true edge to a neighboring contig of this size. However, the graph structure created for contigs of this size is in practice often far from linear due to e.g. small repeated regions and chimeric regions and that is why we want a way to assess edge quality.
The assessment of edges are realized in a score designed to reflect how reads from a read pair library should be placed on contigs if they were in fact correctly assembled close to each other on the genome. It consists of two parts, a link variation score π_{
σ
} and a link dispersity score π_{
ζ
}, which we present in following subsections.
Link variation score (π_{σ}):
Let c_{
i
},c_{
j
} be two correctly assembled contigs at distance d away from each other on the genome (with d small enough for the readpair library to span). Reads linking c_{
i
} and c_{
j
} follow different distributions depending on the size of d, and a good estimate of d can be calculated using MLestimation with the tool GapEst introduced in [17]. Here, we go one step further and answer the question: Given μ,σ and obtained from links observed over c_{
i
} and c_{
j
}, what should the standard deviation of these links be? We denote this quantity with σ_{
od
} (standard deviation of observations given a gap size) to be consistent with the notation used for GapEst. Theorem 1 gives the theoretical expected value of σ_{
od
} which is dependent on the read length r and the length of the longer and shorter contig giving rise to the gap (denoted c_{
min
},c_{
max
}).
Theorem 1.Given μ,σ and d, σ_{
od
}is given by
where g (d) and q (d) are defined in Additional file 1.
This quantity is a measure of how far observed distances are from the theoretical distance. Note that 0≤π_{
σ
}≤1.
Link dispersity score (π_{
ζ
})
The other part of the scoring function is an indicator of how well dispersed links are over the contigs they are connecting, given an estimated gap d between them. This dispersity is scored by the two sample KolmogorovSmirnov statistic [24] that gives a measure of difference in distribution between two independent samples. Letting observations on a contig be a sample, the independence between samples comes from the fact that the aligner has no information about which contigs lie close together, thus links between contigs can be seen as two independent alignment events.
By observing if link distribution is similar for two linked contigs, we can detect abnormal edges that might come from one or several smaller repeats residing within a contig (see Figure 2b). We call this score the dispersity score and it is calculated as follows. Before testing, translation and reflection of the observations are necessary (see Figure 2a).
Let n read pairs link two contigs c_{1} and c_{2} where c_{1} is of length m. Recall that are the ith observations on c_{1} and c_{2} respectively. Let , where is the mean observation on c_{1} and similar for . Samples y^{
i
} are therefore the transformed observations as seen in Figure 2. The empirical distribution of samples is given by
where is the indicator function. Let and be the two empirical distributions for the samples on c_{1} respectively c_{2}, the twosample KS statistic of the observations is then obtained by
It follows that D_{
n
}∈ [0,1] since this quantity measures the largest distance between the two empirical cumulative distributions. The similarity score is defined as π_{
ζ
}= 1  D_{
n
}.
Scoring edges:
To sum up the two previous subsections, we first derived a score for the ratio of the expected to the observed standard deviation of distance for links spanning a gap. Secondly, we gave a similarity score of expected to observed dispersity of links spanning a gap using the twosample KS statistic. The total score of a gapedge is defined as
By definition, we have 0 ≤ π ≤ 2. We have used the heuristic cutoff 0.5 which means that if any of the two quantities deviates more than twice from the assumed value, the score is set to 0, i.e., the edge is discarded from the graph. π is only calculated on edges where c_{
i
},c_{
j
}>μ + 4 σ. That is, any vertex that has more than 2 neighbors in this subgraph is considered to be involved in a region with linking errors since by the constraints of the contig lengths, the library should not be able to link more than one such contig. The score is used to choose the edge with links that best resemble the library. If the two highest scores in such a region are close to each other (their ratio higher than 0.9 set heuristically), we chose to not make a decision. This can for instance occur from larger repeated contigs. This approach finds a mapping ϕ on , representing a scaffolding, in time.
Note about scoring edges
It might be inviting to start using pvalues that can be estimated from the distributions we have defined. This would lead to a statistical test for keeping or discarding edges in the scaffold graph. However, this is not suitable for the problem in hand. Fewer links between contigs gives more uncertainty (leading to volatile pvalues) and can lead to inability to discard many edges with low link support. Edges with many links would also be sensitive to smaller aberrations by increased sensitivity of statistical testing with larger sample sizes. In the case of multiple edges from a contig, we want to compare edges to see which observations have matching distributions. Comparing significance levels of pvalues to make this decision is bad practice since pvalues are nonlinear transforms of data that should only be interpreted under the null hypothesis. That is, the pvalue is a measure of evidence; it is not an estimate of effect size. Looking at the similarity ratio for π_{
σ
} and the KS statistic for π_{
s
} provides a measure that is robust to the number of links and can be used to measure the fit of data when a decision is needed.
Including smaller contigs
Small contigs are defined as having a length less than μ + 4 σ, i.e., all contigs not treated in the previous section. This limit varies with respect to the current library and is used to efficiently create linear scaffolds (as explained in previous section). There are limitations when scaffolding with contigs of size over a particular threshold. Firstly, one will have gaps in these scaffolds where shorter contigs could be placed. Secondly, several small contigs can occur between two large scaffolds making read pairs unable to link them together. We address this issue as follows.
In graph theory, a simple path in is a path without repeated vertices. Let a connection between two large contigs c_{
a
} and c_{
b
} in be defined as a simple path starting at c_{
a
} and ending at c_{
b
} with the rest of its vertices as small contigs. For a connection γ consisting of a sequence of n contigs {c_{1},…,c_{
n
}}, we define g(c_{
i
}) to be the number of links that goes from c_{
i
} to any other contig c_{
j
}∈ γ. Similarly, let b(c_{
i
}) be the number of links that go from c_{
i
} to any other contig c_{
j
}∉ γ. The notation of g and b are chosen to indicate “good” links and “bad” links respectively, for a given connection γ. The score of γ is defined as
i.e., the sum of differences of good and bad links for all contigs belonging in γ. An example region of connections is shown in Figure 3. To find connections between two large contigs c_{
a
},c_{
b
}, we use breadthfirst search in . If more than one connection is found, the highest scoring one is chosen if the score is positive.
If c_{
a
} and c_{
b
} is within an already created scaffold, the algorithm will look for paths with length less than d + 2 σ, where 2 σ allows for uncertainty of the estimate of d. If c_{
a
} and c_{
b
} are not within an already created scaffold, there is no distance constraint on the length of the connection. For dense regions in , there can be an exponential blowup in the number of possible connections. We have set a threshold limiting the breadthfirst search to 1000 iterations. This restricting threshold is motivated by two reasons. Firstly, dense regions are likely to be caused by spurious edges. Thus, paths created in these regions will often have a negative score. Secondly, we have seen from large datasets that correct connections tend to be short. This makes higher ordered layers in the breadth first search contain connections with negative score or paths that does not lead to a valid end vertex (to form a connection).
All connections with a positive score are used to improve the contiguity of the scaffolds. First, gaps within existing scaffolds are filled if a connection with a positive score has been found. In a second step, connections between scaffolds are considered. The extension starts with the highest scoring connection first and proceeds in a descending order of the score. If a contig is found in multiple connections with positive scores, it is only used in the one with the highest score.
Using multiple libraries
If given multiple libraries, BESST uses these libraries in an increasing order of library insert size. Scaffolds created in earlier steps are seen as contigs for the next library.
Implementation
BESST source code is available under the GNU GPL v3 license. It is implemented in Python using Networkx [25] graph library to represent the scaffolding graph and pysam [26] for parsing BAM files. As input BESST takes contigs in a FASTA file and the alignments of the paired reads to the contigs as sorted BAM files. BESST can use several pairedread libraries with different insert sizes. The main output consists of scaffolds in a FASTA file. If several libraries are used, a scaffold FASTA file is given as output in each scaffolding step. The output also consists of AGP and GFFfiles that contain information about the scaffolds, such as position of each contig in the scaffold and length of the gaps. The contigs that were classified as repeats are output in a separate FASTA file.
Preprocessing of
When initializing , BESST computes different statistics of the read library such as mean and standard deviation of insert size (μ,σ) and of the coverage (μ_{
c
},σ_{
c
}). Links with inconsistent insert size defined as o_{1}+ o_{2}>μ + lσ are not considered in the scaffolding since they are likely to be placed on chimeric contigs or misaligned. Here, l is a user defined constant which defaults to 6.
BESST removes the contigs that, based on coverage, behave as repeats. A contig c_{
i
} is classified as a repeat if coverage of c_{
i
} is larger than max {2 μ_{
c
},μ_{
c
}+ tσ_{
c
}}, where t is calculated with the same principle as computing k for π_{
ζ
}.
Data and optional parameters as input to BESST
BESST uses alignments of paired reads to contigs in format of sorted BAM files. A read aligner such as BWA or Bowtie [27, 28] can be used to map the paired reads in forward reverse mode. We use those read pairs whose both ends map to a unique position in the collection of contigs.
Several parameters for BESST can optionally be set on the command line:
μ and σ can be specified instead of being computed internally. This can be good if the assembly is very fragmented.
Minimum number of links needed to create an edge (with 5 as default value).
Coverage cutoff for repeat identification.
Duplicate read remover (based on identical map positions of both fragments in a paired read).
The inclusion of small contigs can be inactivated.
This work was in part funded by the Swedish Research Council (grant 20104634) and the Spruce Genome Project grant from the Knut and Alice Wallenberg Foundation.
We thank Lex Nederbragt and Ole Kristian Tørressen for comments and suggestions that greatly improved the manuscript and Carly Schott for help with editing.
Authors’ Affiliations
(1)
Science for Life Laboratory, School of Computer Science and Communication, KTH Royal Institute of Technology
(2)
Science for Life Laboratory, Department of Biochemistry and Biophysics, Stockholm University
(3)
Science for Life Laboratory, School of Biotechnology, Division of Gene Technology, KTH Royal Institute of Technology
(4)
Swedish eScience Research Centre (SeRC), Department of Numerical Analysis and Computer Science, Stockholm University
References
Zerbino DR, Birney E: Velvet: algorithms for de novo short read assembly using de Bruijn graphs.Genome Res 2008,18(5):821–829.PubMed CentralPubMedView Article
MacCallum I, Przybylski D, Gnerre S, Burton J, Shlyakhter I, Gnirke A, Malek J, McKernan K, Ranade S, Shea TP, Williams L, Young S, Nusbaum C, Jaffe DB: ALLPATHS 2: small genomes assembled accurately and with high continuity from short paired reads.Genome Biol 2009,10(10):103.View Article
Nagarajan N, Read TD, Pop M: Scaffolding and validation of bacterial genome assemblies using optical scaffolding and validation of bacterial genome assemblies using optical restriction maps.Bioinformatics 2008, 24:1229–1235.PubMed CentralPubMedView Article
Mortazavi A, Schwarz E, Williams B, Schaeffer L, Antoshechkin I, Wold B, Sternberg P: Scaffolding a Caenorhabditis nematode genome with RNAseq.Genome Res 2010,20(12):1740–1747.PubMed CentralPubMedView Article
Nystedt B, Street NR, Wetterbom A, Zuccolo A, Lin YC, Scofield DG, Vezzi F, Delhomme N, Giacomello S, Alexeyenko A, Vicedomini R, Sahlin K, Sherwood E, Elfstrand M, Gramzow L, Holmberg K, Hällman J, Keech O, Klasson L, Koriabine M, Kucukoglu M, Käller M, Luthman J, Lysholm F, Rilakovic N, Ritland C, Sena J, Niittylä T, et al.: The Norway spruce genome sequence and conifer genome evolution.Nature 2013,497(7451):579–584.PubMedView Article
Salzberg SL, Phillippy AM, Zimin A, Puiu D, Magoc T, Koren S, Treangen TJ, Schatz MC, Delcher AL, Roberts M, Marcais G, Pop M, Yorke JA: GAGE: a critical evaluation of genome assemblies and assembly algorithms.Genome Res 2012,22(3):557–567.PubMed CentralPubMedView Article
Hunt M, Newbold C, Berriman M, Otto T: A comprehensive evaluation of assembly scaffolding tools.Genome Biol 2014,15(3):42.View Article
Huson DH, Reinert K, Myers EW: The greedy pathmerging algorithm for contig scaffolding.J ACM 2002,49(5):603–615.View Article
Boetzer M, Henkel CV, Jansen HJ, Butler D, Pirovano W: Scaffolding preassembled contigs using SSPACE.Bioinformatics 2011,27(4):578–579.PubMedView Article
Pop M, Kosack DS, Salzberg SL: Hierarchical scaffolding with Bambus.Genome Res 2004,14(1):149–159.PubMed CentralPubMedView Article
Dayarian A, Michael TP, Sengupta AM: SOPRA: scaffolding algorithm for paired reads via statistical optimization.BMC Bioinformatics 2010, 11:345.PubMed CentralPubMedView Article
Roy RS, Chen KC, Sengupta AM, Schliep A: SLIQ: Simple linear inequalities for efficient contig scaffolding.J Comput Biol 2012,19(10):1162–1175.PubMedView Article
Salmela L, Mäkinen V, Välimäki N, Ylinen J, Ukkonen E: Fast scaffolding with small independent mixed integer programs.Bioinformatics 2011,27(23):3259–3265.PubMed CentralPubMedView Article
Gritsenko AA, Nijkamp JF, Reinders MJ, de Ridder D: GRASS: a generic algorithm for scaffolding nextgeneration sequencing assemblies.Bioinformatics 2012,28(11):1429–1437.PubMedView Article
Sahlin K, Street N, Lundeberg J, Arvestad L: Improved gap size estimation for scaffolding algorithms.Bioinformatics 2012,28(17):2215–2222.PubMedView Article
Earl D, Bradnam K, St John J, Darling A, Lin D, Fass J, Yu HO, Buffalo V, Zerbino DR, Diekhans M, Nguyen N, Ariyaratne PN, Sung WK, Ning Z, Haimel M, Simpson JT, Fonseca NA, Birol I, Docking TR, Ho IY, Rokhsar DS, Chikhi R, Lavenier D, Chapuis G, Naquin D, Maillet N, Schatz MC, Kelley DR, Phillippy AM, Koren S, et al.: Assemblathon 1: a competitive assessment of de novo short read assembly methods.Genome Res 2011,21(12):2224–2241.PubMed CentralPubMedView Article
Vezzi F, Narzisi G, Mishra B: Featurebyfeature–evaluating de novo sequence assembly.PLoS ONE 2012,7(2):31002.View Article
Vezzi F, Narzisi G, Mishra B: Reevaluating assembly evaluations with feature response curves: GAGE and assemblathons.PLoS ONE 2012,7(12):52210.View Article
Miller JR, Koren S, Sutton G: Assembly algorithms for nextgeneration sequencing data.Genomics 2010,95(6):315–327.PubMed CentralPubMedView Article
Ribeiro FJ, Przybylski D, Yin S, Sharpe T, Gnerre S, Abouelleil A, Berlin AM, Montmayeur A, Shea TP, Walker BJ, Young SK, Russ C, Nusbaum C, MacCallum I, Jaffe DB: Finished bacterial genomes from shotgun sequence data.Genome Res 2012,22(11):2270–2277.PubMed CentralPubMedView Article
Kolmogorov AN: Sulla determinazione empirica di una legge di distribuzione (On the empirical determination of a distribution law).Giornale dell’Istituto Italiano degli Attuari 1933, 4:83–91.
Li H, Durbin R: Fast and accurate longread alignment with BurrowsWheeler transform.Bioinformatics 2010,26(5):589–595.PubMed CentralPubMedView Article
Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memoryefficient alignment of short DNA sequences to the human genome.Genome Biol 2009,10(3):25.View Article
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.