IMPACT: Interval-based Multi-pass Proteomic Alignment with Constant Traceback

TL;DR

This paper explores adapting Darwin's hardware-accelerated DNA alignment algorithms for protein sequence alignment, proposing a multi-pass GACT variant that improves accuracy by 14% on average.

Contribution

It demonstrates the feasibility of applying Darwin's algorithms to protein sequences and introduces a novel multi-pass GACT algorithm that enhances alignment accuracy.

Findings

Multi-pass GACT improves alignment scores by 14% on average.

Darwin's algorithms can be adapted for protein sequence alignment.

The proposed method maintains efficient hardware implementation.

Abstract

Darwin is a genomics co-processor that achieved a 15000x acceleration on long read assembly through innovative hardware and algorithm co-design. Darwins algorithms and hardware implementation were specifically designed for DNA analysis pipelines. This paper analyzes the feasibility of applying Darwins algorithms to the problem of protein sequence alignment. In addition to a behavioral analysis of Darwin when aligning proteins, we propose an algorithmic improvement to Darwins alignment algorithm, GACT, in the form of a multi-pass variant that increases its accuracy on protein sequence alignment. Concretely, our proposed multi-pass variant of GACT achieves on average 14\% better alignment scores.

Figures

Full text

1 Introduction

Biological sequences are compared through sequence alignment111We use this term to refer to dynamic programming “exact” alignment [2] between two sequences, rather than mapping of reads to a reference, which is sometimes also called alignment.: Given two strings of different lengths, called the query $Q$ and the reference $R$, alignment consists of finding the location of a series of gaps to insert into both sequences to make them align, subject to a scoring function that rewards matches and penalizes mismatches or gaps.

Smith & Waterman (S&W) proposed a dynamic programming algorithm capable of finding an optimal alignment in time and space $O(|Q|\times|R|)$ [2]. The algorithm has two phases: (1) computing a $|Q|\times|R|$ scoring matrix representing the alignment scores of all possible sub-sequences of $Q$ and $R$, and (2) tracing back from the highest-scoring cell in the matrix through a combination of 3 directions ($\leftarrow,\nwarrow,\uparrow$) until a cell with score zero is found. The sequence of cells encountered during the traceback forms the optimal alignment of the two sequences. The S&W algorithm is compute-intensive and has been the focus of many acceleration efforts such as parallelization [3, 4] and hardware acceleration [5, 6].

Recent efforts include Darwin [1], a system that accelerates alignment through both hardware and algorithmic innovations. Crucially, Darwin provides mechanisms for hardware-accelerated alignment with traceback support (providing not only the alignment score, but also the sequence of edits through which the score was achieved) using a constant amount of memory. This allows Darwin to align arbitrary-length sequences using space-limited hardware resources, achieving orders of magnitude speedups on reference guided and de-novo assembly operations. Darwin uses D-SOFT (Diagonal-band Seed Overlapping based Filtration Technique) to reduce the dynamic programming search space, followed by GACT (Genome Alignment using Constant memory Traceback) to perform the alignment. D-SOFT finds good exact-matching seeds or subsequences between the two sequences being aligned, which GACT then uses as its starting point to perform the extension of the alignment, extending from both sides of the seed until an alignment is found. Both component algorithms of Darwin are heuristic, and while Darwin is not guaranteed to find the optimal alignment, it does so in practice most of the time.

Darwin’s algorithms were designed and evaluated with DNA sequences, but fast alignment is also needed for protein sequence analysis. Proteins are different from DNA — they consist of a much larger alphabet, they have highly variable lengths, and their alignment substitution matrices are non-uniform (that is, each alphabet character pair is not a simple match or mismatch, but rather encodes a substitution probability).

These differences lead to some surprising and sometimes counter-intuitive behavior from Darwin that results in suboptimal performance when applied as-is to protein sequences. In this paper, we provide an analysis of this behavior and investigate the underlying causes. In addition, we propose a new multi-pass version of GACT that helps improve alignment scores for proteins. More specifically, this paper makes the following contributions:

• –

A detailed analysis of the performance of Darwin when used to align protein sequences.

• –

