Improving mammalian genome scaffolding using large insert mate-pair next-generation sequencing

Large insert mate-pair library generation

To assess library complexity by determining the maximum number of unique reads obtainable from a MP library, two of the large-insert libraries (20-kb and 25-kb) were sequenced to a higher depth in an additional sequencing run. To assess reproducibility of the adapted MP protocol, three libraries (5-kb, 8-kb, and 15-kb) were generated in duplicate from independently isolated, sheared, and separated genomic DNA samples. Insert size distributions of the individually produced replicates were highly consistent (Figure 1a and Table 1). In total, 192.4 million pairs of MP reads and 160 million pairs of PE reads were generated. A total of 62.3 million non-duplicate MP read pairs and 131 million non-duplicate PE reads were consistently mapped against the rat reference genome, resulting in a genome-wide physical coverage of 228.5× (220× for MP libraries and 8.5× for 170-bp PE). Less than 1% of the paired reads were inverted (one of the reads in other orientation than expected) or everted (both reads in other orientation resulting in wrong order of tags) and approximately 10% were mapped remotely (i.e., to a distant genomic position, significantly deviating from what is expected based on the insert-size distribution). Remote, inverted, or everted events represent a mixture of 1) library construction artefacts due to chimeric molecules, 2) errors in the reference genome assembly (misassemblies) and 3) real structural differences between the reference strain and the substrain tested here. The first category typically involves stochastic events that are supported by a single read pair and that are filtered out by requiring multiple independent supporting read pairs for calling.

Library sequencing and quality assessment

Sequencing libraries may suffer from low complexity due to library amplification steps in the protocol. When the proportion of unique library molecules is low due to inefficient molecular reactions or low amounts of input material, sequencing more reads of that same library would not yield any additional information, but only extra copies of previously sequenced molecules (duplicate reads). We assessed the complexity of each library by plotting the number of read-pairs with unique genome coordinates against the total number of all mapped pairs (Figure 1b). In general, the complexity of the small-insert libraries is higher than that of the large-insert libraries, which more quickly saturate to the level where deeper sequencing delivers predominantly non-informative duplicate reads. Duplicate reads do not necessarily affect the utility of the libraries, because these reads are filtered out as a first step in the analysis procedure; however, low complexity does decrease the capacity to obtain sufficient physical genome coverage. Three sample groups can be distinguished in Figure 1b: (1) high-complexity libraries that deliver approximately 100 million unique pairs (PE, 3kb, 5kb_b and 8kb_b), (2) medium-complexity libraries that result in about 10 million unique pairs (5kb_a, 8kb_a and 20kb), and (3) low-complexity libraries resulting in approximately 1 million unique pairs (15kb_a, 15kb_b, 25kb). Several of the low-complexity libraries show a plateau in the curve, indicating that these have been sequenced to saturation (25kb, 15kb_a). For others (5kb_b, 8kb_b), deeper sequencing would be informative.

Library complexity may be influenced by several experimental conditions. When starting with an equal quantity of genomic DNA, fragmentation for a standard PE library provides approximately 140-fold more unique molecules than for a 25-kb library. Furthermore, MP library preparation involves a circularization step (Additional file 2) that becomes less efficient as the size of the molecule increases. Quantification of DNA before and after circularization (and removal of non-circularized molecules) showed a circularization efficiency of up to 37% for libraries below 10 kb and 5–10% for libraries above 10 kb (Additional file 3). Each of these library generation steps has a negative impact on the recovery of material; for example, an input of 10 μg 25-kb size-selected DNA would result in approximately 6 ng (>4,000-fold reduction) of DNA for adapter ligation and subsequent adapter-mediated PCR. As a consequence, more PCR cycles are required for larger insert libraries to obtain sufficient amounts of library DNA for NGS (Table 1). Although the 3- and 5-kb insert libraries could routinely be generated at high complexity, we observed more technical variation for the large insert libraries. For example, the 20-kb library required only 14 PCR cycles during the library preparation procedure and performed well in the complexity analysis (comparable to 5- and 8-kb libraries). The 15- and 25-kb libraries required 21 and 17 cycles, respectively, and resulted in libraries of lower complexity (Figure 1b). These results indicate that the number of required PCR cycles is a very good predictive parameter for library complexity.

