Integrating Hi-C links with assembly graphs for chromosome-scale assembly

Long-read sequencing and novel long-range assays have revolutionized de novo genome assembly by automating the reconstruction of reference-quality genomes. In particular, Hi-C sequencing is becoming an economical method for generating chromosome-scale scaffolds. Despite its increasing popularity, there are limited open-source tools available. Errors, particularly inversions and fusions across chromosomes, remain higher than alternate scaffolding technologies. We present a novel open-source Hi-C scaffolder that does not require an a priori estimate of chromosome number and minimizes errors by scaffolding with the assistance of an assembly graph. We demonstrate higher accuracy than the state-of-the-art methods across a variety of Hi-C library preparations and input assembly sizes. The Python and C++ code for our method is openly available at https://github.com/machinegun/SALSA.


Introduction
Genome assembly is the process of reconstructing a complete genome sequence from significantly shorter sequencing reads. Most genome projects rely on whole genome shotgun sequencing which yields an oversampling of each genomic locus. Reads originating from the same locus are identified using assembly software, which can use these overlaps to reconstruct the genome sequence [1,2]. Most approaches are based on either a de Bruijn [3] or a string graph [4] formulation. Repetitive sequences exceeding the sequencing read length [5] introduce ambiguity and prevent complete reconstruction. Unambiguous reconstructions of the sequence are output as "unitigs" (or often "contigs"). Ambiguous reconstructions are output as edges linking unitigs. Scaffolding utilizes long-range linking information such as BAC or fosmid clones [6,7], optical maps [8][9][10], linked reads [11][12][13], or chromosomal conformation capture [14] to order and orient contigs. If the linking information spans large distances on the chromosome, the resulting scaffolds can span entire chromosomes or chromosome arms.
Hi-C is a sequencing-based assay originally designed to interrogate the 3D structure of the genome inside a cell nucleus by measuring the contact frequency between all pairs of loci in the genome [15]. The contact frequency between a pair of loci strongly correlates with the one-dimensional distance between them. Hi-C data can provide linkage information across a variety of length scales, spanning tens of megabases. As a result, Hi-C data can be used for genome scaffolding. Shortly after its introduction, Hi-C was used to generate chromosomescale scaffolds [16][17][18][19][20].
LACHESIS [16] is an early method for Hi-C scaffolding that first clusters contigs into a user-specified number of chromosome groups and then orients and orders the contigs in each group independently to generate scaffolds. Thus, the scaffolds inherit any assembly errors present in the contigs. The original SALSA1 [21] method first corrects the input assembly, using a lack of Hi-C coverage as evidence of error. It then orients and orders the corrected contigs to generate scaffolds. Recently, the 3D-DNA [20] method was introduced and demonstrated on a draft assembly of the Aedes aegypti genome. 3D-DNA also corrects the errors in the input assembly and then iteratively orients and orders contigs into a single megascaffold. This megascaffold is then broken, identifying chromosomal ends based on the Hi-C contact map.
There are several shortcomings common across currently available tools. They are sensitive to input assembly contiguity and Hi-C library variations and require tuning of parameters for each dataset. Inversions are common when the input contigs are short, as orientation is determined by maximizing the interaction frequency between contig ends across all possible orientations [16]. When contigs are long, there are few interactions spanning the full length of the contigs, making the true orientation apparent from the higher weight of links. However, in the case of short contigs, there are interactions spanning the full length of the contig, making the true orientation have a similar weight to incorrect orientations. Biological factors, such as topologically associated domains (TADs), also confound this analysis [22]. SALSA1 [21] addressed some of these challenges, such as not requiring the expected number of chromosomes beforehand and correcting assemblies before scaffolding them with Hi-C data. We showed that SALSA1 worked better than the most widely used method, LACHESIS [16]. However, SALSA1 often did not generate chromosome-sized scaffolds. The contiguity and correctness of the scaffolds depended on the coverage of Hi-C data and required manual data-dependent parameter tuning. Building on this work, SALSA2 does not require manual parameter tuning and is able to utilize all the contact information from the Hi-C data to generate near optimal sized scaffolds permitted by the data using a novel iterative scaffolding method. In addition to this, SALSA2 enables the use of an assembly graph to guide scaffolding, thereby minimizing errors, particularly orientation errors.
SALSA2 is an open source software that combines Hi-C linkage information with the ambiguous-edge information from a genome assembly graph to better resolve contig orientations. We propose a novel stopping condition, which does not require an a priori estimate of chromosome count, as it naturally stops when the Hi-C information is exhausted. We show that SALSA2 produces fewer orientation, ordering, and chimeric errors across a wide range of assembly contiguities. We also demonstrate its robustness to different Hi-C libraries with varying levels of intra-chromosomal contact frequencies. When compared to 3D-DNA, SALSA2 generates more accurate scaffolds across most conditions. To our knowledge, this is the first method to leverage assembly graph information for scaffolding Hi-C data. Fig 1(A) shows the overview of the SALSA2 pipeline. SALSA2 begins with a draft assembly generated from long reads such as Pacific Biosciences [23] or Oxford Nanopore [24]. SALSA2 requires the contig sequences and, optionally, a GFA-formatted assembly graph [25] representing the ambiguous reconstructions. Hi-C reads are aligned to the contig sequences, and contigs are optionally split in regions lacking Hi-C coverage. A hybrid scaffold graph is constructed using both ambiguous edges from the GFA and edges from the Hi-C reads, scoring edges according to a "best buddy" scheme. Scaffolds are iteratively constructed from this graph Linkage information obtained from the alignment of Hi-C reads to the assembly. Arrows denote contigs and arcs between arrows denote the inferred linking information from Hi-C reads. Thickness of arcs denote the weight on the Hi-C edge. Thicker edge indicates higher edge weight implied by Hi-C reads (C) Assembly graph obtained from the assembler, where arrows are contigs and arcs denote overlap between contigs(D) Hybrid scaffold graph constructed from the links obtained from the Hi-C read alignments and the overlap graph. Solid edges indicate the linkages between different contigs and dotted edges indicate the links between the ends of the same contig. B and E denote the start and end of contigs, respectively. The E-E edge between blue and red contigs is dashed as this particular orientation between them is not supported by assembly graph, but rather B-E edge is supported. We ignore this dotted edge while computing maximal matching (E) Maximal weighted matching obtained from the graph using a greedy weighted maximum matching algorithm. The numbering of the edges indicates the order in which they were added to the graph. No more solid edges can be added to the matching as it would assign more than one edge to already matched nodes. (F) Edges between the ends of same contigs are added back to the matching to obtain final scaffolds. using a greedy weighted maximum matching. A mis-join detection step is performed after each iteration to check if any of the joins made during this round are incorrect. Incorrect joins are broken and the edges blacklisted during subsequent iterations. This process continues until the majority of joins made in the prior iteration are incorrect. This provides a natural stopping condition, when accurate Hi-C links have been exhausted. Below, we describe each of the steps in detail.

