Efficient phylogenetic tree inference for massive taxonomic datasets: harnessing the power of a server to analyze 1 million taxa

GigaScience , 2024, 13 , 1–12 DOI: 10.1093/gigascience/giae055 Technical Note Efficient phylogenetic tree inference for massi v e taxonomic datasets: harnessing the power of a server to analyze 1 million taxa César Piñeiro 1 , * and Juan C. Pichel 2 , * 1 Information Retrie v al Lab, CITIC, Universidade da Coruña, A Coruña 15008, Spain 2 CiTIUS, Universidade de Santiago de Compostela, Santiago de Compostela 15782, Spain ∗Corr espondence addr ess. César Piñeir o, Information Retrie v al Lab, CITIC, Universidade da Coruña, A Coruña 15008, Spain. E-mail: cesar alfr edo.pineir [email protected] ; Juan C. Pichel, CiTIUS, Universidade de Santiago de Compostela, Santiago de Compostela 15782, Spain. E-mail: juancarlos.pic[email protected] Abstract Bac kgr ound: Ph ylogenies pla y a crucial r ole in biological r esear c h. Unfortunately, the sear c h for the optimal phylogenetic tree incurs significant computational costs, and most of the existing state-of-the-art tools cannot deal with extr emel y large datasets in r easona b le times. Results: In this w ork, w e introduce the new VeryFastTree code (version 4.0), which is able to construct a tree on 1 server using singlepr ecision arithmetic fr om a massi v e 1 million alignment dataset in only 36 hours, which is 3 times and 3.2 times faster than its previous v ersion and FastTr ee-2, r especti v el y. This new v ersion further boosts performance by parallelizing all tr ee trav ersal operations during the tree construction process, including subtree pruning and re gr afting moves. Additionally, it introduces significant new features such as support for new and compressed file formats, enhanced compatibility across a broader range of operating systems, and the inte gr ation of disk computing functionality. The latter feature is particularl y adv anta geous for users without access to high-end servers, as it allows them to manage very large datasets, albeit with an increase in computing time. Conclusions: Experimental r esults esta b lish Ver yFastTr ee as the fastest tool in the state-of-the-art for maximum likelihood ph ylogen y estimation. It is pub licl y av aila b le at https://github.com/citiususc/veryfasttr ee . In addition, Ver yFastTr ee is included as a package in Bioconda, MacPorts, and all Debian-based Linux distributions. Ke yw ords: phylogenetics, v er y large datasets, performance, parallelism Introduction Inferring e volutionary r elationships or phylogenies is a formidable challenge in computational biology. The growth of datasets from next-generation sequencing has made largescale phylogeny estimation crucial. Ho w ever, the computational complexity of inferring phylogenies and performing multiple sequence alignment (MSA) presents a significant obstacle. Established methods like maximum parsimony (MP), maximum likelihood (ML), and Bayesian a ppr oac hes ar e computationall y intensive due to the NP-hard optimization problems they tackle [ 1 ]. As the number of taxa increases, these methods face a common hurdle: an exponential increase in the number of possible trees to explore. The leading heuristics for ML tr ee estimation, suc h as RAxML [ 2 ] and IQ-TREE [ 3 ], employ diverse strategies to search for the tree that maximizes the likelihood score. Although they hav e made consider able performance impr ov ements to handle larger datasets, RAxML, for example, was unable to r eac h conv ergence on a 10,000-sequence dataset e v en after a week [ 4 ]. Note that these tools ar e primaril y optimized for datasets with a limited number of sequences but a significant number of sites (i.e., phylogenomics). Ther efor e, when working with datasets comprising a large number of sequences, users must opt for tools such as FastTree-2 [ 5 ] and VeryFastTree [ 6 ], which are very fast heuristics, but they do not make very substantial attempts to optimize the likelihood score, or they may explore divide-and-conquer strategies [ 7–9 ]. In particular, our tool VeryFastTree was a big step forw ar d in terms of performance, building a tree on a standard server from a large 330k alignment, 3.5 times faster than FastTree-2. Ho w e v er, ther e was still room to impr ov e its speed, scalability, and memory consumption and also to add new functionalities. In this w ork, w e introduce the latest VeryFastTree code, version 4.0, showcasing its potential adv anta ges and ne w featur es compared to both its previous version and FastTree-2. The earlier iteration of VeryFastTree achieved high performance by parallelizing the most time-consuming phase of constructing the tree, specificall y the near est-neighbor interc hanges (NNIs), in comparison to FastTree-2. Ho w ever, the new VeryFastTree further enhances performance by parallelizing all tree traversal operations, including, among others, the subtree pruning and regrafting (SPR) operations, whic h ar e especiall y r ele v ant in terms of computing time when dealing with massive datasets. After a thorough experimental e v aluation, the ne w v ersion pr ov es to be se v er al times faster than the pr e vious one and FastTree-2 on a variety of large datasets. At the same time, VeryFastTr ee-4 incor por ates significant ne w featur es, including support for ne w and compr essed file formats , impro ved compatibility with a wider range of operating systems, and the addition of disk computing functionality. The lat- Recei v ed: October 23, 2023. Revised: April 17, 2024. Accepted: July 11, 2024 ©The Author(s) 2024. Published by Oxford Uni v ersity Pr ess GigaScience. This is an Open Access article distributed under the terms of the Cr eati v e Commons Attribution License ( https://cr eati v ecommons.org/licenses/by/4.0/ ), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is pr operl y cited. Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025 2 | GigaScience , 2024, Vol. 13 ter is especially valuable for users who do not have access to highend servers, as it allows them to process massive datasets, albeit with an increase in computing time. Ne w Fea tures and Optimiza tions VeryFastTree-4 (VFT4) is a big step forw ar d with respect to our first v ersion intr oduced in [ 6 ], fr om now on VFT3, designed to further accelerate the inference of phylogenies for massive alignments. Building upon the strengths of its pr edecessor, VFT4 intr oduces a host of innov ativ e featur es and optimizations aimed at ac hie ving e v en gr eater speed and efficiency. While the cor e principles used by VFT4 remain consistent with its pr e vious v ersions, significant enhancements have been made to push the boundaries of speed. Below are some of the most noticeable impr ov ements, optimizations, and features incorporated into VFT4. Par alleliza tion str a tegy VFT3 ac hie v ed high performance by parallelizing the most timeconsuming phase in the construction of the tree with respect to FastTree-2 (FT2) [ 5 ], the NNIs. Building upon this progress, VFT4 continues this a ppr oac h by parallelizing all operations involving tr ee tr av ersal, including the SPR mov es, further enhancing performance . T hese computations ar e par allelized using 2 strategies: tree partitioning , which divides the tree into multiple nonoverlapping subtrees, and parallel tra verse , a parallel breadth-first traversal . Tree partitioning Once an initial phylogenetic tree is constructed or passed as input argument, it is necessary to split the tree into disjoint subtrees in order to perform computations in parallel using different threads. Note that VFT4 uses rooted trees . T he tree partitioning method is an algorithm that, given an initial tree, prunes nodes from the root to generate several independent subtrees, which are then assigned to different threads to work with them independently. Once all threads have finished their work, those nodes that wer e r emov ed in the pruning pr ocess ar e subsequentl y pr ocessed using just 1 thread. In this way, all nodes in the tree are updated. In the ideal scenario, the computation of a node does not have dependencies on other nodes or these dependencies are contained within the subtree assigned to the same thread. Ho w ever, this is not the case fr equentl y, so dependencies m ust be contr olled to pr e v ent computation err ors . For instance , NNIs and SPR operations perform topology-modifying actions on the tree . T hey requir e tr ee partitioning to ensur e that node exc hanges performed by different threads do not overlap. Tree partitioning must ensure that these operations can be executed in parallel within separate subtrees, allowing them to work independently without interfering with each other. Before getting into the details, it is necessary to understand how the tree is tr av ersed and how data dependencies between nodes are taken into account. We will use the tree of Fig. 1 as example. Each node in the tree is visited before its parents, which means that a depth-first postorder tr av ersal is used [ 5 ]. Let’s assume that the tree is split into 2 subtrees assigned to threads t 1 and t 2 , r espectiv el y. Since the r oot node of the tr ee is not assigned to a thread, it will be labeled as not assigned . Blue and green nodes ar e pr ocessed in par allel by thr eads t 1 and t 2 , while red nodes, whic h hav e dependencies on other thr eads, ar e pr ocessed sequentially after all subtrees have been processed. This ensures that all nodes in the tr ee ar e updated effectiv el y. The dependencies of a node are determined by its relationships with its parents and sib- Root Root Subtree 1 P C N Root Subtree 2 A B Thread t2 Thread t1 Figure 1: Example of tree partitioning and data dependencies with a penalty of 2. The gr a ph r epr esents onl y internal nodes (no leav es). Blue and green nodes are processed in parallel by threads t 1 and t 2 , r espectiv el y. Red nodes (labeled as not assigned ) are processed sequentially after processing all subtrees. lings. If these are not located within the same subtree, the node cannot be computed. This dependency can be quantified as the distance between a node and the farthest one r equir ed for computation. We refer to this distance as the penalty value . T his means that nodes located at le v els lo w er than the penalty v alue fr om the subtr ee’s r oot cannot be pr ocessed in par allel, as they hav e dependencies on nodes outside the subtree. In the example of Fig. 1 , the penalty is 2, so Pcannot be computed since it depends on the root (its grandparent), which lies outside the subtree. On the other hand, C, N, and all their c hildr en hav e their gr andparents within the subtree, allowing them to be computed without conflicts. NNIs and SPR oper ations, whic h ar e the most time-consuming in tree construction, are parallelized using the tree partitioning strategy. In particular, NNIs perform a node exchange with either its parent or its uncle, follo w ed b y a recalculation of the w eights of the involved nodes and their respective parents . T his process incurs a penalty of 2, as the parent or uncle will be 1 le v el abov e the exchanged node, while the new parent will be located at the next le v el. In the case of SPRs, the penalty is dynamic and constrained by the par ameter maxSPRLength , whic h limits the maximum distance a node can move during the regrafting process. Other operations such as computing SH-like supports, updating all br anc h lengths, and optimizing all br anc h lengths will also make use of tree partitioning. VFT4 implements an advanced partitioning algorithm when compared to its predecessor, VFT3. This upgraded algorithm significantl y impr ov es both speed and adaptability across various scenarios. Unlike the pr e vious v ersion, whic h limited partitioning to NNIs, the new algorithm introduces a more versatile approach, allowing its application to other operations, such as those previousl y mentioned. This ne w algorithm considers the minimization of the penalty parameter as an e v aluation criterion. Subsequentl y, the r esulting subtr ees gener ated by this pr ocess, also r eferr ed to as solutions , are assigned to threads with the aim of achieving a good load balance. Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025 Efficient phylogenetic tree inference for massive taxonomic datasets | 3 Partitioning method Our partitioning method is guided by an objective function, whose goal is 2-fold. First, it should balance the workload assigned to eac h thr ead (i.e., the number of internal nodes to be processed by each thread). To model the w orkload, w e define the weight of a node as the number of internal nodes beneath it plus itself. This way, the workload associated to process a subtree is the weight of its root node . T he second goal of the objective function should be to reduce to the minimum the amount of sequential work, which corresponds to processing the not assigned nodes. Note that according to the Amdahl’s law [ 10 ], a small percentage of sequential work harms noticeably the performance and scalability of a parallel application. Let T be the input tree to be partitioned, whic h onl y contains internal nodes. A partition P(T ) is a set of disjoint binary subtrees s i included in T . We must highlight that not all nodes in the tree should be included in a partition. Let | P(T ) | denote the number of subtrees in the partition and n the number of threads. Each subtree s i ∈ P(T ) is assigned to a particular thread in such a way that P j (T ) contains all subtrees in the partition P(T ) assigned to thread j. This w ay, w e can calculate the workload of the subtrees assigned to a particular thread jas wo rklo ad( P j ( T )) =   wo rklo ad(s i ) | s i ∈ P j (T )  (1) We e v aluate a partition P(T ) using the following objectiv e function: sequent ial _ workl oad max { wo rklo ad( P j ( T )) | 1 ≤j ≤n } (2) Note that sequent ial _ workl oad refers to the workload of T when it is processed sequentially. The goal of our partitioning method will be to find a partition P(T ) that maximizes the value of the objective function. The pseudocode of the algorithm to find the best tree partitioning is detailed in Algorithm 1. The objective function aligns with the concept of Speedup , defined as the ratio of the execution time of a task when using multiple threads compared to the execution time when using a single thread. In this way, a value of 1 indicates equivalent speed to sequential processing, while a value of n signifies optimal workload distribution among n threads. The tree partitioning process in VFT4 includes an initialization phase where the weight values for each node in the tr ee ar e computed. As pr e viousl y mentioned, these weights r epr esent the number of descendant nodes for each node, adjusted according to the penalty value . T hen, the initial solution is constructed using the direct child nodes of the tree’s root. In other w or ds, the initial solution has as many subtrees as the root node has children. Nodes in the solution are always sorted by weight for efficiency. Next, an iteration process begins. During each iteration, the current solution is evaluated using Equation 2 . To achieve this, it is necessary to assign subtrees to the different threads using a heuristic method, which will be explained later. Once the solution is e v aluated, a ne w one is cr eated in suc h a way that the subtr ee root node with the highest weight is replaced by its child nodes. The e v aluation r esults of the last Witer ations ar e stor ed to use them as stopping criterion. We refer to Was the tendency window size . Note that the algorithm continues as long as the current solution has nodes to split, there are threads without at least 1 subtree assigned, or the stopping criterion is not r eac hed. Assignment of subtrees to threads The e v aluation of the quality of a solution r equir es ma pping the subtrees of a partition to the considered threads in such a way that the workload is balanced among them. This is a challenge itself and resembles the k-partitioning problem [ 11 ]. The k-partitioning problem is defined as follows: given a set of items { I 1 , I 2 , ..., I n } where item I j is of weight w j > 0 , find a partition S 1 , S 2 , ..., S m of this set with | S i | = k such that the maximum weight of all subsets S i is minimal. Ther efor e, considering our case, items are subtrees and k is the number of threads . T he k-partitioning problem is NP- hard, so algorithms/heuristics to a ppr oximate the solution should be used. Since our pr ocedur e r equir es solving the k-partitioning pr oblem e v ery time a new partition is e v aluated, we hav e opted for a simple gr eedy a ppr oximation based on the first-fit decreasing algorithm whose computational cost is low. The method starts sorting the subtrees in descending order according to their weight. Then, the k first subtrees are distributed among the threads . T he r emaining subtr ees ar e iter ated and assigned to the thread that curr entl y has the minimum workload assigned, that is, the lo w est value of wo rklo ad( P j ( T )) . Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025 4 | GigaScience , 2024, Vol. 13 Figure 2: Comparison of iterations required for achieving near-optimal tree partitioning using our method versus evaluating all possible tree partitionings without a stopping criterion. The bottom figure provides a zoomed-in view of the top one. Experiments were conducted with 30 threads and a dataset consisting of 274,000 unique sequences. Stopping criterion The stopping criterion in the tree partitioning process determines when the algorithm should terminate its execution. This criterion is r equir ed because completing all iter ations, until no mor e nodes are left to split, can sometimes result in a higher computational cost than working on the tree without partitions . T her efor e, the stopping criterion ensures that the algorithm halts its execution when a near-optimal solution has been found. To verify if the stopping criterion is met, we assess the quality of the last Wpartitions using Equation 2 to determine whether the ov er all tr end in r ecent iterations indicates an improvement or deterioration in the solutions (see the StoppingCriterion function in Algorithm 1). In order to ac hie v e this, a counter compar es eac h pair of elements within the tendency list, incrementing it if the earlier element is less than or equal to the later element and decrementing it otherwise. Consequently, the function returns true if the balance is negative, indicating a worsening trend, and false otherwise, signaling that the stopping criterion has not yet been met due to impro ving solutions . By default, Wis set to 50, but users can modify this value through the interface parameter (see the Commands Interface section for details). To demonstrate the advantages of our proposal, Fig. 2 illustrates the comparison between the number of iterations required by our method to ac hie v e a near-optimal tree partitioning and the scenario where all possible tree partitionings are evaluated without any stopping criterion. These results were obtained using 30 threads and the Large dataset that contains 274,000 unique sequences (refer to Table 1 for dataset details). Since 30 threads were considered, according to Equation 2 , the ideal (maximum) quality value of a partitioning is 30. Our method only requires 512 iterations to r eac h a v alue of 29.4, while without the stopping criterion, it is necessary to e v aluate mor e than 200,000 partitionings (iterations). Parallel breadth-first tra ver sal Tree partitioning enables parallel processing of a single tree by assigning disjoint subtrees to each thread, allowing them to compute independently and without interference. Ho w ever, when an operation has no dependencies and solely modifies the node being computed, cr eating subtr ees is inefficient, as eac h node can be dir ectl y assigned to a thread. So, starting from the deepest nodes of the tr ee, eac h thr ead handles a portion of the nodes and computes them sim ultaneousl y, following a br eadth-first tr av ersal pattern, whic h involv es visiting all nodes at a gi ven de pth le v el befor e moving on to the next le v el. Once finished with that le v el, the upper le v els ar e distributed until the r oot of the tr ee is r eac hed. This a ppr oac h is specificall y utilized in the following operations: computing initial pr ofiles, r ecomputing pr ofiles, r ecomputing ML pr ofiles, computing tree length, and computing the likelihood for each site. It is worth noting that FT2 and VFT4 store profiles for the internal nodes in the tree instead of storing a distance matrix, which requires far less memory. Each profile includes a frequency vector for each position and the weighted av er a ge of its c hildr en’s pr ofiles . T hese profiles are used to compute the distances between internal nodes in the tree and also the total distance from a node to all other nodes (see [ 5 , 12 ] for a detailed description). New thread levels To provide users control over parallelism and adaptability to various usage scenarios, VFT4 introduces 5 distinct thread levels. These le v els, including the original le v els fr om VFT3 (now r eferr ed to as le v els 0, 1, and 2), allow users to finely adjust the degree of parallelization based on their specific requirements. At le v el 0, VFT4 employs the same par allelization str ategy as FT2, but with the addition of new parallel blocks to enhance performance. Le v el 1 intr oduces par allel bloc ks that r equir e additional memory for computations, enabling more efficient processing. Le v el 2 utilizes the tree partitioning method to accelerate ML NNI r ounds. Le v el 3 performs in parallel all the computations but SPRs. Last, le v el 4 le v er a ges the tr ee partitioning method to expedite also SPR steps, but it only pays off with datasets containing a high number of alignments. For this r eason, le v el 3 is the default option. Eac h le v el in VFT4 is incr emental with r espect to the pr e vious ones. Ho w e v er, it is important to note that computation at le v el 2 and above follows a different tree traverse order, which may result in differ ent tr ees with r espect to the sequential execution. Ne v ertheless, these results remain strictly correct. By incorporating this m ultile v el a ppr oac h, VFT4 offers users the flexibility to optimize their parallelization strategy according to factors such as dataset size , performance needs , and desired trade-offs. Disk computing As pr e viousl y mentioned, phylogenetic tr ee infer ence r equir es a significant amount of memory, which can become problematic when dealing with numerous sequences and limited resources (e.g., when considering a low-end server with reduced available RAM memory). To tackle this challenge, VFT4 introduces disk computing , a technique that utilizes the hard drive to supplement the memory r equir ements. While disk computing aids in handling large datasets, it does impact performance due to the slo w er access speed of hard disks compared to RAM. Ho w e v er, the adv antages of being able to process larger datasets outweigh the performance trade-off. We classified the dynamically allocated memory used in the tr ee construction pr ocess into 3 types: fixed memory, v ariable memory, and computation memory. In particular: r Fixed memory is allocated at the beginning of the pr ogr am execution once the properties of the input dataset (such as the number of sequences, length, etc.) ar e established and r emains r eserv ed until the pr ogr am’s completion. For instance, fixed memory is used to store sequences and weights in each profile. Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025 Efficient phylogenetic tree inference for massive taxonomic datasets | 5 Table 1: Characteristics of the datasets used in the experimental evaluation. Information obtained using BigSeqKit [ 15 ] Dataset label File name Format Type Sequences Unique sequences Length Large sel03n.masked FASTA AA 331,550 274,401 1,287 Very Large gg_12_10_aligned FASTA DNA 1,075,170 858,234 7,682 Ultra-Large 1-million-taxon-run1 NEXUS DNA 1,000,000 989,109 21,946 r Variable memory , on the other hand, is allocated to store values during specific stages of the computation. For instance, it is the memory used to store the fr equency v ector (conditional probabilities) whose size changes dynamically and cannot be pr ecisel y pr edicted. Note that due to significant v ariations in data sizes , o v er allocation of v ariable memory is often impr actical. r Computation memory refers to private memory required by threads to perform their tasks . T his includes stack memory, which handles function calls; local variables; and temporary data while the pr ogr am is running. Additionally, computation memory includes some temporary variables, which store intermediate results and facilitate complex calculations, as well as data structur es, whic h or ganize and mana ge data for efficient processing. For example, computation memory is used to tempor aril y stor e the ne w computed conditional probabilities for possible exchanges within an NNI iteration. It is important to emphasize that once the temporary values are utilized, they are subsequently released. While fixed and variable memory can be offloaded to the disk, computation memory must be always maintained in RAM. This is because the performance of the computation process is critical. In addition, e v en for huge datasets, computation memory remains significantly smaller in comparison to fixed and variable memories. On the other hand, using swap space as an alternative when a process exceeds RAM capacity has important limitations and dr awbac ks. First, the swa p partition has a fixed size, potentially leading to out-of-memory errors if exceeded. This is especially r ele v ant when dealing with very large datasets, as is our case. Ther efor e, if the dataset exceeds the size of the swa p, onl y an administr ator would hav e the authority to expand the swap partition to pr e v ent the job fr om being aborted. In addition, since all pr ogr ams/a pplications shar e the same swa p space, excessiv e sw apping b y 1 pr ogr am can adv ersel y affect others, leading to contention for disk input/output (I/O) resources and reduced responsiv eness. In contr ast, the mec hanism used by VFT4 in disk computing is based on memory-mapped files, which offers better I/O performance due to the optimized kernel-le v el cac hing mec hanisms and efficient block-based access. Note that when using memory ma pping, an y av ailable disk can be used to stor e data, avoiding the swap partition size limitation. Additionally, memoryma pped files ar e also contr olled by the oper ating system, similar to swap space, and are not manually managed. This ensures efficient handling of data movement between memory and the disk, with the operating system making decisions based on system resources and r equir ements. VFT4 introduces the use of disk computing with 2 parameters: r -disk-computing : it facilitates the transfer of fixed memory to disk. r -disk-dynamic-computing : it facilitates the transfer of variable memory to disk. Employing both parameters allows VFT4 to move data from memory to disk in such a way that very large datasets can be processed e v en on low-end serv ers with limited RAM. Ho w e v er, it is essential to consider that this feature may come at the cost of decreased computational performance. In this way, the performance of VFT4 is dir ectl y influenced by the amount of memory sent to the disk. Hence, it is advisable, whene v er possible, to compute by tr ansferring onl y fixed memory to the disk. This a ppr oac h strikes a delicate balance between memory consumption and computational efficiency, optimizing the ov er all performance of the application. Optimized memory consumption VFT3 has been implemented using C ++ , whic h intr oduces additional memory overhead compared to the C implementation of FT2. Mor eov er, its mor e efficient utilization of thr eads leads to an increase in memory consumption as the number of threads grows. This increase is primarily caused by the replication of data structur es, sync hr onization mec hanisms, stac k space for each thread, and caching effects, among other factors. VFT4, on the other hand, has been redesigned to minimize the use of objects in memory-intensive sections and to release memory as soon as it is no longer needed. Consequently, in sequential execution, the memory r equir ements of VFT4 ar e e v en lo w er than those of FT2. Additionall y, to r educe the ov erhead intr oduced by eac h thr ead, VFT4 optimizes the stor a ge of common temporary data, ensuring that it is stored only once . T his significantly minimizes the overhead associated with m ultithr eading, r esulting in mor e efficient parallel execution compared to VFT3. Other optimizations and functionalities r Support for new and compressed formats : VFT3 and FT2 are limited to supporting the FASTA and Phylip formats, both stored as plain text. In contrast, VFT4 has extended its support to include the widely used Nexus and FASTQ formats . T he Nexus format allows for storing sequences and the initial tree within a single file. Additionally, it is common for datasets downloaded from internet repositories to be compressed in formats such as .gz or .bz . Previously, it was required to manually decompress these files before using them. Ho w ever, VFT4 can dir ectl y r ead compr essed sequences in compatible formats, thanks to the integration of the Zlib and libBZ2 libraries . T his a ppr oac h sav es time by removing the need for manual file decompr ession, r esulting in impr ov ed performance and avoiding the additional effort of reading larger uncompressed files from the disk. r Broader compatibility : VFT4 is a versatile tool that supports Linux, Windows, and macOS, including Windows executables. It is also conv enientl y av ailable in the Bioconda pack- a ge r epository [ 13 ], making it easil y accessible to the bioinformatics comm unity. Furthermor e , for macOS users , it is also available as a MacPorts pac ka ge [ 14 ]. Finall y, for Linux users, VFT4 is included as a pac ka ge in all Debian Linux distribu- Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025 6 | GigaScience , 2024, Vol. 13 tions, simplifying its installation and integration into various computing en vironments . r Better compilation support : The compilation process has been enhanced by adding support for new compilers, such as clang , and incor por ating ne w featur es like AVX512 in Windows builds . T he code has been optimized to comply with the latest compiler standards . Furthermore , parallel compilation has been implemented to accelerate the overall compilation time. Performance Ev alua tion Next, we present experimental results that clearly demonstrate the superior performance of our tool, VFT4, compared to its predecessors, VFT3 and FT2. We have evaluated the following aspects: running time, memory consumption, CPU usage, and topological accuracy. For the e v aluation, we selected 3 v ery lar ge datasets with v arying number of taxa and alignment lengths. Please refer to Table 1 for specific details about these datasets. It is important to note that only unique sequences were considered when constructing the trees. To ensure test reproducibility, we have included information about the parameters used to build the trees for each dataset: r Large : -spr 4 -gamma r Very Large : -nt -gamma -gtr r Ultra-Large : -nt -gamma -gtr The arguments have the following meanings: r -nt : nucleotide alignment used as input. r -spr : sets the number of rounds of SPR moves (default value is 2). r -gtr : utilizes the generalized time-reversible (GTR) model of nucleotide substitution. r -gamma : after optimizing the tree with a fixed rate for each site (the CAT model), VFT (and FT2) will rescale the tree to optimize the Gamma20 likelihood [ 5 ]. Experiments were conducted using 1 server with two 32-core Intel Xeon Ice Lake 8352Y @2.2 GHz processors and 512 GB of RAM. This server is part of a cluster installed at CESGA (Galicia Supercomputing Center, Spain) [ 16 ] running Rocky Linux v8.4 (kernel v4.18.0). We have used in the performance comparison the following tools and v ersions: Fast-Tr ee v.2.1.11 [ 17 ], VeryFastTr ee v3.0 (our pr e vious v ersion), and VeryFastTr ee v4.0.3. Running times First, w e sho w in Fig. 3 the running times when building the trees using single precision. We must highlight that, even if single precision is selected in VFT4 and FT2, many intermediate calculations such as vector reductions are performed with double precision to minimize errors. In the case of the Large dataset, VFT4 outperforms FT2 and VFT3, ac hie ving speed impr ov ements of 2.6 times and 1.6 times r espectiv el y, r esulting in an execution time of just 3.5 hours. When inferring the phylogenetic tr ee fr om the Very Large dataset, VFT4 completes the task in 17.2 hours, whereas FT2 requires 57.7 hours and VFT3 takes 35.4 hours. Finally, dealing with the Ultra-Large dataset, VFT4 is able to build the tree in 35.8 hours (i.e ., 1.5 da ys). Note that the time r equir ed by FT2 and VFT3 incr eases noticeabl y to 4.8 and 4.5 da ys , r espectiv el y. In other w or ds, VFT4 is 3.2 times and 3 times faster than FT2 and VFT3. On the other hand, Fig. 4 presents the running times for tree construction but using double-precision arithmetic. In this case, for the Large dataset, VFT4 is 6.9 times and 2.4 times faster than FT2 and VFT3, r espectiv el y, r educing the execution time to less than 3 hours. Note that VFT4 builds the tree faster when using double precision than considering single precision (2.8 versus 3.5 hours). This is caused by the good behavior of the vectorization strategies used by VFT4 when considering protein alignments. On the other hand, inferring the phylogenetic tree from the Very Large dataset using VFT4 takes 18 hours, while FT2 r equir es 61.7 hours and VFT3 41.3 hours. As a result, VFT4 is again the fastest tool. In conclusion, when dealing with the Ultra-Large dataset, VFT4 can complete the tree-building process within 41.2 hours, or 1.7 da ys . Notably, the time required by FT2 and VFT3 increases significantly to 5.3 and 4.6 da ys , r espectiv el y. To put it simply, VFT outperforms FT2 and VFT3 in speed by 3.1 times and 2.7 times, r espectiv el y. As mentioned in the Introduction, in addition to FT2 and VFT, ther e ar e other state-of-the-art tools for ML tree estimation, with RAxML and IQ-TREE being the most commonly used by the scientific community. Ho w ever, they are limited in their ability to process datasets containing a large number of sequences due to their extensive running times. To validate the observations made in pr e vious w orks [ 4 ], w e assessed both tools, RAxML-NG v.1.2.0 and IQ-TREE v.2.1.3, using our Large dataset. Since RAxML-NG can be executed on a cluster, we conducted experiments with this tool using 4 computing nodes instead of just 1 server. We manually constrained the experiments to a maximum running time of 1 week. Both tools were unable to construct the tree within that time frame. It is worth noting that VFT4 can estimate the tree for the Large dataset in a ppr oximatel y 3 hours (see Figs. 3 and 4 ). This confirms that both RAxML and IQ-TREE are not well suited for datasets containing a very large number of sequences. Ther efor e, to the best of our knowledge, these results establish VFT4 as the fastest tool in the state-of-the-art for ML phylogeny estimation. But mor eov er, VFT4 allows the pr ocessing of massiv e datasets that would otherwise be intractable or would r equir e excessiv el y high computing times. Memory consumption and CPU usage In the pr e vious section, w e pointed out ho w the mor e effectiv e use of thr eads r esults in higher memory consumption as the number of threads increases. VFT4 addresses this issue by optimizing the stor a ge of commonly used temporary data, ensuring it is stored only once. As a result, the overhead associated with multithreading is significantly reduced, leading to a more efficient parallel execution compared to VFT3. Figure 5 displays the CPU usage and memory consumption for all the considered tools when building the tree from the Ultra-Large dataset using single-precision arithmetic. The purpose of this illustration is to showcase the memory optimizations implemented in VFT4. Each time step in the gr a phs corr esponds to 3,600 seconds. It can be observed that the maximum memory consumed by VFT4 is a bit higher with respect to FT2, 272 GB and 228 GB, respectiv el y. This incr ease is attributed to a mor e efficient utilization of threads. It is important to note that VFT4 extensiv el y utilizes parallelism, employing the maximum available number of threads (64) for most of the time. In contrast, FT2 is constrained to using only a few threads, as indicated by the blue lines in the figures. As a consequence, VFT4 r equir es some extra memory for multithreading, but it comes with the benefit of a remarkable decrease in the running time. On the other hand, VFT4 significantly reduces Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025 Efficient phylogenetic tree inference for massive taxonomic datasets | 7 Large Very-Large Ultra-Large 1 5 10 50 100 500 Time (hours) - log scale FastTree-2 VeryFastTree-3 VeryFastTree-4 Figure 3: Running times of FT 2 , VFT 3 , and VFT 4 for building the trees using single-precision arithmetic and different datasets. Large Very-Large Ultra-Large 1 5 10 50 100 500 Time (hours) - log scale FastTree-2 VeryFastTree-3 VeryFastTree-4 Figure 4: Running times of FT 2 , VFT 3 , and VFT 4 for building the trees using double-precision arithmetic and different datasets. CPU Memory (GB) A Memory (GB) B Memory (GB) C Figure 5: CPU usage (# of cores) and memory consumption when building the single-precision trees using as input the Ultra-Large dataset. Time step = 3,600 seconds. Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025 8 | GigaScience , 2024, Vol. 13 memory consumption compared to VFT3. In particular, VFT3 requires twice the maximum memory used by VFT4. Furthermore, it also demonstrates greater efficiency in terms of parallelism. As we pointed out pr e viousl y, VFT4, like FT2, uses sequence profiles of internal nodes of the tree to implement neighbor joining instead of storing a distance matrix. Typically, as stated in [ 12 ], a distance matrix for an alignment with Nsequences, L sites, and a differ ent c har acters in the alphabet r equir es O (N 2 ) space. How- e v er, the upper bound of memory space r equir ed by VFT and FT2 is O (NLa ) . Eac h pr ofile comprises a fr equency v ector for eac h position and the weighted av er a ge of its c hildr en’s pr ofiles. If we take a look in detail, VFT4 r equir es 272 GB of memory to process the Ultra-Large dataset (see Fig. 5 C). Its corresponding par ameters ar e N = 989 , 109 , L = 21 , 946 , and a = 4 (nucleotides). As a consequence, taking into account that the r esulting tr ee contains 2 Nnodes and single precision uses 4 bytes, the memory required by the profiles is the sum of the weights and frequencies for each node in the tree: r Weights r equir e exactl y 2 N ×L ×4 bytes = (2 ×989 , 109) × 21 , 946 ×4 bytes = 173 . 65 GB. r The memory of frequencies can be estimated as 2 N ×a ×g × 4 bytes, with “g” being the number of gaps in the sequence. Since “g” is not a constant value and changes in each node as the topology of the tree evolves, an exact value for the memory space used cannot be pr ovided. In an y case, most of the difference between the total memory used by VFT4 (272 GB) and the memory r equir ed onl y b y the w eights of the profiles (173.65 GB) is due to the frequencies. Disk computing In man y cases, r esearc hers ar e faced with limited computing resources , often ha ving access onl y to low-end serv ers that possess a limited amount of memory. To deal with this issue, as was commented pr e viousl y, we intr oduced disk computing in VFT4, a ne w feature that allows to offload static and dynamic memory to the disk with the aim of handling very large datasets even on small servers . T he obvious drawback is an increase in the running times. An example of the effects of disk computing can be found in Fig. 6 . Each time step in the gr a phs corr esponds to 300 seconds. In particular, Figs. 6 A and 6 B show the CPU usage and the memory footprint in a normal execution of VFT4 when processing the Large dataset considering double-precision arithmetic and 1 and 64 thr eads, r espectiv el y. In these cases, the maximum memory consumed was 58.5 GB (1 thread) and 84.42 GB (64 threads). Note that the spike in the memory consumption is caused by the computation and stor a ge of the profiles at each node of the tree. On the other hand, Figs. 6 C and 6 D also display the CPU usage and the memory footprint but using disk computing with 1 and 64 threads, as well as limiting the memory of the server to just 16 GB. It can be observed that thanks to this new feature, we can successfully build the tree on a small server. Considering 64 threads, for example, the maximum memory used is approximately 5.3 times less compared to a normal execution. It may seem a ppar ent, but without disk computing, the processing of datasets that surpass the available memory would be impossible. It would result in an out of memory err or, whic h would be the case of using both FT2 and VFT3. On the other hand, the VFT4 running time increases, using 1 thr ead fr om 16.8 to 27.5 hours while using 64 thr eads fr om 2.8 to 53.4 hours. In other w or ds, using disk computing on a 16 GB server is from 1.6 times (1 thread) to 19 times (64 threads) slo w er than the standard execution. An interesting observation is that the sequential execution with disk computing is faster than the corresponding execution with 64 threads . T his is due to the limitations of disk and I/O bandwidth, which struggle to handle the extr emel y high number of requests generated when using 64 threads, leading to contention on the bus . T her efor e, we r ecommend r educing the number of threads when utilizing the disk computing feature. Topological accuracy In the liter atur e, se v er al studies hav e examined the performance of leading ML tools, such as RaxML, IQ-TREE, and FT2, regarding accuracy and computational efficiency with large datasets. Liu et al. [ 18 ] compared RAxML and FT2 using simulated datasets for single genes containing 1,000 or more sequences across various estimated alignments . T hey concluded that both methods exhibited similar topological accuracy. Lees et al. [ 19 ], in their exploration of ML heuristics through simulations, found RAxML and IQ-TREE to be comparable in accuracy, both outperforming FT2. Mor eov er, 2 studies [ 20 , 21 ] e v aluated FT2 and RAxML on datasets containing fr a gmentary sequences . T hey observed that FT2 performed less accur atel y than RAxML when handling alignments with a high proportion of fragmentary sequences. Overall, these studies have shown that IQ-TREE and RAxML are both very good ML heuristics with respect to ML scores, while all of the studies have shown that FT2 is indeed very fast, but it is so good at ML score optimization. On the other hand, as we explained in our pr e vious study [ 6 ], VFT is a highly tuned implementation of FT2 that k ee ps the phases , methods , and heuristics used by FT2 to estimate the phylogenetic tree. For this reason, it is only necessary to compare the topological accuracies with respect to FT2 (and our pr e vious version of VFT), since the works commented above already performed a thorough comparison between RAxML, IQ-TREE, and FT2. Similar to [ 5 ], we defined the topological accuracy as the proportion of splits in the true trees that are successfully recovered by eac h r espectiv e tool. T his metric is the in verse of the topological Robinson–Foulds distance [ 22 ], normalized to a range between 0 and 1. Our findings demonstrate that both VFT3 and VFT4 exhibit determinism while maintaining the same le v el of accuracy as FT2. To validate this, we assessed topological accuracy using 5,000- sequences sim ulated pr otein alignments [ 17 ], whic h wer e also emplo y ed in the original FT2 study for the same purpose . T herefore, the true phylogeny is a vailable . All simulations are based on trees from profile alignments of biological sequences and include variable rates across sites . T he simulations are described in more detail in [ 23 ]. Trees were computed using 64 threads, doubleprecision arithmetic, and -gamma -spr 4 parameters. All tools used the same seed for initializing the random number generator. The experimentation script, treecmp.py , is accessible in our repository. We show results for VFT4 considering 2 different levels of parallelism (i.e., thread levels ). As we explained pr e viousl y, le v el 3 performs in parallel all the computations but SPRs, while le v el 4 le v er a ges the tr ee partitioning method to accelerate also SPR mo vements . For VFT3 and VFT4, topological accuracy values were computed by av er a ging 7 measurements (1 for each alignment in the dataset). Conv ersel y, par allel executions of FT2 did not yield deterministic r esults. Consequentl y, we pr esent a r ange of accuracy values . T his range was derived by averaging the minimum and maximum values obtained from 10 executions for each alignment. Accur acy r esults expr essed in percenta ges ar e shown in Table 2 . Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025 Efficient phylogenetic tree inference for massive taxonomic datasets | 9 Time Ste p A Time Ste p B Time Ste p C Time Ste p D Figure 6: Effects of using disk computing when building the double-precision trees using as input the Large dataset. Time step = 300 seconds. Table 2: Topological accuracies (in %) obtained by FT 2 , VFT 3 and VFT 4 using 5,000-sequences simulated protein alignments. Note that the parallel version of FT2 is nondeterministic, so minimum and maximum values are displayed between br ac kets Tool COG438 COG583 COG596 COG642 COG1028 COG1309 COG2814 Average VFT4 (level 3) 85.35 79.11 88.14 82.66 86.19 86.55 79.59 83.94 VFT4 (level 4) 85.43 79.13 87.72 82.19 86.37 87.29 79.92 84.01 VFT3 85.51 79.57 87.82 81.89 86.29 87.09 80.18 84.05 FT2 [85, 85.83] [79.21, 79.91] [88.04, 88.78] [82.11, 82.95] [86.30, 87.03] [86.19, 87.59] [79.06, 80.66] [83.70, 84.68] Based on the results obtained, we conclude that VFT4 produces trees with an accuracy level within the same range as FT2. Additionall y, minor differ ences wer e observ ed when using differ ent thr ead le v els. In particular, the most a ggr essiv e one in terms of par allelism, le v el 4, shows a slightly better behavior than level 3. Commands Interface The new VeryFastTree version 4.0 was designed with extensive cross-platform compatibility, offering support for a variety of operating systems, including Linux, Windows, and macOS. This broad compatibility ensures that users can access and utilize the tool seamlessly across different computing en vironments . Windows users, in particular, benefit from the availability of executable versions of VFT4, simplifying the installation process and widening its user base. Additionally, for those in the bioinformatics field, VFT4 can be effortlessly located within the Bioconda pac ka ge r epository, str eamlining both installation and integr ation into bioinformatics w orkflo ws . VFT4 is also a vailable for macOS users as a MacPorts pac ka ge. Finall y, Linux users will a ppr eciate that VFT4 is included as a pac ka ge in all Debian Linux distributions, making installation straightforw ar d and facilitating its integration into diverse computing setups. Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giae055/7730000 by guest on 14 January 2025