The 5-, 8-, and 15-kb libraries were generated in duplicate using DNA isolated from two different tissues of the same animal. The insert size distribution was found to be highly reproducible (Figure 1a), but the library complexity was much more variable between duplicates (Figure 1b). These differences might have been due to differences in DNA quality (e.g. amount of single strand breaks) or purity (e.g. associated protein or small molecule contaminants) of the DNA and subsequent differences in shearing efficiency. Indeed, DNA yields after size fragmentation were as much as 2.5-fold lower for the duplicate DNA sample (data not shown), which systematically resulted in less complex libraries. Most importantly, however, statistics for the amount of consistently mapped read pairs were comparable for all replicates (Table 1), indicating that the mapped unique read pairs were similar in quality (e.g., low chimaerism) and insert size. Low complexity in libraries could be circumvented by using larger amounts of input DNA and/or by optimization of shearing conditions to concentrate DNA in the desired size range. In our experiments we aimed for a broad size distribution to be able to simultaneously extract DNA for a range of different sizes.

Although the larger insert libraries come with more duplicate reads, far fewer sequencing pairs are required to physically cover the complete genome. It should be noted, that for all MP libraries in the experiments described here more than 10x physical coverage was obtained, including 48x coverage for 20 kb inserts.

To assess the value of the various insert size libraries for genome structure analysis, we determined the ability of each library to (1) physically cover the reference genome and overlap various repeat elements, (2) drive contig scaffolding, and (3) fix contig assembly issues in the current genome assembly (errors in contig order and orientation).

Spanning repeats and physical genome coverage

The ability to physically cover a complete genome by sequencing is not only determined by the length of the read, the insert size of the library, and the number of paired reads, but also depends on genome-specific characteristics, like the composition and distribution of repetitive elements. The rat genome is representative for other mammalian genomes and contains 1.24 Gb of repetitive sequences, which is over 49% of the 2.51 Gb in the current reference genome assembly (RGSC 3.4, v.66 [24]). Retrotransposable LINE (L1) elements are the largest class of repeats with a total length of 474.6 Mb (18.9% of the genome), followed by retrotransposons that are flanked by long terminal repeats (LTRs; 220.9 Mb; Figure 2a and b). To evaluate the effect of library insert size on the degree of physical genome coverage, we merged data from duplicate libraries with the same insert size. Although the MP libraries had far more physical genome coverage, all datasets were normalized to an equal physical genome coverage based on properly mapped and oriented read pairs, which was limited to 8.5x by the available amount of data for the PE library. Next, we determined per library the fraction of bases per contig of the rat reference genome that is physically covered, specifically focusing on repetitive elements. Despite the same physical genome coverage and much higher base coverage, short-insert libraries (PE and 3 kb MP) were much less efficient in spanning long repetitive elements, such as LINEs or LTRs, than larger insert MP libraries (≥ 15 kb) (Figure 2c and d respectively). As expected, PE pairs overlapped hardly any of these elements but also the most widely used 3-kb MP libraries were found to only span approximately half of the 3-kb repeat elements, and only a few elements with sizes above 4 kb. Slightly improved results were observed for the 5-kb and 8-kb MP libraries, where approximately half of the repeats with a matched size could be spanned by at least one mate-pair. The 15-, 20-, and 25-kb libraries spanned over 90% of the repeat elements across the whole size spectrum and all displayed a very similar performance, indicating that there is limited added value for even larger insert sizes.

Contig scaffolding

