Accurate self-correction of errors in long reads using de Bruijn graphs

Abstract Motivation New long read sequencing technologies, like PacBio SMRT and Oxford NanoPore, can produce sequencing reads up to 50 000 bp long but with an error rate of at least 15%. Reducing the error rate is necessary for subsequent utilization of the reads in, e.g. de novo genome assembly. The error correction problem has been tackled either by aligning the long reads against each other or by a hybrid approach that uses the more accurate short reads produced by second generation sequencing technologies to correct the long reads. Results We present an error correction method that uses long reads only. The method consists of two phases: first, we use an iterative alignment-free correction method based on de Bruijn graphs with increasing length of k-mers, and second, the corrected reads are further polished using long-distance dependencies that are found using multiple alignments. According to our experiments, the proposed method is the most accurate one relying on long reads only for read sets with high coverage. Furthermore, when the coverage of the read set is at least 75×, the throughput of the new method is at least 20% higher. Availability and Implementation LoRMA is freely available at http://www.cs.helsinki.fi/u/lmsalmel/LoRMA/.


Introduction
With the diminishing costs, high-throughput DNA sequencing has become a commonplace technology in biological research. Whereas the second generation sequencers produced short but quite accurate reads, new technologies such as Pacific Biosciences and Oxford NanoPore are producing reads up to 50 000 bp long but with an error rate at least 15%. Although the long reads have proven to be very helpful in applications like genome assembly (Koren and Philippy, 2015;Madoui et al., 2015), the error rate poses a challenge for the utilization of this data.
Many methods have been developed for correcting short reads (Laehnemann et al., 2016;Yang et al., 2013) but these methods are not directly applicable to the long reads because of their much higher error rate. Moreover, most research of short read error correction has concentrated on mismatches, the dominant error type in Illumina data, whereas in long reads indels are more common. Recently, several methods for error correction of long reads have also been developed. These methods fall into two categories: either the highly erroneous long reads are self-corrected by aligning them against each other, or a hybrid strategy is adopted in which the long reads are corrected using the accurate short reads that are assumed to be available. Most standalone error correction tools like proovread (Hackl et al., 2014), LoRDEC (Salmela and Rivals, 2014), LSC (Au et al., 2012) and Jabba (Miclotte et al., 2015) are hybrid methods. PBcR (Berlin et al., 2015;Koren et al., 2012) is a tool that can employ either the hybrid or self-correction strategy.
Most hybrid methods like PBcR, LSC and proovread are based on the mapping approach. They first map the short reads on the long reads and then correct the long reads according to a consensus built on the mapped short reads. PBcR extends this strategy to selfcorrection of PacBio reads by computing overlaps between the long reads using probabilistic locality-sensitive hashing and then correcting the reads according to a consensus built on the overlapping reads. As the mapping of short reads is time and memory consuming, LoRDEC avoids the mapping phase by building a de Bruijn graph (DBG) of the short reads and then threading the long reads through this graph to correct them. Jabba is a recent tool that is also based on building a DBG of short reads. While LoRDEC finds matches of complete k-mers in the long reads, Jabba searches for maximal exact matches between the k-mers and the long reads allowing it to use a larger k in the DBG.
In this paper, we present a self-correction method for long reads that is based on DBGs and multiple alignments. First our method performs initial correction that is similar to LoRDEC, but uses only long reads and performs iterative correction rounds with longer and longer k-mers. This phase considers only the local context of errors and hence it misses the long-distance dependency information available in the long reads. To capture such dependencies, the second phase of our method uses multiple alignments between carefully selected reads to further improve the error correction.
Our experiments show that our method is currently the most accurate one relying on long reads only. The error rate of the reads after our error correction is less than half of the error rate of reads corrected by PBcR using long reads only. Furthermore, when the coverage of the read set is at least 75Â, the size of the corrected read set of our method is at least 20% higher than for PBcR.

