The field of phylogenetic inference deals with the problem of inferring hypotheses about the evolutionary relationships of species, mainly based on molecular sequence data. The inferred hypotheses are typically presented in the form of phylogenetic trees. There exist multiple approaches and optimality criteria to address this problem, such as distance based approaches (e.g., Neighbor Joining (NJ) ( 1 )), Maximum Parsimony (MP) ( 2 ), Maximum Likelihood methods (ML) ( 3 ) or Bayesian inference methods (BI) ( 4–8 ). At present, the vast majority of phylogenetic analyses are conducted via ML and BI methods that both rely on the phylogenetic likelihood function that comprises an explicit stochastic model of sequence evolution. The inference of ML optimal trees is computationally hard ( 9 ) due to the super-exponential increase of possible unrooted binary tree topologies as a function of the number of taxa/sequences in the tree. Hence, there exists a plethora of inference tools which deploy different heuristic search strategies to find a tree with a good likelihood. Some of the most widely used tools are FastTree2 ( 10 ), IQTREE2 ( 11, 12 ), and RAxML-NG ( 13 ).

To assess the accuracy of the inferred trees and the performance of the respective heuristics, comparisons against alternative tools are typically based upon a combination of empirical and simulated datasets. The results of such performance assessments may deviate depending on the datasets used. To the best of our knowledge, there does not exist a standardized approach or set of benchmark data to assess and compare these tools. Instead, simulations and empirical datasets are conducted and selected predominantly in an ad hoc manner. Typical informal selection criteria for empirical datasets are the number of sites, the number of sequences (taxa), and the percentage of gaps or missing data in the datasets (see, for example, the respective papers describing FastTree2, IQTREE2, and RAxML-NG).

We examined the papers of the aforementioned analysis tools as well as a study comparing FastTree and RAxML ( 14 ) to determine the criteria deployed to select the respective benchmarking datasets. Then, we queried the RAxMLGrove database ( 15 )(RG) for datasets satisfying these criteria. We found that at most 27% of datasets present in RG satisfy those criteria (see Table 1) and observe that the criteria do not appear to be representative if the data typically being analyzed. Related work assessing the inference accuracy of RAxML ( 16 ) and FastTree ( 17 ) has been conducted by Liu et al. ( 14 ). The authors used 1, 800 simulated and 10 empirical datasets and concluded, that RAxML (the predecessor of the re-designed code RAxML-NG that implements essentially the same algorithm) generally outperforms FastTree in terms of topological accuracy on smaller datasets, but that the differences diminish on large datasets with respect to the number of taxa. We make a similar observation for the updated versions of the two ML tools based upon around 5, 000 empirical datasets from TreeBase. Further, we can characterize and classify these datasets in a systematic manner by predicting their difficulty scores via the Pythia ( 18 ) tool. Pythia assigns values between 0 (easy/strong signal) and 1 (difficult/weak signal). Here, difficulty essentially quantifies the signal strength in the data. Zhou et al. ( 19 ) perform a large-scale analysis of ML-based methods on 19 empirical datasets with up to 200 taxa and partially thousands of genes. They find, that RAxML and IQ-TREE perform similarly in terms of topological accuracy on single-gene datasets, while FastTree performs substantially worse. In our study, we confirm this result for empirical datasets (composed of single- and ca. 10% multi-gene data) with a difficulty score below 0.7. The key contribution of our work is that we deploy empirical and simulated DNA data that representatively covers at least 75% of the data in the two databases (TreeBase and RAxMLGrove). We assess accuracy by means of the inferred tree topologies, the ML scores, and their statistical plausibility (i.e., whether the trees differ significantly from the bestknown tree).

We did not include amino-acid (AA) datasets here to reduce the computational burden and CO2 footprint of our analyses. We conducted our experiments using a Snakemake ( 20 ) pipeline, which is available at https://github.com/angtft/PhyloSmew. This pipeline can be easily extended by additional phylogenetic inference tools and will thus contribute to conducting accuracy inference studies in a more standardized way.

