Panpath-Audit: Sequence-Level Verification of Embedded Pangenome Paths Against Source Assemblies

1. Introduction

Pangenome graphs represent the sequence diversity of a collection while retaining individual genomes or haplotypes as paths through the graph [1]. This model underlies resources such as the draft human pangenome [2] and graph-building systems including Minigraph–Cactus [3] and PGGB [4]. Graphical Fragment Assembly (GFA) files commonly encode segment sequences in S records and embedded traversals in P paths or W walks [5].

Graph topology and source preservation are related but different properties. A topology validator can establish that path steps are present and connected; for example, odgi validate checks path consistency with graph topology [6]. It cannot establish that an independently distributed FASTA collection is the collection represented by the graph, nor that a segment sequence remained unchanged during conversion or transfer. A structurally valid graph may therefore reconstruct a different nucleotide sequence while continuing to satisfy its internal constraints.

panpath-audit addresses the narrow question: did the embedded graph traversals preserve their named source sequences? It is builder-independent, operates directly on plain or gzip-compressed FASTA and GFA files, and reports both complete correspondence and sequence identity. Here we describe its audit contract and apply it to plant and human public pangenomes.

2. Implementation

2.1. Audit Contract

Source identifiers are paired case-sensitively with graph traversal names. Inputs may use exact names, arbitrary prefixes, or the PanSN convention sample#haplotype#contig [7]. Nucleotide comparison is case-insensitive but symbol-exact over the DNA IUPAC alphabet; for example, N does not match A. Reverse-oriented segments use the complete IUPAC complement mapping.

The supported graph contract comprises inline-sequence S records, blunt P paths with unspecified or 0M overlaps, and coordinate-bearing W walks. Non-overlapping walk records with the same PanSN name are combined in coordinate order. Source intervals between explicit walk ranges are classified as not embedded rather than divergent. The program does not validate links, graph topology, external segment references, or non-blunt overlap CIGARs.

Before comparison, an all-or-nothing preflight rejects invalid FASTA symbols, duplicate source names, missing segment sequences, unresolved steps, duplicate traversals, ambiguous walk ranges, and unsupported overlaps. A completed audit classifies every identifier as identical, divergent, missing path, or missing source (Figure 1). Divergence reports the first one-based source and traversal positions and, when present, the segment, orientation, traversal position, and within-segment position. Stable exit codes distinguish identity, divergence, correspondence failure, unsupported GFA, and operational failure.

2.2. Streaming and Reporting

FASTA inputs are streamed twice: once to index mapped identifiers and lengths, and once to supply source sequences to a bounded worker pool. Graph records and segment sequences are indexed in anonymous temporary files instead of retaining the complete graph in memory. A shared admission budget bounds concurrently held sequence bytes; a contig larger than the budget is admitted alone to guarantee progress. This accounting intentionally excludes allocator overhead, compression buffers, parser metadata, operating-system caches, and temporary storage.

JSON and TSV reports include SHA-256 and BLAKE3 digests of complete sources, addressed source ranges, and reconstructed traversals. Comprehensive mode adds separate source and graph base ledgers and bounded exact unit-cost Levenshtein alignment for divergent traversals. First-divergence diagnostics and digests remain available when the alignment work or memory limit is exceeded.

3. Evaluation

3.1. Datasets and Controls

The tomato dataset followed the 23-assembly preparation used to evaluate PGGB [4], originally derived from the tomato graph pangenome resource of Zhou et al. [8]. Chromosome 2 sequences were renamed with PanSN identifiers and audited against the cleaned PGGB GFA containing 23 P paths. Two prepared FASTA controls retained the same graph: one substituted a single base at the midpoint of the 55,739,602 bp MM chromosome, and one omitted the SL5 source entirely.

For the human evaluation, 94 phased Year-1 HPRC assemblies with existing PanSN headers were combined with CHM13 using a PanSN mapping, then audited against the public GRCh38-based Minigraph–Cactus v1.0 clipped and v1.1 full graphs [2,3]. The two graphs name their GRCh38 traversals differently: the clipped graph embeds GRCh38 as P paths named GRCh38.chr1, whereas the full graph embeds it as PanSN walks named GRCh38#0#chr1. GRCh38 was therefore mapped by graph path prefix for the clipped graph and by PanSN sample and haplotype for the full graph. Applying one convention to both graphs reclassifies every GRCh38 contig as missing instead of comparing it, which is itself an instance of the correspondence failure the audit is designed to expose. The clipped graph uses explicit walk ranges and can omit source intervals; the full graph preserves untrimmed walks. Publisher-provided assembly checksums were verified before analysis.

