arXiv:2606.25738v1 [cs.DC] 24 Jun 2026
Endeavor: Efficient PairHMM for Detection of DNA Variants in Genome-Scale Datasets Miguel Graça
Aleksandar Ilic
[email protected] INESC-ID, Instituto Superior Técnico Lisbon, Portugal
[email protected] INESC-ID, Instituto Superior Técnico Lisbon, Portugal
Abstract
1
DNA variant calling represents a key operation in bioinformatics pipelines that aims at identifying genetic variants. Given an evidenced explosion in genomic data availability, there is an urgent need for a high-performant, portable and efficient solution for variant calling, which can further improve our understanding of genomic structure and genetic basis for complex diseases. In its most common formulation, the Pair Hidden Markov Model (PairHMM) algorithm for variant calling stands as the main bottleneck in the pipeline, accounting for up to 70% of the execution time in large-scale genomic datasets. The state-of-the-art approaches for accelerating PairHMM in CPUs and GPUs do not scale to long DNA sequences and only explore very limited anti-diagonal data parallelism, which yields poor performance. In this work, Endeavor is proposed as a new parallelization strategy for PairHMM that redefines its traditional formulation to explore row-level fine-grained parallelism without loss in solution accuracy. Based on this, a novel and portable SIMD-based approach is derived for efficient and highperformance processing of short and long sequences in CPUs and GPUs, leveraging novel levels of parallelism and synchronization to achieve high throughput in sequences up to 100k basepairs for the first time. Evaluation on Intel and AMD CPUs shows that Endeavor outperforms GKL up to 2.14x in peak throughput and GATK HaplotypeCaller by at least 2x in real-world datasets, while NVIDIA and AMD GPUs achieve up to 2.05x speedups in genome-scale datasets when compared to state-of-the-art GPU-based methods.
The Human Genome Project [18], set to sequence more than 90% of the human genome, represented a milestone in genomics that took 13 years to achieve (completed in 2003). During the following 20 years, Next-Generation Sequencing (NGS) technologies [20] have paved the way to sequence genomes faster and cheaper. This rapid technological advance culminated in the development of more ambitious genetic projects, such as the 1000 Genomes Project [12] and the Genome In A Bottle (GIAB) consortium [50], as well as an unprecedented explosion of the available genomic data. The widespread availability of such data led to the emergence of a novel bioinformatics pipeline, known as variant calling [39], which analyses the existence of genetic dissimilarities between the patient’s genome and a reference genome. Variant calling is of the utmost importance to identify genetic differences and understand their correlation with complex diseases [25], which has a significant impact in precision medicine [3], pharmacogenomics [17], and evolutionary biology [23]. As a consequence, specific toolkits have been developed to process genetic data for variant calling, with the Genome Analysis Toolkit (GATK) HaplotypeCaller [13] being one of the most widely used in genomic analysis. The typical workflow of variant calling in GATK follows three steps. Given a dataset with read sequences from a patient, the first step is to generate a set of candidate haplotypes by assembling the reads [41]. The second step is to calculate the probability that a read was generated by a given haplotype. To this end, the Pair Hidden Markov Model (PairHMM) algorithm [51] aligns the read sequence data and the candidate haplotypes, calculating the alignment likelihoods for each candidate and selecting the haplotypes with the highest likelihood. Finally, a Bayesian model is applied to the computed likelihoods to determine the most probable genotype at a given read (if it is equal to a reference genome or a variant). Of these three steps, the PairHMM algorithm is the main bottleneck, responsible for more than 70% of the execution time in genome-scale datasets [26], and therefore is the focus of hardware acceleration efforts. However, the PairHMM algorithm is not trivial to accelerate and the main reasons for this are two-fold. First, PairHMM scales quadratically in time and memory with the sequence size (as it is based on computing matrix cells), making it prohibitive for sequences generated with long-read sequencing technologies. Although memory requirements can be reduced to scale linearly with sequence length [44], the same cannot be done for time requirements. Second, the algorithm uses dynamic programming [15], which hinders the opportunities to parallelize the algorithm due to dependencies in the data flow. Nevertheless, various works in the literature attempted to parallelize PairHMM by leveraging the anti-diagonal independence and
CCS Concepts • Applied computing → Computational biology; • Computer systems organization → Parallel architectures; • Computing methodologies → Parallel algorithms.
Keywords Variant Calling, Parallel Computing, Bioinformatics ACM Reference Format: Miguel Graça and Aleksandar Ilic. 2026. Endeavor: Efficient PairHMM for Detection of DNA Variants in Genome-Scale Datasets. In The 35th International Symposium on High-Performance Parallel and Distributed Computing (HPDC ’26), July 13–16, 2026, Cleveland, OH, USA. ACM, New York, NY, USA, 13 pages. https://doi.org/10.1145/3806645.3807593
This work is licensed under a Creative Commons Attribution 4.0 International License. HPDC ’26, Cleveland, OH, USA © 2026 Copyright held by the owner/author(s). ACM ISBN 979-8-4007-2640-8/2026/07 https://doi.org/10.1145/3806645.3807593
Introduction
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
implementing solutions on FPGAs [4, 21, 44] and GPUs [6, 26, 42]. However, the current solutions suffer from three major drawbacks. First, CPUs are still the de-facto standard in various bioinformatics applications and tools [7, 8, 29, 52]. While highly-optimized PairHMM CPU-based implementations [13] do not offer the performance of massively parallel GPUs, they can evaluate long-read datasets (where the sequence length can exceed 104 to 105 basepairs [22]). Second, the current GPU-based solutions are not scalable with the sequence length, i.e., they cannot deal with long-read datasets, and are not portable. The existent GPU approaches target only NVIDIA GPUs which, given the current heterogeneous computational landscape, lack the necessary portability to target different hardware devices (e.g., the top 5 fastest supercomputers in the world [47] include AMD MI300A GPUs, which cannot be targeted for PairHMM acceleration with the existing solutions from the literature). Finally, the anti-diagonal parallelization strategy is not efficient on both architectures, due to the variable amount of computations that the algorithm performs as it traverses the matrices’ anti-diagonals. Therefore, one of the main research challenges nowadays is how to further improve our understanding of genomes [31] by efficiently processing data from both shortand long-read sequencing technologies and leveraging different hardware architectures to achieve portability. This is a particular gap that this paper intends to close by proposing Endeavor (EfficieNt PairHMM for DEtection of DNA VARiants), a novel parallelization approach for the PairHMM algorithm on CPUs and GPUs that fully leverages the hardware resources and scales to long DNA sequences. The main contributions of this paper are as follows: (1) Redefinition of the processing steps of the PairHMM algorithm to expose previously unexplored levels of parallelism with no loss in solution accuracy; (2) An efficient, portable, and fine-grained parallel approach to achieve high performance on CPUs and GPUs; (3) Scalability to diverse genomic datasets with short- and longread sequences. Contribution (1) pertains to the novel parallelization strategy that is herein proposed to accelerate PairHMM, which leverages row-level parallelism of the matrices to fully exploit performance at the hardware level (with no loss in solution precision when compared to the anti-diagonal approach). Contribution (2) is concerned with the efficient mapping of the proposed approach to extract data parallelism. To benefit from the row-level parallelism on the CPU, the proposed solution focuses on multi-threading and SIMD intrinsics to process read-haplotype pairs at the level of each SIMD lane. On the GPU, Endeavor focuses on shared memory and cross-thread intrinsics to efficiently map the PairHMM computations, which are implemented in CUDA and HIP to target GPUs from major vendors (NVIDIA and AMD). To achieve unprecedented performance for long reads, novel features from recent GPU architectures, namely inter-thread-block operations, are also efficiently leveraged. Contribution (3) focuses on Endeavor’s portability and scalability across different devices. Compared to the Intel Genomics Kernel Library (GKL) [16], Endeavor achieves up to 2.14x improvements in peak throughput in different CPU architectures, while
Miguel Graça and Aleksandar Ilic
evaluation on genome-scale datasets shows that Endeavor outperforms a highly optimized (OpenMP + AVX512) GATK HaplotypeCaller’s PairHMM implementation on CPUs by at least 2x in GIAB datasets. Compared with a state-of-the-art PairHMM GPU-based method (referred herein as gpuPairHMM [45]) that only employs single floating-point precision on vendor-specific GPUs, Endeavor achieves speedups up to 2.05x on large-scale datasets, while providing support for both single- and double-precision, portability across different GPU architectures, and processing sequences up to 100k basepairs, representing an increase by 1 order of magnitude compared to the literature and demonstrating its flexibility to process short- and long-read datasets.
2
Background and Related Work
Variant calling aims to identify genetic variants by comparing a DNA sequence (a string of nucleotides) from a patient (herein referred to as a read sequence) and a reference DNA sequence (herein referred to as a haplotype sequence) [25]. To achieve this, the PairHMM algorithm computes all possible alignments between the read and the haplotype. An alignment describes the relationship between two sequences, where characters on both strings can be aligned to each other or to an empty space (herein referred to as gap) due to inserted or deleted characters on one of the strings. As a probabilistic model, PairHMM attributes a probability for each alignment [14]. The algorithm’s output is the aggregated probability of all alignments, which describes the likelihood that the read sequence is derived from the haplotype sequence. To perform this calculation, PairHMM genomic datasets provide the read-haplotype pairs, and four additional strings (of the same length as the read). These strings provide the read’s quality score (the probability that a given character in the read is erroneous), the insertion and deletion probabilities (the likelihood of an inserted or deleted character in the read when compared to the haplotype), and the gap probabilities in the sequences. Given a read and haplotype sequences, of length 𝑙 and 𝑘, respectively, the PairHMM algorithm operates on three matrices (𝑀, 𝐼 , and 𝐷) and is initialized as follows 𝑀𝑖,0 = 𝐼𝑖,0 = 𝐷𝑖,0 = 𝑀0,𝑗 = 𝐼 0,𝑗 = 0, 𝐷 0,𝑗 = 21020 /𝑘,
(1)
where the first row and column of 𝑀 and 𝐼 are initialized to 0 and 𝐷 is initialized to a large constant to prevent underflow. The recurrence equations for the matrices are given by 𝑀𝑖 𝑗 = 𝑃 (𝑟𝑖 |ℎ 𝑗 )(𝑇𝑀𝑀𝑖 𝑀𝑖 −1,𝑗 −1 + 𝑇𝐼 𝑀 𝐼𝑖 −1,𝑗 −1 + 𝑇𝐷𝑀 𝐷𝑖 −1,𝑗 −1 ), 𝐼𝑖 𝑗 = 𝑇𝑀𝐼𝑖 𝑀𝑖 −1,𝑗 + 𝑇𝐼 𝐼 𝐼𝑖 −1,𝑗 ,
(2)
𝐷𝑖 𝑗 = 𝑇𝑀𝐷𝑖 𝑀𝑖,𝑗 −1 + 𝑇𝐷𝐷 𝐷𝑖,𝑗 −1, where a read position is fixed (𝑖), while all haplotype positions are iterated (𝑗). In detail, 𝑀𝑖 𝑗 is the overall probability of subsequences 𝑟 1, ..., 𝑟𝑖 and ℎ 1, ..., ℎ 𝑗 when 𝑟𝑖 aligns to ℎ 𝑗 , while 𝐼𝑖 𝑗 and 𝐷𝑖 𝑗 are the overall probability when 𝑟𝑖 or ℎ 𝑗 aligns to a gap, respectively. 𝑇𝑀𝑀𝑖 ,𝑇𝑀𝐼𝑖 , and 𝑇𝑀𝐷𝑖 represent transition probabilities that vary with the read position, 𝑟𝑖 . Conversely, the probabilities 𝑇𝐷𝐷 ,𝑇𝐼 𝐼 ,𝑇𝐷𝑀 , and 𝑇𝐼 𝑀 are set to be constant, with 𝑇𝐷𝐷 = 𝑇𝐼 𝐼 = 0.1 and 𝑇𝐷𝑀 = 𝑇𝐼 𝑀 = 0.9. These values come from the input insertion, deletion, and gap probabilities. Lastly, the term 𝑃 (𝑟𝑖 |ℎ 𝑗 ) represents
Endeavor: Efficient PairHMM for Detection of DNA Variants in Genome-Scale Datasets
0
0
0
0
0
0
0
0
0
0
0
0
0 M Matrix
0 I Matrix
D Matrix
Figure 1: Antidiagonal Dependencies of PairHMM.
the conditional probability of the read at position 𝑖, given the haplotype at position 𝑗, calculated as ( − (𝑄𝑖 −33) 10 10 /3 if 𝑟𝑖 ≠ ℎ 𝑗 𝑃 (𝑟𝑖 |ℎ 𝑗 ) = (3) − (𝑄𝑖 −33) 1 − 10 10 if 𝑟𝑖 = ℎ 𝑗 where 𝑄𝑖 is a base quality score for the read at position 𝑖. The final result is a likelihood, 𝐿, which is the cumulative probability of all sequence alignments, calculated as ∑︁ 𝐿= 𝑀𝑙,𝑗 + 𝐼𝑙,𝑗 , (4) 𝑗
which depends only on the last rows of 𝑀 and 𝐼 . Figure 1 displays how the data dependencies of PairHMM imply the parallel processing opportunities only in antidiagonals. When calculating an element of 𝑀, 𝐼 , or 𝐷, the dependencies are on the diagonal, column, and row elements, respectively. Therefore, an anti-diagonal can be completely parallelized, using previous anti-diagonals. To accelerate the PairHMM algorithm, previous works in the literature have leveraged this property for CPU, GPU, and FPGA’s parallelization. One of the most well known implementations in CPUs is included in the Intel Genomics Kernel Library (GKL) [16], which uses AVX and AVX-512 [46] instructions for SIMD parallelization of the matrices’ anti-diagonals, as well as OpenMP to extract the most performance out of multicore processors. While this implementation is used in GATK HaplotypeCaller [13], the anti-diagonal strategy has two major drawbacks: the uncoalesced memory access pattern on the matrices and the variable computations (the first and last anti-diagonals, in the top-left and bottom-right section of the matrices, have significantly fewer elements to process than the antidiagonals in the center). Specialized FPGA-based hardware designs have also been proposed to explore anti-diagonal dependencies [11, 40, 49] and reduce the memory footprint [4, 21, 44]. Although these approaches typically provide better energy-efficiency, their limited computational resources are prohibitive for sequences with long lengths [43], which hinders their scalability. Due to their massively parallel architecture, GPUs are a prominent solution to accelerate PairHMM that also rely on the antidiagonal independence [6, 42], although this approach is not ideal, given the variance of values to calculate on each anti-diagonal and the memory access pattern. In [26], a parallelization strategy is proposed by conceptualizing a processing unit to evaluate antidiagonal cells with an adjustable number of threads. By setting the number of threads to a small number (e.g., 4), the first and last anti-diagonal do not have many idle GPU threads, improving throughput. Additionally, batch parallelization is also proposed to
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
process many sequence pairs simultaneously. However, the implementation assumes that the read lengths do not exceed 300, since it was tailored for Illumina [28] datasets. Therefore, it does not scale for long-read datasets, which have sequence lengths > 104 . Frameworks that are tailored for dynamic programming, such as CUDA Dynamic Parallelism (CDP) [1, 48], have also been explored [27]. CDP allows a CUDA kernel (parent) to launch new kernels (children) without returning control to the host, allowing for nested parallelism. This programming model was used to accelerate the PairHMM algorithm [27], but the kernel launch overhead, as well as the hardware limits on the number of concurrent kernel launches and the depth of nested kernels [10], significantly impact the scalability of CDP for long-read datasets. In the current literature, a novel approach herein referred as gpuPairHMM [45] is the fastest PairHMM implementation, leveraging the anti-diagonal independence and processing read-haplotype pairs at the level of a warp / subwarp by having individual threads calculate multiple elements of 𝑀, 𝐼 , and 𝐷 to achieve unprecedented throughput for PairHMM. While these solutions achieve good results on short-read genomic datasets, the current solutions are vendor-specific and based on a non-ideal parallelism pattern. Therefore, a breakthrough in performance is expected if a PairHMM approach extracts more parallelism from the algorithm while keeping a high throughput and scaling to long reads to process any real genomic dataset in different devices. To this end, Endeavor proposes a row-wise parallelism strategy as a possible solution, which has been applied to the sequence alignment problem [24], but has never been studied for PairHMM. An ideal row-wise definition would do one pass per row, but the existing mathematical proofs for sequence alignment require two passes per row [2] and the key operations are different from PairHMM (maximums instead of sums and products). Although both algorithms share similarities, PairHMM brings its own unique parallelization challenges, which Endeavor addresses in the CPU and GPU methods proposed in this work (Section 3) to deliver high-performance solutions for different hardware architectures to overcome the limitations of the state-of-the-art.
3
Proposed Hardware Acceleration for Variant Calling
In this section, the proposed CPU and GPU-based approaches to accelerate PairHMM is described. In particular, a mathematical derivation starting from the algorithm’s analytical formulation is provided, as it is an essential step to understand how we extract more fine-grained data parallelism. It is worth noting that the proposed algorithm is a novel approach to PairHMM parallel processing, which overcomes the limitations of existing state-ofthe-art methods that only exploit parallelism on anti-diagonals (see Section 2). Afterwards, we elaborate on how the proposed approach is employed to define the proposed methods to process short-read datasets (sequences with a range of 250-800 basepairs), as well as long reads (sequences with > 104 − 105 basepairs).
3.1
Algorithm Design
As mentioned in Section 2 (see Equation 2), there are three main expressions that characterize the PairHMM algorithm and the dependencies between the 𝑀, 𝐼 , and 𝐷 matrices. One aspect worth
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
Miguel Graça and Aleksandar Ilic
noticing about PairHMM is that the final result (see 4) does not require the 𝐷 matrix and is only dependent on 𝑀 and 𝐼 . However, 𝐷 is still necessary for the intermediate calculation of the elements of 𝑀 (see 𝐷𝑖 −1,𝑗 −1 in Equation 2). Based on this observation, we start by expressing 𝐷𝑖 −1,𝑗 −1 as a function of 𝑀 and 𝐼 (as stated in Equation 2) such that 𝐷𝑖 −1,𝑗 −1 = 𝑇𝑀𝐷𝑖 −1 𝑀𝑖 −1,𝑗 −2 + 𝑇𝐷𝐷 𝐷𝑖 −1,𝑗 −2 .
(5)
By relying on the same expression, the first four terms (i.e., 𝑗 ∈ {1, ..., 4}) of the 𝑖 − 1 row of 𝐷 can be expressed as
0
0
0
0
0
0
0
0
0
0 M Matrix
I Matrix
𝐷𝑖 −1,1 = 𝑇𝑀𝐷𝑖 −1 𝑀𝑖 −1,0 + 𝑇𝐷𝐷 𝐷𝑖 −1,0 = 0 Figure 2: Rowwise Dependencies of Endeavor.
𝐷𝑖 −1,2 = 𝑇𝑀𝐷𝑖 −1 𝑀𝑖 −1,1 + 𝑇𝐷𝐷 𝐷𝑖 −1,1 = 𝑇𝑀𝐷𝑖 −1 𝑀𝑖 −1,1 𝐷𝑖 −1,3 = 𝑇𝑀𝐷𝑖 −1 𝑀𝑖 −1,2 + 𝑇𝐷𝐷 𝐷𝑖 −1,2 = 𝑇𝑀𝐷𝑖 −1 𝑀𝑖 −1,2 + 𝑇𝐷𝐷 𝑇𝑀𝐷𝑖 −1 𝑀𝑖 −1,1 = 𝑇𝑀𝐷𝑖 −1 (𝑀𝑖 −1,2 + 𝑇𝐷𝐷 𝑀𝑖 −1,1 ) = 𝑇𝑀𝐷𝑖 −1
2 ∑︁
2−𝑏 𝑀𝑖 −1,𝑏𝑇𝐷𝐷
𝑏=1
𝐷𝑖 −1,4 = 𝑇𝑀𝐷𝑖 −1 𝑀𝑖 −1,3 + 𝑇𝐷𝐷 𝐷𝑖 −1,3 2 = 𝑇𝑀𝐷𝑖 −1 (𝑀𝑖 −1,3 + 𝑇𝐷𝐷 𝑀𝑖 −1,2 + 𝑇𝐷𝐷 𝑀𝑖 −1,1 )
= 𝑇𝑀𝐷𝑖 −1
3 ∑︁
3−𝑏 𝑀𝑖 −1,𝑏𝑇𝐷𝐷 ,
𝑏=1
where 𝐷𝑖 −1,1 = 0 as 𝑀𝑖 −1,0 = 𝐷𝑖 −1,0 = 0 (see Equation 1). If the unrolling procedure is applied to each new term of 𝐷, one obtains a general sum of values in 𝑀𝑖 −1,𝑏 that can be written as follows 𝑗 −3 𝐷𝑖 −1,𝑗 −1 = 𝑇𝑀𝐷𝑖 −1 (𝑀𝑖 −1,𝑗 −2 + 𝑇𝐷𝐷 𝑀𝑖 −1,𝑗 −3 + ... + 𝑇𝐷𝐷 𝑀𝑖 −1,1 )
= 𝑇𝑀𝐷𝑖 −1
𝑗 −2 ∑︁
𝑗 −2−𝑏 𝑀𝑖 −1,𝑏𝑇𝐷𝐷 .
(6)
𝑏=1
As can be observed, Equation 6 shows that 𝐷𝑖 −1,𝑗 −1 depends only on previous values of 𝑀𝑖 −1,𝑏 . Therefore, both 𝑀𝑖 𝑗 and 𝐼𝑖 𝑗 (see Equation 2) depend only on values from the previous row, 𝑀𝑖 −1,𝑏 and 𝐼𝑖 −1,𝑏 . As such, all the elements in a row can be calculated in parallel. In this case, the equations for PairHMM can be expressed as follows 𝑀𝑖 𝑗 = 𝑃 (𝑟𝑖 |ℎ 𝑗 )(𝑇𝑀𝑀𝑖 𝑀𝑖 −1,𝑗 −1 + 𝑇𝐼 𝑀 𝐼𝑖 −1,𝑗 −1 + 𝑇𝐷𝑀 𝑇𝑀𝐷𝑖 −1
𝑗 −2 ∑︁
𝑗 −2−𝑏 𝑀𝑖 −1,𝑏𝑇𝐷𝐷 )
Figure 3: Multithreading+SIMD-Based (CPU) and WarpBased (GPU) Implementations of PairHMM.
(7)
𝑏=1
𝐼𝑖 𝑗 = 𝑇𝑀𝐼𝑖 𝑀𝑖 −1,𝑗 + 𝑇𝐼 𝐼 𝐼𝑖 −1,𝑗 , substituting Equation 2. Defining PairHMM in this fashion provides two major benefits, as is shown in Figure 2. First, 𝐷 is no longer necessary, which reduces the storage requirements for PairHMM. Second, removing 𝐷 eliminates one of the dependencies, which improves parallelism and performance. For a read of length 𝑙 and a haplotype of length 𝑘, the anti-diagonal implementation calculates 𝑀, 𝐼 , and 𝐷, requiring O (𝑙 + 𝑘 − 1) passes (since there are 𝑙 + 𝑘 − 1 antidiagonals). In comparison, the proposed approach calculates only 𝑀 and 𝐼 (which are the necessary matrices for the likelihood), requiring only O (𝑙) passes to do all calculations. Furthermore, the proposed algorithm outputs the same results as GATK HaplotypeCaller [13], without any loss in accuracy.
3.2
Endeavor Approach Overview
To map the redefined PairHMM algorithm to multicore CPUs and highly parallel GPU hardware, a high level overview of Endeavor is provided herein and described in Figure 3. At the CPU level, readhaplotype pairs are read from a dataset and processed at the level of a single SIMD lane (e.g., in Figure 3, using double precision (64-bit) and a SIMD width of 256 bits allows a single thread to process 4 (256/64) read-haplotypes in parallel). As read-haplotype pairs are independent, the use of multiple threads, each one operating on SIMD lanes, allows to maximize the algorithm’s throughput. At the GPU level, the host reads read-haplotype pairs from a dataset, which are sent to the GPU. After performing the PairHMM computations, the GPU sends the likelihoods back to the host. Endeavor’s approach aims at minimizing synchronization in GPUs by
Endeavor: Efficient PairHMM for Detection of DNA Variants in Genome-Scale Datasets
processing the sequences at the level of a warp of size 𝑊 (where 𝑊 = 32 if the GPU is from NVIDIA or Intel and 𝑊 = 64 if the GPU is from AMD to match the warp size on each architecture). All threads iterate over the read’s characters to calculate the rows of 𝑀 and 𝐼 , while each individual thread reads 𝑁 characters from the haplotype and calculates 𝑁 elements of a 𝑀 and 𝐼 row (e.g., in Figure 3, each thread is responsible for calculating 2 elements of a row). After calculating the last row, the values of 𝑀 and 𝐼 are accumulated across all threads to obtain the final likelihood, which is written to global memory and sent back to the host. To achieve high throughput, Endeavor focuses on several different aspects. At the GPU level, read-haplotype sequence pairs can also be processed at the level of a sub-warp of size 𝑊 /𝑃. As an example, if a single thread computes 𝑁 elements of a row of 𝑀 and 𝐼 and haplotypes of size 2𝑁 are to be analyzed, then only 2 threads are necessary for the computations, which allows 16 read-haplotype pairs to be processed if 𝑊 = 32 and 𝑃 = 16 or 32 read-haplotype pairs if 𝑊 = 64 and 𝑃 = 32. It is worth noting that all threads within a warp will still adhere to SIMD processing with no divergency (i.e., performing the same operations on different data). To optimize computations for different sequence lengths, read-haplotype pairs are batched according to their lengths and run on dedicated kernels according to the value of 𝑃, which defines the warp/subwarp size (𝑊 /𝑃). To maximize performance, the value of 𝑁 must be carefully selected to ensure that each thread has enough computations to perform and to avoid register spilling to global memory, which would decrease throughput. In Section 4.2, an analysis of the impact of 𝑁 in different hardware architectures is performed. At the CPU level, the SIMD width impacts the read-haplotype pairs that each thread can process in parallel. Endeavor explores the impact in performance of two different SIMD extensions: AVX (256 bits) and AVX-512 (512 bits). Finally, for both devices, as the computations for PairHMM revolve around probabilities, precision impacts the algorithm’s accuracy and maximum throughput. To this end, Endeavor implements PairHMM in double and single precision. In Section 4.3, Endeavor’s peak throughput is analyzed in these two precisions for various CPU and GPU architectures.
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
Listing 1: CPU kernel for AVX-level PairHMM (FP32). 1 2 3 4 5 6 7
template < unsigned int K > void avx_fp32 ( int * hap , int * read , int * _q , int * _i , int * _d , int seqsize , float * sum , float * q32 ) { # pragma omp parallel { int tid = omp_get_thread_num () ; int nthr = omp_get_num_threads () ;
8
for ( int seq = tid ; seq < seqsize /8; seq += nthr ){ m256 Mp [ K ] , Ip [ K ] , _hap [ K ]; m256 M , y2 , y3 , one , thr , d0 , tmm , pr , qual , tmi , mask ; m256 one0 = m256_set_ps (0.1 f ) , nine = m256_set_ps (0.9 f ) , finalsum = m256_set_ps (0.0 f ) ; m256i rs , q0 , in , d1 ;
9 10 11 12 13 14
rs = m256_load ( read ) ; q0 = m256_load ( q ) ; d1 = m256_load ( d ) ; qual = m256_gather_ps ( q32 , q0 , 4) ; thr = m256_mul_ps ( m256_set_ps (1/3) , qual ) ; one = m256_sub_ps ( m256_set_ps (1) , qual ) ;
15 16 17 18 19 20 21
for ( int h = 0; h < K ; h ++) { _hap [ h ] = m256_load ( hap + h *8) mask = m256_cast_ps ( m256_cmpeq ( rs , _hap [h ]) ); pr = m256_blend_ps ( thr , one , mask ) ; Mp [ h ] = m256_mul_ps ( pr , m256_mul_ps ( nine , m256_div_ps ( m256_set_ps ( ldexpf (1. f , 120. f )) , m256_set_ps ( N ) ) ) ) ; Ip [ h ] = m256_set_ps (0) ; }
22 23 24 25 26
27 28 29
for ( int k = 1; k < K ; k ++) { rs = m256_load ( read + k *8) ; q0 = m256_load ( q + k *8) ; in = m256_load ( i + k *8) ; d0 = m256_mul_ps ( nine , m256_gather_ps ( q32 , d1 , 4) ) ; d1 = m256_load ( d + k *8) ;
30 31 32 33 34 35 36
mask = m256_cast_ps ( m256_cmpeq ( rs , _hap [0]) ); qual = m256_gather_ps ( q32 , q0 , 4) ; tmi = m256_gather_ps ( q32 , in , 4) ; thr = m256_mul_ps ( m256_set_ps (1/3) , qual ); one = m256_sub_ps ( m256_set_ps (1) , qual ); pr = m256_blend_ps ( thr , one , mask ) ; tmm = m256_sub_ps ( m256_sub_ps ( m256_set_ps (1) , m256_gather_ps ( q32 , d1 , 4) ) , tmi ) ; y3 = M = m256_set_ps (0) ;
37 38 39 40 41 42 43 44
3.3
CPU SIMD Implementation
Listing 1 provides an example of the proposed CPU implementation to process read-haplotype pairs of size 𝐾 in single precision, using AVX intrinsics, as described in Figure 3 (note that, since AVX performs computations with 256 bits and single precision values (32-bit) are used, this amounts to processing 8 (256/32) sequence pairs per thread). The main inputs are the read and haplotype sequences, the quality (_𝑞), insertion (_𝑖), and deletion (_𝑑) scores (described as characters, which are associated to powers of 10, as demonstrated in Equation 3), an array with powers of 10 (𝑞32), and the total number of sequence pairs to process (𝑠𝑒𝑞𝑠𝑖𝑧𝑒). The output is an array of sums (each position will hold the sum of the last rows of 𝑀 and 𝐼 for a given read-haplotype pair). For this kernel, it is assumed that all the evaluated sequences have length 𝐾. In line 4, the kernel spawns multiple threads to perform parallel work. Lines 6 to 14 initialize the necessary variables with which each thread will work, namely the 𝑀 and 𝐼 rows (𝑀𝑝, 𝐼𝑝) and haplotype characters (_ℎ𝑎𝑝) in line 10. Note that 𝑀𝑝 and 𝐼𝑝 are
45
for ( int h = 1; h < K ; h ++) { mask = m256_cast_ps ( m256_cmpeq ( rs , _hap [h ]) ); pr = m256_blend_ps ( thr , one , mask ) ; y2 = m256_mul_ps ( pr , m256_add_ps ( m256_add_ps ( m256_mul_ps ( Mp [h -1] , tmm ) , m256_mul_ps ( nine , Ip [h -1]) ) , m256_mul_ps ( d0 , y3 ) ) ) ; y3 = m256_fma_ps ( y3 , one0 , Mp [h -1]) ; Ipr [h -1] = m256_fma_ps ( Mp [h -1] , tmi , m256_mul_ps ( Ip [h -1] , one0 ) ) ; Mp [h -1] = M ; M = y2 ; }
46 47 48 49
50 51 52 53 54 55
Ip [N -1] = m256_fma_ps ( Mp [N -1] , tmi , m256_mul_ps ( Ip [N -1] , one0 ) ) ; Mp [N -1] = M ; }
56 57 58
for ( int h = 0; h < K ; h ++) finalsum = m256_add_ps ( finalsum , m256_add_ps ( Mp [ h ] , Ip [ h ]) ) ;
59 60
m256_store_ps (( sum + seq ) , finalsum ) ;
61
}
62
}
63 64
}
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
Miguel Graça and Aleksandar Ilic
Thread 0 Thread 1 Row M1,1 M1,2 M1,3 M1,4 M1,5 M1,6 M1,7 M1,8 1
Algorithm 1 Warp Inclusive Scan 1: Input: T (value of thread), tid (thread ID) 2: Output: T (result of prefix sum for each thread) 3:
xTDD
V=0 V=
xTDD
4
Ʃ
b=1
xTDD
xTDD
V=0 (TDD4-b) X (M1, b)
V=
xTDD
4
Ʃ
b=1
xTDD
T=
Inclusive Scan (using warp shuffle intrinsics)
(TDD4-b) X (M1, b+4)
TDD-1 x
line 56. In line 59, the results of the last row of 𝑀 and 𝐼 are summed and stored in memory in line 61.
3.4 T=
T
T *= TDD T *= TDD
TDD-1 x
+ M1,1
(T - M1,4)
+ M1,2
4
Ʃ (T b=1 T
DD
4-b
) X (Mi, b)
T *= TDD T *= TDD + M1,5
+ M1,6
Row M2,1 M2,2 M2,3 M2,4 M2,5 M2,6 M2,7 M2,8 2 xTDD
V=0
6:
(TDD4-b) X (M1, b)
xTDD4
T=0
(T - M1,0)
4
Ʃ b=1
_𝑎𝑢𝑥 =__shfl_inst(−1,𝑇 , 𝑜 𝑓 𝑓 𝑠𝑒𝑡); 𝑜 𝑓 𝑓 𝑠𝑒𝑡 ∗𝑁 if 𝑡𝑖𝑑%𝑊 ≥ 𝑜 𝑓 𝑓 𝑠𝑒𝑡: 𝑇 + = _𝑎𝑢𝑥 × 𝑇𝐷𝐷 ; 7: end for 8: return 𝑇 5:
T = _shfl_inst(V, 1) (Thread 0 will have T = 0) T=0
4: for (𝑜 𝑓 𝑓 𝑠𝑒𝑡 = 1; 𝑜 𝑓 𝑓 𝑠𝑒𝑡 < 𝑊 ; 𝑜 𝑓 𝑓 𝑠𝑒𝑡 × = 2) do
xTDD
xTDD
xTDD
xTDD
xTDD
V=0
Figure 4: Endeavor’s GPU Pipeline for 𝑀 (for this example, 𝑁 = 4, with 2 threads per read-haplotype pair). defined as 𝑚256 arrays of size 𝐾 where, for example, the first 32 bits in each position define the 𝑀 and 𝐼 rows for the first read-haplotype pair to process. The first read and quality characters are read in lines 15 to 17. In line 18, a gather operation (𝑚256_𝑔𝑎𝑡ℎ𝑒𝑟 _𝑝𝑠) loads from memory the powers of 10 associated with the quality scores (since different read-haplotype pairs may have different scores and, therefore, require different powers of 10, which will be at different positions in the array). In lines 23 to 27, the haplotype characters are read from memory (𝑚256_𝑙𝑜𝑎𝑑) and the first row of 𝑀 (line 26) and 𝐼 (line 27) is calculated. When the read and haplotype characters are compared (line 24), a mask is created (𝑚256_𝑐𝑚𝑝𝑒𝑞) to select between the two possible values (line 25) that 𝑃 (𝑟𝑖 |ℎ 𝑗 ) can assume (see Equation 3). These two values are calculated on lines 20 and 21 to avoid recalculating 𝑃 (𝑟𝑖 |ℎ 𝑗 ) for each haplotype character, since they depend only on the quality score of the current read character. In lines 30 to 57, the remaining rows are calculated. At each iteration of the outer loop, the next read characters and quality scores are loaded from memory (lines 31 to 35). Next, the new values of 𝑇𝑀𝐼 (line 39) and 𝑇𝑀𝑀 (line 43), as well as the possible values for 𝑃 (𝑟𝑖 |ℎ 𝑗 ) (lines 40 and 41), are calculated. The first element of a 𝑀 row is set to 0 (line 44) due to the initialization in Equation 1 (𝑀1𝑗 depends on 𝑀0,𝑗 −1 and 𝐼 0,𝑗 −1 , which are set to 0). The remaining elements are calculated in the inner loop (lines 46 to 54) according to Equation 7, with the last elements of 𝑀 and 𝐼 being updated in
GPU Cross-Thread Intrinsics
In Section 3.2, Endeavor’s GPU approach was defined to process a read-haplotype pair at the warp-level by leveraging row-wise parallelism. Unlike CPUs, GPU threads cannot efficiently compute a complete row of 𝑀 and 𝐼 for long sequences, as keeping the results in private memory would lead to register spilling and lower performance. Therefore, each thread calculates 𝑁 elements of a row. However, note that the row-wise definition of PairHMM implies that not all elements of 𝑀 require the same amount of computation, which will hinder performance due to poor load balancing. For example, 𝑀𝑖 𝑗 needs 𝑗 − 2 elements to compute the sum derived in Equation 6, 𝑀𝑖,𝑗 −1 , requires 𝑗 − 3 elements, i.e., each subsequent 𝑀𝑖 𝑗 will require one less element. To overcome load balancing issues between different threads within the same warp, a warp shuffle strategy is proposed to ensure the correctness of 𝑀 and illustrated in Figure 4 with an example using 2 threads, each calculating 4 elements of 𝑀 (i.e., 𝑁 = 4). The first two elements (𝑀𝑖,1 and 𝑀𝑖,2 ) do not compute the sum from Equation 6, as it would start from 𝑀𝑖 −1,−1 (which does not exist) and 𝑀𝑖 −1,0 (which is zero), respectively. 𝑀𝑖,3 needs 𝑀𝑖 −1,1 for the sum and the next element, 𝑀𝑖,4 , needs 𝑀𝑖 −1,2 +𝑇𝐷𝐷 𝑀𝑖 −1,1 , where the second term is the sum from the previous calculated element of 𝑀 (𝑀𝑖,3 ) scaled by 𝑇𝐷𝐷 . In general, 𝑀𝑖,𝑘 scales the sum of the previous element, 𝑀𝑖,𝑘 −1 , by 𝑇𝐷𝐷 and adds another element to the sum. While thread 0, which calculates the first 𝑁 elements of 𝑀, can apply the consecutive additions and scaling operations and correctly calculate 𝑀, the same cannot be said for thread 1, which needs the sum from thread 0. In general, thread 𝑇𝐼 𝐷 needs the sums from threads 0 to 𝑇𝐼 𝐷 − 1. The structure of the terms for each thread to calculate 𝑀 according to Equation 6 resembles the inclusive scan [19], apart from 𝑇𝐷𝐷 factors. The first row of 𝑀 (𝑀1,𝑗 ) is easily calculated due to the algorithm’s initialization (because 𝑀0,𝑗 = 𝐼 0,𝑗 = 0 and 𝐷 0,𝑗 = 21020 /|𝐻 |, 𝑀1,𝑗 = 𝑝 1,𝑗 × 𝑇𝐷𝑀 × 21020 /𝑘). As each element is calculated, each thread accumulates the scaled sum (𝑉 ). The following step makes use of intrinsic shuffle operations (__𝑠ℎ𝑓 𝑙_𝑖𝑛𝑠𝑡 ()) so that thread 𝑇𝐼 𝐷 receives 𝑉 from thread 𝑇𝐼 𝐷 − 1 and stores the result in 𝑇 (for CUDA and HIP, this instruction is equivalent to __𝑠ℎ𝑓 𝑙_𝑢𝑝_𝑠𝑦𝑛𝑐 (), and __𝑠ℎ𝑓 𝑙_𝑢𝑝 (), respectively). For thread 0, 𝑇 is set to 0. After setting 𝑇 , warp shuffle operations are executed to perform an inclusive scan, as described in Algorithm 1, in a warp or subwarp of size 𝑊 . In Figure 4, 𝑊 = 2 and only one shuffle instruction is necessary,
Endeavor: Efficient PairHMM for Detection of DNA Variants in Genome-Scale Datasets
where thread 𝑇𝐼 𝐷 receives 𝑇 from thread 𝑇𝐼 𝐷 − 1 and threads with 𝑇𝐼 𝐷 ≥ 1 add it to its own value of 𝑇 by scaling it to a factor of 𝑁 (where 𝑁 = 4 in the case of Figure 4). In general, for a warp 𝑇𝐷𝐷 or subwarp of size 𝑊 , 𝑙𝑜𝑔2 (𝑊 ) shuffle operations are necessary, 𝑁 , as depicted in Algorithm 1. After these scaled by powers of 𝑇𝐷𝐷 operations, each thread has the correct sum for the next row of 𝑀. For example, for thread 1, 𝑇 is given by 𝑇 =
8 ∑︁
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
Listing 2: GPU kernel for warp-level PairHMM (FP64). 1 2
3 4 5 6
8−𝑏 𝑀1,𝑏𝑇𝐷𝐷
(8)
7 8 9
𝑏=1
10
and the sums that thread 1 needs to calculate its elements (𝑀2,5 to 𝑀2,8 ) are given by 𝑀2,9 → 𝑀2,10 → 𝑀2,11 → 𝑀2,12 →
7 ∑︁ 𝑏=1 8 ∑︁ 𝑏=1 9 ∑︁ 𝑏=1 10 ∑︁
template < unsigned int N , unsigned int P > __global__ void fp64_subwarp ( char * hap , char * _rs , char * _q , char * _i , char * _d , double * q64 , double * finalsum ) { short subwarp = 32/ P ; short laneID = threadIdx . x % subwarp , warpID = threadIdx . x / subwarp ; int seqidx = P * blockIdx . x + warpID ; double Mp [ N ] = {0} , Ip [ N ] = {0}; char _hap [ N ] = {0}; double M =0 , I =0 , y1 =0 , y2 =0 , T =0 , V =0 , d0 =0 , tmm =0; char rs , q , in , d1 ;
11
rs = _rs [ tidx ] , q = _q [ tidx ]; for ( unsigned int h = 0; h < N ; h ++) { _hap [ h ] = hap [32* h + tidx ]; Mp [ h ] = ( _hap [ h ] == rs ) ? (1 - q64 [ q ]) : q64 [q ]*0.33; Mp [ h ] *= 0.9* ldexp (1 ,1020) / N ; V = V *0.1 + Mprev [ h ]; } T = __shfl_inst ( -1 , V , 1) ; if ( laneID == 0) T = 0; T = warp_inclusive_scan (T , laneID ) ; for ( int k = 1; k < N ; k ++) { rs = _rs [32* k + tidx ] , q = _q [32* k + tidx ]; in = _i [32* k + tidx ] , d0 = q64 [ _d [32*( k -1) + tidx ]]; d1 = _d [32* k + tidx ]; tmm = 1 - q64 [ in ] - q64 [ d1 ]; V = __shfl_inst ( -1 , Mp [N -1] , 1) ; y1 = __shfl_inst ( -1 , Ip [N -1] , 1) ; y2 = ( _hap [0] == rs ) ? (1 - q64 [ q ]) : q64 [q ]*0.33; y2 *= ( V * tmm + 0.9*( y1 + d0 *10*( T - V ) ) ) ; if ( laneID == 0) y2 = 0; M = y2 , V = y2 ; I = Mp [0]* q64 [ in ] + Ip [0]*0.1; for ( int h = 1; h < N ; h ++) { y2 = ( _hap [ h ] == rs ) ? (1 - q64 [ q ]) : q64 [q ]*0.33; y2 *= ( Mp [h -1]* tmm + 0.9*( Ip [h -1] + d0 * T )); T = T *0.1 + Mp [h -1]; Mp [h -1] = M , Ip [h -1] = I ; M = y2 , I = Mp [ h ]* q64 [ in ] + Ip [ h ]*0.1; V = V *0.1 + y2 ; } Mp [N -1] = M , Ip [N -1] = I ; T = __shfl_inst ( -1 , V , 1) ; if ( laneID == 0) T = 0; T = warp_inclusive_scan (T , laneID ) ; } y1 = 0; for ( int h = 0; h < N ; h ++) y1 += Mp [ h ] + Ip [ h ]; y1 = warp_inclusive_scan ( y1 , laneID ) ; if ( laneID == subwarp - 1) finalsum [ seqidx ] = y1 ;
12 13 14
7−𝑏 −1 𝑀1,𝑏𝑇𝐷𝐷 = 𝑇𝐷𝐷 (𝑇 − 𝑀1,8 )
15 16 17 18
8−𝑏 𝑀1,𝑏𝑇𝐷𝐷 =𝑇
19 20
(9) 9−𝑏 𝑀1,𝑏𝑇𝐷𝐷 = 𝑇𝐷𝐷 𝑇 + 𝑀1,9
21 22 23 24 25 26
10−𝑏 2 𝑀1,𝑏𝑇𝐷𝐷 = 𝑇𝐷𝐷 𝑇 + 𝑇𝐷𝐷 𝑀1,9 + 𝑀1,10
27 28
𝑏=1 29
which, as is shown, can be derived from 𝑇 . For each thread, the first element of 𝑀 is done separately, as it requires values from other threads (in this example, thread 1 needs to retrieve 𝑀1,4 from thread 0 to calculate 𝑀2,5 ), which can be obtained through shuffle intrinsics. Afterwards, the remaining elements can be calculated in a loop, scaling 𝑇 with 𝑇𝐷𝐷 while adding the necessary elements from the previous row and scaling 𝑉 with 𝑇𝐷𝐷 while adding the necessary elements from the current row. As 𝑇𝐷𝐷 is constant (0.1), all its powers are also constant and stored in constant memory. Listing 2 provides the developed GPU kernel code that implements Endeavor’s row-wise PairHMM algorithm. The overall structure is similar to Listing 1, with 𝑀𝑝 and 𝐼𝑝 being defined as arrays of size 𝑁 for each thread, the addition of the warp inclusive scan (lines 21 and 45), which corresponds to the "Inclusive Scan" section in Figure 4), the shuffle intrinsics that allow thread 𝑇𝐼 𝐷 to retrieve the last values for 𝑀 and 𝐼 from thread 𝑇𝐼 𝐷 − 1 (lines 27 and 28), and the shuffle intrinsics to exchange 𝑉 , which has the accumulated sum of 𝑀 elements of each thread (lines 19 and 43). The approach herein devised works for haplotypes of size up to 𝑁 × 𝑊 (where 𝑁 is the number of elements a single thread processes and 𝑊 the warp size). For bigger haplotypes, Endeavor can be scaled to use multiple warps. To achieve this, the last thread of each warp saves in shared memory its last element of 𝑀 and 𝐼 (to be used by the first thread of the next warp to calculate its first element of 𝑀), and the sum of all 𝑀 elements calculated by the warp, given by 𝑇 =
𝑁∑︁ ×𝑊
𝑁 ×𝑊 −𝑏 𝑀𝑖,𝑏𝑇𝐷𝐷
(10)
𝑏=1
as this result will be needed by the next warp to calculate its 𝑀 elements. To ensure correctness when reading from/writing to shared memory, thread synchronization is added.
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
}
Setting 𝑀𝑇 as the maximum number of threads in a block, this approach works for haplotypes up to 𝑁 × 𝑀𝑇 . To process even longer read-haplotype pairs, one can use distributed shared memory to exchange values between different thread blocks within a cluster. To achieve this, the last thread of each block saves in distributed shared memory its last element of 𝑀 and 𝐼 (to be used by the first thread of the next block to calculate its first element of 𝑀), and the sum of all 𝑀 elements calculated by the block, given by 𝑇 =
𝑁∑︁ ×𝑀𝑇
𝑁 ×𝑀𝑇 −𝑏 𝑀𝑖,𝑏𝑇𝐷𝐷
(11)
𝑏=1
as this result will be needed by the next warp to calculate its 𝑀 elements. To ensure correctness when reading from/writing to shared memory, inter-block synchronization is added.
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
Miguel Graça and Aleksandar Ilic
(a) H100
(b) MI300A
Figure 5: TCUPS evolution with the elements processed by each thread (𝑁 ) in H100 (left) and MI300A (right) in FP32 and FP64.
4
Experimental Results
In this section, details on the experimental evaluation of the Endeavor framework are provided. First, we focus on Endeavor’s scalability and portability to answer the following questions: (i) how does the number of elements that each thread calculates (𝑁 ) impact Endeavor’s throughput on GPUs?, and (ii) how does the throughput scale with sequence length on different CPUs and GPUs? Second, Endeavor is tested on large-scale genomic datasets from different sequencing technologies, with varying sequence lengths and sequence pairs to be processed. In doing so, we answer the following questions: (i) how does Endeavor compare with a highly-optimized GATK HaplotypeCaller PairHMM CPU-based method?, and (ii) how does Endeavor compare with state-of-the-art GPU-based methods?
4.1
Experimental Setup
The reported experiments have been performed on 5 GPUs from different major vendors, i.e., NVIDIA V100 (300 W), A100 (300 W), H100 (700 W), and RTX 6000 PRO (600 W) GPUs, using the CUDA 12.8 toolkit, and AMD MI300A (750 W), using ROCM 6.3.4, as well as 2 CPU systems, with AMD EPYC Zen 4 (which is included in MI300A) and Intel Xeon Gold 6438 architectures, both with 96 cores and 192 available threads. To fully examine the capabilities of the proposed approach in datasets generated with short-read and long-read sequencing technologies, datasets from the GIAB [50] consortium are used, namely the NA12878 [32, 35–37], NA24695 [34], NA24631 [38], and NA24149 [33] datasets. For the GIAB datasets, Endeavor is tested on all GPUs and compared with gpuPairHMM [45] (as it is the current best PairHMM implementation for GPU) and GATK HaplotypeCaller, which uses the GKL [16] implementation for PairHMM and is capable of processing long-read datasets, on the Intel Xeon Gold 6438 CPU (with AVX-512 instructions) using all available threads (192). The 10s dataset [9] is also used for benchmarking and comparison with other approaches in the literature. The evaluation of Endeavor and comparison to other state-ofthe-art solutions focuses on throughput, measured as Tera Cell Updates Per Second (TCUPS) and calculated as Í𝑠 (𝑟𝑙𝑒𝑛𝑖 × ℎ𝑙𝑒𝑛𝑖 ) 𝑇𝐶𝑈 𝑃𝑆 = 𝑖 , 𝑡 × 1012
where 𝑡 is the runtime in seconds, 𝑠 is the total number of readhaplotype pairs to process, 𝑟𝑙𝑒𝑛𝑖 and ℎ𝑙𝑒𝑛𝑖 are the length of the 𝑖-th read and haplotype, respectively. For benchmarking, varying sequence lengths are not considered (i.e., all read-haplotype pairs are the same length). For the GPU approach, for fair comparison with the state-of-the-art [45], in Section 4.3, only the kernel runtime is considered, excluding data transfer overheads, while in Section 4.4 the data transfer time is also considered in the execution times.
4.2
GPU Kernel Parameters Analysis
In this section, an analysis on the impact of 𝑁 (the number of elements of 𝑀 and 𝐼 processed per thread) in GPU performance is performed. To this end, experiments on H100 and MI300A are conducted, fixing the read-haplotype size to 512 (this is in the middle of the short-read range, which typically spans 250-800 basepairs [20]) and varying the value of 𝑁 , using the warp size as a reference. Among the tested NVIDIA GPUs, H100 is the best choice to analyze both precisions (while RTX 6000 PRO has higher FP32 peak throughput, the FP64 peak throughput is 1/64 of the single precision, instead of 1/2 as is the case of H100). Figures 5a and 5b show the variations in TCUPS for H100 and MI300A, respectively, on FP64 and FP32, according to 𝑁 . H100 has a warp size of 32 (with 1/32 representing a single thread), while MI300A has a warp size of 64 (with 1/64 representing a single thread). If 𝑁 is too small, the number of threads per block has to be increased, with each thread performing little work, resulting in decreased throughput. On the other hand, if 𝑁 is too big (in the limit, if a single thread calculates all values of 𝑀 and 𝐼 for a single read-haplotype pair), the register pressure results in spilling to the global memory, which, due to its high latency, reduces the kernel’s throughput. Therefore, one must achieve a balance between providing enough elements for threads to process and ensuring that register spilling does not occur. For H100, the optimal value is 𝑁 = 32 (i.e., for the considered length of 512, 16 threads, which is 1/2 of a warp, calculate 32 elements of 𝑀 and 𝐼 rows), followed by 𝑁 = 64, while MI300A achieves the best TCUPS for 𝑁 = 16. The results show that the optimal value of 𝑁 is hardware dependent and can be tuned to achieve maximum performance on different GPU architectures. For this work, V100, A100, H100, and RTX 6000 PRO
Endeavor: Efficient PairHMM for Detection of DNA Variants in Genome-Scale Datasets
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
Figure 6: TCUPS evolution in Intel Xeon Gold 6438 for Endeavor and GKL.
Figure 8: TCUPS evolution in V100, A100, H100, and RTX 6000 for Endeavor and gpuPairHMM.
Figure 7: TCUPS evolution in AMD EPYC Zen 4 for Endeavor and GKL. use 𝑁 = 32 for sequence lengths up to 1024 and 𝑁 = 64 for longer sequences, while MI300A uses 𝑁 = 16 for all sequence lengths.
4.3
Scalability Analysis
To evaluate Endeavor’s scalability, synthetic datasets with sequence lengths ranging from 32 to 131072 basepairs are generated, with the read and haplotype having the same length. In doing so, we focus on answering the following questions: (i) how does Endeavor’s performance scale on CPUs and compares to the state-of-the-art libraries?, and (ii) how does Endeavor perform on GPUs from different vendors and compares to the state-of-the-art approaches? 4.3.1 CPU Analysis. Figures 6 and 7 display the TCUPS evolution in the Intel and AMD CPUs, respectively, as the input readhaplotype length increases. Eight bars are displayed, four for Endeavor combining FP32/FP64 precisions and AVX/AVX-512 SIMD widths, and four for the CPU state-of-the-art approach, GKL [16], for the same precisions and SIMD widths. The results demonstrate that Endeavor outperforms GKL for all tested combinations of SIMD width and floating-point precisions, with the AMD EPYC Zen 4 achieving a peak throughput of 0.625 TCUPS using AVX with FP32 and the Intel Xeon Gold 6438 achieving the best results for AVX-512 and FP32, with 0.415 TCUPS. In contrast, GKL only achieves a peak throughput of 0.323 TCUPS and 0.194 TCUPS, respectively, using AVX-512 with FP32 in both cases. Therefore, Endeavor achieves a
1.93x and 2.14x improvement in PairHMM throughput. Note that Zen 4 does not achieve optimal throughput using AVX-512 as the operations run on a 256-bit data path that requires two passes to process 512 bits [5], which results in lower performance. As the sequence length increases, the throughput for both approaches tends to increase and hit a peak for long reads. As an example, for a sequence length of 16384, Endeavor achieves between 1.69x and 2.83x speedups for the AMD CPU, while Intel CPU achieves speedups between 1.89x and 3.03x. For both devices, the lower speedups occur for AVX-512 with FP64 precision, while the higher speedups occur for AVX with FP32 precision. Finally, note that the proposed approach is scalable to long reads for any precision and any SIMD width, while GKL fails to run in FP64 for the last tested sequence length (131072) in any SIMD width. 4.3.2 GPU Analysis. Figure 8 displays the TCUPS evolution in the V100, A100, H100, RTX 6000 GPUs in FP32 as the input readhaplotype length increases. Eight curves are displayed, four for Endeavor (purple) and four for the GPU state-of-the-art approach, gpuPairHMM [45]. MI300A is not included, as gpuPairHMM only runs in NVIDIA GPUs. As mentioned in Section 4.2, Endeavor uses 𝑁 = 32 for FP32 in NVIDIA GPUs for sequence lengths up to 1024 and 𝑁 = 64 for longer read-haplotype pairs. Given the warp size of 32, a single warp can process read-haplotype pairs of length up to 2048 with FP32 (blue area). For long read-haplotype pairs, processing is done at the level of multiple warps, leveraging shared memory to store necessary values across warps (orange area) and at the level of multiple thread blocks for very long reads (green area), using thread block clusters, which are only available in H100 and RTX 6000. Therefore, for sequences longer than 16384, V100 and A100 still rely on shared memory and multiple warps. As gpuPairHMM’s benchmark only runs for sequence lengths up to 16384, its curve is not represented for longer read-haplotype pairs. For older GPU architectures (V100 and A100), Endeavor is on par with gpuPairHMM, with more significant improvements in
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
Miguel Graça and Aleksandar Ilic
Table 1: GIAB Dataset Description. Sample Technology #Pairs Min Median Max NA12878 [37] SoLiD 113960218 10 50 481 NA24631 [38] BGISEQ-500 (chr 22) 158287379 10 149 498 NA12878 [35] Illumina 1442661880 10 145 897 NA12878 [36] Ion Torrent 1672115436 10 208 499 NA24149 [33] Chromium (chr 22) 131420448 10 151 651 NA12878 [32] PacBio 866041959 10 215 3396 NA24695 [34] ONT (chr 22) 284136274 10 299 12121
Figure 9: TCUPS evolution with sequence length in Endeavor (FP64).
Figure 10: CARM Roofline for Endeavor and gpuPairHMM.
throughput for sequence lengths between 1024 and 4096, where Endeavor achieves at least 1.4x speedup on V100 and 1.16x speedup on A100. For sequence lengths higher than 16384, the throughput is significantly lower due to increased register pressure and spilling to global memory. For more recent GPU architectures (H100 and RTX 6000), Endeavor achieves a peak throughput of 2.85 TCUPS and 4.70 TCUPS, while gpuPairHMM achieves 2.6 TCUPS and 4.38 TCUPS, which represents a 9.6% and 7.3% improvements, respectively. For these two GPUs, speedups are more significant in the 32-128 range, with 1.6x speedups for sequence lengths of 32 basepairs, and in the 2048-8192 range, where both devices achieve at least 1.2x speedups. Finally, Endeavor achieves good performance for sequence lengths up to 131072 basepairs, which is one order of magnitude higher than the largest sequence length that the GPU state-of-the-art approaches can process. Figure 9 provides Endeavor’s FP64 results on all tested GPUs. Note that gpuPairHMM is not represented in this graph as it only supports single precision. The results show that MI300A achieves the best peak throughput in double precision (1.95 TCUPS), followed by H100 (1.42 TCUPS), A100 (0,6 TCUPS), V100 (0.44 TCUPS) and RTX 6000 PRO (0.01 TCUPS). As mentioned in Section 4.2, the FP64 throughput for RTX 6000 PRO is 1/64 of the peak FP32 throughput, which explains the lower TCUPS when compared to other GPU devices. MI300A achieves the best results on sequence lengths up to 16384. For longer sequences, the throughput is significantly lower due to increased register pressure, similarly to V100 and A100. For
sequence lengths between 32768 and 131072, H100 achieves the best throughput, as it leverages thread block clusters to maintain high performance. In FP32, considering the peak throughputs and TDP of each device, RTX 6000 PRO is the most energy-efficient, providing 0.0078 TCUPS/W, followed by H100 (0.004 TCUPS/W), A100 (0.0038 TCUPS/W), and V100 (0.0029 TCUPS/W). In FP64, MI300A is the most energy-efficient (0.0026 TCUPS/W), followed by H100/A100 (0.002 TCUPS/W), V100 (0.0015 TCUPS/W) and RTX 6000 PRO (0.00016 TCUPS/W). Therefore, Endeavor’s portability allows to choose an optimal device considering different metrics (energy efficiency, targeted sequence length, throughput, or precision). Finally, to further demonstrate Endeavor’s efficiency in using the GPU’s resources, Figure 10 displays the RTX 6000’s roofline (as it is the evaluated device with higher FP32 throughput), generated with CARMTool [30], to evaluate Endeavor (for sequence lengths of 32, 2048, 16384, and 131072) and gpuPairHMM (for sequence lengths of 32, 1024, and 16384), with varying sequence lengths to evaluate the performance of both methods as the input size increases. The roofline shows that gpuPairHMM is closer to the memory-bound region (the main limitations come from the memory subsystems, which hinders possible performance improvements), while Endeavor, as a redefined PairHMM approach, is computebound (the main bottleneck is the amount of compute resources in the hardware).
4.4
Evaluation on Real Datasets
In this section, experiments on GIAB datasets from seven different sequencing technologies are performed to evaluate Endeavor’s flexibility and performance on short-read datasets (sequence lengths in the range of 250-800) and long-read datasets (sequence lengths > 104 ). Table 1 provides a description of the GIAB short-read and long-read datasets that are tested. For short-reads, four different datasets are considered, generated by Illumina, SoLiD, BGI-SEQ500, and Ion Torrent, with varying number of sequence pairs to process. For long-reads, three different datasets are evaluated, generated by Chromium Long Ranger, PacBio, and ONT sequencing technologies, with the last one exhibiting the sequences with longest length (12121 basepairs). The BGI-SEQ500, ONT and Chromium datasets refer only to Chromosome 22, while the remaining datasets refer to the complete human genome. The analysis of these datasets is done in three parts. First, Endeavor’s CPU approach is compared to GATK’s PairHMM AVX-512 implementation, which runs on the Intel Xeon Gold 6438 CPU, using 192 threads, in the seven GIAB datasets. Second, Endeavor’s GPU approach is compared to gpuPairHMM in all GIAB datasets
Endeavor: Efficient PairHMM for Detection of DNA Variants in Genome-Scale Datasets
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
Table 2: Execution Times for 10s Dataset (in milliseconds).
Figure 11: Speedup of Endeavor-CPU over GATK AVX-512 (higher is better)
Figure 12: Speedup of Endeavor-GPU over gpuPairHMM (higher is better).
and datasets from the original paper [45]. Finally, the 10s dataset [9], a typical benchmark used in the literature to test novel approaches for PairHMM, is also evaluated with Endeavor and compared to gpuPairHMM and GKL, as well as other well-known PairHMM implementations in the literature (based on CDP [27] and on interand intra-task parallelization [42]). 4.4.1 Comparison with GATK HaplotypeCaller. Figure 11 provides the speedups achieved over GATK HaplotypeCaller’s AVX-512 implementations on all datasets across the tested CPU devices. For a fair comparison, Endeavor runs AVX-512 with single precision and the same number of threads as GATK HaplotypeCaller (192). On the AMD CPU, Endeavor’s speedups range from 1.94x to 20.48x, while on the Intel CPU, the speedups are between 2.29x and 25.32x, with the lowest results for the PacBio dataset and the highest results for the Illumina dataset. The average speedups for both devices are 10.59x (AMD CPU) and 11.98x (Intel CPU). 4.4.2 Comparison with GPU-based methods. Figure 12 provides the speedups achieved over gpuPairHMM with the Endeavor framework on the GIAB, in addition to datasets considered in the gpuPairHMM [45] paper: a subset of the NA12878 dataset from the 1000 Genomes project [12] and the 10s dataset [9] replicated 4096x. The results show that Endeavor provides slight improvements over gpuPairHMM in A100 and H100 on most of the datasets (with the exception of the SoLiD dataset, where Endeavor achieves the highest
Method Xeon Gold 6438 EPYC Zen 4 GKL [16] 15.39 23.95 Endeavor-CPU 11.81 17.24 Method V100 A100 H100 RTX 6000 MI300A Inter-task [42] (Tile = 1) 25.4 19.4 4.60 16.5 Inter-task [42] (Tile = 2) 14.9 11.8 8.81 11.2 Inter-task [42] (Tile = 4) 8.59 6.75 5.20 6.26 Inter-task [42] (Tile = 6) 7.51 5.58 4.94 4.98 Inter-task [42] (Tile = 8) 7.35 6.36 3.65 4.33 Intra-task [42] (Tile = 1) 4.00 3.23 2.32 2.62 Intra-task [42] (Tile = 2) 4.54 3.81 2.83 2.80 Intra-task [42] (Warp) 1.58 1.25 0.91 0.76 CDP [27] 13.2 14.6 7.59 5.94 gpuPairHMM [45] 8.50 7.41 1.77 1.88 Endeavor-GPU 0.37 0.26 0.20 0.15 0.26
speedup, 1.3x). For V100, Endeavor achieves at least 1.22x speedups, while for the newest GPU architecture (RTX 6000 PRO), Endeavor consistently outperforms gpuPairHMM, with speedups between 1.14x (ONT) and 2.05x (SoLiD), and an average 1.47x speedup. Note that, as Endeavor’s results match the results from gpuPairHMM [45] and GATK HaplotypeCaller [13], substituting the latter’s CPU-based implementation with Endeavor’s GPU approach would result in significant speedups with no loss in precision. When compared to GATK HaplotypeCaller, the achieved speedups range between 15.66x-229.96x (V100), 20.88x-305.04x (A100), 48.74x674.75x (H100), 79.58x-1115.91x (RTX 6000), and 48.82x-758.35x (MI300A), with the lowest speedups on the Ion Torrent dataset and the highest on the SoLiD dataset. 4.4.3 Comparison using the 10s dataset. In this section, the 10s dataset [9] is evaluated on CPUs and GPUs. Although 10s is not representative of a real large-scale dataset, it is a benchmark that is widely used in the literature [4, 6, 21, 42, 44, 45] to test novel approaches for PairHMM. To examine Endeavor’s performance on this dataset in CPUs, GKL is employed for comparison, while for GPUs, gpuPairHMM, as well as a benchmark for NVIDIA GPUs known as Genomics-GPU [27], are evaluated. Genomics-GPU is a software suite that includes a CDP implementation of PairHMM, as well as the implementation presented in [42], which exploits inter-task and intra-task parallelism to increase the algorithm’s throughput. Table 2 provides the execution times for the 10s dataset in all the PairHMM codes, including Endeavor. Because the 10s dataset is small (3550 sequence pairs, with the longest sequence having 263 basepairs), it is possible to process all sequence pairs simultaneously. By doing so, Endeavor achieves 11.81 and 22.24 milliseconds on Intel Xeon Gold 6438 and AMD Zen 4 CPUs, respectively, outperforming GKL by 1.3x and 1.39x. On GPUs, Endeavor achieves 0.37, 0.26, 0.20, and 0.15 milliseconds in V100, A100, H100, and RTX 6000, which leads to a speedup of 4.27x, 4.81x, 4.55x, and 5.06x, respectively, over the second fastest solutions presented in the table (which are in bold). While Genomics-GPU and gpuPairHMM only run on NVIDIA GPUs, Endeavor can also be tested on MI300A, which takes 0,26 milliseconds to process 10s, showcasing Endeavor’s portability across different devices.
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
5
Conclusions
Variant calling is a fundamental problem in bioinformatics with a significant impact on precision medicine, pharmacogenomics, and evolutionary biology. Given its importance, the development of efficient solutions to keep up with its computational demand is paramount. To that end, in this work, Endeavor, a novel parallelization strategy, was presented to leverage row-level parallelization of the PairHMM algorithm and achieve better usage of CPU’s and GPU’s hardware resources. The results on Intel and AMD CPUs showed up to 2.14x improvement in peak PairHMM throughput, as well as up to 25.32x speedup when processing genome-scale datasets. On NVIDIA and AMD GPUs, the proposed solution achieves up to 2.05x speedups over state-of-the-art GPU-based methods on the most recent GPU architectures. Finally, the proposed framework can run on short- and long-read datasets, demonstrating its capability to run on datasets from any sequencing technology and any sequence length on CPU and GPU devices.
Acknowledgments This work was supported by national funds through Fundação para a Ciência e a Tecnologia, I.P. (FCT) under projects UID/50021/2025 (DOI: https://doi.org/10.54499/UID/50021/2025) and UID/PRR/50021/ 2025 (DOI: https://doi.org/10.54499/UID/PRR/50021/2025), LISBOA 2030-FEDER-00869000 (2023.18110.ICDT, VERSACOMP, DOI: https: //doi.org/10.54499/2023.18110.ICDT), and the UI/BD/154603/2022 research grant. We also acknowledge the European Union HE Research and Innovation programme under grant agreement No 10109 2877 (SYCLOPS), and the EuroHPC Joint Undertaking for awarding us access to MareNostrum5 at BSC, Spain.
References [1] Andrew Adinetz. 2014. Adaptive parallel computation with CUDA dynamic parallelism. NVIDIA Corporation) Retrieved January 4 (2014), 2016. [2] Srinivas Aluru, Natsuhiko Futamura, and Kishan Mehrotra. 2003. Parallel biological sequence comparison using prefix computations. J. Parallel and Distrib. Comput. 63, 3 (2003), 264–272. [3] Euan A Ashley. 2016. Towards precision medicine. Nature Reviews Genetics 17, 9 (2016), 507–522. [4] Subho S Banerjee, Mohamed El-Hadedy, Ching Y Tan, Zbigniew T Kalbarczyk, Steve Lumetta, and Ravishankar K Iyer. 2017. On accelerating pair-HMM computations in programmable hardware. In 2017 27th International Conference on Field Programmable Logic and Applications (FPL). IEEE, 1–8. [5] Ravi Bhargava and Kai Troester. 2024. AMD next-generation “Zen 4” core and 4th gen AMD EPYC server CPUs. IEEE Micro 44, 3 (2024), 8–17. [6] Beatrice Branchini, Alberto Zeni, and Marco D Santambrogio. 2021. A Methodology for Accelerating Variant Calling on GPU. (2021). [7] Benjamin Buchfink, Klaus Reuter, and Hajk-Georg Drost. 2021. Sensitive protein alignments at tree-of-life scale using DIAMOND. Nature methods 18, 4 (2021), 366–368. [8] Christiam Camacho, Grzegorz M Boratyn, Victor Joukov, Roberto Vera Alvarez, and Thomas L Madden. 2023. ElasticBLAST: accelerating sequence search via cloud computing. BMC bioinformatics 24, 1 (2023), 117. [9] M Carneiro. 2013. Optimization of a Haplotype Pair-HMM class for GPU/FPGA and AVX processing. https://github.com/MauricioCarneiro/PairHMM. [10] Tiago Carneiro Pessoa, Jan Gmys, Francisco Heron de Carvalho Júnior, Nouredine Melab, and Daniel Tuyttens. 2018. GPU-accelerated backtracking using CUDA Dynamic Parallelism. Concurrency and Computation: Practice and Experience 30, 9 (2018), e4374. [11] Ming-Hung Chen, Mao-Jan Lin, Yu-Cheng Li, and Yi-Chang Lu. 2019. Banded PairHMM Algorithm for DNA Variant Calling and Its Hardware Accelerator Design. In 2019 IEEE 19th International Conference on Bioinformatics and Bioengineering (BIBE). IEEE, 563–566. [12] 1000 Genomes Project Consortium et al. 2015. A global reference for human genetic variation. Nature 526, 7571 (2015), 68.
Miguel Graça and Aleksandar Ilic
[13] Mark A DePristo, Eric Banks, Ryan Poplin, Kiran V Garimella, Jared R Maguire, Christopher Hartl, Anthony A Philippakis, Guillermo Del Angel, Manuel A Rivas, Matt Hanna, et al. 2011. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature genetics 43, 5 (2011). [14] Richard Durbin, Sean R Eddy, Anders Krogh, and Graeme Mitchison. 1998. Biological sequence analysis: probabilistic models of proteins and nucleic acids. Cambridge university press. [15] Sean R Eddy. 2004. What is dynamic programming? Nature biotechnology 22, 7 (2004), 909–910. [16] Patrick Foley, Abirami Prabhakaran, Karthik Gururaj, Mishali Naik, Shiva Gopalan, Aleksandr Shargorodskiy, and Ernesto Brau. 2017. Accelerate Genomics Research with the Broad-Intel Genomics Stack. [17] Efstathia Giannopoulou, Theodora Katsila, Christina Mitropoulou, EvangeliaEirini Tsermpini, and George P Patrinos. 2019. Integrating next-generation sequencing in the clinical pharmacogenomics workflow. Frontiers in pharmacology 10 (2019), 384. [18] Richard A Gibbs. 2020. The human genome project changed everything. Nature Reviews Genetics 21, 10 (2020), 575–576. [19] Mark Harris, Shubhabrata Sengupta, and John D Owens. 2007. Parallel prefix sum (scan) with CUDA. GPU gems 3, 39 (2007), 851–876. [20] Taishan Hu, Nilesh Chitnis, Dimitri Monos, and Anh Dinh. 2021. Next-generation sequencing technologies: An overview. Human Immunology 82, 11 (2021). [21] Sitao Huang, Gowthami Jayashri Manikandan, Anand Ramachandran, Kyle Rupnow, Wen-mei W Hwu, and Deming Chen. 2017. Hardware acceleration of the pair-HMM algorithm for DNA variant calling. In Proceedings of the 2017 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays. [22] Miten Jain, Sergey Koren, Karen H Miga, Josh Quick, Arthur C Rand, Thomas A Sasani, John R Tyson, Andrew D Beggs, Alexander T Dilthey, Ian T Fiddes, et al. 2018. Nanopore sequencing and assembly of a human genome with ultra-long reads. Nature biotechnology 36, 4 (2018), 338–345. [23] Hákon Jónsson, Patrick Sulem, Birte Kehr, Snaedis Kristmundsdottir, Florian Zink, Eirikur Hjartarson, Marteinn T Hardarson, Kristjan E Hjorleifsson, Hannes P Eggertsson, Sigurjon Axel Gudjonsson, et al. 2017. Whole genome characterization of sequence diversity of 15,220 Icelanders. Scientific data 4, 1 (2017), 1–9. [24] Ali Khajeh-Saeed, Stephen Poole, and J Blair Perot. 2010. Acceleration of the Smith–Waterman algorithm using single and multiple graphics processors. J. Comput. Phys. 229, 11 (2010), 4247–4258. [25] Daniel C Koboldt. 2020. Best practices for variant calling in clinical sequencing. Genome Medicine 12, 1 (2020), 91. [26] Enliang Li, Subho S Banerjee, Sitao Huang, Ravishankar K Iyer, and Deming Chen. 2021. Improved gpu implementations of the pair-hmm forward algorithm for dna sequence alignment. In 2021 IEEE 39th International Conference on Computer Design (ICCD). IEEE, 299–306. [27] Zhuren Liu, Shouzhe Zhang, Justin Garrigus, and Hui Zhao. 2023. GenomicsGPU: A Benchmark Suite for GPU-accelerated Genome Analysis. In 2023 IEEE International Symposium on Performance Analysis of Systems and Software (ISPASS). IEEE, 178–188. [28] Chengwei Luo, Despina Tsementzi, Nikos Kyrpides, Timothy Read, and Konstantinos T Konstantinidis. 2012. Direct comparisons of Illumina vs. Roche 454 sequencing technologies on the same microbial community DNA sample. PloS one 7, 2 (2012), e30087. [29] Bui Quang Minh, Heiko A Schmidt, Olga Chernomor, Dominik Schrempf, Michael D Woodhams, Arndt Von Haeseler, and Robert Lanfear. 2020. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Molecular biology and evolution 37, 5 (2020), 1530–1534. [30] José Morgado, Leonel Sousa, and Aleksandar Ilic. 2024. CARM Tool: Cache-Aware Roofline Model Automatic Benchmarking and Application Analysis. In 2024 IEEE International Symposium on Workload Characterization (IISWC). IEEE, 68–81. [31] Sergey Nurk, Sergey Koren, Arang Rhie, Mikko Rautiainen, Andrey V Bzikadze, Alla Mikheenko, Mitchell R Vollger, Nicolas Altemose, Lev Uralsky, Ariel Gershman, et al. 2022. The complete sequence of a human genome. Science 376, 6588 (2022), 44–53. [32] National Institute of Health. 2024. NA12878 Pacific Biosciences BAM Dataset. Available at https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/NA12878/PacBio_ SequelII_CCS_11kb/HG001.SequelII.pbmm2.hs37d5.whatshap.haplotag.RTG. trio.bam. NA24149 Chromium Long Ranger [33] National Institute of Health. 2024. BAM Dataset. Available at https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/ AshkenazimTrio/analysis/10XGenomics_ChromiumGenome_LongRanger2.0_ 06202016/HG003_NA24149_father/NA24149_GRCh37.bam. [34] National Institute of Health. 2024. NA24695 Oxford Nanopore Technologies BAM Dataset. Available at https://ftp-trace.ncbi.nlm. nih.gov/giab/ftp/data/ChineseTrio/HG007_NA24695-hu38168_mother/ UCSC_Ultralong_OxfordNanopore_Promethion/HG007_GRCh37_ONTUL_UCSC_20200109.phased.bam. [35] National Institute of Health. 2024. NIH NA12878 Illumina BAM Dataset. Available at https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/NA12878/NIST_ NA12878_HG001_HiSeq_300x/RMNISTHS_30xdownsample.bam.
Endeavor: Efficient PairHMM for Detection of DNA Variants in Genome-Scale Datasets
[36] National Institute of Health. 2024. NIH NA12878 Ion Torrent BAM Dataset. Available at https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/NA12878/ion_exome/ IonXpress_020_rawlib.b37.bam. [37] National Institute of Health. 2024. NIH NA12878 SoLiD BAM Dataset. Available at https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/technical/NA12878_data_other_ projects/alignment/NA12878.SOLID.SRP012400.Xprize_SRR643700.bam. [38] National Institute of Health. 2024. NIH NA24631 BGISEQ500 BAM Dataset. https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/ChineseTrio/HG005_NA24631_ son/NIST_BGIseq_2x150bp_100x/GRCh38/HG005_GRCh38_BGIseq-2x150100x_NIST_20211126.bam [39] Nathan D Olson, Justin Wagner, Nathan Dwarshuis, Karen H Miga, Fritz J Sedlazeck, Marc Salit, and Justin M Zook. 2023. Variant calling and benchmarking in an era of complete human genome sequences. Nature Reviews Genetics 24, 7 (2023), 464–483. [40] Johan Peltenburg, Shanshan Ren, and Zaid Al-Ars. 2016. Maximizing systolic array efficiency to accelerate the PairHMM forward algorithm. In 2016 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). IEEE. [41] Shanshan Ren, Koen Bertels, and Zaid Al-Ars. 2017. GPU-accelerated GATK haplotypecaller with load-balanced multi-process optimization. In 2017 IEEE 17th International Conference on Bioinformatics and Bioengineering (BIBE). IEEE, 497–502. [42] Shanshan Ren, Koen Bertels, and Zaid Al-Ars. 2018. Efficient acceleration of the pair-hmms forward algorithm for gatk haplotypecaller on graphics processing units. Evolutionary Bioinformatics 14 (2018), 1176934318760543. [43] Tony Robinson, Jim Harkin, and Priyank Shukla. 2021. Hardware acceleration of genomics data analysis: challenges and opportunities. Bioinformatics 37, 13 (2021), 1785–1795.
HPDC ’26, July 13–16, 2026, Cleveland, OH, USA
[44] Davide Sampietro, Chiara Crippa, Lorenzo Di Tucci, Emanuele Del Sozzo, and Marco D Santambrogio. 2018. Fpga-based pairhmm forward algorithm for dna variant calling. In 2018 IEEE 29th International Conference on Application-specific Systems, Architectures and Processors (ASAP). IEEE, 1–8. [45] Bertil Schmidt, Felix Kallenborn, Alexander Wichmann, Alejandro Chacon, and Christian Hundt. 2026. gpuPairHMM: High-Speed Pair-HMM Forward Algorithm for DNA Variant Calling on GPUs. IEEE Transactions on Computational Biology and Bioinformatics (2026), 1–8. doi:10.1109/TCBBIO.2026.3657252 [46] Roman Snytsar. 2023. PairHMM Improvements for Modern Instruction Set Architectures. In 2023 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). IEEE, 3328–3331. [47] TOP500.org. [n. d.]. TOP500 June 2025. https://www.top500.org/lists/top500/ 2025/06/. [Online; Jun-2025]. [48] Jin Wang and Sudhakar Yalamanchili. 2014. Characterization and analysis of dynamic parallelism in unstructured GPU applications. In 2014 IEEE International Symposium on Workload Characterization (IISWC). IEEE, 51–60. [49] Rick Wertenbroek and Yann Thoma. 2019. Acceleration of the Pair-HMM forward algorithm on FPGA with cloud integration for GATK. In 2019 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). IEEE, 534–541. [50] Chunlin Xiao, Justin Zook, Shane Trask, Stephen Sherry, and Genome in-a Bottle Consortium. 2014. GIAB: Genome reference material development resources for clinical sequencing. Cancer Research 74, 19_Supplement (2014), 5328–5328. [51] Byung-Jun Yoon. 2009. Hidden Markov models and their applications in biological sequence analysis. Current genomics 10, 6 (2009), 402–415. [52] Zhonghai Zhang, Yewen Li, Ke Meng, Chunming Zhang, and Guangming Tan. 2026. Faster and Cheaper: Pushing the Sequence Alignment Throughput with Commercial CPUs. In Proceedings of the 31st ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming. 466–479.