Skip to content

Using the GBWT in stoat

Jump to bottom
Xian Chang edited this page Sep 30, 2026 · 8 revisions

The GBWT is an index of the haplotype paths through a variation graph. Accessing paths in the GBWT is faster than using other graph representations such as the hg or pg but the interface is a bit more difficult to understand. Here I'll try to explain my understanding of the GBWT interface and how we can use it for STOAT. My usage of the GBWT is based heavily on vg's haplotype_extracter

Background

The GBWT is based on the Burrows Wheeler Transform/FM-index/r-index. This page only explains the GBWT so I recommend watching Ben Langmead's Burrows Wheeler Indexing series, at least the videos on the suffix array (in the suffix indexing playlist), BWT, FM index, RLFM, and r-index. A lot of the documentation about the GBWT interface uses vocabulary from the FM index so it's useful to understand.

Data structures

The GBWT library includes a run-length encoded BWT, FM-index, and (optionally) r-index but I'm just going to explain the parts I think are necessary to understand the interface.

Papers:

GBWT

The GBWT is essentially an FM index of the concatenated haplotypes (encoded as a sequence of node ids, with 0 as the end marker, and actually stored in reverse order), with the BWT split up into nodes. For each node, the GBWT stores the subset of the BWT that encodes the reverse prefixes ending at that node, meaning that for each haplotype traversing a given node, the node record stores the next node in the haplotype.

image

Figure 1 of the GBWT paper shows the GBWT and the node records for each node. Each node v's BWT has a local alphabet Σv , comprising of the nodes w for which there exists an edge (v, w), which is encoded using integers starting from 0. For example, node 4's alphabet would be {5, 6}, which are encoded as 0 and 1. The top part of the node record shows, for node v, for each subsequent node w, the pair (w, BWT.rank(v,w)), where BWT.rank(v,w) is the number of occurrences of the edge (v',w) where v' comes before v in the ordering of nodes. That is, for all nodes before v, how many have edges to w. For node 4, the two outgoing edges go to 5 and 6. Node 5 has one edge entering it that comes from a node with id less than 4.

The bottom part shows the subset of the BWT stored for this node, using the node's local encoding of the alphabet, Σv, rather than the node id itself. In a regular BWT, the BWT is sorted by the right context of each character. Since a GBWT is a BWT of the reverse paths, the BWT is sorted by the left context (the reverse prefix ending at a given node), meaning that the haplotypes in a node are ordered by where they came from (I think). For example, the paths going into node 7 from node 5 are consecutive in 7's BWT, since they have the same reverse prefix (>7>5). Since nodes that come from the same incoming node tend to take the same outgoing node, we will have runs of haplotypes that have the same outgoing node. The BWT for a node is encoded as a pair of outgoing node and the number of haplotypes in the run that take the same next node (internally but maybe that doesn't matter for the interface).

In this way, a position on a haplotype can be identified by its node id, and the offset in the list of haplotypes stored on the node. A range of offsets in the list of haplotypes refers to a group of haplotypes that are consecutive in the BWT. Since the haplotypes are ordered, consecutive haplotypes generally share a similar path leading to the given node.

r-index

The r-index is used to answer the locate() query using a small amount of memory. locate() is used to find all occurrences of a pattern in the original text. This is done by storing a subset of the suffix array, specifically the first (and last?) occurrences of a character in a run of characters in the BWT. This lets us go from an entry in the BWT to its location in the original text.

GBWT interface

The GBWT query interface and search interface are described in the wiki but I'll try to explain them as I understand it and as I used it in stoat

SearchState

A search state represents a node and a set of haplotypes that traverse that node, represented as a range of offsets in the BWT of that node. Since the haplotypes traversing a node are sorted by the incoming edges, a range can represent the set of haplotypes that took the same path to get to this node.

Useful functions:

  • find(): Given a node, get the associated search state covering all haplotypes through the node
  • extend(): A search state can be extend()ed to a next node walking out from the current node, which will return a new search state pointing to the next node, with the range corresponding to the subset of the original range that followed this edge to the next node. Thus running a series of extend()s from a node along a walk will gradually narrow down the range of haplotypes to those that followed this walk.
  • locate(): For a regular GBWT, locate() finds the path identifiers of the search state. Does not find the actual location of the node in a haplotype path, just which haplotype it's on
  • LF(): Given a position (a node on a path), get the next node on the path. Offset i identifies the path as an offset into the node's list of paths

FastLocate (r-index)

This is how I think the r-index works.

Functions:

  • find(): This works as the GBWT's find() to get a search state from a node, and additionally finds the entry in the Suffix Array corresponding to the first occurrence in the range of the search state. This suffix array entry is what is used to go from a location in the BWT back to a location in the original text (the haplotypes).
  • extend(search_state, first): As above, extends the GBWT's extend() to find the first entry in the Suffix Array (first). Note that first is a reference and its value must be the previous suffix array entry found using find() or extend(), and this will update its value to be the new suffix array entry.
  • seqId(): get the path identifier from a suffix array value
  • seqOffset(): get the sequence offset from a suffix array value
    • How it works: The suffix array value is an index into the original text, which is all the path sequences concatinated so since they are forced to have the same length, this seqID() and seqOffset() can get the id and offset along the path using div and mod.
  • locateNext(): Given a value in the suffix array, get the next offset in the interval. As far as I can tell, a contiguous range in the BWT of a node will be contiguous in the suffix array so we can use this to get everything in the search state's range. Running locateNext() for a search state's range is the same as decompressSA() (I think).

Getting path information from an identifier

A path identifier can be found using the locate() function with the GBWT, or with the seqId() function with FastLocate. The metadata for each path identifier is stored in the GBWT's Metadata class. The metadata can be retrieved using functions like get_path_sense() and get_path_sample_name() defined by the gbwtgraph.

Here is how we get the path name from the GBWT's path identifier (based on how they do it in vg):

handlegraph::PathSense sense = gbwtgraph::get_path_sense(gbwt, gbwt_path_id, gbwt_reference_samples);
std::string path_name = handlegraph::PathMetadata::create_path_name(sense,
                                                                    gbwtgraph::get_path_sample_name(gbwt, gbwt_path_id, sense),
                                                                    gbwtgraph::get_path_locus_name(gbwt, gbwt_path_id, sense),
                                                                    gbwtgraph::get_path_haplotype(gbwt, gbwt_path_id, sense),
                                                                    gbwtgraph::get_path_phase_block(gbwt, gbwt_path_id, sense),
                                                                    gbwtgraph::get_path_subrange(gbwt, gbwt_path_id, sense));

Clone this wiki locally