Construction Interface

General

GBWT is built by inserting sequences in batches to a dynamic FM-index (DynamicGBWT). The construction algorithm is based on BCR and RopeBWT2. A single step of the algorithm consists of extending each sequence in the current batch by one node.

Batch size provides a time/space trade-off. Large batches require more memory, as the sequences must be loaded into memory. On the other hand, if the sequences are aligned (to some extent), the number of records modified in each step remains small, and the algorithm is faster.

Each batch must consist of entire sequences. If batch size is specified in number of nodes and the last sequence is not fully contained in the batch, the actual batch size will be smaller than specified.

Some construction options support inserting sequences in both orientations. When using this option, the node identifiers are assumed to be already encoded with Node::encode(), so that the identifier of the reverse node can be determined with Node::reverse(). See Identifiers for further details.

Bidirectional search requires that the index contains both orientations of each sequence. This information is maintained in the header of the index.

Sample interval is a time/space trade-off for locate() queries. Given sample interval N, we store the path identifiers at one out of N positions in each sequence. If the record does not contain a stored identifier for a given offset, we iterate LF-mapping until we find a stored identifier. Sample interval 0 means that we store the identifiers at the end of each sequence.

The construction is single-threaded, though some construction methods do the actual construction in a background worker thread. The memory usage of construction is low enough that it is feasible to run multiple construction jobs (e.g. for multiple chromosomes) in parallel.

If the construction encounters an unrecoverable error, it terminates the program with std::exit(EXIT_FAILURE).

See the data model for input specifications.

Construction tools

Construction

The main construction tool is build_gbwt. It is effectively single-threaded, though it uses separate threads for reading the input and for building the index.

build_gbwt [options] input1 [input2 ...]

This builds an index for all input files output.gbwt. If no output is specified, there is only one input file, and no existing index is loaded, input1 is used as the base name for output.

• -b N: Insert the sequences in batches of N million nodes. Use batch size 0 to insert all sequences in a single batch. Default: 100.
• -c: Check for overlapping variants in generated haplotypes. Use with -p.
• -f: Insert the sequences only in forward orientation. This is the default behavior.
• -i X: Insert the sequences to an existing index X.gbwt.
• -o X: Use X as the base name for output.
• -p: The inputs are parsed VCF files. See Haplotype Generation.
• -P X: Only use phasing information from file X. May repeat. Use with -p.
• -r: Insert the sequences in both orientations. This is required for bidirectional search.
• -s N: Use sample interval N. Default: 1024. Use 0 for no samples.
• -v: Verify the correctness of the index with various queries based on the input.

Example: build_gbwt -b 200 input reads the sequences from input, builds the GBWT in batches of 200 million nodes, and writes the index to input.gbwt.

Example: build_gbwt -r -i index -o output input reads the sequences from input, inserts them in both orientations into index.gbwt, and writes the result to output.gbwt.

Merging

Existing indexes can be merged with merge_gbwt. This tool is single-threaded.

merge_gbwt [options] -o output input1 input2 [input3 ...]

The program reads input1.gbwt, inserts the sequences from input2.gbwt (and further inputs) to it, and writes the merged index to output.gbwt. For best performance, input1 should be the largest input.

• -b N: Use batch size of N sequences. If the batch size is 0, all sequences are inserted as a single batch. Default: 2000.
• -f: Use the fast merging algorithm, assuming that the node ids in the input indexes do not overlap. Ignores -b.
• -o X: Write the output to X.gbwt. Required.
• -s N: Use sample interval N in the inserted sequences. Default: 1024.

Example: merge_gbwt -o merged large small1 small2 inserts the sequences from small1.gbwt and small2.gbwt to large.gbwt and writes the merged index to merged.gbwt.

GBWT construction

DynamicGBWT is defined in dynamic_gbwt.h.

GBWT construction is based on creating an empty index with the default constructor DynamicGBWT() and then inserting sequences into it with one or more DynamicGBWT::insert() calls. There are also ways to merge existing GBWTs.

Single batch

If you can generate the sequences incrementally, this options allows building the GBWT without storing the sequences explicitly. These functions are single-threaded. They are also space-efficient, as they use the text directly as an input without additional buffering.

void insert(const text_type& text, bool has_both_orientations = false, size_type sample_interval = SAMPLE_INTERVAL)
void insert(const text_type& text, size_type text_length, bool has_both_orientations = false, size_type sample_interval = SAMPLE_INTERVAL)
void insert(const vector_type& text, bool has_both_orientations = false, size_type sample_interval = SAMPLE_INTERVAL)

• text: Batch of sequences.
• text_length: Total length of the sequences including endmarkers.
• has_both_orientations: Set to indicate that the batch contains both orientations of each sequence.
• sample_interval: Sample interval for the inserted sequences.

In a common use case, we have a single text_type as a buffer. When the next sequence no longer fits into the buffer, the batch is inserted into the index and the buffer is cleared. Because resizing text_type always causes reallocation, it is more efficient to specify the total length of the sequences instead of resizing the buffer.

Example:

DynamicGBWT gbwt;
for(size_type i = 0; i < source.size(); i += 1000)
{
  text_type batch = source.get(i, std::min(i + 1000, source.size()) - 1);
  gbwt.insert(batch);
}

This builds GBWT for source, generating batches of 1000 sequences and inserting them into an initially empty index.

Construction from disk