3.2. Execution and Results

Version 0.1.2 of panpath-audit [9] was compiled in release mode with Rust 1.98.0 and run with 24 worker threads, a 6144-MiB tracked-memory limit, and comprehensive statistics. Experiments ran on Linux with 24 logical processors and 61.7 GiB host memory; inputs were read from btrfs storage and scratch files were written to btrfs storage. Wall-clock time and maximum resident set size in Table 1 are therefore descriptive telemetry, not a cross-system performance benchmark. The complete evaluation was subsequently rerun on unrelated hardware with a different processor count, operating system, kernel, filesystem, and Rust minor version against byte-identical inputs: every outcome classification, base ledger, and digest was unchanged, and only the timing and memory telemetry differed.

All 23 tomato paths were identical across 1,280,460,312 reconstructed bases. Against the altered source, panpath-audit classified MM#1#chr02 as divergent and the remaining 22 as identical, localizing the substituted base to source position 27,869,801, path position 27,869,801, and offset 69 within segment 1120821. Recovering one altered base from the interior of a 55,739,602 bp chromosome exercises the linear and graph coordinate reports together, including the position within the segment rather than only the segment itself. Against the source set lacking SL5, it reported SL5#1#ch02 as missing source rather than treating the remaining 22 identities as complete correspondence.

In the HPRC clipped graph, 22,923 traversals were identical. The source base ledger partitions 289,594,475,134 bases into 87.99% matched against an embedded traversal, 7.98% belonging to the 17,991 contigs the graph does not embed at all, and 4.03% (11,671,173,035 bases) falling outside the explicit walk ranges of contigs that are embedded. The missing paths are short assembly contigs that Minigraph–Cactus does not place; they are reported as absent correspondence rather than as sequence loss. The full graph embedded 19,798 traversals identically, reported 0 not-embedded bases as expected for untrimmed walks, and left 8.93% of source bases in unplaced contigs. Missing sources numbered 42 and 43 for the clipped and full graphs.

Both audits reported the same 14 divergent traversals, the GRCh38 chromosomes, each carrying 94 substituted bases in total and the same first-divergence source coordinates, even though the two graphs were constructed separately and encode GRCh38 as P paths and W walks respectively.

4. Discussion

These experiments distinguish three properties that are easily conflated: internal graph consistency, complete source–traversal correspondence, and exact sequence recovery. The tomato controls show why reporting only successfully matched paths is insufficient: one changed source or one omitted source must make the complete audit fail. The HPRC analysis additionally shows why coordinate-aware gaps should not be labelled as sequence substitutions.

The 14 divergent GRCh38 chromosomes have a single benign explanation. GRCh38 carries 94 IUPAC ambiguity codes outside the ACGTN alphabet, distributed across exactly those 14 chromosomes and no others, and graph construction replaced each of them with N. Per-chromosome ambiguity-code and substituted-base counts agree exactly in every case, the graph symbol is N at every reported position, and in 13 of the 14 chromosomes the replacement occupies a dedicated single-base segment. Because the comparison is symbol-exact, N does not match R, Y, or M, so the substitution is reported instead of silently absorbed. Recovering the same 94 positions from two independently constructed graphs, through both P paths and W walks, indicates a property of the construction pipeline rather than of one distributed file. Reporting it is the intended behaviour: the audit states what changed and leaves the significance to the reader, which is the distinction a topology validator cannot draw.

The tool is deliberately narrower than a general GFA validator. Its guarantee applies only after the documented input contract has passed, and it does not infer biological homology, validate topology, or support sequence-less segments and non-blunt path overlaps. Exact IUPAC equality is appropriate for artifact integrity but should not be interpreted as variant calling or assembly-quality assessment. Memory telemetry is a tracked admission bound rather than total process memory. These limitations are exposed in reports to prevent an identity result from being interpreted more broadly than warranted.

The resulting binary outcome, structured diagnostics, and stable exit codes make the audit suitable after graph construction or conversion, before data release, and as a continuous-integration gate. Topology validation remains a complementary step.