Nature Biotechnology
(2026) Cite this article
UniFrac measures phylogeny-aware differences between microbiome samples but scales poorly with modern dataset sizes. We introduce an algorithm, DartUniFrac, and a near-optimal implementation with graphics processing unit acceleration that is up to three orders of magnitude faster than UniFrac and scales to millions of samples (pairwise) and billions of taxa. DartUniFrac connects UniFrac with weighted Jaccard similarity and exploits sketching algorithms for fast computation.
UniFrac1,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)4 and American Gut Project (AGP)5, 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 R2 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 datasets2. 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)6. 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)7. Model-based estimates of microbial diversity using the 16S rRNA gene suggest that Earth may harbor more than 1012 microbial species8 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 years7,9,10, for example, faster tree traversal strategies10, better parallel efficiency (Striped UniFrac)7 and hardware acceleration (for example, single instruction, multiple data (SIMD) and graphics processing unit (GPU))9. 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 algorithm1 (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 xj(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)2. 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 scale11,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, respectively15,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 (Jw) 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 practice15, to the best of our knowledge. Overall, the big-O notation for DartUniFrac is O(N × Tactive + N × S × log(S) + N2 × S), where N is the number of samples, Tactive 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(N2 ×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(N2× 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).
a, Schematic overview of the DartUniFrac algorithm and implementation. The phylogenetic tree is represented by optimal balanced parentheses for efficient retrieval of branch lengths relevant for each sample. Then, the input set vectors (each dimension represents a branch of the tree, skipped if not relevant for that sample; there can be millions of branches or more for big trees) are collected and passed to sketching algorithms, DartMinHash (default) or ERS. The output from sketching is a much smaller vector (1,000–2,000, called sketch vectors) storing 64-bit integer hashes (truncated to lower b bits afterward; Supplementary Methods). The number of equal hashes for all slots in the sketching vector out of the total slots equals the weighted Jaccard similarity. CPU multithreading and SIMD or multi-GPU processing can be used to speed up the pairwise computation of equal hashes. DM, distance matrix; Int, integer. b, Scalability of the DartUniFrac-CPU implementation with respect to the number of CPU threads. c, Speedup of DartUniFrac-CPU relative to Striped UniFrac and unifrac-binaries-CPU (hardware-optimized Striped UniFrac, here SIMD for unifrac-binaries-CPU) for weighted UniFrac with respect to the number of samples. d, Speedup of DartUniFrac-GPU relative to unifrac-binaries-GPU for 100,000, 300,000 and 500,000 samples. This benchmark was performed on an Intel Xeon 6972 P 96-thread CPU with NVIDIA A100 (HBM2e) GPUs. e, GPU memory requirements for DartUniFrac-GPU and unifrac-binaries-GPU. f, Procrustes analysis comparing DartUniFrac and Striped UniFrac for the GWMC dataset (unweighted). Procrustes M2 is shown in the plot. P = 0.00038. See Supplementary Fig. 3 for the weighted analysis. g, Mantel correlation coefficients for distance matrices estimated by DartUniFrac and Striped UniFrac (truth) for weighted UniFrac across 5 datasets (all P values < 0.001; Supplementary Table 2). DartMH, DartMinHash. h, RMSE of the DartUniFrac DartMinHash estimator versus the sketch size for the GWMC dataset. Each data point is the mean of n = 10 independent computational replicate runs, each initialized with a different random seed. Error bars represent the s.d. across replicates. i, Mantel correlation coefficients, using Pearson and Spearman correlations, of the DartMinHash estimator versus the sketch size for the GWMC dataset. Each data point is the mean of n = 10 independent computational replicate runs, each initialized with a different random seed. Error bars represent the s.d. across replicates.
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 UniFrac7 and a hardware-optimized version, called unifrac-binaries9; 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) (M2 = 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 M2, PERMANOVA effect size (pseudo-F) and R2 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 datasets17 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)18 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)19, a common container for microbiome feature tables, which is capped at 220 nonzero values in the feature table; with feature tables that extend the limit for nonzero values to 264, 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).
a, Branch set construction time with respect to the number of taxa and sparsity while fixing the number of samples. b, DartMinHash sketching time with respect to the number of samples. c, Running time comparisons between DartUniFrac-GPU in-memory mode and streaming mode. d, Peak RAM comparisons between DartUniFrac-GPU in-memory mode and streaming mode. e–h, PCoA visualization for the EMP dataset for DartUniFrac, unweighted (e) and weighted (f), respectively, and Striped UniFrac, unweighted (g) and weighted (h) (exact UniFrac). The PCoA explanation rate for each principal coordinate is included next to the axis. Data points were grouped by country. A detailed legend for e–h can be found in Supplementary Fig. 11. i, GPU running time for DartUniFrac versus CPU running time for DartUniFrac with respect to the number of samples and the number of GPUs. j, DartUniFrac-GPU memory requirement with respect to the number of samples. VRAM, video RAM. k, Estimation error of our DartMinHash implementation for simulated UniFrac distance between (0, 1), including RMSE, absolute error and theoretical error. A sketch size of 2,048 was used for this analysis. Each data point for RMSE and absolute error is the mean of n = 10 independent computational replicate runs, each initialized with a different random seed. Error bars represent the s.d. across replicates. l, Running time versus sparsity of input weighted sets for DartMinHash and ERS. Sparsity was measured as the fraction of nonzero values considering all the set/branch elements.
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 expectations15. 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 SVD21 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 M2 = 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 profiling22. 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 EMDUniFrac23 (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 disease24 (Supplementary Table 12).
DartUniFrac can greatly accelerate large-scale UniFrac computation for many applications. For example, in the web-enabled microbiome meta-analysis platform Qiita18, 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 UniFrac25 or variance-adjusted UniFrac26 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. EMDUniFrac23,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 algorithms28, 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 EMDUniFrac23 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 addition7. 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 UniFrac29, 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(N2 × T)2, 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(N3) (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(N2 × 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 denoising30 and pool-sample amplicon denoising31,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 resolution20, 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 optimal33 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 (N2 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 implementation21. 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 them37. 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.
