Megascale microbiome analysis with DartUniFrac

UniFrac^1,2,3 is a phylogenetic beta-diversity metric that has been widely used in many microbiome and/or microbial ecology studies (>15,000), including large-scale ones such as the Earth Microbiome Project (EMP)⁴ and American Gut Project (AGP)⁵, due to its capability to incorporate gene/genome evolutionary histories into community dissimilarity metrics. By leveraging branch-length information on the phylogeny, UniFrac frequently yields stronger between-group separation (for example, higher PERMANOVA R² values and clearer ordinations) than nonphylogenetic distances such as Bray–Curtis and Jaccard dissimilarity, particularly when community turnover involves distantly related lineages, as shown in comparative evaluations across diverse datasets². As a fundamental biological community distance metric, UniFrac is not limited to amplicon sequencing-based community profiling techniques, and it is natural to extend it to metagenomic/genomic-based community profiling techniques as long as the unit of interest is clearly defined (for example, amplicon sequence variants/operational taxonomic units (OTUs) or a metagenome-assembled genome-derived species unit based on average nucleotide identity, or ANI)⁶. The computational complexity of UniFrac (weighted) is proportional to the number of taxa in the phylogenetic tree and is quadratic in the number of samples in a given study (pairwise comparisons)⁷. Model-based estimates of microbial diversity using the 16S rRNA gene suggest that Earth may harbor more than 10¹² microbial species⁸ or even more with whole-genome-based approaches. However, we are far from capturing these diverse microbial species. As high-throughput sequencing is becoming readily accessible in standard microbiome studies, we can now sequence thousands of samples with millions of species in a single project. This, however, has created a major challenge—computing all-versus-all UniFrac distances becomes a bottleneck for real-world datasets with even more samples and taxa. Many computational optimizations have been developed for faster UniFrac computation for large numbers of samples (for example, a few thousand) over the past 20 years^7,9,10, for example, faster tree traversal strategies¹⁰, better parallel efficiency (Striped UniFrac)⁷ and hardware acceleration (for example, single instruction, multiple data (SIMD) and graphics processing unit (GPU))⁹. However, these improvements were based on the same UniFrac algorithm and cannot scale further.

