A run-length-compressed skiplist data structure for dynamic GBWTs supports time and space efficient pangenome operations over syncmers

preprint OA: closed CC-BY-4.0
📄 Open PDF Full text JSON View at publisher
AI-generated summary by claude@2026-07, 2026-07-15

This paper introduces a dynamic, run-length-compressed skiplist data structure for graph Burrows-Wheeler transforms, enabling efficient storage and search of pangenome paths and maximal exact matches over syncmer graphs.

One-sentence paraphrase of the abstract; not a substitute for reading it. No clinical advice. How this works

Abstract

Skiplists (Pugh, 1990) are probabilistic data structures over ordered lists supporting 𝒪 (log N ) insertion and search, which share many properties with balanced binary trees. Previously we introduced the graph Burrows-Wheeler transform (GBWT) to support efficient search over pangenome path sets, but current implementations are static and cumbersome to build and use. Here we introduce a doubly-linked skiplist variant over run-length-compressed BWTs that supports 𝒪 (log N ) rank and access operations, and a dynamic version of this that supports 𝒪 (log N + S ) rank and insert operations, where S is the number of symbols in the alphabet. We use these to store and search over paths through a syncmer graph built from Edgar’s closed syncmers, equivalent to a sparse de Bruijn graph. Code is available in rskip.[ch] within the syng package at github.com/richarddurbin/syng . This builds a 5.8 GB lossless GBWT representation of 92 full human genomes, single-threaded in 52 minutes, on top of a 4GB 63bp syncmer set built in 37 minutes. Arbitrarily long maximal exact matches (MEMs) can then be found as seeds for sequence matches to the graph at a search rate of approximately 1Gbp per 10 seconds per thread.
Full text 35,335 characters · extracted from oa-pdf · 3 sections · click to expand

Keywords

skip list · pangenome · GBWT · syncmer graph .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint Rskip data structures for pangenome GBWTs 1 1 Introduction The idea of pangenomes is that rather than use a fixed linear reference genome sequence we should represent the genetic variation known to be present in a species or population when analysing data from a new individual, or set of indi- viduals[23,2]. To do this we combine information from many independently as- sembled genomes, removing redundancy by collapsing together shared segments of sequence. A standard way to do this is in a graph, with shared sequence seg- ments represented by vertices, which are connected by edges that indicate what can follow what. Both the vertices and the edges are bi-directional, like highways or rail lines, which can be travelled in one direction or the other, with the se- quence for a vertex being reverse-complemented when traversed in the opposite direction. Most pangenome graph applications allow all routes across vertices that fol- low edges of the graph with equal likelihood, or perhaps weighted by the fre- quency of a vertex/edge in the set of reference genomes from which the graph was constructed. However there is longer range information in the reference genomes, based on shared haplotypes, that we want to use. In the railway analogy, the graph is the network but to get from A to B you need to take a series of trains that follow established paths through the network. This longer range information is what enables the powerful genotype imputation framework used for genetic association on millions of samples over the last fifteen years[10]. Recent methods for genotype imputation[20] have exploited the positional Burrows-Wheeler transform (PBWT)[3], which provides efficient storage of and matching to haplotype sequence panels. Our goal is to support a similarly ef- ficient and powerful framework for pangenome graphs, where sequences do not just vary at colinearly aligned binary variants. As a step in this direction, I introduce here a computational framework for efficient calculation over graph Burrows-Wheeler transforms (GBWTs), the corresponding data structure to the PBWT for pangenome graphs, and present a practical implementation. Below I present first the use of graph Burrows-Wheeler transforms (GBWTs) to store and search paths through pangenome graphs, and then the use of a vari- ant form of Pugh’s skip list[16], which I call Rskip, for efficient computation over pangenome GBWTs. Then I describe key features of myrskip implementation and its use in thesyng package, which builds pangenome grpahs in which the vertices represent overlappingkmers selected using Edgar’s closed syncmer crite- rion[4].Thisisfollowedbytheresultsofapplying syngtotheHumanPangenome

Reference

