HECTOR: a parallel multistage homopolymer spectrum based error corrector for 454 sequencing data

k-hopo Definition and encoding

HECTOR represents a homopolymer as a pair (R X), where R is the homopolymer length and X ∈ Σ = {A, C, G, T, N}, and encodes a read as a string of homopolymer pairs. For example, the read S = AACCCCCGGGT will be encoded as (2, A) (5, C) (3, G) (1, T). In this case, all of the 3-hopos for S are (2, A) (5, C) (3, G) and (5, C) (3, G) (1, T). In HECTOR, we represent each homopolymer pair as a byte and then calculate the value using the following formula: (128*L + |∑|*R + B - 4), where B is an encoded nucleotide (A = 0, C = 1, G = 2, T = 3, N = 4) and L is the case flag (1 if the nucleotide is lowercase, 0 otherwise). Due to the 1 byte encoding implementation, the value of R has range between 1 and 25.

Homopolymer spectrum construction

HECTOR encodes all reads as strings of homopolymer pairs. The homopolymer spectrum is constructed from all of the encoded reads where the counting of k-hopos is performed using Bloom filters and hash tables and is parallelized using multi-threading. In our implementation, the default value of k is 21. This value is obtained empirically and is adapted from Musket [21].

The first stage filters out as many unique k-hopos as possible using a Bloom filter and stores all non-unique k-hopos in a hash table. However some unique k-hopos are likely to still exist in the hash table due to the false positive probability of a Bloom filter. In this stage, the masters fetch reads in parallel from the input file and then distribute all the k-hopos in the reads to the slaves. The accesses to the file are mutually exclusive and are guaranteed by locks. Hashing is used to distribute each k-hopo to its respective destination slave to ensure a good load balance. Once the destination slave is determined, the canonical k-hopo is transferred to the slave through the corresponding communication channels. The canonical k-hopo is the smaller numerical representation of the k-hopo and its reverse complement. Each slave holds a local Bloom filter and a local hash table to capture all non-unique k-hopos as well as filter out most unique k-hopos. When a k-hopo arrives, the slave performs membership look up in the local Bloom filter for the k-hopos. If the k-hopo queried exists in the Bloom filter, the slave inserts it in the local hash table because it is likely to have more than one occurrence. Otherwise, it is inserted in the Bloom filter. At the end of this stage, all non-unique k-hopos are stored in all local hash tables of all slaves, and the number of unique k-hopos occasionally existing in all hash tables depends on the false positive probability rate of all local Bloom filters. No synchronization between the slave threads is required in this stage due to the independence of the local Bloom filters and hash tables, which greatly benefits efficiency.

The second stage computes the multiplicity of each k-hopo in the hash table to determine the unique hopos that are still stored in the hash table. In this stage the masters follow the same procedure as in the previous stage and all slaves listen to their corresponding communication channels and wait for the arrival of k-hopos. Once a k-hopo arrives, a slave queries the existence of the k-hopo in the local hash table and then increments the multiplicity if it exists. At the end of this stage, each slave holds the multiplicity of each k-hopo stored in its hash table, which is used to determine the uniqueness of a k-hopo.

The third stage removes all unique k-hopos from the hash table leaving all non-unique ones in place. In this stage, each slave deletes all unique k-hopos from its hash table while the masters are idle. At the end of this stage, each slave holds a partition of the set of all non-unique k-hopos in the input reads.