Here, we introduce a UniFrac algorithm, DartUniFrac, that scales to millions of samples with billions of taxa. First, we revisited the original UniFrac algorithm¹ (unweighted; Supplementary Fig. 1a): ({D}_{{{mathrm{UniFrac}}}}(A,B)=frac{sum _{iin E}{ell}_{i}times |{X}_{i}(A)-{X}_{i}(B)|}{sum _{iin E}{ell}_{i}{times} max {{X}_{i}(A),{X}_{i}(B)}}), where ({X}_{i}(S)!:=mathop{rm{max}}limits_{jin text{Desc}(i)}{x}_{!j}(S)in {0,1}) and x_j(S) is defined as ({x}_{!j}(S)=left{begin{array}{ll}1, & {text{if taxon}},{j};{text{is present in sample}},{S},\ 0, & text{otherwise}.end{array}right.). Weighted UniFrac is defined as ({D}_{{{textrm{WUniFrac}}}}=frac{{sum }_{iin E}{{ell}}_{i}left|{sum }_{jin {{textrm{Desc}}}(i)}{a}_{j}-{sum }_{jin {{textrm{Desc}}}(i)}{b}_{!j}right|}{{sum }_{iin E}{{ell}}_{i}({sum }_{jin {{textrm{Desc}}}(i)}{a}_{j}+{sum }_{jin {{textrm{Desc}}}(i)}{b}_{!j})}) (see the detailed definition of each individual variable in Supplementary Fig. 1b and Supplementary Methods, a different but equivalent equation to the original weighted UniFrac)². We prove that unweighted and weighted UniFrac are essentially weighted Jaccard similarity on tree branches (Supplementary Methods and Supplementary Fig. 1a,b). For unweighted UniFrac, ({D}_{{{textrm{UniFrac}}}}=1-{J}_{{textrm{w}}}(mathbf{x},mathbf{y})), where ({J}_{{textrm{w}}}(mathbf{x},mathbf{y})=frac{sum min ({x}_{i},{y}_{i})}{sum max ({x}_{i},{y}_{i})}) (also known as weighted Jaccard similarity), ({x}_{i}={{ell}}_{i}times {X}_{i}(A)) and ({y}_{i}={{ell}}_{i} times {X}_{i}(B)) (Supplementary Fig. 1a). For weighted UniFrac, ({D}_{{{textrm{WUniFrac}}}}=frac{1-{J}_{{textrm{w}}}(mathbf{x},mathbf{y})}{1+{J}_{{textrm{w}}}(mathbf{x},mathbf{y})}), where ({J}_{{textrm{w}}}(mathbf{x},mathbf{y})=frac{sum min {{x}_{i},{y}_{i}}}{sum max {{x}_{i},{y}_{i}}}), ({x}_{i}=)({ell }_{i}times {A}_{i},{y}_{i}=)({ell }_{i} times {B}_{i}) and ({A}_{i}=sum _{jin {{textrm{Desc}}}(i)}{a}_{!j},{B}_{i}=sum _{jin {{textrm{Desc}}}(i)}{b}_{!j}) (Supplementary Fig. 1b). In both unweighted and weighted UniFrac, computing weighted Jaccard similarity ({J}_{{textrm{w}}}(x,y)=frac{sum min {{x}_{i},{y}_{i}}}{sum max {{x}_{i},{y}_{i}}}) is the limiting step because, for real-world datasets, the dimension of x aninid y (equal to the number of branches in the phylogenetic tree) can be millions or even billions. The problem becomes how to compute weighted Jaccard similarity quickly and efficiently. First, we used a succinct balanced-parentheses data structure to represent phylogenetic trees, enabling efficient, pointer-free navigation operations such as parent, child and sibling traversal through auxiliary rank/select and parenthesis-matching indexes (Fig. 1a and Methods). This allows representation and efficient traversal of trees with billions of taxa. Each branch of the tree is now a set element, but the weight of each set element is related to the samples, either the presence or absence of taxa that are descendants of the branch (unweighted) or the sum of the relative abundances of taxa that are descendants of the branch (weighted) (Fig. 1a, Supplementary Fig. 1a,b and Supplementary Methods). After obtaining the weighted sets (Fig. 1a), computing weighted Jaccard similarity is the most expensive step because, for real-world microbiome datasets, there can be billions of taxa/branches and millions of samples. We relied on weighted MinHash (also known as sketching algorithms) to sketch weighted sets for computing weighted Jaccard similarity. MinHash belongs to a category of locality-sensitive hashing algorithms that is widely used in data mining for large-scale web and text comparison. MinHash is also widely used in genomics and metagenomics for computing ANI at large scale^11,12,13,14, and it allows for the computation of weighted Jaccard similarity efficiently with controllable estimation error, which converges to 0 as the sketch size increases (Methods). Among several weighted MinHash algorithms, DartMinHash and efficient rejection sampling (ERS) are the most efficient for sparse and dense sets, respectively^15,16. In the resulting sketch (a low-dimensional representation with length S of original weighted sets from all branches), computing an integer-based Hamming similarity (normalized) for each pair of sketches estimates weighted Jaccard (J_w) and thus also unweighted UniFrac or weighted UniFrac (Fig. 1a and Methods). Real-world microbiome datasets are sparse, such that most branch sets are empty, making the weighted MinHash sketch step extremely fast, with DartMinHash being the fastest weighted MinHash algorithm both in theory and in practice¹⁵, to the best of our knowledge. Overall, the big-O notation for DartUniFrac is O(N × T_active + N × S × log(S) + N² × S), where N is the number of samples, T_active is the number of active branches on average across all samples and S is the weighted MinHash sketch vector length (see detailed big-O analysis in Methods). Big-O notation for the exact UniFrac is O(N² ×T), where T is the total number of branches or taxa and, in practice, S << T (S is less than or equal to 2,048 in practice). The computational bottleneck for DartUniFrac for millions of samples is dominated by the pairwise integer Hamming similarity step (O(N²× S), which computes equal slots for sketch vectors with the same length, where the slots are 16-bit integers) (Fig. 1a and Methods). We first maximized central processing unit (CPU) performance for this step using multithreading and SIMD. However, Hamming similarity computation is bound by memory bandwidth, such that further speedups on a CPU are limited (for example, more CPU threads), but can benefit from hardware with higher memory bandwidth. We therefore offloaded this step to GPUs to further speed up the computation for large numbers of samples (Fig. 1a). A streaming mode (both CPU and GPU implementation), which allows block-by-block computation of the full distance matrix, is also provided when the distance matrix cannot fit into random access memory (RAM) (Fig. 1a).

DartUniFrac-CPU scales well with increasing number of CPU threads (Fig. 1b) and is more than 200 times faster than the state-of-the-art UniFrac algorithms (Striped UniFrac⁷ and a hardware-optimized version, called unifrac-binaries⁹; both are exact UniFrac algorithms) for in-memory mode (Fig. 1c and Table 1). Additionally, DartUniFrac-CPU can compute pairwise UniFrac for a million samples (87,522 taxa) in 1.8 h on a CPU (Supplementary Table 1), provided that enough RAM is available or in approximately 4.5 h in a streaming mode when there is not enough RAM (Supplementary Table 1). Notably, the existing state-of-the-art unifrac-binaries-CPU required more than 20 days to complete the same task on a CPU. In streaming mode, memory requirements can be one order of magnitude smaller for DartUniFrac-CPU (Supplementary Fig. 2). DartUniFrac-GPU is on average ~900 times faster than the unifrac-binaries-GPU implementation (Fig. 1d) while consuming ~24 times less GPU memory for ~87,522 taxa and 500,000 samples (Fig. 1e). Procrustes analysis showed that DartUniFrac aligned almost perfectly with the exact UniFrac for the Global Water Microbiome Consortium (GWMC) dataset (unweighted in Fig. 1f and weighted in Supplementary Fig. 3) (M² = 0.0042 and 0.0039; P < 0.001) and the EMP dataset (Supplementary Fig. 4). Mantel analysis showed strong correlations between DartUniFrac and exact UniFrac distances across all datasets (Mantel r ≥ 0.98; P < 0.001; Fig. 1g and Supplementary Table 2). To better understand the effects of sketch size (the length of vectors storing hashes from DartMinHash and ERS) on DartUniFrac accuracy and running time, we varied the sketch size and conducted a series of experiments on several accuracy metrics for the GWMC dataset. Root-mean-square error (RMSE) converges to 0 as sketch size increases (Fig. 1h and Supplementary Fig. 5a). Mantel correlation coefficients (both Pearson and Spearman), Procrustes M², PERMANOVA effect size (pseudo-F) and R² and nearest-neighbor recall all converge to the ground truth as sketch size increases based on the DartMinHash estimator (Fig. 1i and Supplementary Figs. 5b–e and 6a).

In tests of whether DartUniFrac runtime is sensitive to tree size and number of samples, we showed that, for a fixed number of samples and sparsity, the running time to extract branches via balanced-parentheses representation linearly increases with respect to the number of taxa in the tree and is generally less than a few seconds (Fig. 2a), consistent with the claim above that balanced parentheses can scale to billions of taxa (Supplementary Methods). As sparsity increases without changing the number of samples, running time increases proportionally (Fig. 2a). In practice, increasing the number of samples will increase both the number of taxa and the number of total branches. Therefore, DartMinHash sketching time increases proportionally with the number of samples, but this step is generally fast with multithreaded implementation (Fig. 2b). DartUniFrac-GPU is approximately two times slower in streaming mode than in in-memory mode but requires ~10 times less CPU RAM (Fig. 2c,d). Principal-coordinate analysis (PCoA) on the EMP and GWMC datasets¹⁷ showed that DartUniFrac with the default sketch size provides results nearly identical to exact UniFrac distances on the exact same input data, as measured by Mantel and Procrustes analysis (see Fig. 2e,f for DartUniFrac, unweighted and weighted EMP, and Fig. 2g,h for exact UniFrac, unweighted and weighted, for the EMP dataset; in Supplementary Fig. 7a–d, panels a and c show DartUniFrac and b and d show exact UniFrac for the GWMC dataset). PCoA analysis based on Qiita amplicon studies (279,443 samples)¹⁸ showed similar results (Supplementary Fig. 8a–d, with panels a and c showing unweighted results and b and d showing weighted results). Similar results were also found for large-scale animal gut microbiome datasets and human gut metagenomic datasets (Supplementary Figs. 9a–d and 10a–d). GPU support provides a speedup of another ~20 fold relative to CPU implementation (Fig. 2i and Supplementary Table 3). Integer Hamming similarity computation in DartUniFrac is bound by memory bandwidth, and we observed significantly faster runtimes on recent GPUs relative to CPUs (Fig. 2i and Supplementary Tables 3–6). Because DartUniFrac-GPU only stores only the sketch vector (16 bits per sketch element) for each sample in GPU memory, it can easily scale to millions of samples with moderate GPU memory (for example, 48 GB is enough for 10 million samples with the default sketch size) (Fig. 2j). DartUniFrac-GPU can be even more GPU memory efficient as the number of samples further increases (Methods). DartUniFrac allows computing of UniFrac distances beyond the limitations of BIOM Format (version 2.1.0)¹⁹, a common container for microbiome feature tables, which is capped at 2²⁰ nonzero values in the feature table; with feature tables that extend the limit for nonzero values to 2⁶⁴, DartUniFrac-GPU can finish in 13.8 min for 500,000 samples with 20 million taxa using two GPUs (see detailed hardware information in Supplementary Table 4). For 1,000,000 samples or more, DartUniFrac-GPU is >1,000 times faster than unifrac-binaries-GPU in terms of compute time because both require the same input–output time (Supplementary Table 5).

At the default sketch size of 2,048 for DartMinHash, the observed DartUniFrac RMSE and absolute error across simulated UniFrac distances closely matched the theoretical error (Fig. 2k). ERS also showed similar PCoA results for the EMP dataset (Supplementary Fig. 11a,b) and similar Mantel results with exact UniFrac (Supplementary Table 7) but was ~60 times slower than DartMinHash for the sparse EMP dataset (Supplementary Table 8 and Supplementary Fig. 12a). However, for a denser dataset, the GWMC dataset with an average sparsity of ~10.3% (Supplementary Fig. 12b), ERS was ~2–3 times faster than DartMinHash (Fig. 2l and Supplementary Table 9). For the ERS estimator, the results were consistent with those of the DartMinHash estimator, although ERS was slightly less accurate and was slower than DartMinHash for the same sketch size (Supplementary Figs. 6b, 13a–f and 14a,b). To evaluate the effect of sketch size on running time (both CPU and GPU implementation), we used the AGP and EMP datasets (~50,000 samples) and showed that running time increased linearly with respect to the sketch size (Supplementary Figs. 5f, 13f and 14a,b), dominated by the pairwise integer Hamming step. Therefore, as sketch size increases, DartUniFrac becomes more accurate but slower. The speed–accuracy trade-off is consistent with theoretical expectations¹⁵. Similar results were observed for the AGP dataset (Supplementary Figs. 15a–d, 16a–d, 17a–d and 18a–d). We also studied the effect of sketch size on hierarchical clustering topology to determine how much sketching resolution is needed for DartUniFrac to recover tree structures consistent with the truth. Hierarchical clustering trees from DartUniFrac increasingly matched the truth as sketch size grew: cophenetic correlations rose, and RMSE decreased for both the DartMinHash and ERS estimators, with convergence evident at larger sketch sizes (Supplementary Fig. 19a–d). Because hierarchical clustering is highly sensitive to small distance perturbations, lower accuracy at small sketch sizes such as 2,000–3,000 is expected. Topological agreement also improved with sketch size, as raw Robinson–Foulds distance decreased substantially from ~566 at a sketch size of 1,024 to ~337 at a sketch size of 6,656, indicating more shared internal clades/splits for DartMinHash and ERS (Supplementary Fig. 20a–d). Similar results were observed for cophenetic correlation-based analysis (Supplementary Fig. 21a–d) and Robinson–Foulds distance-based analysis (Supplementary Fig. 22a–d) for the AGP dataset.

To further speed up downstream analysis for a large number of samples (large distance matrix output), we developed a new fast PCoA (fPCoA) algorithm. The key step is to replace the exact but much slower singular value decomposition (SVD) step with a faster approximate SVD—randomized SVD—with controllable error. Among several variants of randomized SVD, we relied on the subspace iteration-type randomized SVD²¹ as subspace iteration was shown to have fewer vectors to orthogonalize (Methods). fPCoA is >100 times faster than exact PCoA and can be ~25% to ~33% faster than the fast PCoA in scikit-bio (Supplementary Table 10). The ordination results showed that fPCoA was nearly identical to exact PCoA for top coordinates (Supplementary Fig. 23; Procrustes M² = 0.00001, P < 0.001). Because DartUniFrac is much faster, it can be used to perform statistical resampling tests on UniFrac-based clustering (for example, unweighted pair group method with arithmetic mean (UPGMA) or neighbor joining), allowing us to quantify how robust the resulting trees are to sampling noise in community composition, a well-recognized issue in PCR-based amplicon community profiling²². For the GWMC dataset (n = 1,185), UPGMA trees obtained from jackknife resampling with DartUniFrac were highly consistent with those derived from the same jackknife procedure using exact UniFrac distances, indicating similar clustering structure (Supplementary Fig. 24). For the full EMP and AGP dataset (n = 50,085), 50 rounds of jackknife resampling with DartUniFrac (CPU) were completed in under 45 min, whereas the same procedure with exact UniFrac (Striped UniFrac) required more than 10 h (Supplementary Table 11). With multi-GPU support, the same procedure can be completed in 6.3 min using DartUniFrac, whereas it requires 2.1 h using unifrac-binaries-GPU (Supplementary Table 11). DartUniFrac also enables large-scale phylogeny-aware differential abundance analysis by providing an unbiased approximation of EMDUniFrac²³ (see detailed proof in Supplementary Methods). For example, we applied DartUniFrac as a component of all-versus-all EMDUniFrac computation to identify key bacterial species from time-series stool samples collected from a patient with inflammatory bowel disease²⁴ (Supplementary Table 12).

DartUniFrac can greatly accelerate large-scale UniFrac computation for many applications. For example, in the web-enabled microbiome meta-analysis platform Qiita¹⁸, users have deposited more than half a million microbiome samples (including metagenomic datasets) and the total number of unique taxon markers is in the tens of millions (including both amplicon sequencing and metagenomic sequencing) (accessed November 2025). It also enables routine cross-study meta-analysis at repository scale and makes sensitivity analyses (for example, jackknifing, bootstrapping, permutation-based significance testing) computationally feasible. When a fixed reference phylogeny is used—as is typical for close-reference community profiling—new samples can be appended by computing only the new rows/columns against the existing sketches, avoiding recomputation of the full DartUniFrac distance matrix. Together, these capabilities lower the computational cost and energy barrier for million-sample phylogenetic beta-diversity analyses. The algorithm can compute only unweighted and weighted UniFrac but not generalized UniFrac²⁵ or variance-adjusted UniFrac²⁶ as the latter two UniFrac algorithms cannot be formulated as weighted Jaccard similarity over tree branches. Because of this, both generalized UniFrac and variance-adjusted UniFrac remain difficult to run beyond a few thousand samples in practice despite their potential advantages in handling abundance-dependent effects, improving power for shifts in rare to moderately abundant taxa and stabilizing noisy abundance differences through variance weighting. EMDUniFrac^23,27, which was proposed recently and is essentially the numerator of weighted UniFrac (Supplementary Fig. 1b), however, cannot be efficiently approximated via other unbiased locality-sensitive hashing algorithms²⁸, whereas weighted MinHash (DartMinHash) for the original UniFrac (both unweighted and weighted) is unbiased. We found that DartUniFrac can also provide an unbiased estimation for EMDUniFrac²³ because the denominator ((sum _{j}{d}_{!j}times({a}_{!j}+{b}_{!j})=sum _{i}{{ell}}_{i}times(sum _{jin {{textrm{Desc}}}(i)}{a}_{!j}+sum _{jin {{textrm{Desc}}}(i)}{b}_{!j}))) in weighted UniFrac (Supplementary Methods) can be computed per sample using a single traversal over branches followed by constant-time pairwise addition⁷. This is much cheaper than computing the numerator exactly for each pair. Based on this unbiased EMDUniFrac estimate, large-scale differential abundance analysis that accounts for evolutionary histories is no longer computationally prohibitive. Notably, the DartUniFrac algorithm can be equally applied to absolute abundance UniFrac²⁹, where we take absolute cell counts of a given species or strain instead of relative abundance as inputs to build the sample set vectors.

Computational complexity analysis showed that DartUniFrac computation time is dominated by the number of samples but not the number of taxa, whereas exact UniFrac algorithms are dominated by O(N² × T)², which becomes impractical when the taxon count T increases with sample size—often approximately proportional to N as additional samples reveal new taxa/strains—such that the worst runtime can approach O(N³) (or even worse) in practice. In fact, none of the packages providing current exact UniFrac implementations support more than 20 million taxa and 500,000 samples (in part due to input file format limitations), which makes DartUniFrac the only option without additional software engineering effort. Crucially, this scalability is about not only a peak throughput but also about feasibility: without sketching, the memory footprint of O(N² × T) methods becomes prohibitive long before computation finishes. By decoupling computation from the raw taxon count and operating on fixed-length signatures, DartUniFrac remains tractable as reference phylogenies and feature catalogs continue to expand (for example, strain-resolved and spatially resolved metagenomics). DartUniFrac can handle real-world large-scale datasets with different sparsity; for example, per-sample amplicon denoising³⁰ and pool-sample amplicon denoising^31,32 methods can have different sparsity (~0.01% and ~5%, respectively)^4,17. The ERS algorithm is especially efficient relative to DartMinHash for datasets with >4% sparsity (Fig. 2l). With the recent advancement in spatial metagenomics technology at fine resolution²⁰, it is not uncommon to obtain even denser datasets with millions or billions of taxa, especially for homogeneous environments such as marine and freshwater environments. The balanced-parentheses data structure for representing trees is theoretically optimal³³ and can scale to trees with billions of taxa. Both DartMinHash and ERS match the current theoretical lower bound for unbiased estimation of weighted Jaccard similarity for sparse and dense datasets, achieving better accuracy and running speed than BagMinHash and Improved Consistent Weighted Sampling (ICWS)^34,35,36. Because these estimators recover similarity purely from hash collision counts, the problem of large-scale UniFrac computation in DartUniFrac essentially reduces to computing integer Hamming similarities (normalized) between fixed-length hash signatures. However, integer Hamming similarity between fixed-length integer hash signatures is purely memory-bandwidth bound, but not compute bound, indicating that hardware with better memory bandwidth is the most efficient way to speed up DartUniFrac computation. This is consistent with our experiments, where DartUniFrac on commodity GPUs (for example, NVIDIA RTX 6000 Pro or A100) was much faster than on CPUs for the integer Hamming similarity step.

The accuracy of DartUniFrac can be further improved by increasing the sketch size, at additional computational cost—for example, by performing integer Hamming similarity on longer vectors/sketches—which in turn places even greater demands on hardware memory bandwidth (both CPU and GPU) to achieve high speed. We also provide both CPU and GPU streaming modes so that a small portion of the large distance matrix (N² entries) can be computed to avoid large RAM requirements for millions of samples. In addition, the main computational bottleneck—computing the pairwise integer Hamming similarity between sketch vectors—is embarrassingly parallel and can be distributed across many computing nodes, with or without GPUs. The pairwise sketch-level computations decompose naturally into independent subproblems (blocks of the distance matrix) that can be assigned to different nodes with minimal communication overhead.

In downstream analysis, for example, PCoA, for the accuracy levels relevant to our applications, we find that power iteration-based randomized SVD achieves the desired precision while operating on a smaller working subspace than a randomized block Krylov implementation²¹. This reduces orthogonalization and memory traffic and yields a modest speedup in practice, consistent with theoretical prediction. Power iteration naturally supports a streaming/out-of-core implementation: each pass over the distance matrix requires only the current sketch, so the matrix can be processed in blocks. In contrast, randomized block Krylov methods must construct and orthogonalize a full block Krylov basis of dimension, which is difficult to maintain in a streaming setting without storing all intermediate blocks or recomputing them³⁷. For an even larger number of samples (for example, more than a million), combined with streaming mode in DartUniFrac, streaming PCoA can reduce memory requirements significantly without sacrificing accuracy.

In summary, DartUniFrac will enable the study of much more diverse environments, such as the soil at fine-grained spatial and temporal resolution and help address microbial ecological and evolutionary questions at unprecedented scales. DartUniFrac paves the way for training deep learning models for microbiomes by providing fast and accurate ground truth at the scale of millions or above.

Balanced-parentheses representation of phylogenetic trees

For a given Newick format tree, we represent the tree using a succinct, bit-level balanced-parentheses encoding instead of a pointer-based tree. During a preorder traversal, we append an open bit (1) when first visiting a node and a close bit (0) after processing its children while storing node labels (for example, taxon IDs) in a compact array in the same order. The resulting bitstring is stored as 64-bit words with a succinct rank/select index and a two-level ‘pioneer’ index for matching parentheses across word boundaries. This supports succinct-index-based tree navigation, including first-child and next-sibling traversal, using packed 64-bit words, auxiliary rank/select structures and a pioneer index for long-distance parentheses matching. In practice, these operations access a small number of contiguous machine words and auxiliary index entries, with no pointer chasing (Fig. 1a). The raw topology is stored as one bit per parenthesis, requiring 2(N + 1) bits for N tree nodes plus a virtual root, with additional overhead for the rank/select and pioneer indexes. Thus, even trees with billions of taxa remain memory feasible and efficiently traversable on modern multicore CPUs. The balanced-parentheses representations are known to be asymptotically succinct for rooted ordered trees and can approach the information-theoretic space bound while supporting efficient navigation operations³³. See Supplementary Methods for a detailed description.

DartMinHash implementation

For any given weighted set (X={(i,{x}_{i})}) with non-negative weights (here i indicates a branch, while x_i is ({{ell}}_{i}times sum _{jin {{textrm{Desc}}}(i)}{a}_{!j}) for weighted UniFrac), our bucket-style DartMinHash sketch produces k independent MinHash sketches as follows. First, we performed a parameterization step: choose k and set the expected dart budget (tapprox {k}times{textrm{ln}} {k}+2times k). Define a rank threshold θ (initially 1.0). Let (W=sum _{i}{x}_{i}) and ({r}_{max }=theta /W). Then we performed dyadic tiling of the weight–rank space (Supplementary Fig. 25). For each element (i, x_i), we iterated over dyadic scales (nu ,rho ge 0). Each pair (v, p) defines a rectangular ‘tile’ with width (Delta w={2}^{nu }/ttimes {2}^{-rho }) and height (Delta r={2}^{rho }times {2}^{-nu }). We traversed tiles in increasing w and r, stopping early when (wge {x}_{i}) or (rge {r}_{max }) (which is (frac{varphi }{parallel mathbf{x} {parallel }_{1}}); (varphi) controls the upper limit on the rank of darts to return, such that the algorithm returns (varphitimes t) darts in expectation) (Supplementary Fig. 25). Next, we started a hash-driven Poisson dart count. For every tile, we formed a 64-bit key by XORing tabular hashes of i, v, p and the tile’s local integer offsets. From this key, we draw a Poisson(1) count via a cumulative distribution function (CDF) table; this gives the number of darts landing inside the tile. For each dart assigned to a tile, we draw two independent pseudorandom values ({U}_{w},{U}_{r}sim {{textrm{uniform}}}(mathrm{0,1})) from the hash and converted them into continuous offsets (w, r) inside the tile: (w={w}_{0}+Delta wtimes {U}_{w}), (r={r}_{0}+Delta rtimes {U}_{r}). This yields a point uniformly distributed over the tile area. We accept the dart if w < x_i and r < r_max. Each accepted dart carries an ID h (another 64-bit tabular hash of the tile key plus a per-dart index) and its rank r. After this, we performed a bucketization step. We maintained k buckets or partitions, each storing the minimum rank seen thus far. A separate 64-bit tabular hasher maps h to a bucket index j = h mod k. If some buckets remain empty after one pass, we increase θ and repeat steps 2 to 5 until all k buckets are filled. For all the hashing-related steps, we used a simple tabulation hashing (32 and 64 bits), which indexes small random tables with byte/word positions and XORs the looked-up entries³⁸. It is fast and stateless and provides strong practical independence when composing multiple keys (IDs, scales and tile offsets) for MinHash-like algorithms. It has been shown that, for simple MinHash, tabulation hashing, especially mixed tabulation hashing, yields concentration essentially as good as truly random hashing for bucket/partition hashing³⁹. DartMinHash produces per-bucket minima via an exponential race; once each bucket has a winner, the collision indicators across buckets behave like independent Bernoulli trials with the mean equal to weighted Jaccard similarity. Under the hashing assumptions (as good as truly random hashing), DartMinHash has variance and concentration that match the unweighted case in practice, although mixed tabulation hashing provides even stronger guarantees³⁹. For any two resulting sketches with length k, the hash collision probability of IDs across positions is an unbiased estimator of weighted Jaccard similarity, and the variance is J(1 – J)/k. Default sketch size k is 2,048 unless specified otherwise. In practice, a sketch size around 2,000 to 3,000 is good enough for drawing biological conclusions from analyses such as PERMANOVA and PCoA based on DartUniFrac. However, users may need to increase the sketch size for analyses such as hierarchical clustering and high-recall (for example, >0.95) nearest-neighbor search, following (kge frac{1}{4{delta }^{2}}) for an expected absolute error ẟ (Supplementary Methods). See the detailed guidelines in Supplementary Table 13.

Big-O notation for DartMinHash and b-bit DartMinHash

The DartMinHash sketch step big-O notation is (Oleft(ktimes{log}(k)+right.)(parallel x{parallel }_{0})({log}left(parallel {mathbf{x}} {parallel }_{1}+right.)(frac{1}{parallel {mathbf{x}} {parallel }_{1}}left.left.right)right))¹⁵, where (parallel mathbf{x}{parallel }_{0}) is the number of nonzero elements in a set (one sample), while (parallel mathbf{x} {parallel }_{1}) is the sum of weights of all nonzero elements. This holds as long as (parallel mathbf{x} {parallel }_{1}) is not extremely small or large, which is the case for real-world microbiome datasets, DartMinHash sketching complexity is simply (O({k}times{mathrm{log}}(k)+ctimes parallel mathbf{x} {parallel }_{0})), where c is a small constant. In the case of DartUniFrac, (parallel mathbf{x} {parallel }_{0}) will be the average number of active branches for each sample, and k is the sketch size. We rewrite the sketching complexity for N samples using T_active as (parallel mathbf{x} {parallel }_{0}) and S as the sketch size: O(N × T_active + N × S × log(S)). Similar to the b-bit MinHash idea⁴⁰, the 64-bit integer hashes in sketch vectors can be reduced to b = 32 or 16 bits (default) by extracting the lower 32 or 16 bits of hashes before computing collision probability⁴⁰. The noise of the new estimator is 2^−b. For b > 14, the noise can be ignored as 2^−b is rather small.

ERS implementation

We first defined caps and red–green indexes. For a weighted set X = {(i, x_i)} (here i indicates a branch, while x_i is ({{ell}}_{i}times sum _{jin {{textrm{Desc}}}(i)}{a}_{j}), where a_j is the abundance of descendant taxa/leaves under branch i) with x_i > 0, choose tight real-valued caps m_i ≥ max x_i(s) across the dataset and build a cumulative array cum with cum[0] = 0, cum[i + 1] = cum[i] + m_i; let (M=sum _{i}{m}_{i}). This defines a line [0, M) partitioned into segments [cum[i], cum[i + 1]]. A draw (rin [0,M)) maps to component i by binary search (‘red–green’ test). Then, for efficiency, we created dense weights for membership tests. Materialize a dense vector w of length D with w_i = x_i (0 otherwise). This allows an O(1) acceptance check once the component index is found. We then generated k independent sequences. For each hash position (jin {0,ldots ,k-1}), generate a fixed-length sequence ({r}_{!j,1},ldots ,{r}_{!j,L}). Each r_j,t is produced by hashing the key (j, t) to a uniform distribution (uin [mathrm{0,1})) and setting r = μ. Next, the red–green acceptance test maps r to component i by upper bounding in cum. The draw is green (accept) if r ≤ cum[i] + x_i; otherwise, it is red (reject). Upon the first green for position j, emit the pair (ID_j, rank_j), where ID_j is a 64-bit hash derived from the accepted draw r and ({text{rank}}_{j}:=t) (the first-hit index). After the red–green test, a densification (or rotation) is performed. If a position j has no green after L attempts, mark it empty and densify draw a per-j offset (oin {1,ldots ,k-1}) and scan circularly to copy the nearest nonempty bucket’s (ID, rank). This procedure is data-independent and preserves the collision-probability estimator. We finally return k pairs (ID_j, rank_j). For any two sketches, the collision probability of IDs across positions is an unbiased estimator of weighted Jaccard similarity, and the variance is J(1 – J)/k. Tight capacities (small M) increase acceptance and allow smaller L. All randomness is via simple tabulation hashing: treat keys as words/bytes, look up in small random tables and XOR results³⁸. Separate 64-/32-bit tabulators drive (1) uniforms u for r; (2) stable IDs from r; and (3) densification offsets—yielding fast, stateless, reproducible pseudorandomness suitable for composing keys (j, t) and r. Let (p=Pr [{{textrm{green}}}]=parallel mathbf{x} {parallel }_{1}/M) under caps m_i (with (M=sum _{i}{m}_{i})). A bucket is empty only if all L trials are rejected, so ({p}_{{{textrm{empty}}}}={(1-p)}^{L}approx {e}^{-{pL}}). We choose L to make empties rare across the sketch, for example, (ktimes{p}_{{{textrm{empty}}}}le epsilon), giving (Lgtrsim mathrm{ln}(k/epsilon )/p). Thus, tighter caps (smaller M) and larger sample mass (parallel mathbf{x} {parallel }_{1}) allow smaller L; in practice, we use a conservative fixed L or adapt L using (p=parallel mathbf{x} {parallel }_{1}/M). In microbiome feature tables that are highly sparse (typically 0.1–10% nonzero), the effective ERS acceptance probability (p=parallel mathbf{x} {parallel }_{1}/M) is often on the order of 10⁻² or higher under tight dataset-wide caps; in this regime, L = 512 makes empty buckets rare for k = 2048 (because (ktimes{(1-p)}^{L}approx ktimes{e}^{-{pL}}ll 1)). Accordingly, we use L = 512 as a conservative default and increase L for datasets where the observed p is smaller. The same b-bit idea can be applied to ERS.

Sketch cost, size and set dimension

DartUniFrac approximates UniFrac distances by estimating the underlying weighted Jaccard similarity J_w between branch-length-weighted feature vectors using weighted MinHash (DartMinHash or ERS). Theoretical results for (weighted) MinHash show that the estimator J_w is unbiased and that its variance scales approximately as ({{textrm{Var}}}({hat{J}}_{{textrm{w}}})propto 1/k), where k is the sketch size and does not depend on the ambient dimensionality (for example, the number of taxa or branches in the tree). Intuitively, the sketch samples from the distribution of nonzero branch masses in each sample; adding branches with zero mass in both samples does not change this distribution and therefore does not affect the approximation error. Consequently, for a given similarity regime and typical sparsity pattern (that is, when the number of nonempty taxa per sample remains stable), keeping k fixed yields approximately constant RMSE in J_w and thus in the resulting UniFrac distances, even as the tree grows to include many more taxa and branches. This is the advantage over exact methods, where if samples contain the taxa, these taxa/branches must be built into the stripe in the Striped UniFrac algorithm and then compared in a pairwise manner, thus substantially increasing the space and computation. For DartUniFrac, due to early-stop conditions in DartMinHash, more taxa/branches will only add to the hashing cost per sample, which is a constant-time operation.

Implementation, parallelization, SIMD Hamming and GPU integer Hamming

DartUniFrac is implemented in pure Rust. All steps, including tree traversal, sample vector parsing, DartMinHash/ERS sketching and integer Hamming distance computation, are fully parallelized via the Rayon library. To further speed up integer Hamming distance computation on the CPU, we used SIMD on modern CPUs via the Rust standard library portable SIMD API. GPU implementation is based on cudarc (https://crates.io/crates/cudarc), a Rust API around the NVIDIA CUDA library, supporting multiple NVIDIA GPUs. Briefly, each sketch vector (representing each sample) is sent from CPU to GPU, and each GPU maintains a copy of all sketches. Because the sketch vector is very small (a few kilobytes), GPU memory consumption is rather small even for millions of samples (~3 GB with default b bits=16) (Fig. 2j). Each GPU is only responsible for computing a small portion of the integer Hamming similarity computation in a multi-GPU setting. The resulting integer Hamming similarity will be sent to the CPU once the computation is finished and stored in CPU RAM. In the streaming setting, a block of the distance matrix will be computed first and then sent to the CPU for writing, while the GPU starts computation for the next block of the distance matrix. Once the writing of the first block is finished, it will be released from CPU memory (RAM). In practice, writing is faster than GPU computing if the block size is not too large compared with the number of samples.

fPCoA based on randomized SVD

Instead of exact SVD in exact PCoA (the speed-limiting step)⁴¹, we use a randomized SVD based on randomized subspace iteration (simultaneous power iteration) for datasets with >5,000 samples rather than the randomized block Krylov method used in scikit-bio (v0.7.0)⁴¹. Subspace iteration typically requires orthogonalizing fewer basis vectors and can be more GEMM heavy (and thus more cache friendly for level 3 BLAS), whereas block Krylov methods offer stronger accuracy guarantees per pass³⁷. We perform the required vector orthogonalization using a LAPACK backend via the Rust lax crate. In our regime—dense N × N distance matrices with up to hundreds of thousands of samples and a modest target rank, randomized subspace iteration and randomized block Krylov have essentially the same dominant cost because both are driven primarily by a small number of passes over the distance matrix (matrix-block multiplies)³⁷. Therefore, only modest speedups are expected from switching between these two approaches for in-memory implementations. The output from fPCoA can be imported directly to scikit-bio and QIIME 2 for visualization. To compare with scikit-bio (v0.7.1) fast PCoA, the ‘fsvd’ option was used to perform randomized block Krylov SVD. k = 10 was used for both fast PCoA in scikit-bio and DartUniFrac.

Benchmarking software packages and dataset generation

We use the original Striped UniFrac (v0.9.3)⁷ and the hardware-optimized Striped UniFrac called unifrac-binaries (v1.6.0)⁹ as the truth. For GPU-based comparisons, the GPU implementation of unifrac-binaries is used. The latter represents the current fastest exact UniFrac implementation with hardware optimization. DartUniFrac v0.3.0 was used for all the benchmarks. EMP data were downloaded from the EMP project⁴. The same processing pipeline was used to obtain OTU feature tables and phylogenetic trees. Briefly, sequence data were error filtered and trimmed to the length of the shortest sequencing run (90 bp) using Deblur (v1.1.1)³⁰. Deblur observation tables were filtered to keep only tag sequences with at least 25 total reads across all samples. For comparison to existing OTU tables, traditional closed-reference OTU picking was done against the 16S rRNA databases Greengenes 13.8 (ref. ⁴²) and SILVA 123 (ref. ⁴³) using SortMeRNA (v6.0.2)⁴⁴ and subsampled open-reference OTU picking was done against Greengenes 13.8. The AGP data were also obtained from the original paper⁵, and a similar processing approach was used to obtain the OTU feature table and the phylogenetic tree. Then, OTU feature tables from EMP and AGP were merged (the same closed-reference 16S rRNA gene database and phylogeny were used) to build a combined EMP and AGP dataset. The GMTOL dataset, which is a subset of EMP with samples that were taken from animal gut, was generated by mapping the metadata from EMP. We then also retrieved all amplicon studies (278,499 samples, ~1 million taxa) from Qiita¹⁸ via Rediom (v0.3.9)⁴⁵ and applied a similar sequence processing pipeline as above to obtain OTU feature tables and phylogeny but against the most recent Greengenes 2 database⁴⁶, a comprehensive 16S rRNA-based taxonomy–phylogeny framework anchored in a whole-genome phylogeny. The three-country THDMI (The Healthy Diet and Microbiome Initiative) dataset was downloaded from the original paper⁴⁷. Briefly, quality-controlled reads from shotgun metagenomic sequencing were mapped to Web of Life version 2 (ref. ⁴⁸) using Bowtie 2 (v2.5.4)⁴⁹ with the SHOGUN (v1.0.0)⁵⁰ parameter set in paired-end mode with a feature table estimated from Woltka (v0.1.7)⁶. Then, the abundance of each genome was obtained by counting the number of reads mapped to the genome with a minimum coverage breadth filter (0.1%)⁵¹. The phylogenetic tree was provided by the abovementioned genome database. For simulated datasets, for example, 100,000, 300,000, 500,000 and 1,000,000 samples, all taxa from EMP and AGP (875,222 taxa) and Greengenes 2 (23,450,268 taxa) were used as input features, and then taxon abundances for each sample were generated from a normal distribution with adjustable sparsity (for example, 0.05–5%). The code repository for simulating feature tables can be found in the Code availability section. All statistical analyses, including PERMANOVA, Procrustes, Mantel and hierarchical clustering topology analyses, were performed using scikit-bio (v0.7.0)⁴¹ and SciPy (v1.16.1)⁵².

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.