Repository navigation
Internals
This page describes implementation details relevant to developers who want to understand or modify BubbleFinder's internals.
BubbleFinder's algorithms are based on the following papers.
The core framework comes from Sena, Politov et al., 2025, which introduces an SPQR-tree approach to find all snarls and superbubbles in linear time.
For ultrabubbles, we rely on a forthcoming preprint (Scalable computation of ultrabubbles in pangenomes by orienting bidirected graphs). It describes how to orient a bidirected graph using DFS, how to resolve conflicts with auxiliary tips, and how to map weak superbubbles back to ultrabubbles.
The weak superbubble detection step uses the linear-time algorithm from Gärtner & Stadler, 2019. BubbleFinder runs this on the oriented skeleton to obtain ultrabubbles.
Finally, the GBZ input format is defined in Sirén & Paten, 2022. BubbleFinder supports both GBZ and GFA as input formats.
If you look at example/tiny1.png, you will notice that the bidirected edge {a+, b+} from the graph image is encoded as:
L a + b - 0M
This is because GFA links are directed. The rule is: for every link
a x b y
in the GFA file (where x, y ∈ {+, -}), we flip the second sign y into ¬y, and build a bidirected edge {a x, b ¬y}. Snarls are then found in this bidirected graph.
Superbubbles are defined on directed graphs, but GFA graphs are bidirected. In superbubbles mode (bidirected, the default), BubbleFinder first converts the bidirected graph into a doubled directed graph: each segment v gets two oriented copies v+ and v-, and superbubbles are detected between oriented endpoints (e.g. (a+, e-) and its mirror (e+, a-)).
To report results at the segment level, without depending on the orientation chosen in the doubled graph, we apply an orientation projection. In practice, this means two things: mirror pairs like (a+, e-) and (e+, a-) count as the same bubble, and we then drop the +/- signs and sort the two segment names.
The final output is a single unordered pair of segment IDs, e.g. a e. The figure below illustrates this process.
In this example, the directed graph on segments has three superbubbles with endpoints (a, b), (b, e) and (e, f). After running the superbubble algorithm on the doubled graph and applying the orientation projection, BubbleFinder reports exactly these three pairs.
--compact-output-chains is an output compaction pass. It looks for linear runs of consecutive endpoint pairs and writes only the two external endpoints of each run.
Example without compaction:
3
a+ b-
b+ c-
c+ d-
With --compact-output-chains:
1
a+ d-
The internal joins are matched by opposite endpoint orientations: b- joins b+, then c- joins c+.
flowchart LR
P1["a+ b-"] -->|"b- / b+"| P2["b+ c-"]
P2 -->|"c- / c+"| P3["c+ d-"]
P1 -. compacted .-> C["a+ d-"]
P3 -. compacted .-> C
Only unambiguous chains are compacted. If an internal endpoint can be joined in more than one way, BubbleFinder keeps the original pairs instead of choosing an arbitrary chain.
flowchart LR
A["a+ b-"]
B["b+ c-"]
X["b+ x-"]
A -->|"b- / b+"| B
A -->|"b- / b+"| X
Amb["ambiguous: b+ appears twice"]
B --> Amb
X --> Amb
Use the normal output if downstream code needs every adjacent pair explicitly.
The snarls solver works on biconnected blocks and their SPQR trees. Large pangenome blocks often contain long degree-2 regions. SPQR compression keeps those regions as macro records and builds the SPQR tree on the remaining core graph.
Here is a schematic block before compression. Dark nodes are kept as core nodes because they are branch points. Pale nodes are degree-2 nodes that can be represented by macros.
graph LR
U((u)) --- X1((x1)) --- X2((x2)) --- V((v))
U --- Y1((y1)) --- V
V --- Z1((z1)) --- Z2((z2)) --- W((w))
U --- W
classDef core fill:#ffffff,stroke:#111111,stroke-width:2px;
classDef deg2 fill:#eeeeee,stroke:#888888;
class U,V,W core;
class X1,X2,Y1,Z1,Z2 deg2;
The core graph only contains graph vertices, drawn as circles. A compressed region is not a new graph vertex: it is an edge label such as m1, pointing to a macro record. The edge u--w is an original edge, so it has no macro record.
graph LR
U((u)) ---|"m1"| V((v))
V ---|"m3"| W((w))
U ---|"original edge"| W
classDef core fill:#ffffff,stroke:#111111,stroke-width:2px;
class U,V,W core;
The macro tree stores what each macro record means. Its boxes are internal records, not extra vertices in the core graph. A series macro is therefore a record in this tree, not a graph node. A parallel macro groups alternative branches between the same poles; a series macro stores one path through degree-2 nodes.
graph TD
M1["m1: parallel macro<br>poles: u, v"]
M1a["series child<br>u - x1 - x2 - v"]
M1b["series child<br>u - y1 - v"]
M3["m3: series macro<br>v - z1 - z2 - w"]
M1 --> M1a
M1 --> M1b
classDef macro fill:#f6f6f6,stroke:#666666;
class M1,M1a,M1b,M3 macro;
The two u-v branches are represented by m1, whose children are two series records. The v - z1 - z2 - w path is independent in this example, so it is represented by m3. In the implementation, macros and core nodes are stored by compact integer IDs; node names are only materialized when output is written.
The main --sp-compress modes are:
| Mode | Meaning |
|---|---|
off |
Build the full SPQR tree without compression. Useful as a baseline. |
on |
Use compressed SPQR construction, then expose a regular SPQR tree to the solver. |
instrument |
Compare baseline and compressed construction, emit timing/statistics data, then keep the baseline tree for solving. |
macro-direct |
Use the macro tree and core SPQR tree directly in the snarl solver, avoiding full reconstruction for compressed blocks. |
In macro-direct, the solver keeps DP state for both the macro nodes and the core SPQR nodes. If a block has no useful compressed structure, BubbleFinder falls back to the regular SPQR path for that block.