Tomtom-lite: Accelerating Tomtom enables large-scale and real-time motif similarity scoring

demonstrating that our re-implementation can achieve over a thousand-fold speedup18 compared with the original tool on reasonable tasks (see Supplementary Note 2).19 Availability and Implementation: Our implementation of Tomtom is freely available as20 a Python package at https://github.com/jmschrei/memesuite-lite, which can be down-21 loaded via pip install memelite.22 1 Introduction23 Calculating the similarity between a pair of sequences has been a central operation in genomics24 since its inception [1]. Initially, these algorithms assumed that both sequences would be25 discrete entities [2] (e.g., composed of nucleotides or amino acids), and were optimized given26 this constraint [3]. Yet, as the field expanded to the discovery of transcription factor (TF)27 binding sites and their regulatory effects, the motifs that were being discovered were not28 themselves discrete. Rather, because TFs do not need to bind to perfect matches, these29 motifs were encoded as position-weight matrices (PWMs) that were frequency-based, where30 each column encoded probabilities for each nucleotide appearing at each position [4]. New31 tools were necessary for learning PWMs from data [5, 6, 7], scanning these PWMs against32 discrete sequences [8, 9], and against each other. Tomtom emerged almost two decades ago33 as a popular algorithm for calculating similarities between a set of query PWMs and a set34 of target PWMs and subsequently converting these scores into p-values that account for the35 length and information content of the motif [10, 11].36 1 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint The core conceptual innovation of Tomtom is the faithful calculation of null distributions.37 Briefly, the challenge is that one cannot assume the query and target PWMs are aligned or38 even the same length. Consequently, most PWM similarity calculations proceed by considering39 all possible ungapped alignments (and potentially reverse complements) and returning the40 maximum similarity score (Fig 1A). Because arriving at this score requires many operations,41 the null distributions for the statistical test must be calculated in a way that represents this42 entire series of operations. In their original work, Gupta et al. [10] propose Tomtom as an43 efficient algorithm for doing so, and a modification of Tomtom by Tanaka et al. [11] adjusts44 the similarity scores based on the number of unaligned positions. See Supplementary Note 145 for a complete description of the Tomtom algorithm, including implementation details.46 Being able to quickly calculate similarities between PWMs has recently taken on newfound47 importance as machine learning (ML) models have been integrated into genomics analyses.48 Two of the most prevalent forms of ML models have been supervised models that directly49 predict genomic activity [12, 13], and language models [14, 15] that capture the distribution50 of the genome [16, 17, 18, 19]. Although these models differ in significant ways, both have51 feature attribution methods that reveal which nucleotides drive model predictions, e.g.in silico52 saturation mutagenesis [20]. Spans of high-attribution nucleotides, called “seqlets”, can then53 be identified [21]. The computational challenge then becomes annotating these seqlets with54 the TF (or TF family) whose binding they most resemble. This can be done either by mapping55 the discrete sequence within seqlets to a motif database, or by mapping a PWM containing56 these attribution values (sometimes called a contribution-weight matrix or CWM [22]) to a57 database of attribution-based motifs. Tomtom is a natural solution to both problems.58 Unfortunately, three practical challenges have limited Tomtom’s adoption. First, Tomtom59 does not scale well to large target databases despite being implemented in C. This is due, in60 large part, to the unnessesary re-calculation of null distributions for each query-target pair.61 Second, Tomtom is not natively multi-threaded and so requires additional command-line tools62 and expertise to utilize modern multicore computing systems. The third is that usage from63 an interactive environment or Jupyter notebook [23] requires writing the queries and targets64 out to disk, running the command-line tool, and then reading the results back from disk into65 memory. These steps can be burdensome and are error prone when trying to quickly do66 interactive analyses.67 Our implementation of Tomtom (referred to hereafter as “tomtom-lite”) overcomes these68 challenges with a numba-accelerated Python function. Upon their first use, numba functions69 are compiled down to machine code that are similarly fast to their C counterparts [24]. These70 functions can use multi-threaded parallelism without significant code changes. Additionally,71 tomtom-lite caches null distributions across query-target pairs, organizes operations to avoid72 cache misses, uses an approximate median in a computational bottleneck, and allows for hash-73 ing of columns in the target database to reduce redundancy in the calculations. Together,74 these improvements make tomtom-lite up to three orders of magnitude faster than Tomtom in75 our evaluations.76 As an initial check of our implementation, we compared the p-values produced by Tomtom77 and tomtom-lite when comparing the entire JASPAR2024 motif database [25] against itself.78 We observed a high correlation between their p-values (r=0.99981, Fig 1B), with tomtom-lite79 being ∼417x faster than Tomtom (∼2420s for Tomtom, ∼5.8s for tomtom-lite). Importantly,80 the p-values have two regimes: r = 0.99981 when p > 1e − 5 containing 99.5% of comparisons,81 and r = 0.97 when p ≤ 1e − 5, where the influence of the approximations is greater. These82