Hi-C library preparation
Hi-C methods first crosslink a sample (cells or tissues) to preserve the genome conformation. The crosslinked DNA is then digested using multiple restriction enzymes (targeting in this case the restriction sites GATC and GANTC). The single-stranded 5'-overhangs are then filled in causing digested ends to be labeled with a biotinylated nucleotide. Next, spatially proximal digested ends of DNA are ligated, preserving both short-and long-range DNA contiguity. The DNA is then purified and sheared to a size appropriate for Illumina short-read sequencing. After shearing, the biotinylated fragments are enriched to assure that only fragments originating from ligation events are sequenced in paired-end mode via Illumina sequencers to inform DNA contiguity.

Read alignment
Hi-C paired end reads are aligned to contigs using the BWA aligner [26](parameters: -t 12 -B 8) as single end reads. First, the reads mapping at multiple locations are ignored as they can cause ambiguities while scaffolding. Reads which align across ligation junctions are chimeric and are trimmed to retain only the start of the read which aligns prior to the ligation junction. After filtering the chimeric reads, the pairing information is restored. Any PCR duplicates in the paired-end alignments are removed using Picard tools [27]. Read pairs aligned to different contigs are used to construct the initial scaffold graph. The suggested mapping pipeline is available at http://github.com/ArimaGenomics/mapping_pipeline.

