Moonshine.jl: a Julia package for genome-scale model-based ancestral recombination graph inference

Moonshine.jl:AJulia PACKAGE FOR GENOME-SCALE MODEL-BASED ANCESTRAL RECOMBINATION GRAPH INFERENCE A PREPRINT Patrick Fournier Département de Mathématiques Université du Québec à Montréal Montréal, Québec, H3C 3P8 [email protected] Fabrice Larribe Département de Mathématiques Université du Québec à Montréal Montréal, Québec, H3C 3P8 [email protected] November 27, 2025 ABSTRACT The ancestral recombination graph (ARG) is the model of choice in statistical genetics to model population ancestries. Software capable of simulating ARGs on a genome scale within a reasonable amount of time are now widely available for most practical use cases. While the inverse problem of inferring ancestries from a sample of haplotypes has seen major progress in the last decade, it does not enjoy the same level of advancement as its counterpart. Up until recently, even moderately sized samples could only be handled using heuristics. In recent years, the possibility of model-based inference for datasets closer to "real world" scenarios has become a reality, largely due to the development of threading-based samplers. This article introduces Moonshine, a Julia package that has the ability, among other things, to infer ARGs for samples of thousands of human haplotypes of sizes on the order of hundreds of megabases within a reasonable amount of time. On recent hardware, our package is able to infer an ARG for samples of densely haplotyped (over one marker/kilobase) human chromosomes of sizes up to 10000 in well under a day on data simulated by msprime. Scaling up simulation on a compute cluster is straightforward thanks to a strictly single-threaded implementation. While model-based, it does not resort to threading but rather places restrictions on probability distributions typically used in simulation software in order to enforce sample consistency. In addition to being efficient, a strong emphasis is placed on ease of use and integration into the biostatistical software ecosystem. Keywords Ancestral Recombination Graph ·ARG Inference ·Coalescent Theory ·Algorithms ·Julia 1 Introduction Coalescent theory is the framework of choice when it comes to inference from genetic data ([1]). A brief introduction and some historical notes on the subject can be found in [2]. Although it originates from the seminal papers of Kingman ([3, 4]) and Hudson ([5]), it didn’t really gain widespread use before the advent of computer software capable of simulating the coalescent process with recombination. ms ([5]), a software written in C, was the first of its kind. It aims at generating a sample of genetic sequences evolving under a Wright-Fisher (WF) model via a coalescent-withrecombination (CWR) approximation. The WF model is known for its simplicity. Somewhat surprisingly, simulating a sample of genetic sequences resulting from the evolution of a population obeying it rapidly becomes a challenging task. This is because it is necessary to keep track of every sequence at each of the numerous discrete steps of the simulation procedure from the original population to the desired sample. The coalescent process works in reverse: it starts from the sample and generates coalescence events, which correspond to sets of sequences finding a common ancestor. Consequently, it is only necessary to keep track of haplotypes associated with these events. The remaining sequences are deemed non-ancestral for the sample and can be disregarded. In addition, exponentially distributed arXiv:2511.21124v1 [q-bio.GN] 26 Nov 2025 Moonshine.jl: ARG Inference A PREPRINT coalescence times are substituted to discrete generation counts. Sampling genetic sequences using the coalescent approximation is a straightforward task. What is not trivial, however, is relaxing the assumptions of the WF model to allow recombination events. The result of simulations performed under a model accounting for these events, along with, somewhat confusingly, the statistical model itself, is known as the ancestral recombination graph (ARG). To avoid ambiguity, we reserve “ARG” for the former and refer to the statistical model through “CWR”. While Hudson’s ms is able to simulate recombination events, it is not efficient enough to meet the requirements of large-scale genomic data analysis. To reduce the computational burden associated with simulations, an approximation to the CWR was developed by Niall Cardin and Gil McVean ([6]). Building on the work of [7], their idea is to approximate the non-Markovian recombination process with a Markovian one called the sequential Markov coalescent (SMC). This gave rise to a plethora of so-called sequential simulators such as MaCS ([8]), fastsimcoal ([9]), SC ([10]) and scrm ([11]). While the sequential approach is not inherently more efficient than the backward-in-time one, it is easier to approximate, leading to substantial improvements in the capacity of CWR samplers. Nonetheless, the idea of backward-in-time simulation was not abandoned, as evidenced by MSMS ([12]) and Discoal ([13]). [14] even developed a forward-in-time simulator able to handle nonWF scenarios. Some years later, the classical backward-in-time CWR approach was brought back into the spotlight by msprime ([15]). By reformulating Hudson’s original ms in terms of a data structure they called sparse trees, their Python package (with performance-critical procedures implemented in C) has the capacity of simulating exactly from the coalescent with recombination more efficiently than the sequential approximations available at the time. Version 1.0 ([16]) brings even more features such as the ability to simulate directly from a WF model. For all these reasons, msprime is nowadays the de facto reference for CWR and even ARG simulation. Programs mentioned above simulate genealogies with the goal of generating a sample of genetic sequences. Throughout this paper, we refer to this process as ARG simulation; the ancestry is viewed as a sample point in the probability space associated with the CWR. The focus is generally on inferring parameters such as recombination/mutation rates or effective population size. The distinctive characteristic of this class of samplers is that the distribution is parametrized by values derived from sequences of markers rather than the markers themselves. This is to be contrasted with software designed to solve the inverse problem of generating likely ancestries, or even, in some cases, a single ancestry directly from the markers. This is known as ARG inference since the ancestry itself is usually the main point of interest. This choice of words is, however, somewhat misleading as so-called inferred genealogies are not necessarily central to the analyses they are involved in. They can be instrumental in the estimation of other parameters, just like their sampled counterparts. It is true that their distribution is not that of the CWR. That being said, this is not to say that alternative likelihood functions enabling maximum likelihood estimation do not exist. Indeed, recent work ([17]) proposes formulations with applicability to broad classes of ancestries in mind. While the problem these samplers solve is by nature more computationally demanding, they provide major benefits in that the ancestries they produce are generally more likely with respect to the sample at hand. In particular, their usefulness cannot be overstated for methodologies that aim at improving the estimation of parameters by treating the ancestry of a sample as a latent variable such as [18] or [19]. Those require integrating these genealogies out, a task involving, in practice, the ability to sample in high probability regions. Early attempts at ARG inference include recom ([20]) and Infs ([21]). More recently, the authors of SC also provides a modified version of their algorithm called SC-sample, capable of solving this inverse problem by generating sampleconsistent graphs. A graph is said to be consistent for a sample with respect to a mutation evolution model if it generates the sample under that model. In [10], they choose to use the popular infinite site model (ISM). Other software have been developed with the goal of inferring the ancestry of a sample, such as tsinfer ([22]) and ARGinfer ([23]). Both assume a sample of sequences of binary markers evolving under the ISM. tsinfer uses a heuristic to partially infer haplotypes of unobserved ancestors. Each site of every sequence in the sample is then traced back to an ancestor by maximum likelihood. This yields a point estimate of the ancestral recombination graph. The goal of ARGinfer is slightly different, as it aims at performing the inference probabilistically. It can compute probability intervals for ARG related quantities. All these simulators are designed to sample ARGs according to a probabilistic model. For that reason, we refer to them as model-based simulators. An extensive review of ARG samplers is available in [24]. A related problem is that of parsimony-based inference, which consists of finding ARGs consistent with a sample using the minimum number of recombination events. This problem has been proved NP-hard ([25]). Nonetheless, attempts to solve it exactly and approximately go back to the early work of Hein ([26]), which was improved by [27]. Other algorithms include Margarita ([28]), ARG4WG([29]) and GAMARG([30]). This paper introduces a Julia package called Moonshine implementing inference of sample-consistent ancestral recombination graphs. The approach is sequential, as it is based on iterative modification of the ARG. A sequence of operations are applied to a coalescent tree to ultimately make it consistent with a sample of haplotypes. The main advantage over back-in-time approaches is the possibility it offers users to specify various levels of approximation 2 Moonshine.jl: ARG Inference A PREPRINT when generating ARGs. The available spectrum ranges from exact simulation without any Markov assumption to first-order approximation à la SMC. The end result of the sampling routine is independent of approximation level. The whole graph, as well as meta-data such as vertex latitudes, associated haplotypes, and intervals of ancestrality for edges, is available to users for subsequent analysis. In fact, Moonshine is designed for easy integration into data analysis workflows. In addition, it is possible to use it solely for inferring a set of ancestries consistent with a given sample of genetic sequences. From that point of view, the software it is most similar to is probably SC-sample. That being said, both packages differ in many respects. At first sight, the most obvious difference is that SC-sample is a C++ standalone software while Moonshine is a Julia package. The latter is much more flexible. A run of SC-sample results in a text file containing sequences of Newick-formated trees representing an ARG which may or may not be easily integrated into an existing pipeline. In contrast, thanks to Julia’s just-in-time compilation, Moonshine can be used interactively without sacrificing performance. It is also fully integrated with Julia’s ecosystem of graph theoretical packages. Another noteworthy difference is at the algorithmic level. SC-sample is a modified version of SC which is itself a modified version of MaCS. Although all sequential algorithms are based on the same idea and, consequently, share many similarities, Moonshine implements its own original algorithm. Furthermore, since it is created with statistical inference in mind, Moonshine treats ARGs as random graphs. It is straightforward to evaluate ARG-related functions such as probability densities for the ARGs themselves or other random variables, such as phenotypes, conditional on an ARG. Implementation of custom functionalities is facilitated by a coherent type hierarchy and thorough documentation of abstract types and interfaces, making the extension of the package’s various components as easy as possible. Finally, interoperability with tskit ([15, 31]) streamlines data management. Integration with the popular genomics package is still a work in progress, but creating samples directly from tskit objects is supported. A convenience constructor for generating random samples from a simple genetic model using msprime ([16]) is provided. This is transparent to the end user, thanks to Moonshine being packaged with its own distribution of msprime. Our objective in developing Moonshine’s is not limited to creating a realistic and convenient sampler; performance is a major priority. We present numerical experiments showing its potential for both coalescent tree construction and ARG inference. Trees for large samples (n= 10000) of long simulated haplotypes (250 Mbp) can be constructed in minutes at high resolution (over one marker per kbp) using Hamming distance between sequences. In the same scenarios, complete ancestries can be inferred in hours. Furthermore, our algorithms are completely single-threaded, enabling us to increase sampling throughput by leveraging concurrency efficiently and easily by launching multiple instances in parallel. 2 Tree Construction Similar to other sequential samplers, the first step of our algorithm is to construct an initial coalescent tree. Consistency with the first marker is not assumed; ARG inference can be carried out even starting from a completely random tree. The sole requirement is that it be a valid coalescent tree, i.e., a full binary tree with a coherent set of latitudes for the vertices. The idea behind this functionality is that it might be of interest to compare the performance of a method under a model without recombination versus one that allows for such events. It would make little sense in that context to give a special status to a single marker, disregarding the remainder of the haplotypes. Consequently, we give the user maximum flexibility when building coalescent trees, which may be of interest both in their own right and as a stepping stone for constructing more complex histories. As will be discussed later, Moonshine is packaged with two haplotype metrics designed for tree building. It is straightforward for the user to implement custom metrics. Within our package, data structures representing genealogies are subtypes of the AbstractGenealogy abstract type. The type of coalescent trees is simply Tree. Given a sample (of type Sample) of phased and polarized genetic sequences (of type Sequence), sampling is controlled by two parameters: the global mutation rate µand a metric on haplotypes d. Assuming a diploid population, sequences of length lwith smarkers, an effective population size Neand a constant per locus mutation rate of µ′, the global mutation rate is µ= 2Neµ′l. These parameters are either computed or explicitly passed to Sample’s constructor. As for the metric, since Moonshine is compatible with binary markers exclusively, a sample is an n-tuple H= (h1,··· , hn)of s-vectors of GF(2), the finite field with two elements. Concretely, for each k, we have hk=h1 k···hs k where h• k∈ {0,1}. Wild and derived alleles are represented by 0 and 1 respectively as is standard in the literature. Let ⊕denote addition modulo 2. Examples of useful metrics, also known as distance functions, include dL(h1, h2) = h1 1⊕h1 2 dH(h1, h2) = s X i=1 hi 1⊕hi 2. 3 Moonshine.jl: ARG Inference A PREPRINT dLis the metric under which the distance between two sequences is zero if and only if the state of their first marker is identical, 1 otherwise. dHis the Hamming distance. Such discrete distances have a natural biological interpretation as the number of mutations between (a subinterval of) sequences. Arbitrary distances can be implemented by the user as subtypes of Distance. Detail-oriented readers might have noticed that the term “metric” is used loosely, as the positivity axiom need not hold; the distance between two distinct haplotypes may be 0. This is necessary to allow inference of trees consistent for a single marker. Technically, the correct mathematical construct is that of a pseudometric. 2.1 Sampling Algorithm Coalescence events are sampled as follows: a vertex vais choosen uniformly among the set of live vertices, that is, the vertices that have not coalesced yet. Another vertex vbis sampled conditional on va. The probability pab of sampling vbgiven vais proportional to pab =µdab Γ(dab + 1) (1) where Γis the gamma function, dab =d(ha, hb)with haand hbthe haplotypes associated with vaand vbrespectively. pab is an unnormalized Poisson probability, where Γis used instead of the usual factorial function to allow non-integer distances. Both vertices then coalesce into vcwith associated haplotype hc=ha⊙hbwhere ⊙is the Hadamard product. The shift in latitude is exponentially distributed with rate parameter equal to the number of live vertices. The latitude of vcis computed with respect to that of the previous event. This procedure replaces vaand vbby vcin the set of live vertices. Repeating it n−1times on a sample of nhaplotypes yields a coalescent tree. The user has the possibility of biasing distance computation by specifying a parameter c0∈R+∪ {0,∞} meaning that, in practice, dab =c0d(ha, hb). c0can be used to reproduce the behavior of other samplers. Setting c0=∞and d=dL, we obtain dab =0if h0 a=h0 b ∞if h0 a6=h0 b resulting in pab =1if h0 a=h0 b 0if h0 a6=h0 b . The resulting sampler will agregate haplotypes with identical status at the first marker, resulting in a tree consistent with it. In practice, dealing with a ratio of such extreme quantities poses a numerical challenge. For instance, if one wishes to use Hamming’s distance, computation of the normalizing constant becomes impossible even for relatively small values of nand s. We address these issues using two tricks. First, we compute pab via Stirling’s approximation. On the logarithmic scale, we obtain log pab ≈dab(log µ−log dab + 1) which is already more manageable. Next, it would be ideal to sample without reverting to the linear scale. Moreover, we would greatly benefit from avoiding the computation of the normalizing constant altogether. It turns out that this is exactly what the so-called Gumbel trick [32] is designed to do. The trick transforms sampling from the target distribution into an optimization problem. For a set of candidate vertices b1,...,bmand a sequence of iid standard Gumbel random variables t1, . . . , tm,arg maxk{log pabk+tk}is distributed as a categorical random variable with the probability of category kbeing equal to pabk. This result gives us a more stable way of sampling the second coalescing vertex vb, at the cost of increased computing time for drawing Gumbel random variables and finding the maximum of the sequence. In many cases, this impact should be minimal compared to the overall execution time, for instance when the tree is to be used for ARG inference. Figure 1 provides a rough estimate of time and memory usage for various scenarios. For cases where time considerations warrant a precision tradeoff, users can sample approximately from the target distribution via a scheme we call secretary sampling, inspired by the famed secretary problem ([33]). It revolves around giving the algorithm a chance of terminating before having traversed the complete set of candidate vertices. Its behavior is controlled by a user-determined threshold parameter t0∈[0,1]. As illustrated in fig. 2, the probability of an early termination decreases with t0. Note that although normalization is not required for sampling, the unnormalized probabilities associated with traversed vertices must be summed in order to evaluate the density of the resulting coalescent tree. This has to be done carefully as departure from the logarithmic scale is unavoidable. Details as well as the complete sampling algorithm are presented in algorithm 1. 4 Moonshine.jl: ARG Inference A PREPRINT Algorithm 1: Tree Sampling 1Fill array ‘live’ with the live edges 2while |live|>1do 3Sample vauniformly from live and remove it from live // log pab +Gumbel, log pab, normalization constant and index of sampled vertex respectively 4Set z← ∞,logp ←0,logc ←0and idx ←0 5Set k←0 6while k<|live|do 7Increment k←k+ 1 8Set b←live[k],logp ←log pab 9if logp <∞then 10 Sample g∼Gumbel(0, 1) 11 Set logc ←logp,z←logc +gand idx ←k 12 Break 13 while k<|live|do 14 Increment k←k+ 1 15 Set b←live[k],logp’ ←log pab 16 if logp’ <∞then // Update logc accurately 17 Set pmin ←min{logp’,logc},pmax ←max{logp’,logc} 18 Update logc ←pmax + log(1 + exp(pmin −pmax)) 19 Sample g∼Gumbel(0, 1) 20 Set z’ ←logp’ +g 21 if z’ ≥zthen 22 Set logp ←logp’,z←z’,idx ←k 23 if k>threshold then 24 Break 25 if idx = 0 then // All probabilities were infinite 26 Sample idx uniformly from 1,...,|live| 27 Set logp ←0,logc ←log |live| 28 Set b←live[idx],vc←2n− |live| 29 Remove vbfrom live, coalesce vaand vbinto vcand add vcto live 30 Add logp −logc −log(1 + |live|)to tree’s log density 5 Moonshine.jl: ARG Inference A PREPRINT 10−4 10−2 100 102 100101102 10−4 10−2 100 102 100101102 10 haplotypes 100 haplotypes 1000 haplotypes 10000 haplotypes Time [s] Haplotype length [Mbp] Haplotype length [Mbp] 100101102 Memory [MByte] 10−2 100 102 Distance Hamming LeftM # Haplotypes 10 100 1000 10000 Figure 1: Time and memory needed for the construction of a coalescent tree as a function of the number of haplotypes, haplotype length, and distance function. Sampling is exact (t0= 1) and no bias is applied (c0= 1). Sampling with respect to either of the two distances yields identical memory usage since their computation does not involve memory allocation. Results are also available in table 1. 3 ARG Inference In Moonshine, ancestral recombination graphs are instances of the ARG type, a subtype of AbstractGenealogy. Since they represent a more realistic model for the ancestry of a sample subject to recombination, ARGs can be viewed as improved versions of coalescent trees. ARGs are constructed sequentially from a Tree by iteratively sampling recombination events until an ancestry consistent with the sample is reached. As is common, consistency is defined through the ISM: a genealogy is consistent with a sample if the number of mutations per marker is at most one. In a context where the recombination rate is several orders of magnitude higher than the mutation rate, a consistent ARG is typically a more realistic genealogy than any other kind of inconsistent ancestry. Moonshine is consequently very well suited to working with single nucleotide polymorphisms (SNPs), which have a low mutation rate. Ancestries are modified by recombination events, which partitions ancestral material into two subintervals: that to the left of an associated point, called a breakpoint, and the material to the right. From a graph-theoretical perspective, a recombination event is generally represented by a vertex of degree 3 with a single child and two parents. As an example, imagine that the child edge of a recombination vertex is ancestral for an interval Iand that the associated breakpoint is b. In that case, one of the parental edges, generally referred to as the left edge, is ancestral for [0, b)∩I 6 Moonshine.jl: ARG Inference A PREPRINT 0 5 10 Threshold 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.5 1.0 Speedup Prob. early terminat on Prob. correct Figure 2: Speedup in tree construction, probability of sampling the correct sequence and probability of early termination as a function of the sampling threshold t0. No bias is applied (c0= 1). The probability of early termination is 1−t0. As the number of candidate vertices increases, the probability of sampling correctly from the target distribution converges to t0(1 −log t0)(depicted here). The exact probability for a single run is given by lemma 1. Results are also available in table 2. 7 Moonshine.jl: ARG Inference A PREPRINT while the other (right) edge is ancestral for [b, ∞)∩I. Although we used the term “interval”, Ican actually be a union of intervals. We will continue to use this terminology when the distinction between the two concepts is not relevant. Recombination vertices are “added”, figuratively speaking, by deleting an edge from the graph and connecting the new vertex with both endpoints of the removed edge. The edge connected to the child vertex is the recombination vertex’s child edge, and the other is its left edge. This procedure leaves the right edge floating. Every recombination is immediately followed by the coalescence of the right edge with the graph. As it is carried out by a different algorithm than the coalescence of two dangling vertices encountered in temporal algorithms or, more trivially, when inferring a coalescent tree, we call those recoalescence events even though they result in an additional coalescence vertex as well. Coalescence and recombination vertices are in many regards mirror images of each other. A coalescence vertex has two child edges and one parental edge. If the child edges are ancestral for two intervals Iland Ir, then so is the parental edge for Il∪Ir. The standard recoalescence procedure begins again by deleting an edge, followed by connecting its incident vertices with the recoalescence vertex. The remaining child edge is then connected to the recombination vertex’s right edge, concluding the procedure. Recoalescence can also occur without edge deletion. In that case, the coalescence vertex becomes the new root of the graph and lacks a parental edge. It is connected downstream to the previous root and the recombination vertex. The type of a recoalescence event depends on its latitude, denoted lc, which itself depends on the recombination’s latitude lr. The root of an ARG corresponds to the sample’s most recent common ancestor (MRCA), sometimes called the grand MRCA (gMRCA). Its latitude is the time to the gMRCA (TgMRCA). Let vrand vcbe the recombination and recoalescence vertices and sr−drand sc−dcthe edges deleted in the recombination-and-recoalescence (RR) steps. Algorithm 2 summarizes the RR procedure. Algorithm 2: Recombination-and-recoalescence 1Add two vertices vrand vcto G 2Delete recombination edge sr—dr 3Add edges sr—vrand vr—dr 4if lc≤TgMRCA then 5Delete recoalescence edge sc—dc 6Add edges sc—vcand vc—dc 7else 8Add edge vc—mrca 9Add edge vc—vr 3.1 Unrestricted Recombination-and-Recoalescence events Moonshine has the capability to sample two kinds of RR events: restricted and unrestricted. Unrestricted events, the subject of this section, have a distribution designed to closely match that of the CWR. Restricted recombination events are sampled to reduce the total number of mutation events on an ARG; these will be discussed at length in the next section. Our package exports a method for sampling an arbitrary number of unrestricted recombination events. Standard theory [7] models the positions (on sequences) and locations (on ancestries) of recombination events as a Poisson point process (PPP). Conditional on their number, both locations and positions are distributed uniformly. When applied directly, this method has the drawback of generating sequences devoided of material ancestral for the sample at hand. Recombination events can be classified depending on whether ancestral material is present on both sides. If it is, the breakpoint can be positioned in ancestral or non-ancestral material. These are referred to as type 1 and type 2 events respectively and exclusively create haplotypes having ancestral material. Events positioned such that only material to their left or right side is ancestral are classified as type 3 and type 4 respectively. Finally, events occurring in entirely non-ancestral sequences are classified as type 5. It is desirable for an algorithm to only generate the first two types of recombination events since the other ones do not contribute to the structure of the sample. Our method follows this approach. We start by drawing a recombination edge with probability proportional to its length. The location of the event is distributed conditional on the selected edge. Then, we sample a position, also known as a breakpoint, uniformly on the mathematical closure of the set of intervals for which the recombination edge is ancestral. This strategy enforces the existence of ancestral material on both sides of the breakpoint and avoids recombination events of type 3, 4 and 5. Additionally, since the closure of ancestral intervals is not, in general, equal to the intervals themselves (it may contain “holes” of non-ancestral material), type 2 events are possible. 8 Moonshine.jl: ARG Inference A PREPRINT 7 6 21 3 4 5 (a) Initial graph. 4−2is the recombination edge in fig. 3c. 5−1is the recombination edge in fig. 3b and recoalescence edge in fig. 3c. 7 6 21 3 4 5 8 9 gMRCA (previous) (b) Recoalescence above the TgMRCA. The resulting ARG is taller than the original. 7 6 21 3 4 5 8 9 (c) Recoalescence below the TgMRCA. Figure 3: Two types of recombination events. Elements removed from the initial graph (fig. 3a) are in red while those added to it are in green. Selecting a conditional distribution for the location of recombination events on branches requires careful consideration. Theory dictates that it should be uniform. However, this results in a problematic frequency of locations close to the branches’ endpoints, contributing to the short branches issue discussed in subsection 3.2.3. Consequently, we use the location-scale family associated with the Beta(2,2) distribution instead. This method preserves the symmetry of the uniform distribution while favoring locations closer to the center of branches, achieving sufficient reduction in the number of short branches to prevent numerical error. The recoalescence process follows standard theory. The latitude is distributed as an inhomogeneous Poisson process with rate equal to the number of branches. The usual time-scale transformation strategy (see subsection 3.2.3) is implemented. A recoalescence edge is sampled uniformly, conditional on the latitude. The process of sampling a single unrestricted recombination event is summarized in algorithm 3. Algorithm 3: Unrestricted Recombination 1Sample a recombination edge in er∈Ewith probability proportional to its length 2Sample a breakpoint r∈(rL, rU)where rLand rUare the leftmost and rightmost positions for which eris ancestral 3Sample a recombination latitude lrdistributed as ldst(er)+ (lsrc(er)−ldst(er))Beta(2,2) where lsrc(er)and ldst(er)are the latitudes of the highest and lowest vertex adjacent to errespectively 4Sample a recoalescence latitude lcdistributed as an inhomogeneous PPP with rate equal to the number of branches via time-scale transformation 5Sample a recoalescence edge ecuniformly among the edges at latitude lc 6Apply algorithm 2 As a sidenote, since a function allowing users to add arbitrary recombination events to a graph is exposed, Moonshine can easily be used to simulate ancestries. We do not recommend doing so, however, especially not using the built-in types Tree and Arg, which are designed to implement ARG inference. Instances of these types store information irrelevant in a simulation context which greatly hinders performance. 9 Moonshine.jl: ARG Inference A PREPRINT 1 1 1 1 1 010 0 0 10 0 0 10 2 10 2 10 0 0 1 (a) Inconsistent graph: marker 2 mutates twice. 1 1 1 1 1 010 0 0 10 0 0 10 1 1 1 1 10 10 2 0 0 1 (b) Wild recoalescence event with a non-ancestral edge leading to a wild branch. Figure 6: Wild RR event leading to a reduction in the total number of mutations with recoalescence on non-ancestral edge. 0101010 0 10 0 0 2 0 0 2 0 0 2 0 0 1 0101010 0 10 0 0 100 1 01 01 0 0 1 2 Figure 7: A single recombination event leads to the elimination of multiple mutations. consideration when sampling the location of a recombination event. Each edge on which a recombination event would result in a reduction in the number of mutations for the current marker has a probability of being chosen proportional to the difference in latitude of its incident vertices, a sampling strategy described earlier as “constrained”. For reasons discussed earlier, the event’s location on the branch is not sampled uniformly but rather assumed to be distributed as the same location-scale family associated with the Beta distribution as the one used for the unconstrained case. Determining the location of a derived recoalescence event is fairly straightforward. According to standard theory ([7]), coalescence occurs at a unit rate with each admissible edge. We simply need to list all available locations and pick one uniformly at random. The set of possible recoalescence edges ERis readily established and does not require graph traversal. The minimum latitude is either that of the recombination event or the smallest latitude among the destination vertices of ER, whichever is greater. Each edge in ERis associated with a probability proportional to its length minus any section outside of admissible latitudes. Once the recoalescence branch is determined, a location is sampled as usual on its admissible portion. Simulating the location of a wild recoalescence event is more involved due to the semi-infinite nature of its support. Its latitude is distributed as the time of the first event of a non-homogeneous Poisson process. A common approach in one-dimensional scenarios such as ours is the time-scale transformation method, which is a form of inverse transform sampling and requires computation of the inverse of the integrated intensity function, also known as the cumulative intensity function. The main issue stems from our decision not to track the number of live edges by latitude. This improves general performance and reduces memory usage, but makes evaluation of the intensity function extremely time-consuming, as it requires graph traversal. We tackle this issue using numerical integration. Although somewhat variable, the intensity function is piecewise constant and therefore a very good candidate for quadrature. We use a logarithmic grid of latitudes to account for the generally decreasing complexity of the ARG topology as height increases. The intensity function is evaluated at quadrature nodes in a single partial traversal of the ARG, reducing overhead to the minimum. By default, our algorithm uses a grid of compile-time constant size 25, but this number can be tuned to balance precision with performance using the preference mechanism. 16 Moonshine.jl: ARG Inference A PREPRINT 10 0 11 11111 01 1 0 0 1 1111 10 0 1 1 01 1 0 0 4, 5 00000 1, 4, 5 2, 3 10 0 11 11111 01 1 0 0 1 1 1 1 1 10 0 1 1 01 1 0 0 00000 1, 4, 5 2, 3 Figure 8: Multiple crossing over. 3.2.4 Multiple Crossing Over In addition to RR events, Moonshine can sample multiple crossing over (MCO) events constrained to reduce the number of mutations. These events arise in the following scenario: assume that the branch on which a recombination event has been sampled is the right parental edge of a recombination vertex. Assume further that a recoalescence with the sibling of the other parental edge’s branch would reduce the number of mutations. When both of these conditions are met, RR events are unnecessary; the same effect can be achieved by modifying the ancestral intervals of the parental edges of the recombination vertex. For this to work, it is necessary to track the partition induced by recombination events. In an ARG, every recombination vertex is associated with two sets of intervals, one for each parental edge. We call these sets recombination masks and denote the mask associated with recombination kby mk={ml k, mr k}. Initially, before any MCO event involving a specific recombination vertex, the mask is rather simple. Let bbe the position of the associated recombination event. ml k= [0, b) mr k= [b, ∞). A MCO event at position b′> b transforms those sets as follows: ml k= [0, b)∪[b′,∞) mr k= [b, ∞)∩[0, b′) = [b, b′). A subsequent event at position b′′ > b′would yield ml k= [0, b)∪([b′,∞)∩[0, b′′)) = [0, b)∪[b′, b′′ ) mr k= [b, b′)∪[b′′,∞). In general, the correct mask can be computed by intersecting the rightmost interval (with respect to the right endpoint) with [0, b′)and taking the union with the other interval and [b′,∞). Although explicit storage is not necessary, coalescence vertices can be thought of as being associated with the mask [0,∞). Recombination masks are used to compute the ancestral intervals of parental edges in the following way: an edge’s ancestral interval is equal to the intersection of its recombination mask with the union of its children’s intervals. Conditions under which a MCO event can occur are very restrictive. Consequently, we expect detection to be underpowered. This might be improved in the future, for example, by allowing the user to bias the sampler towards these events. A MCO is illustrated in fig. 8. 3.3 ARG Update RR and MCO events both modify sequences and ancestral intervals associated with vertices and edges upstream. The affected elements must be updated immediately to ensure the soundness of the remainder of the procedure, which can be achieved by traversing the ARG from the recombination and recoalescence edges toward the root. This seemingly simple operation is, however, computationally demanding. Our first, rather naive implementation was a major bottleneck of the constrained recombination algorithm. To make it as efficient as possible, we limit the update procedure to the elements affected by the event. Recombination events are mostly local, meaning they can only affect vertices and edges located upstream. Consequently, any element not upstream of either the recombination or recoalescence edge can be ignored when updating. In fact, an element must be located upstream of a modified element to be modified itself. This means that the number of updates can be further reduced by keeping track of the state of every vertex and edge before they are updated and stopping when a match between original and updated versions is detected. The algorithm terminates early if both the sequence and the ancestral interval associated with an edge are left unchanged. 17 Moonshine.jl: ARG Inference A PREPRINT The early termination strategy described above dramatically decreases the time dedicated to ARG update. Indeed, coalescences with vertices left untouched limit the spread of changes, often well below the root. It might be conceived that the additional costs associated with storing information about elements before updating them would outweigh the benefits of reducing the number of updated elements. It turns out that a very significant number of recombination events are very local in nature, to such a degree that we have yet to find the point of diminishing return for reasonably large samples. That is not to say, however, that the procedure cannot be improved further. We were able to squeeze even more performance out of it through hashing. Instead of making a copy of the current sequence and ancestral interval before update, we simply compute a hash value for each of those, which we promptly hash together. The procedure is terminated early if the hash of the updated sequence-ancestral interval pair is equal to the original. Since hash functions are not injective, there is a possibility that different pairs may have the same hash, an event known as a collision. Fortunately, this has not been a problem in practice. To put any doubt to rest, the method validate can be applied to the final product of the ARG inference process to ensure that it is exempt, among other things, of the inconsistencies that would emerge from collisions. The complete procedure is presented in algorithm 6. Note that since recombination vertices have a unique child, their associated sequence can be a reference to their child’s. Consequently, line 5 can be performed without additional allocation, dramatically cutting down on memory usage. The early termination strategy is implemented by line 16. Algorithm 6: ARG Update 1Set stack ←[starting_vertex] 2while stack is not empty do 3Set v←pop!(stack) 4if vis a recombination vertex then 5Set the sequence of vto that of its child 6Set the ancestral interval of its parent edges to the intersection of its child’s with the appropriate ancestral mask 7push!(stack,parents(v)) 8else 9Set h←hash(haplotype(v)) 10 Set the sequence of vto the conjunction of its children’s haplotypes 11 if vis the root then 12 Continue 13 Set eto the parent edge of v 14 Update h←hash(h,ancestral_interval(e)) 15 Set the ancestral interval of eto the union of the ancestral intervals of v’s child edges 16 if h6=hash(hash(haplotype(v)),ancestral_interval(e)) then 17 push!(stack,parent(v)) 3.4 Technical Considerations 3.4.1 Markov Approximation The procedure described in this section is exact in the following sense: at any step, recombination and recoalescence edges are free to be sampled anywhere on Gsubject only to the mutation number reduction constraint. It is more similar to the original Wiuf-Hein sequential algorithm than any of its Markovian approximations, such as the SMC or SMC’. This makes our sampler more realistic, but it comes with a performance penalty: as RR events are incorporated into the graph, the number of edges that need to be considered for the next events increases. Consequently, we should expect the number of operations performed by our algorithm to grow as it progresses along sequences. One way to alleviate this computational burden is to make the recombination process artificially Markovian. The procedure we just presented allows for such an approximation, although it was omitted from the description for simplicity. Both restricted and unrestricted RR events may be executed inside a window moving across the sequences, the width of which can be selected by the user. A width of 0 corresponds to a first-order Markovian approximation akin to the SMC, while an infinite width means no approximation at all. In general, assuming distances in base pairs (bp), specifying a width of wfor an RR event occuring at position phas the effect of excluding any edge esuch that ancestral_interval(e)∩[p−w, p +w] = ∅. 18 Moonshine.jl: ARG Inference A PREPRINT 10−6 10−4 10−2 100 10−6 10−4 10−2 100 10 haplotypes 100 haplotypes 1000 haplotypes 10000 haplotypes Time [hr] 1 10 100 1.0 1.5 2.0 1 10 100 Speedup Haplotype Length [Mbp] 0.0 kbp 100.0 kbp Window Width 0.0 kbp 100.0 kbp Inf kbp # Haplotypes 10 100 1000 10000 Figure 9: Time required to sample a single consistent ARG as a function of the number of haplotypes and sequence length, including time spent inferring the initial tree. The speedup provided by various finite window widths is depicted in the lower figure. Technical details are discussed in subsection 5.1. from the set of candidate recombination and recoalescence edges. We emphasize that the window is centered on p, as our sampler is designed for the general task of rendering ARGs consistent rather than building them from a tree in a left-to-right sweep. Moreover, when sampling the recoalescence latitude of an unconstrained event, ignored edges are not taken into account for the computation of the rate of the recombination latitude. Figure 9 shows that the performance impact of reducing the window’s width is significant. For some combinations of sample size and haplotype length, a window width of 0 can lead to over twice as fast computation. Interestingly, the speedup is similar for a width of 100 kbp, suggesting that even relatively modest approximations can lead to substantial reduction in computational resources. Irrespective of window size, fig. 9 suggests that computation time is exponential in both the number of sampled haplotypes and markers. Similarly, fig. 10 suggests exponential growth in memory usage. It appears, however, that reducing the window size increases the size of the resulting ARG. This is largely explained by the increasing number of recombination events sampled as a result. As shown in fig. 11, this number tends to increase as the window width diminishes. Part of this phenomenon might be explained by increased flexibility. Larger window sizes increase the number of candidate edges for recombination and recoalescence events, allowing for constrained coalescence of more similar haplotypes and, ultimately, more parsimonious ancestries. This behavior is all the more interesting as it is entirely spontaneous; 19 Moonshine.jl: ARG Inference A PREPRINT 100 102 104 100 102 104 10 haplotypes 100 haplotypes 1000 haplotypes 10000 haplotypes Size [MByte] 1 10 100 0.85 0.90 0.95 1.00 1 10 100 Spaceup Haplotype Length [Mbp] 0.0 kbp 100.0 kbp Window Width 0.0 kbp 100.0 kbp Inf kbp # Haplotypes 10 100 1000 10000 Figure 10: Size of the object containing the inferred ancestral recombination graph as a function of the number of haplotypes and sequence length. The spaceup is the ratio of the size for infinite window width versus that for the indicated width. Technical details are discussed in subsection 5.1. We emphasize that the total amount of memory required by the sampling procedure exceeds quantities reported here. markers to the right of the recombination breakpoint are not taken into account in the recoalescence edge sampling procedure. 3.4.2 Markers and positions In addition to the ARG update procedure described at the beginning of this section, another major bottleneck in our method was, somewhat unexpectedly, a function called postoidx. It is the inverse of another function called idxtopos, which is designed solely to return the position of a given marker, indexed from the leftmost one. The term “inverse” in that context is used loosely, as only a few positions are associated with a sampled marker. For that reason, we define postoidx formally as a pseudoinverse of idxtopos: postoidx(p) = sup{i∈1,...,s:idxtopos(i)≤p}. Since it is often necessary to mask sequences with respect to an ancestral interval, this function is used extensively within Moonshine’s recoalescence and ARG update procedures. The endpoints of any given interval rarely correspond to the position of a marker. This is compounded by the fact that what we refer to as “ancestral interval” is, in reality, the union of multiple disjoint intervals. Overall, a simple linear search is insufficient; a more efficient approach is 20 Moonshine.jl: ARG Inference A PREPRINT 103.0 104.5 106.0 103.0 104.5 106.0 10 haplotypes 100 haplotypes 1000 haplotypes 10000 haplotypes # Recombinations 1 10 100 0.85 0.90 0.95 1.00 1.05 1 10 100 Recup Haplotype Length [Mbp] 0.0 kbp 100.0 kbp Window Width 0.0 kbp 100.0 kbp Inf kbp # Haplotypes 10 100 1000 10000 Figure 11: Number of recombination events sampled by the ARG inference procedure as a function of the number of haplotypes and sequence length. The recup is the ratio of that number for an infinite window versus that for the indicated width. Technical details are discussed in subsection 5.1. called for. Binary search yields considerable gains, but we can do better. Memoization is a standard strategy for solving similar problems, but it is detrimental to our use case. The cost of computing idxtopos via binary search is small, even in comparison with a lookup in a hash table. Besides, positions are encoded as double-precision floatingpoint values, which limit their usefulness as keys of an associative structure. Perhaps it would be possible to use locality-sensitive hashing to circumvent this issue. We did not explore this avenue, though, because the approach currently implemented, although relatively simple, proved satisfactory. It is essentially a slightly optimized version of interpolation-sequential search ([35]). We take advantage of the fact that the markers’ position vector is not modified by ARG inference. Instead of directly computing the inverse of idxtopos, we do so on a first-order approximation obtained by least-squares fitting, which is straightforward. Even though computing such an approximation is timeconsuming relative to an iteration of bisection search, it only has to be done once, making the overhead negligible. Pseudocode for postoidx is given in algorithm 7. While algorithm 7 is relatively simple, it takes advantage of three features of the positions vector: its static nature, its monotonicity, and the approximate regularity of the distance between markers. A more efficient procedure could likely be designed, but algorithm 7 is at least efficient enough to eliminate the postoidx bottleneck. 21 Moonshine.jl: ARG Inference A PREPRINT Algorithm 7: postoidx // Assume idxtopos(i)≈ai +b 1Function postoidx(p) 2Set i←p−b a,l←1,r←s 3Clamp iin [1,s] 4Update lor rthrough a single binary search iteration 5Perform interpolation-sequential search with i,land r 3.4.3 Memory Allocation Considerable effort has been devoted to using memory as efficiently as possible in the computationally demanding sections of Moonshine. In addition to the obvious benefits in terms of reduced memory requirements, we found dynamic allocation to be a major performance bottleneck in itself. Although Julia is a high-level programming language, fine-grained memory management is straightforward. It is entirely possible to manage allocations directly. Part of the C standard library, comprising functions such as malloc,calloc and free, is exposed to the user via the Libc module included in Julia’s standard library. Our approach, however, is slightly higher level: we make extensive use of slab allocation ([36]) through Bumper.jl ([37]). This package allows us to handle raw pointers when needed while efficiently managing the underlying memory pool automatically. This allows for nearly effortless tight memory management and dramatically reduces the number of allocations needed in core methods. Performancecritical methods accept the buffer keyword argument through which the caller can pass a Bumper.jl buffer. This is how the memory is allocated by the build! method for inferring ARG for instance. The same buffer can be passed around to different function calls for maximum memory usage efficiency. This also makes it extremely easy to seamlessly integrate Moonshine methods in a workflow that already takes advantage of Bumper.jl facilities. 3.5 Complete Algorithm The complete ARG inference procedure is given in algorithm 8. 4 Conclusion We have described a fast and efficient algorithm for ancestral recombination graph inference and the broad outlines of its Julia implementation via the Moonshine package. In addition to exact sampling, an approximate scheme based on user-defined genomic windows is available, which reduces inference time. Inference is grounded in a restricted version of the coalescent-with-recombination distribution rather than heuristics. We have also presented a flexible algorithm for coalescent tree construction that can accommodate a wide range of metrics on haplotypes, exactly or approximately, while avoiding numerical instabilities. As of version 0.3.3, these algorithms are for the most part mature. No significant leap in inference time or memory usage is to be expected in the near future. Our goal for version 1.0.0 is to improve real-world usability. Specifically, the API reference needs to be completed to make implementing alternative models a more pleasant experience. General documentation and guides also need to be improved for the general user. Both are publicly available at https://moonshine.patrickfournier.ca. As of today, version 0.3.4 and later are compatible with files encoded in the variant call format (VCF) through the VCFTools.jl package ([38]). We plan to support FASTA/FASTQ and sequence/binary alignment map format (SAM/BAM) via the FASTX.jl and XAM.jl packages, respectively. Along with compatibility with these input formats, we plan to support missing markers by treating them as non-ancestral material, which is the approach taken by [39, 28, 40]. Moonshine is available on Julia’s general registry and can be easily and quickly obtained using standard facilities. In addition, we plan to package it in two application containers. One will include the Pluto notebook ([41]) package and graphical utilities such as Makie.jl ([42]). Its goal is to reduce the burden associated with ARG inference for practitioners. The other container will be more minimal, including only the minimum required to deploy Moonshine via an orchestration system such as kubernetes. Our hope is to facilitate the creation of high-performance computation clusters in small-to-medium environments. 22 Moonshine.jl: ARG Inference A PREPRINT Algorithm 8: ARG Inference 1Set markeridx ←1,live_edges ←∅ 2Update markeridx and live_edges via algorithm 5 3while markeridx > 0do 4while |live_edges|>1do 5Sample e1&e2uniformly from live_edges and remove them // Consider the possibility of a wild recombination 6Set markerpos ←idxtopos(markeridx) 7if parent(src(e1),markerpos)=src(e2)then 8Set redge ←src(e1)−sibling(dst(e1),markerpos) 9Set vupdate ←src(e2) 10 else if parent(src(e2),markerpos)=src(e1)then 11 Set redge ←src(e2)−sibling(dst(e2),markerpos) 12 Set vupdate ←src(e1) 13 Sample a recombination edge redge with probability proportional to its length among all admissible edges (see subsection 3.2.3 for admissibility criteria) 14 Sample a recombination latitude rlat betwen dst(redge)and src(redge)from a position-scale family associated with the Beta(2, 2) distribution 15 if Recombination is on a wild edge then 16 Sample a recoalescence latitude clat distributed as a non-homogeneous Poisson process on [rlat,∞)via empirical supremum rejection sampling (see subsection 3.2.3 for complete distribution) 17 Sample a recoalescence edge cedge uniformly among those intersecting latitude clat 18 else 19 Fill possible_cedges with derived edges located above redge 20 for k∈1,...,|possible_cedges|do 21 Set ras to the latitude of the lowest valid recoalescence edge associated with eaccording to subsection 3.2.3 22 Set minlatitudes ←max{rlat,r} 23 Set ubound ←max{(latitude ◦src)(e) : e∈possible_cedges} 24 Set lbound ←min minlatitudes 25 Sample recoalescence latitude clat distributed as a non-homogeneous Poisson process on [lbound,ubound]via rejection sampling 26 Sample associated recoalescence edge cedge (see subsection 3.2.3 for complete procedure) 27 Set vupdate to the number of vertices plus 2 28 Sample breakpoint conditional on redge,cedge and markeridx (see subsection 3.2.2) 29 if src(cedge)∈parents(dst(redge)) and dst(redge)is a recombination vertex then 30 Perform an MCO event on vertex dst(redge) 31 else 32 Sample a recombination event on redge at rlat followed by a recoalescence event on cedge at clat at position breakpoint 33 Update ARG via algorithm 6 starting at vupdate 34 Update live_edges 35 Update markeridx and live_edges via algorithm 5 23 Moonshine.jl: ARG Inference A PREPRINT 5 Appendix 5.1 Simulations All simulations were performed with Moonshine version 0.3.3 and Julia 1.11.3 on AMD EPYC 9654 processors. Execution time was measured using the @timed macro. To reduce variability due to the execution environment and differences between ARGs, each reported measurement is the minimum of three runs, averaged over five graphs. Haplotypes are generated by msprime using StandardCoalescent and BinaryMutationModel as the ancestry and mutation models, respectively. Per-locus mutation and recombination rates were both set to 10−8. An effective population size of 104was specified. The calls to msprime.sim_ancestry and msprime.sim_mutations were made through the Moonshine interface to msprime.Moonshine 0.3.3 does not implement the recoalescence latitude sampling procedures described in subsection 3.2.3. In particular, its algorithm for wild events is less performant. Consequently, we expect the performance of versions 0.3.6 onwards, which implement the time-scale transformation approach, to be slightly better when using the default grid. All the code used for the simulation studies, as well as the raw data, is publicly available on Codeberg (https://codeberg.org/ptrk/moonshine.jl-papers). For convenient reproducibility, code to execute each simulation presented and produce related figures is grouped in a single Pluto notebook. As the simulations are computationally intensive and require considerable time to complete, it may be more convenient to run them on remote machines. For that purpose, the notebook can be executed as a standalone Julia script. Figure generation steps will be skipped when doing so; they can be run locally later on. 5.1.1 Results Raw simulation results are presented below. The following abbreviations are used in the column names: •n: number of haplotypes in the sample •L: length of sampled haplotypes •SeU: speedup •SaU: spaceup •RU: recup •W: window width •Recs: recombination events Numbers of markers and recombination events are averaged over the sampled ARGs. Table 1: Tree Building Time n L [Mbp] Distance # Markers Time Size [MByte 10 1 Hamming 588.4 <1s<1 10 1 LeftM 588.4 <1s<1 10 10 Hamming 5725.4 <1s<1 10 10 LeftM 5725.4 <1s<1 10 100 Hamming 56517.4 <1s<1 10 100 LeftM 56517.4 <1s<1 10 250 Hamming 141851 <1s 1 10 250 LeftM 141851 <1s 1 100 1 Hamming 1035.2 <1s<1 100 1 LeftM 1035.2 <1s<1 100 10 Hamming 10435.4 <1s<1 100 10 LeftM 10435.4 <1s<1 100 100 Hamming 103598 <1s 3 100 100 LeftM 103598 <1s 3 100 250 Hamming 259113 <1s 8 100 250 LeftM 259113 <1s 8 1000 1 Hamming 1539.8 <1s<1 1000 1 LeftM 1539.8 <1s<1 24 Moonshine.jl: ARG Inference A PREPRINT 1000 10 Hamming 14950.2 <1s 4 1000 10 LeftM 14950.2 <1s 4 1000 100 Hamming 149768 <1s 38 1000 100 LeftM 149768 <1s 38 1000 250 Hamming 375008 00:00:01 96 1000 250 LeftM 375008 <1s 96 10000 1 Hamming 1928.2 00:00:04 6 10000 1 LeftM 1928.2 00:00:02 6 10000 10 Hamming 19569.2 00:00:15 50 10000 10 LeftM 19569.2 00:00:02 50 10000 100 Hamming 195747 00:02:05 492 10000 100 LeftM 195747 00:00:03 492 10000 250 Hamming 489204 00:05:13 1228 10000 250 LeftM 489204 00:00:03 1228 Table 2: Tree Sampling Threshold Threshold # Markers Time Size [MByte] SeU 0.1 491025 00:00:24 1228 10.1 0.2 491025 00:00:56 1228 4.4 0.3 491025 00:01:30 1228 2.7 0.4 491025 00:02:02 1228 2 0.5 491025 00:02:37 1228 1.6 0.6 491025 00:02:56 1228 1.4 0.7 491025 00:03:41 1228 1.1 0.8 491025 00:04:05 1228 1 0.9 491025 00:04:16 1228 1 1 491025 00:04:06 1228 — Speedups are computed with respect to a threshold of 1. Table 3: ARG Window Width n L [Mbp] W [kbp] # Markers # Recs. Time Size [MByte] SeU SaU RU 10 1 0 549.2 171 <1s<11.2 1 1.1 10 1 100 549.2 176 <1s<11.1 1 1 10 1 ∞549.2 183 <1s<1— — — 10 10 0 5694 1824 <1s 1 1.1 1 1 10 10 100 5694 1875 <1s 1 1.1 1 1 10 10 ∞5694 1811 <1s 1 — — — 10 100 0 56502 19601 00:00:04 145 1.1 1 1 10 100 100 56502 19569 00:00:04 145 1.2 1 1 10 100 ∞56502 19591 00:00:05 144 — — — 10 250 0 141770 49098 00:00:27 893 1 1 1 10 250 100 141770 49210 00:00:25 894 1 1 1 10 250 ∞141770 49005 00:00:26 891 — — — 100 1 0 1064 1389 <1s<11.3 1 1 100 1 100 1064 1501 <1s<11.2 1 1 100 1 ∞1064 1455 <1s<1— — — 100 10 0 10470 15163 <1s 23 1.8 1 1 100 10 100 10470 15108 <1s 22 1.7 1 1 100 10 ∞10470 15054 00:00:01 22 — — — 100 100 0 103530 149880 00:01:04 1989 1.7 1 1 100 100 100 103530 149580 00:01:07 1987 1.7 1 1 100 100 ∞103530 148700 00:01:52 1977 — — — 100 250 0 259080 375660 00:07:29 12403 1.5 1 1 25