New k-mer estimator cracks substitution rates in repetitive DNA

Repetitive DNA has long been a blind spot for one of genomics’ favorite shortcuts. K-mer-based analysis, which compares genomes by counting short strings of length k rather than aligning them base by base, underpins genome assembly, error correction, read mapping, variant calling, genotyping, database search, metagenomic comparison, and alignment-free sequence comparison. But when a genome is riddled with repeats, the counts that these methods rely on become distorted, and the math behind them quietly breaks down.

The new method, described in Algorithms for Molecular Biology by Haonan Wu, Antonio Blanca, and Paul Medvedev, targets a deceptively simple question: given two related sequences, how many substitutions separate them? The Mash tool, a landmark in sketch-based genomics, can answer this accurately for sequences with few repeats by estimating the mutation rate from the Jaccard similarity between their k-mer sketches. That speed allowed the Mash authors to build a phylogeny of 17 primate species in a tiny fraction of the time an alignment-based method would need. Yet the derivation behind Mash’s estimator assumes the genomes contain no repeats, and on highly repetitive regions such as the centromere of a human chromosome it becomes inaccurate. Follow-up work by Blanca and colleagues in 2021 modeled how mutations affect k-mer sets, but only under strong assumptions that the sequence is non-repetitive and that mutations never create spurious k-mer matches. The theoretical foundations for extending such an estimator to repeat-rich sequences have, until now, been lacking.

The core of the problem lies in what happens to a k-mer when mutations strike. Under the substitution process the researchers model, each nucleotide in a source string s stays unchanged with probability 1 minus r and flips to one of the three other nucleotides with probability r divided by 3, where r is the substitution rate. For a k-mer that appears only once in s, the chance it survives intact in the mutated string t is simply (1-r) raised to the power k. The repeat-oblivious estimator, the only previously known approach and the kind used by Mash, builds on this: it assumes s has no repeats and that a mutated k-mer never becomes another k-mer already present in s, then solves the equation that the expected number of shared k-mers equals L times (1-r) to the k, where L is the number of k-mer positions. But for a k-mer occurring multiple times, survival is not that simple. The correct probability that at least one copy persists is 1 minus q raised to the occurrence count, where q is the probability that a single k-mer occurrence contains at least one substitution. Because the occurrence count times (1-q) always exceeds 1 minus q to that count on the interval from 0 to 1, the repeat-oblivious estimator systematically overestimates the rate, and the gap widens as sequences grow more repetitive, as rates fall, and as k shrinks.

Wu, Blanca, and Medvedev’s solution occupies a deliberate middle ground. Rather than fully modeling every nuance of repeat structure, which would be difficult to derive, analyze, and compute, their estimator captures its most essential feature: the abundance histogram of the source string, which records how many distinct k-mers occur exactly once, twice, and so on. Given the observed number of shared k-mers between s and t and this histogram, the estimator is defined as the unique solution to a polynomial equation, found numerically with Newton’s method. The team proved a theorem bounding the estimator’s bias in terms of properties of s alone, so a user can decide in advance whether to trust a particular estimate. A companion theorem supplies a dynamic-programming heuristic that computes, from only the k-mer count L, the k-mer size k, and the rate r, the probability that the k-mer intersection comes back empty, flagging parameter settings where any intersection-based method is doomed to fail.

The empirical tests were deliberately brutal. The researchers extracted four sequences from the human T2T-CHM13v2.0 reference genome to span the repetitiveness spectrum. D-easy, an ordinary substring of chromosome 6, has fewer than 1 percent of its k-mers appearing more than once. D-med is RBMY1A1, a Y chromosome gene built from ALU elements, SINEs, LINEs, and simple repeats, with roughly 3 percent non-singleton k-mers. D-hard is a tandem repeat containing 4.2 similar copies of a 545-nucleotide unit, with more than 40 percent of k-mers non-singleton. D-hardest is a 100,000-nucleotide slice of the active alpha satellite higher-order repeat region of the chromosome 21 centromere, where over 70 percent of k-mers are non-singletons. Alpha satellite DNA, arranged as 171-base-pair monomers stacked into higher-order repeats, is notoriously difficult to assemble or map.