The last stage determines the coverage cut-off for the homopolymer spectrum from the k-hopo coverage histogram. A k-hopo coverage histogram illustrates a mixture of two distributions: one for the coverage of likely correct k-hopos and the other for spurious k-hopos. The k-hopos distributed at the right of the valley are supposed to be true k-hopos (trusted) and the ones on the opposite side to be spurious (untrusted). A homopolymer spectrum is different from a k-mer spectrum i.e. the multiplicity along the x-axis represents a range of actual sequence length. A comparison between a homopolymer spectrum and a k-mer spectrum of SRR000868, SRR000870 and SRR639330 reads is shown in Figure 2. In Illumina datasets where the reads have equal length, the distributions of the k-mer spectrum are normally bimodal. However, due to the irregular read length nature of 454 datasets, the distribution of the k-mer spectrum tends to be unimodal while the distribution of the homopolymer spectrum is normally bimodal. Therefore, by encoding the spectrum construction in homopolymer space, the homopolymer spectrum obtained is analogous to the k-mer spectrum in Musket. As seen in Figure 2a and b the homopolymer spectrums of the SRR000868 and SRR000870 datasets are similar because their datasets were sequenced in the same experiment and have similar coverage (11.7× and 11×), respectively. Figure 2c shows that the SRR639330 dataset has a different homopolymer spectrum because it has a higher coverage (69.8×) and was sequenced in a different experiment. However, in general, all the homopolymer spectrums are bimodal. HECTOR then adopts the coverage cut-off selection method of Musket by choosing the multiplicity corresponding to the smallest density around the valley as the coverage cut-off. As shown in Figure 2, the cut-off values for SRR000868, SRR000870 and SRR639330 datasets are 3, 3 and 6, respectively. These values correspond to the multiplicity associated with the smallest density around the valley in the homopolymer spectrum for each of the datasets, which are 91,923, 100,408 and 38,081, respectively. In addition, HECTOR provides a parameter to allow users to specify the cut-off.

Error correction

Figure 3 shows how sequencing errors i.e. insertions, deletions and substitutions, are represented in the context of k-hopos, given a k-hopo S _k . Due to the representation of the k-hopos in (R, X) format described in the previous section, identifying a homopolymer-length error can then be represented as a difference in R. Therefore, whatever the type of error is, only a single k-hopo pair is actually changed. Thus, explaining why identification of errors in our method can be very efficient. We also ignore error corrections that occur in the first or the last homopolymer of a read. This is because the actual polynucleotide in the reference genome can be much longer than the homopolymer at the beginning or end of a read. Thus, it is impossible to tell if the read should have been longer or shorter by a single base or two.

The examples in Figure 3 are not the only conceivable scenarios of substitution and run-length errors. There are others that could change k-hopos more drastically. We argue in Additional file 1 why these scenarios do generally not occur in 454 data.

The error correction phase of HECTOR consists of three stages i.e. two-sided conservative correction, one-sided aggressive correction and voting-based refinement.

Two-sided conservative correction

In the two-sided correction we assume that there is at most one error in any k-hopo of a read. This assumption is later relaxed in the one-sided aggressive correction. The two-sided correction starts with the classification of trusted and untrusted bases for a read. A base is considered to be trusted if it is covered by any trusted k-hopo. Otherwise, it is considered to be untrusted and is considered to be a potential homopolymer-length error. For a sequencing error occurring at position i of a read of l bases, it causes up to min{k, i, l - i} erroneous k- hopos. HECTOR only evaluates both the leftmost and the rightmost k-hopos that cover position i on the read, instead of all possible k k-hopos that cover position i, thus, significantly improving speed and avoiding high computational overhead. For each base, the correction is made only if an alternative is found to make both the leftmost and the rightmost k-hopos trusted. Otherwise, if more than one alternative is found, the base will be kept unchanged as a result of ambiguity. For a substitution error, the alternative is a homopolymer of length 1 with a different nucleotide. For an indel error, the alternative allows up to 3 bases deletion and up to 3 bases insertion of the same nucleotide. For a read, the two-sided correction will be executed for a fixed number of iterations or until no base change is made.

One-sided aggressive correction

Since the two-sided conservative correction assumes that there is only one error in a single k-hopo it is inadequate in cases where more than one error occurs in a single k-hopo. Therefore, we utilize a one-sided correction to aggressively correct errors in these cases. The core idea is similar to the one described in [21] but we are looking for an insertion or deletion alternative, instead of a substitution alternative. Given a read H, define H _i to denote the base at position i of H, and H _i,k to denote the k-hopo starting at position i. If H _i,k is trusted, but H _i+1,k is untrusted, the homopolymer H _i+k is likely to be a sequencing error. Our one-sided correction aims to find an alternative of H _i+k that yields a trusted H _i+1,k . Unlike the two-sided correction, this correction selects the alternative, which makes the resulting trusted k-hopo H _i+1,k have the largest multiplicity, if more than one alternative is found.

