Billi: Provably Accurate and Scalable Bubble Detection in Pangenome Graphs

Billi: Provably Accurate and Scalable Bubble Detection in Pangenome Graphs | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Billi: Provably Accurate and Scalable Bubble Detection in Pangenome Graphs View ORCID Profile Shreeharsha G Bhat , View ORCID Profile Daanish Mahajan , View ORCID Profile Chirag Jain doi: https://doi.org/10.1101/2025.11.21.689636 Shreeharsha G Bhat 1 Department of Computational and Data Sciences, Indian Institute of Science , KA 560012, India 2 Department of Biotechnology, Indian Institute of Technology Madras , TN 600036, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Shreeharsha G Bhat Daanish Mahajan 3 Department of Interdisciplinary Mathematical Sciences, Indian Institute of Science , KA 560012, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Daanish Mahajan Chirag Jain 1 Department of Computational and Data Sciences, Indian Institute of Science , KA 560012, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Chirag Jain For correspondence: chirag{at}iisc.ac.in Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract A key application of pangenome graphs is the characterization of small and large genomic variants represented as bubbles within the graph. Although bubbles have been extensively studied in directed graphs in the context of genome assembly, there remains a need for a rigorous definition and systematic analysis of bubbles in bidirected graphs, which is the predominant data structure used to represent pangenomes. We show that existing bubble definitions for bidirected graphs do not fully meet the requirements for representing genetic sites and alleles; in particular, overlapping bubbles may not exhibit strict nesting. To address this, we introduce a new sub-graph abstraction called panbubble and prove that it satisfies the desired structural properties for variant characterization. We then present an exact algorithm with runtime 𝒪 (| V | 2 (| V | + | E |)) for detecting all panbubbles in a bidirected graph G = ( V, E ). In addition, we propose a heuristic algorithm that produces identical output as the exact algorithm in practice and scales to large graphs, including both the first and second releases of the Human Pangenome Reference Consortium (HPRC). We implemented our algorithms in the tool Billi ( github.com/at-cg/billi ). On our largest dataset, Billi is more than 15× faster and uses over 5× less memory than VG. 1 Introduction Pangenome graphs provide a flexible framework for representing complex structural variations, including cases where smaller variations are nested within larger ones [ 1 , 3 , 14 , 26 , 27 ]. Unlike conventional variant calling methods, variant discovery in pangenome graphs does not use a single reference genome coordinate system. A key computational challenge in this setting is the identification of bubbles , which are substructures in the graph that likely correspond to genomic variation. Informally, a bubble is a region where multiple paths diverge from a common source vertex and subsequently reconverge at a common sink vertex. The concept of bubbles and efficient algorithms for their detection have been extensively studied in the context of genome assembly [ 2 , 11 , 23 , 28 ]; however, these studies typically assume a directed graph representation. In contrast, pangenome graph construction tools [ 10 , 12 , 19 , 20 ] employ bidirected graphs [ 9 ] to represent double-stranded DNA. Previous efforts [ 8 , 20 , 22 , 24 ] to generalize the concept of bubbles to bidirected graphs build upon the notion of the superbubble , first introduced and formally defined in [ 23 ]. Superbubble is an important class of subgraphs in directed graphs, originally motivated by the goal of simplifying genome assembly graphs. This concept enables the identification of substructures within an assembly graph that arise from sequencing errors, heterozygous variants, and near-identical repeats. In their work, Onodera et al . [ 23 ] precisely defined the conditions under which a pair of distinct vertices s and t can serve as the source and sink of a superbubble, respectively. Specifically, the enclosed superbubble subgraph B must be acyclic; there must be no edge from a vertex outside B to any vertex in B \ { s }; and no edge from a vertex in B \ { t } to a vertex outside B . Superbubbles exhibit elegant structural properties. For instance, the subgraphs of two superbubbles are either disjoint or strictly nested, and every vertex within a superbubble lies on some source-sink path. In pangenome graphs, cyclic subgraphs are of particular interest because genomic rearrangements such as inversions and duplications give rise to paths that loop back. Several generalizations of the superbubble concept have been proposed for bidirected graphs , including snarl [ 24 ], bibubble [ 20 ], and flubble [ 22 ]. These formulations relax the acyclicity constraint, enabling the representation of more complex graph structures. However, they do not fully preserve all properties of superbubbles. For instance, snarls and flubbles may contain vertices that are not part of any source-sink walk (Supplementary Figure S2). Moreover, it is possible that a pair of snarls, bibubbles, or flubbles in a bidirected graph overlap without being strictly nested (Supplementary Figure S3). This can introduce ambiguity in delineating regions corresponding to genomic variation. Although overlapping subgraphs can be filtered using some heuristic to enforce disjointness, this approach risks missing genuine variation. In this paper, we extend the superbubble framework [ 23 ] to bidirected graphs while preserving key structural properties and accommodating cyclic subgraphs. Our main contributions are as follows: We introduce a subgraph abstraction called panbubble . This is inspired by the superbubble concept but adapted to bidirected graphs. We also define hairpin , a related substructure that corresponds to inversions [ 22 ]. We give formal proofs showing that panbubbles and hairpins exhibit the intended nesting behavior. We present an 𝒪 (| V | 2 (| V | + | E |))-time algorithm for finding all panbubbles and hairpins. We also propose a faster 𝒪 (| V |(| V | + | E |))-time heuristic algorithm that yields identical results as our exact algorithm on the test datasets. We implemented both algorithms in our tool Billi . Billi processes the largest publicly available human pangenome graph, containing over a hundred million vertices, in under 75 minutes and with 90 GB of memory. In comparison, computing snarls using VG on this graph required approximately 20 hours and 500 GB of memory. We compare the identified panbubbles with snarls and bibubbles on several pangenome graphs, and observe strong overall agreement. Through Bandage [ 30 ] visualizations, we show that the observed differences arise because prior bubble definitions do not comply with all conditions of panbubble. 2 Notations and problem statement Biedged graph Bidirected graphs are widely used in genomics because they offer a convenient way to represent adjacency relationships between DNA sequences, while allowing traversal of a sequence in both forward and reverse orientations. In this study, we adopt the biedged graph representation of bidirected graphs, following the approach of previous works [ 20 , 24 ]. A biedged graph G b = ( V b , E b , f b ) is an undirected, colored multigraph subject to the following three constraints. Function f b : E b → { gray, black } assigns either gray or black color to each edge. Every vertex of the graph is incident with exactly one black edge. Between any two vertices (which may be identical), there can be at most one gray edge. We call two vertices u, v ∈ V b ( u ≠ v ) adjacent if they are connected by at least one edge of any color. With the above definition, observe that adjacent vertices u and v can have either (i) one gray edge, (ii) one black edge, or (iii) one black and one gray edge between them ( Figures 1 and 2A ). For every vertex v ∈ V b , we use to denote the neighbor vertex of v connected using a black edge. We denote the black edge incident to v as e v . Clearly, for all v ∈ V b . In genome sequence analysis applications where the biedged (or equivalently, bidirected) graph representation is used, each black edge is typically annotated with a DNA sequence. In any walk that traverses a black edge e v , the order in which vertices v and are visited determines the orientation (forward or reverse complement) in which the sequence label of e v is to be read. In this work, we omit these sequence labels and focus solely on the underlying graph topology. Download figure Open in new tab Figure 1: Hairpin ⟨ v 2 ⟩ is highlighted in a blue box. Panbubble ⟨ v 14 , v 23 ⟩ is highlighted in a red box. The blue box encloses all vertices in the set Y ( v 2 ) \ { v 2 }. The red box encloses all vertices in the set B ( v 14 , v 23 ). Note that vertex pairs ( v 0 , v 7 ) and ( v 6 , v 15 ) do not qualify as panbubble entrances because vertex v 2 violates the no-hairpin condition and vertex v 8 violates the contiguity condition, respectively. Download figure Open in new tab Figure 2: (A) Compact biedged graph G b . (B) The corresponding undirected multigraph G U . (C) Partitioning of edges in G U into a set of tree edges and a set of back edges using DFS starting from vertex v src . The notion of a walk in a biedged graph deviates from the standard definition of walks in undirected graphs. A sequence of vertices ( v 0 , v 1 ,…, v n ) with n ≥ 1 is a walk in a biedged graph if (a) for all 1 ≤ i ≤ n , there exists an edge between v i −1 and v i , (b) the walk starts and ends using black edges, that is, and , and (c) the entire walk alternates between black and gray edges. For example, vertex sequence ( v 0 , v 1 , v 2 , v 3 , v 4 , v 5 ) in Figure 1 forms a valid walk. On the other hand, the vertex sequence ( v 0 , v 1 , v 6 , v 2 , v 3 ) is not a walk because it contains two consecutive gray edges. Note that if vertex v appears in a walk, then its neighbor vertex must also appear in that walk. Due to symmetry, ( v n , v n −1 ,…, v 1 , v 0 ) is a walk if ( v 0 , v 1 ,…, v n −1 , v n ) is a walk. Considering any u, v ∈ V b , we say that v is reachable from u if v appears in some walk starting from u . Note that the set of vertices reachable from u always includes u and û because ( u, û ) is a valid walk. We say that a walk ( v 0 , v 1 ,…, v n −1 , v n ) passes through v i for all 1 ≤ i ≤ n — 1. Given any three vertices u, v, w ∈ V b , v is said to be reachable from u without passing through w if v appears in some walk starting from u and that walk does not pass through w . We call a vertex not incident with a gray edge a tip . Two vertex-disjoint subgraphs of a biedged graph are said to be disconnected in the usual sense, that is, if there exists neither black edge nor gray edge that has one endpoint in the first subgraph and the other endpoint in the second subgraph. Unless a biedged graph has two or more disconnected subgraphs, it is said to be connected . A walk ( v 0 , v 1 ,…, v n −1 , v n ) with n ≥ 3 is called a closed walk if v n −1 = v 0 and v n = v 1 . Similarly, a walk ( v 0 , v 1 ,…, v n −1 , v n ) with n ≥ 3 is an inverted closed walk if v n −1 = v 1 and v n = v 0 . For example, in Figure 1 , the vertex sequence ( v 8 , v 9 , v 12 , v 13 , v 11 , v 10 , v 8 , v 9 ) is a closed walk, whereas ( v 2 , v 3 , v 4 , v 5 , v 3 , v 2 ) is an inverted closed walk. Compact biedged graph The process of compacting long non-branching paths into single vertices is a common data reduction step in graph-based sequence analysis. In the context of a biedged graph G b = ( V b , E b , f b ), a compaction operation is defined as follows. Consider two vertices u, v ∈ V b with u ≠ v and . u and v are said to be compactable if (i) they are connected to each other by a gray edge, (ii) the number of gray edges incident to u is one, and (iii) the number of gray edges incident to v is one. The compaction operation on a compactable vertex pair ( u, v ) adds a black edge between û and , and then removes u and v from the graph along with all the edges incident to u or v . The resulting graph is a biedged graph with two fewer edges and two fewer vertices. A compact biedged graph is obtained by applying compaction operations as many times as possible in any order on the given biedged graph. We next introduce the concepts of hairpin and panbubble , which will be useful to characterize genomic variation in pangenome graphs. Definition (Hairpin) Given a compact biedged graph G b = ( V b , E b , f b ), let the set of vertices reachable from vertex s ∈ V b without passing through s be denoted by Y ( s ). A hairpin is said to exist at vertex s if (i) every vertex v ∈ Y ( s ) appears in an inverted closed walk starting from s , (ii) removal of all the edges connecting s and ŝ separates the subgraph induced by Y ( s ) \ { s } from the rest of the graph, and (iii) no vertex in Y ( s ) other than s satisfies these conditions. The separable subgraph induced by Y ( s )\ { s } is called a hairpin . For any vertex s that satisfies the above conditions, we refer to s and e s as the entrance vertex and the entrance edge of the hairpin, respectively. The hairpin is denoted by ⟨ s ⟩. Definition (Panbubble) A panbubble exists between two vertices s and t ( s ≠ t and ) in a compact biedged graph G b = ( V b , E b , f b ) if all the following conditions hold: matching : Let U ( s, t ) denote the set of vertices reachable from s without passing through s or t . Analogously, define U ( t, s ) as the set of vertices reachable from t without passing through t or s . The matching condition requires that U ( s, t ) = U ( t, s ). separable : Define B ( s, t ) as U ( s, t ) \ { s, t }. If black edges e s and e t are removed, the subgraph induced by B ( s, t ) is disconnected from the rest of the graph. contiguity : Every vertex v ∈ U ( s, t ) appears on at least one walk from s to t . no-hairpin : No vertex in is an entrance vertex of a hairpin. minimality : No vertex in U ( s, t ) other than t forms a pair with s and satisfies the above criteria. Similarly, no vertex in U ( s, t ) other than s forms a pair with t and satisfies the above criteria. The separable subgraph induced by B ( s, t ) is called a panbubble . For any set of two vertices { s, t } that satisfies these criteria, we denote the panbubble as ⟨ s, t ⟩ (or equivalently ⟨ t, s ⟩). In a panbubble ⟨ s, t ⟩, s and t are called entrance vertices , and the edges e s and e t are called entrance edges . Note that the entrance vertices s, t and the entrance edges e s , e t do not belong to the panbubble itself, i.e., they are not part of the subgraph induced by B ( s, t ). See Figure 1 for an example. In this paper, we seek to solve the following problem: Problem statement Given a compact biedged graph G b = ( V b , E b , f b ), (i) list the entrance vertex pairs of every panbubble in G b , and (ii) list the entrance vertices of every hairpin in G b . 3 Properties of panbubbles and hairpins In this section, we outline key properties of panbubbles and hairpins, including their possible nesting relationships. These properties will also guide the design of our detection algorithms. For brevity, proofs of all claims in this paper are provided in the Supplementary Document. Lemma 1. In a compact biedged graph G b = ( V b , E b , f b ), the number of distinct panbubbles can be at most | V b / 2|. The number of distinct hairpins can be at most | V b |. A valid walk in a biedged graph cannot traverse two gray edges consecutively, which limits the applicability of standard graph-theoretic techniques. To address this limitation, we construct an auxiliary undirected multigraph derived from G b , in which edge colors are ignored and the walk constraints are relaxed. Construction of an auxiliary undirected multigraph We define our auxiliary undirected multigraph G U = ( V U , E U ) as follows. We assume that G b is connected (if not, each connected component can be handled independently). First, we copy all the edges and the vertices of G b into G U . Note that tips in G b correspond to vertices that have degree one in G U . We assume that G b contains at least one tip * , and hence, G U has at least one degree-one vertex. Next, we introduce one more vertex v src in G U , and we connect v src to all degree-one vertices ( Figures 2A, 2B ). For every vertex v of G b , we denote the corresponding vertex in G U by . Vertex in G b corresponds to in G U . The edge of G U corresponding to a black edge e v in the biedged graph is denoted by . Next, for every vertex v of G b , we ensure that the first vertex in the adjacency list of vertex in G U is . We do this to guarantee that every edge in G U derived from a black edge in G b appears as a tree edge, rather than a back edge, in any depth-first traversal of G U . Unlike biedged graphs, a walk in G U is defined in the usual sense, that is, a sequence of vertices ( v 0 , e 0 , v 1 , e 1 , v 2 , …, e n −1 , v n ) with n ≥ 1 is a walk in G U if edge e i connects the vertices v i and v i +1 for all 0 ≤ i ≤ n —1. Walk ( v 0 , e 0 , v 1 , e 1 , v 2 ,…, e n −1 , v n ) passes through vertices v 1 , v 2 ,…, v n −1 . A vertex is reachable from v ∈ V U if it appears in some walk starting from v . Path and cycle in G U are also defined in the usual sense. A path is a walk whose vertices are distinct. A cycle in G U is a closed path that starts and ends at the same vertex. We do not consider a self-loop as a cycle. Let be the spanning tree obtained by a depth-first traversal of G U starting from v src . Through this traversal, the edge set E U is partitioned into a set of tree edges and a set of back edges ( Figure 2C ). Tree introduces a partial order ≺ over the vertices in G U . For v 1 , v 2 ∈ V U , v 1 ≺ v 2 if and only if (i) v 1 ≠ v 2 , (ii) v 1 and v 2 belong to a root-leaf path in , and (iii) v 1 appears before v 2 in that path from root to leaf. Similarly, for e 1 , e 2 ∈ E U , e 1 ≺ e 2 if and only if (i) e 1 ≠ e 2 , (ii) e 1 and e 2 belong to a root-leaf path in , and (iii) e 1 is traversed before e 2 in that path from root to leaf. A vertex w in V U is called an ancestor of a tree edge e = { u, v } with u ≺ v if either w = u or w ≺ u . Similarly, w is a descendant of edge e if either w = v or v ≺ w . Next, we define the notion of bracket set and cycle equivalence. Definition (Bracket Set) A bracket of a tree edge e in G U is a back edge connecting a descendant of e to an ancestor of e . The bracket set of e is the set of all brackets of e . Definition (Cycle-equivalent, Cycle-equivalence classes) Edges e 1 and e 2 are cycle-equivalent in G U if every cycle in G U contains either both edges or neither edge. The cycle-equivalence class of an edge e in G U is the set of all edges that are cycle-equivalent to e . We now describe key properties of panbubbles. For every panbubble ⟨ s, t ⟩ in G b , the vertices in G U derived from lie along a single root-leaf path of tree in a specific order (Lemma 2). The separable condition of ⟨ s, t ⟩ further implies that edges and must be cycle-equivalent (Lemma 3). If we extract the subgraph of G b corresponding to panbubble ⟨ s, t ⟩ along with its entrance vertices and edges, then every internal black edge connects one endpoint that has a walk to s and the other that has a walk to t (Lemma 4). Furthermore, the no-hairpin condition of ⟨ s, t ⟩ ensures that no vertex v within a panbubble has a corresponding black edge in G U with an empty bracket set, where does not share a common root-leaf path with and (Lemma 5). Lemma 2. For every panbubble ⟨ s, t ⟩ in G b , either or in G U . Lemma 3. For every panbubble ⟨ s, t ⟩ in G b , edges and are cycle equivalent in G U . Lemma 4. Suppose ⟨ s, t ⟩ is a panbubble in G b . Consider walks in G b that do not pass through s or t. For every vertex u ∈ U ( s, t ), there exist such walks from s to û and from u to t, or there exist such walks from s to u and from û to t . Lemma 5. In a panbubble ⟨ s, t ⟩ in G b , there does not exist any vertex such that, in G U , (a) edge has an empty bracket set or contains only a single back edge and (b) neither nor . Next, Theorem 1 states that the above conditions are both necessary and sufficient for a subgraph of G b to represent a valid panbubble. We omit the minimality constraint here because it will be enforced in our algorithm separately by processing candidate vertex pairs in a suitably chosen order. Theorem 1. A vertex pair s, t ∈ V b satisfies all the panbubble conditions (excluding minimality) if and only if all the following conditions hold: Either or in G U , Edges and are cycle equivalent in G U , Consider walks in G b that do not pass through s or t. For all vertices u ∈ U ( s, t ) ∪ U ( t, s ), there exist such walks from s to û and from u to t, or there exist such walks from s to u and from û to t . There is no vertex such that, in G U , (a) edge has an empty bracket set or contains only a single back edge and (b) neither nor . Similarly, we state necessary and sufficient conditions for a subgraph of G b to represent a valid hairpin. Theorem 2. A vertex s is an entrance vertex of a hairpin in G b if and only if all the following hold: The bracket set of edge in G U is either empty or contains a single back edge connecting and . There does not exist any vertex which satisfies the above condition while also satisfying in G U . Consider walks in G b that do not pass through s. For all vertices u ∈ Y ( s ), there exist such walks from s to u and from s to û . Nested relationships between panbubbles and hairpins Pangenome graphs inherently support nested variation, as illustrated in Figure 3 . We characterize the conditions under which the subgraphs of G b representing (i) two distinct hairpins, (ii) two distinct panbubbles, and (iii) a hairpin and a panbubble, may share a vertex. We first prove that if two panbubbles share a vertex, one must be fully contained within the other (Theorem 3). A panbubble and a hairpin can share a vertex only if the panbubble is completely contained in the hairpin (Theorem 4). Finally, two distinct hairpins can never share a vertex (Theorem 5). Download figure Open in new tab Figure 3: An example illustrating various forms of nested relationships between panbubbles and hairpins. Three red boxes and one blue box are used to highlight three panbubbles and one hairpin, respectively. Panbubble ⟨ v 11 , v 16 ⟩ is nested within panbubble ⟨ v 9 , v 20 ⟩. Panbubble ⟨ v 1 , v 8 ⟩ is nested within hairpin ⟨ v 10 ⟩. Theorem 3. Let ⟨ x, y ⟩ and ⟨ m, n ⟩ be two different panbubbles such that B ( x, y ) ∩ B ( m, n ) ≠ ϕ, then either B ( x, y ) is a strict subset of B ( m, n ) or B ( m, n ) is a strict subset of B ( x, y ). Theorem 4. Let ⟨ x, y ⟩ and ⟨ z ⟩ be a panbubble and a hairpin respectively in G b such that B ( x, y ) ∩ ( Y ( z ) \ { z }) ≠ ϕ, then B ( x, y ) is a strict subset of Y ( z ) \ { z }. Theorem 5. Let ⟨ m ⟩ and ⟨ n ⟩ be two distinct hairpins in G b . Then, the two hairpins do not share any vertex, i . e ., ( Y ( m ) \ { m }) ∩ ( Y ( n ) \ { n }) is an empty set . 4 Proposed algorithms to detect panbubbles and hairpins 4.1 Verifying a candidate panbubble entrance vertex pair Notations and assumptions Suppose we have a candidate vertex pair ( s, t ) in G b for evaluation. We assume that preprocessing has already confirmed that ( s, t ) satisfies the first and the second condition in Theorem 1. Details of this preprocessing will be discussed later in Section 4.3 . We also assume that the following arrays have been precomputed: Boolean array A bridge of size | V b | such that A bridge [ v ] = 1 if and only if the bracket set of edge ē v is either empty or it contains only a single back edge . Recall that we obtained spanning tree via a depth-first traversal of G U . Assume that arrays A d and A f record the discovery times and the finishing times of all vertices in G U , respectively. These arrays allow efficient checks of the ancestor-descendant relationships among vertices and edges in . Algorithm We need to check for the third and fourth condition in Theorem 1. While evaluating a vertex pair ( s, t ) in G b , we use two boolean arrays W s and W t , each of size | V b |. Array W s would be used to mark those vertices that have walks originating from them to vertex s without passing through s or t . Similarly, array W t would be used to mark those vertices that have walks originating from them to vertex t without passing through s or t . We initialize arrays W s and W t with zeroes. First, we do a depth-first traversal starting from vertex s in G b . During this traversal, we visit all vertices that are reachable from s without passing through s or t . We customize the standard depth-first traversal slightly to account for the walk restrictions in G b . Consider a single iteration of the depth-first traversal. Suppose the current vertex is u ∈ V b . We perform the following operations: (i) Mark u as visited, (ii) Update W s [ û ] to 1, and (iii) If u ≠ ŝ and , we recursively visit every unvisited neighbor of û . After completion of the traversal, we repeat the same traversal starting from vertex t to compute array W t . Subsequently, we declare vertex pair s, t as satisfying all the conditions of a panbubble (excluding minimality) if both the following conditions are satisfied: (a) For all u ∈ V b , we require W s [ u ] = W t [ û ] = 1 or W s [ û ] = W t [ u ] = 1 or W s [ u ] = W s [ û ] = W t [ u ] = W t [ û ] = 0, and (b) There does not exist any vertex v ∈ V b such that (i) , (ii) A bridge [ v ] = 1, and (iii) tree edges do not share a root-leaf path in . 4.2 Verifying a candidate hairpin entrance vertex Suppose we have a candidate vertex s in G b for evaluation. We will reuse the precomputed arrays A bridge , A d , and A f defined above. We need to check for the three conditions in Theorem 2. While evaluating s , we allocate a boolean array W of size | V b |, initialized with zeroes. Array W will be used to mark those vertices that have walks originating from them to s without passing through s . To compute this array, we follow the depth-first traversal procedure customized for G b starting from vertex s . During this traversal, we visit all vertices that are reachable from s without passing through s . Consider a single iteration of the traversal. Suppose the current vertex is u ∈ V b . We perform the following operations: (i) Mark u as visited, (ii) Update W [ û ] to 1, and (iii) If u ≠ ŝ , we recursively visit every unvisited neighbor of û . After finishing the traversal, we declare vertex s as the entrance vertex of a hairpin if all the following conditions are satisfied: (a) A bridge [ s ] = 1, (b) There does not exist any vertex v ∈ V b \ { ŝ } such that A bridge [ v ] = 1 and in G U , and (c) For all u ∈ V b , we require either W [ u ] = W [ û ] = 1 or W [ u ] = W [ û ] = 0. 4.3 Complete algorithms We first construct the undirected multigraph G U from the compact biedged graph G b ( Section 3 ). We perform a depth-first traversal of G U starting from v src to obtain the spanning tree along with arrays A d and A f . To compute array A bridge , we remove self-loops and parallel edges in G U and use Tarjan’s linear-time bridge-finding algorithm [ 29 ]. Detection of panbubbles (exact algorithm) We partition edges of G U into cycle-equivalent classes by using the linear-time algorithm of Johnson et al . [ 13 ]. Next, suppose C 1 , C 2 ,…, C k are all the cycleequivalent classes. We process these classes one by one. Suppose the class being processed currently is C p , 1 ≤ p ≤ k . Suppose e 1 , e 2 ,…, e n are all the distinct black edges in G b whose corresponding edges in G U belong to class C p . We order these edges such that implies i < j . The ordering of edges can be enforced through a level order traversal of . Using a pair of nested for loops, we test the following for all 1 ≤ i < j ≤ n : (i) whether e i ≺ e j holds, and (ii) whether the vertex pair ( v i , v j ) satisfies the necessary checks discussed in Section 4.1 , assuming with and with . The outer loop iterates over i from 1 to n — 1, and the inner loop iterates over j from i +1 to n . To enforce the minimality condition, we terminate the inner loop early if a panbubble is detected. Detection of panbubbles (heuristic algorithm) We follow that same procedure as the exact algorithm, except that in the inner loop we restrict consideration to the first edge e k such that e i ≺ e k , i.e., the edge with the smallest index k satisfying this condition. The remaining edges are not considered inside the inner loop. This design choice is based on our empirical finding that panbubbles rarely include internal edges that are cycle-equivalent to their entrance edges in G U . However, adversarial cases exist where this assumption fails (Supplementary Figure S4). Detection of hairpins (exact algorithm) We iterate over all vertices v ∈ V b and check whether vertex v satisfies the necessary checks discussed in Section 4.2 . This procedure can be paired with either the exact or heuristic panbubble-detection algorithm. The following theorem states the runtimes of our algorithms. Theorem 6. Given a compact biedged graph G b = ( V b , E b , f b ), the entrance vertices of all panbubbles and hairpins can be enumerated exactly in O(| V b | 2 (| V b |+| E b |)) time. In contrast, the heuristic algorithm performs the same task in O(| V b |(| V b | + | E b |)) time . 5 Experiments Experimental setup We refer to the implementations of our exact and heuristic algorithms as Billi-exact and Billi-heuristic , respectively. Note that both exact and heuristic implementations employ different algorithms for computing panbubbles; however, the computation of hairpins is performed in the same way. We compare the performance of Billi with Pangene [ 20 ] and VG [ 24 ], which compute bibubbles and snarls, respectively. We skip comparing with Povu [ 22 ] because it did not support graphs with non-numeric vertex IDs. All our experiments were done on a dual-socket Intel Xeon Gold 6248R server with 2×24 cores and 750 GB RAM. Among the three tools, only VG supports multi-threading. Accordingly, Billi and Pangene were run using a single thread, and VG was run using 48 threads. Our commands and the tool versions are available in Supplementary Table S1. Datasets We used eight publicly available pangenome graphs of varying sizes. For each graph, we report the number of vertices, edges, and connected components under the biedged graph representation in Table 1 . Since Billi requires compact biedged graphs as input, it performs graph compaction as a preprocessing step. We also list the sizes of the resulting compacted graphs in Table 1 . The first two graphs G 1 and G 2 are pangenome graphs constructed by minigraph [ 19 ] using 90 C4A/C4B haplotypes and 57 MHC haplotypes, respectively. The next two graphs G 3 and G 4 are gene graphs constructed using pangene [ 20 ]. Graph G 3 was constructed using 50 E. coli genomes, whereas graph G 4 was constructed using 100 human haplotypes and 10 great ape haplotypes [ 20 ]. Graphs G 7 and G 8 are the first and second releases of the human pangenome graphs by the Human Pangenome Research Consortium (HPRC) [ 21 ]. These graphs were constructed using 94 and 466 assembled human haplotypes, respectively. Graphs G 5 and G 6 are subgraphs of G 8 corresponding to chromosomes Y and X, respectively. View this table: View inline View popup Download powerpoint Table 1: Sizes of the biedged pangenome graphs used in our benchmark. Results We show the count of bubbles reported by all tools in Table 2 . We also report the counts of hairpins identified by Billi. Billi-exact and Pangene were unable to process the larger graphs G 6- G 8. Nevertheless, on the first five graphs ( G 1- G 5), we observe good agreement among the tools. Notably, Billi-exact and Billi-heuristic produced identical output. This implies panbubbles tend to occur within the closest pair of cycle equivalent edges under the ≺ ordering ( Section 4.3 ). We also observe that the number of snarls reported by VG is either greater than or equal to the count of panbubbles identified by Billi, whereas the bibubble counts produced by Pangene are generally lower. Figure 4 presents Venn diagrams illustrating the concordance among the tools for graphs G 3, G 4, and G 5. View this table: View inline View popup Download powerpoint Table 2: The number of bubbles reported by Billi, VG, and Pangene. ‘–’ symbol indicates that the tool either crashed or did not finish in 24 hours. Download figure Open in new tab Figure 4: A multi-layer Venn diagram to indicate the count of bubbles detected in pangenome graphs G 3 — G 5. The figure quantifies the degree of overlap among the four bubble detection methods. By definition, a snarl is characterized only by separability and minimality [ 24 ]. Figure 5A illustrates a case in which a snarl includes vertices that do not lie on any walk from the source vertex (ID: 208092098) to the sink vertex (ID: 208092203). Because panbubbles satisfy stricter criteria, the set of panbubbles detected by Billi is expected to be a subset of snarls. However, graph G 4 presents an exception. It contains 7 panbubbles that were not reported as snarls by VG ( Figure 4 ). We further checked and found that these subgraphs agree with the formal definition of a snarl. Supplementary Figure S5 provides visualizations of two such instances. On the other hand, Pangene [ 20 ] employs a heuristic implementation that checks every pair of cycle-equivalent edges within a default distance of 100; this occasionally causes it to miss some bibubbles. Furthermore, in graph G 4, Pangene and VG identify three bubbles that Billi does not report. A closer inspection revealed that these subgraphs contain hairpins (see Figure 5B for an example). Billi reports such hairpins separately ( Table 2 ). One key difference between the bibubble and panbubble formulations is that panbubbles, by definition, do not contain hairpins. Allowing hairpins would break the guaranteed nesting structure among overlapping panbubbles (Theorem 3). Download figure Open in new tab Figure 5: Bandage visualization of two subgraphs. (A) VG marked a snarl in graph G 6 that is not a panbubble. (B) Pangene identified a bibubble containing a hairpin in graph G 4. Table 3 summarizes the runtime and memory usage of all tools. Among the evaluated methods, Billiheuristic demonstrates the best performance in both runtime and memory footprint. In contrast, Billi-exact fails to complete on graphs G 6- G 8. The size of the largest cycle-equivalent class increases dramatically from G 5 (125,267) to G 6 (1,390,824). This makes our exact algorithm prohibitively slow. In practice, bubbles are significantly smaller than the full pangenome graph (Supplementary Table S2), which benefits our heuristic algorithm. Billi-heuristic processed the largest graph G 8 within 75 minutes, whereas other methods either took significantly longer or did not finish. View this table: View inline View popup Download powerpoint Table 3: The runtime and peak memory required by various tools. ‘–’ symbol indicates that the tool either crashed or did not finish in 24 hours. Best numbers are highlighted in bold. 6 Conclusions and future work In this work, we presented new subgraph formulations, called panbubble and hairpin, for bidirected graphs that address the key points that must be respected when feeding these bubble boundaries to a variant caller. Building on these definitions, we then developed Billi, a tool that quickly identifies these structures in large pangenome graphs with hundreds of millions of vertices. As pangenome graphs continue to grow, our results offer a practical and robust algorithmic foundation for discovering novel genomic variation at scale. Several research directions emerge from this work. First, it remains an open question whether panbubbles and hairpins can be detected using a faster, linear-time algorithm. Second, our current algorithms assume that the input graph contains a degree-one vertex or tip ( Section 3 ). We leveraged this property in our proofs because such vertices are guaranteed not to lie within a panbubble or hairpin, and these vertices serve as a natural entry point for graph traversal. While the assumption holds for most pangenome graphs, we encountered one graph in the Pangene study [ 20 ] that has a component with no tips. This case suggests that relaxing the no-tip assumption will be important in the future. A third direction is to investigate how the number and structure of panbubbles evolve as additional haplotypes are augmented into a pangenome graph. Lastly, integrating panbubble and hairpin detection with pangenome-based variant calling or genotyping methods may improve accuracy and offer new biological insights. Beyond pangenomics, Billi may also be relevant for de novo assembly workflows because modern assem-blers represent unitig and contig overlaps using bidirected graphs. Detecting panbubbles and hairpins can be useful to isolate subgraphs that require further attention. Acknowledgments We thank Paul Medvedev and Heng Li for providing useful feedback. We also thank the Human Pangenome Reference Consortium (HPRC) for openly sharing their data. This research is supported in part by funding from the DBT/Wellcome Trust India Alliance (IA/I/23/2/506979). We used computing resources provided by the National Energy Research Scientific Computing Center (NERSC), USA. Funder Information Declared DBT/Wellcome Trust India Alliance , IA/I/23/2/506979 Footnotes be21b037{at}smail.iitm.ac.in , daanishm{at}iisc.ac.in ↵ * Should be regarded as joint first-authors. ↵ * In practice, pangenome graphs without tips are rare. Relaxing this assumption is left for future work. References [1]. ↵ Jasmijn A Baaijens , Paola Bonizzoni , Christina Boucher , Gianluca Della Vedova , Yuri Pirola , et al. “ Computational graph pangenomics: a tutorial on data structures and their applications ”. In: Natural computing 21 . 1 ( 2022 ), pp. 81 – 108 . OpenUrl [2]. ↵ Ljiljana Brankovic , Costas S Iliopoulos , Ritu Kundu , Manal Mohamed , Solon P Pissis , et al. “ Linear-time superbubble identification algorithm for genome assembly ”. In: Theoretical Computer Science 609 ( 2016 ), pp. 374 – 383 . OpenUrl [3]. ↵ Sai Chen , Peter Krusche , Egor Dolzhenko , Rachel M Sherman , Roman Petrovski , et al. “ Paragraph: a graph-based structural variant genotyper for short-read sequence data ”. In: Genome biology 20 . 1 ( 2019 ), p. 291 . OpenUrl [4]. HPRC Consortium . HPRC v1 . https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-chm13/hprc-v1.1-mc-chm13.gfa.gz . Online; accessed 1 Oct 2025 . 2025 . [5]. HPRC Consortium . HPRC v2 . https://human-pangenomics.s3.amazonaws.com/pangenomes/scratch/2025_02_28_minigraph_cactus/hprc-v2.0-mc-chm13/hprc-v2.0-mc-chm13.chroms/chrY.vg . Online; accessed 1 Oct 2025 . 2025 . [6]. HPRC Consortium . HPRC v2 . https://human-pangenomics.s3.amazonaws.com/pangenomes/scratch/2025_02_28_minigraph_cactus/hprc-v2.0-mc-chm13/hprc-v2.0-mc-chm13.chroms/chrX.vg . Online; accessed 1 Oct 2025 . 2025 . [7]. HPRC Consortium . HPRC v2 . https://human-pangenomics.s3.amazonaws.com/pangenomes/scratch/2025_02_28_minigraph_cactus/hprc-v2.0-mc-chm13/hprc-v2.0-mc-chm13.gfa.gz . Online; accessed 1 Oct 2025 . 2025 . [8]. ↵ Fawaz Dabbaghie , Jana Ebler , and Tobias Marschall . “ BubbleGun: enumerating bubbles and superbubbles in genome graphs ”. In: Bioinformatics 38 . 17 ( 2022 ), pp. 4217 – 4219 . OpenUrl [9]. ↵ Jack Edmonds and Ellis L Johnson . “ Matching: A well-solved class of integer linear programs ”. In: Combinatorial Optimization—Eureka, You Shrink! Papers Dedicated to Jack Edmonds 5th International Workshop Aussois, France, March 5–9, 2001 Revised Papers . Springer . 2003 , pp. 27 – 30 . [10]. ↵ Erik Garrison , Andrea Guarracino , Simon Heumos , Flavia Villani , Zhigui Bao , et al. “ Building pangenome graphs ”. In: Nature Methods ( 2024 ), pp. 1 – 5 . [11]. ↵ Fabian Gärtner , Lydia Müller , and Peter F Stadler . “ Superbubbles revisited ”. In: Algorithms for Molecular Biology 13 . 1 ( 2018 ), p. 16 . OpenUrl [12]. ↵ Glenn Hickey , Jean Monlong , Jana Ebler , Adam M Novak , Jordan M Eizenga , et al. “ Pangenome graph construction from genome alignments with Minigraph-Cactus ”. In: Nature biotechnology 42 . 4 ( 2024 ), pp. 663 – 673 . OpenUrl [13]. ↵ Richard Johnson , David Pearson , and Keshav Pingali . “ The program structure tree: Computing control regions in linear time ”. In: Proceedings of the ACM SIGPLAN 1994 conference on Programming language design and implementation . 1994 , pp. 171 – 185 . [14]. ↵ Brice Letcher , Martin Hunt , and Zamin Iqbal . “ Gramtools enables multiscale variation analysis with genome graphs ”. In: Genome biology 22 . 1 ( 2021 ), p. 259 . OpenUrl [15]. Heng Li . Pangene graphs . https://zenodo.org/records/14854301/files/Ecoli-v1.1-p0g30r5.gfa.gz . Online; accessed 1 Oct 2025 . 2025 . [16]. Heng Li . Pangene graphs . https://zenodo.org/records/14854301/files/human100p10-v1.1-a1.gfa.gz . Online; accessed 1 Oct 2025 . 2025 . [17]. Heng Li . Sample graphs and sequences for testing sequence-to-graph alignment . https://zenodo.org/records/6617246/files/C4-90.gfa.gz . Online; accessed 1 Oct 2025 . 2022 . [18]. Heng Li . Sample graphs and sequences for testing sequence-to-graph alignment . https://zenodo.org/records/6617246/files/MHC-57.gfa.gz . Online; accessed 1 Oct 2025 . 2022 . [19]. ↵ Heng Li , Xiaowen Feng , and Chong Chu . “ The design and construction of reference pangenome graphs with minigraph ”. In: Genome biology 21 ( 2020 ), pp. 1 – 19 . OpenUrl CrossRef [20]. ↵ Heng Li , Maximillian Marin , and Maha R Farhat . “ Exploring gene content with pangene graphs ”. In: Bioinformatics 40 . 7 ( 2024 ), btae456. [21]. ↵ Wen-Wei Liao , Mobin Asri , Jana Ebler , Daniel Doerr , Marina Haukness , et al. “ A draft human pangenome reference ”. In: Nature 617 . 7960 ( 2023 ), pp. 312 – 324 . OpenUrl [22]. ↵ Njagi Mwaniki , Erik Garrison , and Nadia Pisanti . “ Popping bubbles in pangenome graphs ”. In: arXiv preprint arxiv: 2410.20932 ( 2024 ). [23]. ↵ Taku Onodera , Kunihiko Sadakane , and Tetsuo Shibuya . “ Detecting superbubbles in assembly graphs ”. In: International workshop on algorithms in bioinformatics . Springer . 2013 , pp. 338 – 348 . [24]. ↵ Benedict Paten , Jordan M Eizenga , Yohei M Rosen , Adam M Novak , Erik Garrison , et al. “ Superbubbles, ultrabubbles, and cacti ”. In: Journal of Computational Biology 25 . 7 ( 2018 ), pp. 649 – 663 . OpenUrl [25]. Amatur Rahman and Paul Medvedev . “ Assembler artifacts include misassembly because of unsafe unitigs and underassembly because of bidirected graphs ”. In: Genome Research 32 . 9 ( 2022 ), pp. 1746 – 1753 . OpenUrl [26]. ↵ Jonas Andreas Sibbesen , Lasse Maretty , Danish Pan-Genome Consortium , and Anders Krogh . “ Accurate genotyping across variant classes and lengths using variant graphs ”. In: Nature genetics 50 . 7 ( 2018 ), pp. 1054 – 1059 . OpenUrl [27]. ↵ Jouni Sirén , Jean Monlong , Xian Chang , Adam M Novak , Jordan M Eizenga , et al. “ Pangenomics enables genotyping of known structural variants in 5202 diverse genomes ”. In: Science 374 . 6574 ( 2021 ), abg8871. [28]. ↵ Wing-Kin Sung , Kunihiko Sadakane , Tetsuo Shibuya , Abha Belorkar , and Iana Pyrogova . “ An O(m log m)time algorithm for detecting superbubbles ”. In: IEEE/ACM Transactions on Computational Biology and Bioinformatics 12 . 4 ( 2014 ), pp. 770 – 777 . OpenUrl [29]. ↵ R. Endre Tarjan . A Note on Finding the Bridges of a Graph. Tech. rep. UCB/ERL M427 . Feb . 1974 . url: http://www2.eecs.berkeley.edu/Pubs/TechRpts/1974/29303.html . [30]. ↵ Ryan R Wick , Mark B Schultz , Justin Zobel , and Kathryn E Holt . “ Bandage: interactive visualization of de novo genome assemblies ”. In: Bioinformatics 31 . 20 ( 2015 ), pp. 3350 – 3352 . OpenUrl View the discussion thread. Back to top Previous Next Posted November 22, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Billi: Provably Accurate and Scalable Bubble Detection in Pangenome Graphs Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Billi: Provably Accurate and Scalable Bubble Detection in Pangenome Graphs Shreeharsha G Bhat , Daanish Mahajan , Chirag Jain bioRxiv 2025.11.21.689636; doi: https://doi.org/10.1101/2025.11.21.689636 Share This Article: Copy Citation Tools Billi: Provably Accurate and Scalable Bubble Detection in Pangenome Graphs Shreeharsha G Bhat , Daanish Mahajan , Chirag Jain bioRxiv 2025.11.21.689636; doi: https://doi.org/10.1101/2025.11.21.689636 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7643) Biochemistry (17717) Bioengineering (13910) Bioinformatics (42015) Biophysics (21477) Cancer Biology (18626) Cell Biology (25536) Clinical Trials (138) Developmental Biology (13392) Ecology (19935) Epidemiology (2067) Evolutionary Biology (24356) Genetics (15617) Genomics (22530) Immunology (17755) Microbiology (40437) Molecular Biology (17200) Neuroscience (88703) Paleontology (667) Pathology (2840) Pharmacology and Toxicology (4832) Physiology (7657) Plant Biology (15171) Scientific Communication and Education (2046) Synthetic Biology (4304) Systems Biology (9828) Zoology (2272)