Parallel Smith-Waterman with OpenMP

Isitha Subasinghe University of MelbourneMelbourneVictoriaAustralia isubasinghe@student.unimelb.edu.au

Abstract.

Sequence alignment is the process of comparing and detecting similarities in biological sequences such as RNA and DNA. The idiomatic algorithm presented for this problem is inherently sequential. This report explores how the problem may be parallelised using the ”wavefront” technique along with OpenMP ¹1See https://www.openmp.org . The algorithm implemented using OpenMP was able to offer a speedup of 13 over the sequential algorithm, interestingly, the effeciency of adding more cores to the problem was shown to decrease exponentially.

1. Background

The Smith-Waterman [10] algorithm introduced about 40 years ago, is the only known algorithm that is able to find the optimal local alignment of sequences, however, it is also one of the slowest algorithms for sequence alignment. With the exponential growth of biological databases, far outpacing microprocessor improvements, the use of parallelism to speed up the Smith-Waterman algorithm is of importance to bioinformatics. The algorithm itself is not trivially parallelisable, however, recognising the dependencies for the values in the dynamic programming matrix helps in observing that the algorithm is still however, parallelisable. The dependencies of the values in the dynamic programming matrix are shown in fig. 1.

Figure 1. Dependencies for values in Smith-Waterman matrix

2. Related Work

Searching of DNA and RNA sequences is a fundamental task in bioinformatics and as mentioned previously Smith-Waterman offers optimal local alignment. However, as a result of optimal local alignment being computationally expensive, the most common solutions do not aim to provide this, instead they use heuristics to reduce the computational cost, examples of such tools are FASTA [8] and BLAST [1].

Early work in parallelising Smith-Waterman was based upon diagonal parallelisation, this was initially done through SIMD instructions.SIMD implementations have proven to be popular as a result of their common availability in x86 systems and have demonstrated significant speedups in several published literature [9, 12].

With the rapid accessibility of GPUs over the years following the early work on SIMD parallelism, CUDA compatible GPUs were found to be effective as hardware accelerators for the local optimal alignment problem. The CUDA implementation shown in [5] was demonstrated to be 2 to 30 times faster than other implementations on commodity hardware at the time.

Field Programmable Gate Arrays are integrated circuits that can be programmed through the user, using some hardware description language such as VHDL or Verilog. These implementations are naturally faster than their counterparts in higher level languages as the circuit itself is programmed for the specific problem, however, they are difficult to program and debug. Storaasli [11] claims that their implementation of Smith-Waterman on two Cray XD1 systems was able to obtain a 100x speedup. Interestingly, speedups of 160x have been made through the use of FPGAs as shown in [4].

3. Proposed Approach

3.1. Sequential Algorithm

The sequential algorithm used was near identical to the provided program, however certain optimisations were performed to the program, along with bug fixes. The code for this algorithm is shown below:

Algorithm 1 Sequential Smith-Waterman

1:procedure SeqSmithWaterman

\text{m,n}\leftarrow\text{dimensions of }\textit{matrix}

pgap\leftarrow\textit{gap penalty}

pxy\leftarrow\textit{mismatch penalty}

\text{dp[m][n]}\leftarrow\text{newMatrix(m+1,n+1)}

\text{x}\leftarrow\textit{sequence x}

\text{y}\leftarrow\textit{sequence y}

8: for

\text{i}\leftarrow[0\ldots m]

\text{dp[i][0]}\leftarrow\text{i}*\text{pgap}

10: for

\text{i}\leftarrow[0\ldots n]

11:

\text{dp[0][i]}\leftarrow\text{i}*\text{pgap}

12: for

\text{i}\leftarrow[1\ldots m]

13: for

\text{j}\leftarrow[1\ldots n]

14: if

\text{x[i-1]}=\text{y[j-1]}

then

15:

\text{dp[i][j]}=\text{dp[i-1][j-1]}

16: else

17:

\text{dp[i][j]}=\text{min(dp[j-1][i-1] + pxy, dp[j-1][i] + pgap, dp[j][i-1] + % pgap)}

return traceback(dp)

3.2. Parallel Algorithm

The psuedocode for the parallel algorithm is defined as follows:

Algorithm 2 Parallel Smith-Waterman

1:procedure ParSmithWaterman

\text{m,n}\leftarrow\text{dimensions of }\textit{matrix}

pgap\leftarrow\textit{gap penalty}

pxy\leftarrow\textit{mismatch penalty}

\delta_{i}\leftarrow\textit{tile width}

\delta_{j}\leftarrow\textit{tile height}

outer_{i}\leftarrow 0

outer_{j}\leftarrow 0

diagonals\leftarrow\left\lfloor\frac{m}{\delta_{i}}\right\rfloor+\left\lfloor% \frac{n}{\delta_{j}}\right\rfloor+\begin{cases}1&m\%\delta_{i}+n\%\delta_{j}>0% \\ 0&m\%\delta_{i}+n\%\delta_{j}=0\\ \end{cases}

10:

\text{dp[m][n]}\leftarrow\text{NewMatrix(m+1,n+1)}

11:

\text{x}\leftarrow\textit{sequence x}

12:

\text{y}\leftarrow\textit{sequence y}

13: for

\text{i}\leftarrow[0\ldots m]\textbf{ parallel}

14:

\text{dp[i][0]}\leftarrow\text{i}*\text{pgap}

15: for

\text{i}\leftarrow[0\ldots n]\textbf{ parallel}

16:

\text{dp[0][i]}\leftarrow\text{i}*\text{pgap}

17: for

\text{i}\leftarrow\{0,\ldots,\text{diagonals}-1\}

18:

diff_{i}\leftarrow outer_{i}

19:

diff_{j}\leftarrow n-outer_{j}-1

20:

diag_{i}\leftarrow 1+\left\lfloor\frac{diff_{i}}{\delta_{i}}\right\rfloor

21:

diag_{j}\leftarrow 1+\left\lfloor\frac{diff_{j}}{\delta_{j}}\right\rfloor

22:

length\leftarrow min(diag_{i},diag_{j})

23: for

\text{tile}\leftarrow\{0,\ldots,\text{length}-1\}\textbf{ parallel}

24:

inner_{i}\leftarrow max(1,outer_{i}-(tile*\delta_{i}))

25:

inner_{j}\leftarrow max(1,outer_{j}+(tile*\delta_{j}))

26:

imax\leftarrow max(inner_{i}+\delta_{i},width)

27:

jmax\leftarrow max(inner_{j}+\delta_{j},height)

28: for

\text{i}\leftarrow[inner_{i}\ldots imax]

29: for

\text{j}\leftarrow[inner_{j}\ldots jmax]

30: if

\text{x[i-1]}=\text{y[j-1]}

then

31:

\text{dp[i][j]}=\text{dp[i-1][j-1]}

32: else

33:

\text{dp[i][j]}=\text{min(dp[j-1][i-1] + pxy, dp[j-1][i] + pgap, dp[j][i-1] + % pgap)}

34: if

outer_{i}+\delta_{i}<m

then

35:

outer_{i}\leftarrow outer_{i}+\delta_{i}

36: else

37:

outer_{j}\leftarrow outer_{j}+\delta_{j}

return traceback(dp)

3.2.1. Restrictions

There was a restriction in the implementation of our algorithm, namely the amount of memory available proved constraining in exploiting further opportunities for parrallelism. Computing the sequence alignment among all $\frac{k(k-1)}{2}$ pairs is “embarrassingly parallel”, however, this leads to far greater memory consumption and therefore wasn’t a valid optimisation given our memory restrictions.

3.2.2. NUMA awareness

NUMA (non-uniform memory access) becomes an area of concern when considering programs that use multiple cores. More accurately, when multi-core programs require access to the same memory locations as each other and this sharing is performed over an interconnect. The policy used to help reduce the cost of NUMA in the Parallel implementation is called “First Touch”. “First Touch” ensures that the NUMA node to initialise some memory location, will have that page allocated on a memory device physically closer to the NUMA node.
The steps for memory allocation can be given as follows:

•

dynamic memory allocation call (malloc/new)
•

initial write attempt of dynamically allocated memory
•

page fault triggered when write is attempted
•

physical allocation of memory performed and mapped via MMU

The parallel implementation exploits this property on moderns operating systems to help reduce the cost of non-uniform memory access. This is performed by initialising the matrix used for Smith-Waterman through multiple threads, ensuring that each NUMA node has a slice of the memory and thus a lower number of memory transfers are done over the interconnect, this is shown in algorithm 3.

While this is helpful for distributing the pages among the NUMA nodes, it comes at an additional processing cost. With the GNU 11.1.0 compiler toolchain, we are able to make use of the OpenMP allocator to do this distribution of pages, additionally some further optimisations maybe performed to increase performance, such as removing thread safety from the allocator and ensuring that the allocated memory is aligned, aligning memory is useful when using SIMD as load operations are much faster.

3.2.3. Optimum tile size

The tile size presents a significant influence on the speedup of our algorithm, along with the scheduling algorithm used. The intuition regarding the dimensions of the tile size maybe to minimise it as much as possible, however, this is not always the optimial strategy. Likewise, it may seem that static scheduling maybe the most sensible given that each core would perform relatively the same amount of work, however, this is also untrue. A deeper discussion regarding the scheduling of wavefront parallelism on shared memory multiprocessors is provided in [6]. In our implementation, the tile size was determined by dividing the dimensions by a fixed number and the scheduling was performed dynamically, these details were determined through experimentation.

3.2.4. Cache optimisations

Memory latencies play a significant role in the performance of programs, ensuring false sharing is minimal and maintaining locality of reference are important objectives for any high performance application. The processing of tiles in the wavefront processing of the matrix, helps obtain the previously mentioned objectives as shown by [6]. This is as a result of the layout of memory used in the matrix itself, the entire matrix is allocated as a single contiguous slice of memory. Consider the parallel traversal of the yellow diagonals below in fig. 2. An isolated view of the potential cache lines for each core is also visible in fig. 3, observe how given sufficient width of the matrix, false sharing is extremely unlikely. Additionally, note how the tiles are traversed horizontally and then vertically, this maintains locality of reference as the rows are contiguous in memory which is crucial to high performant code.

Parallel processing across matrix diagonals — Figure 2. Parallel processing of diagonals

Potential cache-line layout across matrix tiles — Figure 3. Potential cache lines

3.2.5. Compiler optimisations

Compiler optimisations were heavily used to further increase the performance of the algorithm and frequent profiling of the program was done via the linux utility “perf“ ²2See https://perf.wiki.kernel.org/index.php/Main_Page. It is important to note that using “std::string” had a significant impact on performance when dynamically linked to our program, this is as a result of symbol lookup performed in a tight loop when traversing the dynamic programming matrix. The remedy was to statically link the C++ standard library, which lead to a significant increase in performance, almost exactly the same performance characteristics as using a simple character array for holding the sequences. Additional optimisation recommendations from [2] were evaluated and used where appropriate. A critical optimisation used was the “-march=native” flag passed onto the GNU compiler, this flag targets the exact architecture that we are compiling the program for, this results in further architecture level optimisations being possible such as usage of AVX-512 ³3See https://www.intel.com.au/content/www/au/en/architecture-and-technology/avx-512-overview.html.

4. Experiments

4.1. Experiment Environment

Two experimentation environments were used throughout the development process. The development environment was used as a quick way to test and filter out ideas that weren’t viable, while valid ideas were tested on the Spartan High Performance Computing cloud [7].

4.1.1. Development Environment

•

CPU - AMD Ryzen™ 7 3700X @ 3.6-4.4GHz
•

Memory - 32GB
•

Compiler - GNU 11.1.0

4.1.2. Spartan Environment