Based on this analysis, a proposal and evaluation of an algorithmic augmentation to GACT to retain near-optimal alignment results for proteins.

The rest of this paper is organized as follows: Section 2 describes the GACT algorithm; Section 3 discusses situations where GACT fails to provide optimal alignments for protein sequences; Section 4 introduces a multi-pass version of GACT that achieves near-optimal alignment of “straight-line” sequences; Section 5 describes future work related to joining intermediate multi-pass GACT results together; Section 6 concludes.

2 The GACT aligner

GACT takes a candidate seed produced by D-SOFT and left- and right-extends it in a tiled fashion using the S&W algorithm. To simplify the following description, our discussion assumes GACT is seeded with the max-scoring cell of the S&W scoring matrix and that only left-extension is performed. GACT has two novelties compared to prior hardware S&W implementations:

• –

GACT computes a single tile of the S&W matrix and traces it back before moving to the next tile, whereas other implementations only start a traceback after computing the complete S&W matrix. The advantage is that the amount of traceback state to be stored in on-chip memory is limited by the tile size $T$.

• –

GACT does not trace back complete tiles, but rather stops at a distance $O$ from one’s sides and uses this position as the starting offset of the next tile to be computed. The end result is that neighboring tiles overlap and some cells are computed twice, possibly to different values.

GACT is not an optimal alignment algorithm since the value of a cell, as defined by S&W, depends on its 3 neighbors ($\leftarrow,\nwarrow,\uparrow$), but these neighbors do not exist for the top and left borders of a tile. However, GACT’s use of overlapping tiles enables recomputed cells to obtain higher values than those computed in their first tile222Cells in the overlapping region of the first tile are close to the tile borders and are therefore very sensitive to border values. The same cells in the second tile are far from the borders and are less sensitive to its values.. Overlapping tiles therefore allow a traceback to follow a path that would otherwise not have existed if a cell value was not recomputed. This behavior can be seen at the interfaces between the tiles in Figure 1.

GACT uses the unmodified S&W algorithm in each tile. It is efficiently implementable in hardware using a systolic architecture of processing elements (PEs) that computes tile cells in a wavefront fashion.

3 Aligning Proteins with GACT

The recurrence relation in the S&W algorithm involves a table lookup into a substitution matrix. This symmetric matrix determines the weight attributed to each character pair in the S&W algorithm and is the primary reason for which GACT does not always achieve optimal alignments on protein sequences. There are several important differences between DNA and protein substitution matrices. First, while DNA is formed from nucleotides with an alphabet size of just 4 characters: {A, C, G, T}, proteins are formed from sequences of amino-acids, building blocks with a much larger alphabet size of 20, leading to much larger matrices. Second, DNA sequences have simple matrices that attribute homogeneous scores to all matches and mismatches (e.g. +2 for matches, -1 for mismatches), whereas protein matrices are heterogeneous and weigh not only matches positively, but also certain mismatches. Finally, entries in protein substitution matrices vary wildly in magnitude, with certain matching amino-acid pairs weighing up to $5\times$ more than others (c.f. Figure 2.)

Put together, these differences make it difficult for GACT to converge back to the path taken by the optimal protein alignment when it makes an erroneous local traceback decision.

3.1 More tile overlaps are not always better

Darwin showed that (1) increasing $T$ while keeping $O$ constant, and (2) increasing $O$ while keeping $T$ constant both result in an increase in the number of observed alignments that are optimal (with best results obtained when both components are increased).

Darwin’s observations, however, do not always hold for protein sequences. Consider the following experiment: Given a $(T,O)$ configuration for which there exists some proteins that do not align optimally, proceed to incrementally increase $O$ (while keeping $T$ constant) until the number of proteins that align optimally increases. Intuitively, one expects that increasing $O$ will always increase the number of relatively-high scoring cells in an overlapped region, and therefore provide more confident values for the traceback to use. So increasing $O$ should enable GACT to find the same or better alignments than at lower overlap margins. However, we actually observe that the total number of optimal alignments obtained by GACT by following this procedure decreases! This implies some of the alignments that were previously optimal with a smaller value of $O$ now diverge with larger values.