To evaluate the utility of the various libraries for guiding genome assembly, we simulated the scaffolding step of such process by using the 137,257 contigs from the current rat genome build and the different MP data sets as input for the SSPACE 2.0 [25] software. To allow for library insert-size comparison, we again used the normalized datasets at 8.5x physical coverage and determined the N50 (at least half of the genome bases are in scaffolds that are equal or exceeding the N50 value) and the number of scaffolds (segments of the genome reference consisting of contigs in known order, separated by gaps) for each individual library and combinations thereof, using the output of the SSPACE software (Figure 3a, Table 2 and Additional file 4). When we consider only the utility of single libraries, the N50 increases from ~38 kb for the PE data to 140–163 kb for the MP libraries of 5 kb and up. PE libraries are not effective in reducing the number of scaffolds as compared to the capillary sequencing-based contigs: a reduction of only 15 scaffolds is obtained (from 137,256 scaffolds in RGSC3.4 to 137,241 scaffolds using the PE data). In contrast, individual MP libraries decreased the number of scaffolds by up to more than 50% (~67,000 for the 5 kb library, which performs best of all individual MP libraries). When considering two insert size libraries, combination of 5 or 8 kb and 20 or 25 kb are most optimal with N50's of ~0.5 Mb. Intriguingly, 3 kb mate-pair libraries, which are most commonly used, showed the worst performance from all MP libraries, also when combined with other libraries. Including all libraries in the scaffolding process results in a further decrease of scaffolds (36,348) with an N50 increase up to 1.3 Mb. Increasing the physical coverage for a single insert library shows to be far less effective than combining libraries with different insert size (Additional file 5 and Additional file 6). For example, when we increase the coverage of the 5 kb insert library to 34x physical coverage, the N50 increases from 141 kb to 262 kb. However, combining 8.5x coverage data for the 5 kb insert library with similar coverage of any other MP library (up to 34x) results in much higher N50 values ranging from ~431 kb up to ~1 Mb.

Reference genome improvement

Finally, we evaluated the value of various MP insert sizes for improving the existing rat genome assembly. To this end, we compared de novo scaffolds constructed using MP data with independently obtained genome-wide optical mapping data. Optical mapping is an integrated system that provides long-range genome structural information by the construction and analysis of genome-wide, ordered restriction maps [9, 10, 26, 27]. We limited the analysis of concordance between sequence scaffolds and optical maps to one of the small rat chromosomes (RNO18), because the fine level optical structural alterations (OSAs) that were automatically called by the optical mapping pipeline [10] were manually curated between sequence scaffolds and optical maps, which required exploration on a case-by-case basis for mediation at the nucleotide-level. We divided the chromosome into 872 100-kb windows and found that 96 out of 872 of such bins harboured structural changes within scaffolds of RNO18 when comparing the MP-updated genome structure with the original genome. The 96 bins contained a total of 199 unique inconsistent connections between contigs within scaffolds. Next, we looked at structural differences between scaffolds of RNO18, based on the comparison of the MP-based scaffolding and the reference genome and observed many more bins to be affected (166/ 872 bins containing a total of 1374 inconsistent links between scaffolds). In total, 236 bins showed one or both types of inconsistent connections. Of these 236 bins, only 106 showed concordance with the reference genome. 130 bins were found to contain OSAs including absence or discordance of alignment between the optical maps and RGSC 3.4 genome assembly (detailed description in Additional file 7 and Additional file 8). Because the optical mapping system constructs ordered restriction maps and does not evaluate genome structure at the nucleotide level, not all discordances detected by the mate-pair analysis are revealed through optical mapping data. For example, small contigs or changes that do not overlap with a SwaI restriction site will not be identified.

We explored two of the largest segments with long-range disagreement between optical maps and RGSC3.4 assembly and conclude that MP-assisted re-scaffolding can recover concordance with the independently generated optical maps (Figure 3b). The complete MP data described here has therefore also been used for building the new genome reference of the rat (Rnor5.0, GenBank ID GCA_000001895.3, unpublished results).

Generation of MP libraries and mapping

