Haplotype-Aware Long-Read Error Correction

To construct MSAs, we adapted code from the Herro software [18]. We briefly summarize the approach here. Herro employs the well-known star alignment heuristic, where the target read is assumed to be the “center of the star”. The MSA is constructed progressively, i.e., adding each overlapping read to the MSA one by one using its precomputed pairwise alignment with the target read. For further details on the star alignment heuristic, see [17].

Once the full MSA is constructed, it is divided into contiguous, non-overlapping windows of fixed length $w$ (Figure 1). Each window represents an MSA of a substring of the target read and its overlapping substrings. Henceforth, each windowed MSA is processed independently. Within each window, alignments from overlapping reads which do not fully span the window are ignored (Figure 1). We use the same window length $w=4096$ as Herro.

If there are less than two alignments from overlapping reads in a window, which can happen occasionally in regions with low coverage, then Herro chops and removes the corresponding substring of the target read instead of leaving it uncorrected. We do the same in HALE. Note that this procedure may split the target read. The rationale is that the uncorrected portions of a read may negatively affect the quality of genome assembly, hence they are removed.

The count of overlapping substrings within a window is influenced by the sequencing coverage of the loci linked to the target read and how repetitive those regions are within the genome. To maintain computational efficiency of next steps, we retain only top-$n$ overlapping substrings in each MSA. Preference is given to those overlapping substrings that have higher sequence identity¹¹1Given a pairwise alignment of two sequences, sequence identity is defined as the number of matching bases over the number of alignment columns. with the target read. If fewer than $n$ overlaps are available, we retain all. In HALE, the default value of parameter $n$ is set to $20$.

Following the approach of Herro [18] and Hifiasm [2], we identify informative alignment columns in the MSA. A column is considered informative if at least two distinct characters from the set $\{A,T,C,G,-\}$ occur in that column with a minimum frequency of $3$. See Figure 3 for an illustration. Only a small fraction of alignment columns are informative in practice. Informative columns typically indicate either biological variation among the aligned substrings or systematic sequencing errors, e.g., homopolymer indel errors. In contrast, non-informative columns are those where alignments largely agree. These columns have limited value for identifying true overlaps in Problem 1. Thus, we remove all non-informative columns from the MSA and use the modified MSA as input to Problem 1.

A polynomial-time algorithm is unlikely for Problem 1 because we have established its hardness in Section 4. We use a direct brute-force approach to find an exact solution. In the MSA input to Problem 1, we have a vector $t$ in ${\{A,T,C,G,-\}}^d$ (representing a subsequence of target read) along with $n$ other vectors in ${\{A,T,C,G,-\}}^d$ (from the overlapping reads). Here, $d$ denotes the number of columns in a given MSA. We solve Problem 1 by exhaustively evaluating all $\left(\genfrac{}{}{0pt}{}{n}{k}\right)$ possible ways of selecting a subset of $k$ vectors from the total of $n$. For a given selection of $k$ vectors, an optimal value of the cluster center $c$ at the $i^{\text{th}}$ coordinate is the character that occurs most frequently at that position among the $k$ vectors and the target vector $t$. Although the runtime of this algorithm is $O\left(\left(\genfrac{}{}{0pt}{}{n}{k}\right)\times kd\right)$, this brute-force solution is still practical due to the reduction in the size of MSA at Steps 2 and 3.

The default value of parameter $k$ in HALE is set to $n/3$. This choice is based on our empirical observation that at least one-third of the overlapping substrings typically originate from the same haplotype as the target read (assuming ploidy $\le 2$). Choosing an appropriate value for $k$ involves a trade-off: a significantly larger $k$ can lead to the selection of reads with false overlaps, while a much smaller $k$ leads to an insufficient number of reads to reliably compute a consensus.

We introduce an additional optimization. Note that Problem 1 treats all MSA columns with equal importance. To further improve the quality of results, we solve a weighted version of Problem 1 defined below. Suppose that weight function $W$ assigns a non-negative weight $W(i)$ to each column $i\in \{1,2,\mathrm{\dots },d\}$ and function $match(\alpha ,\beta )$ returns $1$ if characters $\alpha$ and $\beta$ match, and $0$ otherwise.

