Human Readable Compression of GFA Paths Using Grammar-Based Code

In the past two decades, numerous initiatives have focused on producing reference assemblies for a wide range of species. These reference assemblies play a critical role in tasks such as genotyping, variant calling, and higher-order omics analyses, including epigenomics, metagenomics, and transcriptomics. However, a single reference sequence is insufficient to capture the full genomic diversity within a species and can introduce reference bias into the aforementioned analyses.

In contrast, pangenomes consist of sets of individual haplotype assemblies that collectively represent the genetic variability of a species. Advances in long-read sequencing and assembly technologies have recently enabled the construction of pangenomes from high-quality, reference-grade haplotype assemblies [14, 5] that are particularly suited for identifying structural variation in a species. Although pangenome applications are just emerging, several landmark studies have already demonstrated their advantages in genotyping [22, 3], variant calling [6], and transcriptome analysis [21], outperforming traditional single reference based methods.

In general, as more haplotypes are added to a pangenome, the number of novel genomic variants decreases. A key determinant of the effectiveness of pangenome-based applications is whether the pangenome sufficiently captures the genomic diversity of the species. This is typically assessed by analyzing pangenome growth curves and estimating pangenome openness [24]. These growth curves often plateau, suggesting that a finite set of haplotypes is adequate to represent the genomic variation within a species. Consequently, in the long term, we expect that numerous stable reference pangenomes will be established that will ultimately replace traditional single reference assemblies.

Graph-based data structures have become the predominant approach for representing pangenomes, as they reduce redundancy while preserving properties essential for locating subsequences and visualizing genomic variation. Two widely used graph structures are sequence graphs and de Bruijn graphs. In sequence graphs, non-overlapping genomic segments are represented as vertices, with edges denoting adjacency relationships between them. These segments are typically obtained from sequence alignments. In contrast, de Bruijn graphs are built by decomposing haplotypes into overlapping sequences of fixed length, called $q$-grams. In this model, $q$-grams form the edges of the graph, and vertices represent the $(q-1)$-length overlaps between them. In addition to the graph structure itself, pangenome graphs also retain information on the paths of the original haplotype sequences through the graph.