Overview of LoRDEC
LoRDEC (Salmela and Rivals, 2014) is a hybrid method for the error correction of long reads. It presents the short reads in a DBG and then maps the long reads to the graph. The DBG of a read set is a graph whose nodes are all k-mers occurring in the reads and there is an edge between two nodes if the corresponding k-mers overlap by k À 1 bases. LoRDEC classifies the k-mers of long reads as solid if they are in the DBG and weak otherwise. The correction then proceeds by replacing the weak areas of the long reads by solid ones. This is done by searching paths in the DBG between solid k-mers to bridge the weak areas between them. If several paths are found, the path with the shortest edit distance as compared to the weak region is chosen to be the correct sequence, which replaces the weak region of the long read. The weak heads and tails of the long reads are the extreme regions of the reads that are bordered by just one solid kmer in the beginning (resp. end) of the read. LoRDEC attempts to correct these regions by starting a path search from the solid k-mer and choosing a sequence that is as close as possible to the weak head or tail.
Repetitive regions of the genome can make the DBG tangled. The path search in these areas of the DBG can then become intractable. Therefore, LoRDEC employs a limit on the number of branches it explores during the search. If this limit is exceeded, LoRDEC checks if at least one path within the maximum allowed error rate has been found and then uses the best path found for correction. If no such path has been found, LoRDEC starts a path search similar to the correction of the head and tail of the read, to attempt a partial correction of the weak region.
Some segments of the long reads remain erroneous after the correction. LoRDEC outputs bases in upper case if at least one of the k-mers containing that base is solid, i.e. it occurs in the DBG of the short reads, and in lower case otherwise. For most applications, it is preferable to extract only the upper case regions of the sequences as the lower case bases are likely to contain errors.

Self-correction of long reads
In this section, we will show how an error correction procedure similar to LoRDEC can be used to iteratively correct long reads without short read data. We will use LoRDEC* to refer to LoRDEC in this long reads only mode. Then, we further describe a polishing method to improve the accuracy of correction. Figure 1 shows the workflow of our approach.