If the sequences are stored on disk, this option inserts them in multiple batches. The function is effectively single-threaded, though a it uses separate threads for reading the input and for building the index.

void insert(text_buffer_type& text, size_type batch_size = INSERT_BATCH_SIZE, bool both_orientations = false, size_type sample_interval = SAMPLE_INTERVAL)

• text: Sequences on disk.
• batch_size: Batch size in number of nodes. Use 0 to insert all sequences as a single batch. Default: 100 million.
• both_orientations: Set to true to index the sequences in both orientations.
• sample_interval: Sample interval for the inserted sequences.

Example:

DynamicGBWT gbwt;
text_buffer_type text(input_name);
gbwt.insert(text, 200 * MILLION);

This inserts the sequences from file input_name into an empty GBWT in batches of 200 million nodes.

Merging

This option extracts the sequences from a compressed GBWT and inserts them into the index. It is around 2 times slower than the other construction options. This function is single-threaded.

void merge(const GBWT& source, size_type batch_size = MERGE_BATCH_SIZE, size_type sample_interval = SAMPLE_INTERVAL)

• source: Input GBWT.
• batch_size: Batch size in number of sequences. Use 0 to insert all sequences as a single batch. Default: 2000.
• sample_interval: Sample interval for the inserted sequences.

Example:

DynamicGBWT input1;
sdsl::load_from_file(input1, input1_name);
GBWT input2;
sdsl::load_from_file(input2, input2_name);
input1.insert(input2);

This reads the index from file input1_name as a dynamic GBWT input1 and the index from file input2_name as a compressed GBWT input2. Then it inserts the sequences from input2 into input1.

Fast merging

When the node ids do not overlap (e.g. the indexes are for different chromosomes), the record for a particular node will be identical in the merged GBWT and the source GBWT containing that node. We can then merge the GBWTs quickly by interleaving the records. This constructor is single-threaded.

GBWT(const std::vector<GBWT>& sources)

• sources: Input GBWTs in the order the sequences will be inserted.

Example:

std::vector<GBWT> inputs(argc - 1);
for(int i = 1; i < argc; i++)
{
  std::string input_name = argv[i];
  sdsl::load_from_file(inputs[i - 1], input_name + GBWT::EXTENSION);
}
GBWT merged(inputs);

This reads the GBWT indexes for the base names specified in the command line arguments and merges them into a single index.

Buffered construction

class GBWTBuilder (defined in dynamic_gbwt.h) provides a simple interface for incremental construction. Sequences are inserted into the input buffer one-by-one. When the buffer gets full, it is swapped with the construction buffer. A worker thread is then launched to insert the sequences in the construction buffer into the index.

The build_gbwt tool and the insert(text_buffer_type& text, ...) function are based on GBWTBuilder.

The public interface consists of the following member functions. The index itself can be accessed through the index member variable.

GBWTBuilder(size_type node_width, size_type batch_size = DynamicGBWT::INSERT_BATCH_SIZE, size_type sample_interval = DynamicGBWT::SAMPLE_INTERVAL)

The constructor creates a builder containing an empty index.

• size_type node_width: Number of bits used to represent a node identifier.
• size_type batch_size: Size of the input buffer and the construction buffer. Default: 100 million.
• size_type sample_interval: Sample interval for the inserted sequences.

void swapIndex(DynamicGBWT& another_index)

This function can be used to swap the index stored in the GBWTBuilder with another index. With it, the builder can insert the sequences into an existing index. If sequences have been inserted into the buffer, finish() should be called before calling swapIndex().

• DynamicGBWT& another_index: The index to swap the contents of index with.

void insert(const vector_type& sequence, bool both_orientations = false)

Insert a single sequence into the buffer. The sequence must not contain endmarkers.

• vector_type& sequence: The sequence to be inserted.
• bool both_orientations: Set to true to insert the sequence in both orientations.

void finish()

Finish the construction. Flushes the buffers and recodes the index. Recoding sorts the outgoing edges in each record, making it possible to serialize the index (see File Formats).

Example:

GBWTBuilder builder(bit_length(Node::encode(max_id, true)));
builder.swapIndex(existing_index);
for(auto& sequence : sequences)
{
  builder.insert(sequence, true);
}
builder.finish();
builder.swapIndex(existing_index);

This example first creates a GBWTBuilder, using the width of the largest node identifier max_id in reverse orientation as node_width. It then swaps the index with an existing index and inserts all sequences from sequences in both orientations. Finally it finishes the construction and swaps the contents of the index back to existing_index.

Serialization

GBWT uses the SDSL interface for serializing and loading structures.

size_type serialize(std::ostream& out, sdsl::structure_tree_node* v = nullptr, std::string name = "") const
void load(std::istream& in)

• out: Any output stream.
• v, name: These can be ignored by the user.
• in: Any input stream.

Example: sdsl::store_to_file(index, filename); writes index to filename.

Example: GBWT index; sdsl::load_from_file(index, filename); loads index from filename.

References

Markus J. Bauer, Anthony J. Cox, and Giovanna Rosone: Lightweight algorithms for constructing and inverting the BWT of string collections. Theoretical Computer Science 483:134–148, 2013. DOI: 10.1016/j.tcs.2012.02.002

Heng Li: Fast construction of FM-index for long sequence reads. Bioinformatics 30(22):3274–3275, 2014. DOI: 10.1093/bioinformatics/btu541