The one-sided correction begins with the location of trusted regions for a read. A trusted region means that every base in this region is trusted thus having a minimum length of k. For each trusted region, error corrections are conducted towards each of the two orientations on the read. For each orientation, the error correction process does not stop until it either reaches a neighbouring trusted region or fails to correct the current base.

This correction approach is effective but has the drawback that it strictly relies on the correctness of k-hopo H _i,k . Thus if H _i,k does contain sequencing errors but is deemed to be trusted, this one-sided correction is possible to cause cumulative incorrect corrections, mutating a series of correct bases to incorrect ones. Look-ahead validation and voting-based refinement techniques are used to reduce this adverse effect. The look-ahead validation evaluates the trustworthiness of a predefined maximal number of neighbouring k-hopos that cover the base position at which a sequencing error likely occurs. If all evaluated k-hopos are trusted for a certain alternative on that position, this alternative is reserved as one potential correction. The voting-based refinement is described below. There is also a constraint on the number of corrections that are allowed in any k-hopo of a read. The default value of the constraint is 4. During the correcting process, we track the number of corrections that have been made in any k-hopo. If the number of corrections of a k-hopo exceeds the defined constraint, all corrections made in the k-hopo will be disregarded.

Voting-based refinement

The voting-based refinement method used is similar to the voting algorithm originally used in DecGPU [16]. The voting algorithm attempts to find the correct base by replacing all possible bases at each position of the k-hopo and checking the solidities of the resulting k-hopos. This approach introduces the fewest new errors even though it does not correct as many errors as other correctors. However, the objective of this approach is to reduce of the number of new errors from the one-sided aggressive correction.

Parallelization strategy

Our parallelization strategy uses the master–slave model a typical parallelization paradigm in which masters are dedicated to task distribution, and slaves are assigned to work on individual tasks. Typically in many applications, the master–slave model is implemented using a single master and multiple slaves. In HECTOR, multiple masters are used to avoid the bottleneck caused by the task distribution as the number of slaves grows larger. In general, the model is configured to have more slaves than masters, with a master-to-slave ratio chosen to be 1:3. This type of multiple masters, multiple slaves parallelization has also been implemented by Musket [21]. A hybrid combination of Pthreads and OpenMP parallel programming models is used to implement the master–slave model.

Experimental design

The performance of HECTOR is evaluated in terms of error correction quality, runtime, and parallel scalability. All experiments were conducted on a workstation with two Intel Xeon X5650 hex-core 2.66 GHz CPUs and 96 GB RAM, running Linux (Ubuntu 12.04). Both real and simulated reads are used in the experiments. The real reads are taken from the NCBI Sequence Read Archive (SRA), i.e. SRR000868, SRR000870, SRR639330 and SRR957993. The simulated reads are generated by Mason [25], with default parameters, using E. coli UTI89 and E. coli O104:H4 as reference genomes. The performance of HECTOR is then compared to Coral (v. 1.4) and Acacia (v.1.52) using default parameters.

Correction quality evaluations

In order to evaluate correction quality, we have to distinguish the correct from the erroneous bases in the experimental data. Errors in the sequence reads are identified by using the CUSHAW2 [9] mapping program to align reads to the reference genome. Only uniquely mapped reads with no clippings are considered, with the errors given by the differences to the reference genome. This is common practice since genomic variants such as single nucleotide polymorphisms (SNPs) can be defined as errors, and ambiguously mapped reads can lead to false classifications of bases as well [13, 14, 26]. However, the number of SNPs compared to sequencing errors is insignificant, so that they can be safely ignored for the purposes of experimental validation [27, 28]. In addition, this is a general problem for applications like error correction, assembly, or mapping. Thus, they pose the same problem for any corrector, so that the relative results are not affected.

We have explored SHRiMP2 [29] as a potential mapping algorithm apart from CUSHAW2 and found that the reads mapped by both tools are practically equivalent. We conclude that the evaluation of error correction performance of Coral and HECTOR is independent of the mapper used to establish the gold standard of erroneous bases. For more details, please refer to Additional file 2.