In this section, we first briefly describe the RAxMLGrove database ( 15 ) that we use for simulating Multiple Sequence Alignments (MSAs). Then, we describe the methods for selecting and simulating datasets. We undertake an effort to create realistic simulated MSAs, especially in terms of realistic gap patterns, and we evaluate this attempt in the respective Section on simulation quality. Then, we briefly describe our selection strategy of empirical datasets from the TreeBASE database ( 21 ).

RAxMLGrove Database. Throughout our experiments, we extensively use the RAxMLGrove (RG) database. RG contains anonymized phylogenetic trees inferred by users of RAxML/RAxML-NG web servers at the San Diego Supercomputer Center ( 22 ) and the SIB Swiss Institute of Bioinformatics (https://raxml-ng.vital-it.ch) with their respective inferred model parameters and program execution logs. At present, RG comprises more than 70, 000 datasets and is continously growing. The RAxMLGroveScript repository (https://github.com/angtft/RAxMLGroveScripts) provides functionality for simple data access using an SQLitedatabase with a corresponding Python script (henceforth called RGS). In addition to dataset queries based on a set of dataset characteristics (e.g., the number of taxa, or branch length variance of the tree), RGS provides the generate option to download datasets and generate simulated MSAs based on the estimated model parameters of that dataset. The simulations are conducted by executing Dawg ( 23 ) or AliSim ( 24 ). Since simulation is the main RGS functionality of interest in this work, we briefly list and describe the contents of one such downloaded simulated dataset, after using the generate command, in Table 2.

Simulated Data. We generate simulated data based on the RG database, version 0.7. As our goal is to assess ML tool inference behavior for representative datasets, we proceeded as follows to select the datasets for our analysis: 1. For the sake of simplicity, we only selected trees inferred on DNA alignments under the General Time Reversible (GTR) ( 25 ) model. 2. We selected datasets from the 95-percentiles of the number of taxa and number of site patterns (i.e., the numbers of unique sites in an MSA) distributions respectively, that is, datasets with a number of taxa < 470 and a number of patterns < 18, 764. Thereby, we omitted computations on excessively large datasets, which can be considered as outliers. 3. We initially sorted the datasets by their site patterns per taxa ratio and subsequently divided them into 20, 000 buckets. From every bucket, we randomly selected one dataset as a representative of that bucket. 4. For every representative, we simulated one MSAs using the corresponding RGS functionality, which we describe in more detail below.

Under these criteria, we selected 20, 000 DNA datasets representing more than 80% of overall RG data in terms of signal strength (approximated by patterns/taxa ratio) and simulated MSAs based on the inferred ML trees and their respective inferred model parameters.

During the simulation process, we encountered the problem of simulating MSAs with gaps. To simulate MSA gaps in a realistic way, the RG entry ’OVERALL_GAPS’ that reflects the number of gaps in the original empirical MSA is not sufficient. According to Haag et al. ( 18 ), the fraction of gaps in an MSA, does not have a substantial impact on dataset difficulty. Thus, one possibility is to avoid gaps in simulations and to solely execute the tests on MSAs without gaps. However, this induces another problem. While the mere gap fraction does not constitute a reliable difficulty predictor, the specific occurrence pattern of simulated gaps via an appropriate simulated insertion/deletion process affects the number of distinct site patterns in the MSA. In turn, this affects the resulting difficulty, as the difficulty is correlated with the patterns/taxa ratio (see Figure 1). In our initial experiments, we observed differences exceeding 20% between the number of site patterns in the simulated MSAs without an insertion/deletion process and the number of patterns in the respective RG entries. We used AliSim ( 24 ) to simulate MSAs for our experiments. For simulating gaps, AliSim offers the --indel and --indel-size parameters. The --indel parameter is a user-specified rate for generating insertions and deletions (so-called indels) at sites during the sequence simulation process. The indel lengths are drawn from distributions, which can be specified using the --indel-size parameter. Currently, AliSim supports the Geometric, Negative Binomial (NB), Zipfian, and Lavalette distributions originally proposed in ( 26 ), which can be separately set and selected for the insertion and deletion process. Different distributions also need to Source DNA RG here: RG (sim) IQ-TREE ( 11 ) FastTree2 ( 10 ) (emp) FastTree2 ( 10 ) (sim) RAxML-NG ( 13 ) (emp) RAxML vs FastTree ( 14 ) (sim) RAxML vs FastTree ( 14 ) (emp) Selection criteria < 470 taxa, < 18, 764 patterns 200-800 taxa, #sites / #taxa ≤ 4, #gaps ≤ 70% 500-239, 882 taxa, 65-1, 287 sites 250-5, 000 taxa 36-1, 879 taxa, 18, 328-37, 350, 521 sites ≥ 1, 000 taxa 117-27, 643 taxa be parameterized differently (e.g., mean and variance for indel sizes for the Negative Binomial distribution). Thus, even if we use the same distribution for insertions and deletions, we need to set a total of 6 parameters, that is, the insertion and deletion rates as well as the respective means and variances of the insertion and deletion sizes to parameterize the NB distribution. An additional complication is that the sequence length parameter used by some simulators (e.g., AliSim) only defines the sequence length of the single starting sequence on which the simulation is then carried out along the tree. Thus, when gaps are introduced during the simulation process, the overall MSA length will typically exceed the seed/root sequence length.

Under the assumption, that a simulation under the ‘correct’ settings for these 6 parameters will result in a simulated MSA where the differences between the number of sites, number of patterns, and gap percentage of RG entries and simulated MSAs are minimal, we can define an evaluation function d to quantify how well the simulated MSA matches the original MSA. Let sorig, porig, gorig and ssim, psim, gsim be the number of sites, number of patterns, and fraction of gaps in the RG entries and the simulated MSAs, respectively. Further, let w1, w2, w3 be arbitrary, yet constant weights. We can then define d as follows: ( 1 ) Since it is unclear, how to set these parameters such as to minimize d for a given simulated dataset, we resorted to Bayesian optimization ( 27–29 ), using the Python library skopt ( 30 ), and implemented the Bayesian Optimized iNdel seeKer (BONK) in RGS. Given w1, w2, w3, a scaling factor cseq, insertion and deletion size distributions, and a maximum number of optimization rounds nopt, BONK uses the optimizer to iteratively explore the indel rates and size distributions to minimize d for the simulated MSA. For the sake of simplicity, we arbitrarily chose (preliminary tests indicated no apparent reason to choose one distribution over another) the NB distribution for both, insertion, and deletion sizes. We use cseq to define the sequence length range for the root sequence. As stated above, the length of the root sequence has to be chosen as a function of the gaps being inserted to minimize the difference between sorig and ssim. The length interval is defined as [sorig × (1 − csim), sorig × (1 + csim)]. Preliminary experiments have also shown that, apart from nopt, the minimization of d is strongly affected by the values of w1, w2, w3. These weights penalize some property differences more than others, and since the number of sites, number of patterns, and gaps are correlated, it appears that choosing these weights is not trivial.

Thus, we deployed a second optimizer (results not shown) to optimize the weights with respect to the average d of simsequence length (no gaps) patterns (no gaps) gaps (no gaps) sequence length (after opt) patterns (after opt) gaps (after opt) ulations based on random subsets of datasets from RG. We found that w1 = 9, w2 = 10, w3 = 1, cseq = 0.7 worked sufficiently well, and used these weights for all our simulations with BONK.

After implementing BONK, we discovered SpartaABC ( 31, 32 ), a tool which uses the Approximate Bayesian Computation method ( 33, 34 ) to estimate the indel rates and sizes of a given MSA. The authors use a far more sophisticated distance function than d, overall comprising 27 features. Since the original MSAs are not available in RG, we cannot compute their proposed distance function and include SpartaABC into RGS. Thus, our simulations were conducted with BONK. However, we followed the suggestions regarding SpartaABC (and references to empirical studies) and made adjustments to BONK: We set the explored indel intervals to [0.0, 0.05], switched to the Zipfian distribution for indel sizes, and set the explored interval of the corresponding a shape parameter to [1.001, 2].

In addition to the aforementioned procedures, we also use the presence/absence matrices for partitioned datasets, when available in RG. A presence/absence matrix is a 2dimensional binary matrix denoting the presence of information (sequence data) for taxon k at partition l. If there is no sequence data for taxon k at partition l we set M [k, l] := 0, and M [k, l] := 1 otherwise. After the simulation process, we therefore remove per-partition sequences from the MSA according to M . This is conducted before calculating d, such that the optimizer is aware of the gaps introduced by M . Simulation Quality. We evaluated our BONK method for gapaware simulations by comparing the absolute differences between the properties (i.e., number of sites, number of site patterns, and gap proportion) of the simulated MSA and the respective RG database entry. We show the normalized differences in Table 3. In the worst case, that is, when the optimizer is not able to find parameters improving upon trivial simulations, we obtain an MSA without gaps. On average, we are able to simulate MSAs with properties which are closer to those of the RG entries than under trivial simulations without an indel model in terms of pattern numbers. However, due to the weights in our distance function, the difference in the gap proportions remained roughly the same - with the difference, that we now often insert more gaps than necessary. We further compared the difficulty of the simulated and original MSAs using Pythia ( 18 ). For this, we randomly selected a subset of 989 TreeBASE DNA datasets (note that the original empirical MSAs are required for difficulty prediction, which ( 1 ) no gaps ( 2 ) AliSim mimick ( 3 ) BONK ( 4 ) (shortened) SpartaABC are available in TreeBASE but not in RG), using the same size-based dataset selection criteria as in the preceding experiments, inferred a RAxML-NG tree using 50 parsimony and 50 random starting trees under the GTR+Γ model and used the inferred best-found tree, and model, to conduct the following distinct MSA simulations: ( 1 ) AliSim simulation without any indels; ( 2 ) AliSim simulation with the so-called mimicking function, which uses a specifiable (original) MSA to first infer a tree, then simulate an MSA, and subsequently superimpose the gap pattern from the original MSA to the simulated one; ( 3 ) BONK; ( 4 ) the SpartaABC tool, albeit with modifications to its internal pipeline. In SpartaABC, we removed the realignment step and reduced numbers of default burn-in and optimization rounds (1, 000/10, 000 instead of 10, 000/100, 000 respectively) to reduce the overall computational time, as some of the MSAs required more than 24 hours for a single indel rate inference. After the analysis, we used the suggested indel parameters by SpartaABC to simulate 10 MSAs with AliSim and selected the MSA with the lowest distance (as defined in SpartaABC) to the original MSA. We are fully aware that we do not fully utilize the capabilities of the SpartaABC algorithm, but we include it here to provide an intuition of how this method could perform in terms of difficulty matching, since it might be a valuable addition to the RGS simulation procedure.

In Table 4, we observe the following: Overall, there are small absolute differences between ( 1 ), ( 3 ), and ( 4 ). One might argue that the differences of these three approaches fall within the margin of error of the difficulty prediction. Surprisingly, inserting no gaps at all yields the difficulty prediction that is closest to the original difficulty, albeit the number of patterns more closely matches the original when using ( 3 ) instead of ( 1 ), and the gap proportion errors of ( 1 ) and ( 3 ) with respect to the original are similar. ( 1 ) and ( 3 ) do not require additional information about the underlying MSA, other than the properties already available in RG. Method ( 3 ) requires nopt = 100 steps, which means that it is at least 100 times slower than the single MSA simulation conducted by method ( 1 ). Method ( 4 ) requires an a priori computation of the gap features for the SpartaABC distance function. This could potentially be included into a future RG release, as it requires little additional memory to store these data and does not expose the original MSA to the public (i.e., it can directly be computed on the respective web-servers before being stored in RG). Method ( 2 ) is probably the most intuitive approach. Here, the distance d is negligible on average (data not shown). This method performs worst in terms of difficulty differences. We believe that different sites can have different importance for the inference (this seems at least to be true for empirical data ( 35 )), which a gap insertion by superimposing gaps does not take into account. Method ( 2 ) also requires the complete original MSA (if using the default AliSim mimick function) or at least the gap matrix of the original MSA to be available (if using the already inferred best tree and model present in RG).

Empirical Data. We conducted experiments on 5, 000 MSAs selected from TreeBASE ( 21 ). The dataset filtering criteria and the subsequent selection of representative datasets are analogous to the criteria for selecting simulated data from RG: We selected datasets from the 95-percentiles of number of taxa and site patterns, that is, datasets with a number of taxa < 214 and a number of site patterns < 3, 475. Multiple datasets contained sequences which were completely empty/undetermined. As these sequences cannot be placed onto a tree in any meaningful way, and the ML inference tools handle completely undetermined sequences differently (e.g., IQ-TREE2 terminates instantly, RAxML-NG does not, but inserts such sequences at random), we excluded datasets containing completely undetermined sequences from the analyses. During the experiments, 4 datasets containing ambiguous character encodings were excluded, since IQ-TREE2 incorrectly attempts to analyze them as AA datasets (if our interpretation of the log file is correct) and terminates with an error (this issue has already been reported by a user in the IQ-TREE Google Group). Furthermore, 2 datasets were excluded as they triggered assertions in RAxML-NG to fail during the tree search process. Overall, more than 75% of the datasets diversity (in terms of pattern/taxa ratio) in TreeBASE were covered.

Evaluation Pipeline. We implemented our experimental workflow via a Snakemake ( 20 ) pipeline. The datasets were selected, prepared, and generated as described above. The main steps of the pipeline are the following (also, see Figure 2): 1. Query the database for datasets, sort the results based on the number of taxa by the number of site patterns ratio, and split the results into n buckets 2. Randomly pick a dataset from each bucket 3. Run the FastTree2, IQ-TREE2, and RAxML-NG inferences under the General Time Reversible (GTR) substitution model in combination with the Γ model of rate heterogeneity. To obtain parsimony trees, select the parsimony trees with the highest log likelihood (LnL) score from the set of inferred RAxML-NG parsimony starting trees. Then, use the tree evaluation function of RAxML-NG (which optimizes branch lengths and model parameters while keeping the tree fixed) on all inferred trees to obtain a consistent and comparable LnL score 4. Compute pairwise Robinson-Foulds distances ( 36 ) and LnL-differences between the true trees or best-known ML trees (on empirical data) and the inferred ML trees and run the topological significance tests of IQ-TREE2 on the complete set of trees (including the true tree) 5. Plot accuracy statistics For the tree inferences, we used the respective default parameters to reduce the complexity of the experimental setup: RAxML-NG raxml-ng --msa assembled_sequences.fasta --model GTR+G --prefix [prefix] --seed [seed] --threads 4 --force perf_threads For partitioned datasets we substituted "–model GTR+G" with "–model [partition_file]" as all original partition files were available for empirical datasets. For simulated datasets, RG contains per-partition substitution model and MSA length parameters, based on which RGS simulates the MSAs and generates the partition files.

IQ-TREE2 FastTree2 iqtree2 -s assembled_sequences.fasta -m GTR+G -nt 4 --prefix [prefix] For partitioned datasets, we appended "-p [partition_file]". FastTreeDbl -gtr -gamma

-nt assembled_sequences.fasta To the best of our knowledge, FastTree2 does not support partitioned analyses and therefore all analyses were conducted on unpartitioned MSAs. We could not evaluate some of the FastTree2 trees via RAxML-NG as they contained multifurcations. In that case, we resolved the multifurcations at random prior to LnL evaluation.