Iterative correction
To describe how LoRDEC can be adapted for self-correction of read sets, let Q be a set of long reads to be corrected, and let integer h be the abundancy threshold that is used in choosing the k-mers to the DBG. The correction procedure repeats for an increasing sequence k ¼ k 1 ; . . . ; k t the following steps 1-3: 1. Construct the DBG of set Q using as the nodes the k-mers that occur in Q at least h times; 2. Correct Q using the LoRDEC algorithm with this DBG; 3. Replace Q with the corrected Q.
After the final round, the regions of the reads identified as correct in the last iteration are extracted for further correction with the multiple alignment technique by LoRMA.
As the initial error level is assumed high, the above iterations have to start with a relatively small k ¼ k 1 . With a suitable abundancy threshold h, the DBG should then contain most of the correct k-mers (i.e. the k-mers of the target genome) and a few erroneous ones. Although path search over long weak regions may not be feasible because of strong branching of the DBG, shorter paths are likely to be found and hence, short weak regions can be corrected. After the first round, the correct regions in the reads have become longer because close-by correct regions have been merged whenever a path between them has been found, and thus, we can increase k. Then, with increasing ks, the DBG gets less tangled and the path search over the longer weak regions becomes feasible allowing for the correction of the complete reads. A similar iterative approach has previously been proposed for short read assembly (Bankevich et al., 2012;Peng et al., 2010).
When the path search is abandoned because of excessive branching, the original LoRDEC algorithm still uses the best path found so far to correct the region. Such a greedy strategy improves correction accuracy in a single run, but in the present iterative approach false corrections start to accumulate. Therefore, we make a correction only if it is guaranteed that the correction is the best one available in the DBG, i.e. all branches have been explored.
Abundancy threshold h controls the quality of the k-mers that are used for correction. In our experiments, we used a fixed Fig. 1. Workflow of error correction. LoRDEC* is first applied iteratively to the read set, with an increasing k. The corrected reads are further corrected by LoRMA, which uses multiple alignments to find long-distance dependencies in the reads threshold of h ¼ 4 in all iterations, meaning that the k-mers with less than four occurrences in the read set were considered erroneous.
To justify the value of h, we need to analyse how many times a fixed k-mer of the genome is expected to occur without any error in the reads. Then an h that is about one or two standard deviations below the expected value should give a DBG that contains the majority of the correct k-mers and not too many erroneous ones. We will use an analysis similar to Miclotte et al. (2015).
Let C '!k denote the coverage of a genomic k-mer by exact regions of length at least k. Here exact region refers to a continuous maximal error-free segment of some read in our read set. Figure 2 gives an example of exact regions. Let us add a $ character to the end of each read, and then consider the concatenation of all these reads. In this sequence, an exact region (of length 0 or more) ends either at an error or when encountering the $ character. Let n denote the number of reads, N the length of the concatenation of all reads and p the error rate. Then the probability for an exact region to end at a given position of the concatenated sequence is q ¼ ðpN þ nÞ=ðN þ nÞ. As the reads are long and the error rate is high, we have q % p. The length of the exact regions is distributed according to the geometric distribution GeomðqÞ , and therefore, the probability of an exact region to have length i is PðiÞ ¼ ð1 À qÞ i q. The expected number of exact regions is Nq. An exact region is maximal if it cannot be extended to the left or right. Let R i be the random variable denoting the number of maximal exact regions of length i. Then EðR i Þ ¼ NqPðiÞ ¼ Nq 2 ð1 À qÞ i .
Let C '¼i denote the coverage of a k-mer in the genome by maximal exact regions of length i, and let r i denote the number of maximal exact regions of length i. An exact region of length i, i ! k, covers a fixed genomic k-mer (i.e. the read with that exact region is read from the genomic segment containing that k-mer) if the region starts in the genome from the starting location of the k-mer or from some of the i À k locations before it. Assuming that the reads are randomly sampled from the genome, this happens with probability ði À k þ 1Þ=G, where G is the length of the genome. Therefore, C '¼i is distributed according to the binomial distribution Binðr i ; ði À k þ 1Þ=GÞ (independence of locations of exact regions is assumed), and the expected coverage of a genomic k-mer by maximal exact regions of length i is By the linearity of expectation, the expected coverage of a genomic k-mer by exact regions of length at least k is Because ði À k þ 1Þ=G is small, we can approximate the binomial distribution of C '¼i with the Poisson distribution. Therefore, Assuming that the coverages of a genomic k-mer by maximal exact regions of different lengths are independent, the variance of the coverage by exact regions of length at least k is Figure 3 illustrates EðC '!k Þ for various k and q % p, with 100Â original coverage of the target. Note that original coverage of the target genome by the read set is N/G. For the three datasets in our experiments (Table 1), with coverages 200Â, 208Â and 129Â, the expected coverage EðC '!k Þ has values 9.12, 9.48 and 5.89, respectively, for our initial k ¼ 19 and for our assumed error rate p ¼ 0.15. Hence, our adopted threshold h ¼ 4 is from 0.8 to 1.8 standard deviations below the expected coverage meaning that most of the correct k-mers should be distinguishable from the erroneous ones.