•

CPU - Intel(R) Xeon(R) Gold 6154 CPU @ 3.00-3.7GHz
•

Memory - 32GB
•

Compiler - GNU 10.3.0

4.2. Results

The results after running our parallel algorithm on the “mseq-big13-example” data provided are shown in fig. 4 and presented in below in table 1.

Smith-Waterman runtime results by core count — Figure 4. Results observed

Core(s)	Walltime(microseconds)
1	923129129 $\mu s$
2	491227097 $\mu s$
8	169576852 $\mu s$
16	119517599 $\mu s$
24	112753964 $\mu s$
32	111237046 $\mu s$
36	108483498 $\mu s$

Table 1. Walltime per number of cores

4.3. Determining Efficiency

The efficiency( $\eta(n)$ ) of our algorithm can be determined by the following formula:

(1)

\eta(n)=\frac{\text{Speedup(n)}}{\text{p(n)}}

The efficiency for the results shown in table 1 have been calculated using eq. 1 and displayed in section 4.3.

Core(s)	Efficiency
1	1.572
2	1.478
8	1.070
16	0.759
24	0.537
32	0.408
36	0.372

4.4. Exploring implications of tile size

Core(s)	Tile Size	Walltime(microseconds)
8	$(\frac{m}{8}+1,\frac{n}{8}+1)$	379606318 $\mu s$
16	$(\frac{m}{16}+1,\frac{n}{16}+1)$	382288384 $\mu s$
24	$(\frac{m}{24}+1,\frac{n}{24}+1)$	388430562 $\mu s$
32	$(\frac{m}{32}+1,\frac{n}{32}+1)$	385640454 $\mu s$
36	$(\frac{m}{36}+1,\frac{n}{36}+1)$	381233563 $\mu s$

Table 2. Walltime per number of cores and tile size

4.5. Discussion

Note: The sequential implementation was measured to have taken 1451965419 $\mu s$ .
The results were extremely surprising, note that the single core parallel algorithm is far more performant than the sequential algorithm. While this was unexpected, this is likely due to the substantial amount of micro optimisations and time spent on the implementation of the parallel algorithm itself rather than the complexity properties of the parallel algorithm itself.
The choice of cores to test were not random, they were also made with the intention of determining the effectiveness of our NUMA mitigation strategy mentioned previously. The processor used consisted of 18 cores, this fact ensures that when the core count requested is above 18, we are guaranteed to encounter at least two NUMA nodes. Note that, in practice, this number was in fact much higher at 3-5 NUMA nodes per experiment. The results in table 1 show that our NUMA mitigation strategy was effective, although NUMA certainly would have contributed in adding latencies to our calculations as our core count went up, the increase in parallelism was able to counteract this. The decrease in efficiency as core count increased, demonstrated in section 4.3 is likely a function of NUMA costs and the limit on the parallelisability of the problem itself.

As stated previously, the tile size has a significant impact on the parallelisability of the problem. As shown in table 2, choosing a bad tile size may mean the parallelisability of the problem is artificially limited. The results in table 1, fig. 4 and section 4.3 are a result of dynamically choosing the tile size based on the dimensions of the matrix. To be exact, the dimensions of the matrix were divided by 13, determining this number was through experimentation on the development environment.

5. Further Work

5.1. Generation of tiles

The parallel algorithm implemented does not make full use of the available computing power. When tiles are iterated over, they must wait until all the tiles in the same diagonal have completed their work. This is an unnecessary restriction, as it results in a large amount of time spent on the barrier. Further work could explore how to reduce the level of granularity in our dependencies, new tiles could simply be generated after a single tiles work unit has completed. A simple algorithm would be to start processing tile $(i+\delta_{i},j)$ and $(i,j+\delta_{j})$ when $(i,j)$ has completed its unit of work, this does however introduce the problem of tiles being visited more than once but this could simply be avoided through using a xorfilter to determine if a pair of indices have been visited before.

5.2. Inter-sequence SIMD parallelism