The results were stark. On the two hardest datasets, the new estimator stayed within a few percent of the true substitution rate across rates from 0.1 to 24 percent, while the repeat-oblivious estimator’s reliable range shrank dramatically, covering only about 10 to 24 percent on D-hardest. At a true rate of 1.1 percent on that dataset, the repeat-oblivious estimator averaged 10.4 percent, nearly a tenfold overestimate, while the new method returned 1.1 percent. The failure mode is algebraic: the old estimator approximates a repeated k-mer’s survival probability as the occurrence count times (1-q), which inflates the estimate. On the easier datasets the old method performed nearly as well as the new one, except at very low rates, where it posted a 230 percent relative error at a rate of 0.1 percent on D-med. Above roughly 24 percent on the hard datasets, the observed intersection was frequently empty, and both estimators defaulted to reporting 100 percent, an instability the authors attribute to a fundamental constraint shared by all intersection-based estimators.

The team also mapped the interplay between k and r across a wide grid, measuring average relative absolute error over 100 replicates per cell. Smaller k increases the number of repeats and makes estimation harder, while larger k or r raises q, the per-occurrence mutation probability, pushing toward empty intersections. Since users typically choose k but only know a rough range for r, and since substitution rates above 25 percent are unlikely for biologically functional sequences, the heatmap serves as a practical guide. Across all four datasets, choosing k between 10 and 20 worked well. Simulations also confirmed the theoretical bias bounds: the empirical mean of the estimator usually fell within the predicted bounds, the upper bound proved nearly tight because overlapping k-mer occurrences are rare in these data, and where the lower bound misbehaved on D-hardest at low rates, the culprit was a loose variance bound rather than the estimator itself.

Sketching, the trick that makes k-mer methods scale to massive datasets, slots neatly into the new framework. Using a FracMinHash sketch, which keeps only k-mers whose hash falls below a threshold theta, the researchers simply divided the sketched intersection size by theta before feeding it into the same polynomial equation. On D-hardest, sketching introduced no systematic bias, only increased variance that grows as theta shrinks, letting users trade sketch size for precision. In a case study with real evolutionary stakes, the team turned to the X chromosome centromeres of five great apes, chimpanzee, bonobo, gorilla, Bornean orangutan, and Sumatran orangutan, extracted from recent telomere-to-telomere assemblies, with lengths ranging from 1.7 megabases in chimpanzee to 8.2 megabases in Sumatran orangutan. Using a FracMinHash sketch with theta of 0.01 and k of 31, their method sampled orders of magnitude fewer k-mers than Mash with its million-k-mer sketch. Because centromeres evolve through mechanisms more complex than simple substitution, the computed distances were not expected to track real evolutionary distance linearly, but they proved accurate enough for tree building: both estimators recovered the correct species tree with Neighbor Joining, yet only the new estimator produced the correct topology with UPGMA, a method more sensitive to systematic distance biases, while Mash’s UPGMA tree misplaced Gorilla.

The method does have limits. It requires knowing the abundance histogram of the source string, though tools like GenomeScope can estimate it from sequencing data, or a related genome’s histogram can serve as a stand-in. Adapting the estimator to unassembled sequencing reads, which its design makes amenable since it ignores k-mer multiplicities in the intersection, remains future work, as does tightening the bias bounds and deriving confidence intervals, a harder problem requiring deeper understanding of the estimator’s distribution. Runtime is not a concern: the estimator completes in under a second, since Newton’s method converges in on the order of log log iterations for this polynomial. The software is available on GitHub, and the authors have released all simulated datasets and scripts to reproduce every figure.

Subject of Research: A repeat-aware k-mer estimator for substitution rates between repetitive genomic sequences

Article Title: A k-mer-based estimator of the substitution rate between repetitive sequences

Article References: Wu, H., Blanca, A., & Medvedev, P. (2026). A k-mer-based estimator of the substitution rate between repetitive sequences. Algorithms for Molecular Biology, 21(1), Article 14. https://doi.org/10.1186/s13015-026-00304-1

Image Credits: AI Generated

DOI: 10.1186/s13015-026-00304-1

Keywords: k-mer, substitution rate, repetitive sequences, centromere, alpha satellite, Mash, sketching, FracMinHash, phylogenetics, bias bounds, genome assembly, T2T-CHM13v2.0