The quality of the error correction is then evaluated in terms of recall, specificity, gain, precision, and F-score. We define TP (True Positive) as the number of erroneous bases that are successfully corrected, FP (False Positive) as the number of newly introduced errors, i.e. the number of correct bases that are changed to be erroneous, TN (True Negative) as the number of correct bases that remain unchanged and FN (False Negative) as the number of erroneous bases that remain undetected. Please note that homopolymers that fall in the start and end of reads are not counted to obtain the statistics. Recall is calculated as TP/(TP + FN), specificity as TN/(FP + TN), gain as (TP - FP)/(TP + FN), precision as TP/(TP + FP) and F-score as 2 × precision × recall/(precision + recall), respectively. All the recall, specificity, gain, precision, and F-score values in related tables have been multiplied by 100 and all best values have been highlighted in bold.

Evaluation on simulated datasets

Evaluation on real datasets

Discussion of performance on simulated and real datasets

There are several factors that cause differences in performance between real and simulated data. Firstly, the coverage in simulated data is more uniform, whereas real data tends to follow a more variable distribution. This leads to areas in the genome that are now correctable; it can also cause false corrections, because correct sequence gets erroneously classified as errors. Secondly, the error profiles of real and simulated data can be different leading to further differences in correction performance. Overall, for the simulated datasets, HECTOR shows a consistent performance regardless of coverage while the performance of Coral improves as the coverage increases. This is consistent with the performance of the real datasets, in which HECTOR is superior for the dataset with the least coverage (SRR000868) while the performance of Coral gets better as the coverage of the datasets increases. The biggest performance difference between Coral and HECTOR occurs on the SRR639330 dataset, which has an extreme coverage of almost 70 × .

Performance analysis based on homopolymer length

The correction quality performance is then further analysed based on the length of homopolymers. The sensitivity and specificity of both HECTOR and Coral for homopolymer lengths of 1, 2, 3, 4, 5, 6, 7 and larger than 7 on the SRR000868 dataset are shown in Figure 4. The homopolymers are classified based on the length of the homopolymers in the original dataset. For example, if a homopolymer of length 3 gets corrected to any length, then this is classified as a correction for a homopolymer of length 3.

As shown in Figure 4a, Coral shows better sensitivity for homopolymers of length one and shows an irregular trend while the sensitivity of HECTOR improves as the homopolymer length increases. Figure 4b shows that in terms of specificity, both Coral and HECTOR show comparable performance for homopolymer length of 1, 2 and 3. As the homopolymer length gets longer, Coral shows better specificity than HECTOR. It should be noted that typically, the longer the length of the homopolymer, the less occurrence it has in the sequences. Thus, there is a higher probability of homopolymers with very long length which occur less than the required cut-off value of the homopolymer spectrum. These homopolymers will be classified as untrusted and therefore will be falsely corrected, which is one of the reasons why the specificity of HECTOR goes down as the length of the homopolymer increases.

Run time and parallel scalability evaluations

In addition to correction quality, runtime and parallel scalability are important factors that must be taken into account, especially for large-scale datasets. We have evaluated the runtime and parallel scalability of both HECTOR and Coral using the aforementioned real datasets using 2, 4, 8 and 12 threads. Runtimes are measured in wall clock time. Figure 5 shows the runtime comparison between HECTOR and Coral. HECTOR is superior to Coral for all cases. On average, HECTOR is about 3.7× faster than Coral.

The run times of HECTOR are vastly superior compared to Acacia. For the SRR000868, SRR000870, SRR639330 and SRR957993 datasets, Acacia has a run time of 25, 24, 185 and 22 hours, respectively. In comparison, HECTOR has a run time of 136.9, 114.6, 459.7 and 231.3 seconds, respectively, using 2 threads. The run time of Acacia is not shown in Figure 5 because of the huge run time difference compared to Coral and HECTOR. In addition, Acacia is also only a sequential program, so it is not possible to measure the parallel run time with multiple threads.

The parallel scalability of HECTOR and Coral are measured by varying the number of CPU threads used. As HECTOR employs a master–slave model, all of its speed-ups are calculated against the runtime with two threads. Figure 6 illustrates the speedups of HECTOR and Coral with different number of threads. HECTOR showed a better parallel scalability compared to Coral.