Figure 3 shows one such example where an optimal alignment is achieved with $(T,O)=(32,6)$, but where the traceback ends prematurely with $(T,O)=(32,7)$. This phenomenon can be explained by analyzing the values of the upper-left cells encompassed by the 2nd tile of Figure 3(a) that are absent in the same tile of Figure 3(b). The highly variable values in the substitution matrix, coupled with the relatively low scores in GACT tiles, means that displacing a tile by just 1 cell can cause it to miss a high-scoring substitution, and therefore trigger a sharp drop in tile cell values that ultimately cause a traceback to end prematurely at a zero cell.

The conclusion is that a GACT suffers excessively for small local alignment mistakes on protein sequences, a situation that can also occur for DNA sequences, but which is far less probable due to the low amplitude scores in DNA substitution matrices that smooths error propagation.

3.2 Hardware resource requirements

The impressive alignment throughput of Darwin relies on two mechanisms: (1) each GACT module is an array of PEs that compute S&W cells in a wavefront progression; and (2) GACT modules are replicated 64 times to handle incoming seed hits provided by D-SOFT in parallel.

The traditional systolic architecture used to implement S&W is economical in hardware resources when used for DNA sequences as the substitution “matrix” can simply be replaced with a comparator, two registers to store the match reward or mismatch/gap penalty, and finally a multiplexor that selects one of the two values.

This optimized implementation cannot be used when protein sequences are being aligned, as the substitution matrix is too large to be modeled with simple logic/registers and must be stored in an on-chip memory per-PE for fast lookup. It is therefore not possible to put as many PEs in each GACT module given the increased on-chip memory requirement, and thus processing large tiles becomes linearly slower for protein sequences.

4 Improving the GACT heuristic

Based on our observations, we need to reduce:

• –

Premature traceback termination due to sensitivity to GACT’s placement of the current tile;

• –

Traceback divergence due to an erroneous traceback in a previous tile that causes the current tile from being misplaced.

The last point stems from the fact that GACT uses its current traceback to help place the next tile, but this assumes that the traceback leading up to the tile borders is correct. This assumption is one of the causes of “forward” traceback divergence, but we show that the phenomenon can be attenuated with some additional computation.

All issues with GACT arise from low scores being present in the S&W scoring matrix. A natural solution to low-scoring cells would be to increase their values when we detect they are falling too low, however we cannot perform such an operation since one does not know what “too low” is333Long alignments would have a high “low threshold”, whereas short alignments would have a low “low threshold”, but we do not know in advance whether the alignment is long or short.. Instead we propose the following 6-step multi-pass algorithm that increases confidence in the placement of the next tile by increasing confidence that the traceback of the previous tile is correct, while simultaneously retaining GACT’s key contribution (constant traceback memory usage):

Compute cells of tile 1; 2. 2.

Perform traceback of tile 1 until distance $O$ from border; 3. 3.

Compute cells of tile 2 (no traceback); 4. 4.

Recompute cells of tile 1, taking into account elevated values from overlapped region with tile 2; 5. 5.

Perform final traceback of tile 1 until distance $O$ from border; 6. 6.

Advance to the next tile.

Figure 4(a) shows an example where the original GACT algorithm suffers from an erroneous traceback in a previous tile, which then causes a future tile to be misplaced (by just 1 cell) and ultimately lead to premature traceback termination. In turn, Figure 4(b) shows how the proposed multi-pass variant of GACT handles such a situation by raising cell scores enough such that new traceback “routes” that were previously unseen become available.

Note that the algorithm is completely determined by the dataset on which it operates and does not rely on any pre-determined “low threshold” for cell boosting. Furthermore, one does not need to recompute the values of cells that were overlapped with tile 2 as they can be forwarded from tile 2’s on-chip memory.

When evaluated on a small bacterial protein dataset containing roughly 25K sequences, we found that our multi-pass variant obtains on average 14% better alignment scores when used under identical parameters as the original GACT algorithm, with relative score improvements ranging from 0% (for identical alignments as GACT) up to 29000% (for alignments where the multi-pass variant allows long sequences to align much better).