To allow for the construction of large insert MP libraries, we modified the standard SOLiD 4 mate-pair library preparation protocol (Additional file 2). In short, genomic DNA was isolated from Brown Norway (BN/RijHsd) rat brain and testis tissue. DNA (100 μg) was sheared under mild conditions using HydroShear (JHSH204007, 20 cycles, SC15) and subsequently end-repaired (Epicentre End-It^TM DNA-end repair kit) in 1 mL End-It mix per 100 μg input DNA. CAP adapters were ligated in 500-μL reaction volumes of New England Biolabs Quick Ligase reaction mix. The amount of ligated adapter was determined based on the DNA content (100 pmol CAP adapter/pmol DNA). Following CAP-adapter ligation, DNA fragments were purified with phenol:chloroform:isoamylalcohol (PCI, pH7.9) by gentle mixing and centrifugation in MaXtract high-density tubes (QIAgen, 1.5mL, #129046). The fragmented DNA was separated via pulsed-field gel electrophoresis (PFGE; Bio-Rad CHEF Mapper XA system). PFGE conditions and settings were: 1% low melt agarose gel (Invitrogen, #16520100), 0.5× TBE, 14°C, 19 hours, forward current: 9.0 V/cm, switch time 0.08 s–0.46 s, reverse current 6.0 V/cm, switch time 0.08 s–0.46 s. Multiple size ranges (<7 kb, 7–10 kb, 10–14 kb, 14–18 kb, 18–24 kb, 24–33 kb, and >33 kb) were selected from the gel using a 1 kb extension ladder (Invitrogen, #10511-012). For unknown reasons, actual library insert sizes after library construction and mapping tend to be lower than their initial appearance on gel. A probable explanation for this is that most library construction steps may show a small bias towards smaller molecules (e.g. during circularization). Further on in the manuscript, each library will be referred to as the actual insert size as determined after data analysis. Following gel excision, DNA was carefully recovered using GELase^TM (Epicentre, #G09200). DNA fragments were circularized with a biotinylated internal adapter at a final DNA concentration of 1 nanogram per microliter (ng/μl). For every 40-μl reaction volume, 1 μl T4 DNA ligase was used. The reaction mixture was purified, and non-circularized fragments were removed by a plasmid-safe DNase treatment (Epicentre, #E3101K). Following linear DNA removal, DNA polymerase I-directed nick translation “pushed” the nick from the adapters into the circularized target DNA to generate sufficient tag length for sequencing (~100 bp for each tag, 13 minutes on ice water [0°C], inactivated with PCI). T7 exonuclease and S1 nuclease treatment were used to digest the circles at the position of the nick. The digested fragments of approximately 300 bp in size were end-repaired and bound to MyOne C1 streptavidin beads via the biotinylated internal adapter. Standard SOLiD P1 and P2 adapters were ligated to the blunt ends of the library molecules (a-tailing and alternative adapters should be used at this step to make the library compatible with the other sequencing platforms like Illumina, see Additional file 10), followed by another round of nick translation to remove the nick introduced by adapter ligation. Mate-pair libraries were amplified by PCR for 13–21 cycles, depending on the library. Amplification of the 14–18-kb size range samples (with an estimated final insert size of 10–12 kb) did not result in sufficient material (most likely because of unsuccessful adaptor ligation) and were not further included in the process. For all other insert size libraries we continued with templated bead preparation and libraries were successively sequenced on the SOLiD 4 system. For all libraries together, 192.4 million pairs of MP reads (AB/SOLiD V4, two slides) were sequenced. Paired reads were mapped against rat reference genome RGSC 3.4 using BWA v0.5.9 [37], and non-unique (based on identical read start sites for the forward and reverse read) and ambiguously mapped read pairs were removed from the data set.

Generation of the paired-end library (170-bp) and mapping

For construction of the paired-end library (PE; 170-bp insert), the SOLiD 3 protocol for fragment library preparation was used (SOLiD™ 3 System Library Preparation Guide; Section 2.1). DNA (3 μg) derived from the Brown Norway rat was used for shearing using Covaris S2 (10 cycles of 60 s, intensity 5, 100 cycles/burst, 4°C). Sheared DNA was end-repaired using the Epicentre End-It^TM DNA-end repair kit. P1 and P2 adapters were ligated to the DNA fragments, and the library molecules were selected based on size (220–300 bp, including 90 bp for both adapters). PCR amplification was done for 5 cycles with primers specific to the adapters to obtain sufficient library molecules for ePCR and sequencing. Sequencing was done in paired-end mode (forward and reverse tag; 50 bp and 35 bp, respectively) on the SOLiD 3 system (1 slide AB/SOLiD V3). Approximately 1.6 × 10⁸ paired-end reads were sequenced (95% non-clonal) and mapped with BWA, resulting in a data set with a median insert size of 170 bp.

Calculation of library insert size and contig scaffolding

Forward and reverse reads were mapped independently against contigs of rat RGSC3.4 genome assembly using BWA 0.5.9. Only read pairs with a single best hit for each tag (X0 flag equal to 1) were taken into consideration for estimate of insert size distribution. Analysis of library complexity was done by randomly sampling of reads from a library and determining the number of non-clonal pairs. Next, read pairs with exactly the same mapping coordinates of forward and reverse tags were marked as clonal and excluded from further analysis. Distribution of insert sizes was estimated from read pairs with proper orientation and distance between tags (below 100 kb).

To allow comparison of different library insert sizes for contig scaffolding, we created a subset of data for each library that corresponds to 8.5x non-clonal physical coverage of rat genome. We randomly sampled read pairs from each library, computing physical coverage represented by non-clonal pairs with expected orientation and distance between tags on chromosomal level (skipping pairs corresponding to first and last percentiles of insert size distribution). Read pairs from the normalized datasets with forward and reverse tags mapped to different contigs were selected for scaffolding analysis. Scaffolding was performed using SSPACE v2.0 software [25] with default parameters. The order in which the libraries were used by SSPACE was as recommended by the SSPACE manual - always from the smallest to largest insert size.

Repeat analysis

Repeat annotation of the rat genome reference was obtained from Ensembl database [24] (v.66) and was used for calculation of size distribution and abundance of different repeat types. This annotation was used to determine the percentage of repeat elements spanned by mapped fragments from each mate-paired library. We used the normalized dataset where every library had 8.5x physical coverage (as described above). We considered only unambiguously mapped read pairs that had a proper orientation of tags and did not exceed 99^th percentile of fragment size distribution. Since individual copies of a mobile element or repeat class differ in size, we used 500 bp windows for calculation of percentage of mobile elements overlapped by fragments from each library.

Comparison to optical maps

Original RGSC3.4 scaffolds and those obtained after NGS-assisted re-scaffolding were compared to each other by nucleotide BLAST search. Nearly identical (>99%) super-kb segments were plotted as Harr-plot visualization graphs. Most evident discordant regions were manually selected and the corresponding genomic segments were compared to optical maps [10] of Brown Norway rats. Optical maps are generated from large, randomly sheared high-molecular weight genomic DNA molecules that are stretched on a microscope slide. After stretching, the DNA molecule is digested with a SwaI restriction enzyme. While the DNA remains attached to the slide, the cuts become visible under the microscope as small gaps and the sizes of the stained DNA fragments can be measured. Multiple optical maps are then combined to form a comprehensive reference optical map for the Brown Norway rat genome. The optical maps used in this study are produced in the lab of David C. Schwartz and are available upon request. To compare our assembly with existing Brown Norway rat optical maps, nucleotide sequences of MP-enhanced scaffolds were digested in silico with a SwaI restriction enzyme. The in silico digested fragments were plotted next to optical maps, originally aligned to the RGSC 3.4 assembly and visually inspected for concordance.

Availability of data

The data sets supporting the results of this article are available in the ArrayExpress repository under accession number E-MTAB-1082.