Enhanced whole exome sequencing by higher DNA insert lengths

Background Whole exome sequencing (WES) has been proven to serve as a valuable basis for various applications such as variant calling and copy number variation (CNV) analyses. For those analyses the read coverage should be optimally balanced throughout protein coding regions at sufficient read depth. Unfortunately, WES is known for its uneven coverage within coding regions due to GC-rich regions or off-target enrichment. Results In order to examine the irregularities of WES within genes, we applied Agilent SureSelectXT exome capture on human samples and sequenced these via Illumina in 2 × 101 paired-end mode. As we suspected the sequenced insert length to be crucial in the uneven coverage of exome captured samples, we sheared 12 genomic DNA samples to two different DNA insert size lengths, namely 130 and 170 bp. Interestingly, although mean coverages of target regions were clearly higher in samples of 130 bp insert length, the level of evenness was more pronounced in 170 bp samples. Moreover, merging overlapping paired-end reads revealed a positive effect on evenness indicating overlapping reads as another reason for the unevenness. In addition, mutation analysis on a subset of the samples was performed. In these isogenic subclones, the false negative rate in the 130 bp samples was almost double to that in the 170 bp samples. Visual inspection of the discarded mutation sites exposed low coverages at the sites flanked by high amplitudes of coverage depth. Conclusions Producing longer insert reads could be a good strategy to achieve better uniform read coverage in coding regions and hereby enhancing the effective sequencing yield to provide an improved basis for further variant calling and CNV analyses.


Background
During the past years whole exome sequencing (WES) has gained much popularity in research and diagnostics, as focussing on protein-coding regions reduces sequencing costs compared to whole genome sequencing (WGS) [1][2][3][4]. Concentrating on exonic regions minimises the sequencing target area of the human genome with about 3 Gbp to less than 2 % [3,4]. There is a broad application area of WES such as variant calling [3][4][5] and analysis of copy number variations (CNV) [6,7] demonstrating its general usefulness in the genomic field. In terms of quality validation, many research efforts were focussing on the comparison of different exome enrichment platforms [8][9][10][11] and their performance to WGS techniques [12,13]. Apparently, one recurrent bias appearing with WES is its inhomogeneous coverage across targeted protein-coding regions, which is suggested to be resolved by increasing read depth. However, increasing depth for WES would also place the economic costs for the alternative WGS in similar range to WES [12,13]. For example, in one study the exome coverage of ≥ 20× needed an average of 160× with WES, whereas WGS was sufficient with 44-56× [12].
In the midst of this debate, we would like to propose a new aspect to the technical side of WES. During standard WES (2 × 101 paired-end sequencing) of our cancer cell lines, we observed a strong irregular distribution of read coverage along exons, which had a size of ≥1 kbp. In the following, the genomic DNA fraction sequenced as paired-end reads and flanked by illumina adapters is denominated as insert, in order to confine these to sequencing library fragments including illumina adapters. Calculation of the original DNA insert size revealed that the genomic DNA insert had a peak size of 130 bp. This prompted us to consider, whether longer purified DNA inserts might improve the evenness of read coverage. In the past, DNA insert lengths of 200-250 bp for 2 × 90 paired-end reads were mentioned to contribute to library optimisation, nonetheless, results were not shown [8]. Additionally, different fragmentation techniques improved coverage, yet DNA insert length were indicated for sonication (161 bp) but not for enzymatic fragmentation [10]. In another study short inserts of 100-200 bp are suggested instead of 500-800 bp [1], however, at that time the development of exome capture design just started and more importantly, very short single reads of 26 bp were analysed. Similarly, short reads of 35 and 50 bp were sequenced for a further publication, in which short insert sizes of 120 bp are recommended given the short median length of 120 bp of human exons [14]. Apparently, no specific study on the effect of DNA insert sizes to paired-end sequencing has been published to date. Therefore, we tested in this study, whether the genomic DNA insert length influences the uniformity of read coverage within targeted regions. Beside samples with 130 bp peak insert length, a sample group with 170 bp peak inserts was produced and validated. The evenness score [14] was applied as metric for assessing the effectiveness of target region coverage. Altogether, we would like to provide a short technical note on the effect of DNA insert length on the evenness of coverage for paired-end sequenced reads.