A popular format for storing pangenome graphs is the Graphical Fragment Assembly (GFA) format (https://github.com/GFA-spec/GFA-spec), a tabular representation in which each line constitutes a record encoding a specific type of information, such as a node (S line), an edge (L line), or a haplotype path (P or W line). The popularity of the format stems largely from its simplicity, human readability, and scripting-friendliness. However, since P and W lines store path information in an explicit, uncompressed form, as illustrated in Figures 2(a) and 2(d), file size is typically dominated by these entries that consume excessive amounts of storage, undermining the original advantages of the format. This has led to the ironic situation where, although the size of the pangenome graph itself follows the expected behavior of a saturating growth curve, the corresponding file size increases linearly with the number of haplotypes, as illustrated in Figure 1. For example, the GFA file for a recently released pangenome of human chromosome 19 comprising 1,000 haplotypes occupies 14 GB, whereas the portion encoding the graph structure accounts only for 178 MB. This discrepancy severely hampers the usability of pangenome graphs, as the resulting data volumes become increasingly unmanageable, and it poses a genuine risk to the broader adoption of pangenome-based applications.

Currently, to our knowledge only one tool is specifically designed for compressing pangenome graphs: gbz [23]. This tool constructs a variant of the positional Burrows–Wheeler Transform (PBWT) from the haplotypes encoded in a pangenome graph and also encodes the topology of the sequence graph using succinct data structures. Beyond storage, the gbz format provides efficient in-memory representations that support pangenome analysis; for instance, it underpins the giraffe read aligner for mapping short reads to pangenome graphs [22]. However, gbz only partially supports incremental updates–while modifications such as adding nodes or new haplotypes are possible, refining existing haplotypes, e.g., by replacing a short sequences of nodes with another or by splitting nodes into smaller ones, requires the rebuilding of the entire data structure. The adoption of the gbz format within the pangenome graph community has remained limited, largely due to its sophisticated data structures, binary file format, and the consequential dependence on external libraries for construction and usage.

In this work, we propose an extension of the more popular GFA file format that allows haplotype paths to be specified in a compressed form, while preserving human readability. Our approach is based on context-free grammars. To this end, we introduce a new line type, the Q line, to define meta-nodes representing paths in the graph. This enables the representation of repetitive substrings of haplotype paths of length two or more with a single meta-node, resulting in a more compact path encoding. To reference these grammar-defined subpaths within haplotype paths, we introduce an additional line type, the Z line. We also present sqz, a proof-of-concept implementation, written in Rust, available at https://github.com/codialab/sqz.

There is a long-standing body of research on the concept of grammar-based code [9]. In grammar-based coding, nonterminal symbols are introduced to replace subsequences of the original text that occur multiple times. Each nonterminal is associated with a production rule that maps it to a sequence of symbols–potentially including other nonterminals–that ultimately spells the original subsequence. A prominent grammar-based coding scheme is byte pair encoding (BPE), in which each production rule corresponds to a digram, i.e., a pair of adjacent symbols, that appears more than once in the text. To construct the grammar, BPE iteratively identifies the most frequent digram, replaces all of its occurrences with a new nonterminal symbol, and records the corresponding production rule. This process continues until no digram appears more than once in the transformed text. BPE and its variants, such as SentencePiece [12], have become standard in machine learning, particularly in large language model (LLM) applications, where compact and consistent tokenization is crucial for managing vocabulary size and handling rare or unseen words efficiently [19, 18].

We implemented our method, sqz, in Rust and released the source code under the MIT license at https://github.com/codialab/sqz.

We evaluated sqz on three human pangenome graphs:

• $\mathrm{■}$

  Chr19, comprising 1,000 haplotypes [8],

• $\mathrm{■}$

  HPRC v2.0 MC GRCh38 full and

• $\mathrm{■}$

  HPRC v2.0 MC GRCh38 clipped, both comprising 464 haplotypes [14].

The full version of the HPRC v2.0 MC GRCh38 graph includes complete assemblies, but alignment quality is poor in certain regions. Particularly centromeres contain essentially large chunks of unaligned sequence. In contrast, the clipped version removes these segments, limiting the pangenome graph to regions with reliable alignments.

We benchmarked our method against three compression tools: bgzip [1], gbz [23], and sequitur [15]. bgzip produces gzip-compatible compressed binary files and is widely used for compressing pangenome graphs in the GFA format, primarily because its output can be decompressed with gzip, which is pre-installed on most Unix systems. Its key advantage over gzip is support for indexing, allowing selective decompression of file segments. gbz, in contrast, is specifically designed for compressing pangenome graphs and integrates with tools in the GBWT ecosystem, as discussed in Section 1. sequitur is a classical grammar-based compression algorithm in natural language processing (NLP).

Among existing BPE based algorithms, Re-Pair [13] is most similar to our approach. Both algorithms iteratively replace the most frequent digram with a production rule or meta-node, proceeding in descending order of frequency. However, Re-Pair operates on linear text rather than graph structures, and it does not account for inverted digrams, unlike sqz and gbz. Although several Re-Pair implementations are available, many operate solely on text or binary files and are not designed to handle characters spanning multiple bytes. Therefore, we were unable to successfully apply them to haplotype paths, which are composed of node identifiers that cannot be represented by a single byte.

In contrast, we were able to evaluate sequitur [15], another BPE algorithm. sequitur constructs and applies production rules as soon as it reads a repeated digram. As a result, sequitur’s memory usage scales with the size of the compressed output rather than the input. However, its compression performance depends on the input order, since it processes digrams in the order they appear rather than by frequency. sequitur allows the user to set the size of its core data structure, a hash table. This affects both memory usage and runtime: larger tables reduce collisions and improve performance at the cost of increased memory. For all evaluations, we set the hash table size to 2,000 MB, which was the maximum value supported by the program.

As an initial experiment, we used the Chr19 pangenome to evaluate how compression scales with the number of haplotypes in a pangenome graph. To this end, we randomly removed haplotypes until only a specified number remained, and removed all S- and L-lines corresponding to nodes and edges that are not covered by the remaining subset. We then applied sqz, bgzip, sequitur, and gbz to the resulting trimmed graphs. Additionally, we tested a combined approach using sqz followed by bgzip compression, as bgzip is commonly used in transferring GFA files. The same thing we also did for gbz.

The resulting file sizes and compression ratios are shown in Figure 5. The compression ratio is defined as $\text{uncompressed}/\text{compressed}$. bgzip exhibits constant compression behavior: file size increases linearly with the number of haplotypes, and the compression ratio remains constant at about a rate of $4$. For fewer than 100 haplotypes, bgzip achieves better compression than the purely grammar-based approaches sqz and sequitur. For all other methods, the compression ratio improves as more haplotypes are added, benefiting from increased redundancy. Up to 100 haplotypes, the difference between sqz and sequitur remains small (e.g., $0.19$ GB vs. $0.24$ GB), but it gets larger with increasing number of haplotypes. Among all methods, gbz+bgzip delivers the best compression results. sqz+bgzip slightly outperforms gbz for fewer haplotypes.

A second benchmark was conducted on the clipped and full HPRC v2.0 MC GRCh38 graphs, that have been split up into their individual chromosome graphs. This setup enables an evaluation of how compression performance is affected by the similarity among haplotypes, due to the full graphs containing less similar haplotype paths than their clipped counterparts.

Due to the large input sizes (up to 36 GB), we were unable to run sequitur for this benchmark, as it required excessive runtimes. We evaluated compression performance in terms of both file size and compression ratio. For the clipped graphs, sqz+bgzip and gbz perform similarly again, achieving compression ratios between 60 and 80, depending on the specific graph. sqz alone achieves a compression ratio of around 20, while bgzip performs worst overall. The best result achieves gbz+bgzip, reaching compression ratios of up to 150. Except for bgzip, compression results on the full graphs vary considerably, as the amount of heterochromatic sequence differs across chromosomes. For chromosomes 1 through 11, bgzip typically outperforms sqz, but from chromosome 12 onward, sqz yields smaller files, with the exception of chromosomes 1 and 9. gbz consistently produces smaller files than bgzip, while the combined gbz+bgzip approach generally results in the smallest file sizes. With the exception of chromosome 9, sqz+bgzip is always very close in terms of resulting file size.

To assess how effectively haplotypes are compressed using sqz, we examined the number of nodes per haplotype for each chromosome in the clipped graphs, illustrated in Figure 7. sqz achieves a reduction in haplotype length by more than two orders of magnitude. This analysis also reveals a non-negligible variance across haplotypes, with some chromosomes, e.g., chromosome 15, containing haplotypes that are much less compressed than others.

Another key consideration when working with large files is the runtime and memory consumption of the compression algorithm. We compared the performance of sequitur, sqz, gbz and bgzip on the Chr19 pangenome using the set of trimmed graphs containing a varying numbers of haplotypes that we constructed for the first experiment.

In terms of runtime, sequitur outperforms sqz on graphs with up to 800 haplotypes. However, for the two largest graphs with 900 and 1,000 haplotypes, sqz is faster and, for the latter graph, requires more than an hour less runtime. The fastest method, however, is gbz (with the exception of the smallest two graphs, where bgzip is equally as fast or faster). On memory usage, bgzip performs the best by having a constantly very low memory usage. Still, sequitur and gbz outperform sqz by a wide margin. Even for the largest graph, sequitur stays below 5 GB of memory usage and gbz below 20 GB, while sqz requires up to 176 GB.

To evaluate the performance of computing coverage tables using the in-memory representation of the grammar, we implemented two Python scripts that make use of libraries networkx and numpy. One does the calculation on an uncompressed GFA file of the Chr19 pangenome. The second one follows the approach described in the previous section and takes a GFA file compressed by sqz. The script for the uncompressed GFA file takes 7:49 min, while the one on the sqz-compressed graph only takes 2:39 min. The memory usage is reduced over ten-fold when using the script with the compressed graph, requiring only 14.8 GB instead of 159.6 GB. We also ran the two approaches over all of the HPRC v2.0 graphs (clipped and full), the results are shown in Figure 8. Similar to the Chr19 graph, working with the grammar-based strategy is faster and uses less memory for all chromosomes. Particularly, the difference in terms of runtime is large, with the original method needing 2:08 hrs, while the compressed files only need 12 min for the clipped chromosome 1.

By introducing Q- and Z-records, we propose an extension to the GFA format that supports human-readable, grammar-based compression. Grammar-based compression is simple and intuitive, yet is as effective as any finite-state compression scheme [10]. It is highly versatile and has given rise to a broad body of research, resulting in a variety of compression schemes [9]. Our proposed GFA extension fully supports incremental updates: nodes and edges can be added to the graph without interfering with the grammar captured in Q-records. Likewise, new haplotypes can be encoded using the existing grammar, although this may lead to less efficient compression. In addition, common pangenome graph modifications, such as “unchopping”, i.e., the merging of a universal digram into a single node or the decomposition of a node into a series of smaller ones, affect the grammar only locally such that at most a few low-level meta-nodes need to be updated. Most importantly, we argue that grammar-encoded haplotypes need not be decompressed for analysis. Instead, they enable faster and more memory-efficient haplotype analysis, without the need to use dedicated external libraries. The simplicity of the grammar allows it to be represented by data structures available in widely used scripting languages such as Python.

Our BPE-based implementation, sqz, serves as a proof of concept for the proposed GFA extension. sqz often achieves higher compression ratios than bgzip and sequitur, and performs on par with gbz when combined with bgzip. The latter combination is only outperformed by gbz+bgzip that provides an even higher compression ratio. However, sqz offers an advantage over gbz in that it provides a simple, human-readable, and easily modifiable encoding.

We envision further variants of grammar based encodings that can be realized with our proposed GFA extension, such as grammars based on $q$-grams for $q>2$, or grammars derived from maximal exact matches (MEMs) which are widely used in pangenomics [20]. Moreover, Q-records can carry biological meaning by annotating shared sequences associated with specific alleles. This enables not only compression, but also allows to encode haplotypes as combinations of alleles, embodying a natural, biologically motivated representation.

Conceptually, our proposal represents a dichotomy to existing approaches in pangenomics that encode differences between haplotypes, such as snarl decomposition [17, 2] or reference-based variant formats like the variant call format (VCF). We argue that capturing commonalities rather than differences offers several advantages. First, decompression does not require an external reference or computationally intensive edit operations to reconstruct the original sequence. Second, there is often no unique or universally accepted way to represent differences, leading to arbitrary decisions that impede comparability across datasets, a challenge exemplified by the persistent difficulties in VCF merging [11, 4]. Nevertheless, some of these representations could potentially be reconciled with our grammar-based approach. For example, a snarl decomposition can naturally give rise to a simple grammar, where each tip and each shared block within a chain, as well as each distinct variant within a snarl, is represented as a meta-node.