phastSim

Scenario 1: simulating MSAs that mimic a user-provided MSA (Figure 1B) A common use-case for alignment simulation software is that users want to simulate an MSA that mimics the evolutionary history of a given MSA, for example, because this is needed for parametric bootstrap analysis. Until now, this required at least a two-step process whereby users first inferred the tree and model in one piece of software, then used these as input to the MSA simulation tool. The resulting MSA often failed to capture many characteristics of the original MSA, such as the position of gaps and the sitespecific evolutionary rates. AliSim improves this process by first running IQ-TREE to infer an evolutionary model and a tree from the input MSA, and then immediately generating any number of simulated MSAs from the inferred tree and model with the same gap patterns of the original MSA. For simulations under a mixture model, AliSim randomly assigns a model component of the mixture to each site according to the site posterior probability distribution of the mixture. For site-frequency mixture models, AliSim applies the posterior mean site frequencies (Wang et al. 2018) . Similarly, AliSim employs the posterior mean site rates to better reflect the underlying evolutionary rate variation across sites. All these mechanisms help produce simulated MSAs that better reflect the relevant features of the original MSAs.

Scenario 2: simulating MSAs from a random tree and/ or random parameters from empirical/user-defined distributions When using Seq-Gen or Dawg, users need to provide as input a tree with branch lengths. To avoid this sometimes cumbersome step, AliSim allows users to generate a random tree under biologically plausible models such as the Birth-Death model (Kendall 1948) and the Yule-Harding model (Yule 1925; Harding 1971) . For the Yule-Harding model, users only need to specify the number of leaves of the tree. For the birth-death model, users need to additionally provide the speciation and extinction rate. Branch lengths are randomly generated from an exponential distribution with a mean of 0.1 (by default) or from a userdefined distribution specified by a list of numbers.

Where users wish to simulate alignments that mimic empirical MSAs in the absence of a set of input alignments, AliSim can generate a stationary distribution of nucleotides from empirical distributions that were previously estimated from a large collection of empirical datasets (Naser-Khdour et al. 2021). Other parameters, such as substitution rates, nonsynonymous/synonymous rate ratios, transition, and transversion rates, can be drawn from user-defined lists of numbers, allowing AliSim to incorporate arbitrary distributions for all simulation parameters.

Existing simulators typically employ either the rate matrix approach or the probability matrix approach to evolve sequences along a tree (see Methods). However, their performance varies with different sequence lengths ( ) and branch lengths ( ). Therefore, AliSim automatically switches between the rate matrix and probability matrix approaches to optimize the computing time. To determine when to use each approach, we compared the runtime of the rate matrix approach with the probability matrix approach on simulations using different combinations of

and (see Methods).

The simulation results showed that the rate matrix approach is generally faster than the probability matrix approach when 7 <2 .226 and 7 <17.307 for the discrete and continuous rate heterogeneity models, respectively. Therefore, the adaptive approach will apply the rate matrix approach for those branches that satisfy this inequality; otherwise, it will apply the probability matrix approach.

of common phylogenomic simulation conditions. Specifically, we simulated 8deep9 data with 30K sites and 10K to 1M sequences, and we simulated 8long9 data with 30K sequences and 10K to 1M sites (see Methods). We note before presenting these results that phastSim is designed specifically to simulate data along trees with a large proportion of extremely short branches. The trees in these simulations do not match these conditions, and so one might expect phastSim to perform poorly here. We include phastSim here for completeness, and present a comparison of phastSim and AliSim under the conditions for which phastSim was designed later.

As expected, when the datasets were small, runtimes and memory usage were similar across all pieces of software. However, AliSim shows increasing advantages in runtimes and memory usage as the datasets get bigger. For example, for the deepest dataset (30K sites and 1M sequences), Seq-Gen, Dawg, INDELible, phastSim, and AliSim required 1.9, 1.7, 4.8, >24, and 1.3 hours of runtime, respectively. A nd for the longest dataset (30K sequences and 1M sites), Seq-Gen, Dawg, INDELible, phastSim, and AliSim required 5.1, 16.4, 3.6, >24, and 1.3 hours respectively. Thus, AliSim is a fast sequence simulator under a range of conditions. Runtimes for (A) 8deep9-data (varying number of sequences and 30K sites) and (B) 8long9-data simulations (varying number of sites and 30K sequences) and peak memory consumption for (C) 8deep9data and (D) 8long9-data simulations.