Our motivation to introduce the weights is as follows. A majority of errors in long reads are insertions and deletions, especially in homopolymers regions [11]. Accordingly, informative columns where the most frequent and the second most frequent characters are from $\{A,T,C,G\}$ are likely to represent true biological variation. We assign a weight $1$ to such columns. In contrast, other informative columns where either the most frequent or the second most frequent character is “$-$” may correspond to either true indel variation or sequencing artifacts. We assign a weight $1/10$ to these columns. The choice of this parameter was based on our empirical observations. Developing a more principled approach for determining the weight parameters is left for future work.

Recall that before sampling alignment columns in the MSA (Step 3), we had a substring of the target read aligned with overlapping substrings from other reads, forming an MSA of length $w$, where $w$ is the window length parameter. The alignment of the substring from the target read corresponds to a vector in ${\{A,T,C,G,-\}}^w$. We update this vector as follows. For each non-informative column, where biological variation is not expected, we replace the character with the most frequent character in that column, which is equivalent to taking a majority vote. In all informative columns, we use the the optimal cluster center $c$, as determined by solving Problem 1, to update the characters. Finally, we generate the corrected version of the target read substring by removing any gap symbols (“$-$”) from the updated vector.

We utilized two sequencing datasets, one from a diploid HG002 human sample²²2https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/scratch/HG002/sequencing/hifirevio/, and the other from a haploid Col-0 Arabidopsis thaliana sample³³3https://ngdc.cncb.ac.cn/gsa/browse/CRA004538/CRR302668. Following the benchmarking approach used in Herro, we restricted our evaluation to reads that are at least $10$ kbp in length. Reads shorter than $10$ kbp were removed.

To reduce runtime of the experiments involving the human dataset, we sampled reads from chromosome 9 of the HG002 genome. Chromosome 9 contains a high abundance of satellite DNA sequences (long arrays of tandem repeats) which makes it one of the more challenging chromosomes to assemble [14]. To isolate reads from chromosome 9, we aligned all reads to HG002v1.1 assembly⁴⁴4https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/HG002/assemblies/hg002v1.1.fasta.gz (GCA_018852605.3, GCA_018852615.3) and retained those reads whose primary alignments mapped to either the paternal or maternal haplotype of chromosome 9. To evaluate the impact of sequencing coverage on error correction performance, we created two versions of HG002 and A. thaliana datasets, one with $40\times$ coverage and the other with $60\times$ coverage by randomly sampling reads.

We compared HALE against two state-of-the-art read correction tools, Herro [18] and Hifiasm [2]. Herro is a deep learning method that was primarily developed for correcting nanopore sequencing reads. Since a pretrained model for HiFi reads is not available in Herro, we used its existing model trained on Oxford Nanopore (R10) data. We assumed that the ONT-trained model generalizes to HiFi reads because HiFi reads are relatively more accurate and exhibit fewer systematic errors.

In addition to these baseline methods, we also compared HALE with two simplified correction methods. As a reminder, HALE performs a majority vote on non-informative alignment columns while using a more sophisticated technique based on solving Problem 1 in the informative alignment columns (Section 5). The naive methods employ the same heuristics described in Section 5, but bypass solving Problem 1 (Step 4).

In the first naive approach (Naive-1), a majority vote is performed across all alignment columns, including informative ones, during Step 5. In the second naive approach (Naive-2), the original values in the aligned target read are preserved at informative columns, while a majority vote is applied only to non-informative columns. This comparison allows us to isolate the benefit of solving Problem 1 as an optimization step in HALE.

For real sequencing datasets, there is no ground truth available on the errors made by the sequencing instrument. Thus, for both HG002 and A. thaliana datasets, we rely on high-quality references on which the reads can be mapped. For HG002 datasets, we use HG002v1.1 assembly (GCA_018852605.3, GCA_018852615.3) released by the Telomere-to-telomere consortium. For A. thaliana datasets, we use a high-quality assembly⁵⁵5https://ngdc.cncb.ac.cn/gwh/Assembly/21820/show from [21].