We could apply inter-sequence SIMD parallelism into our current thread-parallel implementation to obtain instruction level parallelism. While this is plausible, an important consideration that warrants investigation is determining an efficient method to move data from the heap and into the SIMD register, this technique will likely have significant speedups with AVX-512.

6. Conclusion

We have demonstrated how OpenMP may be used to develop a parallel algorithm for optimal local alignment of sequences. The algorithm developed showed effective use of parallelism to achieve a significant speedup of 13 over the sequential algorithm.

References

S. F. Altschul, T. L. Madden, A. A. Schäffer, J. Zhang, Z. Zhang, W. Miller, and D. J. Lipman (1997) Gapped blast and psi-blast: a new generation of protein database search programs. Nucleic acids research 25 (17), pp. 3389–3402. Cited by: §2.
A. Fog (2020) Optimizing software in c++. Self-publishing. Cited by: §3.2.5.
C. Lameter (2013) An overview of non-uniform memory access. Communications of the ACM 56 (9), pp. 59–54. Cited by: Appendix A.
I. T. Li, W. Shum, and K. Truong (2007) 160-fold acceleration of the smith-waterman algorithm using a field programmable gate array (fpga). BMC bioinformatics 8 (1), pp. 1–7. Cited by: §2.
S. A. Manavski and G. Valle (2008) CUDA compatible gpu cards as efficient hardware accelerators for smith-waterman sequence alignment. BMC bioinformatics 9 (2), pp. 1–9. Cited by: §2.
N. Manjikian and T. S. Abdelrahman (1996) Scheduling of wavefront parallelism on scalable shared-memory multiprocessors. In Proceedings of the 1996 ICPP Workshop on Challenges for Parallel Processing, Vol. 3, pp. 122–131. Cited by: §3.2.3, §3.2.4.
B. MEADE, L. LAFAYETTE, G. Sauter, and D. TOSELLO (2017) Spartan hpc-cloud hybrid: delivering performance and flexibility. University of Melbourne. External Links: Link, Document Cited by: §4.1.
W. R. Pearson and D. J. Lipman (1988) Improved tools for biological sequence comparison. Proceedings of the National Academy of Sciences 85 (8), pp. 2444–2448. Cited by: §2.
T. Rognes and E. Seeberg (2000) Six-fold speed-up of smith–waterman sequence database searches using parallel processing on common microprocessors. Bioinformatics 16 (8), pp. 699–706. Cited by: §2.
T. F. Smith, M. S. Waterman, et al. (1981) Identification of common molecular subsequences. Journal of molecular biology 147 (1), pp. 195–197. Cited by: §1.
O. O. Storaasli and D. Strenski (2007) Accelerating genome sequencing 100x with fpgas. Technical report Oak Ridge National Lab.(ORNL), Oak Ridge, TN (United States). Cited by: §2.
M. Zhao, W. Lee, E. P. Garrison, and G. T. Marth (2013) SSW library: an simd smith-waterman c/c++ library for use in genomic applications. PloS one 8 (12), pp. e82138. Cited by: §2.

Appendix A Appendix

Algorithm 3 Parallel Matrix initialisation

procedure NewMatrix

\text{m,n}\leftarrow\text{dimensions of }\textit{matrix}

\text{dp}\leftarrow\text{malloc(sizeof(int *)*m)}

\text{size}\leftarrow\textit{m*n}

\text{dp0}\leftarrow\text{malloc(sizeof(int)*size)}

6: for

\text{i}\leftarrow\{0,\ldots,\text{size}-1\}\textbf{ parallel}

\text{dp0[0]}\leftarrow 0

for

\text{i}\leftarrow\{0,\ldots,m-1\}\textbf{ parallel}

\text{dp[0]}\leftarrow 0

\text{dp[0]}\leftarrow\text{dp0}

for

\text{i}\leftarrow\{1,\ldots,m-1\}

12:

\text{dp[0]}\leftarrow\text{dp[i-1] + n}

return dp