Production of WES libraries with two different insert sizes
For exploring the effect of different genomic DNA insert lengths on the uniformity of coverage in whole exome sequencing (WES), we varied the DNA shearing by acoustic fragmentation (Covaris). Two DNA inserts groups of six samples each were yielded with approximately 130 and 170 bp length (Fig. 1a+b). Exome enrichment was conducted with Agilent SureSelectXT All Exon v5+UTR/v5 and sequencing produced 18-54 million reads per sample (see Table 1). After trimming, removing PCR duplicates and mapping the reads to the human GENCODE genome v21, manual inspection of alignments to many exons exceeding 1 kbp indicated large amplitudes of coverage for 130 bp compared to 170 bp samples (Fig. 1c). This observation, that longer inserts might compensate for "mountain-valley" profiles in coverage, encouraged us to examine this in detail.
Since the mean DNA insert peak for each sample group with 130 and 170 bp, respectively, was smaller than the resulting paired-end sequenced reads of 2 × 101 bp, a high percentage of paired-end reads contained overlapping sequences. As these overlaps did not contain further information for e.g. variant calling and CNV analysis, pairedend reads were joined where overlapping sequences were found and aligned to the human genome. Trimming, read mapping and joining statistics are summarised in Table 2 and Phred quality scores for the sequencing cycles are visualised in Fig. 1d. Nearly all sequenced bases for further analyses had quality scores ≥30.
In the following sections unmerged and merged sequences are compared along with contrasting 130 versus 170 bp insert results.

Differences in exome capture and read depth
We applied two different exome enrichment kits in this study, namely Agilent SureSelectXT v5+UTR and v5 on the 130 bp genomic DNA insert and 170 bp group, respectively (see Table 1). The genomic target region both kits are covering were nearly identical except for the UTR stretches which were exclusively contained in v5+UTR: v5 target regions constitute 68 % of v5+UTR (75 Mb) and 99.9996 % of v5 (50 Mb) were included in v5+UTR (Fig. 2a). Hence, for further comparison the specific target region needed to be considered.
Although aiming to adapt the read numbers to the respective target regions during the sequencing process, the mean coverage to the target region sizes 75 and 50 Mbp was increased for 130 bp compared to 170 bp inserts (Fig. 2b). However, this difference attenuated when computing coverage means for merged paired-end sequences. Similar results were obtained when calculating the mean coverage of 130 and 170 bp inserts on their respective target regions and on overlapping target regions of v5 and v5+UTR (Fig. 2c). Intriguingly, concentrating on the fraction in the respective target region at ≥10× read depth revealed a smaller fraction of 130 bp samples covered at that mininum read depth than for 170 bp (Fig. 2d). Low coverages in turn mean impaired mutation detection in these regions, whereas at the same time excessive coverage of target regions in 130 bp samples (Fig. 1c) seems dispensible.

