Source-linked AI summary
The design and construction of reference pangenome graphs
Heng Li, Xiaowen Feng, Chong Chu
TL;DR
Integrating high-quality assemblies into a reference pangenome while preserving linear-reference coordinates remains challenging, and existing graph approaches have important scalability and representation limitations. The paper introduces a graph-based model and formats implemented in minigraph, which constructs pangenome graphs from multiple assemblies and encodes structural variation.
Problem
Single linear references cannot represent human population diversity, while existing pangenome-graph approaches lack satisfactory whole-genome construction, scalability, or long-insertion representation.
Method
The authors introduce the rGFA reference-pangenome graph format, associated GAF alignment format, and an incremental minigraph procedure for constructing graphs from multiple assemblies.
Results
Minigraph constructs a pangenome graph from twenty human assemblies in three hours and encodes structural variations of length 100bp or longer.
Takeaways & Limitations
The resulting reference pangenome graphs preserve linear-reference coordinates while compactly representing structural variation across multiple genomes.
Takeaways & Limitations
VCF-based reference-graph construction cannot properly encode variations on long insertions and is therefore limited to simple variations.
Abstract
from arXiv · showhide
The recent advances in sequencing technologies enables the assembly of individual genomes to the reference quality. How to integrate multiple genomes from the same species and to make the integrated representation accessible to biologists remain an open challenge. Here we propose a graph-based data model and associated formats to represent multiple genomes while preserving the coordinate of the linear reference genome. We implemented our ideas in the minigraph toolkit and demonstrate that we can efficiently construct a pangenome graph and compactly encode tens of thousands of structural variants missing from the current reference genome.
Background
Reference pangenome graphs address the inability of a single linear reference to represent human genetic diversity, but existing approaches have important limitations. The paper introduces rGFA and minigraph to preserve reference coordinates while incrementally constructing graphs from multiple assemblies.
- Motivation: Human reference genome GRCh38 cannot represent population diversity, leaving individuals with thousands of large germline variants absent from the reference.These variants may be associated with phenotype and can be missed or misinterpreted when short reads are mapped to GRCh38.
- Motivation: Long-read sequencing enables reference-quality individual assemblies, creating an opportunity to integrate multiple genomes into a reference pangenome.The paper describes a reference pangenome as a collection of genomes intended to address issues with a single linear reference.
- Existing approaches: Full-text indexes can retrieve genomes and test k-mers, but redundant mappings to multiple genomes can confuse downstream analyses.The text identifies unresolved multiple mappings as a practical problem for alignment-based analysis.
- Existing approaches: Colored compacted de Bruijn graphs preserve sample or population colors but discard chromosomal coordinates and may encode connections absent from the input genomes.These properties prevent colored cDBGs from serving as reference pangenome graphs.
- Existing approaches: Given hundreds of human genomes, deBGA and Pufferfish may face scaling challenges, while GBWT lacks practical tools for constructing cDBGs for many human genomes.The paper presents scalability and construction practicality as open questions for these representations.
- Paper contribution: The paper introduces rGFA and minigraph, which encode structural variations of length 100bp or longer and construct a graph from twenty human assemblies in three hours.The approach is demonstrated with human and primate graphs generated from multiple assemblies or haplotypes.
Sequence graphs
Sequence graphs represent DNA sequences on vertices and connections on bidirected edges, while GFA supplies a segment-based coordinate system. The paper extends GFA with stable coordinates so graph representations remain linked to linear reference annotations.
- Sequence graph model: A sequence graph G(V, E) is a bidirected graph whose vertices carry DNA sequences and whose edges have orientations at both endpoints.The four edge types determine how sequences are spelled along graph walks or paths.
- GFA: GFA represents vertices with S-lines and bidirected edges with L-lines, using segments and links as its terminology.The GFA grammar consists of one or more segment or link records.
- Coordinate limitation: GFA assigns each graph base a segment coordinate (segId, segOffset), but splitting a segment changes coordinates without changing graph-spelled sequences.This instability decouples the coordinate system from linear annotations and makes basic sequence graphs inadequate as reference graphs.
- rGFA: rGFA extends GFA with three tags recording each segment’s origin from a linear genome, creating a stable coordinate system anchored to the linear reference.Stable coordinates allow positions and existing annotations to be mapped onto the graph and paths to be reported in reference coordinates.
- rGFA: In rGFA, each segment has one origin, preventing collapse of different regions from one sequence and preserving a linear flavor.The format also disallows edge overlaps and multiple edges between the same vertex pair to reduce implementation complexity.
The Graphical mApping Format (GAF)
The paper introduces GAF for representing sequence-to-graph alignments and describes minigraph’s mapping and incremental graph-construction procedures. GAF preserves both graph-segment and stable-reference coordinates while minigraph prioritizes approximate mapping locations and efficient graph updates.
- The Graphical mApping Format (GAF): GAF extends PAF with a tab-delimited format whose sixth column encodes a path through the graph.Paths may use segment coordinates or equivalent stable coordinates, with optional SAM-like tags and base alignment in the cg tag.
- The Graphical mApping Format (GAF): Minigraph seeds alignments with minimizers, chains compatible hits, connects them through the graph, and selects the most weighted graph path.The mapping diagram describes a graph-chain construction that begins with minimizer-based seeds and linear chains.
- The Graphical mApping Format (GAF): Incremental construction maps each assembly to an existing graph and augments it with long sequences that map poorly.The graph is built iteratively rather than assembled independently for every genome.
- The Graphical mApping Format (GAF): Stable-coordinate paths can merge adjacent intervals from the same reference sequence and reduce to PAF for forward mappings to one reference sequence.The format supports equivalent segment- and stable-coordinate paths and recommends replacing a path derived from one reference with the entire forward reference path.
- Sequence-to-graph mapping: Minigraph adapts a minimap2-like algorithm to chromosome-long genome-to-graph alignment by finding local segment hits while ignoring graph topology.It chains local hits when they connect on the graph, possibly through cycles, to obtain approximate mapping locations.
- Sequence-to-graph mapping: Minigraph does not perform base-level alignment, which contributes to lower runtime and memory use than GraphAligner on simulated reads.The comparison used 68,857 simulated reads mapped over 8 CPU threads.
- Sequence-to-graph mapping: Minigraph was more accurate than GraphAligner, while its errors usually preserved genomic loci and GraphAligner’s errors more often selected wrong loci.The authors characterize minigraph as better at approximate locations and GraphAligner as better at disambiguating similar graph paths.
- Sequence-to-graph mapping: vg indexed the same graph in 14.7 wall-clock hours and mapped reads in 1.8 hours over 8 threads, but produced no mapped reads in the output.The authors report that vg was tens of times slower than minigraph and GraphAligner and did not work with their data.
Generating pangenome graphs
The study constructs pangenome graphs incrementally from assemblies, focusing on structural variations while enforcing filters needed for the rGFA model. Human-graph analyses reveal repeat-rich and frequently multiallelic variation, alongside assembly and alignment difficulties in complex VNTR regions.
- Generating pangenome graphs: The construction filter removes assembly regions covered by two or more primary alignments longer than 20kb.This avoids paralogous regions and ensures that minigraph-generated graphs can be modeled by rGFA.
- Generating pangenome graphs: 2.7% of 14,792 dipcall SVs were missed by minigraph, while 6.0% of 14,932 minigraph SVs were missed by dipcall.The comparison used SVs of at least 100bp and treated calls as missed when no counterpart occurred within 1000bp.
- Generating pangenome graphs: A human graph built from 20 haplotypes took 2.7 wall-clock hours over 24 CPU threads and contained 37,332 variations.The graph had 148,618 segments, 214,995 links, and a peak memory use of 98.1GB.
- A human pangenome graph: Human graph variations are enriched with Alus and VNTRs, with repeat-less variations also enriched near chromosome ends.The study additionally identified 85 processed pseudogenes among the variations.
- A human pangenome graph: Over half of VNTR variations are multiallelic, and the figure describes their distribution by repeat class, allele count, genomic location, and longest-allele length.VNTRs are defined here as tandem repeats with unit motif length ≥7bp.
- A human pangenome graph: Multiallelic VNTR regions can show insertions of different lengths and disagreement between assemblies, which the authors tend to attribute to local misassemblies.Collapsed assemblies may have more difficulty because the two haplotypes in an individual often differ at these loci.
- A human pangenome graph: VNTRs are difficult to align and call, and minigraph’s lack of base-level alignment may make optimal alignments harder in these complex regions.The authors call for improved data quality, assembly algorithms, and graph-mapping algorithms for detailed VNTR analysis.
A great ape pangenome graph
The paper presents rGFA and minigraph as a reference-pangenome framework that preserves linear-reference coordinates while representing structural variation across genomes. Applications include population-scale SV summarization, improved interpretation of short-read errors, and comparative analysis across great apes.
- A great ape pangenome graph: 206,452 variations were identified in the great ape graph, over four times more than in the human graph.The graph combined GRCh38, one chimpanzee, two gorillas, and one orangutan.
- A great ape pangenome graph: About half of great ape graph variations originated from orangutan, the species most distant from human.
- A great ape pangenome graph: The great ape graph had an L1-to-Alu ratio close to 1:1, much higher than the human graph.The paper suggests this may correlate with elevated L1 activity in great apes.
- Applications to variant calling: 21% of 32,879 false-positive SNPs fell in blacklist regions covering less than 1% of the human genome.The authors conclude that a noticeable fraction of short-read false SNP calls resulted from misalignment involving structural variants.
- Reference graph design: rGFA uses a linear reference genome as its backbone and preserves the conceptual linearity of input genomes.Unlike collapsed graph models, rGFA maintains a stable relationship between graph segments and reference coordinates.
- Implementation boundaries: Minigraph constructs graphs from 20 genomes in three hours but currently lacks sample information and base-level alignment.The method also cannot align against graphs encoding all small variants, whose many short segments impede initial seeding.
Seeding and linear chaining
Minigraph seeds alignments with minimizers and chains anchors using a gap model that permits long gaps. A modified AVL-tree RMQ accelerates chaining, while a heuristic combines thorough local searches with approximate long-range searches.
- Minigraph uses minimizers on graph segments as seeds and applies a minimap2-like chaining algorithm with different scoring.
- An anchor is a tuple (x, y, w) describing matching intervals on a reference-graph segment and the query sequence.
- The gap cost assigns infinity to gaps above G, a length-and-distance penalty to intermediate gaps, and zero cost to gapless anchors.
- G = 100000 in graph construction, compared with 5000 in minimap2; minigraph therefore permits larger gaps and penalizes them more heavily.
- Solving the chaining recurrence directly takes O(n^2), so minigraph introduces a heuristic because minimap2 can miss optimal alignments across long segmental duplications.
- Dynamic 1-dimension Range-Min-Query: RMQ stores tuples in a modified balanced AVL tree, supports dynamic insertion and deletion, and answers queries in O(m log m) time.
- Minigraph searches a small window thoroughly and a large window approximately using a faster gap-cost function.
Graph chaining
Minigraph first creates linear chains without graph topology, then performs graph-aware dynamic programming over chains connected by walks in the graph. The procedure largely uses distances rather than base sequences, producing minor effects in the resulting graphs but limiting optimal path choice.
- Minigraph generates linear chains while ignoring graph topology, then applies a second chaining round that accounts for topology.
- One chain precedes another when query coordinates increase and a graph walk connects them; among multiple walks, minigraph examines the shortest 16.
- Graph chaining computes each chain’s optimal score with dynamic programming using inter-chain weights and the linear-chain score.
- The method largely ignores base sequences when weighting chains, so similar-length walks may not yield the best base alignment; the authors report minor effects in produced graphs.
- The authors propose adding more expensive base alignment for cases where the current algorithm may fail.
- A graph chain is an ordered list of anchors whose query coordinates increase and whose graph segments are adjacent when consecutive anchors differ by segment.
The minigraph graph generation algorithm
Minigraph identifies divergent query regions from primary chains, aligns those regions at base level, and adds sufficiently large or divergent variations to the graph.
- Minigraph identifies primary chains and defines orthogonal query regions using chains longer than 100kb and non-overlap with other chains longer than 20kb.
- It scores adjacent anchor pairs by graph distance and query distance to locate subregions that differ substantially from the reference.
- The cumulative score H(i, j) becomes large for highly divergent query–graph regions, which are detected as maximal scoring intervals.
- For each divergent region, minigraph performs base alignment and retains variations involving an INDEL ≥100bp or a ≥100bp region with identity below 80%.
- Retained variations are used to augment the existing graph.
Annotating variations
Variation annotation combines RepeatMasker, SDUST, and tandem-repeat detection to classify low-complexity regions and distinguish VNTRs, STRs, and other regions.
- RepeatMasker classifies interspersed repeats in the longest allele sequence of each variation, but does not directly provide reliable VNTR or STR annotations.
- Minigraph combines RepeatMasker and SDUST results to collect low-complexity regions and separately detects pure tandem repeats with motifs occurring at least twice.
- A low-complexity region is classified as VNTR or STR when 70% of it has the corresponding annotation; remaining regions are classified as Other-LCR.
Creating blacklist regions
Blacklist regions are generated by extending each graph variation, aligning inserted sequences to GRCh38, filtering low-quality alignments, and combining interval sets defined by alignment identity.
- Each graph variation is expanded into a GRCh38 interval by 50 bp at both ends, forming I0.
- Inserted sequences are aligned to GRCh38 with minimap2 using the asm20 preset and a 2 kb alignment parameter.
- Alignments with mapping quality below 5 are filtered before defining intervals by sequence-identity ranges.
- The blacklist is computed as I0 union I(0, 0.99) excluding I(0.998, 1).