Parallel Computing 87 (2019) 11–24 Contents lists available at ScienceDirect Parallel Computing journal homepage: www.elsevier.com/locate/parco Exploring efficient data parallelism for genome read mapping on multicore and manycore architecturesR Shaolong Chen∗, Miquel Angel Senar Universidad Autonoma de Barcelona, Computer Architecture and Operating System, Barcelona, Spain a r t i c l e i n f o a b s t r a c t Article history: Nowadays heterogeneous architectures formed by multicore and manycore systems have become attrac- Received 1 September 2018 tive solutions to cope with the data booming in genomic-based studies. Our work explores the efficient Revised 28 March 2019 usage of heterogeneous architectures in such area. In particular, we have studied the use of manycore Accepted 28 April 2019 components like the Xeon Phi accelerator, which has proved to be a convenient choice because it allows Available online 15 May 2019 an easy migration of applications developed for multicore servers based on the × 86 architecture. Our Keywords: study also focuses on the problem of sequence alignment, which is one of the fundamental and most Sequence alignment costly computational stages in most genome variant studies. We concentrate our attention on BWA, one Heterogeneous architecture of the most popular sequence aligners, and we have focused our attention on three types of heteroge- Intel xeon neous systems, one containing Intel multi-core CPUs and accelerators, one that are made up of several Intel xeon phi multi-core servers, and one large-scale system. Each with different characteristics in terms of number NUMA Node of CPUs, number of cores and system organization memory. Although the problem of alignment of se- quences fits in the embarrassingly parallel pattern, achieving good performance and good scalability in heterogeneous environments can be complex. We have analyzed different strategies based on the dis- tribution of data and the replication of certain data structures and we found that MDPR (Multi-level Data Parallelization and Replication) strategy has shown the best results in all the heterogeneous plat- forms tested. Its results have surpassed other strategies proposed in the literature and have shown its malleability to be used in different heterogeneous environments without the need to apply specific ad- justments according to the underlying architecture. In the design of MDPR, different static and dynamic data distribution strategies have also been evaluated. The best results were obtained by the static strat- egy, which has a significant preprocessing cost. However, the dynamic strategy of data distribution using a round-robin mechanism obtained similar times without the need for the preprocessing stage. Although our proposal was applied to BWA using human genome data samples, this strategy can be easily applied to other sequence datasets and alignment tools that have similar operating principles with those of BWA aligner. © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/) 1. Introduction coprocessor, which has became an attractive solution for existing programs on multicore systems that can be easily ported because In recent years, heterogeneous systems consisting of different Xeon Phi is also based on the × 86 architecture. architectures have become a significant trend in the high perfor- Compared to computing evolution (doubling the number of mance computing arena [1]. Manycore architectures, such as GPUs, transistors every 18 months), the amount of sequence genome data have obtained considerable popularity, despite the fact that the only doubles every 12 months [4]. However, this is an spectacular specific programming languages such as OpenCL and CUDA require growth of data that constitutes an important challenge that still significant coding skills [2]. Intel Xeon Phi [3] is another manycore requires efficient applications to process them. Sequence alignment, variant calling and variant annotation R comprises three fundamental operations in genome data studies Conflict of interest. None. ∗ Corresponding author. [5]. Sequence alignment is a crucial step that can provide pri- E-mail addresses: shaolong.chen@uab.cat (S. Chen), miquelangel.senar@uab.cat mary consequences for the other two operations. It involves map- (M.A. Senar). ping short reads to a genome sequence reference as accurately as https://doi.org/10.1016/j.parco.2019.04.014 0167-8191/© 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/) 12 S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 possible [6]. Although many aligners have been developed in re- such as Blast [10], Bowtie [11], SOAP2 [12] and BWA [13], but few cent years, such as MAQ, SOAP and BLAST, they all exhibit a signif- works have been applied to Intel Xeon Phi manycore systems or icant execution time and a large memory footprint. Although, its other heterogeneous systems. execution time is similar to other existing aligners, BWA (Burrows- Several researchers have studied the problem of read alignment Wheeler Aligner) [7] takes advantage of the BWT indexing tech- in multicore architectures with shared memory. Zhang et al. [14] il- nique that decreases memory requirements of the alignment pro- lustrated that sequence alignment is memory-sensitive and tried cess. This reduction in memory needs combined with its mapping to optimize the usage of BWT index by decreasing their irregu- accuracy have made BWA one of the most popular aligners in the lar memory accesses. Kathiresan et al. [15] presented several par- scientific community. allelization strategies, including memory cache improvements, that In this paper, we explore different strategies that allow BWA to could be useful on multicore system. Macedo et al. [16] introduced run on heterogeneous systems that combine multicore CPUs with a hybrid programming schema using OpenMP and MPI to facilitate manycore accelerators or different multicore systems. We analyze multiple sequence alignment application on heterogeneous multi- some proposed strategies in the literature and propose enhance- core clusters, but they do not explore their approach to sequence ments and new strategies to improve global performance. This alignment in manycore architectures. Olivier et al. [17] compared paper extends our previous studies in which we proposed some different task scheduling algorithms on multicore systems to eval- data parallelization/replication strategies and made an evaluation uate their impact in overall performance. Lenis et al. [18] inves- of BWA-based aligners in manycore-based systems. In those works tigated memory interleaving and data parallelization strategies for our study focused only on heterogeneous systems that incorporate NUMA multicore systems. Chen et al. [19] studied an offloading ac- accelerators. In one case, we analyzed improvements when the ap- celeration scheme applied to BWA aligner in FPGA architectures. plications were executed in native mode [8], and in [9] we studied Houtgast et al. [20,21] made an effort to adopt alignment appli- different variants to execute BWA in symmetric mode. The MDPR cations on GPU- and FPGA-based architectures. Several proposals strategy analyzed in the current paper derives from these previ- of alignment performance improvement on FPGA are discussed in ous works and it has been generalized and extended so that it can [22]. Additionally, many papers [23] studied sequence alignment be executed indistinctly not only in accelerator-based servers but improvements on manycore architecture based on GPUs. in other types of heterogeneous environments formed by multiple Particular proposals applied to BWA have also been produced CPUs organized around different NUMA architectures. The strategy in recent years. Many of them describe different approaches to MDPR utilizes multiple instances in order to reduce memory bot- improve the performance of BWA aligner on some particular ar- tlenecks by decreasing memory congestion and improving memory chitecture. Zhang et al. [14] found aligner BWA is a memory- locality. We present a study where we analyze the performance sensitive application by optimizing cache mechanisms on multi- of different BWA variants on multicore and manycore system. We core architectures. Houtgast et al. [24] accelerated BWA by offload- have augmented also MDPR with different mechanisms of dynamic ing some highly parallel kernels on a GPU system. Houtgast et al. data distribution and we have evaluated them with two represen- [25] utilized a simple adaptive load balance algorithm on GPUs and tative datasets and three heterogeneous architectures. The main achieved a considerable performance improvement. Other existing contributions of this paper can be summarized as follows: approaches are using Hadoop technology, such as BigBWA [26] that achieved a significant improvement of performance on multiple • We propose a strategy (named MDPR) that has been applied nodes. Or SparkBWA [27] that tried to use another big data tech- to BWA by adding new mechanisms of data parallelism and nology, Spark, to boost BWA aligner. MICA [28] explored the extra data replication. The strategy exhibits a hierarchical struc- parallelism available in each Xeon-Phi core. Cui et al. [6] proposed ture based on different instances of computation that in- mBWA which works on Xeon-Phi systems using the offload mode. clude replication of certain data structures and the distribu- mBWA utilizes a multi-level parallelization strategy that includes tion of input data among all parallel instances. data IO parallelization and a parallel pipeline in reads alignment. • We analyze three data distribution strategies for MDPR, one Tian et al. [29] presented several practical SIMD vectorization tech- static and two dynamic (even and round-robin). niques on Xeon-Phi architectures. However, such SIMD techniques • We have carried out an experimental study comparing dif- require complex code rewriting in order to incorporate them with ferent versions of BWA on multicore and manycore sys- the existing genome aligners. Memeti et al. [30] and Wang et al. tems, including the strategies proposed in this work and [31] provided solutions for large scale DNA analysis which exploits other alternatives proposed in the literature. This study has thread affinity level and SIMD techniques in the Xeon Phi, but their shown us the best alternatives to run BWA on manycore sys- solution does not deal with data parallelism among threads or in- tems based on Intel Xeon Phi and on heterogeneous systems stances. pBWA [32] is an efficient parallel version of BWA based based on Intel and AMD multicore processors. on the Open MPI library. It not only preserves the multi-thread The rest of the paper is organized as follows. Section 2 de- capability provided by BWA but also adds efficient parallelization scribes related work in performance improvement of sequence for core alignment functions on the heterogeneous cluster. Unfor- alignment algorithms, including BWA, on parallel architectures. tunately, it uses Pthread library and cannot be evaluated for a large Some general background in multicore and manycore architectures number of instances of pBWA on the Xeon Phi architecture due is presented in Section 3 and BWA aligner is described in Section 4. to its memory limitations. Herzeel et al. [33] evaluated an Intel Section 5 discusses DP, DR, DPR and MDPR strategies. Section 6 de- Cilk implementation of BWA on heterogeneous systems that cor- scribes our experimental environment and Section 7 shows the re- rected data load imbalance produced by differences in data reads. sults we obtained in the experiments. Finally, Section 8 discusses We compared this BWA Cilk version with the original BWA Pthread the main conclusions of this research and highlights some future version and an BWA OpenMP-based version proposed in our pre- work. vious work [9]. 2. Related work 3. Multicore-manycore architectures Our work explores effective data parallelism solutions for the Nowadays, systems based on NUMA (Non Uniform Memory Ac- BWA aligner that are suitable for heterogeneous architectures. cess) architectures are the most used in the field of high perfor- Many sequence aligners have been studied on multicore system, mance computing. These systems are often used in conjunction S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 13 Fig. 1. Heterogeneous multicore-manycore architecture. with manycore accelerators that increase the number of cores The Intel Xeon Phi has an × 86 architecture that allows three available in the system at a reasonable cost. In this section we execution modes. An application can be executed in three different provide a brief description of the basic architectures that we have ways when submitted to a Xeon Phi system [35]: used in our study, namely, multicore systems with NUMA nodes • Native mode: the application is executed on each compo- and manycore accelerators. nent (processor or accelerator) independently. Fig. 1 shows the main architectures of the systems used in this • Symmetric mode: the accelerator is seen as a regular node work. Two multicore systems are depicted at the left side and cen- in the system and the application is executed on all compo- ter of the figure (an AMD-Opteron based server and an Intel Xeon- nents (main processor and an accelerator). based server). Intel Xeon Phi is at the right side. Our AMD Opteron • Offload mode: the application is executed on the main pro- and Intel Xeon systems have a NUMA node with different link lev- cessor but selected highly parallel sections pass to the accel- els across processors. erator. AMD-Opteron and Intel Xeon systems have both a NUMA archi- tecture with four independent sockets, each socket has a one or 4. Genome aligner: BWA two multicore processors and each processor is attached to a local memory bank. Each processor and its corresponding memory bank Genome read aligners provide the relative position of short is referred as a NUMA node. One of the most significant charac- reads within a reference genome. Despite the particular differences teristics of NUMA systems is that accesses from threads running of each aligner, they all share a similar mode of operation that can in local cores (located in the same NUMA node) have lower la- be summarized as follows: tency than accesses that go to physical memory located in different NUMA nodes (as shown in Fig. 2). This asymmetry in memory ac- 1. There is a set of reads that can be mapped to the reference cesses involves additional complexity when executing parallel ap- genome independently. plications because access to remote banks increases the execution 2. There is a reference genome data structure (genome index) time and, in the case of genome alignment, introduces also con- that is read-only data and is used to map each individual gestion problems. read. Intel Xeon Phi accelerator (seen at right side of Fig. 1) has a 3. The results consist in populating a shared data structure that bi-ring architecture with 57 to 61 cores inside. Each core in the would be written on an output file at the end of the align- Xeon Phi has a SPU, a VPU with 512-bit SIMD capability, and a pri- ment. vate 512KB L2 cache that is maintained fully coherent by a global- Read aligners need to construct an index for short reads or for distributed tag directory (TD) [34], which forms a unified shared L2 genome reference or for both. Depending on the index construc- cache of 30.5MB. Xeon Phi can achieve 1 TFlops in double precision tion, read aligners can be classified into three groups [23]: hash or 2 TFlops in single precision. In addition, memory bandwidth can table-based, FM-index based and merge sorting based. The former reach 352GB per second theoretically thanks to 16 memory chan- two groups contain the majority of aligners used by the scientific nels. The memory controllers (GDDR MC) and the PCIe Client Logic community. BLAST, MAQ, SOAP and ZOOM are examples of hash provide a direct interface to the GDDR5 memory and the PCI ex- table-based aligners. Their main drawback is the large memory press bus, respectively. consumption required for the index. In contrast, the memory re- quirements of aligners based on the FM-index mechanism are sig- nificantly lower, which explains their increasing popularity. We focus our study on BWA, written by Li [7], which is one of the most popular sequence alignment tools from the FM-index based family. BWA consists of three algorithms, BWA-backtrack, BWA-SW and BWA-MEM. Our study focuses on BWA-backtrack that has already been used by other works in the literature, such as mBWA [6], pBWA [32] and BMIC [36]. Thus, we will try to pro- vide a view as complete as possible of the performance that each of the different strategies can obtain when applied to the same ba- sic application. BWA creates the index of genome reference beforehand. Later, Fig. 2. NUMA memory architecture. short read data is processed using fixed sized chunks, in a round- 14 S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 the system and, 2) the memory congestion exhibited by the bank where the reference genome is allocated. Moreover, the library and the organization used to handle ap- plication threads can cause certain variations in the execution times obtained. In our paper [9], it was observed that the par- allelization using Pthread, OpenMP and Cilk presented noticeable differences in performance, with OpenMP and Cilk being the best alternatives to Pthread, which was the library originally used in BWA. Besides the thread library, the main performance problems in genome read aligners arise from the use of memory in NUMA systems. Memory allocation is done according to a first-touch pol- icy and this means that the genome reference index is allocated in a single bank, that will become a bottleneck when multiple threads try to access it [18]. Fig. 3. Scalability on 2-socket 2-NUMA 24-thread system. 5. Strategies based on data parallelization and replication The scalability problems of BWA are mainly due to the con- tention generated by the need for multiple threads to access the same memory bank where the reference genome is located. To mitigate this problem we propose a series of alternatives based on a basic common scheme. In this scheme, there is a first level of parallelism based on groups of threads that handle groups of reads. On the other hand, the reference genome will be replicated in different instances, each of which will be accessible by different thread groups. The number of groups of threads and the number of replicas of the reference genome can vary and can be organized in different ways. The final objective is to analyze which variant ob- tains the best results in multicore or manycore architectures. In the manycore systems we are interested in finding strategies that are Fig. 4. Scalability on 4-socket 4-NUMA 64-thread system. efficient using all the resources of the system (that is, executing in symmetric mode). And we are also interested in solutions for manycore systems that can be easily transferred to multicore sys- robin pattern among threads. Each thread in BWA processes one tems composed of servers with different characteristics in terms of chunk of data, 256 K by default, that is independently mapped to performance. The variants that have been implemented and tested the genome reference. Like many other aligners, BWA has multi- are the ones described below. threading capability (using the Pthread library) to solve read map- ping operations in parallel since there is no data dependency in 5.1. Data management on homogeneous architectures: DP, DR and this action. The addition of multi-threading to an aligner does not DPR strategies guarantee that processing resources are used efficiently. In fact, it is quite common to observe that efficiency decreases when the We are exploring data management that involves both paral- thread count grows large enough [8]. lelism and replication. As shown in Fig. 5, two datasets are in- This phenomenon is shown in Figs. 3 and 4. These figures illus- volved in the alignment process. N corresponds to the set of all trate the results obtained by three sequence aligners (BWA, SOAP2 short reads, and n is the number of read subsets in which N is di- and BOWTIE2) that are all FM-index-based application when map- vided. M represents the reference genome, and m is the total num- ping two large short read genome examples of 17.6GB against hu- ber of replicas made of such genome. Therefore, read subsets are man genome. Each figure illustrates the results obtained when referred as N1, .., Ni, . . . , Nn and reference replicas are referred as aligners are executed on two different servers. Both servers are M1, ..M j, . . . Mm according to Fig. 5. based on NUMA architecture but they differ in the total number of processors, cores and NUMA nodes available. The first system (Fig. 3) has 2-sockets with 2-NUMA Intel Xeon CPU with 12 cores each). Twenty-four threads can be executed simultaneously using Hyper Thread technology. The second system (Fig. 4) is a 4-socket 4-NUMA Intel Xeon CPU. Each socket contains 8-core processor, so the total number of cores is 32 and 64 threads can be executed si- multaneously. In both cases, each socket is connected to a memory bank, which constitutes a NUMA node. NUMA nodes are connected by QuickPath Interconnect link. The performance degradation experimented by all three strate- gies can be observed in Fig. 3 as the parallelization degree in- creases. Execution times are reduced until half of the available cores in each architecture are used. From that point, the execu- tion times hardly improve or even get worse slightly. This abnor- mal phenomena in performance appears when memory becomes the application bottleneck and it could be explained as follows, 1) the increasing number of accesses to remote memory bancks in Fig. 5. Sequence alignment procedure. S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 15 Fig. 6. Implementation of DP, DR and DPR. Although the reference genome could also be divided, and thus 3. Collecting all alignment sub-results and merging them into  increasing the degree of parallelism available in the application final output R = Ri j; the time of this step time is T(R); and reducing memory consumption, this approach has been dis- i=n, j=m Total execution time is T = T (M ) + T (R ) + i=1, j=1 T (Ri j ). Ac- carded in our study for three main reasons. First, the index of each part of the reference genome would have to be reconstructed at tually, the mapping consumes most of the execution time because run-time, increasing the total execution time (in the traditional ap- loading the index and collecting subresults mainly involve linear operations of reading and writing files. Hence, we can simplify the proach, the index is built only once offline and that index is used i=n, j=m in all the executions of the aligner). Second, the final results of expression to T = i=1, j=1 T (Ri j ), being this stage the one that each subset (or blocks) have to be recombined to be assigned the can benefit from parallelization. correct position within the whole reference genome. Finally, a cer- There are different variations that can be applied to the values tain loss of precision in the final results is experienced due to reads for n and m. that might map in the cuts between successive blocks of the refer- 1. Original execution: n = 1, m = 1. ence genome. This case corresponds to the original structure of BWA. As In our approach, BWA is organized as a set of thread groups shown in Fig. 6a, multiple OpenMP threads are processing (Ng) that are responsible for mapping some reads. In practice, a all short reads using a single reference genome. thread group represents a group of threads in OpenMP or a pro- 2. DP strategy: n > 1, m = 1. cess in MPI, as shown in Fig. 6. Each thread group is executed In this case (see Fig. 6b), short reads are evenly divided into using a different number of threads depending on the hardware to n parts based on the number of thread groups that are cores available in the system. One thread group in the alignment used (Ng), and a single copy of the reference genome is is logically organized as a collection of threads that work on in- used by all thread groups. We refer to this strategy as data dependent datasets of short reads and share the genome index parallelization (DP). Alignment process runs in every thread among all the threads. For a total number of potential threads Nt group with the shared genome index and its relative subset and a number of thread groups Ng, each thread group consists of of short read. Subset results are finally gathered and merged Nt/Ng threads. For example, Nt is 24 and 240 in the cases of In- from all thread groups. tel Xeon and Xeon Phi, respectively. When applying two thread 3. DR strategy: n = 1, m > 1. groups (Ng = 2), 12(24/2) and 120(240/2) threads are included in In this case, the reference genome is replicated m times, each thread group on the Intel Xeon and on the Xeon Phi, respec- while short reads are shared among all thread groups. We tively. refer to this strategy as data replication (DR). As shown in The process of alignment has three main phases: Fig. 6c, the reference genome is replicated for every thread group in memory. The alignment process in each thread 1. Loading the index for the genome reference Mj; the time of group accesses to the shared set of short reads. Finally, the this phase is T (M j ) = T (M ); results are merged as in the previous case. 2. Mapping short read subset Ni to genome reference Mj, which 4. DPR strategy: n > 1, m > 1. achieves mapping sub-result Ri j = Map(Ni, M j ); the time of In this case, both short reads and reference genome have this step is T(Rij); several parts or copies (the set of reads is divided in n parts 16 S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 Fig. 7. MDPR strategy. and reference genome is replicated m times, being n and m The short read dataset (S) is partitioned into several subsets: the same value). As shown in Fig. 6d, each group of threads one subset for each node in the system (Fig. 7 shows the case of has one reference genome replica and short read dataset is two nodes: node X(S ) and node Y(S )). We have several thread evenly divided into as many parts as thread groups. This groups or instances in node X (denoted as S 1, . . . , S i, . . . , S x, strategy is referred as DPR as it combines the previous two 1 <= i <= x) and several thread groups in node Y (denoted as strategies. S 1, . . . , S j, . . . , S y, 1 <= j <= y). Thread groups are generated by taking into account the memory requirements of the aligner (in particular, the size of the reference genome) and the amount of available space on the memory banks in each system. Reference 5.2. Data management on heterogeneous architecture: MDPR strategy genome replicas are denoted as R 1, . . . , R x, and R 1, . . . , R y. Each BWA instance executes several threads to align short reads as- The strategies described above can be applied in homogeneous signed to it. It is worth mentioning that the number of instances multicore systems with shared memory. In the case of systems on a heterogeneous architecture does not require to be identical. In based on accelerators (executing in symmetric mode) or in the fact, in the experimental part we will show results obtained with case of independent multicore systems, we propose a new strat- various combinations that we tested. egy that generalizes DPR. It is referred as MDPR (Multi-level Data It is important to note that the division of the data set between Parallelization and Replication) and is illustrated in Fig. 7. The de- all the nodes of the system does not imply a significant increase sign of MDPR involves several challenges related to the use of an in communications because each group of threads is only provided environment based on heterogeneous and distributed computing with the appropriate pointers to the file where all the reads are resources. located. That is, the distribution of data is not done by a master On the one hand, a work distribution mechanism must be es- process that reads the data file and distributes the blocks to the tablished using message passing primitives that do not introduce other processes by using message passing mechanisms. In our case, excessive overhead in terms of communication and synchroniza- the master process only communicates the start and end points tion. MDPR minimizes communication and synchronization needs of each block to each process; afterwards, each process will read because there is no central process in charge of reading and dis- directly that part of the file that is accessible through the shared tributing the input data, and collecting output results. Each par- file system that is attached to all the computing nodes (Lustre, in ticipant process in the computation reads its data directly and our case). only needs to be informed of the file positions where this data According to this schema, the best performance should be ob- is located. Additionally, results are produced independently by tained if all thread groups consume a similar time. However, the each process and they are finally combined by the central process heterogeneity of the whole system may lead to different execution through a simple append operation. times if short reads are distributed evenly. This load imbalance sit- On the other hand, the configuration of the application and the uation can be ameliorated by partitioning the short read dataset distribution of the work must be done in such a way that a simi- in an uneven way. This partitioning can be done in different ways lar execution time is guaranteed in each of the nodes that make and we have evaluated three possibilities: a static distribution that up the system. MDPR has been designed in such a way that it divides the dataset off-line and two dynamic distribution mecha- can be dynamically configured according to the characteristics of nisms that divide the dataset at run-time. Details of each mecha- the NUMA architecture (adjusting the number of instances used) nism follows: and incorporates data distribution mechanisms (both static and dy- namic) that allow the distribution of data in an unbalanced way in 1. Static distribution. order to adapt the whole application to the relative performance In this mode, size of datasets S and S are computed be- exhibited by each node of the system. These mechanisms and their forehand and the initial dataset of short read is divided stat- subsequent evaluation are described below. ically at the start of execution. Each node reads then its S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 17 process a block of size (size(S )/(x + y )). The ratio between node X and node Y could set up as x/y. In node X we have size(S ) = size(S ) ∗ x/(x + y ), and in node Y we have size(S ) = size(S ) ∗ y/(x + y ). Similarly to what happens in the static distribution case, the original dataset is divided at the beginning of the execution and no significant overhead is introduced. However, this scheme does not lead in general to a good load balance if relative performance of nodes is significantly different. As shown in the central part of Fig. 8, master instance generates equal size divisions of the original dataset and each instance (slave) will process it. 3. Round-robin distribution. This distribution mechanism implies that a master instance distributes chunks of reads to slave instances when they have finished previous chunks. This mechanism should guar- antee a better load balance than the previous one at the ex- pense of a larger overhead at run time. Originally, each BWA thread processes a set of reads of the same size (256 K, by default). That size has shown a good trade-off between the computation time for mapping all reads in the set and the I/O time consumed with that set. In our implementation, this set size has been maintained and, therefore, chunck size is calculated by multiplying 256 K by the number of threads that are executed in each instance. For instance, if we execute BWA aligner in a system that has 64 available Fig. 8. Data distribution mechanisms. threads, chunk size should be at least 16 M (256 K∗ 64) if one instance is applied or 8 M (256 K∗ 32) if two instances block at the beginning without the need to perform new are applied. This mechanism is illustrated at bottom part in distribution operations later. The ratio between S and S is Fig. 8, where chunks of short read are distributed following computed according to the relative performance achieved by a round-robin way. An instance receives a new chunck when the best result achieved by DPR in each architecture. For in- it has completed the previous one. As mentioned previously, stance, in the case of two nodes (X and Y), for a certain chunk distribution simply implies the communication of the number of instances on the node X and on the node Y (x start and end points of the corresponding chunk. and y instances, respectively), the time-cost Tx and Ty are measured. Then the Ty/Tx ratio between node X and node 6. Experimental environment Y is used to divide the input data S. Namely, the ratio be- tween S and S is computed according to the relative per- To evaluate the different versions of BWA that have been de- formance achieved in each node. Hence, for the node X we scribed in this article and compare them with other versions pro- have size(S ) = size(S ) ∗ T y/(T x + T y ) and for the node Y we posed in the literature, we have used two datasets and three dif- have size(S ) = size(S ) ∗ T x/(T x + T y ). Obviously, this static ferent architectures. Below we show the details of these data sets distributions does not introduce any overhead at run-time and these architectures. but it requires the previous execution of several examples Table 1 shows the sample datasets we used in our experiments that to calculate the relative performance of each node of (SRR766060 and YH110112). The first one was obtained from the the system. This mechanism is illustrated at the upper part 10 0 0 Genomes Project [37] and the second one corresponds to a of Fig. 8. genome sequence of a Han Chinese individual [38]. HS37D5 corre- 2. Even distribution. sponds to the human reference genome. In this case, all instances of MDPR receive a block of data Table 2 summarizes the main characteristics of the three sys- of the same size. If the total number of instances is (x + y) tems used in our experiments (referred as S1, S2 and S3). Sandman and the total size of the dataset (S) is size(S), each instance and Sandman-Xeon Phi (S1) constitutes an heterogeneous system Table 1 Experiment dataset. Name Data set Data size Read length Read number – Short Read SRR766060(SRR) 17.4GB 100 34.2 M Paired-end – YH110112(YH) 4.4GB 100 14.8 M Segment Genome Reference HS37D5 3GB – – Human Genome Table 2 Experiment system. System Hostname Type Processor NUMA Core Thread Memory Data Width S1 Sandman Intel Xeon E5 2620 2 2 12 24 64GB 64bits Sandman-Xeon Phi Intel Xeon Phi 7120 1 1 60 240 16GB 512bits S2 Penguin Intel Xeon E5 4620 4 4 32 64 128GB 64bits Batman AMD Opteron Processor 6376 4 8 32 64 128GB 64bits S3 AMD+Intel AMD Opteron+Intel Xeon 17 25 126 252 576GB 64bits 18 S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 made of 2 Intel Xeon E5-2620 processors and one Intel Xeon Phi Table 3 BWA aligners on basic execution modes. 7120 accelerator. Sandman has 2 processors, 2-NUMA nodes (with 32GB of memory in each node), 12 cores and 24 available threads Architecture Native mode Symmetric mode Offload mode with hyper threading enabled. Sandman-Xeon Phi consists of 60 Xeon BWA-Xeon sBWA cores, 240 available threads and 16GB memory on card. Each core Xeon Phi BWA-Xeon Phi mBWA runs at 1333 Mhz and can process 4 hardware threads. The second system S2 that we used in our experiments consists of two multi- core servers: Penguin (equipped with Intel processors) and Batman of BWA running in symmetric mode with two instances of BWA (equipped with AMD processors). Penguin has 4 Xeon processors, aligner (one in Xeon and one in Xeon Phi) connected by MPI ser- 4-NUMA nodes, 32 cores, 64 threads and 32GB memory size per vices. Finally, mBWA [6], runs BWA in an offload mode. It starts in processor, and Batman has 4 sockets (with 2 processors per socket the Xeon part and the main time consuming parts are off-loaded to which implies a total of 8-NUMA nodes). Furthermore, a large- the Xeon Phi part. Fig. 9 shows average execution times obtained scale system (S3) that was used in MDPR experiments consisting by each strategy for each dataset. The longer column in one color of 6 nodes (combining AMD Opteron and Intel Xeon processors), represents the results from dataset YH, while the shorter column with a total of 126 physical cores and the capability of executing in the same color stands for the results from dataset SRR. 252 threads. The architectures of these systems are illustrated in Maximal parallelism in our system is 24 threads in Xeon Fig. 1. (Sandman) and 240 threads in Xeon Phi (Sandman-Xeon Phi). All experiments were executed using Intel Parallel Studio XE. Hence, BWA-Xeon case illustrates scalability of 6, 12, 18 and 24 And a minimum of 6 executions of each independent case with threads. BWA-Xeon Phi and mBWA cases use 60, 120, 180 and 240 the same environment and configuration were carried out for the threads. Execution of sBWA was done with the following number purpose of getting mean values with high confidence intervals. of threads: 6 (on the Xeon) and 60 (on the Xeon Phi), 12 and 120, 18 and 180, and 24 and 240, respectively. All strategies show 7. Experimental results and evaluation a similar tendency in each dataset. Despite the high number of threads used, BWA-Xeon Phi exhibits the worst execution times, In this section we show several experiments comparing differ- with sBWA being the best performing version. sBWA always out- ent versions of BWA in order to find the one that provides better performs native and offload cases. But, although sBWA simultane- performance. The first set of experiment shows a comparison be- ously takes advantage of available cores both on the Xeon and on tween the basic mode of BWA (native mode) and strategies from the Xeon Phi, the performance improvement is not linear with the the literature that run in offload mode and in symmetric mode number of threads. on Xeon Phi systems. The next experiments analyze the behavior These results can be explained because sequence alignment has of the strategies based on data parallelization or data replication large data I/O operations and is very sensitive to the time costs of described in Section 5.1. Finally, experiments carried out with the such operations. Reading the reference genome and the datasets, MDPR strategy are shown to evaluate the performance of the three and writing the final mappings are more expensive in terms of ex- mechanisms of data distribution both on our manycore system and ecution time in the Xeon Phi than in the Xeon. In the case of of- on our multicore systems. fload mode, mBWA performs main data I/O operation in the Xeon part instead of in the Xeon Phi part, and offloads highly parallel 7.1. Evaluation of parallelism on basic execution modes loops into the Xeon Phi part. Thus when comparing the perfor- mance of the same amount of threads in mBWA with that of BWA- These experiments evaluate the scalability of three basic ex- Xeon Phi, mBWA spends approximately two thirds of the execution ecution modes of BWA according to the description provided in time spent by the BWA-Xeon Phi. The same situation occurs with Section 3 (i.e., native, offload and symmetric modes). We ran the BWA-Xeon and BWA-Xeon Phi, the former case reaches better per- applications on our manycore system (Sandman and Sandman- formance than the latter one, despite a lower amount of threads Xeon Phi, as described in Table 2) and we used two datasets. We applied in BWA-Xeon. In the case of sBWA, the asymmetric distri- tested the BWA versions shown in Table 3. BWA-Xeon and BWA- bution of initial data compensates the loss in performance incurred Xeon Phi correspond to the original version of BWA running in na- by data I/O operations on the Xeon Phi. tive mode either in the Xeon part of the system or in the Xeon Fig. 9 also illustrates another problem that affects all four align- Phi part of the system. sBWA is a straightforward implementation ers: scalability does not improve significantly as a growing number Fig. 9. Scalability of BWA aligners (Identical color in short column: SRR and in long column: YH). S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 19 Fig. 10. Scalability comparison on Sandman. Fig. 11. Scalability comparison on Sandman-Xeon Phi. of threads is used. As mentioned in Section 4, BWA uses a shared is not so significant because there are no clearly-separated NUMA data structure containing the reference genome. This data struc- nodes and the performance is improved until a point is reached ture is loaded at one particular memory bank of the system by where the ring connecting all the memory banks is saturated. the master thread. Thus this memory bank turns into a congestion bottleneck of performance as the number of threads grows. This 7.2.2. Analysis of thread groups and data replicas in DP, DR and DPR memory bank needs to support more concurrent accesses coming One of the factors associated with the design of the DP, DR and from different threads to populate the corresponding caches. DPR strategies is the number of thread groups or the number of replicas used by them. With the experiments shown in this subsec- 7.2. Evaluation of data parallelism strategies tion we want to answer questions like what is the ideal number of groups or replicas that should be used in a particular architecture? In this subsection, we show the results obtained in the evalua- or can we use any practical rule that can be applied to deduce such tion of the three different data parallelism strategies described in number? Section 5.1. All strategies were executed in system S1. We executed DP, DR and DPR varying the number of thread groups and replicas that were used (x and y). In all cases, we used 7.2.1. Scalability of DP, DR and DPR the maximal number of available cores. This means that 24 and We set up two thread groups for execution in the strategy of 240 threads were used in all the cases on the Xeon part (Sand- DP, DR and DPR, namely n = 2 in DP, m = 2 in DR and n = m = 2 man) and on the Xeon Phi part (Sandman-Xeon Phi), respectively. in DPR. Figs. 10 and 11 show relative speedups on Intel Xeon and The number of OpenMP threads used at each instance (tx and ty) Xeon Phi, respectively. The speed-up is calculated under dataset on the Xeon and Xeon Phi are tx = 24/x and ty = 240/y. For in- SRR using as a baseline case the execution time of original BWA stance, if 2 thread groups are applied at the whole system, each with one thread running on Sandman. Xeon thread group utilizes 12 OpenMP threads and each Xeon Phi In general, the three strategies have a consistent behavior in thread group utilizes 120 OpenMP threads. Thus, the number of all cases and it can be concluded that DPR obtains the best per- thread groups 2, 3, 4, 6 and 12, could stand for the number of formance, followed by DP and DR. Results in the Xeon exceed by OpenMP threads 12, 8, 6, 4 and 2 in each thread group in Fig. 12 a factor of 10 approximately those obtained in the Xeon Phi. The and 120, 80, 60, 40 and 20 in each thread group in Fig. 13. For overhead of I/O explains part of these lack of performance, as it DR and DPR, the number of thread groups is limited by memory was mentioned in the previous section. This loss in performance is consumption of each group. Basically, one thread group of BWA also explained because BWA does not take advantage of the 512-bit aligner consumes about 5GB memory, which is mainly required to SIMD units available in the Xeon Phi. On the other hand, the effect store the human genome sequence reference (4.4GB). According to of memory contention is attenuated thanks to the use of replicas the memory size available on the Xeon part (64GB) and on the of the reference genome. This effect is greater in the case of the Xeon Phi part (16GB), 12 and 3 should be the maximum number Xeon system because it has two clearly separated NUMA nodes. of thread groups that can be launched at each part. We tested DP, In the case of Xeon Phi, the use of replicas allows increasing the DR and DPR using thread groups of 2, 3, 4, 6 and 12 thread groups number of concurrent accesses. But in this case, the improvement on both systems, but in the cases of DPR-BWA and DR-BWA only 3 20 S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 Fig. 12. Performance comparison on Sandman. Fig. 13. Performance comparison on Sandman-Xeon Phi. Fig. 14. BWA memory access. Fig. 15. DPR-BWA memory access. thread groups were applied in the experiment of Xeon Phi due to each memory bank. Threads that belong to each group also take the memory limitations. advantage of those cache entries that contain the index positions Results are shown in Figs. 12 and 13. It can be observed that which are read in the first stages of the alignment. By increasing the best performance is achieved when a small number of groups the number of thread groups, the usage of index entries in the are used (2 or 3). Beyond these values, adding more groups wors- cache memory gets worse and there is an increasing number of ens performance. These results seem coherent if we take into ac- cache misses that does not compensate for the reduction of mem- count the number of NUMA nodes in the system (2 in the case of ory congestion. Sandman and 1 in the case of Sandman Phi). Having more repli- The performance improvements observed in these experiments cas helps distribute traffic between different memory banks and are based on the idea that each group of threads will make ac- prevents a single bank from becoming congested. However, adding cesses to the memory node of its NUMA domain, thus taking ad- more replicas increases the overhead and reduces the performance vantage of the greatest available bandwidth. Using the NUMATOP of the cache because there is more competition between different tool, we obtained measurements on the distribution of memory ac- groups of threads without increasing the transference capacity of cesses between the different memory nodes. Figs. 14 and 15 show S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 21 the behavior of the original version of BWA and our DPR strategy previous sections). The second one (referred as S2) is made of an when they were executed in Sandman, a system that has 2 NUMA Intel-based server and an AMD-based server. Furthermore, we also nodes formed by a processor and a memory bank. The two figures tested MDPR in a large system (S3) that comprises 6 nodes (made show, on the one hand, the proportion of local memory accesses of Intel Xeon and AMD Opteron processors), 17 sockets, 25 NUMA (LMA) and the speedup of the application compared to the sequen- domains, capable of running 252 threads. tial execution. As can be seen in Fig. 14, BWA exhibits a proportion According to the results obtained in the evaluation of DPR with of local memory access greater than 80% if it is executed in a sin- different number of replicas the following configurations of MDPR gle processor (12 threads). were tested in S1: 1+1(2), 2+1(3), 2+2(4), 2+3(5) and 4+2(6). When increasing the parallelism, the additional threads are ex- The first value corresponds to the number of replicas in the Xeon ecuted in the second processor and all of them will generate re- part and the second value corresponds to the number of replicas mote accesses of memory because they have to access the node in the Xeon Phi. System S2 has 4+4 sockets and 4+8 NUMA do- where the reference index is located. The LMA index decreases to mains, while system S3 has a total of 17 sockets and 25 NUMA do- just over 40%. As a result, remote accesses have a negative effect on mains. Thus, the configurations evaluated in S2 are 1+1(2), 2+2(4), the execution time of the application because latency is higher for 4+4(8) and 4+8(12); in S3 the configurations were 6, 9, 17 and 25. remote accesses and they increase congestion in the node where All available hardware threads were used in all cases of MDPR ex- the index is located. This negative impact on the execution time periments (i.e., Intel Xeon and Xeon Phi could use 24 threads plus explains that when doubling from 12 to 24 threads the speed-up 240 threads, the S2 system could support 64 threads each, and S3 only improves by a factor of just over 2. could provide 252 threads). On the contrary, the use of instances (as shown in Fig. 15) in As shown in Figs. 16–18, the static distribution of MDPR DPR always maintains a proportion of local accesses above 80% al- achieves the best execution times in all cases. The even distribu- though the number of instances used is increased. There is no neg- tion is the worst strategy since it only makes an initial distribution ative impact on memory access times nor does congestion occur. of the data without subsequent corrections according the relative Anyway, as we saw earlier, the best speedup (factor 12, approxi- performance of the nodes. As an intermediate strategy, the round- mately) is obtained when the number of instances is low (2 or 4, robin method obtains results close to the static one without the which means one or two instances executed in each NUMA node). need for preliminary setup computations. On the other hand, if we If the number of instances is too high (6 or more) the groups of analyze the number of replicas used, it is observed that the best threads compete and interfere with the use of the cache memory results are obtained when the number of replicas is adjusted to without improving the access times to the main memory. the number of NUMA domains available in the system. This be- havior has already been observed in the DPR experiments on the 7.2.3. Performance evaluation of MDPR manycore system and is confirmed in the execution of MDPR in This last set of experiments analyzes the performance of the systems S2 and S3. MDPR strategy, evaluating the three data distribution policies de- The results also show that, although the even distribution is scribed in Section 5.2. the one that achieves worst results in all scenarios, this phe- MDPR-BWA allows the execution of BWA aligner in the sym- nomenon is more significant when applied in heterogeneous sys- metric mode and it can also be used on distributed memory sys- tems where the differences in performance between its compo- tems. The best execution time should be achieved if reads are dis- nents is very significant. In our case, this is observed in the S1 tributed so that all the threads finish at the same time. The distri- system, where the Xeon Phi part is significantly slower than the bution of reads can be done in different ways taking into account Xeon part (as seen in previous experiments). Therefore, when the that there is a trade-off between the overhead incurred and the same amount of data is distributed to each instance, the execu- load balance that can be achieved. In our case, we have evaluated tion in the instances running in the slowest component will de- a static mechanism that requires some benchmarking executions termine the overall execution time. As the degree of heterogene- beforehand but does not incur any overhead at run time. We have ity decreases and the number of instances increases, the perfor- also tested two dynamic methods (even distribution and round- mance differences of both dynamic data distribution strategies are robin), which incur in more overhead in execution time but do not close to the performance of the static strategy. If we take into require any prior benchmarking. account the cost incurred by the static strategy to compute the In these experiments, we used the same datasets used in previ- relative performance of each computing node, the use of the dy- ous experiments (SRR+YH) and three different heterogeneous sys- namoc round-robin strategy constitutes a good trade-off between tems. The first one (referred as S1) is the one equipped with an the reduction in execution times and the overheads added at Intel Xeon and Xeon Phi (Sandman and Sandman-Xeon Phi used in run-time. Fig. 16. MDPR comparison on Xeon and Xeon Phi(S1). 22 S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 Fig. 17. MDPR comparison on Intel and AMD NUMA nodes(S2). Fig. 18. MDPR comparison on Intel and AMD NUMA nodes(S3). 8. Discussion and conclusion static strategy that initially distributes the entire data set taking into account the performance of each node of the system, which In this paper, we have evaluated different data parallelism were previously evaluated through a benchmarking process. How- strategies to execute BWA aligner in heterogeneous architectures ever, we believe that the use of a dynamic strategy that distributes made of multicore and/or manycore nodes. BWA has scalability data chunks in a round-robind fashion is the best practical alterna- limitations when large numbers of threads are used. The use of tive because the loss of performance is not significant and its use a shared reference genome generates problems of memory con- does not require the need to perform any previous benchmarking tention that explain such scalability issues. of nodes. We have analized several strategies in order to reduce the nega- First, we evaluate different versions of BWA that run on many- tive impact incurred by BWA due to memory contention. The basic core systems in native mode, offload mode and symmetric mode. idea that inspires our strategies is to divide the application into Second, we analyze different versions of BWA that using different a series of groups (or instances) that will cooperate in the tasks strategies to parallelize the distribution of short reads or to repli- of read alignment; each group is formed by a certain number of cate the reference genome. threads that will process each chunk of data in which the original Tables 4–7 summarize the best results achieved in all our ex- dataset is divided. In addition, the application can use data repli- periments by each BWA-based aligner presented in this paper cation. In our case, this serves to replicate the reference genome in when they were executed with the maximum number of resources. each group or instance. The strategy that incorporates all these ele- Each entry in the tables shows the aligner that was used, the sys- ments is the one that we referred as MDPR and is the one that has tem on which it was executed, the dataset, the average execution obtained the best results, in general. MDPR can be used in any het- time that was obtained, the corresponding speed-up, the number erogeneous systems based on manycore and/or multicore nodes. of groups or instances used and the total number of threads. The From our experiments, we have observed that the best MDPR speed up is obtained by comparison with the original implemen- results are obtained when the number of instances is equal to tation of BWA (BWA-Pthread) using one thread with SRR dataset the number of NUMA domains of the underlying architecture. In (first entry in Table 4). this way, the structure of the application is being adjusted to the Tables 4 and 5 show the results obtained by different BWA ver- characteristics of the different nodes that build the whole system: sions executed using a native o a offload method in the Xeon part there are as many groups as there are NUMA domains and each or in the Xeon-Phi part, respectively, of our S1 system. In both group executes with as many threads as cores the processor of that cases we have also compared an official BWA implementation (re- domain has. ferred as BWA-ALN-Xeon-Phi) from Intel Corporation [39]. BWA- In our study we have also analyzed three strategies to perform ALN-Xeon-Phi was tested in native mode (Tables 4 and 5) and the distribution of data among the different groups that are cre- it was also tested in symmetric mode in S1 cluster, as shown in ated in MDPR. The best execution times were obtained by the Table 6. S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 23 Table 4 BWA performance summary on Intel Xeon. Aligner System Dataset Time-cost(s) Speed-up Group num Thread num BWA-Pthread S1-Xeon SRR 18,658 1X 1 1 BWA-ALN-Xeon-Phi S1-Xeon SRR 1873 10X 1 24 YH 2985 6.3X 1 24 DP S1-Xeon SRR 1580 11.8X 2 24 YH 2591 7.2X 2 24 DR S1-Xeon SRR 1740 10.7X 2 24 YH 2784 6.7X 2 24 DPR S1-Xeon SRR 1392 13.4X 3 24 YH 1777 10.5X 3 24 Table 5 BWA Performance Summary on Intel Xeon Phi. Aligner System Dataset Time-cost(s) Speed-up Group num Thread num BWA-ALN-Xeon-Phi S1-Xeon Phi SRR 2895 6.4X 1 240 YH 4623 4X 1 240 mBWA S1-Xeon Phi SRR 2427 7.7X 1 240 YH 3817 4.9X 1 240 DP S1-Xeon Phi SRR 2630 7.1X 3 240 YH 3395 5.5X 3 240 DR S1-Xeon Phi SRR 3347 5.6X 2 240 YH 4118 4.5X 2 240 DPR S1-Xeon Phi SRR 2555 7.3X 2 240 YH 2955 6.3X 3 240 Table 6 BWA performance summary on S1. Aligner System Dataset Time-cost(s) Speed-up Group num Thread num BWA-ALN-Xeon-Phi S1 SRR 1156 16.1X 2(1+1) 24+240 YH 1733 10.8X 2(1+1) 24+240 sBWA S1 SRR 1301 14.3X 2(1+1) 24+240 YH 2047 9.1X 2(1+1) 24+240 MDPR-static S1 SRR 966 19.3X 5(2+3) 24+240 YH 1477 12.6X 5(2+3) 24+240 MDPR-even S1 SRR 1501 12.4X 5(2+3) 24+240 YH 2443 7.6X 5(2+3) 24+240 MDPR-roundrobin S1 SRR 1086 17.2X 5(2+3) 24+240 YH 1584 11.8X 5(2+3) 24+240 Table 7 BWA performance summary on S2 and S3. Aligner System Dataset Time-cost(s) Speed-up Group num Thread num MDPR-static S2 SRR 204 91.5X 12(4+8) 64+64 YH 385 48.5X 12(4+8) 64+64 S3 SRR 109 171.2X 25 252 YH 175 106.6X 25 252 MDPR-even S2 SRR 260 71.8X 8(4+4) 64+64 YH 535 34.9X 8(4+4) 64+64 S3 SRR 167 111.7X 25 252 YH 212 88X 17 252 MDPR-roundrobin S2 SRR 245 76.2X 12(4+8) 64+64 YH 411 45.4X 12(4+8) 64+64 S3 SRR 145 128.7X 17 252 YH 195 95.7X 17 252 As seen in Tables 4 and 5, our strategy DPR that use data repli- MDPR-roundrobin achieves values of 17.2 and 11.8. It is worth not- cation and data parallelization obtain better execution times in ing that although our MDPR approach allows a significant improve- general. DPR achieved a speed-up of 13.4 and 10.5 for SRR and YH, ment of performance in manycore systems based on Intel Xeon respectively, in the Xeon case, and a speed-up of 7.3 and 6.3 in the Phi, the characteristics of this accelerator (limited I/O capacity and Xeon-Phi case. Only mBWA achieved a slightly better result in the main memory interconnected by a ring bus) limit its potential per- case of SRR dataset with a seep-up of 7.7. formance in comparison with multicore systems with NUMA ar- Table 6 shows the best results obtained by the strategies run- chitecture, if we take into account the total number of cores used. ning in symmetric mode on the entire S1 system. In this case, two Speed-up values obtained at the S1 system are significantly lower versions of MDPR (with static distribution and with round-robin) that those obtained at S2 and S3 systems (see Table 7), but the to- obtain better results than the other two strategies that do not use tal number of cores in S1 was greater than that used in S2 and S3. data replication (BWA-ALN-Xeon-Phi) and sBWA. The best strategy Finally, Table 7 shows the results of the different MDPR variants is MDPR-static with speed-ups of 19.3 and 12.6, respectively, while when running on larger systems formed by several heterogeneous 24 S. Chen and M.A. Senar / Parallel Computing 87 (2019) 11–24 multicore nodes. MDPR with static distribution achieves the best [15] N. Kathiresan, M.R. Temanni, R. Al-Ali, Performance improvement of BWA MEM results, followed by the round-robin version. algorithm using data-parallel with concurrent parallelization, in: Parallel, Dis- tributed and Grid Computing (PDGC), 2014 International Conference on, IEEE, From all the results obtained, we can conclude that the MDPR 2014, pp. 406–411. strategy is a very attractive alternative to reduce the execution [16] E. de Araujo Macedo, A. Boukerche, Hybrid MPI/OpenMP strategy for biolog- time of BWA in the distributed systems, either by using a static ical multiple sequence alignment with DIALIGN-TX in heterogeneous multi- core clusters, in: Parallel and Distributed Processing Workshops and Phd Fo- distribution or a dynamic round-robin pattern. rum (IPDPSW), 2011 IEEE International Symposium on, IEEE, 2011, pp. 418– Our study also reveals an idea that can be applied to other 425. problems in the field of bioinformatics. Although in this field there [17] S.L. Olivier, A.K. Porterfield, K.B. Wheeler, J.F. Prins, Scheduling task parallelism on multi-socket multicore systems, in: Proceedings of the 1st International is a great growth in the quantity and quality of existing data sets, Workshop on Runtime and Operating Systems for Supercomputers, ACM, 2011, in practice, the resolution of problems such as sequence alignment pp. 49–56. are limited mainly by the size of reference genomes (usually be- [18] J. Lenis, M.A. Senar, A performance comparison of data and memory allocation strategies for sequence aligners on NUMA architectures, Cluster Comput. 20 (3) tween 3GB and 30GB) which are the structure that should remain (2017) 1909–1924. in main memory all the time. Data files with sequences are read [19] Y.-T. Chen, J. Cong, J. Lei, P. Wei, A novel high-throughput acceleration en- gradually in the form of chunks and, therefore, they do not sup- gine for read alignment, in: Field-Programmable Custom Computing Machines pose any significant limitation in terms of memory requirements in (FCCM), 2015 IEEE 23rd Annual International Symposium on, IEEE, 2015, pp. 199–202. current systems. Consequently, the design of parallel applications [20] E.J. Houtgast, V.-M. Sima, K. Bertels, Z. Al-Ars, Gpu-accelerated BWA-MEM by using an architecture based on instances that contain replicas genomic mapping algorithm using adaptive load balancing, in: International of certain data structures constitutes a solution with low algorith- Conference on Architecture of Computing Systems, Springer, 2016, pp. 130– 142. mic complexity that can provide significant improvements in sys- [21] E.J. Houtgast, V.-M. Sima, K. Bertels, Z. Al-Ars, An FPGA-based systolic array tems whose memory architecture is of the NUMA type, as shown to accelerate the BWA-MEM genomic mapping algorithm, in: Embedded Com- by MDPR in the case of BWA. puter Systems: Architectures, Modeling, and Simulation (SAMOS), 2015 Inter- national Conference on, IEEE, 2015, pp. 221–227. [22] H.-C. Ng, S. Liu, W. Luk, Reconfigurable acceleration of genetic sequence align- Acknowledgements ment: a survey of two decades of efforts, in: Field Programmable Logic and Applications (FPL), 2017 27th International Conference on, IEEE, 2017, pp. 1– 8. This work is supported by the Spanish Government [23] H. Li, N. Homer, A survey of sequence alignment algorithms for next-genera- Project TIN201784553C21R] and China Scholarship Council tion sequencing, Brief. Bioinformatics 11 (5) (2010) 473–483. [201406890 0 07]. [24] E.J. Houtgast, V.-M. Sima, G. Marchiori, K. Bertels, Z. Al-Ars, Power-efficiency analysis of accelerated BWA-MEM implementations on heterogeneous comput- ing platforms, in: ReConFigurable Computing and FPGAs (ReConFig), 2016 In- References ternational Conference on, IEEE, 2016a, pp. 1–8. [25] E.J. Houtgast, V.-M. Sima, K. Bertels, Z. Al-Ars, GPU-accelerated BWA-MEM [1] C.F. Baes, M.A. Dolezal, J.E. Koltes, B. Bapst, E. Fritz-Waters, S. Jansen, C. Flury, genomic mapping algorithm using adaptive load balancing, in: International H. Signer-Hasler, C. Stricker, R. Fernando, et al., Evaluation of variant identifi- Conference on Architecture of Computing Systems, Springer, 2016b, pp. 130– cation methods for whole genome sequencing data in dairy cattle, BMC Ge- 142. nomics 15 (1) (2014) 948. [26] J.M. Abuín, J.C. Pichel, T.F. Pena, J. Amigo, BigBWA: approaching the bur- [2] I. Buck, Gpubench: evaluating GPU performance for numerical and scien- rows–wheeler aligner to big data technologies, Bioinformatics 31 (24) (2015) tific application, in: Proc. 1st ACM Workshop General-Purpose Computing on 40 03–40 05. Graphics Processors (GP^2’04), 2004. [27] J.M. Abuín, J.C. Pichel, T.F. Pena, J. Amigo, Sparkbwa: speeding up the alignment [3] L. Stein, Genome annotation: from sequence to biology, Nat. Rev. Genet. 2 (7) of high-throughput dna sequencing data, PLoS ONE 11 (5) (2016) e0155461. (2001) 493. [28] S.-H. Chan, J. Cheung, E. Wu, H. Wang, C.-M. Liu, X. Zhu, S. Peng, R. Luo, T.-W. [4] J.G. Reid, A. Carroll, N. Veeraraghavan, M. Dahdouli, A. Sundquist, A. English, Lam, MICA: a fast short-read aligner that takes full advantage of intel many M. Bainbridge, S. White, W. Salerno, C. Buhay, et al., Launching genomics integrated core architecture (mic), arXiv:1402.4876 (2014). into the cloud: deployment of mercury, a next generation sequence analysis [29] X. Tian, H. Saito, S.V. Preis, E.N. Garcia, S.S. Kozhukhov, M. Masten, pipeline, BMC Bioinformatics 15 (1) (2014) 30. A.G. Cherkasov, N. Panchenko, Practical SIMD vectorization techniques for in- [5] S. Pabinger, A. Dander, M. Fischer, R. Snajder, M. Sperk, M. Efremova, B. Kra- tel® xeon phi coprocessors, in: Parallel and Distributed Processing Symposium bichler, M.R. Speicher, J. Zschocke, Z. Trajanoski, A survey of tools for variant Workshops & PhD Forum (IPDPSW), 2013 IEEE 27th International, IEEE, 2013, analysis of next-generation genome sequencing data, Brief. Bioinformatics 15 pp. 1149–1158. (2) (2014) 256–278. [30] S. Memeti, S. Pllana, Accelerating dna sequence analysis using intel (r) xeon [6] Y. Cui, X. Liao, X. Zhu, B. Wang, S. Peng, mBWA: a massively parallel se- phi (tm), in: Trustcom/BigDataSE/ISPA, 2015 IEEE, volume 3, IEEE, 2015, quence reads aligner, in: 8th International Conference on Practical Applica- pp. 222–227. tions of Computational Biology & Bioinformatics (PACBB 2014), Springer, 2014, [31] L. Wang, Y. Chan, X. Duan, H. Lan, X. Meng, W. Liu, XSW: accelerating biolog- pp. 113–120. ical database search on xeon phi, in: Parallel & Distributed Processing Sym- [7] H. Li, R. Durbin, Fast and accurate short read alignment with burrows–wheeler posium Workshops (IPDPSW), 2014 IEEE International, IEEE, 2014, pp. 950– transform, Bioinformatics 25 (14) (2009) 1754–1760. 957. [8] S. Chen, M.A. Senar, Accelerating BWA aligner using multistage data paral- [32] L. Ping, Speeding up large-scale next generation sequencing data analysis with lelization on multicore and manycore architectures, Procedia Comput. Sci. 80 pBWA, J. Appl. Bioinf. Comput. Biol. (2012). (2016) 2438–2442. [33] C. Herzeel, T.J. Ashby, P. Costanza, W. De Meuter, Resolving load balancing is- [9] S. Chen, M.A. Senar, Improving performance of genomic aligners on intel xeon sues in BWA on NUMA multicore architectures, in: International Conference on phi-based architectures, in: 2018 IEEE International Parallel and Distributed Parallel Processing and Applied Mathematics, Springer, 2013, pp. 227–236. Processing Symposium Workshops (IPDPSW), IEEE, 2018, pp. 570–578. [34] T. Cramer, D. Schmidl, M. Klemm, D. an Mey, Openmp programming on intel [10] S.F. Altschul, T.L. Madden, A.A. Schäffer, J. Zhang, Z. Zhang, W. Miller, D.J. Lip- r xeon phi tm coprocessors: an early performance comparison, in: Proc. Many man, Gapped blast and psi-blast: a new generation of protein database search Core Appl. Res. Community (MARC) Symp, 2012, pp. 38–44. programs, Nucleic Acids Res. 25 (17) (1997) 3389–3402. [35] J. Reinders, An overview of programming for intel xeon processors and intel [11] B. Langmead, C. Trapnell, M. Pop, S.L. Salzberg, Ultrafast and memory-efficient xeon phi coprocessors, Intel Corporation, Santa Clara, 2012. alignment of short dna sequences to the human genome, Genome Biol. 10 (3) [36] Y. Cui, X. Liao, X. Zhu, B. Wang, S. Peng, B-Mic: an ultrafast three-level parallel (2009) R25. sequence aligner using mic, Interdiscip. Sci. 8 (1) (2016) 28–34. [12] R. Li, C. Yu, Y. Li, T.-W. Lam, S.-M. Yiu, K. Kristiansen, J. Wang, Soap2: an im- [37] .G.P. Consortium, et al., A global reference for human genetic variation, Nature proved ultrafast tool for short read alignment, Bioinformatics 25 (15) (2009) 526 (7571) (2015) 68. 1966–1967. [38] J. Wang, Y. Li, R. Luo, B. Liu, Y. Xie, Z. Li, X. Fang, H. Zheng, J. Qin, B. Yang, [13] H. Li, R. Durbin, Fast and accurate short read alignment with burrows–wheeler et al., Updated genome assembly of yh: the first diploid genome sequence of transform, Bioinformatics 25 (14) (2009) 1754–1760. a han chinese individual (version 2, 07/2012), GigaScience Database (2012). [14] J. Zhang, H. Lin, P. Balaji, W.-c. Feng, Optimizing burrows-wheeler transfor- [39] C.I. You, Liang (Intel); Congdon, Building and optimizing BWA-ALN 0.5.10 m-based sequence alignment on multicore architectures, in: Cluster, Cloud and for intel xeon phi coprocessors, URL https://software.intel.com/en-us/ Grid Computing (CCGrid), 2013 13th IEEE/ACM International Symposium on, articles/recipe- building- and- optimizing- bwa- aln- 0510- for- intel- xeon- phi- IEEE, 2013, pp. 377–384. coprocessors (2014).