Comparing uniformity of coverage
As a measure for the skewed distribution of exome captured sequences in the coding regions, the evenness score was calculated allowing for target region and library size correction [14] (Fig. 3). Hereby, after normalising the coverage of each sample to the integral for each coverage curve to 1 (Fig. 3a+b), the area below the curve between 0 and 1 is defined as the evenness score [14]. Bearing in mind that the higher the evenness score, the a d b c Fig. 1 DNA shearing to 130 and 170 bp fractions before Illumina adapter ligation; sequencing base quality. a DNA insert length distribution per sample. b Peak insert lengths for the two different sample groups. c Alignment histograms for 130 bp insert samples (red) exhibited high amplitudes of coverage within the exon in comparison to 170 bp (blue) as exemplified by this gene BMP4 via IGV. Target regions of Agilent v5 and v5+UTR are given in the last two lines. Please note the 3× fold higher maximum coverage of 130 bp samples. d High Phred score quality values for mapped paired-end reads. Base calling quality was high after trimming and mapping to the human genome. As expected, for both reads in forward and reverse direction (1-100 and 101-200 bases) read quality increased during the first 10 cycles and dropped gradually due to de-phasing errors of Illumina's sequencing pipeline. After joining paired-end reads, quality scores improved between 75-125 cycles, as the best scores were kept while merging. Quality scores were ≥30 throughout nearly all cycles and similar between 130 and 170 bp samples more even the coverage, the evenness score for 170 bp insert samples was clearly higher than 130 bp and hence its coverage more even than the 130 bp group (Fig. 3c) despite the higher average coverage of 130 bp (Fig. 2c).
Excluding overlapping sequences within the paired-end reads yielded a substantial rise in evenness for both insert groups.
Since using different enrichment kits for 130 and 170 bp samples respectively, the additional UTR target region fraction might be the cause of the unevenness observed in 130 bp inserts. However, the differences in the evenness scores were comparable to the results above when calculated on the common target region of v5 for both 130 and 170 bp inserts (Fig. 3c). The augmented evenness score for 130 bp inserts in the common coding target regions compared to its v5+UTR target region might even hint on inferior read coverage in the UTRs or inversely enhanced coverage (Fig. 2c) and uniformity for the coding regions.
For WES applying paired-end sequencing, libraries with small genomic DNA insert length produce overlapping sequences. These overlapping bases within one pairedend read provide no extra information on an alternate DNA strand or another allele, since they stem from the identical original genomic DNA sequence. The more overlapping bases within a paired-end read, the more bases remain unused, hence the effectivity to gain coverage shrinks with low insert sizes. Moreover, joining pairedend reads yielded in higher evenness scores particularly for 130 bp inserts (Fig. 3c), thereby showing another negative impact of redundant overlapping sequences.
On the other side, the median size of human exons is 120 bp, thus many bases will map off-target with longer inserts [14], however, the coverage in longer exons would reach higher uniformity and higher minimum depth (Fig. 2) instead of distinct "mountain-valley" coverage silhouette. Increasing coverage depth unfortunately would not yield in proportionally homogeneous coverage ( Fig. 2 and [12]). Higher evenness in turn is prerequisite to effectively detect mutations, which is achieved with longer inserts (Fig. 3c). These longer inserts can be produced by a plain technical adjustment in the DNA sample fragmentation step. It may be speculated that peak library sizes of >200 bp will result in even better evenness and simultaneously minimise costly overlapping read bases and reduce excess coverages.

Missing mutations in shorter insert samples
To demonstrate that shorter DNA inserts do increase the false negative rate for mutation analysis, we included four isogenic subclones of the human cell line HG-3. All four HG-3 subclones were sequenced at a comparable depth to   (Fig. 4a), although coverage was slightly higher in 130 bp inserts for common target coding regions (Fig. 2c). Manual inspection of 44 selected mutations with a minimum depth ≥20× in both 130 bp peak insert samples and simultaneously <10× in 170 bp samples and vice versa delivered mutations in four genes, that were discarded by longer insert samples, but mutations in 20 genes were missed by shorter inserts (e.g. Fig. 4b). Many of the failed mutations by the shorter inserts were found in regions, where amplitudes between maximum and minimum depth were high in the target regions forming a "mountain-valley" pattern as for the mutation in OR5H15 (Fig. 4c). Although this mutation dropped out in the shorter insert samples due to low depth in a "valley", the biallelic nature seen in the longer insert was also true for shorter inserts. Here, higher coverage depth in target regions for 130 bp samples could not improve mutation detection in this gene. Expanding the insert length resulted in mapping more reads to off-target regions of OR5H15, which was more appropriate for this bait design.
The question arose whether the bait-balancing differences between the two exome enrichment kits may account for the skewed coverage of 130 bp compared to 170 bp samples. Bait-balancing is applied for adjusting the number of probes to the binding efficiency of the targeted regions. For OR5H15 as an example all coverage histogram peaks reflected the target regions, which were identical to v5 and v5+UTR (Fig. 4c). An unoptimised bait-balance between the given oligos within this single exon for both insert groups was apparent and the bait-balance of the v5 and v5+UTR kit was similar if not identical. It seems, that the manufacturer's recommended 130-150 bp peak insert fragmentation was insufficient to span the entire annotated exon of OR5H15 on the basis of the target region rather than due to bait-balancing effects of the different capture kits used.
To overcome this technical shortcoming, the baits of the capture assay could be redesigned in closer proximity, or the potential of the current platform could be exploited by applying longer insert sizes. This would be a trade-off between too low as well as excessive read coverages and capturing off-target regions. Since flanking sequences to target regions often reside within annotated gene regions, achieving higher uniform coverage should be prioritised.