Contig correction
As any assembly is likely to contain mis-assembled sequences, SALSA2 uses the physical coverage of Hi-C pairs to identify suspicious regions and break the sequence at the likely point of mis-assembly. We define the physical coverage of a Hi-C read pair as the region on the contig spanned by the start of the leftmost fragment and the end of the rightmost fragment. A drop in physical coverage indicates a likely assembly error. In SALSA1, contigs are split when a fixed minimum coverage threshold is not met. A drawback of this approach is that coverage can vary, both due to sequencing depth and variation in Hi-C link density. Fig 2 sketches the new contig correction algorithm implemented in SALSA2. Instead of the single coverage threshold used in SALSA1, a set of suspicious intervals is found with a sweep of thresholds. For a sweep of thresholds, we find the continuous stretches of regions which have lower physical coverage. Note that there can be multiple intervals for a particular threshold that have multiple stretches of low coverage. In such case, we only consider the interval of the maximum size. These intervals denote the regions of potential misassembly on the contig. Using the collection of these intervals as an interval graph, we find the maximal clique, which is the maximal set of intervals intersecting at any location along the contig. This maximal clique represents the region of the contig which had low coverage for the majority of the tested cutoffs. This can be done in O(NlogN) time, where N is the number of intervals. For a maximal clique, the region between the start and end of the smallest interval in the clique is flagged as a mis-assembly and the contig is split into three pieces-the sequence to the left of the region, the junction region itself, and the sequence to the right of the region. The intuition behind choosing the smallest interval is to accurately pinpoint the location of assembly error. Note that this algorithm finds only one misassembly per contig. For more rigorous misassembly detection, same algorithm can be run multiple times on each contig until no more drops in physical coverage are found.