suggest that the approximations have little effect when the goal is more discriminative83 in nature, e.g., identifying which targets might be good matches, but might influence results84 more when one cares about precisely ranking strong hits against each other.85 We then evaluated how the two implementations scale with respect to various data prop-86 2 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint Calculate Similarity Between o and o Alignment Similarity Scores Score PDF TOMTOM Background Distributions Query T arget (Aligned) T arget (Unaligned) p-value 0 20 40 60 80 Genomic Position 0.04 0.02 0.00 0.02 0.04 Attributions HAP3 (0.266) NFYB (0.172) ZNF140 (-0.0252) 10 13 10 10 10 7 10 4 10 1 TOMTOM p-value 10 14 1 0 11 10 8 10 5 10 2 tomtom-lite p-value x = y 10 0 10 1 10 2 10 3 10 4 # Queries 10 1 10 0 10 1 10 2 Time (s) tomtom-lite (1 thread) tomtom-lite (2 threads) tomtom-lite (4 threads) tomtom-lite (8 threads) tomtom-lite (16 threads) tomtom-lite (32 threads) TOMTOM (1 thread) 10 0 10 1 10 2 # Queries 10 1 10 2 Time (s) 10 0 10 1 10 2 10 3 10 4 # Targets 10 1 10 0 10 1 10 2 Time (s) Seqlet Count 0 5000 10000 15000 20000 25000 Arnt PAX6 RORA RXRA::VDR REL RELA NR1H2::RXRA Znf423 NFIC::TLX1 ZNF354C Pou5f1::Sox2 EWSR1-FLI1 Arid3a INSM1 RARA::RXRA JASPAR2024 Top Motif Hit Count G. A. C. F . D. E. B. Figure 1: A schematic and results of our Tomtom implementation. (A) A schematic demonstrating that for each possible alignment (each row of purple/pink boxes) a score is calculated as the sum of all aligned columns between the target and the query (in pink), and the maximum score is kept and converted to a p-value. (B) P-values from Tomtom and tomtom-lite when comparing JASPAR2024 against itself. Only 1% of the >5M points are displayed. (C) Timings for Tomtom (in gray) and tomtom-lite with an increasing number of queries. (D) Timings when using a database of 1M random PWMs. (E) Timings when considering a target database of increasing size. (F) An example attribution track with seqlets called and annotated using Tomtom and JASPAR. (G) Count of the number of seqlets mapped to each motif in JASPAR according to tomtom-lite. 3 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint erties. See Supplementary Note 2 for details on our exact evaluations. Both implementations87 scaled linearly with the number of queries but, due to differing scaling constants, tomtom-lite88 was 43x faster with 1 query and 265x faster with 128 queries (Fig 1C). Running tomtom-lite89 with 16,384 queries and 32 threads was faster than running Tomtom with 16 queries. We90 noticed that additional threads diminished in improvement past 8 and hypothesized that this91 was because the communication cost of managing that number of threads became higher than92 each task itself. Accordingly, we considered a more compute-intensive task where the target93 database was 1 million randomly generated PWMs that were longer, on average, than those in94 JASPAR. In this evaluation, we observed better scaling up to 32 threads (Figure 1D), suggest-95 ing that more threads will not always help when the target database is too small. Next, we96 found that both implementations scaled linearly with the number of targets but tomtom-lite97 scales significantly better because it caches null distributions. With only 1 target, tomtom-lite98 was only 2.1x faster with 1 thread and 2.89x faster with 8 threads, but with 16,384 targets99 tomtom-lite is 528x faster with 1 thread and 1,015x faster with 8 threads (Fig 1E).100 We demonstrate these speed improvements in two practical settings: real-time interactive101 analysis, and large-scale motif similarity scoring. The first setting is invaluable for the explo-102 ration of data and model predictions, the generation of hypotheses, and the prototyping of103 new methods. The second setting is more representative of the computational burdens faced104 by analysis tools aiming to scale genome-wide or further. In both settings, we consider down-105 stream applications of a ML model, ChromBPNet [12], that makes predictions for chromatin106 accessibility (specifically, ATAC-seq) in K562 from nucleotide sequence.107 In our interactive example, we consider applying ChromBPNet to a specific locus of interest108 to identify motifs driving accessibility. Here, we follow the standard protocol for calculating109 DeepLIFT/SHAP attributions [26, 27] (see Supplementary Note 2 for details), identify seqlets,110 and then use tomtom-lite to annotate these seqlets using JASPAR2024. Our procedure identi-111 fied 4 seqlets (Fig 1F, only 3 seqlets visualized), which tomtom-lite annotated in ∼0.1 seconds112 and Tomtom annotated in ∼2.5 seconds.113 In our larger scale example, we consider the goal of counting the number of seqlets genome-114 wide that map to each motif in JASPAR2024. Specifically, we consider ChromBPNet attribu-115 tions across all 203,804 ATAC-seq peaks in K562, which result in 634,776 seqlets. Annotating116 these seqlets using the JASPAR2024 database, which has 2,346 entries, takes only ∼1,260117 seconds ( ∼20 minutes) using 8 threads. This speed makes it possible for one to perform118 genome-wide seqlet annotation on a commercially available laptop. When we consider the top119 motifs by count, we find that Arnt, PAX6, and RORA are in the top 15 motifs (Fig 1G).120 Although our evaluations here focus on traditional PWMs where each column contains a121 value for each nucleotide or amino acid, Tomtom and tomtom-lite do not require this. Rather,122 they only requires that each column represent a position and that the features are semantically123 consistent across sequences. This means that one could natively apply Tomtom to the setting124 where PWMs contain embeddings of each nucleotide, such as those derived from a language125 model. Given the popularity of language models across domains and their demonstrated126 usefulness when applied to the genome, we anticipate that motif representations derived in127 this manner may augment frequency- or attribution-based PWMs.128 Our results demonstrate that tomtom-lite is significantly faster than Tomtom and that129 these speed improvements have important practical benefits. These speed improvements are130 most beneficial at both ends of the compute spectrum: when the number of queries is small,131 potentially because one is doing interactive analyses, and when the number of queries and132 targets is very large. Given these benefits, we anticipate that tomtom-lite will be invaluable133 for tools involving genomic annotation and discovery.134 4 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint

We thank Austin Wang, Justin Sanders, and Nikolaus Mandlburger135 for their feedback on the manuscript, and William Noble and Timothy Bailey for feedback and136 discussions on the original Tomtom implementation. Research at the Institute of Molecular137 Pathology (IMP) is supported by Boehringer Ingelheim GmbH and the Austrian Research138 Promotion Agency (FFG, FO999902549). For the purpose of Open Access, the authors have139 applied a CC BY public copyright license to any Author Accepted Manuscript (AAM) version140 arising from this submission. The computational results presented were obtained using the141 CLIP cluster (https://clip.science).142 References143 [1] T F Smith and M S Waterman. Identification of common molecular subsequences. J.144 Mol. Biol., 147(1):195–197, March 1981.145 [2] Waqar Haque, Alex Aravind, and Bharath Reddy. Pairwise sequence alignment algo-146 rithms: a survey. In Proceedings of the 2009 conference on Information Science, Technol-147 ogy and Applications, New York, NY, USA, March 2009. ACM.148 [3] F Chiaromonte, V B Yap, and W Miller. Scoring pairwise genomic sequence alignments.149 Pac. Symp. Biocomput., pages 115–126, 2002.150 [4] Arttu Jolma, Jian Yan, Thomas Whitington, Jarkko Toivonen, Kazuhiro R Nitta, Pasi151 Rastas, Ekaterina Morgunova, Martin Enge, Mikko Taipale, Gonghong Wei, Kimmo Palin,152 Juan M Vaquerizas, Renaud Vincentelli, Nicholas M Luscombe, Timothy R Hughes,153 Patrick Lemaire, Esko Ukkonen, Teemu Kivioja, and Jussi Taipale. DNA-binding speci-154 ficities of human transcription factors. Cell, 152(1-2):327–339, January 2013.155 [5] T L Bailey and C Elkan. Fitting a mixture model by expectation maximization to discover156 motifs in biopolymers. Proc. Int. Conf. Intell. Syst. Mol. Biol. , 2:28–36, 1994.157 [6] Xiaotu Ma, Ashwinikumar Kulkarni, Zhihua Zhang, Zhenyu Xuan, Robert Serfling, and158 Michael Q Zhang. A highly efficient and effective motif discovery method for ChIP-159 seq/ChIP-chip data using positional information. Nucleic Acids Res., 40(7):e50, April160 2012.161 [7] Timothy L Bailey. STREME: accurate and versatile sequence motif discovery. Bioinfor-162 matics, 37(18):2834–2840, September 2021.163 [8] Janne Korhonen, Petri Martinm¨ aki, Cinzia Pizzi, Pasi Rastas, and Esko Ukkonen.164 MOODS: fast search for position weight matrix matches in DNA sequences. Bioinformat-165 ics, 25(23):3181–3182, December 2009.166 [9] Charles E Grant, Timothy L Bailey, and William Stafford Noble. FIMO: scanning for167 occurrences of a given motif. Bioinformatics, 27(7):1017–1018, April 2011.168 [10] Shobhit Gupta, John A Stamatoyannopoulos, Timothy L Bailey, and William Stafford169 Noble. Quantifying similarity between motifs. Genome Biol. , 8(2):R24, February 2007.170 [11] Emi Tanaka, Timothy Bailey, Charles E Grant, William Stafford Noble, and Uri Keich.171 Improved similarity scores for comparing motifs. Bioinformatics, 27(12):1603–1609, June172 2011.173 [12] Anusri Pampari, Anna Shcherbina, Evgeny Z Kvon, Michael Kosicki, Surag Nair,174 Soumya Kundu, Arwa S Kathiria, Viviana I Risca, Kristiina Kuningas, Kaur Alasoo,175 William James Greenleaf, Len A Pennacchio, and Anshul Kundaje. ChromBPNet: bias176 factorized, base-resolution deep learning models of chromatin accessibility reveal cis-177 regulatory sequence syntax, transcription factor footprints and regulatory variants.bioRx-178 ivorg, page 2024.12.25.630221, January 2025.179 5 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint [13] ˇZiga Avsec, Vikram Agarwal, Daniel Visentin, Joseph R Ledsam, Agnieszka Grabska-180 Barwinska, Kyle R Taylor, Yannis Assael, John Jumper, Pushmeet Kohli, and David R181 Kelley. Effective gene expression prediction from sequence by integrating long-range in-182 teractions. Nat. Methods, 18(10):1196–1203, October 2021.183 [14] Eric Nguyen, Michael Poli, Marjan Faizi, Armin Thomas, Callum Birch-Sykes, Michael184 Wornow, Aman Patel, Clayton Rabideau, Stefano Massaroli, Yoshua Bengio, Stefano185 Ermon, Stephen A Baccus, and Chris R´ e. HyenaDNA: Long-range genomic sequence186 modeling at single nucleotide resolution. arXiv [cs.LG] , June 2023.187 [15] Hugo Dalla-Torre, Liam Gonzalez, Javier Mendoza-Revilla, Nicolas Lopez Car-188 ranza, Adam Henryk Grzywaczewski, Francesco Oteri, Christian Dallago, Evan Trop,189 Bernardo P de Almeida, Hassan Sirelkhatim, Guillaume Richard, Marcin Skwark, Karim190 Beguir, Marie Lopez, and Thomas Pierrot. Nucleotide transformer: building and evaluat-191 ing robust foundation models for human genomics. Nat. Methods, pages 1–11, November192 2024.193 [16] Maxwell W Libbrecht and William Stafford Noble. Machine learning applications in194 genetics and genomics. Nat. Rev. Genet., 16(6):321–332, June 2015.195 [17] Jacob Schreiber and Ritambhara Singh. Machine learning for profile prediction in ge-196 nomics. Curr. Opin. Chem. Biol. , 65:35–41, December 2021.197 [18] Xuehai Hu, Alisdair R Fernie, and Jianbing Yan. Deep learning in regulatory genomics:198 from identification to design. Current Opinion in Biotechnology, 79:102887, February199 2023.200 [19] G¨ okcen Eraslan,ˇZiga Avsec, Julien Gagneur, and Fabian J Theis. Deep learning: new201 computational modelling techniques for genomics. Nat. Rev. Genet., 20(7):389–403, July202 2019.203 [20] Gherman Novakovsky, Nick Dexter, Maxwell W Libbrecht, Wyeth W Wasserman, and204 Sara Mostafavi. Obtaining genetics insights from deep learning via explainable artificial205 intelligence. Nat. Rev. Genet., 24(2):125–137, February 2023.206 [21] Avanti Shrikumar, Katherine Tian, ˇZiga Avsec, Anna Shcherbina, Abhimanyu Banerjee,207 Mahfuza Sharmin, Surag Nair, and Anshul Kundaje. Technical note on transcription fac-208 tor motif discovery from importance scores (TF-MoDISco) version 0.5.6.5. arXiv [cs.LG],209 October 2018.210 [22] ˇZiga Avsec, Melanie Weilert, Avanti Shrikumar, Sabrina Krueger, Amr Alexandari, Khy-211 ati Dalal, Robin Fropf, Charles McAnany, Julien Gagneur, Anshul Kundaje, and Julia212 Zeitlinger. Base-resolution models of transcription-factor binding reveal soft motif syntax.213 Nat. Genet., 53(3):354–366, March 2021.214 [23] Thomas Kluyver, Benjamin Ragan-Kelley, Fernando P´ erez, Brian Granger, Matthias Bus-215 sonnier, Jonathan Frederic, Kyle Kelley, Jessica Hamrick, Jason Grout, Sylvain Corlay,216 Paul Ivanov, Dami´ an Avila, Safia Abdalla, Carol Willing, and Jupyter development team.217 Jupyter notebooks - a publishing format for reproducible computational workflows. In218 Fernando Loizides and Birgit Scmidt, editors, Positioning and Power in Academic Pub-219 lishing: Players, Agents and Agendas , pages 87–90, Netherlands, 2016. IOS Press.220 [24] Siu Kwan Lam, Antoine Pitrou, and Stanley Seibert. Numba: a LLVM-based python JIT221 compiler. In Proceedings of the Second Workshop on the LL VM Compiler Infrastructure222 in HPC, New York, NY, USA, November 2015. ACM.223 [25] Ieva Rauluseviciute, Rafael Riudavets-Puig, Romain Blanc-Mathieu, Jaime A Castro-224 Mondragon, Katalin Ferenc, Vipin Kumar, Roza Berhanu Lemma, J´ er´ emy Lucas, Jeanne225 6 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint Ch` eneby, Damir Baranasic, Aziz Khan, Oriol Fornes, Sveinung Gundersen, Morten Jo-226 hansen, Eivind Hovig, Boris Lenhard, Albin Sandelin, Wyeth W Wasserman, Fran¸ cois227 Parcy, and Anthony Mathelier. JASPAR 2024: 20th anniversary of the open-access228 database of transcription factor binding profiles. Nucleic Acids Res. , 52(D1):D174–D182,229 January 2024.230 [26] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje. Learning important features231 through propagating activation differences. arXiv [cs.CV] , April 2017.232 [27] Scott Lundberg and Su-In Lee. A unified approach to interpreting model predictions.233 arXiv [cs.AI] , May 2017.234 7 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint Supplementary Note 1: Tomtom Algorithm235 Tomtom is a statistical test that converts similarity scores between two motifs into p-values. A236 key challenge in calculating these p-values is that these motifs are not assumed to be aligned,237 and so the similarity is derived as the maximum similarity across all potential ungapped238 alignments (excluding insertions and deletions). Consequently, the primary contribution of239 Tomtom is the algorithm for faithfully calculating null distributions in a manner that accounts240 for this maximum operation across all possible alignments. Importantly, Tomtom is not itself241 a similarity score and so, conceptually, can be used with any similarity or distance score S242 between individual positions, though the original implementation only allows users to choose243 between a few hard-coded ones. In practice, additively decomposable score functions are244 significantly faster to compute and negative Euclidean distance works well. Here, we will245 provide a description for how to recreate the original Tomtom implementation, with notes on246 what tomtom-lite changed to improve speed.247 Notationally, we will follow the original Tomtom work. Q ∈ Rwq,d is a single “query” of248 length wq and with an alphabet of size d. T ∈ RwT ,d is the concatenation of all n sequences in249 a “target database” across the length dimension, with Ti ∈ Rwi,d being the i-th entry in this250 database and wT = nP i=1 wi. Here, the target database simply means a collection of sequences251 that are being scanned against the query, potentially a motif database like JASPAR2024.252 Finally, to be consistent with the original Tomtom work we will refer to each of the vectors of253 size d as “columns”, meaning that there are wq columns in the query and wT columns in the254 target database.255 1.1 Integerized Similarity Matrix256 The first step of Tomtom is to calculate a matrix of similarity scores between the columns in the257 query and those in the target database. We calculate Θ ∈ Rwq,wT where each row represents258 one column in the query and each column represents one column in the target database and259 Θi,j = S(Qi, T j) where T j represents one column in the concatenation, as opposed to an260 entire target. As described by Tanaka et al., normalizing these scores accounts for differences261 in the number of unaligned columns when considering all potential alignments. This step262 is important to prevent single-column alignments on the edges that happen to match from263 scoring higher than an imperfect match between the core of the motifs. To do so, we calculate264 the median mi for each of the wq rows in Θ. In general, this normalization is robust to265 approximations of the median, and so tomtom-lite approximates this median for efficiency266 reasons (see Supplementary Note 2). Finally, these raw similarity scores are converted into267 integers ranging from 1 to t where t is the maximal similarity, which can be set by the user268 but defaults to 100. As a complete description, this procedure involves269 α = ⌊min i,j Θi,j − mi⌋ (1) where α is the floor of the minimum median-normalized value in Θ,270 β = max i,j Θi,j − mi (2) where β is the maximum median-normalized value (without taking the floor) in Θ, and271 ζ = ⌊ t β − α ⌋ (3) where ζ is the bin width. Given these values, we can directly convert the original (not median-272 normalized) similarity scores into integerized scores via273 1 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint Γi,j = ⌊(Θi,j − mi)ζ + 0.5⌋ + αζ (4) In the original Tomtom work by Gupta et al., Γ is used to represent both the initial matrix274 of continuous similarity scores and also the subsequent matrix of processed integer scores. For275 conceptual simplicity, we have separated them.276 Admittedly, these steps are more complicated than traditional methods for converting a277 known range of number into integer bins, and it is unclear the value that they add, but they278 are necessary for perfectly reproducing the original Tomtom implementation.279 1.2 Null distributions280 The next step is to calculate the null distributions of similarity scores to use in the statistical281 test. We use the plural because one distribution is needed for each target length from 1 to282 max i wi, because differing target lengths means a different number of potential alignments we283 are taking the maximum operation over. We begin this step by calculating the marginal PDF284 of scores for each query column f ∈ Rwq,t. More formally,285 fi,j = 1 wT wTX k=1 δ(Θi,k = j) (5) where δ is the Kronecker delta function which is 1 when the term is true and 0 otherwise. This286 corresponds to Eq. 2 in the original Tomtom paper and results in an f where the sum of each287 row is equal to 1 and the values are the fraction of entries in Γi that correspond to that integer288 value.289 We then use these marginal PDFs to calculate PDFs for each span of query columns, with290 the crucial assumption necessary for calculating a null distribution that these marginal PDFs291 are independent from each other. These span score PDFs encode the probabilities of observing292 each integer score when the target is aligned to only this span of columns in the query, assuming293 independence between the query columns. They will be important for quickly calculating the294 null distributions because they will be used in the simulated alignment procedure. Because Ti295 can be shorter than Q, we must consider all internal spans in addition to those that originate296 or end at the edges of Q. Mathematically, each span score PDF is297 P (yi:j = k) = tX l=1 P (yi:j−1 = k − l)P (xj = l) (6) where yi:j refers to the span PDF between columns i and j in the query and xj refers to298 marginal distribution over the j-th column. Notably, P (xj) is already calculated in fj and299 this equation can be solved recursively.300 Accordingly, in our implementation we define A ∈ Rwq,wq,twq as the set of these span score301 PDFs. Here, the first dimension is the start of the span, the second dimension is the end of302 the span (inclusive), the third dimension is the size of the maximum integer score possible in303 any span, and the values are the probability of that score given that span. The lower triangle304 formed by the first two dimensions will be unfilled because spans cannot be negative. Because305 each Ai,j is the probability distribution P (yi:j), summing this tensor along the last dimension306 will produce a matrix where the upper triangle is entirely 1s.307 Our first step in filling in A is to set the diagonals to be marginal distributions such that308 Ai,i = fi. Then, we proceed to fill out each row recursively as309 Ai,j,k = tX l=1 Ai,j−1,k−lfi,l (7) 2 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint Essentially, the probability of observing the integer score k is the probability of observing310 k − l when the span were one column shorter multiplied by the probability of taking a step of311 size l, summed over all l between 1 and t. This is recursive because the distribution at index312 j can be directly calculated from the one at index j − 1.313 Finally, we can use A to calculate our null distributions by simulating the alignment pro-314 cedure. Because the returned score is the maximum of scores across all alignments, we need315 to calculate the distribution of the maximum integer that would be emitted by all relevant316 PDFs in A. Note that we are not simply taking the maximum probability at each position317 in the set of PDFs, but rather are calculating a new distribution describing the probability318 of the maximum sampled value across a set of PDFs being equal to a certain value. Because319 calculating this given a large set of distributions can be tricky to implement, we instead choose320 to iteratively calculate the maximum of a pair of distributions until the entire set has been321 exhausted. Specifically, we use322 zi = xiyc i + yixc i − xiyi (8) where xc and yc are the CDFs of x and y respectively and z is the PDF of the maximums.323 Although applying Eq. 8 iteratively to a set of PDFs is straightforward, determining the324 set of PDFs to use requires some thought. Essentially, the number and composition of PDFs325 depends on the lengths of Q and Ti. There are two situations. When wq ≤ wi, this set contains326 all PDFs in the first row of A except for the last one, A0,0...wq−2, and all PDFs in the last327 column of A except the first one, A1...wq−1,wq−1 (the index A0,wq−1 being excluded is the same328 in both cases, and is the top right index in the matrix). Then, this top-right index A0,wq−1 is329 included a total of wi − wq + 1 times; one time for each alignment of the entire query against330 a target that is longer than it. Because we are taking the max of a set of PDFs, including the331 same PDF multiple times does have an effect. When wq > w i, we begin to use the internal332 spans. Specifically, if the canonical diagonal of a matrix is the 0-th diagonal, we use all PDFs333 in the wi − 1st diagonal of A, all the PDFs in the first row leading up to the that diagonal,334 and all the PDFs in the final column leading to the bottom right index. For example, when335 wq = 4 and wi = 2, we would use A0,0, A0,1, A1,2, A2,3, A3,3.336 The above process is repeated for each target length and results in a PDF. Because we337 want to quickly calculate p-values, we convert each PDF into 1 - CDF so that indexing into338 the distribution directly yields a p-value. Because our PDFs are over discrete integers, this339 can be achieved simply by taking the cumulative sum from the smallest to the largest integer340 in the distribution, and subtracting each value from 1.341 A critical point is that these null distributions depend only onwi and not on the composition342 of the target. In the original Tomtom implementation, this null distribution is recalculated for343 each pair of Q and Ti, despite null distributions being identical for any target in a given target344 database of the same length. In tomtom-lite, we precalculate null distributions for all lengths345 from 1 to max i wi and re-use them across targets. Because calculating these distributions takes346 a significant amount of time, this change is the primary speed improvement in tomtom-lite.347 1.3 Calculating p-values348 At this point, all that is left is to calculate the similarity scores for each possible ungapped349 alignment between Q and Ti, take the maximum, and convert this score into a p-value using350 the null distribution corresponding to the length ofTi. Fortunately, this process is easy because351 we have already calculated the integerized distances between all columns in the query and all352 columns in the target database. Considering the slice of Γ corresponding only to the columns353 of Ti, denoted γ (remember that Γ has wq rows and wT columns, spanning the entire target354 database, meaning γ has wq rows and wi columns), the similarity scores for the alignments can355 3 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint be easily calculated as the sum of the diagonals of γ. The maximal such diagonal-sum is then356 kept and converted to a p-value by indexing into the appropriate 1 - CDF. When considering357 reverse complements, the query and it’s reverse complement are both processed independently,358 the maximal score between the two is kept, and the p-value becomes 1 − (1 − min(p1, p2))2359 where p1 and p2 are the relevant p-values.360 1.4 Implementation Details for tomtom-lite361 Despite tomtom-lite being a Python function, we have demonstrated that it is significantly362 faster than Tomtom, which is implemented in C. There are several factors that contribute to363 tomtom-lite’s speed, and they can be divided broadly into those that are beneficial without364 reducing precision, and those that are beneficial but do cause a loss of precision. In practice,365 we have found that the loss of precision for the default settings is minimal, but may become366 meaningful in non-standard settings or when changing the defaults.367 1.4.1 Speed Improvements Without Precision Loss368 First, tomtom-lite is implemented using numba. Numba is a just-in-time (JIT) compiler that369 compiles functions into machine code the first time they are run. This offers C-like speed370 for numeric operations with Python-like syntax. Built-in caching functionality enables these371 compiled functions to be saved to disk, so the compilation cost only occurs the first time the372 function is run. Combined, this means that despite being available as a Python function,373 tomtom-lite would be as fast as a C implementation.374 Second, numba enables multi-threaded parallelism without significant code changes and this375 parallelism is done across the processing of queries. Importantly, this is truly multi-threading376 with shared memory and not simply multi-processing, where multiple instances of Python are377 started and communication costs are high. This is possible because the compiled code is not378 bound by the global interpreter lock (GIL) that restrains Python code. Simply by usingprange379 and using the parallel=True decorator appropriately, for loops can be converted into being380 parallel so long as the entire operation remains within numba (e.g., no calls to pure Python381 functions or globals) and calls within the loop do not depend on each other. This is satisfied in382 our implementation because each query can be processed independently. To avoid having each383 thread allocating small amounts of memory for each array needed for internal computation,384 tomtom-lite allocates a scratchpad before entering the loop whose first dimension is equal to385 the number of threads being used. Each thread uses the portion of the scratchpad allocated386 to its thread ID to avoid interfering with each other’s computations. The final output of each387 thread is written out to a pre-allocated results block.388 Third, through careful implementation, tomtom-lite dramatically reduce the number of389 cache misses when calculating alignment scores. At a high level, when memory is accessed on390 modern systems, the elements surrounding a desired index are usually also loaded. Because391 the cost of starting an I/O operation is significantly higher than continuing it or performing392 compute on those elements, implementations that can make use of these adjacent elements393 can be several times faster than those that do not, even though, conceptually, the same set of394 operations are being performed. A “cache miss” is when elements are loaded because they are395 adjacent to earlier elements, then discarded, then directly accessed later on, and represent a396 wasted opportunity to eliminate an I/O call.397 The potential for cache misses arises when summing the diagonals of γ. If we access an398 element and load the next elements in the row but then need to access the next element in399 the diagonal, we get an almost maximal number of cache misses because our operations are400 structured in a way to almost never make use of adjacent data. tomtom-lite eliminates these401 cache misses by, instead of proceeding along diagonals, proceeding along rows in the data and402 4 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint adding the element to an array of shape wq +wi −1. By ordering the operations in this manner,403 we observe a ∼2-4x speed improvement in the step at the cost of having to keep the scores404 from all alignments in memory and taking the max at the end, instead of being able to keep405 a running maximum across diagonals.406 1.4.2 Speed Improvements with Precision Loss407 In the process of turning Θ into Γ, the median of each row is subtracted from elements in408 that row. This has the effect of normalizing the score to account for the number of un-409 aligned columns. However, calculating the exact median is time-consuming as doing so takes410 O(wT log w T ) time, which can grow to be quite large in some applications. Consequently, we411 replace this exact median with an approximate median that takes only O(wT + m) time for a412 user-defined m, which is usually small with respect to wT . In this approximate formulation,413 the range between the minimum and maximum value (which is are calculated on-the-fly dur-414 ing the construction of Θ) is divided into m bins. After scanning through the wT items in415 Θi to count the number of instances of each, we can scan through the m bins until we have416 encountered 50% of elements. For additional precision, we can store not only the count but417 the sum of elements within each bin and so, at the end, return the average value of elements418 within the bin. Although our empirical results have shown that Tomtom is robust to the use419 of an approximate median, the returned p-values may be different.420 A final speed improvement involves reducing redundant calculations when calculating Θ.421 Essentially, many columns in T will be redundant due to the same character appearing in422 many motifs (e.g., a position in the PWM almost exclusively being an A, C, G, or T) and423 uninformative positions being fairly similar because they are almost the uniform distribution.424 Rather than re-calculating the similarity between the query column and these repeat columns425 in T , we can identify sets of similar columns in T as a preprocessing step and then calculate426 Θ over this set of non-redundant columns.427 However, a key challenge is identifying these redundant columns in T quickly. The most428 conceptually straightforward way is to calculate the similarity between all columns inT and set429 a threshold for clustering, but this would likely take longer than the original task of calculating430 the similarity between Q and T . Accordingly, we use a hashing algorithm to identify similar431 columns. Briefly, for each of the d elements in Tj, we divide the range between the minimum432 and maximum value into a user-defined number of bins b. Each column can then be assigned a433 hash index between 0 and bd where elements with the same hash are similar, without needing434 to know anything more than the values within the column. Note that the number of hashes435 will likely far exceed the number of examples, but that we do not need to allocate memory for436 all bd potential hashes.437 When using hashed columns, minor algorithmic modifications need to be made to the other438 steps. When calculating the marginal distributions f we need to modify the Kronecker delta439 function to not return 1 when the statement is true, but rather return a count corresponding440 to the number of items in the hash index being considered. Once f is calculated, the null441 distributions can be calculated without modification.442 A more complicated modification needs to be made when calculating the alignment scores443 while still minimizing cache misses. Specifically, instead of being able to extract γ directly444 as a series of contiguous columns, we have to look up the hash indexes for each column in445 Ti and iteratively load those columns. We proceed with a similar process to before, where446 the entirety of the column is used to update part of the vector of similarity scores. However,447 because we are operating on columns and the pre-loading mechanism works best on rows,448 we have to transpose the entire Γ matrix such that we are loading rows instead of columns.449 Because we are no longer loading a series of memory-adjacent rows (as there is no guarantee450 5 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint and it is actually quite unlikely that adjacent columns in Ti are adjacent in hash index) there451 is a small cost to doing it this way. However, despite slightly slowing down the calculation of452 similarity scores, this significantly speeds up the calculation of Γ to the point where it is no453 longer a compute bottleneck.454 6 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint Supplementary Note 2: Methods and Evaluations455 In this work, we performed several experiments to demonstrate that tomtom-lite produces456 similar results to the original Tomtom implementation and to time the two implementations.457 The code and Jupyter notebooks for these comparisons can be found in the memesuite-458 lite GitHub repository: https://github.com/jmschrei/memesuite-lite. These compar-459 isons were performed on a compute server with 76 Intel(R) Xeon(R) Gold 6138 CPU @460 2.00GHz cores. Many of these comparisons involved the JASPAR2024 motif database (called461 “JASPAR” for the rest of this note). We used the non-redundant core version, which can462 be be found at https://jaspar.elixir.no/download/data/2024/CORE/JASPAR2024_CORE_463 non-redundant_pfms_meme.txt464 In these comparisons, each call to the Tomtom command-line tool was as follows:465 tomtom q.meme t.meme -motif-pseudo 0 -thresh 1 -text -verbosity 1 > o.tomtom466 where q.meme is a MEME-formatted file containing the queries, t.meme is a MEME-formatted467 file containing the targets, and o.tomtom is a TSV containing the outputs of Tomtom. Our468 timings for Tomtom include only running this command, and do not include the creation of469 q.meme and t.meme or subsequent reading of o.Tomtom into memory. Similarly, when running470 our Python function we only time the running of the function itself.471 1.1.1 Comparing p-values between Tomtom and tomtom-lite472 Our first evaluation involved comparing the p-values produced by Tomtom and by tomtom-473 lite to ensure that they were similar. This comparison was done by using JASPAR for both474 the queries and the targets. After running and timing the two commands, we constructed a475 p-value similarity matrix that was square because the number of queries and targets was the476 same. Because self-comparisons are definitionally very similar and because Tomtom uses 64-bit477 operations by default and tomtom-lite uses 32-bit operations by default, we set the diagonal of478 this matrix to 1 to reduce the influence of machine precision. We report the Pearson correlation479 in two regimes: when the Tomtom p-values are above 1 e − 5 and when they are below that.480 1.1.2 Increasing Query Size481 Our next evaluation was to consider timings when running an increasing number of queries482 against JASPAR. These queries are derived from the JASPAR itself to ensure a biologically483 plausible comparison. Since our comparisons ultimately end up considering many more queries484 than are in the JASPAR motif, we end up cycling through it multiple times. Each query is485 independent from the others, and so running multiple identical queries through a comparison486 does not influence our timing (or statistical) results.487 Having shown that the benefits of multithreading diminish with increasing threads, we488 considered a second target database that is comprised entirely of random PWMs. These489 PWMs were generated using the following procedure: first, generate motifs of alphabet size 4490 with width 15 by independently randomly drawing from N(0, 1), second, exponentiate these491 values to get strictly positive values, and third, divide each column by the sum of the columns492 to convert each value into something resembling a probability. A benefit of this exponentiation493 operation is that it results in columns that are strongly biased towards one character at each494 position, which is similar to real motifs. We then follow the same procedure for generating495 queries as when we used JASPAR, except without the need to cycle through the target database496 multiple times because it is so big.497 7 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint 1.1.3 Increase Target Size498 We then turned to considering timings when the query size is held constant and the target499 size is fixed. Because multithreading was not shown to be particularly effecive when target500 databases were small, we considered only up to 8 threads and used only 20 queries. We501 constructed these target databases as the first n entries in JASPAR, looping back over the502 JASPAR once we have included all of the motifs, in a similar manner to how we constructed503 large query sets previously.504 2 ChromBPNet and Attributions505 Finally, we considered two evaluations involving feature attributions from a ChromBPNet506 model. We downloaded the ChromBPNet model for ATAC-seq signal in K562 from the EN-507 CODE Project (https://www.encodeproject.org, accession ENCSR868FGK). This acces-508 sion corresponds to an archive of five ChromBPNet models that have been trained and eval-509 uated on different folds, and we use only the one that is trained on fold 0. . We use the bp-510 netlite repository ( https://github.com/jmschrei/bpnet-lite) for loading this ChromBP-511 Net model into PyTorch and the tangermeme repository ( https://github.com/jmschrei/512 tangermeme) for subsequent analyses, specifically, the implementation of DeepLIFT/SHAP,513 the recursive seqlet caller, the seqlet annotation function that uses tomtom-lite’s Tomtom514 implementation, and the attribution plotting utilities.515 In the interactive example, we considered an enhancer of GATA2 whose coordinates are516 at chr3:128481848-128493691 on hg38. We extracted a 2,114 bp window centered at the517 middle of these coordinates and calculated feature attributions using DeepLIFT/SHAP with518 default settings (20 dinucleotide shuffles as the null distribution). To call seqlets, we used519 an unpublished recursive seqlet calling algorithm based on a statistical test and a recursive520 definition that each span within a seqlet must also be called as a seqlet. These seqlets were521 then annotated by taking the underlying discrete sequences and mapping them to JASPAR522 using tomtom-lite’s implementation of Tomtom.523 In our larger-scale example, we repeated the above procedure across all ATAC-seq peaks in524 K562. We began by downloading the peak calls from ENCODE (accession ENCFF558BLC).525 We then calculated attributions for each peak and called seqlets using the exact same procedure526 as above except that we only used 5 dinucleotide shuffles for DeepLIFT/SHAP. This was for527 computational efficiency, as calculating attributions with even only 5 shuffles took several528 hours. Finally, we used tomtom-lite to map the identified seqlets to JASPAR using 8 threads.529 Code reproducing these analyses can be found in the tutorials folder of the memesuite530 repository, in the notebook denoted as belonging to the application note.531 8 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted May 31, 2025. ; https://doi.org/10.1101/2025.05.27.656386doi: bioRxiv preprint