AliSim shows dramatic improvements over other software in peak memory usage. Figures 2C and 2D show that these improvements become large even for fairly modestly-sized datasets. For the deepest dataset (30K sites and 1M sequences), Seq-Gen, Dawg, INDELible, and AliSim required 56.4, 56.3, 502.3, and 1.3 GB RAM, respectively (Figure 2C; phastSim peak memory usage was not recorded as it took >24hrs to run). For the longest dataset (30K sequences and 1M sites), Seq-Gen, Dawg, INDELible, and AliSim consumed 56, 537, 499, and 0.2 GB RAM, r espectively (as above, phastSim was excluded because it took >24hrs to run). Importantly, the memory usage of AliSim only grows sub-linearly with respect to the dataset size. For the deep simulations, when increasing the number of sequences from 10K to 1M (100-fold), the RAM consumption of AliSim only increased from 137MB to 1.3GB (~10-fold increase, Figure 2C). For the long simulations, increasing the sequence length from 10K to 1M sites (100fold) only increased the RAM usage marginally from 156MB to 222MB (less than a 2-fold increase, memory consumption is now dominated by storing the tree data structure and not by the simulated sequences.

We also tested the performance of these tools on simulating MSAs from COVID-19-like trees, which differ from commonly-simulated trees because they contain a large proportion of extremely short branch lengths, and for which phastSim was explicitly designed. In these conditions (Suppl. Figure S1A), phastSim and AliSim were the two fastest pieces of software, requiring 6 and 9 minutes respectively, whereas Seq-Gen, Dawg, and INDELible took 1.7, 1.7, and 3.2 hours respectively. In terms of RAM consumption, phastSim and AliSim only needed 1.4 and 1.3 GB RAM, respectively, while Seq-Gen, Dawg, and INDELible required 56, 51, and 502 GB RAM (Suppl. Figure S1B). We note that the performance of phastSim in these conditions is particularly remarkable because it is written in Python. Because of this, it may be that the language itself rather than the sequence simulation algorithms of phastSim are what limit its performance, and we speculate that phastSim may be able to be even faster if re-written in C or C++.

The adaptive approach helps AliSim achieve high performance by selecting the most efficient simulation approach for each branch. For example, in simulations under trees where branch lengths were generated from an exponential distribution with a mean of 0.1, the adaptive method will apply the probability matrix approach rather than the rate matrix approach on most branches, simply because most branches are longer than the switching parameters. The benefits of the adaptive approach can be measured in our simulations by forcing AliSim to use one method. For example, using the adaptive approach, AliSim took only 1.3 hours to simulate 1M sequences of 30K sites (Figure 2A). However, if we force AliSim to employ only the rate matrix approach, it takes more than 5 hours to simulate the same dataset. Similarly, the adaptive approach took only 9 minutes to simulate a dataset on a COVID-19-like tree (Suppl. Figure S1A), but if we force AliSim to use the probability matrix approach, the same simulation takes 1.3 hours.

To validate the AliSim, we simulated 287 MSAs with 100 sequences across a wide range of substitution models and insertion-deletion rates of 0.0, 0.02, 0.04, 0.06, 0.08 and 0.1. These choices of indel rat es follow empirical studies (Cartwright 2009). We then ran IQ-TREE to determine the best-fit model using ModelFinder (Kalyaanamoorthy et al. 2017) and reconstructed phylogenetic trees under the best-fit model. We compared the topology between the true trees and the inferred trees using the RobinsonFoulds distance (Robinson and Foulds 1981).

Supplementary Table S1 showed that in 147 tests (51.22%), the true model was recovered as the best-fit model. In 243 tests (84.67%), 246 tests (85.71%), a nd 267 tests (93.03%), the true models appear in th e top-2, top-3, and top-4 best models, respectively. The average Robinson-Foulds distance between the true trees and the inferred trees across all test cases was 1.99 (s.e. 0.133). That means the inferred tree s differed from the true trees by only 0.995 out of 9 7 (1.03%) internal branches. The tree lengths (sum of branch lengths) of the inferred trees differed from the true trees by only 1.9%.

For simulations with non-zero indel rates, the average differences in the alignment length and proportion of gaps between MSAs simulated by AliSim and those by INDELible were 0.5% and 0.22%, respectively (Suppl. Table S2).

In conclusion, AliSim is a fast and memory-efficient simulation tool, which simplifies and speeds up many common workflows in phylogenetics. Because AliSim has a small memory footprint, it can be used to simulate even very large alignments on personal computers.

Naively, when simulating sequences along a bifurcating tree with tips, we need to store up to (2 3 1) sequences, consisting of ( 2 1) internal nodes and tips. To save memory, we release the memory allocated to the sequence of an internal node if the sequences of its children nodes are already generated. Additionally, in simulations without partitions and indels, AliSim writes out the tip sequences to the output file immediately after simulating them, then frees the memory. This approach considerably reduces the maximum number of sequences that need to be stored in memory from (2 2 1) to the maximum depth of the tree. For a balanced bifurcating tree, this maximal depth is 6( ) + 1, leading to a substantial reduction in memory usage.

We benchmarked AliSim against Seq-Gen, INDELible, Dawg, and phastSim. For INDELible, we used Method 2 which is more efficient than Method 1 in simulations with continuous rate heterogeneity across sites (Fletcher and Yang 2009). The benchmark was r un on a Linux server with 2.0 GHz AMD EPYC 7501 32-Core Processor and 1-TB RAM. Inspired by the newly emerged COVID-19 data, we tested the ability of all tools to simulate alignments with 30K sites and 10K, 100K, 500K, and 1M DNA sequences. For the model of sequence evolution, we applied the GTR model with a 0.2 proportion of invariant sites and a continuous Gamma model of rate heterogeneity across sites (shape parameter of 0.5). We ran different software to simulate sequences along a random Yule-Harding tree with exponentially distributed branch lengths with a mean of 0.1. We call this the 8deep9-data simulation. Moreover, to mimic the size of real phylogenomic datasets, we simulated MSAs with 30K sequences and increased sequence length from 10K to 1M sites. This is called 8long9-data simulation.