Approaching Long Genomic Regions and Large Recombination Rates with msParSm as an Alternative to MaCS Carlos Montemuiño1, Antonio Espinosa1, Juan C. Moure1, Gonzalo Vera2, Porfidio Hernández1 and Sebastián Ramos-Onsins2 1 Computer Architecture and Operating Systems Department (CAOS), Universitat Autònoma de Barcelona, Bellaterra, Spain. 2Centre for Research in Agricultural Genomics (CRAG) Consortium CSIC-IRTA-UAB-UB Edifici CRAG, Campus UAB, Bellaterra, Spain. Abstract: The msParSm application is an evolution of msPar, the parallel version of the coalescent simulation program ms, which removes the limitation for simulating long stretches of DNA sequences with large recombination rates, without compromising the accuracy of the standard coalescence. This work introduces msParSm, describes its significant performance improvements over msPar and its shared memory parallelization details, and shows how it can get better, if not similar, execution times than MaCS. Two case studies with different mutation rates were analyzed, one approximating the human average and the other approximating the Drosophila melanogaster average. Source code is available at https://github.com/cmontemuino/msparsm. Keywords: coalescence, recombination, sequential Markov coalescent, HPC, MPI Citation: Montemuiño et al. Approaching Long Genomic Regions and Large Correspondence: cmontemu@acm.org Recombination Rates with msParSm as an Alternative to MaCS. Evolutionary Bioinformatics 2016:12 223–228 doi: 10.4137/EBO.S40268. Copyright: © the authors, publisher and licensee Libertas Academica Limited. This is an open-access article distributed under the terms of the Creative Commons CC-BY-NC TYPE: Software or Database Review 3.0 License. Received: May 31, 2016. ReSubmitted: July 19, 2016. Accepted for  aper subject to independent expert blind peer review. All editorial decisions made P publication: July 21, 2016. by independent academic editor. Upon submission manuscript was subject to anti- plagiarism scanning. Prior to publication all authors have given signed confirmation of Academic editor: Liuyang Wang, Associate Editor agreement to article publication and compliance with all applicable ethical and legal Peer Review: Two peer reviewers contributed to the peer review report. Reviewers’ requirements, including the accuracy of author and contributor information, disclosure of reports totaled 450 words, excluding any confidential comments to the academic editor. competing interests and funding sources, compliance with ethical requirements relating to human and animal study participants, and compliance with any copyright requirements Funding: This work has been supported by projects (number: AGL2013-41834-R and of third parties. This journal is a member of the Committee on Publication Ethics (COPE). TIN2014-53234-C2-1-R) of Spanish Ministerio de Ciencia y Tecnologia. The authors confirm that the funder had no influence over the study design, content of the article, or Published by Libertas Academica. Learn more about this journal. selection of this journal. Competing Interests: Authors disclose no potential conflicts of interest. Introduction be extremely efficient for short sequences (less than few mega The field of population genetics endeavors to ascertain how base pairs), but it becomes computationally demanding for basic evolutionary forces, such as natural selection, recom­ simulating long genome regions with large recombination bination, and mutation, shape the patterns of genetic varia­ rate.13 New approaches have emerged to overcome this issue, tion within and between populations.1,2 Computer simulation such as the sequential Markov coalescent and Markov Chain software has traditionally been used to explore analytically Monte Carlo methods5,8,14–16 and alternative implementations intractable genetic models.3,4 Continuous advances in numeri­ of the exact coalescent, including algorithmic and data struc­ cal simulation and the wide availability of computational ture optimizations, such as scrm and msprime.17,18 resources allow researchers to use numerical simulation to test In this article, we focus on coalescent simulators, which mathematical models in virtual populations, and even to ana­ are shown to be computationally intractable when working lyze genetic data.5,6 As a direct consequence, there is a plethora at the genome scale,8,11,19 and more specifically on msPar, 20 of simulators available, each one tailored to a specific scenario. which is the parallel version of Hudson’s ms coalescent simu­ This forces geneticists to choose a simulator depending on the lator, 21 ie, the most classical and widely used coalescent simu­ research being conducted.3 lator.3,7,22 msPar addresses both the problem of sampling long The challenge is to provide flexible simulation programs genomic regions with large recombination rates and running with the capacity for dealing with complex and realistic analysis requiring a large number of samples.20 We have previ­ demographic/evolutionary models, using efficient algorithms ously shown that msPar execution time is orders of magnitude to deal with genome-wide data sets processing complexity.4,7,8 faster than Hudson’s ms when simulating large samples of The availability of cost-effective next-generation sequencing long genomic sequences and large recombination rates using technologies has democratized whole-genome sequencing, a High Performance Computing (HPC) cluster, 20 but it is making large genomic data sets available to most research­ less efficient than the other approximated standard coalescent ers.2,9,10 Among the two approaches to simulation algorithms, simulators such as MaCS and fastsimcoal.13–15 forward-in-time and coalescent-based, the latter is most Our interest in further developing the standard coales­ widely used because of its efficiency and flexibility.3,11,12 In cent simulators resides in the underlying side effects shown by particular, the standard coalescent approach was shown to approximated methods, like the one implemented by MaCS: Evolutionary Bioinformatics 2016:12 223 Montemuiño et al 1. We can suffer loss of accuracy for certain evolutionary both MPI processes are in the same node, thereby saving one models as when sampling from populations separated by copy operation compared against the previous transport layer reduced gene flow.23 used included in OpenMPI (ie, sm BTL). 2. Ignoring type 2 recombination events (from Marjoram To take the advantage of Vader BTL, we have changed and Wall’s classification 24) may impact on the mean and the master–worker topology used by msPar. Given N as the variance of most recent common ancestor times when total number of available processes for parallel computation, long sequences are simulated.25 we define one single global master process coordinating the N−1 remaining worker processes. In msParSm, we have one Regarding the implementation details, we have consid­ master process per compute node (henceforth referred to as erably improved msPar’s memory management, improved its node master) and one global master process that coordinates master–worker implementation, removed the data structure the M node masters (where M is the number of compute dependencies inherited from Hudson’s ms, and made use of nodes). Each master node coordinates the worker processes the Vader shared memory byte transport layer (BTL) pro­ located in its node, returning a consolidated message with vided by OpenMPI v1.8.26 By applying a message-passing/ all of the generated samples back to the global master when shared-memory parallelization strategy, msParSm is able to computation is performed. It is important to note that there is achieve similar execution times to those of MaCS, simulating going to be one node hosting two master processes: the global long stretches of DNA sequences with large recombination master and the node master. rates. Thus, we significantly reduce the computational burden By using non-blocking MPI collectives, we enabled node of the standard coalescent method. For the purpose of com­ masters to also generate samples while waiting for local worker parison, we have used two case studies, with fixed mutation communications; therefore, the loss of worker processes com­ rates approximating the human and Drosophila melanogaster pared with msPar is compensated as much as possible. In a average values. setup with M compute nodes and P cores per node, in msPar, Accuracy of simulations is maintained to be the same we have a total number of M*P − 1 dedicated worker processes, as the one from Hudson’s ms. It is ensured by the fact that while in msParSm the number is M*(P − 1) − 1. the code in charge of sample generation was maintained For the use case of the scientist using a single node for unchanged from original Hudson’s ms code. the computation, eg, a single fat node or even a workstation, It is also important to note that msParSm, as same as the application is not going to use a global master but a single msPar and Hudson’s ms, only generates samples under a master that is also going to generate replica samples. Wright–Fisher neutral model of genetic variation. Another improvement we made is to allow MPI pro­ cesses located in the same compute node as the global master, Methods to directly output the generated samples avoiding the use of In this section, we describe the parallelization approach MPI point-to-point communications with a master process. taken over msPar, the experimental setup we designed for Summarizing the major improvements, we evidence gathering, and the hardware/software used to run the experiments. • greatly decreased the number of memory allocation calls Parallelization approach. We present a new master– by favoring memory reusing and reallocation as much as worker application design that includes the following: remov­ possible; ing the data structure dependencies inherited from Hudson’s • refactored the master–worker strategy to leverage the use ms, refactoring msPar’s memory management, and improv­ of the Vader BTL; ing the communication patterns used in the msPar’s master– • included the master process in the replica generation worker implementation. process; First, we have taken the advantage of the OpenMPI’s • made worker processes more autonomous to reduce the new transport layer (Vader shared-memory BTL) for trans­ number of MPI communications; ferring data between communication end points, which • removed all global structures inherited from Hudson's ms, was introduced in version 1.8. The Message Passing Inter­ making it possible to explore a fine-grained paralleliza­ face (MPI) transport layer used by msPar is known as sm tion approach (eg, using OpenMP and/or CUDA). BTL (shared-memory BTL), which follows a copy-in/ copy-out pattern: when an MPI process X sends a mes­ Experimental setup and case studies. We compared our sage to a process Y, the message is first copied from the msParSm with Hudson’s ms and MaCS in terms of execution X’s buffer to the shared memory, and then the receiver time, using two case studies with fixed population mutation (ie, the process Y) copies the message from the shared memory rate values (ie, θ = 4Nµ, where N is the current population into a buffer.27 The Vader shared memory BTL provides sup­ size and µ is the mutation rate per site), one of them approxi­ port for XPMEM Linux kernel module and allows to directly mating the estimated human average (case 1) and the other transfer messages between sender and receiver buffers when one approximating the estimated D. melanogaster average 224 Evolutionary Bioinformatics 2016:12 Approaching with msParSm as an alternative to MaCS (case 2; θ = 0.001 and 0.005 values, meaning that 1 and 5 to sequential Hudson’s ms version. Nevertheless, for the case are differences between two random individuals for each of maximal recombination (Fig. 1D), we observed a slight 1000 sites, respectively). performance degradation. This happened because the node The number of generated replicas was set at 300, each one memory was mostly exhausted by the MPI processes (90 GB with 100 chromosomes, a genetic region of 1 × 106 bp, and in use out of 96 GB available), implying a penalty because of population recombination per site (4Nρ, where ρ is the muta­ the page swapping. tion rate between two contiguous sites) ranging from 0.0005 In Figure 2, the same graphic data as before are shown, to 0.08 values. Both evolutionary models exhibit no popula­ but for the case study 2. Compared with the case study 1, we tion genetic structure. observed a worse performance for lower recombination rate The parameters used in each of the simulations performed (Fig. 2A), but there was a similar performance for upper values in this work are given in Supplementary Table 1. (ie, 2000, 4000, and 8000). An explanation for this behavior Hardware and software. We implemented the paral­ is associated with the message payload exchanged between lel application in ANSI C and compiled it using the Linux MPI processes. Average message with sample data is 4.5 times GNU Compiler 4.9.1. Inter-process communication was higher for the case study 2, suggesting that the time spent in implemented using OpenMPI 1.10.1.28 We ran the binaries communication is less relevant as long as the recombination on a cluster using up to eight nodes each with two Intel Xeon rate increases. X5660 six-core processors using hyper-threading technology In Figure 3, the average consumed memory in function running at 2.80 GHz, 12 MB L3 cache, and 96 GB of Dou­ of the recombination rate and cores used for the computation ble data rate synchronous dynamic-random access memory is shown, which is registered in case study 1. In Figure 4, the (DDR-RAM). This configuration provided us 12 physical same exact data are shown, but for the case study 2. From processors per compute node. both figures, we can observe how the node’s main memory is exhausted when the maximal recombination rate is reached Results and Discussion (Fig. 4D), irrespective of the number of involved cores, We evaluated the performance of msParSm by analyzing the explaining why the efficiency is impacted given the current speedup and efficiency as a function of increasing numbers of experimental setup. compute nodes. The speedup is defined as Sp = T1/Tp, where The most important contribution for the geneticist is the p is the number of processors, T1 is the execution time of the reduction in the overall running time of the simulation pro­ sequential application, and Tp is the execution time of the cess. In Table 1, the average time Hudson’s ms, MaCS, msPar, parallel application with p processors. Efficiency is defined as and msParSm spent running both case studies using differ­ Ep = Sp /p. ent problem sizes (ie, recombination rate values) is given. We In Figure 1, the speedup of msParSm compared with observe that how the execution time of msParSm is orders of the sequential Hudson’s ms version is shown, which was magnitude faster than Hudson’s ms, faster by a factor than registered in case study 1. We observed a high parallel effi­ msPar, and significantly better than MaCS in most configura­ ciency of msParSm and a strong increase in speed in relation tions, mainly when mutation rate is higher. A Speed-up (case study 1) B Speed-up (case study 1) 90.0 90.0 80.0 ρ = 1000 80.0 ρ = 2000 65.9 70.0 70.0 Speed-up Speed-up 60.0 60.0 47.3 50.0 50.0 38.9 40.0 31.0 40.0 30.0 18.1 30.0 20.2 20.0 9.7 20.0 9.8 10.0 10.0 0.0 0.0 12 cores 24 cores 48 cores 96 cores 12 cores 24 cores 48 cores 96 cores C Speed-up (case study 1) D Speed-up (case study 1) 90.0 90.0 80.0 ρ = 4000 76.4 80.0 ρ = 8000 74.1 70.0 70.0 Speed-up Speed-up 60.0 60.0 50.0 42.5 50.0 41.5 40.0 40.0 30.0 22.2 30.0 21.7 20.0 10.9 20.0 10.9 10.0 10.0 0.0 0.0 12 cores 24 cores 48 cores 96 cores 12 cores 24 cores 48 cores 96 cores Figure 1. Speedup analysis of case study 1 in function of the number of cores used for computation, grouped by a scaled recombination rate (A to D), registered in case study 1. Subfigures show speedup of the application increases adding more computational cores, slowly decaying when considering recombination rates of 4Nρ = 8000. Evolutionary Bioinformatics 2016:12 225 Montemuiño et al A Speed-up (case study 2) B Speed-up (case study 2) 90.0 90.0 80.0 ρ = 1000 80.0 ρ = 2000 68.8 70.0 70.0 Speed-up Speed-up 60.0 60.0 50.0 42.7 50.0 44.6 40.0 40.0 26.5 30.0 30.0 22.9 20.0 15.3 20.0 8.7 11.3 10.0 10.0 0.0 0.0 12 cores 24 cores 48 cores 96 cores 12 cores 24 cores 48 cores 96 cores C Speed-up (case study 2) D Speed-up (case study 2) 90.0 90.0 80.0 ρ = 4000 74.5 80.0 ρ = 8000 74.2 70.0 Speed-up 70.0 60.0 Speed-up 60.0 50.0 41.8 50.0 41.5 40.0 40.0 30.0 21.7 30.0 21.9 20.0 10.6 20.0 10.9 10.0 10.0 0.0 0.0 12 cores 24 cores 48 cores 96 cores 12 cores 24 cores 48 cores 96 cores Figure 2. Obtained speedup of case study 2 in function of the number of cores used for computation, grouped by used recombination rate (A to D). As in the previous case study 1, this figure shows increasing speedups when adding more computational resources. A Memory (ρ = 1000) B Memory (ρ = 2000) 100 100 Gigabytes Gigabytes 80 80 60 60 40 40 20 20 9.34 8.44 7.39 7.27 3.71 2.88 2.51 2.35 0 0 12 cores 24 cores 48 cores 96 cores 12 cores 24 cores 48 cores 96 cores C Memory (ρ = 4000) D Memory (ρ = 8000) 96.27 92.80 91.88 100 100 89.57 80 80 Gigabytes Gigabytes 60 60 40 28.59 26.37 40 25.77 25.18 20 20 0 0 12 cores 24 cores 48 cores 96 cores 12 cores 24 cores 48 cores 96 cores Figure 3. Average consumed memory for case study 1 in function increasing the number of computational cores and recombination values, registered in case study 1. Each subfigure shows the consumed memory values when running the simulation using 12, 24, 48, and 96 cores, starting with a scaled recombination of 4Nρ = 1000 in (A), and doubling it through the remaining subfigures (C to D) until reaching a scaled recombination of 4Nρ = 8000. A Memory (ρ = 1000) B Memory (ρ = 2000) 100 100 Gigabytes 80 80 Gigabytes 60 60 40 40 20 20 13.60 10.71 8.61 7.51 7.37 4.56 3.35 2.82 0 0 12 cores 24 cores 48 cores 96 cores 12 cores 24 cores 48 cores 96 cores C Memory (ρ = 4000) D Memory (ρ = 8000) 101.29 94.93 100 100 92.13 90.65 Gigabytes 80 80 Gigabytes 60 60 40 33.05 29.14 26.77 40 26.16 20 20 0 0 12 cores 24 cores 48 cores 96 cores 12 cores 24 cores 48 cores 96 cores Figure 4. Average consumed memory for case study 2 in function increasing the number of cores used for computation. Each subfigure (A to D) shows the same information as in Figure 3, but with data registered when running the case study 2. 226 Evolutionary Bioinformatics 2016:12 Approaching with msParSm as an alternative to MaCS Table 1. Comparison of average time cost (in minutes) between ms, MaCS, msPar, and msParSm. The execution times of both msParSm and msPar are grouped in function of the number of cores used for computation. A combination of font styles and background color is used to facilitate the reading: values in italics are related to msPar, while bold style is used for msParSm; shaded table cells indicate cases that the execution time of either msParSm or msPar is worse than MaCS. Rho Ms MaCS msParSm/mspar 96 cores 48 cores 24 cores 12 cores Case study 1 1000 3.81 3.82 0.08 0.44 0.12 0.44 0.21 0.64 0.40 0.59 2000 29.77 7.16 0.45 2.75 0.76 2.75 1.47 3.91 3.05 3.42 4000 322.80 13.76 4.23 19.57 7.60 26.07 14.57 28.48 29.64 30.52 8000 2605.58 26.86 35.15 149.48 62.78 214.36 120.02 223.98 240.12 249.5 Case study 2 1000 3.95 5.80 0.22 1.14 0.21 1.22 0.26 1.61 0.44 1.37 2000 34.85 9.14 0.51 3.33 0.78 3.53 1.52 4.89 3.09 4.23 4000 320.30 15.88 4.30 20.02 7.67 27.61 14.78 29.43 30.12 32.03 8000 2602.00 29.32 35.05 151.35 62.75 222.60 119.03 224.18 239.47 249.27 Although the parallel execution of msParSm is behind Supplementary Material MaCS, the scalability of presented msParSm is generally good Supplementary Table 1. Parameters used in each of the when adding more computational resources with the possible simulations performed. limitation of the available memory of the nodes. References Conclusion 1. Hudson R. Gene genealogies and the coalescent process. In: Futuyama D, In this work, we developed a new application for parallelizing the Antonovics J, eds. Oxford Survey in Evolutionary Biology. Vol 7. Oxford: Oxford University Press; 1991:1–44. Hudson’s ms, a standard coalescent simulator, coined as msParSm, 2. Akey JM, Shriver MD. A grand challenge in evolutionary and population genet­ addressing the sampling of long stretches of DNA sequences with ics: new paradigms for exploring the past and charting the future in the post- large recombination rates, and we showed that this new applica­ genomic era. Front Genet. 2011;2:47. 3. Hoban S, Bertorelle G, Gaggiotti OE. Computer simulations: tools for popula­ tion has comparable performance with MaCS when working with tion and evolutionary genetics. Nat Rev Genet. 2011;13(2):110–22. large recombination rate, without the side effects of MaCS. 4. Arenas M. Simulation of molecular data under diverse evolutionary scenarios. PLoS Comput Biol. 2012;8(5):e1002495. Interestingly, the application achieves high parallel effi­ 5. Sanford J, Baumgardner J, Brewer W, Gibson P, ReMine W. Mendel’s ciency figures (.70%) when working with large recombination accountant: a biologically realistic forward-time population genetics program. Scalable Comput Pract Exp. 2001;8(2):147–65. rate, suggesting that it could be run on more nodes than used 6. Sanford J, Nelson C. The next step in understanding population dynamics: for this work and get better execution times. In addition, our comprehensive numerical simulation. In: Fust MC, ed. Studies in Population Genetics. InTech; 2012. Available at: http://www.intechopen.com/books/studies-in- application can outperform MaCS in most of the cases, even population-genetics/the-next-step-in-understanding-population-dynamics- when running the application in one single node. comprehensive-numerical-simulation. Accessed March 29, 2016. We believe that this new application will facilitate simu­ 7. Carvajal-Rodríguez A. Simulation of genomes: a review. Curr Genomics. 2008;9(3):155–9. lating coalescent processes with long genomic regions and 8. Carvajal-Rodríguez A. Simulation of genes and genomes forward in time. Curr large recombination rates in population genomics and evo­ Genomics. 2010;11(1):58–61. 9. Liu DJ, Leal SM. Replication strategies for rare variant complex trait association lutionary biology, without compromising the accuracy of the studies via next-generation sequencing. Am J Hum Genet. 2010;87(6):790–801. standard coalescence. 10. Pool JE, Hellmann I, Jensen JD, Nielsen R. Population genetic inference from genomic sequence variation. Genome Res. 2010;20(3):291–300. 11. Kim Y, Wiehe T. Simulation of DNA sequence evolution under models of recent Author Contributions directional selection. Brief Bioinform. 2009;10(1):84–96. 12. Peng B, Kimmel M. simuPOP: a forward-time population genetics simulation Conceived and designed the experiments: SRO, CM. environment. Bioinformatics. 2005;21(18):3686–7. Analyzed the data: CM. Wrote the first draft of the manu­ 13. Yang T, Deng H-W, Niu T. Critical assessment of coalescent simulators in script: CM. Contributed to the writing of the manuscript: modeling recombination hotspots in genomic sequences. BMC Bioinformatics. 2014;15:3. PH, SRO. Agreed with the manuscript results and conclu­ 14. Chen GK, Marjoram P, Wall JD. Fast and flexible simulation of DNA sequence sions: SRO, PH, AE, GV, JCM, CM. Jointly developed the data. Genome Res. 2009;19(1):136–42. 15. Excoffier L, Foll M. Fastsimcoal: a continuous-time coalescent simulator of structure and arguments for the article: SRO, PH, CM. Made genomic diversity under arbitrarily complex evolutionary scenarios. Bioinformatics. critical revisions and approved the final version: SRO, PH, 2011;27(9):1332–4. 16. Grünwald NJ, Goss EM. Evolution and population genetics of exotic and AE, GV, JCM. All the authors reviewed and approved the re-emerging pathogens: novel tools and approaches. Annu Rev Phytopathol. final manuscript. 2011;49(1):249–67. Evolutionary Bioinformatics 2016:12 227 Montemuiño et al 17. Staab PR, Zhu S, Metzler D, Lunter G. scrm: efficiently simulating long 23. Eriksson A, Mahjani B, Mehlig B. Sequential Markov coalescent algorithms sequences using the approximated coalescent with recombination. Bioinformatics. for population models with demographic structure. Theor Popul Biol. 2015;31(10):1680–2. 2009;76(2):84–91. 18. Kelleher J, Etheridge AM, McVean G. Efficient coalescent simulation and genea­ 24. Marjoram P, Wall JD. Fast “coalescent” simulation. BMC Genet. 2006;7:16. logical analysis for large sample sizes. PLoS Comput Biol. 2016;12(5):e1004842. 25. Wang Y, Zhou Y, Li L, et al. A new method for modeling coalescent processes 19. Liang L, Zöllner S, Abecasis GR. GENOME: a rapid coalescent-based whole with recombination. BMC Bioinformatics. 2014;15:273. genome simulator. Bioinformatics. 2007;23(12):1565–7. 26. FAQ: Tuning the Run-Time Characteristics of MPI Sm Communications. 20. Montemuiño C, Espinosa A, Moure JC, et al., eds. MsPar: a parallel coalescent Available at: https://www.open-mpi.org/faq/?category=sm. Accessed March 25, simulator. Euro-Par 2013: Parallel Processing Workshops. Lecture Notes in Computer 2016. Science. Berlin, Heidelberg: Springer; 2013:321–30. Available at: http://link. 27. Hjelm NT, Gutierrez SK, Gorentla Venkata M. On the Current State of Open MPI springer.com/chapter/10.1007/978-3-642-54420-0_32. Accessed February 24, on Cray Systems. Oak Ridge National Laboratory (ORNL); Oak Ridge Leader­ 2016. ship Computing Facility (OLCF); 2014. Available at: http://www.osti.gov/ 21. Hudson RR. Generating samples under a Wright-Fisher neutral model of genetic scitech/biblio/1150898-current-state-open-mpi-cray-systems. Accessed April 4, variation. Bioinformatics. 2002;18(2):337–8. 2016. 22. Ewing G, Hermisson J. MSMS: a coalescent simulation program including 28. Open MPI: Open Source High Performance Computing. Available at: https:// recombination, demographic structure and selection at a single locus. Bioinfor- www.open-mpi.org/. Accessed March 29, 2016. matics. 2010;26(16):2064–5. 228 Evolutionary Bioinformatics 2016:12