Polishing with multiple alignments
The error correction performed by LoRDEC* does not make use of long range information contained in the reads. In particular, approximate repeats of the target are collapsed in the DBG into a path with alternative branches. In practice, such repeat regions are corrected towards a copy of the repeat but not necessarily towards the correct copy. However, the correct copy is more likely uncovered because we choose the path that minimizes the edit distance between the weak region to be corrected and the sequence spelled out by the path. Therefore, if we have several reads from the same location, the majority of them are likely corrected towards the correct copy.
Our multiple alignment error correction exploits the long range similarity of reads by identifying the reads that are likely to originate from the same genomic location. If the reads contain a repeat area, the most abundant copy of the repeat present in the reads is likely the correct one. Then by aligning the reads with each other we can correct them towards this most abundant copy. The approach we use here bears some similarity to the method used in Coral (Salmela and Schrö der, 2011).
As preprocessing phase for the method, we build a DBG of all the reads using abundancy threshold h ¼ 1 to ensure that all k-mers present in the reads are indexed. Then we enumerate the simple paths of the DBG and find for each read the unique path that spells it out. Each such path is composed of non-overlapping unitig segments that have no branches. We call such segments the parts of a path. We associate to each path segment (i.e. a unitig path of the DBG) a set of triples describing the reads traversing that segment. Each triple consists of read id, part id and the direction of the read Fig. 2. Division of a read into maximal exact regions, shown as boxed areas. The shaded boxes give the regions that could cover a 4-mer Fig. 3. Expected coverage of a genomic k-mer by exact regions of length at least k for a read set with coverage 100Â for different error rates p on this path. Hence, the path for a read i consists of segments who have a triplet with i as the read id and with part id values 1, 2,. . ., the path being composed of these segments in the order of the part id value (Fig. 4). Using this information, it is now possible to reconstruct each read from the DBG except that the reads will be prefixed (suffixed) by the complete simple path that starts (ends) the read.
In the second phase of our method, we take the reads one by one and use the DBG to select reads that are similar to the current read. We follow the path for the current read and gather the set of reads sharing k-mers with it, which can be done using the triplets of the augmented DBG. Out of these reads, we then first select each read R such that the shared k-mers span at least 80% of the shorter one of the read R and the current read. Furthermore, out of these reads, we select those that share the most k-mers with the current read. We call this read set the friends of the current read. The number of selected reads is a parameter of our method (by default 7).
We then proceed to compute a multiple alignment of the current read and its friends. To keep the running time feasible, we use the same simple method as in Coral (Salmela and Schrö der, 2011). First, the current read is set to be the initial consensus. Then we take each friend of the current read one by one, align them against the current consensus using banded alignment, and finally update the consensus according to the alignment. Finally, we inspect every column of the multiple alignment and correct the current read towards the consensus if the consensus is supported by at least two reads.
We implemented the above procedure in a tool called Long Read Multiple Aligner (LoRMA) using the GATB library (Drezen et al., 2014) for the implementation of the DBG.

Experimental results
We ran experiments on three datasets that are detailed in Table 1. The simulated Escherichia coli dataset was generated with PBSIM (Ono et al., 2013) using the following parameters: mean accuracy 85%, average read length 10 000, and minimum read length 1000. The other two datasets are real data. Although our method works solely on the PacBio reads, the table also includes statistics of complementary Illumina reads that were used to compare our method against hybrid methods that need also short reads. All experiments were run on 32 GB RAM machines equipped with 8 cores.

Evaluation of the quality of error correction
In the simulated dataset, the genomic position where each read derives from is known. Therefore, the quality of error correction on the simulated dataset is evaluated by aligning the corrected read against the corresponding correct genomic sequence. We allow free deletions in the flanks of the corrected read because the tools trim regions they are not able to correct. To check if the corrected reads align to the correct genomic position, we aligned the corrected reads  on the reference genome with BLASR (Chaisson and Tesler, 2012) keeping only a single best alignment for each read. The following statistics were computed: • Size: The relative size of the corrected read set as compared to the original one. • Error rate: The number of substitutions, insertions and deletions divided by the length of the correct genomic sequence. • Correctly aligned: The relative number of reads that align to the same genomic position where the read derives from.
To evaluate the quality of error correction on the real datasets, we used BLASR (Chaisson and Tesler, 2012) to align the original and corrected reads on the reference genome. For each read, we used only a single best alignment because a correct read should only have one continuous alignment against the reference. Thus, chimeric reads will be only partially aligned. We computed the following statistics: • Size: The relative size of the corrected read set as compared to the original one. Together, these statistics measure three aspects of the quality of error correction. Size measures the throughput of the method. Aligned and error rate together measure the accuracy of correction. Finally genome coverage estimates if reads deriving from all regions of the genome are corrected.