Consortium phase 1 release of 92 full human genomes[9]. 2 Theory Here we discuss using GBWTs to represent a set of paths on a pangenome graph and provide efficient search over these. We will use the term ‘vertex’ for the elemnents of the graph associated with sequences, in our case syncmers, and ‘node’ for nodes in the internal Rskip data structure that we use to manage the routing tables for the paths, which happens to also be a form of graph. .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint 2 Richard Durbin 2.1 Path traversal and search in a sparse graph via local GBWTs Fig. 1.LF mapping to follow a path entering as offset 1 in the path bundle from C and exiting as offset 2 in the path bundle to X. The local GBWTG provides the routing information, with the index inG being the sum of the incoming offset and the total count of preceding symbols in the input directoryD. All offsets and indexes are 0- based. Adapted from [21] with permission. Given a graph vertex with a set of paths running through it, and an ordering on each of the bundles of paths coming into the vertex from other vertices, it is natural to also list the vertices from which the paths come, and hence derive a global ordering of incoming paths (left side of Figure 1). Then a simple routing table of corresponding output nodes, as highlighted in purple in Figure 1), is all that is needed to specify the output vertex for each incoming path. Given that the paths within each outgoing bundle are ordered colinearly with incoming paths, as for the X’s and Y’s in Figure 1, adjacent paths in a bundle will share recent sequence, and by construction the routing table will correspond to the subsection of the Burrows-Wheeler (BWT) transform over sequences of vertices where the last symbol corresponded to the current vertex. This graph vertex BWT is called the Graph Burrows-Wheeler transform[21,5]. Because sharing past sequence is predictive of sharing future sequence and specifically the next symbol, BWTs in general and the GBWT in particular tend to contain runs of the same symbol. They are therefore efficiently compressed by run-length compression. The run-length compressed GBWT or rGBWT for Figure 1 would be((Y,1),(X,1),(Z,1),(X,2),(Y,2)). When following a path using a GBWT we need operationsocc(D, s) on the input directory D to give the total count of input paths from symbols listed before s, access(G, i) to give the symbol at positioni in the GBWT G and .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint Rskip data structures for pangenome GBWTs 3 rank(G, i, s) to give the number of timess appears before position i in G, as shown in Figure 1. The same operations support efficient matching. There are standard algorithms linear in match length in terms ofocc() and rank() op- erations to find the number of paths matching a substring or to find maximal exact matches (MEMs) of a longer search string (this is a little imprecise for MEMs, but let’s ignore that here). Note that this construction shows that for the update operations we do not require counts for all vertex symbols, only the ones incident to the current vertex, nor a global ordering of the symbols, just a local ordering of those that are incident. This is good because there are many millions of vertices, and different orderings may be preferable for different vertices. Typically BWTs are built statically by one process, and then an index struc- ture is built over them by a different process to support the required operations. Here we introduce a novel data structure that naturally supports dynamic build- ing of run-length compressed BWTs for efficient storage, traversal and search. 2.2 Skip lists for dynamic run-length-compressed arrays If we don’t want to move large memory blocks when we insert into a dynamic array the simplest structure would be a linked list, but this incurs linear time traversal cost for allaccess() and rank(). There are many succinct data struc- tures that can reduce these toO(log R) (where R is the number of runs) or for some operations even constant time[14], but typically these assume relatively small alphabets, whereas our alphabet size is potentially unbounded, reaching in practice at least tens of thousands due to the repetitive structure of the genome (see results). Here we take advantage of the elegant skip list data structure introduced by Pugh in 1990[16], which augments a linked list with additional layers generated by a random process in which each node has an upward parent with probability p (Figure 2A). Skip lists provide a lightweight alternative to balanced trees, sup- porting access() in expected O(log R) time as in Figure 2B, and alsoO(log R) insert() (and if wanted delete()) by keeping track of the last node visited at each level on a stack and updating their pointers and counts as necessary for a new column. Storing run lengths is natural in a skip list because the nodes at upper levels already store counts. To supportrank() once we are located at a node, we need to find the sum of counts of the same symbol in preceding nodes. We could count forward along the list how many are in succeeding nodes, then subtract from a total count that we maintain for each symbol, but this would be a linear time operation. Instead I add a second set ofright pointers that point to the next node at the same level with the same symbol, and counts for how many of that symbol are present up to that node. Effectively this is em- bedding a skip list for each symbol within the same set of columns as in the original primary skip list. We also need to addup pointers in reverse to the down pointers. See Figure 2C. Now we can find the number of copies of symbol s by moving up and right along theup and sRight pointers, going up whenever possible, and accumulating sCounts as we go. This is exactly reciprocal to the .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint 4 Richard Durbin Fig. 2.Skip list and Rskip data structures. A: a standard skip list over a run-length array, with columns up to (in this case) two levels. B: the access() operation is performed in expectedO(log R) time by starting at the top of the first column, moving right accumulating a partial sum, and dropping levels when the sum for the next right edge would exceed the target; only the partial sums coloured in red need to be calculated. C: Rskip structure with additional pointers and counts between nodes for the same symbol. initial access() search, so has expected timeO(log Rs) where there areRs runs for symbol s. Analogous to theinsert() operation for a standard skip list, it is possible to also insert new columns into the Rskip structure: once the rank is known the relevant node is found again usingsAcccess() on the symbol’s own skip list, and the stack of nodes up to that point can have theirsRight pointers and sCount values updated. Although the previous paragraph satisfies the requirement ofrank() for LF mappingtotraverseapath,wherewestartfromanodewiththerequiredsymbol, for general rank(s,i) for an arbitrary symbols we would need to find the next node with that symbol, which will require a linear search on the bottom level. The expected time for this isO(S) where S is the alphabet size, since if symbol s has frequencyfs then the expected time for this search is1/fs so the expected time for a random symbol isP s fs1/fs = P s 1 = S. Although this does not happen when traversing a path, it can happen when matching a new sequence, or inserting a new column. Fortunately this is worst case. First, entries with the same symbol tend to be clustered - indeed this is what facilitates run-length compression and indeed inserting into the structure by incrementing a run is by far the most common insertion, requiring no new nodes and just updating counts not pointers. Second, symbols that occur in only one run (singletons) can be special-cased and due to the skewed distribution of pangenome graph edges these are frequent. .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint Rskip data structures for pangenome GBWTs 5 In practice, I use two variants of the canonical Rskip structure described above. When building and requiring support forinsert() I add back pointers left and sLeft to simplify the update code at the cost of some additional memory. In particular, when doing this there is no need to maintain a special start column for each symbol. For static use while searching, it is possible to replace the count and sCount in the nodes in Figure 2C by partial sums up to that node, in a single linear sweep through the structure. Then therank() operation that followsaccess() during path following is a simple lookup, since the partialsSum of the symbol of a node is stored within the node. Furthermore, general rank() can be executed inO(log Rs + S) time by following thesRight and down pointers while looking at thesum entries then scanning left for a node with the required symbol. 3 Implementation The Rskip data structure and algorithms were implemented in C filerskip.c with public interface inrskip.h. These are available within the syng package at github.com/richarddurbin/syng, which implements pangenome construc- tion and, for now, basic search operations based on rskip and Edgar’s closed syncmers[4]. 3.1 Rskip implementation To avoid memory fragmentation and ensure cache locality the entire Rskip data structure is held in a single contiguous memory block made of an array ofmax nodes. In static search mode I use five 32-bit integers per node of typeFixed to record right, sRight and down pointers as offsets into the array, andsum, sSum partial sums. On the bottom level we set a top bit of sSum as a flag and use the down slot to hold the symbol identifiersym. The first node of the block is used as a header, recordingmax, the number of symbolsnSym and the startingnode.Thisisfollowedby nSymnodesformingadirectoryforthesymbols, which hold external symbol information and the total count for the symbol in the whole array. The nextnSym nodes following the directory are the tops of the first columns for each symbol. Column height for new nodes is sampled from a geometric distribution with mean height 1.6 by iterative sampling from a lightweight random number generator so long as the returned value is under 3/8, which is close to the optimum1/e. The dynamic Rskip variant uses eleven integers per node: six for bidirectional pointers right, left, sRight, sLeft, down, up, four for counts before and after before, count, sBefore, sCount andonefor sym.The before/sBefore counts are redundant with[left].count/[sLeft].sCount but simplify code at the left boundary which would otherwise require additional storage. There is again a header node and directory, in this case with the directory entries holding a pointer to the top of the left-most highest column for the relevant symbol. To .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint 6 Richard Durbin manage new node allocation, afreepointer is kept in the header, and array size is doubled when more space is required Because the Rskip structures are relatively heavy, I use a much simplified Linearnode array for small run lists, in both static and dynamic modes. Specif- ically these contain a maximum of 128 two-byte nodes, with nodes holding either a one-byte symbol identifier and count, or if the count is 255 then the subsequent nodeholds atwo-bytecount. Theheaderis nowtwonodes,and thedirectory uses four nodes per symbol. Again there is afree pointer in the header in dynamic mode. Implicit in the description above is that internal to the Rskip I use symbol values denoted sym from the range0..nSym-1. Calls to therskip package are made with a general integer, or in the case ofsyng with a general (symbol,offset) pair, and the correspondingsym is looked up/stored in the directory. Currently this lookup is done with a linear search, which in principle makes the time complexity O(S) where S is the local alphabet sizenSym. It would be possible to use a hash table or equivalent to make this constant time, but in practice the frequencies of symbols in pangenome graphs are very skewed (most genetic variants are rare). As a consequence, since the directory is built progressively as paths are added, high frequency symbols are at or close to the start of the directory list, meaning that the expected time for linear searches is short. Actual expected search times for a human pangenome are given in the results section. 3.2 syng implementation Since the focus of this paper is on the Rskip data structure and algorithms, only brief details of the syng implementation are given here. Syng pangenome graphs are bi-directional, with each vertex representing a syncmer and its reverse complement. Syncmers are fixed lengthkmers with the property that one of their terminal smers has less than or equal hash value comparedtoallinternal smers[4].Thisgivesthemaguaranteedwindowproperty, like minimizers to which they are related[19], with consecutive syncmers in a sequence overlapping by at leasts bases, though unlike minimizers syncmers are not context dependent: whether akmer is a syncmer is an intrinsic property of the kmer sequence. Hence syncmers cover any sequence (except for its ends), and a list of syncmers together with offsets between them, plus terminal ends of length <= k − s, are sufficient to reconstruct the sequence. For syng we require that k is odd (defaultk = 63, s = 8) and use a positive syncmer index value to indicate traversal in the direction that lexicographically comes before its reverse complement, and minus that value for the reverse-complement direction. Syng interconverts sequences and syncmer lists via a rapid scanning algorithm and a hash table. For each vertex syng keeps a 16-byte data structure with two sides,in and out, together with an 8-bit set of flags. In many cases there is only one edge for a side, termed asimple side; then the next syncmer index is stored directly in the side object, together with the offset to it and the edge count. Otherwise, a .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint Rskip data structures for pangenome GBWTs 7 pointer to an Rskip is stored, which maintains the list of edges in its directory together with a BWT supporting path reconstruction. To follow a path using LF-mapping with a BWT standardly requires a rank operation on the BWT and a sorted occurrence array. In our GBWT case only a small subset of the large total alphabet is relevant to each step, corresponding to the incoming and outgoing edges for each vertex in what is overall a very sparse graph. The rank operation is provided byrsRank() on the Rskip BWT for the out side, and we support the occurrence operation via our standard linear search of the directory for the in side, as when converting converting (syncmer,offset) pairs to rawsyms. During building this means that we need separate counts for incoming and outgoing edges for each side, but after inserting paths in both directions these become equal, and we only keep a single value in the static representation. 3.3 Storage in ONEcode files The syng graph structure is stored in a .1gbwt ONEcode file. Drawing on 30 years of bioinformatics experience, ONEcode provides a flexible, efficient framework for bioinformatics data storage, implemented with Gene Myers at github.com/thegenemyers/ONEcodeandusedforexamplebyPathPhynder[11] and FastGA[13]. It supports strongly-typed record-based binary files with built- in schemas, indices and data-specific compression (separate codec adaptively trained per record type), and a simple dependency-free interface (single C source file and header), thread support and standard text interconvertibility. Thegbwt ONEcode schema support serialisation of the vertex, edge and path structure stored in syng vertex and Rskip structures, effectively equivalent to the infor- mation stored in a GFA file[1]. Although syng operates on syncmer graphs, a .1gbwt file can store an arbitrary bidirectional sequence graph with paths over it. In parallel, syng uses a .1khash file type to store the syncmer sequences, enabling fast reconstruction of its syncmer data structures. 4 Results I demonstrate the performance of the Rskip data structure using 92 human genomes comprising the Human Pangenome Reference Consortium (HPRC) re- lease 1[9] minus the two HG002 genomes which will be used as test sequences for alignment. These contain 37,269 sequences totalling 277.4 Gbp, which generated 234 billion instances of 193 million syncmers, with average coverage 120.8, in 37 minutes starting from FASTA sequences. The syncmer hash table and index take 10.4GB in memory, but can be stored in 4.0GB on disk as a .1khash file. It took syng 52 minutes on a Linux server, single-threaded, to build a full bi-directional GBWT from the syncmer lists, using 15.7GB max memory. Time to add genomes increased from ∼22 seconds per genome to ∼40 seconds per genome (Figure 3). The resulting graph contained 339.8 million simple vertex sides connected to only one edge, and 46.2 million with multiple edges that .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint 8 Richard Durbin required Rskip objects. Of these 46 million were stored asLinear arrays with on average 2.4 symbols and 5.4 runs taking 2.0 GB memory, and 175 thousand as Dynamic arrays with on average 53.5 symbols and 410 runs taking 7.8 GB (9.8GB total). The total BWT length (sum of runs) was 4.2 billion forLinear nodesand1.1billionfor Dynamicnodes,withanaveragerunlengthof16.4across both types. In total there were 229 million edges, but the incidence distribution has a very long tail due to complex vertices from repetitive regions such as centromeres: the vertex side with the most edges had 13,506 and the maximum number of runs for one edge was 191,212 with maximum total count 1,274,525. There were 151,991 sides with more than 20 edges. The expected directory list search length was 6.3. On disk, the resulting .1gbwt file used 5.8GB. Fig. 3.Construction time pergenome and cumulative memory used when building a syng pangenome from 92 HPRCv1 genome sequences For initial mapping experiments I used the Pacific Biosciences HiFi readset from individual HG002 available froms3-us-west-2.amazonaws.com/human- pangenomics/T2T/HG002/assemblies/polishing/HG002/v1.1/mapping/ hifi_revio_pbmay24/hg002v1.1_hifi_revio_pbmay24.bam,comprising12.8 millionreadstotalling205Gbp(averagelength16.0kb).Instaticmodeforsearch- ing, the Rskip data structures takes less memory, 1.4GB forLinear and 2.6GB for Fixed (4.0GB total), both because there is no free space for expansion and because Fixed Rskip nodes are smaller thanDynamic nodes. With eight threads it takes 14.0 seconds to load the syncmer table and 36.0 seconds to load the ver- tices, edges and GBWT data using 8 threads. A single forward pass search of the 205Gbp took 468 seconds or 2.3 seconds per Gbp, coming to 8.6 minutes total .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint Rskip data structures for pangenome GBWTs 9 for the search including loading tables. Parallelisation performance did not scale with thread number, presumably because of cache contention accessing both the syncmer and GBWT tables: single-threaded times were 125.6 seconds to load the graph (3.5x) and 1890 seconds to search (9.2 seconds per Gb, 4x). This scan found 204 million maximal exact matches (MEMs) of average length 1304bp, with only 249 reads failing to find matches. MEMs are termi- nated either by a syncmer not found in the pangenome reference, which may be caused either by a sequencing error or a true genetic variant, or by ancestral re- combination generating a new haplotype not present in any reference sequence. In this case we believe the primary reason is sequencing error. We know that the PacBio HiFi sequencing error rate is around 0.1%, and that most errors are in homopolymer run lengths. If we homopolymer compress both the references and the reads, resulting in a 30% length reduction, then average match length increases to 4421bp compressed, equivalent to approximately 6300bp in original bases. 5 Discussion I have presented a GBWT computational model and implementation built on O(logR)core operations, which operates efficiently at the scale of a hundred human genomes. Empirical time complexity growth as seen in Figure 3 depends on many factors beyond run lengths, but is increasing relatively slowly and apparently sublinearly, so it is reasonable to expect this implementation to scale to the thousand haplotype level envisioned in current pangenome projects. Notably,syngconstruction is much faster than the time taken by Minigraph- Cactus[7] or PGGB[6] to construct the HPRCv1vg graphs[9]. That said, the syng graphs are fundamentally different in that the Minigraph-Cactus and pggb graphs are based on multiple sequence alignments that aim to identify large scalecolinearchromosomalhomologywithfewifanycycles,whereasthesyncmer graph built bysyngcontains many cycles through repetitive sequences. Similarly the GBWT described here is fundamentally different to the one implemented in the vg package that supports search using giraffe[22]. The DNA sequences of the vertices invg graphs are not unique, and so thegiraffe GBWT must index the {A,C,G,T} strings, rather than sequences of graph vertices, which in our case are syncmers each with their own distinct sequence so indexable by sequence. Our construction is much closer to a de Bruijn graph as used in sequence assemblers[15,8], but only using a covering set of sparsekmers as in MBG[17] and Verkko[18]. Indeed, the syncmer code insyng (though notrskip) hasalreadybeenusedinthe syncasmassemblerthatformspartoftheOraganelle Assembly TookKitoatk[24]. The MEM finding results presented above are only an initial step towards sequence searching. I note that the current code does not even report all MEMs, but this can be fixed given that our BWT is bidirectional by iterative for- ward/backward matching, and in any case enough matches were found to each read to act as seeds for more complete alignment that accepts sequencing er- .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint 10 Richard Durbin rors and mutations, perhaps using principles from Myers’ wave alignment algo- rithm[12]. In the long run, as indicated in the introduction, the goal is to support impu- tation of genome sequences from partial data such as low coverage or short read data sets together with a pangenome reference panel. This will require pulling out the most likely pair of haplotype sequences through the graph, based on coverage of read alignments together with the longer range haplotype structure represented by the paths stored in the GBWT. Acknowledgments.I thank Chenxi Zhou and Gene Myers for ideas and comments relating to this work, and more generally many other colleagues over the years for stim- ulating relevant discussions. This work was supported by Wellcome Discovery award 317408/Z/24/Z and a grant on Quantum Pangenomics from Wellcome Leap’s Q4Bio program. Disclosure of Interests.R.D. is a scientific advisory board member of Dovetail Inc. with a small financial interest in its holding company EdenRoc Inc.