To evaluate read correction performance, we aligned raw reads as well as reads corrected by different methods to their respective references using Minimap2. If a read aligns to multiple regions, then only its primary alignment was considered. We used BamConcordance [22] to examine the alignments. BamConcordance generates an empirical log-scaled measure of base error probabilities (Q-concordance) for each read. The formula used for calculating Q-concordance is $-10{\log }_{10}(1-\text{identity})$, where identity refers to the sequence identity observed in the read alignment. As an example, Q-concordance $20$ corresponds to a sequence identity of $99\%$, indicating that $1\%$ of the bases are erroneous. In addition to the Q-concordance statistics, we also report the percentage of reads containing at least $x$ mismatch (similarly, $x$ indel) errors, for $x\in \{1,2,3\}$, before and after correction.

As noted earlier in Section 5, the approach used in HALE and Herro discards a segment of a read if it cannot be reliably corrected due to insufficient overlapping read support. For a fair comparison, we evaluated only those reads that were fully corrected (i.e., not chopped) by all methods. Without this, HALE and Herro would have an unfair advantage over Hifiasm.

We show a comparison of the accuracy of raw reads and those corrected by HALE, Herro, and Hifiasm, respectively, in Figure 4a. We achieved comparable accuracy as Herro and Hifiasm. All three correction methods led to more than a 10-fold improvement in read accuracy relative to the raw, uncorrected reads.

Next, we separately assessed mismatch and indel errors in the corrected reads. In Figure 4b, we highlight that the fraction of HG002 chr9 reads with one or more mismatch errors significantly drops after error correction. This result suggests that all three methods are effective in preserving biological variations in the reads. At $40\times$ sequencing depth, HALE exhibits a marginally higher percentage of reads containing at least $x$ mismatch errors ($x\in 1,2,3$) compared to Herro and Hifiasm. However, this difference diminishes at $60\times$ depth. We also tested HALE using reads from HG002 chromosomes 13 and 18, observing a consistent trend (Appendix Figure 6).

In the case of A. thaliana dataset, which is a haploid genome sample, HALE is again comparable to the best performing method (Figures 4d, 4e). Herro lags slightly behind HALE and Hifiasm here, possibly because its model was trained using a diploid sample and nanopore reads [18].

We also manually checked and visualized a few reads in IGV [16] where HALE failed to correct one or more errors. In some cases, we found that uneven sequencing coverage in some genomic regions resulted in MSA windows where fewer than $k$ true overlaps were present. As a result, our algorithm, constrained by the fixed cutoff of $k$, selected substrings from false overlaps for correction.

Solving Problem 1 (Section 3) is an important component of HALE as it enables us to preserve the haplotype-specific variation in reads. We compared HALE with two simplified versions of our algorithm, Naive-1 and Naive-2, respectively, which omit this optimization step.

Runtime comparisons between HALE, Herro, and Hifiasm are challenging for two main reasons. First, Hifiasm uses its own custom all-vs-all read overlapping algorithm that is tightly integrated inside its code, whereas HALE and Herro use an external tool, Minimap2 [7]. Second, Herro leverages GPU acceleration, whereas both HALE and Hifiasm run on CPUs. In the following, we report the wall-clock runtimes of all-vs-all read overlapping and error correction steps for HALE and Herro. For Hifiasm, we report the total end-to-end assembly time. We report the results using chr9 HG002 $60\times$ dataset.

Minimap2 used $26$ minutes to compute all-vs-all overlaps using $64$ CPU threads. HALE used $9$ minutes for error correction using the same number of threads. In comparison, Herro required $10$ minutes on a server with four NVIDIA A100 GPUs. Hifiasm completed the full assembly in $16$ minutes using $64$ CPU threads, with approximately $90\%$ of that time spent on read overlapping and error correction steps. Based on these results, we conclude that HALE is computationally more efficient than Herro. Hifiasm is the fastest overall due to its highly optimized read overlapping implementation. As all-vs-all overlap computation is a runtime bottleneck while using HALE, using a faster read overlapper would reduce the runtime.