Parameters of our method
We ran experiments on the real E.coli dataset to test the effect of parameters on the performance of our method. First, we tried several progressions of k in the first phase where LoRDEC* is run iteratively. We started all iterations with k ¼ 19 because given the high error rate of the data k must be small for correct k-mers to occur in the read data. The results of these experiments are presented in Table 2. With more iterations, the size of the corrected read set and the aligned proportion of reads decrease, but the aligned regions are more accurate. The decrease in the size of the corrected read set may be a result of better correction because PacBio reads have more insertions than deletions. However, the decrease in the aligned proportion of the reads may indicate some accumulation of false corrections. The runtime of the method increases with the number of iterations but later iterations take less time as the reads have already been partially corrected during the previous rounds. To balance out these effects, we chose to use a moderate number of iterations, i.e. k ¼ 19, 40, 61, by default, which also optimizes the error rate of the aligned regions.
LoRMA also builds a DBG of the reads and thus we need to specify k. We investigated the effect of the value of k on the E.coli dataset. Table 3 shows the effect of k on the performance of LoRMA. Because the DBG is only used to detect similar reads in LoRMA, the performance is not greatly affected by the choice of k. There is a slight decrease in the throughput of the method as k increases as well as a slight increase in runtime but these effects are very modest. For the rest of the experiments, we set k ¼ 19.
Another parameter of the method is the size of the set of friends of the current read (-friends parameter). We tested also the effect of this parameter on the E.coli dataset. As the optimal value of this  97 19,22,25,28,31,34,37,40,43,46 66.630 99.311 0.274 20.65 19,22,25,28,31,34,37,40,43,46, 49,52,55,58  parameter might depend on the coverage of the dataset, we created several subsets of this dataset with different coverage to investigate this. Table 4 shows the results of these experiments. We can see that the accuracy of the correction increases as the size of the friends set increases. However, for the dataset with the lowest coverage, 75Â, the coverage of the genome by the corrected reads decreases when the size of the friends set is increased indicating that lower coverage areas are not well corrected. We can also see that increasing the size of the friends set increases the running time of the method. To keep the running time reasonable, we decided to set the default value of the parameter at a fairly low value, 7.

Comparison against previous methods
We compared our new method against PBcR (Berlin et al., 2015;Koren et al., 2012) which is to the best of our knowledge, the only previous self-correction method for long reads, and LoRDEC (Salmela and Rivals, 2014), proovread (Hackl et al., 2014) and Jabba (Miclotte et al., 2015) which also use short complementary reads. Table 5 shows the results on the simulated dataset comparing our new method to PBcR using long reads only. Table 6 shows the results of the comparison of our new method against previous methods on the real datasets. In the following, we will use LoRDEC to refer to the hybrid correction method using also short reads and LoRDEC*þLoRMA for our new method in which LoRDEC* is run in long reads selfcorrection mode followed by LoRMA. PBcR pipeline from Celera Assembler version 8.3rc2 was run without the assembly phase and memory limited to 16 GB. PBcR was run both only using PacBio reads and by utilizing also the short read data. For PBcR utilizing also short read data, the PacBio reads were divided into three subsets each of which was corrected in its own run. Proovread v2.12 was run with the sequence/fastq files chunked to 20M as per the usage manual and used 16 mapping threads. LoRDEC used an abundancy threshold of 3 and k-mer size was set to 19 similar to the experiments by Salmela and Rivals (2014). Jabba 1.1.0 used k-mer size 31 and short output mode. LoRMA was run with 6 threads. The k-mer sizes for LoRDEC*þLoRMA iteration steps were chosen 19, 40 and 61. For proovread and LoRDEC, we present results for trimmed and split reads. Table 5 shows that on the simulated data both PBcR and LoRDEC*þLoRMA are able to correct most of the data. Our new method achieves a lower error rate and higher throughput. We see that the fraction of corrected reads aligning to the correct genomic position is lower for LoRDEC*þLoRMA than for PBcR when all reads are considered, which suggests that LoRDEC*þLoRMA tends to overcorrect some reads. However, for corrected reads longer than 2000 bp this difference disappears, and thus, we can conclude that the overcorrected reads are short. When compared to the other selfcorrection method, PBcR, our new tool has a higher throughput and produces more accurate results on both real datasets as shown in Table 6. Out of the hybrid methods, Jabba has a lower error rate than LoRDEC*þLoRMA but its throughput is lower. When compared to the other hybrid methods, LoRDEC*þLoRMA has comparable accuracy and throughput. All hybrid methods produce corrected reads that do not cover the whole E.coli reference, which could be a result of coverage bias in the Illumina data. On the yeast data proovread produced few corrected reads and thus the coverage of the corrected reads is very low. Table 6 shows that our method is slower and uses more memory than PBcR in self-correction mode but its disk usage is lower. On the E.coli dataset our new method is faster than proovread and PBcR utilising short read data but slower than LoRDEC, Jabba or PBcR using only PacBio data. On the yeast dataset, we are faster than PBcR in hybrid mode but slower than the others. On the E.coli and yeast datasets, LoRDEC*þLoRMA uses 45% and 37%, respectively, of its running time on LoRDEC* iterations. On both datasets, the error rate of the reads after LoRDEC* iterations and trimming was 0.5%.