References

1. Graphical fragment assembly (GFA) format specification, https://github.com/GFA-spec/GFA-spec 2. Bao, Z., Weigel, D.: Complexity welcome: Pangenome graphs for comprehensive population genomics. Quant. Plant Biol.6(e43), e43 (27 Oct 2025) 3. Durbin, R.: Efficient haplotype matching and storage using the positional burrows- wheeler transform (PBWT). Bioinformatics30(9), 1266–1272 (1 May 2014) 4. Edgar, R.: Syncmers are more sensitive than minimizers for selecting conserved k-mers in biological sequences. PeerJ9, e10805 (5 Feb 2021) 5. Gagie, T., Manzini, G., Sirén, J.: Wheeler graphs: A framework for BWT-based data structures. Theor. Comput. Sci.698, 67–78 (25 Oct 2017) 6. Garrison, E., Guarracino, A., Heumos, S., Villani, F., Bao, Z., Tattini, L., Hag- mann, J., Vorbrugg, S., Marco-Sola, S., Kubica, C., Ashbrook, D.G., Thorell, K., Rusholme-Pilcher, R.L., Liti, G., Rudbeck, E., Golicz, A.A., Nahnsen, S., Yang, Z., Mwaniki, M.N., Nobrega, F.L., Wu, Y., Chen, H., de Ligt, J., Sudmant, P.H., Huang, S., Weigel, D., Soranzo, N., Colonna, V., Williams, R.W., Prins, P.: Build- ing pangenome graphs. Nat. Methods pp. 1–5 (21 Oct 2024) 7. Hickey, G., Monlong, J., Ebler, J., Novak, A.M., Eizenga, J.M., Gao, Y., Human Pangenome Reference Consortium, Marschall, T., Li, H., Paten, B.: Pangenome graph construction from genome alignments with minigraph-cactus. Nat. Biotech- nol. 42(4), 663–673 (Apr 2024) 8. Li, H., Durbin, R.: Genome assembly in the telomere-to-telomere era. Nat. Rev. Genet. (22 Apr 2024) 9. Liao, W.W., Asri, M., Ebler, J., others: A draft human pangenome reference. Na- ture 617(7960), 312–324 (May 2023) 10. Marchini, J., Howie, B.: Genotype imputation for genome-wide association studies. Nat. Rev. Genet.11(7), 499–511 (2 Jul 2010) 11. Martiniano, R., De Sanctis, B., Hallast, P., Durbin, R.: Placing ancient DNA se- quences into reference phylogenies. Mol. Biol. Evol.39(2) (3 Feb 2022) .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint Rskip data structures for pangenome GBWTs 11 12. Myers, E.W.: AnO(ND) difference algorithm and its variations. Algorithmica1(1), 251–266 (1 Nov 1986) 13. Myers, G., Durbin, R., Zhou, C.: FastGA: fast genome alignment. Bioinform. Adv. 5(1), vbaf238 (6 Oct 2025) 14. Navarro, G.: Compact Data Structures: A Practical Approach. Cambridge Univer- sity Press, Cambridge, England (5 Sep 2016) 15. Pevzner, P.A., Tang, H., Waterman, M.S.: An eulerian path approach to DNA fragment assembly. Proc. Natl. Acad. Sci. U. S. A.98(17), 9748–9753 (14 Aug 2001) 16. Pugh, W.: Skip lists: a probabilistic alternative to balanced trees. Commun. ACM 33(6), 668–676 (1 Jun 1990) 17. Rautiainen, M., Marschall, T.: MBG: Minimizer-based sparse de bruijn graph con- struction. Bioinformatics37(16), 2476–2478 (15 Aug 2021) 18. Rautiainen, M., Nurk, S., Walenz, B.P., Logsdon, G.A., Porubsky, D., Rhie, A., Eichler, E.E., Phillippy, A.M., Koren, S.: Telomere-to-telomere assembly of diploid chromosomes with verkko. Nat. Biotechnol. (16 Feb 2023) 19. Roberts, M., Hayes, W., Hunt, B.R., Mount, S.M., Yorke, J.A.: Reducing storage requirements for biological sequence comparison. Bioinformatics20(18), 3363–3369 (12 Dec 2004) 20. Rubinacci, S., Delaneau, O., Marchini, J.: Genotype imputation using the posi- tional burrows wheeler transform. PLoS Genet.16(11), e1009049 (16 Nov 2020) 21. Sirén, J., Garrison, E., Novak, A.M., Paten, B., Durbin, R.: Haplotype-aware graph indexes. Bioinformatics36(2), 400–407 (15 Jan 2020) 22. Sirén, J., Monlong, J., Chang, X., Novak, A.M., Eizenga, J.M., Markello, C., Sibbe- sen, J.A., Hickey, G., Chang, P.C., Carroll, A., Gupta, N., Gabriel, S., Blackwell, T.W., Ratan, A., Taylor, K.D., Rich, S.S., Rotter, J.I., Haussler, D., Garrison, E., Paten, B.: Pangenomics enables genotyping of known structural variants in 5202 diverse genomes. Science374(6574), abg8871 (17 Dec 2021) 23. Taylor, D.J., Eizenga, J.M., Li, Q., Das, A., Jenike, K.M., Kenny, E.E., Miga, K.H., Monlong, J., McCoy, R.C., Paten, B., Schatz, M.C.: Beyond the human genome project: The age of complete human genome sequences and pangenome references. Annu. Rev. Genomics Hum. Genet.25(1), 77–104 (27 Aug 2024) 24. Zhou, C., Brown, M., Blaxter, M., Darwin Tree of Life Project Consortium, Mc- Carthy, S.A., Durbin, R.: Oatk: a de novo assembly tool for complex plant organelle genomes. Genome Biol.26(1), 235 (7 Aug 2025) .CC-BY 4.0 International licenseavailable under a (which 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 preprintthis version posted March 29, 2026. ; https://doi.org/10.64898/2026.03.26.714584doi: bioRxiv preprint

Text is read by the "Ask this paper" AI Q&A widget below. Extraction quality varies by source — PMC NXML preserves structure cleanly, OA-HTML may include some navigation residue, and OA-PDF can have broken hyphenation. The publisher copy (via DOI) is the canonical version.

My notes (saved in your browser only)

Ask this paper AI returns verbatim quotes from the full text · source: oa-pdf

Answers must be backed by verbatim quotes from this paper's full text. Hallucinated quotes are dropped automatically; if no verbatim passage answers the question, we say so. How this works

Citation neighborhood (no data yet)

We don't have any in-corpus citations linked to this paper yet. This is a recent paper (2026) — citers typically take a year or two to land, and the OpenAlex reference graph may still be filling in.

Source provenance

europepmc
last seen: 2026-05-20T01:45:00.602351+00:00
unpaywall
last seen: 2026-05-22T02:00:06.705733+00:00
License: CC-BY-4.0