Conclusions
Although widely used, one major drawback of WES is its skewed coverage distribution within the targeted exome. By simply enlarging genomic DNA fragment sizes before exome capturing, the evenness of coverage can be augmented. We think that WES with an average coverage of 80x in contrast to WGS will remain feasible for smaller studies with limited budget in the next years, therefore any optimisation of this technology is assumed to affect a broad community.
Hence, increasing the DNA insert length maybe even longer than 170 bp will gain better uniform read coverage for WES and thus provide an improved basis for variant calling and CNV analyses at minimised sequencing costs. a c b Fig. 3 Evenness between different insert groups and unmerged/merged sequences. Before (a) and after (b) normalisation of coverage to the fraction of respective target regions for unmerged sequences. The complete integral of normalised coverage to the target region is summing up to 1. c The evenness score computed from the area under the curve of unmerged (Fig. 3b) and merged sequences between 0 and 1 normalised coverage. The closer the evenness score is to 1, the better the uniformity of base coverage. The impact of higher insert length was evident; merged inserts gained top evenness scores regardless of relating to the specific corresponding target region or to overlapping target regions

Samples, exome enrichment and sequencing
A selection of 2 × 6 human cancer cell lines (see Table 1) was prepared for WES all held by the DSMZ cell line bank (http://www.dsmz.de) and cultured as described previously [15].
Fragmentation of 100 ng purified genomic DNA (gDNA) in 55 μl Tris-EDTA buffer was done on Covaris S2 and the procedure adjusted to obtain fragments with a peak length of 130 and 170 bp, respectively. After library preparation from 100 ng of fragmented gDNA using Agilent SureSelectXT v5 (50 Mb) and v5+UTR (75 Mb), libraries were purified, size validated and prepared for sequencing according to the manufacturer's protocols. Libraries were sequenced on Illumina HiSeq2500 using TruSeq SBS Kit v3-HS (2 × 101, paired-end run). Concentration, quality, fragment sizes of purified genomic DNA (gDNA) and libraries were controlled by Agilent Technologies 2100 Bioanalyzer (Agilent Technologies; Waldbronn, Germany).

Sequence processing, mapping and data analysis
Before mapping raw reads in fastq format, sequences were trimmed at the ends for low quality (<Q20) or adapter contamination by fastq-mcf of the ea-utils toolbox (version 1.  Reads were sorted in IGV to bases at the mutation site, hence all detected Ts for 130 bp samples (3 and 1, respectively) are indicated. c Specific target regions of v5+UTR and v5 in gene OR5H15 were identical to which the coverage histogram peaks map. Target regions are given in the last two lines. The amplitudes were higher for the 130 bp samples as well as the maximum read depth in the visible region compared to 170 bp samples. At the same time pronounced amplitudes were also obvious for 170 bp within the gene region of OR5H15 implying that sequencing longer DNA fragments would have gained even smaller amplitudes and a higher coverage across the gene region. The four subclones carried four further mutations (grey arrows) beside the failed mutation of 130 bp samples (red arrow) indicating sequence similarity of the subclones Agilent SureSelect All Exon v5 (S04380110) and Agilent SureSelect All Exon v5+UTR (S04380219). These files were converted to gencode v21 coordinates by the UCSC liftover tool and files (https://genome-store.ucsc.edu/). Insert sizes were determined by Picard tools. Graphs were created in the R/Bioconductor environment (http://www. bioconductor.org/) in particular applying ggplot2 [18]. For comparability the evenness score [14] served as metrics for the uniform coverage in the target regions.
Variant calling was conducted by VarScan 2 [19] to identify mutations in four subclones of the HG-3 cell line with a minimum depth of 10× and 2× for an alternative allele. The four subclones were HG3CD5n_cl1 and HG3CD5p-cl7 for 130 bp peak insert sizes and HG3CD5n_cl41 and HG3CD5n_cl48 for 170 bp (see Table 1). Afterwards, the Variant Effect Predictor (release 77, http://www.ensembl.org/info/docs/tools/vep/ index.html) helped to concentrate on mutations in coding regions with severe consequences such as missense and stop gained effects.
Data were deposited at ArrayExpress under the accession number E-MTAB-4527.