Note that tile-based alignment can be accelerated if (1) tiles are smaller, i.e. quadratically fewer cells to compute, and (2) tiles are laid out as far as possible from each other, reducing cell value recomputation. Though it is true that this multi-pass variant of GACT roughly triples the compute requirements of the original algorithm in terms of number of tiles computed, we believe that the increased confidence in the tracebacks produced allows us to minimize tile sizes and overlap margins while obtaining similar results, therefore enabling a quadratic reduction in the number of cells that need to be computed.

Ultimately, GACT remains a heuristic algorithm and while our proposed enhancement alleviates some of its issues, we do not eliminate them. We discuss possible solutions in Section 5.

5 Future Work

5.1 Problematic sequences

One notable issue that neither the original GACT algorithm nor our multi-pass variant solves is that of traceback diverging due to the presence of a single long deletion or insertion in the optimal alignment. Indeed, it is impossible for any tile-based alignment algorithm to find the optimal traceback if a deletion or insertion has a length longer than the chosen tile size as scores would obviously reach zero. Furthermore, increasing overlap margins is no solution here as successively placing a tile where it goes through such “vertical” or “horizontal” traceback strips will ultimately cause the tile to no longer cover the region of interest as tiles are always placed diagonally from each other by at least one cell. The only way for a tiled alignment to work on such problematic sequences is if the tile size is larger than the longest deletion or insertion, but one cannot know this value in advance. Increasing tile sizes also has a negative effect on alignment speed.

GACT and our multi-pass variant are therefore well-suited for finding “diagonal” alignments in a S&W score matrix, but not for long “vertical” or “horizontal” ones. This brings us to our discussion of D-SOFT, Darwin’s filtering algorithm, and how we can improve it to help solve the long deletion/insertion problem.

5.2 Interval-based alignment

Myers’ proposed an algorithm for finding seed hits that have a high probability of being “close” to the optimal alignment in the presence of noisy long reads [7]. Roughly, the idea is to (1) find exact matches of seeds of size $k$ from the query in the reference, then (2) count the number of non-overlapping seed hits in diagonal bands, and finally (3) select a band that passes a certain threshold. Darwin’s D-SOFT implementation translates Myers’ hardware-unfriendly algorithm into one more suitable for hardware acceleration and is able to provide quality seeds to the GACT arrays for extension.

However, D-SOFT is unsuitable for use with protein sequences for one main reason: While the probability of finding sequences of 10-15 identical nucleotides in long DNA sequences is very high, the probability of finding a sequence of even three common amino-acids between two proteins is almost zero. Therefore any seeding algorithm that relies on exact matches to find seeds will experience seed starvation.

Since exact string matching does not work, we envision a filtering algorithm that uses the same subsequence matching method as S&W. Concretely, we propose to select multiple small samples from the query and the reference sequences, and to perform independent S&W computations on these regions to obtain a single value for each, effectively that of the max-scoring cell. We can then find diagonal sequences of high-scoring “blocks” in this “down-sampled” space and issue them as seeds to the GACT arrays for alignment. Finally, these aligned segments would need to be stitched together to form a complete alignment without affecting their alignment score.

What still remains to be determined is how big the sampled blocks must be and how this decision affects the choice of tile sizes that GACT can use afterwards to guarantee it can handle all (smaller) vertical/horizontal jumps that may be encountered during its alignments.

6 Conclusion

This paper analyzed the limitations of Darwin’s design for protein sequence alignment, a more difficult problem due to the increased search space and hardware requirements. It proposes a multi-pass aligner variant that helps increase the aligner accuracy of Darwin on protein sequences. Our proposed multi-pass variant of GACT achieves on average 14% better alignment scores. We also discuss a refinement to Darwin’s filtering algorithm that can be used in conjunction with our multi-pass algorithm to bridge a gap that is currently yet unhandled by both Darwin’s aligner and our multi-pass variant: obtaining good alignments under the presence of large insertions and deletions.

Acknowledgements

We thank Yatish Turakhia for sharing the code for his D-SOFT filter and GACT aligner upon which this work is based, as well as for his valuable assistance.