Assembly graph construction
For our experiments, we use the unitig assembly graph produced by Canu [28] (Fig 1(C)), as this is a more conservative assembly output than contig sequences that represent various traversals of this graph. SALSA2 requires only a GFA format [  Assembly graph based scaffolding

Scaffold graph construction
The scaffold graph is defined as G(V, E), where nodes V are the ends of contigs and edges E are derived from the Hi-C read mapping ( Fig 1B). The idea of using contig ends as nodes is similar to that used by the string graph formulation [4].
Modeling each contigs as two nodes allows a pair of contigs to have multiple edges in any of the four possible orientations (forward-forward, forward-reverse, reverse-forward, and reverse-reverse). The graph then contains two edge types-one explicitly connects two different contigs based on Hi-C data, while the other implicitly connects the two ends of the same contig.
As in SALSA1, we normalize the Hi-C read counts by the frequency of restriction enzyme cut sites in each contig. This normalization reduces the bias in the number of shared read pairs due to the contig length as the number of Hi-C reads sequenced from a particular region are proportional to the number of restriction enzyme cut sites in that region. For each contig, we denote the number of times a cut site appears as C(V). We define edges weights of G as: is the number of Hi-C read pairs mapped to the ends of the contigs u and v. By the ends, we mean the first and second half of the contig if divided at the midpoint along its length.
We observed that the globally highest edge weight does not always capture the correct orientation and ordering information due to variations in Hi-C interaction frequencies within a genome. To address this, we defined a modified edge ratio, similar to the one described in [20], which captures the relative weights of all the neighboring edges for a particular node.
The best buddy weight BB(u, v) is the weight W(u, v) divided by the maximal weight of any edge incident upon nodes u or v, excluding the (u, v) edge itself. Computing best buddy weight naively would take O(|E| 2 ) time. This is computationally prohibitive since the graph, G, is usually dense. If the maximum weighted edge incident on each node is stored with the node, the running time for the computation becomes O(|E|). We retain only edges where BB(u, v) > 1. This keeps only the edges that are the best incident edge on both u and v. Once used, the edges are removed from subsequent iterations. Thus, the most confident edges are used first but initially low-scoring edges can become best in subsequent iterations.
For the assembly graph, we define a similar ratio. Since the edge weights are optional in the GFA specification and do not directly relate to the proximity of two contigs on the chromosome, we use the graph topology to establish this relationship. Let � u denote the reverse complement of the contig u. Let σ(u, v) denote the length of shortest path between u and v. For each edge (u, v) in the scaffold graph, we find the shortest path between contigs u and v in every possible orientation, that is, σ(u, v), sðu; � vÞ, sð� u; vÞ and sð� u; � vÞ. With this, the score for a pair of contigs is defined as follows: where x and y are the orientations in which u and v are connected by a shortest path in the assembly graph. Essentially, Score(u, v) is the ratio of the length of the second shortest path to the length of the shortest path in all possible orientations. Once again, we retain edges where Score(u, v) > 1. If the orientation implied by the assembly graph differs from the orientation implied by the Hi-C data, we remove the Hi-C edge and retain the assembly graph edge ( Fig 1D). Computing the score graph requires |E| shortest path queries, yielding total runtime of O(|E| � (|V| + |E|)) since we do not use the edge weights.

Contig layout
Once we have the hybrid graph, we lay out the contigs to generate scaffolds. Since there are implicit edges in the graph G between the beginning and end of each contig, the problem of computing a scaffold layout can be modeled as finding a weighted maximum matching in a general graph, with edge weights being our ratio weights. In a weighted maximum matching, a set of edges from a graph is chosen in such a way that they have no endpoints common and the sum of edge weights is maximized. In the case of scaffolding, a maximum weighted matching implies a layout of contigs, where no end can be used twice, that is maximally consistent with the data being used for scaffolding (Hi-C in our case). If we find the weighted maximum matching of the non-implicit edges (that is, edges between different contigs) in the graph, adding the implicit edges to this matching would yield a complete traversal. However, adding implicit edges to the matching can introduce a cycle. Such cycles are prevented by removing the lowest-weight non-implicit edge. Computing a maximal matching takes O(|E||V| 2 ) time [31]. We iteratively find a maximum matching in the graph by removing nodes found in the previous iteration. Using the optimal maximum matching algorithm this would take O(|E||V| 3 ) time, which would be extremely slow for large graphs. Instead, we use a greedy maximal matching algorithm which is guaranteed to find a matching within 1/2-approximation of the optimum [32]. The greedy matching algorithm takes O(|E|) time, thereby making the total runtime O(|V||E|). The algorithm for contig layout is sketched in Algorithm 1.

Iterative mis-join correction
Since the contig layout is greedy, it can introduce errors by selecting a false Hi-C link which was not eliminated by our ratio scoring. These errors propagate downstream, causing large chimeric scaffolds and chromosomal fusions. We examine each join made within all the scaffolds in the last iteration for correctness. Any join with low spanning Hi-C support relative to the rest of the scaffold is broken and the links are blacklisted for further iterations.
We compute the physical coverage spanned by all read pairs aligned in a window of size w around each join. For each window, w, we create an auxiliary array, which stores −1 at position i if the physical coverage is greater than some cutoff δ and 1, otherwise. We then find the maximum sum subarray in this auxiliary array, since it captures the longest stretch of low physical coverage. If the position being tested for a mis-join lies within the region spanned by the maximal clique generated with the maximum sum subarray intervals for different cutoffs (Fig 2), the join is marked as incorrect. The physical coverage can be computed in O(w + N) time, where N is the number of read pairs aligned in window w. The maximum sum subarray computation takes O(w) time. If K is the number of cutoffs(δ) tested for the suspicious join finding, then the total runtime of mis-assembly detection becomes O(K(N + 2 � w)). The parameter K controls the specificity of the mis-assembly detection, thereby avoiding false positives. The algorithm for mis-join detection is sketched in Algorithm 2. When the majority of joins made in a particular iteration are flagged as incorrect by the algorithm, SASLA2 stops scaffolding and reports the scaffolds generated in the penultimate iteration as the final result.

Algorithm 2 Misjoin detection and correction algorithm
Cov: Physical coverage array for a window size w around a scaffold join at position p on a scaffold A: Auxiliary array I: Maximum sum subarray intervals Break the scaffold at position p end if

Dataset description
We created artificial assemblies, each containing contigs of same size, by splitting the GRCh38 [33] reference into fixed-sized contigs of 200 to 900 kbp. This gave us eight assemblies. The assembly graph for each input was built by adding edges for any adjacent contigs in the genome. So the simulated assembly graph was linear with edges between two adjacent contigs for each contig in the graph.
For real data, we use the recently published NA12878 human dataset sequenced with Oxford Nanopore [34] and assembled with Canu [28]. We use a Hi-C library from Arima Genomics (Arima Genomics, San Diego, CA) sequenced to 40x Illumina coverage (SRX3651893). Table 1 shows the statistics for this library. We compare results with the Table 1. Hi-C library statistics for different datasets used in this paper. Mapped read pairs denote the total number of Hi-C read pairs aligned before mapping quality filtering. Intra-contig read pairs account for the read pairs where both the reads align to same contig and inter-contig read pairs account for the read pairs where both reads align to different contigs.

Library
Total read pairs Mapped read pairs Intra-contig pairs Inter-contig pairs Assembly graph based scaffolding original SALSA(commit-833fb11), SALSA2 with and without the assembly graph input(commit-68a65b4), and 3D-DNA (commit-3f18163). We did not compare our results with LACHESIS because it is no longer supported and is outperformed by 3D-DNA [20]. SALSA2 was run using default parameters, with the exception of graph incorporation, as listed. For 3D-DNA, alignments were generated using the Juicer alignment pipeline [35] with defaults (-m haploid -t 15000 -s 2), except for mis-assembly detection, as listed. A genome size of 3.2 Gbp was used for contiguity statistics for all assemblies.
For evaluation, we also used the GRCh38 reference to define a set of true and false links from the Hi-C graph. We mapped the assembly to the reference with MUMmer3.23 (nucmer -c 500 -l 20) [36] and generated a tiling using MUMmer's show-tiling utility. For this "true link" dataset, any link joining contigs in the same chromosome in the correct orientation was marked as true. This also gave the true contig position, orientation, and chromosome assignment. We masked sequences in GRCh38 that matched known structural variants from a previous assembly of NA12878 [37] to avoid counting true variations as scaffolding errors.

Evaluation on simulated contigs
Assembly correction. We simulated assembly error by randomly joining 200 pairs of contigs from each simulated assembly. All erroneous joins were made between contigs that were more than 10 Mbp apart or were assigned to different chromosomes in the reference. The remaining contigs were unaltered. We then aligned the Arima-HiC data and ran our assembly correction algorithm. When the algorithm marked a mis-join within 20 kbp of a true error we called it a true positive, otherwise we called it a false positive. Any unmarked error was called a false negative. The average sensitivity over all simulated assemblies was 77.62% and the specificity was 86.13%. The sensitivity was highest for larger contigs (50% for 200 kbp versus more than 90% for contigs greater than 500 kbp, S1 Table) implying that our algorithm is able to accurately identify errors in large contigs, which can have a negative impact on the final scaffolds if not corrected. Although we used a cutoff of 20 kbp to evaluate sensitivity and specificity, most of the predicted locations of misassembly were within 5 kbp from the true misassembly location (S2 Fig). Scaffold mis-join validation. As before, we simulated erroneous scaffolds by joining contigs which were not within 10 Mbp in the reference or were assigned to different chromosomes. Rather than pairs of contigs, each erroneous scaffold joined 10 contigs and we generated 200 such erroneous scaffolds. The remaining contigs were correctly scaffolded (ten contigs per scaffold) based on their location in the reference. The average sensitivity was 67.66% and specificity was 100% (S2 Table)(no correct scaffolds were broken). Most of the unflagged joins occurred near the ends of scaffolds and could be captured by decreasing the window size. Similar to assembly correction, we observed that sensitivity was highest with larger input contigs. Most of the misjoins missed by the algorithm were near the ends of scaffolds. The issue in detecting mis-assemblies in these regions is the low Hi-C physical coverage. Also, the other missed joins were between the small contigs which are hard to capture with Hi-C data alone. This evaluation highlights the accuracy of the mis-join detection algorithm to avoid over-scaffolding and provide a suitable stopping condition.
Scaffold accuracy. We evaluated scaffolds across three categories of error: orientation, order, and chimera. An orientation error occurs whenever the orientation of a contig in a scaffold differs from that of the scaffold in the reference. An ordering error occurs when a set of three contigs adjacent in a scaffold have non-monotonic coordinates in the reference. A chimera error occurs when any pair of contigs adjacent in a scaffold align to different chromosomes in the reference. We broke the assembly at these errors and computed corrected scaffold lengths and NGA50 (analogous to the corrected NG50 metric defined by Salzberg et al. [38]). This statistic corrects for large but erroneous scaffolds which have an artificially high NG50. We did not include SALSA1 in the comparison because for small contig sizes (200 kbp to 500 kbp), none of the scaffolds contained more than 2 contigs. For larger sizes (600 kbp to 900 kbp), the contiguity widely varied depending upon the minimum confidence parameter for accepting links between contigs.
Hi-C scaffolding errors, particularly orientation errors, increased with decreasing assembly contiguity. We evaluated scaffolding methods across a variety of simulated contig sizes. Fig 3 shows the comparison of these errors for 3D-DNA, SALSA2 without the assembly graph, and SALSA2 with the graph. SALSA2 produced fewer errors than 3D-DNA across all error types and input sizes. The number of correctly oriented contigs increased significantly when assembly graph information was integrated with the scaffolding, particularly for lower input contig sizes (Fig 3). For example, at 400 kbp, the orientation errors with the graph were comparable to the orientation errors of the graph-less approach at 900 kbp. The NGA50 for SALSA2 also increased when assembly graph information was included (Fig 4). This highlights the power of the assembly graph to improve scaffolding and correct errors, especially on lower contiguity assemblies. This also indicates that generating a conservative assembly, rather than maximizing contiguity, can be preferable for input to Hi-C scaffolding. All the assemblies described in this paper are available online and can be found at https://obj.umiacs.umd.edu/paper_ assemblies/index.html. Table 2 lists the metrics for NA12878 scaffolds. We include an idealized scenario, using only reference-filtered Hi-C edges for comparison. As expected, the scaffolds generated using only true links had the highest NGA50 value and longest error-free scaffold block. SALSA2 scaffolds were generally more accurate and contiguous than the scaffolds generated by SALSA1 and 3D-DNA, even without use of the assembly graph. The addition of the graph further improved the NGA50 and longest error-free scaffold length.

Evaluation on NA12878
We also evaluated the assemblies using Feature Response Curves (FRC) based on scaffolding errors [40]. An assembly can have a high raw error count but still be of high quality if the Assembly graph based scaffolding errors are restricted to only short scaffolds. FRC captures this by showing how quickly error is accumulated, starting from the largest scaffolds. Fig 5(D) shows the FRC for different assemblies, where the X-axis denotes the cumulative % of assembly errors and the Y-axis denotes the cumulative assembly size. The assemblies with more area under the curve accumulate fewer errors in larger scaffolds and hence are more accurate. SALSA2 scaffolds with and without the graph have similar areas under the curve and closely match the curve of the assembly using only true links. The 3D-DNA scaffolds have the lowest area under the curve, implying that most errors in the assembly occur in the long scaffolds. This is confirmed by the lower NGA50 value for the 3D-DNA assembly ( Table 2).
Apart from the correctness, SALSA2 scaffolds were highly contiguous and reached an NG50 of 112.8 Mbp (cf. GRCh38 NG50 of 145 Mbp). Fig 6 shows the alignment ideogram for the input contigs as well as the SALSA2 assembly. Every color change indicates an alignment break, either due to error or due to the end of a sequence. The input contigs are fragmented with multiple contigs aligning to the same chromosome, while the SALSA2 scaffolds are highly contiguous and span entire chromosomes in many cases. Fig 7(A) shows the contiguity plot with corrected NG stats. As expected, the assembly generated with only true links has the highest values for all NGA stats. The curve for SALSA2 assemblies with and without the assembly graph closely matches this curve, implying that the scaffolds generated with SALSA2 are approaching the optimal assembly of this Arima-HiC data.
We also evaluated the ability of scaffolding short-read assemblies for both 3D-DNA and SALSA2. We did not include SALSA1 in this comparison because it is not designed to scaffold short-read assemblies. We observed that use of the assembly graph when scaffolding significantly reduced the number of orientation errors for SALSA2, increasing the scaffold NGA50 and largest chunk almost two-fold. When compared to 3D-DNA without input assembly correction, SALSA2 with the assembly graph generates scaffolds of much higher NGA50 (7.99 Mbp vs. 1.00 Mbp). The number of scaffolding errors across all the categories was much lower in SALSA2 compared to 3D-DNA. We computed the CPU runtime and memory usage for both the methods while scaffolding long and short read assemblies with Arima-HiC data. We excluded the time required to map reads in both cases. While scaffolding the long-read assembly SALSA2 was 30-fold faster and required 3-fold less memory than 3D-DNA (11.44 CPU hours and 21.43 Gb peak memory versus 3D-DNA with assembly correction at 318 CPU hours and 64.66 Gb peak memory). For the short-read assembly, the difference in runtime was even more pronounced. SALSA2 required 36.8 CPU hours and 61.8 Gb peak memory compared to 2980 CPU hours and 48.04 Gb peak memory needed by 3D-DNA without assembly correction. When run with assembly correction, 3D-DNA ran over 14 days on a 16-core machine without completing so we could not report any results. Table 2. Scaffold and correctness statistics for NA12878 assemblies scaffolded with different Hi-C libraries. "True links" is an idealized case where the Hi-C links have been filtered in advance. The NG50 of human reference GRCh38 is 145 Mbp. The ratio between NG50 and NGA50 represents how many erroneous joins affect large scaffolds in the assembly. The bigger the difference between these values, the more aggressive the scaffolding was at the expense of accuracy. Longest chunk represents the longest error-free portion of the scaffolds. We observed that the 3D-DNA mis-assembly detection was overly aggressive in some cases, and so we ran some assemblies both with and without this feature. For the Illumina assembly as an input, 3D-DNA w correction did not finish within two weeks and is omitted. An evaluation of a previously published [20] 3D-DNA assembly from short-read contigs is included in S3 Table but did not exceed SALSA2's NGA50.  Assembly graph based scaffolding

Robustness to input library
We next tested scaffolding using two libraries with different Hi-C contact patterns. The first, from [42], is sequenced during mitosis. This removes the topological domains and generates fewer off-diagonal interactions. The other library was from [43], are in vitro chromatin sequencing library (Chicago) generated by Dovetail Genomics (L1). It also removes off-diagonal matches but has shorter-range interactions, limited by the size of the input molecules. As seen from the contact map in Fig 8, both the mitotic Hi-C and Chicago libraries follow different interaction distributions than the standard Hi-C (Arima-HiC in this case). Table 1 shows the mapping statistics for these libraries. We ran SALSA2 with defaults and 3D-DNA with both the assembly correction turned on and off.  The X-axis denotes the NGAX statistic and the Y-axis denotes the corrected block length to reach the NGAX value. SALSA2 results were generated using the assembly graph, unless otherwise noted. https://doi.org/10.1371/journal.pcbi.1007273.g007 For mitotic Hi-C data, we observed that the 3D-DNA mis-assembly correction algorithm sheared the input assembly into small pieces, which resulted in more than 25,000 errors and more than half of the contigs incorrectly oriented or ordered. Without mis-assembly correction, the 3D-DNA assembly has a higher number of orientation (250 vs. 81) and ordering (215 vs. 54) errors compared to SALSA2. The feature response curve for the 3D-DNA assembly with breaking is almost a diagonal (Fig 5(B)) because the sheared contigs appeared to be randomly joined. SALSA2 scaffolds contain longer stretches of correct scaffolds compared to 3D-DNA with and without mis-assembly correction (Fig 7(B)). SALSA1 scaffolds had a similar error count to SALSA2 but were less contiguous.
For the Chicago libraries, 3D-DNA without correction had the best NGA50 and largest correct chunk. However, the scaffolds had more chimeric join errors than SALSA2. SALSA2 outperformed 3D-DNA in terms of NG50, NGA50, and longest chunk when 3D-DNA was run with assembly correction. 3D-DNA uses signatures of chromosome ends [20] to identify break positions which are not prominently present in Chicago data. As a result, it generated more chimeric joins compared to SALSA2. However, the number of order and orientation errors was similar across the methods. Since Chicago libraries do not provide chromosome-spanning contact information for scaffolding, the NG50 value for SALSA2 is 5.8 Mbp, comparable to the equivalent coverage assembly (50% L1+L2) in [43] but much smaller than Hi-C libraries. Interestingly, SALSA1 was able to generate scaffolds of similar contiguity to SALSA2, which can be attributed to the lack of long range contact information in the library. SALSA2 is robust to changing contact distributions. In the case of Chicago data it produced a less contiguous assembly due to the shorter interaction distance. However, it avoids introducing false chromosome joins, unlike 3D-DNA, which appears tuned for a specific contact model.

Scaffolding non-model organisms
To evaluate the effectiveness of SALSA2 on a non-model organism, we used Hi-C data from recently published Anopheles funestus genome assembly which was scaffolded using an independent method (Phase Genomics or LACHESIS [16]) and manually curated using Illumina mate-pair support as well as FISH information [44]. This genome had high heterozygosity as the data was sequenced from a colony of mosquitoes rather than a single individual. Due to this, the assembly had a high duplication rate and was almost double the expected genome size. We scaffolded both the full assembly and the assembly after running purge haplotigs [45] using SALSA2 and 3D-DNA. For the post purge assembly, 3D-DNA generated an assembly with higher continuity but with more errors and a similar NA50 to SALSA2. (S3 Table). However, neither method performed well for the full assembly. SALSA2 was more contiguous than 3D-DNA (S4 Table) but was still very fragmented and much larger than the expected genome size. We conclude that heterozygous genome scaffolding remains a challenge and assemblies must either be de-duplicated beforehand or improved algorithms for scaffolding, such as [46] are needed.

Discussion
In this work, we present the first Hi-C scaffolding method that integrates an assembly graph to produce high-accuracy, chromosome-scale assemblies. Our experiments on both simulated and real sequencing data for the human genome demonstrate the benefits of using an assembly graph to guide scaffolding. We also show that SALSA2 outperforms alternative Hi-C scaffolding tools on assemblies of varied contiguity, using multiple Hi-C library preparations.
SALSA2's misassembly correction and scaffold misjoin validation can be improved in several ways. The current implementation does not detect a misjoin between two small contigs with high accuracy, mainly because Hi-C data does not have enough resolution to correct such errors. Also, we do not account for any GC bias correction when using the Hi-C coverage for detecting misjoins. Addressing these challenges in misjoin detection and misassembly correction is the immediate next step to improve the SALSA2 software.
The human genome is relatively homozygous compared to many other species. Assembly of many species is further complicated by DNA input requirements which necessitates pooling multiple individuals. SALSA2 does not remove duplication present in an input assembly and thus requires pre-processing by another tool, such as Purge Haplotigs [45] or haplomerger [47]. Once contigs are classified into the "primary" and "haplotig" sets, SALSA2 could be run on each of the sets independently.
Hi-C scaffolding has been historically prone to inversion errors when the input assembly is highly fragmented. The integration of the assembly graph with the scaffolding process can overcome this limitation. Orientation errors introduced in the assembly and scaffolding process can lead to incorrect identification of structural variations. On simulated data, more than 50% of errors were due to inversions, and integrating the assembly graph reduced these by as much as 3 to 4 fold. We did not observe as much improvement with the NA12878 test dataset because the contig NG50 was much higher than in the simulation. However, it is not always possible to assemble multi-megabase contigs. In such cases, the assembly graph is useful for limiting Hi-C errors.
Most existing Hi-C scaffolding methods also require an estimate for the number of chromosomes for a genome. This is implicitly taken to be the desired number of scaffolds to output. As demonstrated by the Chicago, mitotic, and replicate [48] Hi-C libraries, the library as well as the genome influences the maximum correct scaffold size. It can be impractical to sweep over hundreds of chromosome values to select a "best" assembly. Since SALSA2's mis-join correction algorithm stops scaffolding after the useful linking information in a dataset is exhausted, no chromosome count is needed as input.
Obtaining the chromosome-scale picture of the genome is important and there is a tradeoff between accuracy and continuity of the assembly. However, we believe that manual curation to remove assembly errors is an expensive and involved process that can often outpace the cost of the rest of the project. Most of the assembly projects using Hi-C data have had a significant curation effort to date [19,49]. Thus, we believe that not introducing errors in the first place is a better strategy to avoid the later burden of manual curation of small errors in chromosomes. The Hi-C data can be used with other independent technologies, such as optical mapping or linked-reads to reach accurate chromosome-scale scaffolds. 3D-DNA was recently updated to not require the chromosome count as input but the algorithm used has not been described. Interestingly, it no longer generates single-chromosome scaffolds in our experiments, a major claim in [20], supporting a conservative scaffolding approach. Even while scaffolding short-read assemblies, we observed that SALSA2 generated more accurate scaffolds than 3D-DNA, indicating the utility of SALSA2 in scaffolding existing short-read assemblies of different genomes with the newly generated Hi-C data.
As the Genome10K consortium [50] and independent scientists begin to sequence novel lineages in the tree of life, it may be impractical to generate physical or genetics maps for every organism. Thus, Hi-C sequencing combined with SALSA2 presents an economical alternative for the reconstruction of chromosome-scale assemblies.