Exact sketch read mapping algorithm: recall 0.88 vs 0.76 for minimap2 on simulated chromosome Y reads

Researchers have formulated the read mapping problem at the level of sequence sketches and described an exact dynamic programming algorithm that finds all hits above a given similarity threshold, according to a paper in Algorithms for Molecular Biology. On simulated long reads from the T2T assembly of human chromosome Y, at a precision equivalent to minimap2, the recall of the algorithm is 0.88, compared with 0.76 for minimap2 biomedcentral.comAlgorithms for Molecular Biologybiomedcentral.comOpen source ↗.

Traditionally, a similar sequence meant one with a high alignment score, and read mappers were treated as heuristic solutions to that problem. For sketch-based mappers there was no problem formulation describing what an exact algorithm should solve, and no sketch-based method could find all mapping positions for a read above a score threshold.

Problem formulation and algorithm

The authors define the task on sketches: given a read sketch and a reference sketch, report all candidate intervals of the reference whose score exceeds a threshold. A candidate must begin and end with k-mers that appear in the read sketch, and it must not be contained in a longer interval with a strictly greater score. The authors also introduce a parameterized linear score function and prove that it separates candidate mappings at least as well as the weighted Jaccard index.

The algorithm runs in O(|t| + |p| + ℓ²) time and O(ℓ log ℓ) space. Here |t| is the number of k-mers in the reference sketch, |p| is the number of k-mers in the read sketch, and ℓ is the number of times k-mers from the read sketch occur in the reference sketch biomedcentral.comAlgorithms for Molecular Biologybiomedcentral.comOpen source ↗. Memory is kept low by computing the dynamic programming matrix column by column and holding only two columns at a time.

Test on human chromosome Y

Chromosome Y contains ampliconic regions, where it is desirable to find all good mapping positions rather than a single best one. The reads were simulated to mimic Pacific Biosciences HiFi data: lengths followed a gamma distribution with a mean of 9,000 base pairs, and the overall mutation rate was 0.2 percent across substitutions, insertions and deletions. After reads from repetitive centromeric and telomeric arrays that aligned to more than twenty locations were removed, 32,295 reads remained.

The algorithm was compared with minimap2, Winnowmap2 and an exact alignment-based tool built on edlib. It used the same minimizer parameters as minimap2: k-mer length 15 and window size 10, with minimizers occurring more than 100 times in the reference sketch masked. Ground truth was established with BLAST alignments. At precision equivalent to minimap2, the recall was 0.88 against 0.76.