The effect of coverage
Especially for larger genomes, it is of interest to know how much coverage is needed for the error correction to succeed. We investigated this by creating random subsets of the E.coli dataset with coverages 25Â, 50Â, 100Â and 150Â. We then ran our method and PBcR (Berlin et al., 2015;Koren et al., 2012) on these subsets to investigate the effect of coverage on the error correction performance. Table 7 shows the results of these experiments. The other tools, LoRDEC, Jabba and proovread, use also the complementary Illumina reads and the coverage of PacBio reads does not affect their performance.
When the coverage is high, the new method retains a larger proportion of the reads than PBcR and is more accurate, whereas when the coverage is low, PBcR retains more of the data and a larger proportion of it can be aligned. However, the error rate remains much lower for our new tool. The reads corrected by PBcR also cover a larger part of the reference when the coverage is low.

Conclusions
We have presented a new method for correcting long and highly erroneous sequencing reads. Our method shows that efficient alignment free methods can be applied to highly erroneous long read data. The current approach needs alignments to take into account the global context of errors. Reads corrected by the new method have an error rate less than half of the error rate of reads corrected by previous self-correction methods. Furthermore, the throughput of the new method is 20% higher than previous self-correction methods with read sets having coverage at least 75Â.
Recently several algorithms for updating the DBG instead of constructing it from scratch when k changes have been proposed (Boucher et al., 2015;Cazaux et al., 2014). However, these methods are not directly applicable to our method because also the read set changes when we run LoRDEC* iteratively on the long reads.
Our method works solely on the long reads, whereas many previous methods require also short accurate reads produced by e.g. Illumina sequencing, which can incorporate sequencing biases in PacBio reads. This could have very negative effect on sequence quality, especially since Illumina suffers from GC content bias and some context-dependent errors (Nakamura et al., 2011;Schirmer et al., 2015).
As further work, we plan to improve the method to scale up to mammalian size genomes. We will investigate a more compact representation of the path labels in the augmented DBG to replace the simple hash tables currently used. Construction of multiple alignment could also be improved by exploiting partial order alignments (Lee et al., 2002) which have been shown to work well with PacBio reads (Chin et al., 2013).
Another direction of further work is to investigate the applicability of the new method on long reads produced by the Oxford NanoPore MinION platform. Laver et al. (2015) have reported an error rate of 38.2% for this platform and they also observed some GC content bias. Both of these factors make the error correction problem more challenging, and therefore, it will be interesting to see a comparison of the methods on this data.