GAPiM: Discovering Genetic Variations on a Real Processing-in-Memory System

GAPiM: Discovering Genetic Variations on a Real Processing-in-Memory System | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article GAPiM: Discovering Genetic Variations on a Real Processing-in-Memory System Naomie Abecassis, Juan Gómez-Luna, Onur Mutlu, Ran Ginosar, Aphélie Moisson-Franckhauser, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-3245756/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Variant calling is a fundamental stage in genome analysis that identifies mutations (variations) in a sequenced genome relative to a known reference genome. Pair-HMM is a key part of the variant calling algorithm and its most compute-intensive part. In recent years, Processing-in-Memory (PiM) solutions, which consist of placing compute capabilities near/inside memory, have been proposed to speed up the genome analysis pipeline. We implement the Pair-HMM algorithm on a commercial PiM platform developed by UPMEM. We modify the Pair-HMM algorithm to make it more suitable for PiM execution with acceptable loss of accuracy. We evaluate our implementation on single chromosomes and whole genome sequencing datasets, demonstrating up to 2x speedup compared to existing CPU accelerations and up to 3x speedup compared to FPGA accelerations. Processing in Memory Pair-HMM Variant calling DNA sequence analysis High Performance Computing Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Introduction The completion of the Human Genome Project in 2003 and the advent of next-generation sequencing (NGS) technology in the early 2000s marked the beginning of a new era in DNA sequencing [1,2], enabling faster and more efficient sequencing of large amounts of DNA, and driving the growth in the field of genomics. Over the years, the progress made in the field has led to a rapid increase in the amount of available genomic data and a decrease in the cost of sequencing [3]. Genomic data is analyzed to understand the genetic information carried by a DNA molecule and to identify DNA alterations that can cause diseases, such as genetic disease and cancer. The process of identifying variations in the sequenced data compared to a reference genome is known as variant calling and can take up to 40% of the total time spent in the genome sequencing pipeline [4]. Variant calling is an essential step in DNA sequencing pipelines because it enables the identification of differences, or variations, between the DNA sequence being analyzed and a reference genome or other sequences. These variations can include single nucleotide polymorphisms (SNPs, i.e., a DNA sequence variation that occurs when a single nucleotide (adenine, thymine, cytosine, or guanine) in the genome sequence is altered), small insertions or deletions (indels 2 ), or larger structural variations such as copy number variations (CNVs, i.e., refer to the genetic trait involving the number of copies of a particular gene present in the genome of an individual) or genomic rearrangements. Variant calling is important for a variety of applications, including disease diagnosis [18,20], personalized medicine, and understanding evolutionary relationships between organisms [19]. By identifying genetic variations, variant calling can help researchers and clinicians identify disease-causing mutations, develop targeted therapies, and better understand the genetic basis of disease. In addition, variant calling can be used to study the genetic diversity of populations or track the spread of infectious diseases. With the increased need for genomic data and continuous advancements in technology, it is predicted that the amount of genomic data generated will continue to rise [5]. Most of the computational algorithms developed for genome analysis are designed for use in a CPU system. The analysis of the sequenced data requires intensive data movement between DRAM and processing elements which are expensive, both in terms of time and energy [6,7,8,9]. In addition, the variant calling task is a highly parallelizable task as different regions in the genome can be analyzed independently of the others. While CPUs, FPGAs and GPUs are commonly used in variant calling to take advantage of parallel processing [30,31,32,33,35,36], the speedup achieved by these hardware architectures is limited as explained in the evaluation section, emphasizing the need for a new architecture paradigm. In this paper we focus on Processing-in-Memory computer architecture paradigm whose main purpose is alleviating the memory wall , which limits the performance and energy inefficiency of a conventional von Neumann architecture for data intensive applications such as genome analysis. Processing in Memory (PiM) architecture is a promising solution for genomic data analysis considering the expansion of the available data [10,11,12,13,14,15]. PiM architectures can solve data access bottleneck and lower the energy consumption. Pair-HMM is the most time-consuming part of the variant caller in the GATK pipeline. Thus, we investigate the PiM architecture as a potential solution of acceleration. PiM architecture provides high bandwidth and low latency memory access. Additionally, PiM provides a large amount of compute parallelism and enables a throughput that scales with memory size. In this work, we show that a PiM architecture can efficiently implement an end-to-end solution for genomic analysis, even when certain components of the genome analysis pipeline are not highly data intensive. The main contributions of this work are as follows: We rethink and adapt the Pair-HMM algorithm to make it suitable for a real Processing-in-Memory architecture. We implement the Pair-HMM algorithm on a commercial PiM system with more than 2500 PiM cores (DPUs). We comparatively evaluate the PiM Pair-HMM implementation and demonstrate up to 2x time performance improvement compared to high-performance CPU and up to 3x performance improvement compared to SOTA FPGA solutions. Background A. The Genome Analysis ToolKit (GATK) GATK (Genome Analysis Toolkit) is a collection of tools developed by the Broad Institute to identify SNPs and indels in genomic data [16]. SNPs are the most common sort of genetic variation among human beings [17]. SNPs and indels play an important role in the study of human cancer genes [18,19,20]. The GATK best practices pipeline is one of the most widely used for variant calling, due to its high accuracy [21]. This pipeline is composed of three essential parts. The first one is data pre-processing, the second is variant discovery, and the third is variant evaluation. GATK is an open-source software available on GitHub [22]. B. HaplotypeCaller algorithm The second part of the GATK pipeline, variant discovery, is implemented by the HaplotypeCaller algorithm. This algorithm gets as inputs a processed BAM file (a set of aligned DNA reads) and a reference sequence file, and outputs a VCF file (variant calling file). The HaplotypeCaller \(1-\varvec{\epsilon }\) algorithm can take up to 35% of the complete genome analysis pipeline time as reported in [23]. The HaplotypeCaller algorithm has four main steps: Active region determination based on the presence of a sufficient number of variations that exceed a predetermined threshold [16]. Active regions, i.e., with sufficient variations, are the target of the next three steps. • De-Bruijn like graph assembly is used to assemble the reads of the region and determine a list of haplotypes [24]. A haplotype is a sequence that represents a certain combination of alleles (alternative gene forms) present in a region. Then the Smith-Waterman algorithm is used to align each haplotype to the reference haplotype to identify potential genetic variations. • Pair-HMM forward algorithm is used to perform a pairwise alignment of each read against each haplotype [25]. The algorithm calculates the likelihood of each haplotype given each read by summing the likelihoods of all possible alignments of the read to the haplotype. This likelihood reflects the probability that the read and the haplotype are related, based on their sequence similarity. • Genotype sampling estimates the genotypes (combination of alleles) at each variant site in the region using Bayes’ theorem in combination with the previously obtained likelihoods. C. Algorithm analysis The GATK HaplotypeCaller software includes a feature that tracks the time required to complete the Smith-Waterman and the Pair-HMM steps of the workflow. Additionally, we used profiling tools [55,56] and present the time repartition of the program in Fig. 1 (running on an Intel Core I7-5820K x86_64, 6 cores (12 threads), 3.30GHz, 64GB). The figure represents the average relative execution time for chromosomes 1–3 of sample NA12878 (commonly used sample from a female individual of European descent). As observed, the Pair-HMM component of the algorithm accounts for a significant portion of 49% of the total time, prompting our decision to accelerate this section. Moreover, we check the memory boundedness of the HaplotypeCaller algorithm by using hardware performance counters [52]. As reported in [52], the memory boundedness of the program, which is a component of the non-execution of the CPU ( cycle_activity.cycles_no_execute ), can be computed by the following formula: $$memory bound= \text{m}\text{a}\text{x}(RESOURCE\_STALLS:SB , CYCLE\_ACTIVITY:STALLS\_L1D\_PENDING) \left(1\right)$$ We present in Table 1 the results of this analysis for a single chromosome (chromosome 1 randomly chosen) and for a WGS running. As can be observed, the memory boundedness of the HaplotypeCaller algorithm is low, which means that memory access is unlikely to become a bottleneck of a HaplotypeCaller. The intuitive conclusion is that PiM would not be the most efficient solution for HaplotypeCaller. Nevertheless, there are three main advantages to PiM architecture even this non-memory-bound program can benefit from: PiM architectures have a large number of execution units that provide large amounts of parallelism. PiM compute throughput scales with memory capacity. PiM systems are not limited in terms of memory capacity (or not as much) as other accelerators (e.g., GPUs or FPGAs with HBM memory). Table 1 Memory boundedness of HaplotypeCaller algorithm Chromosome 1 WGS resource_stalls.sb 146,149,728,421 7,680,287,475,816 cycle_activity.stalls_l1d_pending 791,177,991,761 65,333,865,516,882 cycle_activity.cycles_no_execute 7,590,369,180,122 510,501,141,773,479 Memory boundedness 10% 12.8% D. Pair-HMM forward algorithm In this work we focus on the Pair-HMM step which is the most resource-intensive part of the HaplotypeCaller algorithm. In this part we compute the likelihood of each read against each haplotype. The likelihood calculation is based on a pairwise alignment using a Hidden Markov Model (represented by a state machine) and the qualities of the read DNA bases (qualities produced by a sequencer). Given a DNA read R of length n and a haplotype H of length m we define the transitions of the state machine as the probabilities of the possible changes in state that occur during the process of aligning the read to the haplotype. In other words, the state machine represents the different states that the alignment can go through, such as a match (state M), an insertion (state I) or a deletion (state D) in the read. The transitions between the states are defined by the following probabilities: $${\delta } ={10}^{-\frac{{\text{Q}}_{\text{d}\text{e}\text{l}}}{10}}, {\iota } ={10}^{-\frac{{\text{Q}}_{\text{i}\text{n}\text{s}}}{10}}, {\epsilon } ={10}^{-\frac{{\text{Q}}_{\text{g}\text{c}\text{p}}}{10}}$$ $$prio{r}_{i,j}= \left\{\begin{array}{c}1-{10}^{-\frac{{Q}_{b}}{10}} if rea{d}_{i}=haplotyp{e}_{j}\\ \frac{{10}^{-\frac{{Q}_{b}}{10}}}{3} if rea{d}_{i}\ne haplotyp{e}_{j}\end{array}\right.$$ 2 Currently, \({\text{Q}}_{\text{d}\text{e}\text{l}}\) , \({\text{Q}}_{\text{i}\text{n}\text{s}}\) and \({\text{Q}}_{\text{g}\text{c}\text{p}}\) (respectively the deletion, insertion, and gap continuation qualities) are included in the PHMM model but are not produced by sequencers and thus default values determined by GATK [22] are used for each base. However, \({\text{Q}}_{\text{b}}\) is referred to as the quality of each base in the read and is provided by the sequencer machine. The finite state machine of Fig. 2 models the transitions between the three different states ( M , the match state, I , the insertion state and D , the deletion state). Using this state machine with dynamic programming, the algorithm fills up three matrices, M , I and D as shown in Equations (3;1) and (3;2) . Each element \(i,j\) of the matrix represents the total likelihood of all paths from the beginning of the read to position \(i\) such that the path ends with the corresponding state (match state for matrix M , deletion state for matrix D , and insertion state for matrix I ). Initialization : \(\forall 0\le i<n, 0\le j<m,\) $$\left\{\begin{array}{c}{M}_{i,0}=0 , {I}_{i,0}=0, {D}_{i,0}=\frac{1}{m}\\ {M}_{0,j}=0 , {I}_{0,j}=0, {D}_{0,j}=0\\ {M}_{\text{0,0}}=0 , {I}_{\text{0,0}}=0, {D}_{\text{0,0}}=0\end{array}\right. (3;1)$$ Recursion step : \(\forall 1\le i\le n, 1\le j\le m,\) \(\left\{\begin{array}{c}{M}_{i,j}=prio{r}_{i,j}*[\left(1-\left({\iota }+{\delta }\right)\right){ M}_{i-1,j-1}+\\ +\left(1-{\epsilon }\right) {(\text{I}}_{\text{i}-1,\text{j}-1} + {\text{D}}_{\text{i}-1,\text{j}-1}\left)\right]\\ {I}_{i,j}= {{\epsilon }\text{*}\text{I}}_{\text{i}-1,\text{j}} + {{\iota }\text{*}\text{M}}_{\text{i}-1,\text{j}}\\ {D}_{i,j}={{\epsilon }\text{*}\text{D}}_{\text{i},\text{j}-1} + {{\delta }\text{*}\text{M}}_{\text{i},\text{j}-1}\end{array}\right.(3;2\) ) Termination step : $$L=\sum _{j=1}^{m}{M}_{n,j}+{I}_{n,j} (3;3)$$ The terminations step computes the final likelihood for the current read-haplotype. Related work Genome analysis requires high-performance computing to process and store large and continuously growing genomic datasets. Several processing in memory (PiM) solutions have been proposed to accelerate various parts of genome analysis pipeline. For example, the development of UpVC [13] on the UPMEM architecture has sped up the read alignment and variant calling process and yielded a time and energy improvement compared to GPU or FPGA implementations, using a simpler and less accurate algorithm than PHMM. Other works have focused on accelerating the read alignment part of the pipeline by using in-memory databases [26,27], or by using resistive content addressable memory (BioSEAL [11]). RASSA [28] accelerates read pre-alignment using approximate search-capable resistive content addressable memory. GRIM filter [29] integrates computation within a logic layer optimized to exploit 3D stacked memory layers to filter seed location. More recently GenStore [51] has implemented filters for in-memory reads for genome analysis tasks, and research conducted by [15] demonstrates that the wavefront algorithm achieves greater read alignment throughput on PiM architecture. In contrast to the majority of prior research, it is crucial to emphasize that our work is conducted on a real Processing-in-Memory architecture. GAPiM is the first GATK pipeline implementation on a commercial PiM platform. To our knowledge, GATK has never been implemented on an actual PIM platform before. Typical GATK accelerators are based on CPU, FPGA, and GPU. Most of these accelerations implement the Pair-HMM step. CPU acceleration solutions use a combination of SIMD extensions such as AVX and AVX-512 and multithreading [30]. FPGA-based solutions mainly focus on accelerating the floating-point calculations of Pair-HMM. Most FPGA solutions use systolic arrays [31,32,33]. Another solution for FPGA is to create a specialized unit to execute the PairHMM command and replicate it to provide efficient parallelization [34]. But these FPGA solutions suffer from limited on-chip memory. Our solution, GAPiM, alleviates the memory wall which may limit the CPU and FPGA implementations and hence enables better performance and scalability. We believe PiM approaches are also beneficial to investigate for emerging deep learning based variant calling tools [58]. We leave the investigation of the evaluation, design, and analysis of such tools to future work. GPU accelerations are also used since the Pair-HMM algorithm is easily parallelizable, focusing on intra-task and inter-task parallelism [35] and allowing to retrieve the optimal alignment [36]. Recently, an ASIC that replaces floating point multiplication by 20-bit addition in log domain employs bound checks to maintain correctness has been proposed in [37]. GAPiM provides a better scalability and flexibility when compared with GPU or ASIC implementations. UPMEM PiM architecture In this work we utilize the PiM platform developed by UPMEM [39]. It is the first commercially available PiM solutions. It has been shown to be effective in speeding up memory-intensive workloads [13,40,41]. UPMEM platform is a standard DDR4-2400 DIMM with several PiM chips, each one composed of 8 parallel processors, called DPUs (DRAM Processing Units). Although this paper we refer to DPUs as PiM cores. The PiM cores are operated by a host CPU that sends the program to be executed and collects the processing results when PiM cores complete their processing tasks [39,42,43]. A. Chip level organization Within each memory chip (illustrated in Figure 3 ), there are 8 asynchronous PiM cores that operate independently of one another. Each PiM core is associated with a 64MB DRAM bank called MRAM (Main RAM) that can also be accessed by the host CPU. Different PiM cores cannot communicate with each other directly. The host CPU can transfer data to/from MRAM banks of individual PiM cores. This design allows each computing unit to efficiently process its own fragment of the dataset. Calculations are performed in-situ, within each unit at a memory bandwidth of almost 1GB/s. B. The PiM core A PiM core is a 32-bit RISC-based processor (with 64-bit capabilities) running at 400MHz. It is a multithreaded processor with up to 24 hardware threads, called PiM threads. The PiM core has a pipeline depth of 14 stages, however, only the last three stages of the pipeline can execute in parallel with stages of the next instruction in the same thread. Therefore, instructions from the same thread must be dispatched 11 cycles apart, requiring at least 11 threads to fully utilize the pipeline [39,42,43]. Each PiM core has an exclusive access to a 24KB IRAM (Instruction RAM), and to a 64KB WRAM (Working RAM) shared by all the threads, as well as the 64MB MRAM as illustrated in Fig. 4 . The WRAM has a lower access latency than MRAM and can be accessible in one cycle. To allow several threads to access a same resource we use a mutual exclusion (mutex) mechanism, i.e., a synchronization object is used to control access to a shared resource and ensure that only one PiM thread can access that resource at a time, avoiding conflicts and data inconsistencies. C. The PiM core A PiM core is a 32-bit RISC-based processor (with 64-bit capabilities) running at 400MHz. It is a multithreaded processor with up to 24 hardware threads, called PiM threads. The PiM core has a pipeline depth of 14 stages, however, only the last three stages of the pipeline can execute in parallel with stages of the next instruction in the same thread. Therefore, instructions from the same thread must be dispatched 11 cycles apart, requiring at least 11 threads to fully utilize the pipeline [39,42,43]. Each PiM core has an exclusive access to a 24KB IRAM (Instruction RAM), and to a 64KB WRAM (Working RAM) shared by all the threads, as well as the 64MB MRAM as illustrated in Figure 4 . The WRAM has a lower access latency than MRAM and can be accessible in one cycle. To allow several threads to access a same resource we use a mutual exclusion (mutex) mechanism, i.e., a synchronization object is used to control access to a shared resource and ensure that only one PiM thread can access that resource at a time, avoiding conflicts and data inconsistencies. In this setup, programmers explicitly transfer data between MRAM and WRAM using DMA calls which must be 8-byte aligned and execute load and store operations on WRAM. Adapting HaplotypeCaller algo-rithm to UPMEM PiM architecture We modified the Pair-HMM algorithm to enable its efficient implementation on UPMEM PiM architecture while optimizing the hardware utilization. A. Proposed log domain fixed point model The need to redesign the Pair-HMM for PiM core stems from two observations, as follows: first, multiplications in PiM core are much slower than additions (approximately 10 MOPS for multiplication vs. 60 MOPS for 32-bit integer addition) [40]. Second, floating point operations are significantly less efficient on PiM cores compared to integer operations (approximately 5 MFLOPS for additions) [40]. This is because the PiM core does not implement floating point ALUs and these operations are emulated by the UPMEM runtime library in software [44, 45]. The Pair-HMM algorithm requires a high number of multiplications (Equations (3;1) and (3;2) ) and uses floating point calculations (see priors and transition probabilities). These calculations are significantly slower than integer addition/subtraction on the wimpy PiM cores: multiplication takes six times longer and floating-point addition takes twelve times longer than integer addition. To overcome this limitation, we first convert all equations to log domain, replacing multiplications with additions. For example, the first equation of Eq. (3; 2 ) is transformed as follows, computing the value of \(\text{l}\text{o}\text{g}\left({M}_{i,j}\right)\) rather than \({M}_{i,j}\) : \(\text{log}\left({M}_{i,j}\right)=\text{l}\text{o}\text{g}\left(prio{r}_{i,j}*\right[\left(1-\left({\iota }+{\delta }\right)\right){ M}_{i-1,j-1}+\left(1-{\epsilon }\right) {(\text{I}}_{\text{i}-1,\text{j}-1} + {\text{D}}_{\text{i}-1,\text{j}-1}\left)\right] )\) \(=\text{log}\left(prio{r}_{i,j}\right)+\text{l}\text{o}\text{g}(\left(1-\left({\iota }+{\delta }\right)\right){ M}_{i-1,j-1}+\left(1-{\epsilon }\right) {(\text{I}}_{\text{i}-1,\text{j}-1} + {\text{D}}_{\text{i}-1,\text{j}-1}))\) \(=\text{log}\left(prio{r}_{i,j}\right)+\text{log}\left({10}^{\text{log}\left(\left(1-\left({\iota }+{\delta }\right)\right){ M}_{i-1,j-1}\right)}+{10}^{\text{l}\text{o}\text{g}(\left(1-{\epsilon }\right) {(\text{I}}_{\text{i}-1,\text{j}-1} + {\text{D}}_{\text{i}-1,\text{j}-1}))}\right)\) \(=\text{log}\left(prio{r}_{i,j}\right)+\text{log}\left({10}^{\text{log}\left(1-\left({\iota }+{\delta }\right)\right)+\text{log}\left({ M}_{i-1,j-1}\right)}+{10}^{\text{log}\left(1-{\epsilon }\right)+\text{l}\text{o}\text{g}{(\text{I}}_{\text{i}-1,\text{j}-1} + {\text{D}}_{\text{i}-1,\text{j}-1})}\right)\) \(=\text{log}\left(prio{r}_{i,j}\right)+\text{log}\left({10}^{a}+{10}^{b}\right)=\) \(=\text{log}\left(prio{r}_{i,j}\right)+\text{log}\left({10}^{a}\left(1+{10}^{b-a}\right)\right)=\) \(=\text{log}\left(prio{r}_{i,j}\right)+\text{log}\left({10}^{a}\right)+\text{log}\left(1+{10}^{b-a}\right)=\) \(=\text{log}\left(prio{r}_{i,j}\right)+a+\text{l}\text{o}\text{g}(1+{10}^{b-a})\) $$\text{W}\text{h}\text{e}\text{r}\text{e} a=\text{log}\left(1-\left({\iota }+{\delta }\right)\right)+\text{log}\left({ M}_{i-1,j-1}\right) \text{a}\text{n}\text{d} b=\text{log}\left(1-{\epsilon }\right)+\text{l}\text{o}\text{g}{(\text{I}}_{\text{i}-1,\text{j}-1} + {\text{D}}_{\text{i}-1,\text{j}-1}).$$ Several methods for PiM execution of log and other transcendental functions have been proposed ‎[57]. To implement log-domain additions, we use a look-up table (LUT) to quickly calculate the value of log(1 + 10 X ). This allows us to compute the value of the logarithmic \(sum \left(log\right(x+y\left)\right)\) according to the above equations. Such LUTs are implemented in PiM core memory prior to the operation. The log values of constants and known-in-advance initial values are precalculated by the host CPU and stored in the PiM core memory prior to the operation. Our second contribution is replacing the floating-point calculations by fixed point ones. This reduces the number of slow multiplications and eliminates the need for emulated floating-point operations, considerably improving the performance and efficiency of the algorithm implementation on PiM cores. We modify the Pair-HMM recursion step as shown in Eq. (4) below, where \({+}_{F}\) represents a fixed-point regular addition and \({+}_{L}\) represents a logarithmic addition (i.e., log(a + b)) performed using the LUT approach. $$\left\{\begin{array}{c}{M}_{i,j}={log}\left(prio{r}_{i,j}\right){+}_{\varvec{F}}[\left(log\left(1-\left(\iota +\delta \right)\right){+}_{\varvec{F}}{ M}_{i-1,j-1}\right)\\ {+}_{L}\left({log}\left(1-\epsilon \right){+}_{\varvec{F}} {(I}_{i-1,j-1} {+}_{L} {D}_{i-1,j-1}\right)\left)\right]\\ {I}_{i,j}= {{log}\left(\epsilon \right){+}_{\varvec{F} }I}_{i-1,j} {+}_{L} log({\iota ){+}_{\varvec{F}} M}_{i-1,j}\\ {D}_{i,j}={{log}\left(\epsilon \right){+}_{\varvec{F}} D}_{i,j-1} {+}_{L} {{log}\left(\delta \right){+}_{\varvec{F}} M}_{i,j-1}\end{array}\left(4\right)\right.$$ B. Accuracy analysis Replacing the floating-point operations by fixed-point ones affects the computation precision and therefore, impacts the accuracy of the algorithm (i.e., the ability of the algorithm to correctly identify all the variants). We investigate this accuracy impact by comparing our modified GATK implementation with the original GATK HaplotypeCaller algorithm. We use two samples from the International Genome Sample Resource (IGSR) [46], NA19685 and NA12878 , for each one WES (Whole Exome Sequencing, i.e., protein-coding regions of genes in the genome [47]) with 30x coverage) and WGS (Whole Genome Sequencing) with low coverage. These genomes are aligned to the Hg38 reference and the GATK best practices pipeline is applied on each one. We optimize the fixed-point format of our calculations using the Matlab Fixed-Point Designer tool [48]. This tool suggests the optimal fixed-point configuration for the variables involved, ensuring accurate and reliable results. It enables us to select the appropriate wordlength for each variable, minimizing the accuracy loss in the Pair-HMM execution. To compare the performances of the variant callers we first define the following: True positive (TP): variants found by both the original GATK and by fixed-point precision GATK. • False Negative (FN): variants found by original GATK but not found by fixed-point precision GATK. • False Positive (FP): non-existing variants found by fixed-point precision GATK. The metrics of precision (P) and sensitivity (S) are defined as follow, the reference being the original GATK program: $$P=\frac{TP}{TP+FP}, S=\frac{TP}{TP+FN} \left(5\right)$$ We use the vcfeval tool [50] to compute the different metrics testing four samples (Whole Exome Sequencing (WES) and Whole Genome Sequencing (WGS)) of two different subjects (NA19685 and NA12878). Our analysis shows that fixed-point utilization affects the VCF output, as can be seen in Table 2 below. However, the extent to which this impact is acceptable or not may depend on application requirements and should be evaluated accordingly. Table 2 Fixed-point impact on sensitivity and precision Sample NA19685 WES NA19685 WGS NA19685 CHR18 NA12878 WES Coverage High Low High High Bits 32 32 12 12 Decimal bits 25 25 7 8 P 0.9998 0.9996 0.9717 0.9993 S 0.9998 0.9996 0.9996 0.9996 C. PiM core initialization We pre-calculate the log of the initial values and constants in the host CPU since the log function does not exist in the PiM core runtime library. These log values are then transferred in a fixed-point format directly to the PiM core. For instance, the initialization value for \({D}_{i,0}\) is sent to the PiM core as \(fixed\_point\left(log\right(\frac{1}{m}\left)\right)\) , and the transition values are sent as the fixed point of their log values. Implementation on PiM cores A. CPU-PiM data transfers The data of all the regions assigned to a particular PiM core is transferred to the PiM core’s MRAM. It includes the reads of the regions and their auxiliary information (read bases, length, and log-domain transitions probabilities), as well as the haplotypes of the regions and their information (read sequence, length and log-domain initial value ( fixed_point(log(1/m)) as detailed in Section V.C. )). B. PiM kernel Each PiM thread is assigned a read; then each PiM thread calculates the likelihood of such read against each haplotype in the region. Whenever a read from a new region is assigned to the PiM core, the haplotypes for that region are transferred to the WRAM and remain until the processing for the region is completed. This approach allows us to effectively exploit data reuse, as the haplotypes can be accessed multiple times during the region processing. Each running PiM thread fetches (function reserve_read in Algorithm 1 ) a read from a pool of reads (either at the start of the program or when it completes handling its current read) and then proceeds to transfer the relevant data associated with that read to the WRAM to compute the likelihood result of the associated read against each haplotype of the region. To manage the read allocation to PiM threads, we maintain a global counter that is shared among all PiM threads, and use a mutex to prevent simultaneous allocation of reads to multiple PiM threads, as explained above (part IV.C). A simplified pseudo-code run by each PiM thread is presented in Algorithm 1. Algorithm 1. Pseudo code for PiM core 1: function allocate_read_for_thread 2: mram_read(read length) 3: mram_read(read sequence) 4: mram_read(read transitions) 5: function reserve_read 6: mutex_lock() 7: result = free_read_index++ 8: if region is new then 9: mram_read(region haplotypes) 10: mutex_unlock 11: return result 12: function main 13: while reserve_read < reads number do 14: allocate_read_for_thread 15: for each haplotype in region do 16: initialize matrices 17: for each base in haplotype do 18: for each base in read do 19: update matrices 20: sum likelihood 21: mram_write(result likelihood) 22: return C. Efficient handling of on-chip memory To optimize the memory usage, the haplotypes present in the WRAM (that have been transferred from the MRAM) are stored in a circular buffer. The utilization of the circular buffer allows for a continuous reuse of memory by replacing the oldest data with new data, preventing wastage, and making efficient use of the available memory. Next, we transfer the read information from the MRAM to the WRAM, including the length, sequence, and transitions probabilities (respectively lines 2,3,4 in Algorithm 1 ). In order to fill the matrices efficiently and optimize memory usage, we only keep two lines of each matrix (for each PiM thread). This is because, as presented in Eq. ( 2 ), the value in each matrix cell depends only on its neighbors as illustrated in Fig. 5 . The intermediate likelihood result for each read-haplotype pair is stored in the WRAM. After completing the pairwise alignment of the read and the haplotype, the resulting likelihood is written to the MRAM memory in a result matrix (line 21 in Algorithm 1 ). In order to maximize the bandwidth of the WRAM-MRAM transfer, we coalesce writes from several PiM threads. D. Load balancing over PiM cores To optimize the performance of the parallel PiM cores, it is important to ensure that the workload is evenly distributed among them. One alternative approach to the method we employ in this work involves dividing genome regions using a SNP database, as reported in [38], assuming that the number of variants in the sample is comparable to the number of variants in the database. Our approach is to allocate several regions to each PiM core based on their complexity, which is determined by the total length of reads and haplotypes in the region. The load balancing is done by the host CPU that evaluates the complexity of each region and potentially divides regions into smaller fragments to meet a specified complexity target. E. Asynchronous implementation E. As presented in Section V, the pair-HMM algorithm requires redesign to run efficiently on PiM cores, including pre-processing and post processing steps running on the host CPU, such as splitting the regions. PiM cores can be launched asynchronously, i.e., we can launch them by rank of 64 PiM cores and each rank can be launched independantly of the others. We leverage this to minimize the load imbalance among PiM cores ranks. To efficiently utilize available ressources and optimize PiM cores work we construct a pipeline as illustrated in Fig. 6 below. To link between the different stages of this pipeline we use buffer queues (one input queue and one output queue). This pipeline is executed on the host CPU and contains three stages: Pre-processing stage handled by one single thread that writes to the input queue. PiM cores populating and launching stage . In this stage, a thread is created for each PiM core rank and reads from the input queue to transfer data to PiM cores and writes to the output queue the final result produced by the PiM cores. • Post processing stage handled by a single thread that reads from the output queue and writes the result to an output file. By constructing the pipeline in this manner, we are able to achieve latency hiding during the pre-processing, data transfer, and post-processing steps. This is because the time spent on PiM core computation is much greater than the time spent on these tasks. Thus, when computing the running time on future architectures (as discussed in VII) we disregard the impact of non PiM core time, as explained in [53]. The profiler tool provided by UPMEM [59] illustrates the pipeline parallelization as shown in Fig. 7 , with two levels of parallelism: the first is inter-rank parallelism, meaning that several ranks run in parallel (green segments represents PiM core time), and the second is thread parallelism on the host side, meaning that pre-processing task (blue segments), data transfer task (red segments), and post-processing tasks (orange segments) are performed concurrently. Evaluation A. Evaluation Methodology The algorithm is implemented on servers provided by UPMEM. The architecture is as follows: Host: Intel Xeon Silver 4215 CPU, 2.5GHz 2560 PiM cores divided into 40 ranks. Each PiM core operates at a frequency of 400MHz. The benchmarks used in this paper are obtained using a publicly available Whole Genome Sequencing dataset of NA12878 from the International Genome Sample Ressource (IGSR) [46]. We use real datasets in order to ensure transparency and accurately reflect the performance of the Pair-HMM algorithm in real-world scenarios. We present the results for chromosomes 1, 2 and 18 as well as for WGS dataset. We compare the processing time of our solution with different CPU implementations, including the original Java implementation, the AVX-accelerated version, and the OpenMP multithreaded version (with 4 threads and 8 threads). Additionally, we present the time results of the open-source FPGA implementation from [34], which includes two versions, one with 24 “workers” and the second one with 96 “workers”, where a worker is an acceleration unit (respectively represented as 24wk and 96wk). This Pair-HMM implementation is run on a F1 instance of AWS [54]. B. Results The GATK execution time figures are presented in Table 3 below along with the speedup figures, using the original Java implementation [22] of the algorithm as the baseline. The results indicate that the existing UPMEM architecture (2560 PiM cores operating at 400MHz, indicated by bold numbers) can provide speedups similar to that achieved by CPU or FPGA. For example, the UPMEM platform is able to process chromosome 1 in 117 seconds, which is slightly longer than 100 seconds it took the CPU with OpenMP 4 threads, and slightly faster than the 124 seconds needed for the FPGA implementation. The performance achieved by this existing UPMEM architecture represents up to 10.7x speedup compared to the original Java implementation of the algorithm. UPMEM has laid plans for continued improvement and innovation of its product offerings, and these results serve as a prediction of what can be expected from these future architectures. We present in Figure 8 below the performance roadmap [60] of the UPMEM architecture as inspired by [49]. Each percent figure shows the relative improvement of DPU’s frequency in the related time period. A scaled-up UPMEM platform with 2560 PiM cores running at 550MHz (and above) surpasses both the CPU with OpenMP 4 threads and the FPGA based platforms. This can be seen in the speedup achieved for chromosomes 1 and 2, which are 14.8x and 12.7x respectively compared to 12.5x and 10.6x for CPU (OpenMP 4t) and to 10x and 11.6x for the FPGA speedup respectively. This is a promising result as UPMEM plans on releasing a version of their architecture with PiM cores working at 550MHz in Q2 2023. Table 3 Running time in seconds and acceleration factor of the Pair-HMM algorithm on different architectures CHR1 CHR2 CHR18 WGS Java (reference) 1254.8 1x 638.7 1x 162 1x 15159.6 1x AVX 317.1 4x 163 3.9x 40.2 4x 3759 4x OpenMP 4 threads 100.5 12.5x 60 10.6x 15.8 10.2x 1189.4 12.7x OpenMP 8 threads 75.3 16.7x 46.9 13.6x 12.9 12.6x 886.7 17x OpenMP 4t FPGA 24wk 140.8 8.9x 57.6 11x 12.3 13.2x 1298 11.7x OpenMP 4t FPGA 96wk 124.4 10x 55.2 11.6x 12.2 13.3x 1200 12.6x UPMEM 2554 PiM cores 350MHz 134 9.4x 79 8x 23 7x 1637 9.2x UPMEM 2560 PiM cores 400MHz 117 10.7x 69 9.3x 20.1 8x 1429 10.6x UPMEM 2560 PiM cores 550MHz 85 14.8x 50.2 12.7x 14.6 11.1x 1039 14.6x UPMEM 3584 PiM cores 550MHz 60.8 20.6x 35.8 17.8x 10.4 15.6x 742.3 20.4x UPMEM 3584 PiM cores 600MHz 55.7 22.5x 32.8 19.5x 9.6 16.9x 680.5 22.3x UPMEM 5120 PiM cores 600MHz 39 32.2x 23 27.8x 6.7 24.2x 476.3 31.8x Figure 9 presents (a) the comparison between CPUs execution time (OpenMP four threads and OpenMP eight threads) and PiM cores execution time and (b) the comparison between FPGAs execution time and PiM cores execution time over different architectures for the chromosome 1: C. Power consumption and energy The power consumption of the PiM core depends on its version and clock speed. The actual consumption of a PiM core is not yet available on the current chip version but the PiM core’s TDP (Thermal Design Power, i.e., the measure of the maximum amount of power that the PiM core consumes) is provided by the manufacturer. When running at 350MHz, the TDP of the new version of the PiM core is 112mW, whereas at 400MHz it consumes 220mW [60]. PiM cores that run at 550MHz are designed to consume less energy, and their TDP is expected to be around 200mW. We compute the overall TDP of the PiM system used for testing, (according to the values provided by the manufacturer [60]) composed of 2560 PiM cores operating at 350MHz. We recall that the UPMEM chip is composed of eight PiM cores and thus this system contains 320 chips. \(TDP\left( PiM core \right)=112mW\) \(TDP\left( One UPMEM chip\right)=8\text{*}TDP\left( PiM core \right)=0.9W\) \(TDP\left( System \right)=320\text{*}0.9=288W\) D. Summary analysis To summarize this analysis, we compare the speedup achieved by each studied architecture as a function of the number of acceleration units (i.e., the level of parallelism), as presented in Fig. 10 below. The X axis of the graph is the number of acceleration units (as specifically defined in each architecture) and the Y axis of the graph represents the speedup factor for the corresponding accelerator. The PiM UPMEM x-axis represents a projection of performances on future architectures, having as a basis the current existing PiM system (2560 DPUs working at 400MHz). The speedups of OpenMP and FPGA accelerators increase very mildly as the number of acceleration units grows. In contrast, the UPMEM accelerator’s performance scales practically linearly. Indeed, the OpenMP multithreading implementation is limited and suffers from high Cycles Per Instruction (CPI) leading to a low performance increase. Similarly, FPGAs suffer from limited on-chip memory and their high hardware resource consumption makes it impractical to add more workers units. E. Discussion UPMEM PiM platform based Pair-HMM algorithm outperforms state-of-the-art solutions in terms of performance and power consumption. Based on UPMEM performance roadmap, such performance and energy efficiency gaps are expected to grow in the near future. This is because our implementation operates UPMEM DPU at very modest operating frequency which leaves an ample opportunity for further improvements. At the same time, we cannot expect the operating frequency of high-performance CPU and FPGA which are already fully optimized, to improve at similar pace. Implementing the GATK on UPMEM platforms carries an additional unique advantage. Genome analysis pipeline, which GATK is a major part of, comprises indexing, read mapping, read alignment and in some cases, genome assembly (reference guided or de-novo). Most of these operations are data-intensive and therefore strongly benefit from PiM implementation [10,11,12,13,14,15]. If indexing, read mapping and other data-intensive operations were to be on a PiM platform, while maintaining the conventional CPU-centric (CPU or FPGA based) GATK implementation, it would lead to several significant drawbacks. First, the overall performance would be adversely affected due to the necessity to transfer extremely large amount of sequenced data from PiM platform to a CPU or a FPGA (genomic read file can reach more than 100GB for a single human genome). Second, the total cost of ownership would be significantly higher. Indeed, a heterogeneous solution with two different processing clusters (PiM and FPGA for example) induces higher capital expenditure (CAPEX) due to the cost of the different hardware, as well as higher operational expenditure (OPEX) due to the need for maintenance and increased energy consumption. The ability to implement the entire genome analysis pipeline, from the read mapping to the GATK on the same PiM platform resolves both these inefficiencies. As demonstrated in this paper, to fully realize the benefits of PiM architecture in genome analysis, it is necessary to develop algorithms that are optimized for these architectures, including the parallelization approach, genome region assignment, computing in log domain and so on. Conclusions This work implements and evaluates Pair-HMM, the main part of the GATK variant calling algorithm and the most computationally intensive part of genome analysis pipeline, on a commercial PiM platform developed by UPMEM. We redesign the original algorithm to optimize it for the UPMEM PiM platform. The suggested modifications include moving the computation to log domain and replacing floating point calculations by fixed point one with very insignificant loss of the variant calling accuracy. To optimize the UPMEM platform utilization and performance, we develop a scheme for optimal allocation of genome fragments to individual PiM cores. This work further demonstrates the potential of PiM architectures as a future solution for end-to-end genome analysis. Despite some parts of the process not being memory-bound, the results show that even complex algorithms such as Pair-HMM can be effectively implemented on PiM architecture and provide speedup over state-of-the-art FPGA or CPU solution. Declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Availability of data and materials The code of this research is available at https://github.com/naomieab/UPMEM_HAPLOTYPECALLER Competing interests None to declare. Funding This work was supported by European Union’s Horizon Europe programme for research and innovation (grant Number 101047160). Authors’ contributions The code of this research has been written by Naomie Abecassis and Aphelie Moisson-Franckhauser. The research has been supervised and led by Dr. Leonid Yavits and Prof. Ran Ginosar and has also been led in partnership with Dr. Juan Gómez-Luna and Prof. Onur Mutlu. All authors reviewed the manuscript. References Pettersson E, Lundeberg J, Ahmadian A. Generations of sequencing technologies. Genomics. 2009 Feb 1;93(2):105-11. Reis-Filho JS. Next-generation sequencing. Breast cancer research. 2009 Dec;11(3):1-7. Mardis ER. DNA sequencing technologies: 2006–2016. Nature protocols. 2017 Feb;12(2):213-8. Roodi M, Moshovos A. Gene sequencing: where time goes. In2018 IEEE International Symposium on Workload Characterization (IISWC) 2018 Sep 30 (pp. 84-85). IEEE. Stephens ZD, Lee SY, Faghri F, Campbell RH, Zhai C, Efron MJ, Iyer R, Schatz MC, Sinha S, Robinson GE. Big data: astronomical or genomical?. PLoS biology. 2015 Jul 7;13(7):e1002195. Mutlu O, Ghose S, Gómez-Luna J, Ausavarungnirun R. Processing data where it makes sense: Enabling in-memory computation. Microprocessors and Microsystems. 2019 Jun 1;67:28-41. Ghose S, Boroumand A, Kim JS, Gómez-Luna J, Mutlu O. Processing-in-memory: A workload-driven perspective. IBM Journal of Research and Development. 2019 Aug 8;63(6):3-1. Kim VS, Fallin C, Lee D, Ausavarungnirun R, Luo GP, Mutlu O, Gibbons PB, Kozuch MA, Mowry TC. RowClone: Fast and efficient In-DRAM copy and initialization of bulk data. Microarchitecture. Mutlu O. Memory scaling: A systems architecture perspective. In2013 5th IEEE International Memory Workshop 2013 May 26 (pp. 21-25). IEEE. Khalifa M, Ben-Hur R, Ronen R, Leitersdorf O, Yavits L, Kvatinsky S. FiltPIM: In-memory filter for DNA sequencing. In2021 28th IEEE International Conference on Electronics, Circuits, and Systems (ICECS) 2021 Nov 28 (pp. 1-4). IEEE. Kaplan R, Yavits L, Ginosasr R. BioSEAL: In-memory biological sequence alignment accelerator for large-scale genomic data. InProceedings of the 13th ACM International Systems and Storage Conference 2020 May 30 (pp. 36-48). Angizi S, Fahmi NA, Zhang W, Fan D. Pim-assembler: A processing-in-memory platform for genome assembly. In2020 57th ACM/IEEE design automation conference (DAC) 2020 Jul 20 (pp. 1-6). IEEE. Lavenier D, Cimadomo R, Jodin R. Variant calling parallelization on processor-in-memory architecture. In2020 IEEE International Conference on Bioinformatics and Biomedicine (BIBM) 2020 Dec 16 (pp. 204-207). IEEE. Diab S, Nassereldine A, Alser M, Gómez Luna J, Mutlu O, El Hajj I. A framework for high-throughput sequence alignment using real processing-in-memory systems. Bioinformatics. 2023 May 1;39(5):btad155. Diab S, Nassereldine A, Alser M, Luna JG, Mutlu O, Hajj IE. High-throughput pairwise alignment with the wavefront algorithm using processing-in-memory. arXiv preprint arXiv:2204.02085. 2022 Apr 5. https://gatk.broadinstitute.org/ (2023) Shastry BS. SNP alleles in human disease and evolution. Journal of human genetics. 2002 Nov;47(11):561-6. Engle LJ, Simpson CL, Landers JE. Using high-throughput SNP technologies to study cancer. Oncogene. 2006 Mar;25(11):1594-601. Mullaney JM, Mills RE, Pittard WS, Devine SE. Small insertions and deletions (INDELs) in human genomes. Human molecular genetics. 2010 Oct 15;19(R2):R131-6. Bernig T, Chanock SJ. Challenges of SNP genotyping and genetic variation: its future role in diagnosis and treatment of cancer. Expert review of molecular diagnostics. 2006 May 1;6(3):319-31. Zhao S, Agafonov O, Azab A, Stokowy T, Hovig E. Accuracy and efficiency of germline variant calling pipelines for human genome data. Scientific reports. 2020 Nov 19;10(1):20222. https://github.com/broadinstitute/gatk (2023) Al-Ars Z, Wang S, Mushtaq H. SparkRA: enabling big data scalability for the GATK RNA-seq Pipeline with Apache Spark. Genes. 2020 Jan 3;11(1):53. Benjamin, D. (2018). Local Assembly in HaplotypeCaller and Mutect. Available online: https://github.com/broadinstitute/gatk/blob/master/docs/local_assembly.pdf Benjamin, D. (2018). Pair HMM probabilistic realignment in HaplotypeCaller and Mutect. Broad Institute. Available online: https://github.com/broadinstitute/gatk/blob/master/docs/pair_hmm.pdf Schapranow MP, Plattner H. HIG—An in-memory database platform enabling real-time analyses of genome data. In2013 IEEE International Conference on Big Data 2013 Oct 6 (pp. 691-696). IEEE. Firnkorn D, Knaup-Gregori P, Lorenzo Bermejo J, Ganzinger M. Alignment of high-throughput sequencing data inside in-memory databases. Ine-Health–For Continuity of Care 2014 (pp. 476-480). IOS Press. Kaplan R, Yavits L, Ginosar R. RASSA: resistive prealignment accelerator for approximate DNA long read mapping. IEEE Micro. 2018 Dec 28;39(4):44-54. Kim JS, Senol Cali D, Xin H, Lee D, Ghose S, Alser M, Hassan H, Ergin O, Alkan C, Mutlu O. GRIM-Filter: Fast seed location filtering in DNA read mapping using processing-in-memory technologies. BMC genomics. 2018 May;19(2):23-40. Foley P, Prabhakaran A, Gururaj K, Naik M, Gopalan S, Shargorodskiy A, Brau E. Accelerate Genomics Research with the Broad-Intel Genomics Stack. Rauer C, Finamore N. Accelerating genomics research with opencl and fpgas. Altera, Now Part of Intel, Tech. Rep. 2016 Mar. Ren S, Sima VM, Al-Ars Z. FPGA acceleration of the pair-HMMs forward algorithm for DNA sequence analysis. In2015 IEEE International Conference on Bioinformatics and Biomedicine (BIBM) 2015 Nov 9 (pp. 1465-1470). IEEE. Sampietro D, Crippa C, Di Tucci L, Del Sozzo E, Santambrogio MD. Fpga-based pairhmm forward algorithm for dna variant calling. In2018 IEEE 29th International Conference on Application-specific Systems, Architectures and Processors (ASAP) 2018 Jul 10 (pp. 1-8). IEEE. Wertenbroek R, Thoma Y. Acceleration of the Pair-HMM forward algorithm on FPGA with cloud integration for GATK. In2019 IEEE International Conference on Bioinformatics and Biomedicine (BIBM) 2019 Nov 18 (pp. 534-541). IEEE. Code available at: https://github.com/rick-heig/PHMM-F1 Ren S, Bertels K, Al-Ars Z. Efficient acceleration of the pair-hmms forward algorithm for gatk haplotypecaller on graphics processing units. Evolutionary Bioinformatics. 2018 Mar;14:1176934318760543. Ren S, Ahmed N, Bertels K, Al-Ars Z. GPU accelerated sequence alignment with traceback for GATK HaplotypeCaller. BMC genomics. 2019 Apr;20(2):103-16. Wu X, Subramaniyan A, Wang Z, Narayanasamy S, Das R, Blaauw D. A High-Throughput Pruning-Based Pair-Hidden-Markov-Model Hardware Accelerator for Next-Generation DNA Sequencing. IEEE Solid-State Circuits Letters. 2020 Dec 16;4:31-5. Ren S, Bertels K, Al-Ars Z. GPU-accelerated GATK haplotypecaller with load-balanced multi-process optimization. In2017 IEEE 17th International Conference on Bioinformatics and Bioengineering (BIBE) 2017 Oct 23 (pp. 497-502). IEEE. UPMEM Processing In-Memory (PIM), technology paper, unpublished results Gómez-Luna J, El Hajj I, Fernandez I, Giannoula C, Oliveira GF, Mutlu O. Benchmarking a new paradigm: Experimental analysis and characterization of a real processing-in-memory system. IEEE Access. 2022 May 10;10:52565-608. Gómez-Luna J, El Hajj I, Fernandez I, Giannoula C, Oliveira GF, Mutlu O. Benchmarking memory-centric computing systems: Analysis of real processing-in-memory hardware. In2021 12th International Green and Sustainable Computing Conference (IGSC) 2021 Oct 18 (pp. 1-7). IEEE. UPMEM, Introduction to UPMEM PIM. Processing-in-memory (PIM) on DRAM Accelerator (White Paper). 2018. Devaux F. The true processing in memory accelerator. In2019 IEEE Hot Chips 31 Symposium (HCS) 2019 Aug 1 (pp. 1-24). IEEE Computer Society. UPMEM, UPMEM User Manual. Version 2021.1.0. 2021. LLVM. Compiler-RT, LLVM project. https://github.com/llvm/llvm-project/tree/main/compiler-rt/lib/builtins, 2021. https://www.internationalgenome.org/ (2023) Rabbani B, Tekin M, Mahdieh N. The promise of whole-exome sequencing in medical genetics. Journal of human genetics. 2014 Jan;59(1):5-15. https://www.mathworks.com/products/fixed-point-designer.html (2023) Hennessy JL, Patterson DA. A new golden age for computer architecture. Communications of the ACM. 2019 Jan 28;62(2):48-60. Cleary JG, Braithwaite R, Gaastra K, Hilbush BS, Inglis S, Irvine SA, Jackson A, Littin R, Rathod M, Ware D, Zook JM. Comparing variant call files for performance benchmarking of next-generation sequencing variant calling pipelines. BioRxiv. 2015 Aug 2:023754. Mansouri Ghiasi N, Park J, Mustafa H, Kim J, Olgun A, Gollwitzer A, Senol Cali D, Firtina C, Mao H, Almadhoun Alserr N, Ausavarungnirun R. GenStore: a high-performance in-storage processing system for genome sequence analysis. InProceedings of the 27th ACM International Conference on Architectural Support for Programming Languages and Operating Systems 2022 Feb 28 (pp. 635-654). Molka D, Schöne R, Hackenberg D, Nagel WE. Detecting memory-boundedness with hardware performance counters. InProceedings of the 8th ACM/SPEC on International Conference on Performance Engineering 2017 Apr 17 (pp. 27-38). Yavits L, Morad A, Ginosar R. The effect of communication and synchronization on Amdahl’s law in multicore systems. Parallel Computing. 2014 Jan 1;40(1):1-6. Amazon EC2 F1 Instances. Available: https://aws.amazon.com/ec2/instance-types/f1/ (2023) https://github.com/aayasin/perf-tools (2023) Intel Corporation (December, 2017), Intel® VTune™ Amplifier. https://software.intel.com/en-us/intel-vtune-amplifier-xe (2023) Gómez-Luna J, Guo Y, Oliveira GF, Sadrosadati M, Mutlu O. TransPimLib: A Library for Efficient Transcendental Functions on Processing-in-Memory Systems. arXiv preprint arXiv:2304.01951. 2023 Apr 3. Poplin R, Chang PC, Alexander D, Schwartz S, Colthurst T, Ku A, Newburger D, Dijamco J, Nguyen N, Afshar PT, Gross SS. A universal SNP and small-indel variant caller using deep neural networks. Nature biotechnology. 2018 Nov;36(10):983-7. Available at https://github.com/google/deepvariant https://sdk.upmem.com/2023.2.0/261_Application_profiling.html#application-profiling (2023) Y. Falevoz. UPMEM. Personal Communication. October 2022. Additional Declarations No competing interests reported. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-3245756","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":225395708,"identity":"b895bd59-9c34-4120-a9b4-02fad2ecc9d3","order_by":0,"name":"Naomie Abecassis","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABH0lEQVRIiWNgGAWjYLCCBAYGOSQuD1gEDA48wK3FGMaWQNWSgE09BCQ2oGhBNg0dmLMfPrrhQYVden/76bTHPH8Y6vjZzx578KDiDgN/+wFGbLZY9qSl3Ug4k5w740zudmPeNgYJyZ68dIOEM88YJM4kYHWYwYEcsxuJbcy5Gxhyt0nzNjBIGNzgMZNIbDvMwHADu18Mzr8BavlXn27A/3abNNBhEvYwLfK4tNwA2dJwOMFAAmgLDxvQFgmoFgOcWp4B/XLsuOGMG2+3G85tk5CccSbHTCLhzGEewzOJDdgdlnzs5o+aann+/txtD978seHnbz9jJvmj4rCc3PHDhz98wNSCDNgYwNECBcDYYWzArwGiZRSMglEwCkYBJgAANptr04b2GX0AAAAASUVORK5CYII=","orcid":"","institution":"Technion Israel Institute of Technology","correspondingAuthor":true,"submittingAuthor":false,"prefix":"","firstName":"Naomie","middleName":"","lastName":"Abecassis","suffix":""},{"id":225395709,"identity":"ad9653f6-aac0-4f31-a93d-24083ed1b7b9","order_by":1,"name":"Juan Gómez-Luna","email":"","orcid":"","institution":"ETH Zurich Zurich","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Juan","middleName":"","lastName":"Gómez-Luna","suffix":""},{"id":225395710,"identity":"3a165394-6c9f-4b8e-9698-55776570d647","order_by":2,"name":"Onur Mutlu","email":"","orcid":"","institution":"ETH Zurich Zurich","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Onur","middleName":"","lastName":"Mutlu","suffix":""},{"id":225395711,"identity":"6a984e0c-77a2-48f1-9a73-4d6326defd47","order_by":3,"name":"Ran Ginosar","email":"","orcid":"","institution":"Technion Israel Institute of Technology","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Ran","middleName":"","lastName":"Ginosar","suffix":""},{"id":225395712,"identity":"3b9dba44-d410-413a-b3d3-9e3d3cb65838","order_by":4,"name":"Aphélie Moisson-Franckhauser","email":"","orcid":"","institution":"UPMEM","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Aphélie","middleName":"","lastName":"Moisson-Franckhauser","suffix":""},{"id":225395713,"identity":"0c9de710-675c-44e1-8609-c4bb5aeae462","order_by":5,"name":"Leonid Yavits","email":"","orcid":"","institution":"Bar-Ilan University","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Leonid","middleName":"","lastName":"Yavits","suffix":""}],"badges":[],"createdAt":"2023-08-08 13:44:23","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-3245756/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-3245756/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":41785960,"identity":"4135ef13-ae1b-456e-a55d-07184ce7e1b8","added_by":"auto","created_at":"2023-08-18 20:42:28","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":32083,"visible":true,"origin":"","legend":"\u003cp\u003eExecution time breakdown of GATK HaplotypeCaller\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/4d23ed18029ad3fca25e8ae6.png"},{"id":41787574,"identity":"38dc6630-7445-47b1-a72f-d65a10d1414b","added_by":"auto","created_at":"2023-08-18 20:58:28","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":21352,"visible":true,"origin":"","legend":"\u003cp\u003eFinite state machine model for HaplotypeCaller\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/4a5693ad6cee4fbe265318c1.png"},{"id":41786679,"identity":"ab491c8e-1677-41a8-a954-27f860499c1a","added_by":"auto","created_at":"2023-08-18 20:50:28","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":55516,"visible":true,"origin":"","legend":"\u003cp\u003eUPMEM PiM architecture [39]\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/d47e54b963450c35df331541.png"},{"id":41785968,"identity":"d8553412-cf8f-4425-859d-e859fb7c1b77","added_by":"auto","created_at":"2023-08-18 20:42:28","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":84061,"visible":true,"origin":"","legend":"\u003cp\u003eUPMEM chip diagram (reproduced from [38])\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/21069ffdd2173593db13b305.png"},{"id":41785961,"identity":"c8ff90e1-9500-48d7-9f3a-b021cf6bde24","added_by":"auto","created_at":"2023-08-18 20:42:28","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":31336,"visible":true,"origin":"","legend":"\u003cp\u003eCell computation dependencies in the matrix\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/c678f3a14a1d2fb3136cf3ea.png"},{"id":41786680,"identity":"a871de4b-08df-466c-b169-fc412abb5b4d","added_by":"auto","created_at":"2023-08-18 20:50:28","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":12993,"visible":true,"origin":"","legend":"\u003cp\u003eHost CPU pipeline stages\u003c/p\u003e","description":"","filename":"6.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/78c2d204e73783f7a394a2f6.png"},{"id":41787575,"identity":"8bcc377c-c521-4dc2-b3ba-a67b317b481c","added_by":"auto","created_at":"2023-08-18 20:58:28","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":50830,"visible":true,"origin":"","legend":"\u003cp\u003eTime diagram of the system pipeline\u003c/p\u003e","description":"","filename":"7.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/249a3ea4306836ad587e718d.png"},{"id":41785969,"identity":"d4b2adc4-17f9-44fc-8092-f3a21716ad1f","added_by":"auto","created_at":"2023-08-18 20:42:28","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":15328,"visible":true,"origin":"","legend":"\u003cp\u003ePerformance roadmap of UPMEM architecture\u003c/p\u003e","description":"","filename":"8.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/37f45ffcad15f1e1fa368582.png"},{"id":41785962,"identity":"8dd0744c-f2d1-4f8b-8685-f193c67cc660","added_by":"auto","created_at":"2023-08-18 20:42:28","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":19077,"visible":true,"origin":"","legend":"\u003cp\u003eRunning time comparison between architectures\u003c/p\u003e","description":"","filename":"9.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/e87286ec4a6fcaefe3eaf1f1.png"},{"id":41785965,"identity":"7d8b61c0-1cd3-4a88-8060-ddf2bf636ea3","added_by":"auto","created_at":"2023-08-18 20:42:28","extension":"png","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":40090,"visible":true,"origin":"","legend":"\u003cp\u003eSpeedup factor graph for different architectures\u003c/p\u003e","description":"","filename":"10.png","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/dfe893261e4c9ef3c326fcbe.png"},{"id":44346185,"identity":"de8cd29a-680e-45ed-9bd1-bd35273b94c8","added_by":"auto","created_at":"2023-10-10 08:37:36","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":786439,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3245756/v1/7ee56c9c-9253-4c72-a7a8-c50b67ddf94b.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"GAPiM: Discovering Genetic Variations on a Real Processing-in-Memory System","fulltext":[{"header":"Introduction","content":"\u003cp\u003eThe completion of the Human Genome Project in 2003 and the advent of next-generation sequencing (NGS) technology in the early 2000s marked the beginning of a new era in DNA sequencing [1,2], enabling faster and more efficient sequencing of large amounts of DNA, and driving the growth in the field of genomics. Over the years, the progress made in the field has led to a rapid increase in the amount of available genomic data and a decrease in the cost of sequencing [3].\u003c/p\u003e\n\u003cp\u003eGenomic data is analyzed to understand the genetic information carried by a DNA molecule and to identify DNA alterations that can cause diseases, such as genetic disease and cancer. The process of identifying variations in the sequenced data compared to a reference genome is known as variant calling and can take up to 40% of the total time spent in the genome sequencing pipeline [4]. Variant calling is an essential step in DNA sequencing pipelines because it enables the identification of differences, or variations, between the DNA sequence being analyzed and a reference genome or other sequences. These variations can include single nucleotide polymorphisms (SNPs, i.e., a DNA sequence variation that occurs when a single nucleotide (adenine, thymine, cytosine, or guanine) in the genome sequence is altered), small insertions or deletions (indels\u003csup\u003e2\u003c/sup\u003e), or larger structural variations such as copy number variations (CNVs, i.e., refer to the genetic trait involving the number of copies of a particular gene present in the genome of an individual) or genomic rearrangements.\u003c/p\u003e\n\u003cp\u003eVariant calling is important for a variety of applications, including disease diagnosis [18,20], personalized medicine, and understanding evolutionary relationships between organisms [19]. By identifying genetic variations, variant calling can help researchers and clinicians identify disease-causing mutations, develop targeted therapies, and better understand the genetic basis of disease. In addition, variant calling can be used to study the genetic diversity of populations or track the spread of infectious diseases.\u003c/p\u003e\n\u003cp\u003eWith the increased need for genomic data and continuous advancements in technology, it is predicted that the amount of genomic data generated will continue to rise [5].\u0026nbsp;Most of the computational algorithms developed for genome analysis are designed for use in a CPU system. The analysis of the sequenced data requires intensive data movement between DRAM and processing elements which are expensive, both in terms of time and energy [6,7,8,9].\u0026nbsp;In addition, the variant calling task is a highly parallelizable task as different regions in the genome can be analyzed independently of the others. \u0026nbsp;While CPUs,\u0026nbsp;FPGAs\u0026nbsp;and GPUs\u0026nbsp;are commonly used in variant calling to take advantage of parallel processing\u0026nbsp;[30,31,32,33,35,36], the speedup achieved by these hardware architectures is limited as explained in\u0026nbsp;the evaluation section, emphasizing the need for a new architecture paradigm.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eIn this paper we focus on Processing-in-Memory computer architecture\u0026nbsp;paradigm whose main purpose is alleviating\u0026nbsp;the \u003cem\u003ememory wall\u003c/em\u003e, which limits the performance and energy inefficiency of a conventional von Neumann architecture\u0026nbsp;for data intensive applications such as genome analysis.\u0026nbsp;Processing in Memory (PiM) architecture is a promising solution for genomic data analysis considering the expansion of the available data [10,11,12,13,14,15]. PiM architectures can solve data access bottleneck and lower the energy consumption.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003ePair-HMM is the most time-consuming part of the variant caller in the GATK pipeline. Thus, we investigate the PiM architecture as a potential solution of acceleration. PiM architecture provides high bandwidth and low latency memory access. Additionally, PiM provides a large amount of compute parallelism and enables a throughput that scales with memory size.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eIn this work, we show that a PiM architecture can efficiently implement an end-to-end solution for genomic analysis, even when certain components of the genome analysis pipeline are not highly data intensive.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eThe main contributions of this work are as follows:\u003c/p\u003e\n\u003col\u003e\n \u003cli\u003eWe rethink and adapt the Pair-HMM algorithm to make it suitable for a real Processing-in-Memory architecture.\u003c/li\u003e\n \u003cli\u003eWe implement the Pair-HMM algorithm on a commercial PiM system with more than 2500 PiM cores (DPUs).\u003c/li\u003e\n \u003cli\u003eWe comparatively evaluate the PiM Pair-HMM implementation and demonstrate up to 2x time performance improvement compared to high-performance CPU and up to 3x performance improvement compared to SOTA FPGA solutions.\u0026nbsp;\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Background","content":"\u003cp\u003eA. \u003cem\u003eThe Genome Analysis ToolKit (GATK)\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eGATK (Genome Analysis Toolkit) is a collection of tools developed by the Broad Institute to identify SNPs and indels in genomic data [16]. SNPs are the most common sort of genetic variation among human beings [17]. SNPs and indels play an important role in the study of human cancer genes [18,19,20]. The GATK best practices pipeline is one of the most widely used for variant calling, due to its high accuracy [21]. This pipeline is composed of three essential parts. The first one is data pre-processing, the second is variant discovery, and the third is variant evaluation. GATK is an open-source software available on GitHub [22].\u003c/p\u003e\n\u003cp\u003eB. \u003cem\u003eHaplotypeCaller algorithm\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe second part of the GATK pipeline, variant discovery, is implemented by the HaplotypeCaller algorithm. This algorithm gets as inputs a processed BAM file (a set of aligned DNA reads) and a reference sequence file, and outputs a VCF file (variant calling file). The HaplotypeCaller \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(1-\\varvec{\\epsilon }\\)\u003c/span\u003e\u003c/span\u003ealgorithm can take up to 35% of the complete genome analysis pipeline time as reported in [23]. The HaplotypeCaller algorithm has four main steps:\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eActive region determination\u003c/em\u003e based on the presence of a sufficient number of variations that exceed a predetermined threshold [16]. Active regions, i.e., with sufficient variations, are the target of the next three steps.\u003c/p\u003e\n\u003cp\u003e\u0026bull; \u003cem\u003eDe-Bruijn like graph assembly\u003c/em\u003e is used to assemble the reads of the region and determine a list of haplotypes [24]. A haplotype is a sequence that represents a certain combination of alleles (alternative gene forms) present in a region. Then the Smith-Waterman algorithm is used to align each haplotype to the reference haplotype to identify potential genetic variations.\u003c/p\u003e\n\u003cp\u003e\u0026bull; \u003cem\u003ePair-HMM forward algorithm\u003c/em\u003e is used to perform a pairwise alignment of each read against each haplotype [25]. The algorithm calculates the likelihood of each haplotype given each read by summing the likelihoods of all possible alignments of the read to the haplotype. This likelihood reflects the probability that the read and the haplotype are related, based on their sequence similarity.\u003c/p\u003e\n\u003cp\u003e\u0026bull; \u003cem\u003eGenotype sampling\u003c/em\u003e estimates the genotypes (combination of alleles) at each variant site in the region using Bayes\u0026rsquo; theorem in combination with the previously obtained likelihoods.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eC. Algorithm analysis\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe GATK HaplotypeCaller software includes a feature that tracks the time required to complete the Smith-Waterman and the Pair-HMM steps of the workflow. Additionally, we used profiling tools [55,56] and present the time repartition of the program in Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e (running on an Intel Core I7-5820K x86_64, 6 cores (12 threads), 3.30GHz, 64GB). The figure represents the average relative execution time for chromosomes 1\u0026ndash;3 of sample NA12878 (commonly used sample from a female individual of European descent). As observed, the Pair-HMM component of the algorithm accounts for a significant portion of 49% of the total time, prompting our decision to accelerate this section.\u003c/p\u003e\n\u003cp\u003eMoreover, we check the memory boundedness of the HaplotypeCaller algorithm by using hardware performance counters [52]. As reported in [52], the memory boundedness of the program, which is a component of the non-execution of the CPU (\u003cem\u003ecycle_activity.cycles_no_execute\u003c/em\u003e), can be computed by the following formula:\u003c/p\u003e\n\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e$$memory bound= \\text{m}\\text{a}\\text{x}(RESOURCE\\_STALLS:SB , CYCLE\\_ACTIVITY:STALLS\\_L1D\\_PENDING) \\left(1\\right)$$\u003c/div\u003e\n\u003c/div\u003e\n\u003cp\u003eWe present in Table \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e the results of this analysis for a single chromosome (chromosome 1 randomly chosen) and for a WGS running. As can be observed, the memory boundedness of the HaplotypeCaller algorithm is low, which means that memory access is unlikely to become a bottleneck of a HaplotypeCaller. The intuitive conclusion is that PiM would not be the most efficient solution for HaplotypeCaller. Nevertheless, there are three main advantages to PiM architecture even this non-memory-bound program can benefit from:\u003c/p\u003e\n\u003col\u003e\n \u003cli\u003ePiM architectures have a large number of execution units that provide large amounts of parallelism.\u003c/li\u003e\n \u003cli\u003ePiM compute throughput scales with memory capacity.\u003c/li\u003e\n \u003cli\u003ePiM systems are not limited in terms of memory capacity (or not as much) as other accelerators (e.g., GPUs or FPGAs with HBM memory).\u003c/li\u003e\n\u003c/ol\u003e\n\u003cp\u003e\u003c/p\u003e\u0026nbsp;\u003ctable id=\"Tab1\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\n \u003cdiv class=\"CaptionContent\"\u003e\n \u003cp\u003eMemory boundedness of HaplotypeCaller algorithm\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\u0026nbsp;\u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eChromosome 1\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eWGS\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eresource_stalls.sb\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e146,149,728,421\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e7,680,287,475,816\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003ecycle_activity.stalls_l1d_pending\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e791,177,991,761\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e65,333,865,516,882\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003ecycle_activity.cycles_no_execute\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e7,590,369,180,122\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e510,501,141,773,479\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eMemory boundedness\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e10%\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e12.8%\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n\u003c/table\u003e\n\u003cp\u003e\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eD. Pair-HMM forward algorithm\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eIn this work we focus on the Pair-HMM step which is the most resource-intensive part of the HaplotypeCaller algorithm. In this part we compute the likelihood of each read against each haplotype.\u003c/p\u003e\n\u003cp\u003eThe likelihood calculation is based on a pairwise alignment using a Hidden Markov Model (represented by a state machine) and the qualities of the read DNA bases (qualities produced by a sequencer). Given a DNA read R of length \u003cem\u003en\u003c/em\u003e and a haplotype H of length \u003cem\u003em\u003c/em\u003e we define the transitions of the state machine as the probabilities of the possible changes in state that occur during the process of aligning the read to the haplotype. In other words, the state machine represents the different states that the alignment can go through, such as a match (state M), an insertion (state I) or a deletion (state D) in the read. The transitions between the states are defined by the following probabilities:\u003c/p\u003e\n\u003cdiv id=\"Equb\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e$${\\delta } ={10}^{-\\frac{{\\text{Q}}_{\\text{d}\\text{e}\\text{l}}}{10}}, {\\iota } ={10}^{-\\frac{{\\text{Q}}_{\\text{i}\\text{n}\\text{s}}}{10}}, {\\epsilon } ={10}^{-\\frac{{\\text{Q}}_{\\text{g}\\text{c}\\text{p}}}{10}}$$\u003c/div\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Equ1\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equ1\" name=\"EquationSource\"\u003e$$prio{r}_{i,j}= \\left\\{\\begin{array}{c}1-{10}^{-\\frac{{Q}_{b}}{10}} if rea{d}_{i}=haplotyp{e}_{j}\\\\ \\frac{{10}^{-\\frac{{Q}_{b}}{10}}}{3} if rea{d}_{i}\\ne haplotyp{e}_{j}\\end{array}\\right.$$\u003c/div\u003e\n \u003cdiv class=\"EquationNumber\"\u003e2\u003c/div\u003e\n\u003c/div\u003e\n\u003cp\u003eCurrently, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\text{Q}}_{\\text{d}\\text{e}\\text{l}}\\)\u003c/span\u003e\u003c/span\u003e, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\text{Q}}_{\\text{i}\\text{n}\\text{s}}\\)\u003c/span\u003e\u003c/span\u003e and\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\text{Q}}_{\\text{g}\\text{c}\\text{p}}\\)\u003c/span\u003e\u003c/span\u003e (respectively the deletion, insertion, and gap continuation qualities) are included in the PHMM model but are not produced by sequencers and thus default values determined by GATK [22] are used for each base. However, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\text{Q}}_{\\text{b}}\\)\u003c/span\u003e\u003c/span\u003e is referred to as the quality of each base in the read and is provided by the sequencer machine.\u003c/p\u003e\n\u003cp\u003eThe finite state machine of \u003cem\u003eFig.\u0026nbsp;2\u003c/em\u003e models the transitions between the three different states (\u003cem\u003eM\u003c/em\u003e, the match state, \u003cem\u003eI\u003c/em\u003e, the insertion state and \u003cem\u003eD\u003c/em\u003e, the deletion state).\u003c/p\u003e\n\u003cp\u003eUsing this state machine with dynamic programming, the algorithm fills up three matrices, \u003cem\u003eM\u003c/em\u003e, \u003cem\u003eI\u003c/em\u003e and \u003cem\u003eD\u003c/em\u003e as shown in Equations \u003cem\u003e(3;1)\u003c/em\u003e and \u003cem\u003e(3;2)\u003c/em\u003e. Each element \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(i,j\\)\u003c/span\u003e\u003c/span\u003e of the matrix represents the total likelihood of all paths from the beginning of the read to position \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(i\\)\u003c/span\u003e\u003c/span\u003e such that the path ends with the corresponding state (match state for matrix \u003cem\u003eM\u003c/em\u003e, deletion state for matrix \u003cem\u003eD\u003c/em\u003e, and insertion state for matrix \u003cem\u003eI\u003c/em\u003e).\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eInitialization\u003c/em\u003e:\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\forall 0\\le i\u0026lt;n, 0\\le j\u0026lt;m,\\)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n\u003cdiv id=\"Equc\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equc\" name=\"EquationSource\"\u003e$$\\left\\{\\begin{array}{c}{M}_{i,0}=0 , {I}_{i,0}=0, {D}_{i,0}=\\frac{1}{m}\\\\ {M}_{0,j}=0 , {I}_{0,j}=0, {D}_{0,j}=0\\\\ {M}_{\\text{0,0}}=0 , {I}_{\\text{0,0}}=0, {D}_{\\text{0,0}}=0\\end{array}\\right. (3;1)$$\u003c/div\u003e\n\u003c/div\u003e\n\u003cp\u003e\u003cem\u003eRecursion step\u003c/em\u003e:\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\forall 1\\le i\\le n, 1\\le j\\le m,\\)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e\n\u003cp\u003e\u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(\\left\\{\\begin{array}{c}{M}_{i,j}=prio{r}_{i,j}*[\\left(1-\\left({\\iota }+{\\delta }\\right)\\right){ M}_{i-1,j-1}+\\\\ +\\left(1-{\\epsilon }\\right) {(\\text{I}}_{\\text{i}-1,\\text{j}-1} + {\\text{D}}_{\\text{i}-1,\\text{j}-1}\\left)\\right]\\\\ {I}_{i,j}= {{\\epsilon }\\text{*}\\text{I}}_{\\text{i}-1,\\text{j}} + {{\\iota }\\text{*}\\text{M}}_{\\text{i}-1,\\text{j}}\\\\ {D}_{i,j}={{\\epsilon }\\text{*}\\text{D}}_{\\text{i},\\text{j}-1} + {{\\delta }\\text{*}\\text{M}}_{\\text{i},\\text{j}-1}\\end{array}\\right.(3;2\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e)\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eTermination step\u003c/em\u003e:\u003c/p\u003e\n\u003cdiv id=\"Equd\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equd\" name=\"EquationSource\"\u003e$$L=\\sum _{j=1}^{m}{M}_{n,j}+{I}_{n,j} (3;3)$$\u003c/div\u003e\n\u003c/div\u003e\n\u003cp\u003eThe terminations step computes the final likelihood for the current read-haplotype.\u003c/p\u003e"},{"header":"Related work","content":"\u003cp\u003eGenome analysis requires high-performance computing to process and store large and continuously growing genomic datasets. Several processing in memory (PiM) solutions have been proposed to accelerate various parts of genome analysis pipeline. For example, the development of UpVC [13] on the UPMEM architecture has sped up the read alignment and variant calling process and yielded a time and energy improvement compared to GPU or FPGA implementations, using a simpler and less accurate algorithm than PHMM. Other works have focused on accelerating the read alignment part of the pipeline by using in-memory databases [26,27], or by using resistive content addressable memory (BioSEAL [11]). RASSA [28] accelerates read pre-alignment using approximate search-capable resistive content addressable memory. GRIM filter [29] integrates computation within a logic layer optimized to exploit 3D stacked memory layers to filter seed location. More recently GenStore [51] has implemented filters for in-memory reads for genome analysis tasks, and research conducted by [15] demonstrates that the wavefront algorithm achieves greater read alignment throughput on PiM architecture. In contrast to the majority of prior research, it is crucial to emphasize that our work is conducted on a real Processing-in-Memory architecture.\u003c/p\u003e \u003cp\u003eGAPiM is the first GATK pipeline implementation on a commercial PiM platform. To our knowledge, GATK has never been implemented on an actual PIM platform before. Typical GATK accelerators are based on CPU, FPGA, and GPU. Most of these accelerations implement the Pair-HMM step. CPU acceleration solutions use a combination of SIMD extensions such as AVX and AVX-512 and multithreading [30]. FPGA-based solutions mainly focus on accelerating the floating-point calculations of Pair-HMM. Most FPGA solutions use systolic arrays [31,32,33]. Another solution for FPGA is to create a specialized unit to execute the PairHMM command and replicate it to provide efficient parallelization [34]. But these FPGA solutions suffer from limited on-chip memory. Our solution, GAPiM, alleviates the memory wall which may limit the CPU and FPGA implementations and hence enables better performance and scalability. We believe PiM approaches are also beneficial to investigate for emerging deep learning based variant calling tools [58]. We leave the investigation of the evaluation, design, and analysis of such tools to future work.\u003c/p\u003e \u003cp\u003eGPU accelerations are also used since the Pair-HMM algorithm is easily parallelizable, focusing on intra-task and inter-task parallelism [35] and allowing to retrieve the optimal alignment [36]. Recently, an ASIC that replaces floating point multiplication by 20-bit addition in log domain employs bound checks to maintain correctness has been proposed in [37]. GAPiM provides a better scalability and flexibility when compared with GPU or ASIC implementations.\u003c/p\u003e"},{"header":"UPMEM PiM architecture","content":"\u003cp\u003eIn this work we utilize the PiM platform developed by UPMEM [39]. It is the first commercially available PiM solutions. It has been shown to be effective in speeding up memory-intensive workloads [13,40,41]. UPMEM platform is a standard DDR4-2400 DIMM with several PiM chips, each one composed of 8 parallel processors, called DPUs (DRAM Processing Units). Although this paper we refer to DPUs as PiM cores. The PiM cores are operated by a host CPU that sends the program to be executed and collects the processing results when PiM cores complete their processing tasks [39,42,43].\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eA. Chip level organization\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eWithin each memory chip (illustrated in \u003cem\u003eFigure 3\u003c/em\u003e), there are 8 asynchronous PiM cores that operate independently of one another. Each PiM core is associated with a 64MB DRAM bank called MRAM (Main RAM) that can also be accessed by the host CPU. Different PiM cores cannot communicate with each other directly. The host CPU can transfer data to/from MRAM banks of individual PiM cores. This design allows each computing unit to efficiently process its own fragment of the dataset. Calculations are performed in-situ, within each unit at a memory bandwidth of almost 1GB/s.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eB. The PiM core\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eA PiM core is a 32-bit RISC-based processor (with 64-bit capabilities) running at 400MHz. It is a multithreaded processor with up to 24 hardware threads, called PiM threads. The PiM core has a pipeline depth of 14 stages, however, only the last three stages of the pipeline can execute in parallel with stages of the next instruction in the same thread. Therefore, instructions from the same thread must be dispatched 11 cycles apart, requiring at least 11 threads to fully utilize the pipeline [39,42,43]. Each PiM core has an exclusive access to a 24KB IRAM (Instruction RAM), and to a 64KB WRAM (Working RAM) shared by all the threads, as well as the 64MB MRAM as illustrated in \u003cem\u003eFig.\u0026nbsp;4\u003c/em\u003e. The WRAM has a lower access latency than MRAM and can be accessible in one cycle. To allow several threads to access a same resource we use a mutual exclusion (mutex) mechanism, i.e., a synchronization object is used to control access to a shared resource and ensure that only one PiM thread can access that resource at a time, avoiding conflicts and data inconsistencies.\u003c/p\u003e\n\u003ch2\u003eC. The PiM core\u003c/h2\u003e\n\u003cp\u003eA PiM core is a 32-bit RISC-based processor (with 64-bit capabilities) running at 400MHz. It is a multithreaded processor with up to 24 hardware threads, called PiM threads.\u0026nbsp;The\u0026nbsp;PiM core\u0026nbsp;has a pipeline depth of 14 stages, however, only the last three stages of the pipeline can execute in parallel with stages of the next instruction in the same thread. Therefore, instructions from the same thread must be dispatched 11 cycles apart, requiring at least 11 threads to fully utilize the pipeline\u0026nbsp;[39,42,43]. Each PiM core has an exclusive access to a 24KB IRAM (Instruction RAM), and to a 64KB WRAM (Working RAM) shared by all the threads, as well as the 64MB MRAM as illustrated in \u003cem\u003eFigure 4\u003c/em\u003e. The WRAM has a lower access latency than MRAM and can be accessible in one cycle. To allow several threads to access a same resource we use a mutual exclusion (mutex) mechanism, i.e., a synchronization object is used to control access to a shared resource and ensure that only one PiM thread can access that resource at a time, avoiding conflicts and data inconsistencies.\u003c/p\u003e\n\u003cp\u003eIn this setup, programmers explicitly transfer data between MRAM and WRAM using DMA calls which must be 8-byte aligned and execute load and store operations on WRAM.\u003c/p\u003e"},{"header":"Adapting HaplotypeCaller algo-rithm to UPMEM PiM architecture","content":"\u003cp\u003e\u003cstrong\u003eWe modified the Pair-HMM algorithm to enable its efficient implementation on UPMEM PiM architecture while optimizing the hardware utilization.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eA. Proposed log domain fixed point model\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe need to redesign the Pair-HMM for PiM core stems from two observations, as follows: first, multiplications in PiM core are much slower than additions (approximately 10 MOPS for multiplication vs. 60 MOPS for 32-bit integer addition) [40]. Second, floating point operations are significantly less efficient on PiM cores compared to integer operations (approximately 5 MFLOPS for additions) [40]. This is because the PiM core does not implement floating point ALUs and these operations are emulated by the UPMEM runtime library in software [44, 45].\u003c/p\u003e\n\u003cp\u003eThe Pair-HMM algorithm requires a high number of multiplications (Equations \u003cem\u003e(3;1)\u003c/em\u003e and \u003cem\u003e(3;2)\u003c/em\u003e) and uses floating point calculations (see priors and transition probabilities). These calculations are significantly slower than integer addition/subtraction on the wimpy PiM cores: multiplication takes six times longer and floating-point addition takes twelve times longer than integer addition. To overcome this limitation, we first convert all equations to log domain, replacing multiplications with additions. For example, the first equation of Eq.\u0026nbsp;(3;\u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e) is transformed as follows, computing the value of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\text{l}\\text{o}\\text{g}\\left({M}_{i,j}\\right)\\)\u003c/span\u003e\u003c/span\u003e rather than \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({M}_{i,j}\\)\u003c/span\u003e\u003c/span\u003e:\u003c/p\u003e\n\u003cp\u003e\u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(\\text{log}\\left({M}_{i,j}\\right)=\\text{l}\\text{o}\\text{g}\\left(prio{r}_{i,j}*\\right[\\left(1-\\left({\\iota }+{\\delta }\\right)\\right){ M}_{i-1,j-1}+\\left(1-{\\epsilon }\\right) {(\\text{I}}_{\\text{i}-1,\\text{j}-1} + {\\text{D}}_{\\text{i}-1,\\text{j}-1}\\left)\\right] )\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(=\\text{log}\\left(prio{r}_{i,j}\\right)+\\text{l}\\text{o}\\text{g}(\\left(1-\\left({\\iota }+{\\delta }\\right)\\right){ M}_{i-1,j-1}+\\left(1-{\\epsilon }\\right) {(\\text{I}}_{\\text{i}-1,\\text{j}-1} + {\\text{D}}_{\\text{i}-1,\\text{j}-1}))\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(=\\text{log}\\left(prio{r}_{i,j}\\right)+\\text{log}\\left({10}^{\\text{log}\\left(\\left(1-\\left({\\iota }+{\\delta }\\right)\\right){ M}_{i-1,j-1}\\right)}+{10}^{\\text{l}\\text{o}\\text{g}(\\left(1-{\\epsilon }\\right) {(\\text{I}}_{\\text{i}-1,\\text{j}-1} + {\\text{D}}_{\\text{i}-1,\\text{j}-1}))}\\right)\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(=\\text{log}\\left(prio{r}_{i,j}\\right)+\\text{log}\\left({10}^{\\text{log}\\left(1-\\left({\\iota }+{\\delta }\\right)\\right)+\\text{log}\\left({ M}_{i-1,j-1}\\right)}+{10}^{\\text{log}\\left(1-{\\epsilon }\\right)+\\text{l}\\text{o}\\text{g}{(\\text{I}}_{\\text{i}-1,\\text{j}-1} + {\\text{D}}_{\\text{i}-1,\\text{j}-1})}\\right)\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(=\\text{log}\\left(prio{r}_{i,j}\\right)+\\text{log}\\left({10}^{a}+{10}^{b}\\right)=\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(=\\text{log}\\left(prio{r}_{i,j}\\right)+\\text{log}\\left({10}^{a}\\left(1+{10}^{b-a}\\right)\\right)=\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(=\\text{log}\\left(prio{r}_{i,j}\\right)+\\text{log}\\left({10}^{a}\\right)+\\text{log}\\left(1+{10}^{b-a}\\right)=\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(=\\text{log}\\left(prio{r}_{i,j}\\right)+a+\\text{l}\\text{o}\\text{g}(1+{10}^{b-a})\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e\u003c/p\u003e\n\u003cdiv id=\"Eque\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Eque\" name=\"EquationSource\"\u003e$$\\text{W}\\text{h}\\text{e}\\text{r}\\text{e} a=\\text{log}\\left(1-\\left({\\iota }+{\\delta }\\right)\\right)+\\text{log}\\left({ M}_{i-1,j-1}\\right) \\text{a}\\text{n}\\text{d} b=\\text{log}\\left(1-{\\epsilon }\\right)+\\text{l}\\text{o}\\text{g}{(\\text{I}}_{\\text{i}-1,\\text{j}-1} + {\\text{D}}_{\\text{i}-1,\\text{j}-1}).$$\u003c/div\u003e\n\u003c/div\u003e\n\u003cp\u003eSeveral methods for PiM execution of log and other transcendental functions have been proposed \u0026lrm;[57]. To implement log-domain additions, we use a look-up table (LUT) to quickly calculate the value of log(1\u0026thinsp;+\u0026thinsp;10\u003csup\u003eX\u003c/sup\u003e). This allows us to compute the value of the logarithmic\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(sum \\left(log\\right(x+y\\left)\\right)\\)\u003c/span\u003e\u003c/span\u003e according to the above equations. Such LUTs are implemented in PiM core memory prior to the operation. The log values of constants and known-in-advance initial values are precalculated by the host CPU and stored in the PiM core memory prior to the operation.\u003c/p\u003e\n\u003cp\u003eOur second contribution is replacing the floating-point calculations by fixed point ones. This reduces the number of slow multiplications and eliminates the need for emulated floating-point operations, considerably improving the performance and efficiency of the algorithm implementation on PiM cores. We modify the Pair-HMM recursion step as shown in \u003cem\u003eEq.\u0026nbsp;(4)\u003c/em\u003e below, where \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({+}_{F}\\)\u003c/span\u003e\u003c/span\u003e represents a fixed-point regular addition and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({+}_{L}\\)\u003c/span\u003e\u003c/span\u003e represents a logarithmic addition (i.e., log(a\u0026thinsp;+\u0026thinsp;b)) performed using the LUT approach.\u003c/p\u003e\n\u003cdiv id=\"Equf\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equf\" name=\"EquationSource\"\u003e$$\\left\\{\\begin{array}{c}{M}_{i,j}={log}\\left(prio{r}_{i,j}\\right){+}_{\\varvec{F}}[\\left(log\\left(1-\\left(\\iota +\\delta \\right)\\right){+}_{\\varvec{F}}{ M}_{i-1,j-1}\\right)\\\\ {+}_{L}\\left({log}\\left(1-\\epsilon \\right){+}_{\\varvec{F}} {(I}_{i-1,j-1} {+}_{L} {D}_{i-1,j-1}\\right)\\left)\\right]\\\\ {I}_{i,j}= {{log}\\left(\\epsilon \\right){+}_{\\varvec{F} }I}_{i-1,j} {+}_{L} log({\\iota ){+}_{\\varvec{F}} M}_{i-1,j}\\\\ {D}_{i,j}={{log}\\left(\\epsilon \\right){+}_{\\varvec{F}} D}_{i,j-1} {+}_{L} {{log}\\left(\\delta \\right){+}_{\\varvec{F}} M}_{i,j-1}\\end{array}\\left(4\\right)\\right.$$\u003c/div\u003e\n\u003c/div\u003e\n\u003cp\u003e\u003cem\u003eB. Accuracy analysis\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eReplacing the floating-point operations by fixed-point ones affects the computation precision and therefore, impacts the accuracy of the algorithm (i.e., the ability of the algorithm to correctly identify all the variants). We investigate this accuracy impact by comparing our modified GATK implementation with the original GATK HaplotypeCaller algorithm.\u003c/p\u003e\n\u003cp\u003eWe use two samples from the International Genome Sample Resource (IGSR) [46], \u003cem\u003eNA19685\u003c/em\u003e and \u003cem\u003eNA12878\u003c/em\u003e, for each one WES (Whole Exome Sequencing, i.e., protein-coding regions of genes in the genome [47]) with 30x coverage) and WGS (Whole Genome Sequencing) with low coverage. These genomes are aligned to the Hg38 reference and the GATK best practices pipeline is applied on each one.\u003c/p\u003e\n\u003cp\u003eWe optimize the fixed-point format of our calculations using the Matlab Fixed-Point Designer tool [48]. This tool suggests the optimal fixed-point configuration for the variables involved, ensuring accurate and reliable results. It enables us to select the appropriate wordlength for each variable, minimizing the accuracy loss in the Pair-HMM execution.\u003c/p\u003e\n\u003cp\u003eTo compare the performances of the variant callers we first define the following:\u003c/p\u003e\n\u003cp\u003eTrue positive (TP): variants found by both the original GATK and by fixed-point precision GATK.\u003c/p\u003e\n\u003cp\u003e\u0026bull; False Negative (FN): variants found by original GATK but not found by fixed-point precision GATK.\u003c/p\u003e\n\u003cp\u003e\u0026bull; False Positive (FP): non-existing variants found by fixed-point precision GATK.\u003c/p\u003e\n\u003cp\u003eThe metrics of precision (P) and sensitivity (S) are defined as follow, the reference being the original GATK program:\u003c/p\u003e\n\u003cdiv id=\"Equg\" class=\"Equation\"\u003e\n \u003cdiv class=\"mathdisplay\" id=\"FileID_Equg\" name=\"EquationSource\"\u003e$$P=\\frac{TP}{TP+FP}, S=\\frac{TP}{TP+FN} \\left(5\\right)$$\u003c/div\u003e\n\u003c/div\u003e\n\u003cp\u003eWe use the \u003cem\u003evcfeval\u003c/em\u003e tool [50] to compute the different metrics testing four samples (Whole Exome Sequencing (WES) and Whole Genome Sequencing (WGS)) of two different subjects (NA19685 and NA12878). Our analysis shows that fixed-point utilization affects the VCF output, as can be seen in Table \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e below. However, the extent to which this impact is acceptable or not may depend on application requirements and should be evaluated accordingly.\u003c/p\u003e\u0026nbsp;\u003ctable id=\"Tab2\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv class=\"CaptionNumber\"\u003eTable 2\u003c/div\u003e\n \u003cdiv class=\"CaptionContent\"\u003e\n \u003cp\u003eFixed-point impact on sensitivity and precision\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eSample\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eNA19685\u003c/em\u003e\u003c/p\u003e\n \u003cp\u003eWES\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eNA19685\u003c/em\u003e\u003c/p\u003e\n \u003cp\u003eWGS\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eNA19685\u003c/em\u003e\u003c/p\u003e\n \u003cp\u003eCHR18\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eNA12878\u003c/em\u003e\u003c/p\u003e\n \u003cp\u003eWES\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eCoverage\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eHigh\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eLow\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eHigh\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eHigh\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eBits\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e32\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e32\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e12\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e12\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eDecimal bits\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e25\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e25\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e7\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e8\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eP\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.9998\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.9996\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.9717\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.9993\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eS\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.9998\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.9996\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.9996\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.9996\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n\u003c/table\u003e\n\u003cp\u003e\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eC. PiM core initialization\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eWe pre-calculate the log of the initial values and constants in the host CPU since the log function does not exist in the PiM core runtime library. These log values are then transferred in a fixed-point format directly to the PiM core. For instance, the initialization value for \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({D}_{i,0}\\)\u003c/span\u003e\u003c/span\u003e is sent to the PiM core as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(fixed\\_point\\left(log\\right(\\frac{1}{m}\\left)\\right)\\)\u003c/span\u003e\u003c/span\u003e, and the transition values are sent as the fixed point of their log values.\u003c/p\u003e"},{"header":"Implementation on PiM cores","content":"\u003cp\u003eA. \u003cem\u003eCPU-PiM data transfers\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe data of all the regions assigned to a particular PiM core is transferred to the PiM core\u0026rsquo;s MRAM. It includes the reads of the regions and their auxiliary information (read bases, length, and log-domain transitions probabilities), as well as the haplotypes of the regions and their information (read sequence, length and log-domain initial value (\u003cem\u003efixed_point(log(1/m))\u003c/em\u003e as detailed in Section \u003cem\u003eV.C.\u003c/em\u003e)).\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eB. PiM kernel\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eEach PiM thread is assigned a read; then each PiM thread calculates the likelihood of such read against each haplotype in the region. Whenever a read from a new region is assigned to the PiM core, the haplotypes for that region are transferred to the WRAM and remain until the processing for the region is completed. This approach allows us to effectively exploit data reuse, as the haplotypes can be accessed multiple times during the region processing.\u003c/p\u003e\n\u003cp\u003eEach running PiM thread fetches (function \u003cem\u003ereserve_read\u003c/em\u003e in \u003cem\u003eAlgorithm 1\u003c/em\u003e) a read from a pool of reads (either at the start of the program or when it completes handling its current read) and then proceeds to transfer the relevant data associated with that read to the WRAM to compute the likelihood result of the associated read against each haplotype of the region. To manage the read allocation to PiM threads, we maintain a global counter that is shared among all PiM threads, and use a mutex to prevent simultaneous allocation of reads to multiple PiM threads, as explained above (part IV.C).\u003c/p\u003e\n\u003cp\u003eA simplified pseudo-code run by each PiM thread is presented in Algorithm 1.\u003c/p\u003e\n\u003cdiv class=\"gridtable\"\u003e\n \u003cdiv align=\"left\" class=\"colspec\"\u003e\u003cbr\u003e\u003c/div\u003e\u0026nbsp;\u003ctable id=\"Taba\" border=\"1\"\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cem\u003eAlgorithm 1. Pseudo code for PiM core\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003ctfoot\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e1: \u003cstrong\u003efunction\u003c/strong\u003e allocate_read_for_thread\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e2: mram_read(read length)\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e3: mram_read(read sequence)\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e4: mram_read(read transitions)\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e5: \u003cstrong\u003efunction\u003c/strong\u003e reserve_read\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e6: mutex_lock()\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e7: result\u0026thinsp;=\u0026thinsp;free_read_index++\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e8: \u003cstrong\u003eif\u003c/strong\u003e region is new \u003cstrong\u003ethen\u003c/strong\u003e\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e9: mram_read(region haplotypes)\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e10: mutex_unlock\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e11: \u003cstrong\u003ereturn\u003c/strong\u003e result\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e12: \u003cstrong\u003efunction\u003c/strong\u003e main\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e13: \u003cstrong\u003ewhile\u003c/strong\u003e reserve_read\u0026thinsp;\u0026lt;\u0026thinsp;reads number \u003cstrong\u003edo\u003c/strong\u003e\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e14: allocate_read_for_thread\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e15: \u003cstrong\u003efor\u003c/strong\u003e each haplotype in region \u003cstrong\u003edo\u003c/strong\u003e\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e16: initialize matrices\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e17: \u003cstrong\u003efor\u003c/strong\u003e each base in haplotype \u003cstrong\u003edo\u003c/strong\u003e\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e18: \u003cstrong\u003efor\u003c/strong\u003e each base in read \u003cstrong\u003edo\u003c/strong\u003e\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e19: update matrices\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e20: sum likelihood\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e21: mram_write(result likelihood)\u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"1\"\u003e22: \u003cstrong\u003ereturn\u003c/strong\u003e\u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tfoot\u003e\n \u003c/table\u003e\n\u003c/div\u003e\n\u003ch2\u003eC. Efficient handling of on-chip memory\u003c/h2\u003e\n\u003cp\u003eTo optimize the memory usage, the haplotypes present in the WRAM (that have been transferred from the MRAM) are stored in a circular buffer. The utilization of the circular buffer allows for a continuous reuse of memory by replacing the oldest data with new data, preventing wastage, and making efficient use of the available memory. Next, we transfer the read information from the MRAM to the WRAM, including the length, sequence, and transitions probabilities (respectively lines 2,3,4 in \u003cem\u003eAlgorithm 1\u003c/em\u003e).\u003c/p\u003e\n\u003cp\u003eIn order to fill the matrices efficiently and optimize memory usage, we only keep two lines of each matrix (for each PiM thread). This is because, as presented in Eq.\u0026nbsp;(\u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e), the value in each matrix cell depends only on its neighbors as illustrated in Fig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003e. The intermediate likelihood result for each read-haplotype pair is stored in the WRAM. After completing the pairwise alignment of the read and the haplotype, the resulting likelihood is written to the MRAM memory in a result matrix (line 21 in \u003cem\u003eAlgorithm 1\u003c/em\u003e). In order to maximize the bandwidth of the WRAM-MRAM transfer, we coalesce writes from several PiM threads.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eD. Load balancing over PiM cores\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo optimize the performance of the parallel PiM cores, it is important to ensure that the workload is evenly distributed among them. One alternative approach to the method we employ in this work involves dividing genome regions using a SNP database, as reported in [38], assuming that the number of variants in the sample is comparable to the number of variants in the database.\u003c/p\u003e\n\u003cp\u003eOur approach is to allocate several regions to each PiM core based on their complexity, which is determined by the total length of reads and haplotypes in the region. The load balancing is done by the host CPU that evaluates the complexity of each region and potentially divides regions into smaller fragments to meet a specified complexity target.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eE. Asynchronous implementation\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eE. As presented in Section V, the pair-HMM algorithm requires redesign to run efficiently on PiM cores, including pre-processing and post processing steps running on the host CPU, such as splitting the regions. PiM cores can be launched asynchronously, i.e., we can launch them by rank of 64 PiM cores and each rank can be launched independantly of the others. We leverage this to minimize the load imbalance among PiM cores ranks. To efficiently utilize available ressources and optimize PiM cores work we construct a pipeline as illustrated in \u003cem\u003eFig.\u0026nbsp;6\u003c/em\u003e below. To link between the different stages of this pipeline we use buffer queues (one input queue and one output queue). This pipeline is executed on the host CPU and contains three stages:\u003c/p\u003e\n\u003cul\u003e\n \u003cli\u003e\u003cem\u003ePre-processing stage\u003c/em\u003e handled by one single thread that writes to the input queue.\u003c/li\u003e\n \u003cli\u003e\u003cem\u003ePiM cores populating and launching stage\u003c/em\u003e. In this stage, a thread is created for each PiM core rank and reads from the input queue to transfer data to PiM cores and writes to the output queue the final result produced by the PiM cores.\u003c/li\u003e\n \u003cli\u003e\u0026bull; \u003cem\u003ePost processing stage\u003c/em\u003e handled by a single thread that reads from the output queue and writes the result to an output file.\u003c/li\u003e\n\u003c/ul\u003e\n\n\u003cp\u003eBy constructing the pipeline in this manner, we are able to achieve latency hiding during the pre-processing, data transfer, and post-processing steps. This is because the time spent on PiM core computation is much greater than the time spent on these tasks. Thus, when computing the running time on future architectures (as discussed in VII) we disregard the impact of non PiM core time, as explained in [53].\u003c/p\u003e\n\u003cp\u003eThe profiler tool provided by UPMEM [59] illustrates the pipeline parallelization as shown in Fig. \u003cspan class=\"InternalRef\"\u003e7\u003c/span\u003e, with two levels of parallelism: the first is inter-rank parallelism, meaning that several ranks run in parallel (green segments represents PiM core time), and the second is thread parallelism on the host side, meaning that pre-processing task (blue segments), data transfer task (red segments), and post-processing tasks (orange segments) are performed concurrently.\u003c/p\u003e\n"},{"header":"Evaluation","content":"\u003cp\u003eA. \u003cem\u003eEvaluation Methodology\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe algorithm is implemented on servers provided by UPMEM. The architecture is as follows:\u003c/p\u003e\n\u003cul\u003e\n \u003cli\u003eHost: Intel Xeon Silver 4215 CPU, 2.5GHz\u003c/li\u003e\n \u003cli\u003e2560 PiM cores divided into 40 ranks.\u003c/li\u003e\n \u003cli\u003eEach PiM core operates at a frequency of 400MHz.\u003c/li\u003e\n\u003c/ul\u003e\n\u003cp\u003e\u0026nbsp;The benchmarks used in this paper are obtained using a publicly available Whole Genome Sequencing dataset of \u003cem\u003eNA12878\u003c/em\u003e from the International Genome Sample Ressource (IGSR) [46]. We use real datasets in order to ensure transparency and accurately reflect the performance of the Pair-HMM algorithm in real-world scenarios. We present the results for chromosomes 1, 2 and 18 as well as for WGS dataset.\u003c/p\u003e\n\u003cp\u003eWe compare the processing time of our solution with different CPU implementations, including the original Java implementation, the AVX-accelerated version, and the OpenMP multithreaded version (with 4 threads and 8 threads). Additionally, we present the time results of the open-source FPGA implementation from [34], which includes two versions, one with 24 \u0026ldquo;workers\u0026rdquo; and the second one with 96 \u0026ldquo;workers\u0026rdquo;, where a worker is an acceleration unit (respectively represented as 24wk and 96wk). This Pair-HMM implementation is run on a F1 instance of AWS [54].\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eB. Results\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe GATK execution time figures are presented in Table \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003e below along with the speedup figures, using the original Java implementation [22] of the algorithm as the baseline.\u003c/p\u003e\n\u003cp\u003eThe results indicate that the existing UPMEM architecture (2560 PiM cores operating at 400MHz, indicated by bold numbers) can provide speedups similar to that achieved by CPU or FPGA. For example, the UPMEM platform is able to process chromosome 1 in 117 seconds, which is slightly longer than 100 seconds it took the CPU with OpenMP 4 threads, and slightly faster than the 124 seconds needed for the FPGA implementation. The performance achieved by this existing UPMEM architecture represents up to 10.7x speedup compared to the original Java implementation of the algorithm.\u003c/p\u003e\n\u003cp\u003eUPMEM has laid plans for continued improvement and innovation of its product offerings, and these results serve as a prediction of what can be expected from these future architectures. We present in \u003cem\u003eFigure 8\u003c/em\u003e below the performance roadmap [60] of the UPMEM architecture as inspired by [49]. Each percent figure shows the relative improvement of DPU\u0026rsquo;s frequency in the related time period.\u003c/p\u003e\n\u003cp\u003eA scaled-up UPMEM platform with 2560 PiM cores running at 550MHz (and above) surpasses both the CPU with OpenMP 4 threads and the FPGA based platforms. This can be seen in the speedup achieved for chromosomes 1 and 2, which are 14.8x and 12.7x respectively compared to 12.5x and 10.6x for CPU (OpenMP 4t) and to 10x and 11.6x for the FPGA speedup respectively. This is a promising result as UPMEM plans on releasing a version of their architecture with PiM cores working at 550MHz in Q2 2023.\u003c/p\u003e\n\u003cp\u003e\u003c/p\u003e\n\u003ctable id=\"Tab3\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv class=\"CaptionNumber\"\u003eTable 3\u003c/div\u003e\n \u003cdiv class=\"CaptionContent\"\u003e\n \u003cp\u003eRunning time in seconds and acceleration factor of the Pair-HMM algorithm on different architectures\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\u0026nbsp;\u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eCHR1\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eCHR2\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eCHR18\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eWGS\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eJava\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e(reference)\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e1254.8\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e1x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e638.7\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e1x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e162\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e1x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e15159.6\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e1x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eAVX\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e317.1\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e4x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e163\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e3.9x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e40.2\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e4x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e3759\u003c/p\u003e\n \u003cp\u003e4x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eOpenMP\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e4 threads\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e100.5\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e12.5x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e60\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e10.6x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e15.8\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e10.2x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e1189.4\u003c/p\u003e\n \u003cp\u003e12.7x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eOpenMP\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e8 threads\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e75.3\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e16.7x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e46.9\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e13.6x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e12.9\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e12.6x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e886.7\u003c/p\u003e\n \u003cp\u003e17x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eOpenMP 4t\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003eFPGA 24wk\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e140.8\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e8.9x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e57.6\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e11x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e12.3\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e13.2x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e1298\u003c/p\u003e\n \u003cp\u003e11.7x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eOpenMP 4t\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003eFPGA 96wk\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e124.4\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e10x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e55.2\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e11.6x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e12.2\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e13.3x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e1200\u003c/p\u003e\n \u003cp\u003e12.6x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eUPMEM\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e2554 PiM cores\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e350MHz\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e134\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e9.4x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e79\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e8x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e23\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e7x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e1637\u003c/p\u003e\n \u003cp\u003e9.2x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eUPMEM\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e2560 PiM cores\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e400MHz\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003e117\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e10.7x\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003e69\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e9.3x\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003e20.1\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e8x\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003e1429\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e10.6x\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eUPMEM\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e2560 PiM cores\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e550MHz\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e85\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e14.8x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e50.2\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e12.7x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e14.6\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e11.1x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e1039\u003c/p\u003e\n \u003cp\u003e14.6x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eUPMEM\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e3584 PiM cores\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e550MHz\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e60.8\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e20.6x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e35.8\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e17.8x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e10.4\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e15.6x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e742.3\u003c/p\u003e\n \u003cp\u003e20.4x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eUPMEM\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e3584 PiM cores\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e600MHz\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e55.7\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e22.5x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e32.8\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e19.5x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e9.6\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e16.9x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e680.5\u003c/p\u003e\n \u003cp\u003e22.3x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e\u003cstrong\u003eUPMEM\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e5120 PiM cores\u003c/strong\u003e\u003c/p\u003e\n \u003cp\u003e\u003cstrong\u003e600MHz\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e39\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e32.2x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e23\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e27.8x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e6.7\u003c/p\u003e\n \u003cp\u003e\u003cem\u003e24.2x\u003c/em\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e476.3\u003c/p\u003e\n \u003cp\u003e31.8x\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n\u003c/table\u003e\n\u003cp\u003e\u003c/p\u003e\n\u003cp\u003eFigure \u003cspan class=\"InternalRef\"\u003e9\u003c/span\u003e presents (a) the comparison between CPUs execution time (OpenMP four threads and OpenMP eight threads) and PiM cores execution time and (b) the comparison between FPGAs execution time and PiM cores execution time over different architectures for the chromosome 1:\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eC. Power consumption and energy\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe power consumption of the PiM core depends on its version and clock speed. The actual consumption of a PiM core is not yet available on the current chip version but the PiM core\u0026rsquo;s TDP (Thermal Design Power, i.e., the measure of the maximum amount of power that the PiM core consumes) is provided by the manufacturer. When running at 350MHz, the TDP of the new version of the PiM core is 112mW, whereas at 400MHz it consumes 220mW [60]. PiM cores that run at 550MHz are designed to consume less energy, and their TDP is expected to be around 200mW.\u003c/p\u003e\n\u003cp\u003eWe compute the overall TDP of the PiM system used for testing, (according to the values provided by the manufacturer [60]) composed of 2560 PiM cores operating at 350MHz. We recall that the UPMEM chip is composed of eight PiM cores and thus this system contains 320 chips.\u003c/p\u003e\n\u003cp\u003e\u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(TDP\\left( PiM core \\right)=112mW\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(TDP\\left( One UPMEM chip\\right)=8\\text{*}TDP\\left( PiM core \\right)=0.9W\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e \u003cspan class=\"InlineEquation\"\u003e\u0026nbsp;\u003cspan class=\"mathinline\"\u003e\\(TDP\\left( System \\right)=320\\text{*}0.9=288W\\)\u003c/span\u003e\u0026nbsp;\u003c/span\u003e\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eD. Summary analysis\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo summarize this analysis, we compare the speedup achieved by each studied architecture as a function of the number of acceleration units (i.e., the level of parallelism), as presented in Fig. \u003cspan class=\"InternalRef\"\u003e10\u003c/span\u003e below. The X axis of the graph is the number of acceleration units (as specifically defined in each architecture) and the Y axis of the graph represents the speedup factor for the corresponding accelerator. The PiM UPMEM x-axis represents a projection of performances on future architectures, having as a basis the current existing PiM system (2560 DPUs working at 400MHz).\u003c/p\u003e\n\u003cp\u003eThe speedups of OpenMP and FPGA accelerators increase very mildly as the number of acceleration units grows. In contrast, the UPMEM accelerator\u0026rsquo;s performance scales practically linearly. Indeed, the OpenMP multithreading implementation is limited and suffers from high Cycles Per Instruction (CPI) leading to a low performance increase. Similarly, FPGAs suffer from limited on-chip memory and their high hardware resource consumption makes it impractical to add more workers units.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eE. Discussion\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eUPMEM PiM platform based Pair-HMM algorithm outperforms state-of-the-art solutions in terms of performance and power consumption. Based on UPMEM performance roadmap, such performance and energy efficiency gaps are expected to grow in the near future. This is because our implementation operates UPMEM DPU at very modest operating frequency which leaves an ample opportunity for further improvements. At the same time, we cannot expect the operating frequency of high-performance CPU and FPGA which are already fully optimized, to improve at similar pace.\u003c/p\u003e\n\u003cp\u003eImplementing the GATK on UPMEM platforms carries an additional unique advantage. Genome analysis pipeline, which GATK is a major part of, comprises indexing, read mapping, read alignment and in some cases, genome assembly (reference guided or de-novo). Most of these operations are data-intensive and therefore strongly benefit from PiM implementation [10,11,12,13,14,15]. If indexing, read mapping and other data-intensive operations were to be on a PiM platform, while maintaining the conventional CPU-centric (CPU or FPGA based) GATK implementation, it would lead to several significant drawbacks. First, the overall performance would be adversely affected due to the necessity to transfer extremely large amount of sequenced data from PiM platform to a CPU or a FPGA (genomic read file can reach more than 100GB for a single human genome). Second, the total cost of ownership would be significantly higher. Indeed, a heterogeneous solution with two different processing clusters (PiM and FPGA for example) induces higher capital expenditure (CAPEX) due to the cost of the different hardware, as well as higher operational expenditure (OPEX) due to the need for maintenance and increased energy consumption. The ability to implement the entire genome analysis pipeline, from the read mapping to the GATK on the same PiM platform resolves both these inefficiencies. As demonstrated in this paper, to fully realize the benefits of PiM architecture in genome analysis, it is necessary to develop algorithms that are optimized for these architectures, including the parallelization approach, genome region assignment, computing in log domain and so on.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eThis work implements and evaluates Pair-HMM, the main part of the GATK variant calling algorithm and the most computationally intensive part of genome analysis pipeline, on a commercial PiM platform developed by UPMEM. We redesign the original algorithm to optimize it for the UPMEM PiM platform. The suggested modifications include moving the computation to log domain and replacing floating point calculations by fixed point one with very insignificant loss of the variant calling accuracy. To optimize the UPMEM platform utilization and performance, we develop a scheme for optimal allocation of genome fragments to individual PiM cores. This work further demonstrates the potential of PiM architectures as a future solution for end-to-end genome analysis. Despite some parts of the process not being memory-bound, the results show that even complex algorithms such as Pair-HMM can be effectively implemented on PiM architecture and provide speedup over state-of-the-art FPGA or CPU solution.\u003c/p\u003e"},{"header":"Declarations","content":"\u003ch2\u003e\u003cstrong\u003eEthics approval and consent to participate\u0026nbsp;\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003ch2\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003ch2\u003e\u003cstrong\u003eAvailability of data and materials\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eThe code of this research is available at https://github.com/naomieab/UPMEM_HAPLOTYPECALLER\u003c/p\u003e\n\u003ch2\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eNone to declare.\u003c/p\u003e\n\u003ch2\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eThis work was supported by European Union\u0026rsquo;s Horizon Europe programme for research and innovation (grant Number 101047160).\u003c/p\u003e\n\u003ch2\u003e\u003cstrong\u003eAuthors\u0026rsquo; contributions\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eThe code of this research has been written by Naomie Abecassis and Aphelie Moisson-Franckhauser. The research has been supervised and led by Dr. Leonid Yavits and Prof. Ran Ginosar and has also been led in partnership with Dr. Juan G\u0026oacute;mez-Luna and Prof. Onur Mutlu. All authors reviewed the manuscript.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003ePettersson E, Lundeberg J, Ahmadian A. Generations of sequencing technologies. Genomics. 2009 Feb 1;93(2):105-11.\u003c/li\u003e\n\u003cli\u003eReis-Filho JS. Next-generation sequencing. Breast cancer research. 2009 Dec;11(3):1-7.\u003c/li\u003e\n\u003cli\u003eMardis ER. DNA sequencing technologies: 2006\u0026ndash;2016. Nature protocols. 2017 Feb;12(2):213-8.\u003c/li\u003e\n\u003cli\u003eRoodi M, Moshovos A. Gene sequencing: where time goes. In2018 IEEE International Symposium on Workload Characterization (IISWC) 2018 Sep 30 (pp. 84-85). IEEE.\u003c/li\u003e\n\u003cli\u003eStephens ZD, Lee SY, Faghri F, Campbell RH, Zhai C, Efron MJ, Iyer R, Schatz MC, Sinha S, Robinson GE. Big data: astronomical or genomical?. PLoS biology. 2015 Jul 7;13(7):e1002195.\u003c/li\u003e\n\u003cli\u003eMutlu O, Ghose S, G\u0026oacute;mez-Luna J, Ausavarungnirun R. Processing data where it makes sense: Enabling in-memory computation. Microprocessors and Microsystems. 2019 Jun 1;67:28-41.\u003c/li\u003e\n\u003cli\u003eGhose S, Boroumand A, Kim JS, G\u0026oacute;mez-Luna J, Mutlu O. Processing-in-memory: A workload-driven perspective. IBM Journal of Research and Development. 2019 Aug 8;63(6):3-1.\u003c/li\u003e\n\u003cli\u003eKim VS, Fallin C, Lee D, Ausavarungnirun R, Luo GP, Mutlu O, Gibbons PB, Kozuch MA, Mowry TC. RowClone: Fast and efficient In-DRAM copy and initialization of bulk data. Microarchitecture.\u003c/li\u003e\n\u003cli\u003eMutlu O. Memory scaling: A systems architecture perspective. In2013 5th IEEE International Memory Workshop 2013 May 26 (pp. 21-25). IEEE.\u003c/li\u003e\n\u003cli\u003eKhalifa M, Ben-Hur R, Ronen R, Leitersdorf O, Yavits L, Kvatinsky S. FiltPIM: In-memory filter for DNA sequencing. In2021 28th IEEE International Conference on Electronics, Circuits, and Systems (ICECS) 2021 Nov 28 (pp. 1-4). IEEE.\u003c/li\u003e\n\u003cli\u003eKaplan R, Yavits L, Ginosasr R. BioSEAL: In-memory biological sequence alignment accelerator for large-scale genomic data. InProceedings of the 13th ACM International Systems and Storage Conference 2020 May 30 (pp. 36-48).\u003c/li\u003e\n\u003cli\u003eAngizi S, Fahmi NA, Zhang W, Fan D. Pim-assembler: A processing-in-memory platform for genome assembly. In2020 57th ACM/IEEE design automation conference (DAC) 2020 Jul 20 (pp. 1-6). IEEE.\u003c/li\u003e\n\u003cli\u003eLavenier D, Cimadomo R, Jodin R. Variant calling parallelization on processor-in-memory architecture. In2020 IEEE International Conference on Bioinformatics and Biomedicine (BIBM) 2020 Dec 16 (pp. 204-207). IEEE.\u003c/li\u003e\n\u003cli\u003eDiab S, Nassereldine A, Alser M, G\u0026oacute;mez Luna J, Mutlu O, El Hajj I. A framework for high-throughput sequence alignment using real processing-in-memory systems. Bioinformatics. 2023 May 1;39(5):btad155.\u003c/li\u003e\n\u003cli\u003eDiab S, Nassereldine A, Alser M, Luna JG, Mutlu O, Hajj IE. High-throughput pairwise alignment with the wavefront algorithm using processing-in-memory. arXiv preprint arXiv:2204.02085. 2022 Apr 5.\u003c/li\u003e\n\u003cli\u003ehttps://gatk.broadinstitute.org/ (2023)\u003c/li\u003e\n\u003cli\u003eShastry BS. SNP alleles in human disease and evolution. Journal of human genetics. 2002 Nov;47(11):561-6.\u003c/li\u003e\n\u003cli\u003eEngle LJ, Simpson CL, Landers JE. Using high-throughput SNP technologies to study cancer. Oncogene. 2006 Mar;25(11):1594-601.\u003c/li\u003e\n\u003cli\u003eMullaney JM, Mills RE, Pittard WS, Devine SE. Small insertions and deletions (INDELs) in human genomes. Human molecular genetics. 2010 Oct 15;19(R2):R131-6.\u003c/li\u003e\n\u003cli\u003eBernig T, Chanock SJ. Challenges of SNP genotyping and genetic variation: its future role in diagnosis and treatment of cancer. Expert review of molecular diagnostics. 2006 May 1;6(3):319-31.\u003c/li\u003e\n\u003cli\u003eZhao S, Agafonov O, Azab A, Stokowy T, Hovig E. Accuracy and efficiency of germline variant calling pipelines for human genome data. Scientific reports. 2020 Nov 19;10(1):20222.\u003c/li\u003e\n\u003cli\u003ehttps://github.com/broadinstitute/gatk (2023)\u003c/li\u003e\n\u003cli\u003eAl-Ars Z, Wang S, Mushtaq H. SparkRA: enabling big data scalability for the GATK RNA-seq Pipeline with Apache Spark. Genes. 2020 Jan 3;11(1):53.\u003c/li\u003e\n\u003cli\u003eBenjamin, D. (2018). Local Assembly in HaplotypeCaller and Mutect. Available online: https://github.com/broadinstitute/gatk/blob/master/docs/local_assembly.pdf\u003c/li\u003e\n\u003cli\u003eBenjamin, D. (2018). Pair HMM probabilistic realignment in HaplotypeCaller and Mutect. Broad Institute. Available online: https://github.com/broadinstitute/gatk/blob/master/docs/pair_hmm.pdf\u003c/li\u003e\n\u003cli\u003eSchapranow MP, Plattner H. HIG\u0026mdash;An in-memory database platform enabling real-time analyses of genome data. In2013 IEEE International Conference on Big Data 2013 Oct 6 (pp. 691-696). IEEE.\u003c/li\u003e\n\u003cli\u003eFirnkorn D, Knaup-Gregori P, Lorenzo Bermejo J, Ganzinger M. Alignment of high-throughput sequencing data inside in-memory databases. Ine-Health\u0026ndash;For Continuity of Care 2014 (pp. 476-480). IOS Press.\u003c/li\u003e\n\u003cli\u003eKaplan R, Yavits L, Ginosar R. RASSA: resistive prealignment accelerator for approximate DNA long read mapping. IEEE Micro. 2018 Dec 28;39(4):44-54.\u003c/li\u003e\n\u003cli\u003eKim JS, Senol Cali D, Xin H, Lee D, Ghose S, Alser M, Hassan H, Ergin O, Alkan C, Mutlu O. GRIM-Filter: Fast seed location filtering in DNA read mapping using processing-in-memory technologies. BMC genomics. 2018 May;19(2):23-40.\u003c/li\u003e\n\u003cli\u003eFoley P, Prabhakaran A, Gururaj K, Naik M, Gopalan S, Shargorodskiy A, Brau E. Accelerate Genomics Research with the Broad-Intel Genomics Stack.\u003c/li\u003e\n\u003cli\u003eRauer C, Finamore N. Accelerating genomics research with opencl and fpgas. Altera, Now Part of Intel, Tech. Rep. 2016 Mar.\u003c/li\u003e\n\u003cli\u003eRen S, Sima VM, Al-Ars Z. FPGA acceleration of the pair-HMMs forward algorithm for DNA sequence analysis. In2015 IEEE International Conference on Bioinformatics and Biomedicine (BIBM) 2015 Nov 9 (pp. 1465-1470). IEEE.\u003c/li\u003e\n\u003cli\u003eSampietro D, Crippa C, Di Tucci L, Del Sozzo E, Santambrogio MD. Fpga-based pairhmm forward algorithm for dna variant calling. In2018 IEEE 29th International Conference on Application-specific Systems, Architectures and Processors (ASAP) 2018 Jul 10 (pp. 1-8). IEEE.\u003c/li\u003e\n\u003cli\u003eWertenbroek R, Thoma Y. Acceleration of the Pair-HMM forward algorithm on FPGA with cloud integration for GATK. In2019 IEEE International Conference on Bioinformatics and Biomedicine (BIBM) 2019 Nov 18 (pp. 534-541). IEEE. Code available at: https://github.com/rick-heig/PHMM-F1\u003c/li\u003e\n\u003cli\u003eRen S, Bertels K, Al-Ars Z. Efficient acceleration of the pair-hmms forward algorithm for gatk haplotypecaller on graphics processing units. Evolutionary Bioinformatics. 2018 Mar;14:1176934318760543.\u003c/li\u003e\n\u003cli\u003eRen S, Ahmed N, Bertels K, Al-Ars Z. GPU accelerated sequence alignment with traceback for GATK HaplotypeCaller. BMC genomics. 2019 Apr;20(2):103-16.\u003c/li\u003e\n\u003cli\u003eWu X, Subramaniyan A, Wang Z, Narayanasamy S, Das R, Blaauw D. A High-Throughput Pruning-Based Pair-Hidden-Markov-Model Hardware Accelerator for Next-Generation DNA Sequencing. IEEE Solid-State Circuits Letters. 2020 Dec 16;4:31-5.\u003c/li\u003e\n\u003cli\u003eRen S, Bertels K, Al-Ars Z. GPU-accelerated GATK haplotypecaller with load-balanced multi-process optimization. In2017 IEEE 17th International Conference on Bioinformatics and Bioengineering (BIBE) 2017 Oct 23 (pp. 497-502). IEEE.\u003c/li\u003e\n\u003cli\u003eUPMEM Processing In-Memory (PIM), technology paper, unpublished results\u003c/li\u003e\n\u003cli\u003eG\u0026oacute;mez-Luna J, El Hajj I, Fernandez I, Giannoula C, Oliveira GF, Mutlu O. Benchmarking a new paradigm: Experimental analysis and characterization of a real processing-in-memory system. IEEE Access. 2022 May 10;10:52565-608.\u003c/li\u003e\n\u003cli\u003eG\u0026oacute;mez-Luna J, El Hajj I, Fernandez I, Giannoula C, Oliveira GF, Mutlu O. Benchmarking memory-centric computing systems: Analysis of real processing-in-memory hardware. In2021 12th International Green and Sustainable Computing Conference (IGSC) 2021 Oct 18 (pp. 1-7). IEEE.\u003c/li\u003e\n\u003cli\u003eUPMEM, Introduction to UPMEM PIM. Processing-in-memory (PIM) on DRAM Accelerator (White Paper). 2018.\u003c/li\u003e\n\u003cli\u003eDevaux F. The true processing in memory accelerator. In2019 IEEE Hot Chips 31 Symposium (HCS) 2019 Aug 1 (pp. 1-24). IEEE Computer Society.\u003c/li\u003e\n\u003cli\u003eUPMEM, UPMEM User Manual. Version 2021.1.0. 2021.\u003c/li\u003e\n\u003cli\u003eLLVM. Compiler-RT, LLVM project. https://github.com/llvm/llvm-project/tree/main/compiler-rt/lib/builtins, 2021.\u003c/li\u003e\n\u003cli\u003ehttps://www.internationalgenome.org/ (2023)\u003c/li\u003e\n\u003cli\u003eRabbani B, Tekin M, Mahdieh N. The promise of whole-exome sequencing in medical genetics. Journal of human genetics. 2014 Jan;59(1):5-15.\u003c/li\u003e\n\u003cli\u003ehttps://www.mathworks.com/products/fixed-point-designer.html (2023)\u003c/li\u003e\n\u003cli\u003eHennessy JL, Patterson DA. A new golden age for computer architecture. Communications of the ACM. 2019 Jan 28;62(2):48-60.\u003c/li\u003e\n\u003cli\u003eCleary JG, Braithwaite R, Gaastra K, Hilbush BS, Inglis S, Irvine SA, Jackson A, Littin R, Rathod M, Ware D, Zook JM. Comparing variant call files for performance benchmarking of next-generation sequencing variant calling pipelines. BioRxiv. 2015 Aug 2:023754.\u003c/li\u003e\n\u003cli\u003eMansouri Ghiasi N, Park J, Mustafa H, Kim J, Olgun A, Gollwitzer A, Senol Cali D, Firtina C, Mao H, Almadhoun Alserr N, Ausavarungnirun R. GenStore: a high-performance in-storage processing system for genome sequence analysis. InProceedings of the 27th ACM International Conference on Architectural Support for Programming Languages and Operating Systems 2022 Feb 28 (pp. 635-654).\u003c/li\u003e\n\u003cli\u003eMolka D, Sch\u0026ouml;ne R, Hackenberg D, Nagel WE. Detecting memory-boundedness with hardware performance counters. InProceedings of the 8th ACM/SPEC on International Conference on Performance Engineering 2017 Apr 17 (pp. 27-38).\u003c/li\u003e\n\u003cli\u003eYavits L, Morad A, Ginosar R. The effect of communication and synchronization on Amdahl\u0026rsquo;s law in multicore systems. Parallel Computing. 2014 Jan 1;40(1):1-6.\u003c/li\u003e\n\u003cli\u003eAmazon EC2 F1 Instances. Available: https://aws.amazon.com/ec2/instance-types/f1/ (2023)\u003c/li\u003e\n\u003cli\u003ehttps://github.com/aayasin/perf-tools (2023)\u003c/li\u003e\n\u003cli\u003eIntel Corporation (December, 2017), Intel\u0026reg; VTune\u0026trade; Amplifier. https://software.intel.com/en-us/intel-vtune-amplifier-xe (2023)\u003c/li\u003e\n\u003cli\u003eG\u0026oacute;mez-Luna J, Guo Y, Oliveira GF, Sadrosadati M, Mutlu O. TransPimLib: A Library for Efficient Transcendental Functions on Processing-in-Memory Systems. arXiv preprint arXiv:2304.01951. 2023 Apr 3.\u003c/li\u003e\n\u003cli\u003ePoplin R, Chang PC, Alexander D, Schwartz S, Colthurst T, Ku A, Newburger D, Dijamco J, Nguyen N, Afshar PT, Gross SS. A universal SNP and small-indel variant caller using deep neural networks. Nature biotechnology. 2018 Nov;36(10):983-7. Available at https://github.com/google/deepvariant\u003c/li\u003e\n\u003cli\u003ehttps://sdk.upmem.com/2023.2.0/261_Application_profiling.html#application-profiling (2023)\u003c/li\u003e\n\u003cli\u003eY. Falevoz. UPMEM. Personal Communication. October 2022.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Processing in Memory, Pair-HMM, Variant calling, DNA sequence analysis, High Performance Computing","lastPublishedDoi":"10.21203/rs.3.rs-3245756/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-3245756/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eVariant calling is a fundamental stage in genome analysis that identifies mutations (variations) in a sequenced genome relative to a known reference genome. Pair-HMM is a key part of the variant calling algorithm and its most compute-intensive part. In recent years, Processing-in-Memory (PiM) solutions, which consist of placing compute capabilities near/inside memory, have been proposed to speed up the genome analysis pipeline. We implement the Pair-HMM algorithm on a commercial PiM platform developed by UPMEM. We modify the Pair-HMM algorithm to make it more suitable for PiM execution with acceptable loss of accuracy. We evaluate our implementation on single chromosomes and whole genome sequencing datasets, demonstrating up to 2x speedup compared to existing CPU accelerations and up to 3x speedup compared to FPGA accelerations.\u003c/p\u003e","manuscriptTitle":"GAPiM: Discovering Genetic Variations on a Real Processing-in-Memory System","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2023-08-18 20:42:23","doi":"10.21203/rs.3.rs-3245756/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"d3edbd6c-7eea-4ad8-bb94-700d2bfe8538","owner":[],"postedDate":"August 18th, 2023","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2023-10-10T08:29:26+00:00","versionOfRecord":[],"versionCreatedAt":"2023-08-18 20:42:23","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-3245756","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-3245756","identity":"rs-3245756","version":["v1"]},"buildId":"ehx78VzkSd0WSzXnipQa-","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}